PROGRAM CCTN IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind PARAMETER (MX=4000,MXB=140,IOX=1500,n30=170,nel=118) c MX maximum number of atoms c MXB maximum number of overlaped atoms for each constant c IOX max # of overlaps c n30- dimension for the cubic and quartic constants c number of polarizability components for sos real*8 R(3,MX),RS(3,MX),AM(MX,3,3),alphf,QS(MX,3,3,3), 3PS(MX,3,3),AS(MX,3,3),PB(MX,3,3),AB(MX,3,3),ALPHA(3*MX,3,3), 1ALPHAS(3*MX,3,3),AA(3*MX,3,3,3),AAS(3*MX,3,3,3),a33(3*MX), 1G(3*MX,3,3),GS(3*MX,3,3),AAT(3*MX,3,3,3),GT(3*MX,3,3), 1QB(MX,3,3,3) integer*4 IND(MX),IB(MX,MXB),ND(MX),NDS(MX),IS(MX,MXB),iweights, 2INDBIG(IOX,MX),IBONDSB(4,MX),NBO(MX),NBOB(MX),j33(3*MX) LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM,LXX3,LXX4, 2L36,L37,lxfile,lpol,lpolg,lpola,lq,LDATOM,lsos,lmut,lmpt,lgrad, 3sgrad,fgrad,lnmr,lspin,LQUAD,SOSNM CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname,pname CHARACTER*8 sss,ssc CHARACTER*6 ssb DIMENSION NDNAME(IOX,MX),coo(3) c DIMENSION NDNAME(IOX,MX),cb(81*n30**3),cs(81*n30**3) REAL*8,allocatable:: cb(:),cs(:),ct(:),FFS(:),alsos(:),alsosb(:), 1weights(:),FF(:,:),FT(:,:),GG(:),NMRNS(:),SPINNS(:),SPINT(:,:), 1NMRT(:),nmrb(:),spinb(:,:),qq(:,:) integer*4 iffs common/cname/asname(IOX,MX,3,3),polname(IOX,3,3),po(3,3), 1psname(IOX,MX,3,3),aaname(IOX,3*MX,3,3,3),gname(IOX,3*MX,3,3), 2alphaname(IOX,3*MX,3,3),rname(IOX,3,MX),qsname(IOX,MX,3,3,3), 3gpo(3,3),apo(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3), 2gradname(IOX,3*MX),natname(IOX),NAME c WRITE(6,3000) 3000 FORMAT(/, 1 ' Transfer of Molecular Tensors in Cartesian Coordinates' 8,/,/, ' By Petr Bour, Prague 1994-2009',/,/, 9 ' INPUT FILES:',/, 9 ' Target: BIG.X coordinates',/, 1 ' BIG. FC force field',/, 1 ' BIG.TEN APT and AAT tensors',/, 3 ' BIG.TTT ROA tensors',/, 3 ' Source: SMALL.X coordinates',/, 2 ' SMALL. FC force field',/, 2 ' SMALL.TEN APTs and AATs',/, 2 ' SMALL.QEN qudrupole derivativess',/, 4 ' SMALL.TTT ROA tensors',/, 4 ' (or named .x,.fc, .ten and .ttt)',/, 4 ' Options: CTT.INP atom assignment',/,/, 5 ' OUTPUT:',/, 5 ' FILE.X same as BIG.X',/, 7 ' FILE. FC improved BIG.FC',/, 7 ' FILE.TEN improved BIG.TEN',/, 7 ' FILE.TTT improved BIG.TTT',/, 8 ' CTT.OUT control list',/) OPEN(3,FILE='CCT.OUT') WRITE(3,3000) call inipol(po,gpo,apo) NATS_max=0 bohr=0.529177D0 iffs=-1 lpol=.false. lpolg=.false. lpola=.false. c c Input geometry of the target molecule: inquire(FILE='BIG.X',EXIST=lxfile) if(lxfile)then OPEN(2,FILE='BIG.X',STATUS='OLD') CALL READRX(MX,NAT,R,ND,IBONDSB,NBOB) else call report('BIG.X not found') endif CLOSE(2) WRITE(3,*)NAT, ' atoms in the big molecule' c c Open file with option and overlaps and read it in: OPEN(2,FILE='CCT.INP',STATUS='OLD') CALL readopt(NS,IB,IS,IND,MX,MXB,INDBIG,IOX,NO,NAT, 1LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT,LAPT, 2LALF,LRDO,LNAMES,NAME,LPOLYMER,LHALTONERROR,CUTOFF,n1,n2,okind, 3LXX3,LXX4,ALIM,L36,L37,lq,LDATOM,IOPOL,lsos,lmut,lmpt,IFC,EFC, 4iweights,lgrad,sgrad,fgrad,alphf,lnmr,lspin,LQUAD,SOSNM) CLOSE(2) c c if tensors constructed from named fragments, load them all: if(LNAMES)then NATS_max=0 DO 11 I=1,NO call getname(NAME(I),basicname,NE) write(6,*) write(6,*)' Fragment ',basicname(1:NE) xname=basicname(1:NE)//'.x' write(6,*)xname(1:NE+2) inquire(file=xname,exist=lxfile) if(.not.lxfile)call report('Coordinates not found') open(2,file=xname) CALL READRXS(MX,NATS,RS,NDS,NBO) CLOSE(2) NATNAME(I)=NATS do 13 ia=1,NATS NDNAME(I,ia)=NDS(ia) do 13 ix=1,3 13 rname(I,ix,ia)=RS(ix,ia) if(NATNAME(I).gt.NATS_max)NATS_max=NATNAME(I) 11 WRITE(3,*)NATS, ' atoms in the fragment' else OPEN(2,FILE='SMALL.X',STATUS='OLD') CALL READRXS(MX,NATS,RS,NDS,NBO) CLOSE(2) WRITE(3,*)NATS, ' atoms in the small molecule' WRITE(6,*)NATS, ' atoms in the small molecule' NATS_max=NATS NATNAME(1)=NATS do 131 ia=1,NATS NDNAME(1,ia)=NDS(ia) do 131 ix=1,3 131 rname(1,ix,ia)=RS(ix,ia) endif c names c c check consistency in atomic types: do 22 i=1,NO if(LNAMES)then NATS=natname(i) do 23 ia=1,NATS 23 NDS(ia)=ndname(i,ia) endif DO 22 IAB=1,NAT IF(INDBIG(I,IAB).NE.0)THEN IF(ND(IAB).NE.NDS(ABS(INDBIG(I,IAB))))THEN WRITE(3,3300)I,IAB,ND(IAB),NDS(ABS(INDBIG(I,IAB))) 3300 FORMAT('Overlap ',I2,', atom ',I3,': ',2I4) CALL REPORT(' Different atomic types !') ENDIF IF(IABS(INDBIG(I,IAB)).gt.NATS)then write(3,*)'overlap ',i write(3,*)' NATS ',NATS write(3,*)'IND ',IAB write(3,*)'INDB ',INDBIG(I,IAB) write(6,*)'overlap ',i write(6,*)' NATS ',NATS write(6,*)'IND ',IAB write(6,*)'INDB ',INDBIG(I,IAB) call report('Overflow NATS') ENDIF ENDIF 22 CONTINUE c if(lgrad)then allocate(GG(3*NAT)) INQUIRE(FILE='BIG.GR',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3007) WRITE(6,3007) 3007 FORMAT('File BIG.GR not found, constants set to zero') DO 20 I=1,3*NAT 20 GG(I)=0.0d0 ELSE OPEN(2,FILE='BIG.GR') DO 19 I=1,3*NAT 19 read(2,*)(GG(ix+3*(I-1)),ix=1,3) CLOSE(2) WRITE(3,3008)NAT WRITE(6,3008)NAT 3008 FORMAT(' Gradient for the big molecule found,',/,I10,' atoms') ENDIF ENDIF if(lnmr)then allocate(nmrb(NAT)) INQUIRE(FILE='BIG.NMR',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3012) WRITE(6,3012) 3012 FORMAT('File BIG.NMR not found, constants set to zero') DO 25 I=1,NAT 25 nmrb(I)=0.0d0 ELSE call rdnmr(nmrb,'BIG.NMR') ENDIF endif if(lspin)then allocate(spinb(NAT,NAT)) INQUIRE(FILE='BIG.SPI',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3013) WRITE(6,3013) 3013 FORMAT('File BIG.SPI not found, constants set to zero') DO 26 I=1,NAT DO 26 J=1,NAT 26 spinb(I,J)=0.0d0 ELSE call rdspin(spinb,NAT,'BIG.SPI') ENDIF endif if(LDIA.or.LOFF)then allocate(FF(3*NAT,3*NAT)) INQUIRE(FILE='BIG.FC',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3003) WRITE(6,3003) 3003 FORMAT('File BIG.FC not found, constants set to zero') DO 1 I=1,3*NAT DO 1 J=1,3*NAT 1 FF(I,J)=0.0d0 ELSE OPEN(2,FILE='BIG.FC',STATUS='OLD') CALL READFF(NAT,FF) CLOSE(2) WRITE(3,3001)NAT WRITE(6,3001)NAT 3001 FORMAT(' Force field of the big molecule found,',/,I10,' atoms') ENDIF ENDIF c c Special section for named fragments: c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn if(lnames)then allocate(FFS(NO*(3*NATS_max)**2)) write(6,*)' allocated named FFS' if(lnmr)allocate(NMRNS(NO*NATS_max)) if(lspin)allocate(SPINNS(NO*NATS_max**2)) do 14 i=1,NO NATS=natname(i) call getname(NAME(i),basicname,NE) write(6,*) write(6,*)' Fragment ',basicname(1:NE) c c Re-create RS for ROA and VCD tensor transformations: do 1311 ia=1,NATS do 1311 ix=1,3 1311 RS(ix,ia)=rname(i,ix,ia) c IF(lgrad)THEN fname=basicname(1:NE)//'.gr' write(6,*)fname(1:NE+3) INQUIRE(FILE=fname,EXIST=OPT) IF(.NOT.OPT)THEN CALL REPORT(basicname(1:NE+3)//' not found.') ELSE call rgr(fname(1:NE+3),NATS,IOX,MX,i,gradname) WRITE(3,3009)NATS WRITE(6,3009)NATS 3009 FORMAT(' Gradient found,',I10,' atoms') ENDIF ENDIF IF(lnmr)then fname=basicname(1:NE)//'.nmr' write(6,*)fname(1:NE+4) INQUIRE(FILE=fname,EXIST=OPT) IF(.NOT.OPT)THEN CALL REPORT(basicname(1:NE)//'.nmr not found.') ELSE allocate(NMRT(NATS)) call rdnmr(NMRT,fname) do 24 ii=1,NATS 24 NMRNS((i-1)*NATS_max+ii)=NMRT(ii) deallocate(NMRT) ENDIF ENDIF IF(lspin)then fname=basicname(1:NE)//'.spi' write(6,*)fname(1:NE+4) INQUIRE(FILE=fname,EXIST=OPT) IF(.NOT.OPT)THEN CALL REPORT(basicname(1:NE)//'.spi not found.') ELSE allocate(SPINT(NATS,NATS)) call rdspin(SPINT,NATS,fname) do 29 ii=1,NATS do 29 jj=1,NATS 29 SPINNS((i-1)*NATS_max**2+(ii-1)*NATS_max+jj)=SPINT(ii,jj) deallocate(SPINT) ENDIF ENDIF IF(LDIA.OR.LOFF)THEN fname=basicname(1:NE)//'.fc' write(6,*)fname(1:NE+3) INQUIRE(FILE=fname,EXIST=OPT) IF(.NOT.OPT)THEN CALL REPORT(basicname(1:NE)//'.fc not found.') ELSE allocate(FT(3*NATS,3*NATS)) OPEN(2,FILE=fname,STATUS='OLD') CALL READFF(NATS,FT) CLOSE(2) do 21 ii=1,3*NATS do 21 jj=1,3*NATS 21 FFS((i-1)*(3*NATS_max)**2+(ii-1)*3*NATS_max+jj)=FT(ii,jj) deallocate(FT) WRITE(3,3002)NATS WRITE(6,3002)NATS 3002 FORMAT(' Force field found,',/,I10,' atoms') ENDIF ENDIF c IF(LABS.OR.LVCD)THEN nname=basicname(1:NE)//'.ten' write(6,*)nname(1:NE+4) OPEN(15,FILE=nname,STATUS='OLD') CALL READTEN(MX,PS,AS,LAPT,RS,NATS) CLOSE(15) do 18 ia=1,NATS do 18 ix=1,3 do 18 jx=1,3 asname(i,ia,ix,jx)=AS(ia,ix,jx) 18 psname(i,ia,ix,jx)=PS(ia,ix,jx) IF(LQUAD)THEN nname=basicname(1:NE)//'.qen' write(6,*)nname(1:NE+4) allocate(qq(3*NATS,6)) call reqd(nname,qq,NATS,PS,MX,RS) call asnqn(qsname,IOX,MX,i,qq,NATS) deallocate(qq) ENDIF ENDIF C IF(LROA.or.LRAM)THEN tname=basicname(1:NE)//'.ttt' write(6,*)tname(1:NE+4) OPEN(15,FILE=tname,STATUS='OLD') CALL READTTT(MX,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO' CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,1,RS) ENDIF c write the local tensors: if(LWR) 1 CALL WRITETTT(MX,NATS,ALPHAS,AAS,GS,tname(1:NE)//'.l.ttt') do 28 iax=1,3*NATS do 28 ix=1,3 do 28 jx=1,3 alphaname(i,iax,ix,jx)=ALPHAS(iax,ix,jx) gname(i,iax,ix,jx)=GS(iax,ix,jx) do 28 kx=1,3 28 aaname(i,iax,ix,jx,kx)=AAS(iax,ix,jx,kx) C Now small tensors are in atomic origins c c read static polarizability: pname=basicname(1:NE)//'.pol' inquire(file=pname,exist=lpol) if(lpol)then call w36(pname//' found') call readpol(po,pname) endif c c read G tensor pname=basicname(1:NE)//'.g.pol' inquire(file=pname,exist=lpolg) if(lpolg)then call w36(pname//' found') call readpolg(gpo,pname) endif c c read A tensor pname=basicname(1:NE)//'.a.pol' inquire(file=pname,exist=lpola) if(lpola)then call w36(pname//' found') call readpola(apo,pname) endif c transform G, A to small molecule center if(lpola.or.lpolg)then do 301 ix=1,3 coo(ix)=0.0d0 zsum=0.0d0 do 3011 ia=1,NATS if(lq)then sf=dble(NDNAME(i,ia)) else sf=1.0d0 endif zsum=zsum+sf 3011 coo(ix)=coo(ix)+rname(i,ix,ia)*sf 301 coo(ix)=coo(ix)/bohr/zsum call TTT1(po,apo,gpo,1,coo,1) endif c c transcript alpha,G,A to named strings: do 302 ix=1,3 do 302 jx=1,3 gpolname(I,ix,jx)=gpo(ix,jx) polname(I,ix,jx)=po(ix,jx) do 302 kx=1,3 302 apolname(I,ix,jx,kx)=apo(ix,jx,kx) c ENDIF c LROA LRAM c c dynamic polarizabilities if(lsos)then open(44,file=basicname(1:NE)//'.sos') read(44,*)np,wmin,wmax if(I.eq.1)then allocate(alsos(NO*np*nel)) write(6,*)' allocated named alsos' endif c read polarizabilities: call readsos(alsos,np,I,nel) c for each frequency, rewrite to smaller arrays and transform: c get the center: do 1301 ix=1,3 coo(ix)=0.0d0 zsum=0.0d0 do 13011 ia=1,NATS if(lq)then sf=dble(NDNAME(i,ia)) else sf=1.0d0 endif zsum=zsum+sf 13011 coo(ix)=coo(ix)+rname(i,ix,ia)*sf 1301 coo(ix)=coo(ix)/bohr/zsum call trdpol(np,alsos,coo,I,nel) if(I.eq.1)then if(iweights.ne.0)then if(iweights.lt.0)then c atomic weights for all frequencies allocate(weights(NATS_max*NO*np)) write(6,*)' allocated weights - all freqs' else c atomic weights independent on frequencies allocate(weights(NATS_max*NO)) write(6,*)' allocated weights - constant' endif endif endif if(iweights.ne.0)then call rweights(NO,i,iweights,weights,NATS_max,NATS, 1 basicname(1:NE)//'.wei',np) write(6,*)'atomic weights read from ', 1 basicname(1:NE)//'.wei' endif endif C 14 continue c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn else if(lgrad)then INQUIRE(FILE='SMALL.GR',EXIST=OPT) IF(.NOT.OPT)CALL REPORT('SMALL.GR with gradient not found') call rgr('SMALL.GR',NATS,IOX,MX,1,gradname) endif IF(lnmr)then INQUIRE(FILE='SMALL.NMR',EXIST=OPT) IF(.NOT.OPT)THEN CALL REPORT('SMALL.NMR not found.') ELSE allocate(NMRNS(NATS)) call rdnmr(NMRNS,'SMALL.NMR') ENDIF ENDIF IF(lspin)then INQUIRE(FILE='SMALL.SPI',EXIST=OPT) IF(.NOT.OPT)THEN CALL REPORT('SMALL.SPI not found.') ELSE allocate(SPINT(NATS,NATS)) call rdspin(SPINT,NATS,'SMALL.SPI') do 27 ii=1,NATS do 27 jj=1,NATS 27 SPINNS((ii-1)*NATS+jj)=SPINT(ii,jj) deallocate(SPINT) ENDIF ENDIF if(LOFF.or.LDIA)then INQUIRE(FILE='SMALL.FC',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,*)'Force field of small not found' WRITE(*,*)'Force field of small not found' IF(LDIA.OR.LOFF)CALL REPORT('Cannot do force field') ELSE allocate(FFS(1:(3*NATS)**2)) allocate(FT(3*NATS,3*NATS)) write(6,*)' allocated FFS' OPEN(2,FILE='SMALL.FC',STATUS='OLD') CALL READFF(NATS,FT) CLOSE(2) WRITE(3,3002)NATS WRITE(6,3002)NATS do 30 ii=1,3*NATS do 30 jj=1,3*NATS 30 FFS((ii-1)*3*NATS+jj)=FT(ii,jj) deallocate(FT) ENDIF endif IF(LABS.OR.LVCD)THEN OPEN(15,FILE='SMALL.TEN',STATUS='OLD') CALL READTEN(MX,PS,AS,LAPT,RS,NATS) CLOSE(15) IF(LQUAD)THEN allocate(qq(3*NATS,6)) call reqd('SMALL.QEN',qq,NATS,PS,MX,RS) call asnq(QS,qq,MX,NATS) deallocate(qq) ENDIF ENDIF C C IF(LROA.OR.LRAM)THEN OPEN(15,FILE='SMALL.TTT',STATUS='OLD') CALL READTTT(MX,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO for small' CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,1,RS) ENDIF C Now small tensors are in atomic origins C IF(LRDO)THEN WRITE(3,*)' Alpha from SMALL_ALPHA.TTT for small' OPEN(15,FILE='SMALL_ALPHA.TTT',STATUS='OLD') CALL READTTT(MX,ALPHAS,AAT,GT,NATS) CLOSE(15) ENDIF c c write the local tensors: if(LWR) 1 CALL WRITETTT(MX,NATS,ALPHAS,AAS,GS,'SMALL.L.TTT') c c read static polarizability: inquire(file='SMALL.POL',exist=lpol) if(lpol)then call w36('SMALL.POL found') call readpol(po,'SMALL.POL') endif c c read G inquire(file='SMALL.G.POL',exist=lpolg) if(lpolg)then call w36('SMALL.G.POL found') call readpolg(gpo,'SMALL.G.POL') endif c c read A inquire(file='SMALL.A.POL',exist=lpola) if(lpola)then call w36('SMALL.A.POL found') call readpola(apo,'SMALL.A.POL') endif c transform G, A to small molecule center if(lpola.or.lpolg)then do 313 ix=1,3 coo(ix)=0.0d0 zsum=0.0d0 do 3131 ia=1,NATS if(lq)then sf=dble(NDS(ia)) else sf=1.0d0 endif zsum=zsum+sf 3131 coo(ix)=coo(ix)+RS(ix,ia)*sf 313 coo(ix)=coo(ix)/bohr/zsum call TTT1(po,apo,gpo,1,coo,1) endif c c record as first item to the named strings: do 312 ix=1,3 do 312 jx=1,3 polname(1,ix,jx)=po(ix,jx) gpolname(1,ix,jx)=gpo(ix,jx) do 312 kx=1,3 312 apolname(1,ix,jx,kx)=apo(ix,jx,kx) endif c lroa/ram c dynamic polarizabilities if(lsos)then open(44,file='SMALL.SOS') read(44,*)np,wmin,wmax allocate(alsos(np*nel)) write(6,*)' allocated alsos' c read polarizabilities: call readsos(alsos,np,1,nel) c for each frequency, rewrite to smaller arrays and transform: c get the center: do 2301 ix=1,3 coo(ix)=0.0d0 zsum=0.0d0 do 23011 ia=1,NATS if(lq)then sf=dble(NDS(ia)) else sf=1.0d0 endif zsum=zsum+sf 23011 coo(ix)=coo(ix)+RS(ix,ia)*sf 2301 coo(ix)=coo(ix)/bohr/zsum call trdpol(np,alsos,coo,1,nel) if(iweights.ne.0)then if(iweights.lt.0)then allocate(weights(NATS*np)) else allocate(weights(NATS)) endif call rweights(1,1,iweights,weights,NATS,NATS,'WEIGHTS.TXT',np) write(6,*)'WEIGHTS.TXT read in' endif endif ENDIF c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn c lnames c IF(LABS.OR.LVCD)THEN INQUIRE(FILE='BIG.TEN',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3004) WRITE(6,3004) 3004 FORMAT('File BIG.TEN not found, constants set to zero') DO 2 I=1,NAT DO 2 J=1,3 DO 2 K=1,3 PB(I,J,K)=0.0d0 2 AB(I,J,K)=0.0d0 ELSE OPEN(15,FILE='BIG.TEN') CALL READTEN(MX,PB,AB,LAPT,R,NAT) write(6,*)'Tensors BIG.TEN found' write(3,*)'Tensors BIG.TEN found' CLOSE(15) ENDIF IF(LQUAD)THEN INQUIRE(FILE='BIG.QEN',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3014) WRITE(6,3014) 3014 FORMAT('File BIG.QEN not found, constants set to zero') DO 34 I=1,NAT DO 34 J=1,3 DO 34 K=1,3 DO 34 L=1,3 34 QB(I,J,K,L)=0.0d0 ELSE allocate(qq(3*NAT,6)) call reqd('BIG.QEN',qq,NAT,PB,MX,R) call asnq(QB,qq,MX,NAT) deallocate(qq) write(6,*)'Tensors BIG.QEN found' write(3,*)'Tensors BIG.QEN found' ENDIF ENDIF ENDIF C INQUIRE(FILE='BIG.TTT',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3005) WRITE(6,3005) 3005 FORMAT('File BIG.TTT not found, tensors set to zero') DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX DO 3 I=1,3 DO 3 J=1,3 ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 DO 3 K=1,3 3 AA(IIND,I,J,K)=0.0d0 ELSE OPEN(15,FILE='BIG.TTT') CALL READTTT(MX,ALPHA,AA,G,NAT) CLOSE(15) ENDIF IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO for big' CALL TTTT(3*MX,ALPHA,AA,G,NAT,3,R) ELSE CALL TTTT(3*MX,ALPHA,AA,G,NAT,1,R) ENDIF c write the local tensors: if(LWR) 1 CALL WRITETTT(MX,NAT,ALPHA,AA,G,'BIG.L.TTT') C if(lsos)then WRITE(6,*)' Transferring vibrational tensors omitted' else WRITE(6,*)' Transferring the tensors ...' CALL IMPROVE(FF,FFS,R,RS,NAT,NATS,NS,IND,IB,IS, 1 AM,INDBIG,NO,IBONDSB,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA, 2 IWG,LROA,LRAM,ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT, 3 LNAMES,LHALTONERROR,CUTOFF,n1,n2,okind,NATS_max,lpol, 1 lpolg,lpola,lq,ND,LDATOM,IOPOL,GG,lgrad,sgrad,fgrad,alphf, 1 lnmr,lspin,NMRNS,SPINNS,nmrb,spinb,LQUAD,QS,QB) endif if(iffs.ne.-1)deallocate(FFS) write(6,*) if(lnmr)call wrnmr(NAT,nmrb,ND) if(lspin)call wspin(spinb,NAT) c if(lgrad)call svgr(GG,NAT) if(LDIA.or.LOFF)then if(IFC.eq.0)then OPEN(2,FILE='FILE.FC') CALL WRITEFF(NAT,FF) CLOSE(2) call w36(' FILE.FC written ') else CALL WRITEAFF(MX,NAT,FF,EFC,CUTOFF,R,LWR,j33,a33) call w36(' AFILE.FC written ') endif endif IF(LABS.OR.LVCD)then CALL WRITETEN(MX,PB,AB,R,NAT) if(LQUAD)call wrqq(MX,PB,QB,R,NAT) endif IF(LROA.or.LRAM)THEN c write the local tensors: if(LWR) 1 CALL WRITETTT(MX,NAT,ALPHA,AA,G,'FILE.L.TTT') CALL TTTT(3*MX,ALPHA,AA,G,NAT,2,R) if(.not.lroa)then do 31 i=1,3*NAT do 31 ix=1,3 do 31 jx=1,3 G(i,ix,jx)=0.0d0 do 31 kx=1,3 31 AA(i,ix,jx,kx)=0.0d0 write(6,*)' ROA tensors set to zero' write(3,*)' ROA tensors set to zero' endif c write the common origin tensors: CALL WRITETTT(MX,NAT,ALPHA,AA,G,'FILE.TTT') ENDIF OPEN(2,FILE='FILE.X') CALL WRITERX(MX,NAT,R,ND,IBONDSB) CLOSE(2) WRITE(3,*) WRITE(6,*) WRITE(3,*)' FILE.X written ' WRITE(6,*)' FILE.X written ' c if(lsos)then allocate(alsosb(np*nel)) write(6,7666)lmpt,lmut,iweights 7666 format(' Frequency-dependent polarizability transfer',/,/, 2 ' lmpt: ',l3,/, 3 ' lmut: ',l3,/, 1 ' iweights: ',i3,/) if(lmpt)then call IMPROVEMPT(np,INDBIG,IOX,MX,IO,LNAMES,NATNAME, 1 RS,R,IND,rname,NATS,LINV,NO,NAT,lq,ND,alsosb,alsos, 1 lmut,wmin,wmax,nel) else if(iweights.ne.0)then call IMPROVESOSw(np,INDBIG,IOX,MX,IO,LNAMES,NATNAME, 1 RS,R,IND,rname,NATS,LINV,NO,NAT,alsosb,alsos, 1 lmut,wmin,wmax,weights,nel,NATS_max,iweights,CUTOFF,SOSNM) else call IMPROVESOS(np,INDBIG,IOX,MX,IO,LNAMES,NATNAME, 1 RS,R,IND,rname,NATS,LINV,NO,NAT,lq,ND,alsosb,alsos, 1 lmut,wmin,wmax,nel,CUTOFF,SOSNM) endif endif call writesos(alsosb,np,wmin,wmax,nel,SOSNM) endif if(LXX3)then if(NAT.gt.n30)call report('too many atoms for LXX3') if(.not.lnames)then INQUIRE(FILE='SMALL.33',EXIST=OPT) IF(.NOT.OPT)CALL REPORT('SMALL.33 not found') allocate(cs((3*NATS)**3)) CALL READFC('SMALL.33',NATS,cs,ncs) write(6,*)ncs,' constants in SMALL.33' write(3,*)ncs,' constants in SMALL.33' else allocate(cs(NO*81*NATS_max**3)) do 17 I=1,NO NE=80 DO 1500 NE=80,1,-1 1500 if(NAME(I)(NE:NE).ne.' ')goto 1600 1600 NATS=NATNAME(I) allocate(ct((3*NATS)**3)) CALL READFC(NAME(I)(1:NE)//'.33',NATS,ct,ncs) write(6,6017)I,ncs 6017 format(' Overlap ',i3,',',I15,' Cijk constants') do 171 j=1,(NATS*3)**3 171 cs((I-1)*81*NATS_max**3+j)=ct(j) 17 deallocate(ct) endif INQUIRE(FILE='BIG.33',EXIST=OPT) ncs=0 allocate(cb((3*NAT)**3)) IF(OPT)then CALL READFC('BIG.33',NAT,cb,ncs) else CALL vz0((3*NAT)**3,cb) endif write(6,*)ncs,' constants in BIG.33' write(3,*)ncs,' constants in BIG.33' CALL IMPROVE34(3,cb,cs,R,RS,NAT,NATS,IND,INDBIG,NO, 1 IBONDSB,NBOB,LWR,IWG,LSTRICT,LINV,LNAMES,LHALTONERROR,CUTOFF, 2 okind,NATS_max) open(20,file='SCR.33') call WRITEFC(20,NAT,cb,ALIM,ncs) write(6,*)'SCR.33 written' deallocate(cs) deallocate(cb) endif c if(LXX4)then if(NAT.gt.n30)call report('too many atoms for LXX4') if(.not.lnames)then INQUIRE(FILE='SMALL.44',EXIST=OPT) IF(.NOT.OPT)CALL REPORT('SMALL.44 not found') allocate(cs((3*NATS)**3)) CALL READFC('SMALL.44',NATS,cs,ncs) write(6,*)ncs,' constants in SMALL.44' write(3,*)ncs,' constants in SMALL.44' else allocate(cs(NO*81*NATS_max**3)) do 172 I=1,NO NE=80 DO 1501 NE=80,1,-1 1501 if(NAME(I)(NE:NE).ne.' ')goto 1601 1601 NATS=NATNAME(I) allocate(ct((3*NATS)**3)) CALL READFC(NAME(I)(1:NE)//'.44',NATS,ct,ncs) write(6,6027)I,ncs 6027 format(' Overlap ',i3,',',I15,' Diijk constants') do 173 j=1,(NATS*3)**3 173 cs((I-1)*81*NATS_max**3+j)=ct(j) 172 deallocate(ct) endif INQUIRE(FILE='BIG.44',EXIST=OPT) ncs=0 allocate(cb((3*NAT)**3)) IF(OPT)then CALL vz0((3*NAT)**3,cb) CALL READFC('BIG.44',NAT,cb,ncs) endif write(6,*)ncs,' constants in BIG.44' write(3,*)ncs,' constants in BIG.44' CALL IMPROVE34(4,cb,cs,R,RS,NAT,NATS,IND,INDBIG,NO, 1 IBONDSB,NBOB,LWR,IWG,LSTRICT,LINV,LNAMES,LHALTONERROR,CUTOFF, 2 okind,NATS_max) open(20,file='SCR.44') call WRITEFC(20,NAT,cb,ALIM,ncs) write(6,*)'SCR.44 written' deallocate(cs) deallocate(cb) endif c c transfer of two-atomic quartic constants: do 36 ii=1,2 if((ii.eq.1.and.L36).or.(ii.eq.2.and.L37))then if(ii.eq.1)then sss='SMALL.36' ssb='BIG.36' ssc='SCR.36' else sss='SMALL.37' ssb='BIG.37' ssc='SCR.37' endif if(NAT.gt.n30)call report('too many atoms for 36/7') if(.not.lnames)then INQUIRE(FILE=sss,EXIST=OPT) IF(.NOT.OPT)CALL REPORT(sss//' not found') allocate(cs(81*NATS**3)) open(20,file=sss) CALL READ36(20,NATS,cs,ncs) close(20) write(6,*)ncs,' constants in '//sss write(3,*)ncs,' constants in '//sss else allocate(cs(NO*81*NATS_max**3)) do 174 I=1,NO NE=80 DO 1502 NE=80,1,-1 1502 if(NAME(I)(NE:NE).ne.' ')goto 1602 1602 NATS=NATNAME(I) allocate(ct(81*NATS**3)) if(ii.eq.1)then write(6,*)NAME(I)(1:NE)//'.36' OPEN(20,FILE=NAME(I)(1:NE)//'.36') else OPEN(20,FILE=NAME(I)(1:NE)//'.37') write(6,*)NAME(I)(1:NE)//'.37' endif CALL READ36(20,NATS,ct,ncs) close(20) write(6,6037)I,ncs write(3,6037)I,ncs 6037 format(' Overlap ',i3,',',I15,' 37iijk constants') do 175 j=1,81*NATS**3 175 cs((I-1)*81*NATS_max**3+j)=ct(j) 174 deallocate(ct) endif ncs=0 INQUIRE(FILE=ssb,EXIST=OPT) allocate(cb(81*NAT**3)) IF(OPT)then open(20,file=ssb) CALL READ36(20,NAT,cb,ncs) close(20) else do 176 i=1,81*NAT**3 176 cb(i)=0.0d0 endif write(6,*)ncs,' constants in '//ssb write(3,*)ncs,' constants in '//ssb if(ii.eq.2)then CALL IMPROVE34(5,cb,cs,R,RS,NAT,NATS,IND,INDBIG,NO, 1 IBONDSB,NBOB,LWR,IWG,LSTRICT,LINV,LNAMES,LHALTONERROR,CUTOFF, 2 okind,NATS_max) else CALL IMPROVE36(cb,cs,R,RS,NAT,NATS,IND, 1 INDBIG,NO,IBONDSB,NBOB,LWR,IWG,LSTRICT,LINV,LNAMES, 3 LHALTONERROR,CUTOFF,okind,NATS_max) endif open(20,file=ssc) if(ii.eq.1)then call WRITE36(20,NAT,cb,ALIM,ncs,.true.) else call WRITE36(20,NAT,cb,ALIM,ncs,.false.) endif close(20) deallocate(cs) deallocate(cb) write(6,*)ssc//' written,',ncs,' constants >',ALIM write(3,*)ssc//' written,',ncs,' constants >',ALIM endif 36 enddo c CLOSE(3) WRITE(6,*)' Program finished OK.' STOP END SUBROUTINE IMPROVE36(db,ds,RB,RS,NAT,NATS,IND, 1INDBIG,NO,IBONDS,NBOB,LWR,IWG,LSTRICT,LINV,LNAMES, 3LHALTONERROR,CUTOFF,okind,NATS_max) implicit none INTEGER*4 MX,MXB,IOX PARAMETER (MX=4000,MXB=140,IOX=1500) INTEGER*4 okind,IND(MX),INDBIG(IOX,MX),LLX,IBONDS(4,MX), 1NBOB(MX),IOL(IOX),ISUM,I,IA,II,IIOII,IIX,IO,i1,i2,lx, 2IOIJ,IP,IPASS,IPO,IWG,IX,J,JA,IERC,JX,NAT,NATIO,NATS,NE,NF,NO, 3JJX,ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,MX),nats_o(IOX), 5KKX,KX,is,IAS,JAS,IPOOLD,NATS_max,ido real*8 RB(3,MX),RS(3,NATS),TIJ(3),OTIJ(3),WF(IOX),db(*),ds(*), 1DISTM,SUM,FNEW,DIST,FOLD,CTF,CUTOFF,RM(IOX,3,3),xi,yi,zi,xj,yj, 2zj,SIA,SJB,SKC,SLD LOGICAL LWR,LSTRICT,LINV,LNAMES,LHALTONERROR real*8 A,XS,XB,TOL,RMS integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS real*8 aname,pname,aaname,gname,alphaname,rname,pol,polname, 1gpo,apo,apolname,gpolname,gradname,qsname integer*4 natname character*80 NAME(IOX) common/cname/aname(IOX,MX,3,3),polname(IOX,3,3),pol(3,3), 1pname(IOX,MX,3,3),aaname(IOX,3*MX,3,3,3),gname(IOX,3*MX,3,3), 2alphaname(IOX,3*MX,3,3),rname(IOX,3,MX),qsname(IOX,MX,3,3,3), 3gpo(3,3),apo(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3), 2gradname(IOX,3*MX),natname(IOX),NAME c IF(NO.EQ.0)call report('NO=0') WRITE(3,*)' Two-atomic quartic constants requested' IF(LSTRICT)WRITE(3,*)' Only transformed atoms in the beds' IF(IWG.EQ.1)WRITE(3,*)' Multiple overlaps will be weighted' IF(LINV)WRITE(3,*)' Smaller fragment as optical antipode' IF(LNAMES)WRITE(3,*)' Named small fragments' IF(LHALTONERROR)WRITE(3,*)' Will stop on error' if(okind.eq.0)write(3,*)'selection on centers of masses' if(okind.eq.1)write(3,*)'selection on RMS distances' if(okind.eq.2)write(3,*)'selection on molecular RMS distances' WRITE(3,3999)CUTOFF 3999 format('Distance cutoff ',f12.2,' A') ctf=CUTOFF**2 IOIJ=0 IPOOLD=0 C C Look at all atom diads i-j: DO 30 IA=1,NAT write(6,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) xi=RB(1,IA) yi=RB(2,IA) zi=RB(3,IA) DO 30 JA=IA,NAT xj=RB(1,JA) yj=RB(2,JA) zj=RB(3,JA) if((xi-xj)**2+(yi-yj)**2+(zi-zj)**2.lt.ctf)then if(LWR)write(3,3008)IA,JA 3008 format(/,' Atoms ',2I6) C C Look at all overlaps io, find which fits best ijk c c for okind=2, look at the overlaps only once c cccccccccccccccccccccccccccccccccccccccccccccccccccc IF(okind.ne.2.or.(IA.eq.1.and.JA.EQ.1))THEN C C Calculate center of i/j radiusvectors: DO 22 IX=1,3 22 TIJ(IX)=(RB(IX,IA)+RB(IX,JA))/3.0d0 DISTM=1.0d5 NF=0 DO 21 IO=1,NO C C Check, if io contains ij: IF(INDBIG(IO,IA).GT.0.AND.INDBIG(IO,JA).GT.0)THEN NF=NF+1 IOL(NF)=IO C C Calculate center of the overlap and distance from i-j center c i.........TIJ.............j c ----------OTIJ--------------------- DO 23 IX=1,3 23 OTIJ(IX)=0.0D0 NATIO=0 DO 24 II=1,NAT IIOII=INDBIG(IO,II) C If Lstrict, take only transformed atoms IF(LSTRICT.AND.IIOII.LT.0)IIOII=0 IIOII=ABS(IIOII) IF(IIOII.NE.0)THEN NATIO=NATIO+1 DO 25 IX=1,3 25 OTIJ(IX)=OTIJ(IX)+RB(IX,II) ENDIF 24 CONTINUE c c If required, calculate RMS overlap and use this criterium instead c of the distances,remember rotation matrix (RM) if(okind.gt.0)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IO, 1 IA,JA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) do 55 IX=1,3 do 55 JX=1,3 55 RM(NF,IX,JX)=A(IX,JX) ipass_o(NF)=IPASS WF(NF)=DSQRT(RMS) ierc_o(NF)=IERC itu_o(NF)=ITU do 58 I=1,ITU 58 ind_o(NF,I)=IND(I) do 59 I=1,3 do 59 J=1,3 59 RM(NF,I,J)=A(I,J) nats_o(NF)=NATS IF(RMS.LT.DISTM)THEN DISTM=RMS IOIJ=IO ENDIF else c Use distances: DO 26 IX=1,3 26 OTIJ(IX)=OTIJ(IX)/DBLE(NATIO) DIST=(OTIJ(1)-TIJ(1))**2+(OTIJ(2)-TIJ(2))**2+(OTIJ(3)-TIJ(3))**2 IF(DIST.LT.DISTM)THEN DISTM=DIST IOIJ=IO ENDIF WF(NF)=DSQRT(DIST) endif ENDIF 21 CONTINUE c WRITE(3,3003)NF,IOIJ 3003 FORMAT(I5,' overlaps found; best is # ',I4) C C Calculate the weighting factors: SUM=0.0d0 DO 43 I=1,NF 43 IF(WF(I).GT.1.0d-9)SUM=SUM+1.0d0/WF(I) DO 42 I=1,NF IF(WF(I).LE.1.0d-9)THEN WF(1)=1.0d0 NF=1 GOTO 44 ENDIF WF(I)=1.0d0/WF(I)/SUM 42 CONTINUE 44 CONTINUE C c cccccccccccccccccccccccccccccccccccccccccccccccccc ENDIF c C Loop over possible overlaps: ISUM=0 DO 45 IP=1,NF IPO=IOL(IP) IF(IWG.EQ.0)THEN IF(IPO.NE.IOIJ)GOTO 45 WF(IP)=1.0d0 ENDIF ISUM=ISUM+1 C C Find the atomic groups to be placed on together for IA-JA: c and the rotation matrix if(okind.eq.0)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1 IA,JA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) else IPASS=ipass_o(IP) IERC=ierc_o(IP) ITU=itu_o(IP) do 56 I=1,ITU 56 IND(I)=ind_o(IP,I) do 57 I=1,3 do 57 J=1,3 57 A(I,J)=RM(IP,I,J) NATS=nats_o(IP) endif IF(LWR)WRITE(3,*)IPASS,' passes' if(IERC.ne.0)goto 20 WRITE(3,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU) WRITE(3,3005)(INDBIG(IPO,IND(I)),I=1,ITU) 3004 FORMAT('Overlap ',I2,': weighting factor: ',F7.4,';', 1 I4,' atom bed:',/, 2 ' Big : ',20I3) 3005 FORMAT( ' Small: ',20I3,/) if(LNAMES)write(3,3399)NAME(IPO) 3399 FORMAT(A80) C IF(LINV)THEN DO 50 I=1,3 DO 50 J=1,3 50 A(I,J)=-A(I,J) ENDIF C C Transform the force constant if(LNAMES)then c retrieve FC of fragment IPO c c load only if not in memory already if(IPOOLD.EQ.0.or.IPOOLD.NE.IPO)then NE=80 DO 1500 NE=80,1,-1 1500 if(NAME(IPO)(NE:NE).ne.' ')goto 1600 1600 NATS=NATNAME(IPO) IPOOLD=IPO endif ido=(IPO-1)*81*NATS_max**3 ELSE ido=0 ENDIF IAS=INDBIG(IPO,IA) JAS=INDBIG(IPO,JA) DO 361 IX=1,3 DO 361 JX=1,3 DO 361 KX=1,3 DO 361 LX=1,3 c c TYPE AABB: i1=NAT*((9*IA+3*IX+JX-13)*NAT*9+3*(3*JA-4+KX))+3*JA-3+LX FOLD=db(i1) C Zero out at the start of the summation IF(ISUM.EQ.1)db(i1)=0.0d0 FNEW=0.0d0 DO 371 IIX=1,3 SIA=A(IX,IIX) DO 371 JJX=1,3 SJB=SIA*A(JX,JJX) DO 371 KKX=1,3 SKC=SJB*A(KX,KKX) DO 371 LLX=1,3 SLD=SKC*A(LX,LLX) is=NATS*((9*IAS+3*IIX+JJX-13)*NATS*9+3*(3*JAS-4+KKX))+3*JAS-3+LLX 371 FNEW=FNEW+SLD*ds(ido+is) IF(LWR)WRITE(3,3002)IX,JX,KX,LX,FOLD,FNEW,IIT,TOL,WF(IP) 3002 FORMAT('d',4I2,2F11.5,I5,F10.7,F7.4) db(i1)=FNEW*WF(IP)+db(i1) IF(IA.NE.JA)then c TYPE BBAA: i2=NAT*((9*JA+3*KX+LX-13)*NAT*9+3*(3*IA-4+IX))+3*IA-3+JX db(i2)=db(i1) c c TYPE AAAB: i1=NAT*((9*IA+3*IX+JX-13)*NAT*9+3*(3*IA-4+KX))+3*JA-3+LX FOLD=db(i1) C Zero out at the start of the summation IF(ISUM.EQ.1)db(i1)=0.0d0 FNEW=0.0d0 DO 372 IIX=1,3 SIA=A(IX,IIX) DO 372 JJX=1,3 SJB=SIA*A(JX,JJX) DO 372 KKX=1,3 SKC=SJB*A(KX,KKX) DO 372 LLX=1,3 SLD=SKC*A(LX,LLX) is=NATS*((9*IAS+3*IIX+JJX-13)*NATS*9+3*(3*IAS-4+KKX))+3*JAS-3+LLX 372 FNEW=FNEW+SLD*ds(ido+is) IF(LWR)WRITE(3,3002)IX,JX,KX,LX,FOLD,FNEW,IIT,TOL,WF(IP) db(i1)=FNEW*WF(IP)+db(i1) i2=NAT*((9*IA+3*IX+JX-13)*NAT*9+3*(3*JA-4+LX))+3*IA-3+KX db(i2)=db(i1) c c TYPE BBBA: i1=NAT*((9*JA+3*IX+JX-13)*NAT*9+3*(3*JA-4+KX))+3*IA-3+LX FOLD=db(i1) C Zero out at the start of the summation IF(ISUM.EQ.1)db(i1)=0.0d0 FNEW=0.0d0 DO 373 IIX=1,3 SIA=A(IX,IIX) DO 373 JJX=1,3 SJB=SIA*A(JX,JJX) DO 373 KKX=1,3 SKC=SJB*A(KX,KKX) DO 373 LLX=1,3 SLD=SKC*A(LX,LLX) is=NATS*((9*JAS+3*IIX+JJX-13)*NATS*9+3*(3*JAS-4+KKX))+3*IAS-3+LLX 373 FNEW=FNEW+SLD*ds(ido+is) IF(LWR)WRITE(3,3002)IX,JX,KX,LX,FOLD,FNEW,IIT,TOL,WF(IP) db(i1)=FNEW*WF(IP)+db(i1) i2=NAT*((9*JA+3*IX+JX-13)*NAT*9+3*(3*IA-4+LX))+3*JA-3+KX db(i2)=db(i1) endif 361 continue 45 CONTINUE c loop over overlaps C 20 continue endif 30 CONTINUE c IA,JA RETURN END SUBROUTINE IMPROVE34(i34,cb,cs,RB,RS,NAT,NATS,IND,INDBIG,NO, 1IBONDS,NBOB,LWR,IWG,LSTRICT,LINV,LNAMES,LHALTONERROR,CUTOFF, 2okind,NATS_max) implicit none INTEGER*4 MX,MXB,IOX PARAMETER (MX=4000,MXB=140,IOX=1500) INTEGER*4 okind,IND(MX),INDBIG(IOX,MX),INDLS,LLX,NS,N, 4IBONDS(4,MX),NBOB(MX),IOL(IOX),ISUM,I,IA,II,i1,i1b,i2b, 1IIOII,IIX,INDIB,INDIS,INDJS,IO,i34,ix1,ix2,IOIJ,IP,IPASS, 1IPO,IWG,IX,J,JA,IERC,I2,J2,k,IPOOLD,JSTART,JX,NAT,NATIO,NATS, 1NE,NF,NO,INDJB,JJX,ipass_o(IOX),ierc_o(IOX),itu_o(IOX), 1ind_o(IOX,MX),nats_o(IOX),INDKB,INDKS,KA,KKX,KSTART,KX,LX, 1NATS_max,ido real*8 RB(3,MX),RS(3,NATS),t1(81),t2(81),WF(IOX), 1TIJX,OTIJX,TIJY,OTIJY,TIJZ,OTIJZ, 1cb(*),cs(*),DISTM,SUM,FNEW,DIST,FOLD,CTF,CUTOFF, 6RM(IOX,3,3),xi,yi,zi,xj,yj,zj,xk,yk,zk,SIA,SIB,SJC,SKD LOGICAL LWR,LSTRICT,LINV,LNAMES,LHALTONERROR real*8 A,XS,XB,TOL,RMS integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS real*8 aname,pname,aaname,gname,alphaname,rname,polname,pol, 1gpo,apo,apolname,gpolname,gradname,qsname integer*4 natname character*80 NAME(IOX) common/cname/aname(IOX,MX,3,3),polname(IOX,3,3),pol(3,3), 1pname(IOX,MX,3,3),aaname(IOX,3*MX,3,3,3),gname(IOX,3*MX,3,3), 2alphaname(IOX,3*MX,3,3),rname(IOX,3,MX),qsname(IOX,MX,3,3,3), 3gpo(3,3),apo(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3), 2gradname(IOX,3*MX),natname(IOX),NAME c IF(NO.EQ.0)call report('NO=0') if(i34.eq.3)then WRITE(3,*)' Cubic constants requested' else if(i34.eq.4)then WRITE(3,*)' Quartic constants requested' else if(i34.eq.5)then WRITE(3,*)' Three-atomic quartic constants requested' else write(3,*)i34 call report('unknown option') endif endif endif IF(LSTRICT)WRITE(3,*)' Only transformed atoms in the beds' IF(IWG.EQ.1)WRITE(3,*)' Multiple overlaps will be weighted' IF(LINV)WRITE(3,*)' Smaller fragment as optical antipode' IF(LNAMES)WRITE(3,*)' Named small fragments' IF(LHALTONERROR)WRITE(3,*)' Will stop on error' if(okind.eq.0)write(3,*)'selection on centers of masses' if(okind.eq.1)write(3,*)'selection on RMS distances' if(okind.eq.2)write(3,*)'selection on molecular RMS distances' WRITE(3,3999)CUTOFF 3999 format('Distance cutoff ',f12.2,' A') ctf=CUTOFF**2 IOIJ=0 IPOOLD=0 NS=3*NATS N=3*NAT C C Look at all atom triads i-j-k: DO 30 IA=1,NAT write(6,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) if(i34.eq.3)then jstart=IA else jstart=1 endif xi=RB(1,IA) yi=RB(2,IA) zi=RB(3,IA) DO 30 JA=jstart,NAT xj=RB(1,JA) yj=RB(2,JA) zj=RB(3,JA) if((xi-xj)**2+(yi-yj)**2+(zi-zj)**2.lt.ctf)then DO 20 KA=JA,NAT xk=RB(1,KA) yk=RB(2,KA) zk=RB(3,KA) if((xi-xk)**2+(yi-yk)**2+(zi-zk)**2.lt.ctf)then if((xj-xk)**2+(yj-yk)**2+(zj-zk)**2.lt.ctf)then if(LWR)write(3,3008)IA,JA,KA 3008 format(/,' Atoms ',3I6) C C Look at all overlaps io, find which fits best ijk c c for okind=2, look at the overlaps only once c cccccccccccccccccccccccccccccccccccccccccccccccccccc IF(okind.ne.2.or.(IA.eq.1.and.JA.EQ.1))THEN C C Calculate center of i/j/k radiusvectors: TIJX=(xi+xj+xk)/3.0d0 TIJY=(yi+yj+yk)/3.0d0 TIJZ=(zi+zj+zk)/3.0d0 DISTM=1.0d5 c c NO ... total number of overlaps c NF ... number of overlaps containing IA JA KA NF=0 DO 21 IO=1,NO C C Check, if io contains ijk: IF(INDBIG(IO,IA).GT.0.AND. 1 INDBIG(IO,JA).GT.0.AND.INDBIG(IO,KA).GT.0)THEN NF=NF+1 IOL(NF)=IO C C Calculate center of the overlap and distance from i-j center c i.........TIJ.............j c ----------OTIJ--------------------- OTIJX=0.0D0 OTIJY=0.0D0 OTIJZ=0.0D0 NATIO=0 DO 24 II=1,NAT IIOII=INDBIG(IO,II) C If Lstrict, take only transformed atoms IF(LSTRICT.AND.IIOII.LT.0)IIOII=0 IIOII=ABS(IIOII) IF(IIOII.NE.0)THEN NATIO=NATIO+1 OTIJX=OTIJX+RB(1,II) OTIJY=OTIJY+RB(2,II) OTIJZ=OTIJZ+RB(3,II) ENDIF 24 CONTINUE OTIJX=OTIJX/DBLE(NATIO) OTIJY=OTIJY/DBLE(NATIO) OTIJZ=OTIJZ/DBLE(NATIO) c c If required, calculate RMS overlap and use this criterium instead c of the distances,remember rotation matrix (RM) if(okind.gt.0)then call CRMS34(IPASS,NAT,INDBIG,IOX,MX,IO,IA,JA,KA,NBOB,okind, 1 IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) do 55 IX=1,3 do 55 JX=1,3 55 RM(NF,IX,JX)=A(IX,JX) ipass_o(NF)=IPASS WF(NF)=DSQRT(RMS) ierc_o(NF)=IERC itu_o(NF)=ITU do 58 I=1,ITU 58 ind_o(NF,I)=IND(I) do 59 I=1,3 do 59 J=1,3 59 RM(NF,I,J)=A(I,J) nats_o(NF)=NATS IF(RMS.LT.DISTM)THEN DISTM=RMS IOIJ=IO ENDIF else c Use distances: DIST=(OTIJX-TIJX)**2+(OTIJY-TIJY)**2+(OTIJZ-TIJZ)**2 IF(DIST.LT.DISTM)THEN DISTM=DIST IOIJ=IO ENDIF WF(NF)=DSQRT(DIST) endif ENDIF 21 CONTINUE c WRITE(3,3003)IA,JA,KA,NF,IOIJ 3003 FORMAT(3I4,':',I4,' overlaps found; best is # ',I4) C C Calculate the weighting factors: SUM=0.0d0 DO 43 I=1,NF 43 IF(WF(I).GT.1.0d-9)SUM=SUM+1.0d0/WF(I) DO 42 I=1,NF IF(WF(I).LE.1.0d-9)THEN WF(1)=1.0d0 NF=1 GOTO 44 ENDIF WF(I)=1.0d0/WF(I)/SUM 42 CONTINUE 44 CONTINUE C c cccccccccccccccccccccccccccccccccccccccccccccccccc ENDIF c C Loop over possible overlaps: ISUM=0 DO 45 IP=1,NF IPO=IOL(IP) IF(IWG.EQ.0)THEN IF(IPO.NE.IOIJ)GOTO 45 WF(IP)=1.0d0 ENDIF ISUM=ISUM+1 C C Find the atomic groups to be placed on together for IA-JA: c and the rotation matrix if(okind.eq.0)then call CRMS34(IPASS,NAT,INDBIG,IOX,MX,IPO, 1 IA,JA,KA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) else IPASS=ipass_o(IP) IERC=ierc_o(IP) ITU=itu_o(IP) do 56 I=1,ITU 56 IND(I)=ind_o(IP,I) do 57 I=1,3 do 57 J=1,3 57 A(I,J)=RM(IP,I,J) NATS=nats_o(IP) endif IF(LWR)WRITE(3,*)IPASS,' passes' if(IERC.ne.0)goto 20 WRITE(3,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU) WRITE(3,3005)(INDBIG(IPO,IND(I)),I=1,ITU) 3004 FORMAT('Overlap ',I2,': weighting factor: ',F7.4,';', 1 I4,' atom bed:',/, 2 ' Big : ',20I3) 3005 FORMAT( ' Small: ',20I3,/) if(LNAMES)write(3,3399)NAME(IPO) 3399 FORMAT(A80) C IF(LINV)THEN DO 50 I=1,3 DO 50 J=1,3 50 A(I,J)=-A(I,J) ENDIF C C Transform the force constant if(LNAMES)then c retrieve FC of fragment IPO c c load only if not in memory already if(IPOOLD.EQ.0.or.IPOOLD.NE.IPO)then NE=80 DO 1500 NE=80,1,-1 1500 if(NAME(IPO)(NE:NE).ne.' ')goto 1600 1600 NATS=NATNAME(IPO) IPOOLD=IPO endif ido=(IPO-1)*81*NATS_max**3 ELSE ido=0 ENDIF NS=3*NATS if(i34.eq.3)then c c rewrite cs into temporary array do 60 IIX=1,3 INDIS=3*INDBIG(IPO,IA)-3+IIX do 60 JJX=1,3 INDJS=3*INDBIG(IPO,JA)-3+JJX do 60 KKX=1,3 INDKS=3*INDBIG(IPO,KA)-3+KKX 60 t1(KKX+3*(JJX-1+3*(IIX-1)))=cs(ido+ 1 INDIS+NS*(INDJS-1+NS*(INDKS-1))) c c transform first index do 61 IX=1,3 do 61 JJX=1,3 do 61 KKX=1,3 t2(KKX+3*(JJX-1+3*(IX-1)))=0.0d0 do 61 IIX=1,3 61 t2(KKX+3*(JJX-1+3*(IX-1)))=t2(KKX+3*(JJX-1+3*(IX-1)))+ 1 t1(KKX+3*(JJX-1+3*(IIX-1)))*A(IX,IIX) c transform second index do 62 IX=1,3 do 62 JX=1,3 do 62 KKX=1,3 t1(KKX+3*(JX-1+3*(IX-1)))=0.0d0 do 62 JJX=1,3 62 t1(KKX+3*(JX-1+3*(IX-1)))=t1(KKX+3*(JX-1+3*(IX-1)))+ 1 t2(KKX+3*(JJX-1+3*(IX-1)))*A(JX,JJX) c transform third index do 63 IX=1,3 do 63 JX=1,3 do 63 KX=1,3 t2(KX+3*(JX-1+3*(IX-1)))=0.0d0 do 63 KKX=1,3 63 t2(KX+3*(JX-1+3*(IX-1)))=t2(KX+3*(JX-1+3*(IX-1)))+ 1 t1(KKX+3*(JX-1+3*(IX-1)))*A(KX,KKX) c rewrite into new field DO 36 IX=1,3 INDIB=3*IA-3+IX JSTART=1 IF(IA.EQ.JA)JSTART=IX DO 36 JX=JSTART,3 INDJB=3*JA-3+JX KSTART=1 IF(JA.EQ.KA)KSTART=JX DO 36 KX=KSTART,3 INDKB=3*KA-3+KX i1=INDIB+N*(INDJB-1+N*(INDKB-1)) FOLD=cb(i1) C Zero out at the start of the summation IF(ISUM.EQ.1)cb(i1)=0.0d0 FNEW=t2(KX+3*(JX-1+3*(IX-1))) IF(LWR)WRITE(3,3001)IX,JX,KX,FOLD,FNEW,IIT,TOL,WF(IP) 3001 FORMAT('c',3I5,2F11.5,I5,F10.7,F7.4) cb(i1 )=cb(i1)+FNEW*WF(IP) cb(INDIB+N*(INDKB-1+N*(INDJB-1)))=cb(i1) cb(INDJB+N*(INDKB-1+N*(INDIB-1)))=cb(i1) cb(INDJB+N*(INDIB-1+N*(INDKB-1)))=cb(i1) cb(INDKB+N*(INDJB-1+N*(INDIB-1)))=cb(i1) 36 cb(INDKB+N*(INDIB-1+N*(INDJB-1)))=cb(i1) endif if(i34.eq.4)then DO 361 IX=1,3 c double index IX INDIB=3*IA-3+IX DO 361 JX=1,3 INDJB=3*JA-3+JX KSTART=1 IF(JA.EQ.KA)KSTART=JX DO 361 KX=KSTART,3 INDKB=3*KA-3+KX i1=INDIB+N*(INDJB-1+N*(INDKB-1)) FOLD=cb(i1) C Zero out at the start of the summation IF(ISUM.EQ.1)cb(i1)=0.0d0 FNEW=0.0d0 DO 371 IIX=1,3 INDIS=3*INDBIG(IPO,IA)-3+IIX SIA=A(IX,IIX) DO 371 JJX=1,3 INDJS=3*INDBIG(IPO,IA)-3+JJX SIB=SIA*A(IX,JJX) DO 371 KKX=1,3 INDKS=3*INDBIG(IPO,JA)-3+KKX SJC=SIB*A(JX,KKX) DO 371 LLX=1,3 INDLS=3*INDBIG(IPO,KA)-3+LLX SKD=SJC*A(KX,LLX) c FNEW=FNEW+SKD*D(INDIS,INDJS,INDKS,INDLS) if(INDIS.eq.INDJS)then FNEW=FNEW+SKD*cs(ido+INDIS+NS*(INDKS-1+NS*(INDLS-1))) else if(INDIS.eq.INDKS)then FNEW=FNEW+SKD*cs(ido+INDIS+NS*(INDJS-1+NS*(INDLS-1))) else if(INDIS.eq.INDLS)then FNEW=FNEW+SKD*cs(ido+INDIS+NS*(INDJS-1+NS*(INDKS-1))) else if(INDJS.eq.INDKS)then FNEW=FNEW+SKD*cs(ido+INDJS+NS*(INDIS-1+NS*(INDLS-1))) else if(INDJS.eq.INDLS)then FNEW=FNEW+SKD*cs(ido+INDJS+NS*(INDIS-1+NS*(INDKS-1))) else if(INDKS.eq.INDLS)then FNEW=FNEW+SKD*cs(ido+INDKS+NS*(INDIS-1+NS*(INDJS-1))) endif endif endif endif endif endif 371 continue IF(LWR)WRITE(3,3002)IX,JX,KX,FOLD,FNEW,IIT,TOL,WF(IP) 3002 FORMAT('d',3I5,2F11.5,I5,F10.7,F7.4) cb(i1 )=cb(i1)+FNEW*WF(IP) 361 cb(INDIB+N*(INDKB-1+N*(INDJB-1)))=cb(i1) endif if(i34.eq.5)then c c rewrite cs into temporary array do 601 IIX=1,3 do 601 JJX=1,3 I2=(9*INDBIG(IPO,IA)+3*IIX+JJX-13)*NATS*9 do 601 KKX=1,3 J2=3*(INDBIG(IPO,JA)-1)+KKX do 601 LLX=1,3 i1=NATS*(I2+3*(J2-1))+3*(INDBIG(IPO,KA)-1)+LLX 601 t1(LLX+3*(KKX-1+3*(JJX-1+3*(IIX-1))))=cs(ido+i1) c c transform first index do 611 IX=1,3 do 611 JJX=1,3 do 611 KKX=1,3 do 611 LLX=1,3 611 t2(LLX+3*(KKX-1+3*(JJX-1+3*(IX -1))))= 1 t1(LLX+3*(KKX-1+3*(JJX-1+3*(1 -1))))*A(IX,1)+ 1 t1(LLX+3*(KKX-1+3*(JJX-1+3*(2 -1))))*A(IX,2)+ 1 t1(LLX+3*(KKX-1+3*(JJX-1+3*(3 -1))))*A(IX,3) c transform second index do 621 IX=1,3 do 621 JX=1,3 do 621 KKX=1,3 do 621 LLX=1,3 621 t1(LLX+3*(KKX-1+3*(JX -1+3*(IX-1))))= 1 t2(LLX+3*(KKX-1+3*(1 -1+3*(IX-1))))*A(JX,1)+ 1 t2(LLX+3*(KKX-1+3*(2 -1+3*(IX-1))))*A(JX,2)+ 1 t2(LLX+3*(KKX-1+3*(3 -1+3*(IX-1))))*A(JX,3) c transform third index do 631 IX=1,3 do 631 JX=1,3 do 631 KX=1,3 do 631 LLX=1,3 631 t2(LLX+3*(KX -1+3*(JX -1+3*(IX -1))))= 1 t1(LLX+3*(1 -1+3*(JX -1+3*(IX -1))))*A(KX,1)+ 1 t1(LLX+3*(2 -1+3*(JX -1+3*(IX -1))))*A(KX,2)+ 1 t1(LLX+3*(3 -1+3*(JX -1+3*(IX -1))))*A(KX,3) c transform fourth index do 632 IX=1,3 do 632 JX=1,3 do 632 KX=1,3 do 632 LX=1,3 632 t1(LX +3*(KX -1+3*(JX -1+3*(IX -1))))= 1 t2(1 +3*(KX -1+3*(JX -1+3*(IX -1))))*A(LX,1)+ 1 t2(2 +3*(KX -1+3*(JX -1+3*(IX -1))))*A(LX,2)+ 1 t2(3 +3*(KX -1+3*(JX -1+3*(IX -1))))*A(LX,3) c rewrite into new field c indices XI1,IX2 to IA: DO 362 IX1=1,3 DO 362 IX2=1,3 DO 362 JX=1,3 J=3*(JA-1)+JX DO 362 KX=1,3 K=3*(KA-1)+KX c i1b- index for the BIG system i1b=NAT*((9*IA+3*IX1+IX2-13)*nat*9+3*(J-1))+K FOLD=cb(i1b) C Zero out at the start of the summation IF(ISUM.EQ.1)cb(i1b)=0.0d0 FNEW=t1(KX+3*(JX-1+3*(IX2-1+3*(IX1-1)))) IF(LWR)WRITE(3,3042)IX1,IX2,JX,KX,FOLD,FNEW,IIT,TOL,WF(IP) 3042 FORMAT(4I5,2F11.5,I5,F10.7,F7.4) cb(i1b)=FNEW*WF(IP)+cb(i1b) c i2b symmetry-related BIG index JA <-> KA i2b=NAT*((9*IA+3*IX1+IX2-13)*nat*9+3*(K-1))+J 362 cb(i2b)=cb(i1b) endif 45 CONTINUE c loop over overlaps C endif endif 20 CONTINUE c KA endif 30 CONTINUE c IA,JA RETURN END SUBROUTINE IMPROVEMPT(np,INDBIG,IOX,MX,IO,LNAMES,NATNAME, 1RS,RB,IND,rname,NATS,LINV,NO,NAT,lq,ND,alt,al,lmut,wmin,wmax, 1nel) implicit none INTEGER*4 np,iw,ix,jx,kx,IO,IOP,IOX,MX,NATNAME(*), 1INDBIG(IOX,MX),IND(*),NO,noo,nu,iao,NAT,ND(*),ic,IiO, 1ii,jj,i,j,NATS,MXB,istart,nel parameter (MXB=140) real*8 bohr,bigpolr(3,3),bigpoli(3,3),biggpolr(3,3), 1biggpoli(3,3),bigapolr(3,3,3),bigapoli(3,3,3),RS(3,MX), 1RB(3,MX),rname(IOX,3,MX),su,alt(*), 1polr(3,3),poli(3,3),gpolr(3,3),gpoli(3,3), 1apolr(3,3,3),apoli(3,3,3),al(*), 1wmin,wmax c alt resultant total big polarizability c al input overlap polarizabilities real*8, allocatable::coo(:,:),zsum(:),opolr(:,:,:), 1opoli(:,:,:),gopolr(:,:,:),gopoli(:,:,:),Pi(:,:), 2aopolr(:,:,:,:),aopoli(:,:,:,:),am(:,:,:),X(:,:), 3Pr(:,:) logical LNAMES,LINV,lq,lmut real*8 A,XS,XB,TOL,RMS integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS c allocate(coo(NO,3),opolr(NO,3,3),opoli(NO,3,3), 1zsum(NO),gopolr(NO,3,3),gopoli(NO,3,3), 1aopolr(NO,3,3,3),aopoli(NO,3,3,3),am(NO,3,3)) c bohr=0.529177D0 write(6,8000)wmin,wmax,np,NO,LMUT,NAT 8000 format(/,/,' Dynamic polarizability transfer-MPT',/,/, 1 ' wmin:',f10.1,' wmax:',f10.1,' np :',i10,/, 2 ' NO :',i10, ' LMUT:',L10, ' NAT:',i10,/) c c loop over overlaps: c oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo DO 615 IO=1,NO write(6,6010)IO 6010 format(' overlap',i6,', transformation matrix:') c make the transformation matrix for whole overlap: call CRMSo(NAT,INDBIG,IOX,MX,IO,LNAMES,natname,RS,RB,IND, 2rname,NATS,LINV) do 610 ix=1,3 do 610 jx=1,3 610 am(IO,ix,jx)=A(ix,jx) write(6,6019)A 6019 format(3(3f10.3,/)) c c calculate overlap center: do 53 ix=1,3 53 coo(IO,ix)=0.0d0 noo=0 nu=0 zsum(IO)=0.0d0 do 52 iao=1,NAT if(INDBIG(IO,iao).GT.0)then if(lq)then su=dble(ND(iao)) else su=1.0d0 endif zsum(IO)=zsum(IO)+su noo=noo+1 do 54 ix=1,3 54 coo(IO,ix)=coo(IO,ix)+RB(ix,iao)*su c c look at other overlaps, and find out if this atoms is part of them: ic=0 do 552 IOP=1,NO 552 if(IOP.ne.IO.and.INDBIG(IOP,iao).ne.0)ic=ic+1 if(ic.eq.0)nu=nu+1 endif 52 continue if(noo.gt.0)then do 55 ix=1,3 55 coo(IO,ix)=coo(IO,ix)/zsum(IO)/bohr endif 615 write(6,6009)(coo(IO,ix)*bohr,ix=1,3),noo,NATS,nu 6009 format(' Center: ',3f12.3,' A',/, 1i5,' atoms transferred of',i5,',',i5,' unique') c oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo allocate(X(NO*24,NO*24),Pr(NO*24,NO*24),Pi(NO*24,NO*24)) do 63 i=1,NO*24 do 63 j=1,NO*24 Pr(i,j)=0.0d0 Pi(i,j)=0.0d0 63 X(i,j)=0.0d0 c ************* Calculate general distance tensor: ********* call doX(NO,coo,X) c c loop over frequencies do 1 iw=1,np write(6,*)iw c total polarizability do 607 ix=1,3 do 607 jx=1,3 biggpolr(ix,jx)=0.0d0 bigpolr(ix,jx)=0.0d0 biggpoli(ix,jx)=0.0d0 bigpoli(ix,jx)=0.0d0 do 607 kx=1,3 bigapolr(ix,jx,kx)=0.0d0 607 bigapoli(ix,jx,kx)=0.0d0 DO 551 IO=1,NO do 613 ix=1,3 do 613 jx=1,3 613 A(ix,jx)=am(IO,ix,jx) c c if named fragments, fill alpha, G, A for each overlap: if(lnames)then IiO=IO else IiO=1 endif ii=0 jj=0 do 2 ix=1,3 do 2 jx=1,3 ii=ii+1 polr(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii ) poli(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii+ 9) gpolr(ix,jx)=al(nel*((IiO-1)*np+iw-1)+ii+18) gpoli(ix,jx)=al(nel*((IiO-1)*np+iw-1)+ii+27) do 2 kx=1,3 jj=jj+1 apolr(jx,kx,ix)=al(nel*((IiO-1)*np+iw-1)+jj+36) 2 apoli(jx,kx,ix)=al(nel*((IiO-1)*np+iw-1)+jj+63) c c rotate polarizabilities and write to new strings: call rpp( NO,IO, opolr, polr,A) call rpp( NO,IO,gopolr,gpolr,A) call rppp(NO,IO,aopolr,apolr,A) call rpp( NO,IO, opoli, poli,A) call rpp( NO,IO,gopoli,gpoli,A) 551 call rppp(NO,IO,aopoli,apoli,A) c c ************* Calculate general polarizability tensor: ********* do 7 i=1,NO istart=24*(i-1)+1 do 7 ix=1,3 ii=0 do 7 jx=1,3 Pr(istart+ix-1,istart+jx-1)=opolr(i,ix,jx) Pr(istart+ix+2,istart+jx+2)=opolr(i,ix,jx) Pr(istart+ix-1,istart+jx+8)=gopolr(i,ix,jx) Pr(istart+jx+8,istart+ix-1)=gopolr(i,ix,jx) Pr(istart+ix+2,istart+jx+5)=-gopolr(i,ix,jx) Pr(istart+jx+5,istart+ix+2)=-gopolr(i,ix,jx) Pi(istart+ix-1,istart+jx-1)=opoli(i,ix,jx) Pi(istart+ix+2,istart+jx+2)=opoli(i,ix,jx) Pi(istart+ix-1,istart+jx+8)=gopoli(i,ix,jx) Pi(istart+jx+8,istart+ix-1)=gopoli(i,ix,jx) Pi(istart+ix+2,istart+jx+5)=-gopoli(i,ix,jx) Pi(istart+jx+5,istart+ix+2)=-gopoli(i,ix,jx) do 7 kx=jx,3 ii=ii+1 Pr(istart+ix- 1,istart+ii+11)=aopolr(i,ix,jx,kx)/3.0d0 Pr(istart+ii+11,istart+ix- 1)=aopolr(i,ix,jx,kx)/3.0d0 Pr(istart+ix+ 2,istart+ii+17)=aopolr(i,ix,jx,kx)/3.0d0 Pr(istart+ii+17,istart+ix+ 2)=aopolr(i,ix,jx,kx)/3.0d0 Pi(istart+ix- 1,istart+ii+11)=aopoli(i,ix,jx,kx)/3.0d0 Pi(istart+ii+11,istart+ix- 1)=aopoli(i,ix,jx,kx)/3.0d0 Pi(istart+ix+ 2,istart+ii+17)=aopoli(i,ix,jx,kx)/3.0d0 7 Pi(istart+ii+17,istart+ix+ 2)=aopoli(i,ix,jx,kx)/3.0d0 c a 0 0 G A/3 0 c 0 a -G 0 0 A/3 c 0 -Gt 0 0 0 0 c Gt 0 0 0 0 0 c A/3 0 0 0 0 0 c 0 A/3 0 0 0 0 c c get redressed polarizabilies, common origin call MDP(coo,x,Pr,Pi,NO,bigpolr,bigpoli,biggpolr,biggpoli, 1bigapolr,bigapoli) c record to total polarizability at this frequency: ii=0 jj=0 do 3 ix=1,3 do 3 jx=1,3 ii=ii+1 alt(nel*(iw-1)+ii )=bigpolr(ix,jx) alt(nel*(iw-1)+ii+ 9)=bigpoli(ix,jx) alt(nel*(iw-1)+ii+18)=biggpolr(ix,jx) alt(nel*(iw-1)+ii+27)=biggpoli(ix,jx) do 3 kx=1,3 jj=jj+1 alt(nel*(iw-1)+jj+36)=bigapolr(jx,kx,ix) 3 alt(nel*(iw-1)+jj+63)=bigapoli(jx,kx,ix) 1 continue c iw-frequency RETURN END SUBROUTINE IMPROVESOS(np,INDBIG,IOX,MX,IO,LNAMES,NATNAME, 1RS,RB,IND,rname,NATS,LINV,NO,NAT,lq,ND,alt,al,lmut,wmin,wmax, 1nel,CUTOFF,SOSNM) implicit none INTEGER*4 np,iw,ix,jx,kx,IO,IOP,IOX,MX,NATNAME(*),nel, 1INDBIG(IOX,MX),IND(*),NO,noo,nu,iao,NAT,ND(*),ic,IiO, 1ii,jj,i,j,NATS,MXB parameter (MXB=140) real*8 bohr,bigpolr(3,3),bigpoli(3,3),biggpolr(3,3), 1bigpolvr(3,3),bigpolvi(3,3),CUTOFF,cut2, 1biggpoli(3,3),bigapolr(3,3,3),bigapoli(3,3,3),RS(3,MX), 1RB(3,MX),rname(IOX,3,MX),su,temr(3,3),temi(3,3),alt(*), 1polr(3,3),poli(3,3),gpolr(3,3),gpoli(3,3), 1apolr(3,3,3),apoli(3,3,3),al(*),rt(3),rt2,tt(3,3), 1ttt(3,3,3),tar(3,3),tai(3,3),tcr1(3,3),tcr2(3,3), 1tci1(3,3),tci2(3,3),tcr(3,3),tci(3,3),pi,sf,taii(3,3), 1tair(3,3),tajr(3,3),taji(3,3),atar1(3,3),atar2(3,3), 1atai1(3,3),atai2(3,3),ataii(3,3),ataji(3,3),atajr(3,3), 1tvi(3),tvr(3),sui,sur,atair(3,3),wmin,wmax,polvr(3,3),polvi(3,3), 1x,y,z,xx,yy,zz,w,wau,wnm,dw,CM,wcm c alt resultant total big polarizability c al input overlap polarizabilities real*8, allocatable::coo(:,:),zsum(:),opolr(:,:,:), 1opoli(:,:,:),gopolr(:,:,:),gopoli(:,:,:),opolvr(:,:,:), 2aopolr(:,:,:,:),aopoli(:,:,:,:),am(:,:,:),opolvi(:,:,:) logical LNAMES,LINV,lq,lmut,SOSNM real*8 A,XS,XB,TOL,RMS integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS c allocate(coo(NO,3),opolr(NO,3,3),opoli(NO,3,3), 1opolvr(NO,3,3),opolvi(NO,3,3), 1zsum(NO),gopolr(NO,3,3),gopoli(NO,3,3), 1aopolr(NO,3,3,3),aopoli(NO,3,3,3),am(NO,3,3)) c bohr=0.529177D0 cut2=CUTOFF**2 write(6,8000)wmin,wmax,np,NO,LMUT,NAT,CUTOFF 8000 format(/,/,' Dynamic polarizability transfer',/,/, 1 ' wmin:',f10.1,' wmax:',f10.1,' np :',i10,/, 2 ' NO :',i10, ' LMUT:',L10, ' NAT:',i10,/, 2 ' CUT :',f10.1,/) c c loop over overlaps: c oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo DO 615 IO=1,NO write(6,6010)IO 6010 format(' overlap',i6,', transformation matrix:') c make the transformation matrix for whole overlap: call CRMSo(NAT,INDBIG,IOX,MX,IO,LNAMES,natname,RS,RB,IND, 2rname,NATS,LINV) do 610 ix=1,3 do 610 jx=1,3 610 am(IO,ix,jx)=A(ix,jx) write(6,6019)A 6019 format(3(3f10.3,/)) c c calculate overlap center: do 53 ix=1,3 53 coo(IO,ix)=0.0d0 noo=0 nu=0 zsum(IO)=0.0d0 do 52 iao=1,NAT if(INDBIG(IO,iao).GT.0)then if(lq)then su=dble(ND(iao)) else su=1.0d0 endif zsum(IO)=zsum(IO)+su noo=noo+1 do 54 ix=1,3 54 coo(IO,ix)=coo(IO,ix)+RB(ix,iao)*su c c look at other overlaps, and find out if this atoms is part of them: ic=0 do 552 IOP=1,NO 552 if(IOP.ne.IO.and.INDBIG(IOP,iao).ne.0)ic=ic+1 if(ic.eq.0)nu=nu+1 endif 52 continue if(noo.gt.0)then do 55 ix=1,3 55 coo(IO,ix)=coo(IO,ix)/zsum(IO)/bohr endif 615 write(6,6009)(coo(IO,ix)*bohr,ix=1,3),noo,NATS,nu 6009 format(' Center: ',3f12.3,' A',/, 1i5,' atoms transferred of',i5,',',i5,' unique') c oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo c CM=219474.63d0 dw=(wmax-wmin)/dble(np-1) w=wmin-dw c loop over frequencies do 1 iw=1,np w=w+dw if(SOSNM)then wau=(1.0d7/w)/CM else wau=w/CM endif wcm=wau*CM wnm=1.0d7/wcm c total polarizability-zero out: call zp(biggpolr,biggpoli,bigpolr,bigpoli,bigpolvr,bigpolvi, 1bigapolr,bigapoli) DO 551 IO=1,NO do 613 ix=1,3 do 613 jx=1,3 613 A(ix,jx)=am(IO,ix,jx) c c if named fragments, fill alpha, G, A for each overlap: if(lnames)then IiO=IO else IiO=1 endif c transcript polarizabilities: call trpol(polr,poli,polvr,polvi,gpolr,gpoli, 1apolr,apoli,al,nel,np,IiO,iw) c c rotate polarizabilities and write to new strings: call rpp( NO,IO, opolr, polr,A) call rpp( NO,IO, opolvr, polvr,A) call rpp( NO,IO,gopolr,gpolr,A) call rppp(NO,IO,aopolr,apolr,A) call rpp( NO,IO, opoli, poli,A) call rpp( NO,IO, opolvi, polvi,A) call rpp( NO,IO,gopoli,gpoli,A) call rppp(NO,IO,aopoli,apoli,A) c c construct alpha,G,A of all system: c plain sum: do 5519 IX=1,3 do 5519 JX=1,3 bigpolr( IX,JX) =bigpolr( IX,JX)+opolr( IO,IX,JX) bigpoli( IX,JX) =bigpoli( IX,JX)+opoli( IO,IX,JX) bigpolvr(IX,JX) =bigpolvr(IX,JX)+opolvr(IO,IX,JX) bigpolvi(IX,JX) =bigpolvi(IX,JX)+opolvi(IO,IX,JX) biggpolr(IX,JX) =biggpolr(IX,JX)+gopolr(IO,IX,JX) biggpoli(IX,JX) =biggpoli(IX,JX)+gopoli(IO,IX,JX) do 5519 KX=1,3 bigapoli(IX,JX,KX)=bigapoli(IX,JX,KX)+aopoli(IO,IX,JX,KX) 5519 bigapolr(IX,JX,KX)=bigapolr(IX,JX,KX)+aopolr(IO,IX,JX,KX) c c polarizability part for origin dependence of G call odg(biggpolr,coo(IO,1),coo(IO,2),coo(IO,3),opolr,IO,NO) call odg(biggpoli,coo(IO,1),coo(IO,2),coo(IO,3),opoli,IO,NO) c c polarizability part for origin dependence of A: call oda(bigapolr,coo(IO,1),coo(IO,2),coo(IO,3),opolr,IO,NO) call oda(bigapoli,coo(IO,1),coo(IO,2),coo(IO,3),opoli,IO,NO) 551 continue c IO, loop over overlaps c correction of polarizability for mutual interaction of groups if(lmut)then do 608 IO=1,NO x=coo(IO,1) y=coo(IO,2) z=coo(IO,3) do 608 IOP=1,NO if(IOP.ne.IO)then xx=coo(IOP,1) yy=coo(IOP,2) zz=coo(IOP,3) rt(1)=x-xx rt(2)=y-yy rt(3)=z-zz rt2=rt(1)**2+rt(2)**2+rt(3)**2 if(rt2.gt.0.0d0)then if(CUTOFF.le.0.0d0.or.rt2.lt.cut2)then c c T-tensor call mktt(rt,tt) c c TTT-tensor call mkttt(rt,ttt) c c aTa to a: call ctd(IO,opolr,temr,NO) call ctd(IO,opoli,temi,NO) call ct(tar,tt,temr) call ct(tai,tt,temi) call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(tcr1,temr,tar) call ct(tcr2,temi,tai) call ct(tci1,temi,tar) call ct(tci2,temr,tai) do 612 ix=1,3 do 612 jx=1,3 bigpolvi(ix,jx)=bigpolvi(ix,jx)+tci1(ix,jx)+tci2(ix,jx) bigpolvr(ix,jx)=bigpolvr(ix,jx)+tcr1(ix,jx)-tcr2(ix,jx) bigpoli(ix,jx)=bigpoli(ix,jx)+tci1(ix,jx)+tci2(ix,jx) 612 bigpolr(ix,jx)=bigpolr(ix,jx)+tcr1(ix,jx)-tcr2(ix,jx) c ATa-aTA to a: do 614 i=1,3 do 614 j=1,3 tcr(i,j)=0.0d0 tci(i,j)=0.0d0 do 614 ix=1,3 do 614 jx=1,3 do 614 kx=1,3 tcr(i,j)=tcr(i,j)+ 1 (aopolr(IO ,ix,jx,i)*ttt(ix,jx,kx)*opolr(IOP,kx,j)- 1 aopolr(IOP,ix,jx,j)*ttt(ix,jx,kx)*opolr(IO ,kx,i))/3.0d0+ 1 (-aopoli(IO ,ix,jx,i)*ttt(ix,jx,kx)*opoli(IOP,kx,j)+ 1 aopoli(IOP,ix,jx,j)*ttt(ix,jx,kx)*opoli(IO ,kx,i))/3.0d0 614 tci(i,j)=tci(i,j)+ 1 (aopolr(IO ,ix,jx,i)*ttt(ix,jx,kx)*opoli(IOP,kx,j)- 1 aopolr(IOP,ix,jx,j)*ttt(ix,jx,kx)*opoli(IO ,kx,i))/3.0d0+ 1 (aopoli(IO ,ix,jx,i)*ttt(ix,jx,kx)*opolr(IOP,kx,j)- 1 aopoli(IOP,ix,jx,j)*ttt(ix,jx,kx)*opolr(IO ,kx,i))/3.0d0 call tadd(bigpolr,tcr) call tadd(bigpoli,tci) call tadd(bigpolvr,tcr) call tadd(bigpolvi,tci) c c GTG to alpha: pi=4.0d0*atan(1.0d0) c sf=(2.0d0*pi/(532.0d0*10.0d0/bohr))**2 sf=(2.0d0*pi/(wnm *10.0d0/bohr))**2 do 616 i=1,3 do 616 j=1,3 tcr(i,j)=0.0d0 tci(i,j)=0.0d0 do 616 ix=1,3 do 616 jx=1,3 tci(i,j)=tci(i,j) 1 +gopoli(IO ,i,jx)*tt(ix,jx)*gopolr(IOP,j,jx)*sf 1 +gopolr(IO ,i,jx)*tt(ix,jx)*gopoli(IOP,j,jx)*sf 616 tcr(i,j)=tcr(i,j) 1 +gopolr(IO ,i,jx)*tt(ix,jx)*gopolr(IOP,j,jx)*sf 1 -gopoli(IO ,i,jx)*tt(ix,jx)*gopoli(IOP,j,jx)*sf call tadd(bigpolvr,tcr) call tadd(bigpolvi,tci) call tadd(bigpolr,tcr) call tadd(bigpoli,tci) c c contribution of alpha...alpha coupling to G: call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(tajr,tt,temr) call ct(taji,tt,temi) call ctd(IO,opolr,temr,NO) call ctd(IO,opoli,temi,NO) call ct(atar1,temr,tajr) call ct(atar2,temi,taji) call ct(atai1,temr,taji) call ct(atai2,temi,tajr) do 620 ix=1,3 call xvp(tci(ix,1),tci(ix,2),tci(ix,3),-coo(IO,1)/2.0d0, 1 -coo(IO,2)/2.0d0,-coo(IO,3)/2.0d0,atai1(1,ix)+atai2(1,ix), 1 atai1(2,ix)+atai2(2,ix),atai1(3,ix)+atai2(3,ix)) 620 call xvp(tcr(ix,1),tcr(ix,2),tcr(ix,3),-coo(IO,1)/2.0d0, 1 -coo(IO,2)/2.0d0,-coo(IO,3)/2.0d0,atar1(1,ix)-atar2(1,ix), 1 atar1(2,ix)-atar2(2,ix),atar1(3,ix)-atar2(3,ix)) call tadd(biggpolr,tcr) call tadd(biggpoli,tci) c c contribution of alpha...A coupling to G do 624 ix=1,3 do 624 jx=1,3 taii(ix,jx)=0.0d0 taji(ix,jx)=0.0d0 tair(ix,jx)=0.0d0 tajr(ix,jx)=0.0d0 do 624 i=1,3 do 624 j=1,3 taii(ix,jx)=taii(ix,jx)+ttt(ix,i,j)*aopoli(IO ,i,j,jx) taji(ix,jx)=taji(ix,jx)+ttt(ix,i,j)*aopoli(IOP,i,j,jx) tair(ix,jx)=tair(ix,jx)+ttt(ix,i,j)*aopolr(IO ,i,j,jx) 624 tajr(ix,jx)=tajr(ix,jx)+ttt(ix,i,j)*aopolr(IOP,i,j,jx) do 625 ix=1,3 do 625 jx=1,3 ataii(ix,jx)=0.0d0 ataji(ix,jx)=0.0d0 atair(ix,jx)=0.0d0 atajr(ix,jx)=0.0d0 do 625 i=1,3 ataii(ix,jx)=ataii(ix,jx)+opolr(IOP,ix,i)*taii(i,jx) 1 +opoli(IOP,ix,i)*tair(i,jx) ataji(ix,jx)=ataji(ix,jx)+opolr(IO ,ix,i)*taji(i,jx) 1 +opoli(IO ,ix,i)*tajr(i,jx) atair(ix,jx)=atair(ix,jx)+opolr(IOP,ix,i)*tair(i,jx) 1 -opoli(IOP,ix,i)*taii(i,jx) 625 atajr(ix,jx)=atajr(ix,jx)+opolr(IO ,ix,i)*tajr(i,jx) 1 -opoli(IO ,ix,i)*taji(i,jx) do 622 ix=1,3 call xvp(tcr(ix,1),tcr(ix,2),tcr(ix,3),-coo(IO,1)/6.0d0, 1 -coo(IO,2)/6.0d0,-coo(IO,3)/6.0d0,atair(ix,1)-atajr(1,ix), 1 atair(ix,2)-atajr(2,ix),atair(ix,3)-atajr(3,ix)) 622 call xvp(tci(ix,1),tci(ix,2),tci(ix,3),-coo(IO,1)/6.0d0, 1 -coo(IO,2)/6.0d0,-coo(IO,3)/6.0d0,ataii(ix,1)-ataji(1,ix), 1 ataii(ix,2)-ataji(2,ix),ataii(ix,3)-ataji(3,ix)) call tadd(biggpoli,tci) call tadd(biggpolr,tcr) c c contribution of Alpha..G to G: call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(taji,tt,temi) call ct(tajr,tt,temr) do 623 ix=1,3 do 623 jx=1,3 tci(jx,ix)=-( 1 gopolr(IO,jx,1)*taji(1,ix)+gopoli(IO,jx,1)*tajr(1,ix)+ 1 gopolr(IO,jx,2)*taji(2,ix)+gopolr(IO,jx,3)*taji(3,ix)+ 1 gopoli(IO,jx,2)*tajr(2,ix)+gopoli(IO,jx,3)*tajr(3,ix)) 623 tcr(jx,ix)=-( 1 gopolr(IO,jx,1)*tajr(1,ix)-gopoli(IO,jx,1)*taji(1,ix)+ 1 gopolr(IO,jx,2)*tajr(2,ix)+gopolr(IO,jx,3)*tajr(3,ix)- 1 gopoli(IO,jx,2)*taji(2,ix)-gopoli(IO,jx,3)*taji(3,ix)) call tadd(biggpolr,tcr) call tadd(biggpoli,tci) c contribution of alpha..alpha to A call ctd(IO,opolr,temr,NO) call ctd(IO,opoli,temi,NO) tvi(1)=coo(IO,1)*temi(1,1)+coo(IO,2)*temi(2,1) 1 +coo(IO,3)*temi(3,1) tvi(2)=coo(IO,1)*temi(1,2)+coo(IO,2)*temi(2,2) 1 +coo(IO,3)*temi(3,2) tvi(3)=coo(IO,1)*temi(1,3)+coo(IO,2)*temi(2,3) 1 +coo(IO,3)*temi(3,3) tvr(1)=coo(IO,1)*temr(1,1)+coo(IO,2)*temr(2,1) 1 +coo(IO,3)*temr(3,1) tvr(2)=coo(IO,1)*temr(1,2)+coo(IO,2)*temr(2,2) 1 +coo(IO,3)*temr(3,2) tvr(3)=coo(IO,1)*temr(1,3)+coo(IO,2)*temr(2,3) 1 +coo(IO,3)*temr(3,3) do 626 ix=1,3 do 626 jx=1,3 do 626 kx=1,3 sui=1.5d0*( 1 (coo(IO,ix)*temi(jx,1)+coo(IO,jx)*temi(ix,1))*tar(1,kx)+ 1 (coo(IO,ix)*temi(jx,2)+coo(IO,jx)*temi(ix,2))*tar(2,kx)+ 1 (coo(IO,ix)*temi(jx,3)+coo(IO,jx)*temi(ix,3))*tar(3,kx)+ 1 (coo(IO,ix)*temr(jx,1)+coo(IO,jx)*temr(ix,1))*tai(1,kx)+ 1 (coo(IO,ix)*temr(jx,2)+coo(IO,jx)*temr(ix,2))*tai(2,kx)+ 1 (coo(IO,ix)*temr(jx,3)+coo(IO,jx)*temr(ix,3))*tai(3,kx)) if(ix.eq.jx)sui=sui 1 -tvr(1)*tai(1,kx)-tvr(2)*tai(2,kx)-tvr(3)*tai(3,kx) 1 -tvi(1)*tar(1,kx)-tvi(2)*tar(2,kx)-tvi(3)*tar(3,kx) sur=1.5d0*( 1 (coo(IO,ix)*temr(jx,1)+coo(IO,jx)*temr(ix,1))*tar(1,kx)+ 1 (coo(IO,ix)*temr(jx,2)+coo(IO,jx)*temr(ix,2))*tar(2,kx)+ 1 (coo(IO,ix)*temr(jx,3)+coo(IO,jx)*temr(ix,3))*tar(3,kx)- 1 (coo(IO,ix)*temi(jx,1)+coo(IO,jx)*temi(ix,1))*tai(1,kx)- 1 (coo(IO,ix)*temi(jx,2)+coo(IO,jx)*temi(ix,2))*tai(2,kx)- 1 (coo(IO,ix)*temi(jx,3)+coo(IO,jx)*temi(ix,3))*tai(3,kx)) if(ix.eq.jx)sur=sur 1 -tvr(1)*tar(1,kx)-tvr(2)*tar(2,kx)-tvr(3)*tar(3,kx) 1 +tvi(1)*tai(1,kx)+tvi(2)*tai(2,kx)+tvi(3)*tai(3,kx) bigapoli(ix,jx,kx)=bigapoli(ix,jx,kx)+sui 626 bigapolr(ix,jx,kx)=bigapolr(ix,jx,kx)+sur c contribution of alpha...A to A do 627 ix=1,3 do 627 jx=1,3 do 627 kx=1,3 sur=aopolr(IO,ix,jx,1)*tajr(1,kx) 1 +aopolr(IO,ix,jx,2)*tajr(2,kx)+aopolr(IO,ix,jx,3)*tajr(3,kx) 1 -aopoli(IO,ix,jx,1)*taji(1,kx) 1 -aopoli(IO,ix,jx,2)*taji(2,kx)-aopoli(IO,ix,jx,3)*taji(3,kx) sui=aopoli(IO,ix,jx,1)*tajr(1,kx) 1 +aopoli(IO,ix,jx,2)*tajr(2,kx)+aopoli(IO,ix,jx,3)*tajr(3,kx) 1 +aopolr(IO,ix,jx,1)*taji(1,kx) 1 +aopolr(IO,ix,jx,2)*taji(2,kx)+aopolr(IO,ix,jx,3)*taji(3,kx) bigapoli(ix,jx,kx)=bigapoli(ix,jx,kx)+sui 627 bigapolr(ix,jx,kx)=bigapolr(ix,jx,kx)+sur c endif c if(cut2.eq.0.or.rt2.lt.cut2)then endif c if(rt2.gt.0.0d0)then endif c IO<>IOP 608 continue endif c record to total polarizability: ii=0 jj=0 do 3 ix=1,3 do 3 jx=1,3 ii=ii+1 alt(nel*(iw-1)+ii )=bigpolr(ix,jx) alt(nel*(iw-1)+ii+ 9)=bigpoli(ix,jx) alt(nel*(iw-1)+ii+90)=bigpolvr(ix,jx) alt(nel*(iw-1)+ii+99)=bigpolvi(ix,jx) alt(nel*(iw-1)+ii+18)=biggpolr(ix,jx) alt(nel*(iw-1)+ii+27)=biggpoli(ix,jx) do 3 kx=1,3 jj=jj+1 alt(nel*(iw-1)+jj+36)=bigapolr(jx,kx,ix) 3 alt(nel*(iw-1)+jj+63)=bigapoli(jx,kx,ix) 1 continue c iw-frequency RETURN END SUBROUTINE IMPROVE(FFB,FFS,RB,RS,NAT,NATS,NS,IND,IB,IS,AM,INDBIG, 1NO,IBONDS,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,IWG,LROA,LRAM, 2ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT,LNAMES, 3LHALTONERROR,CUTOFF,n1,n2,okind,NATS_max,lpol, 1lpolg,lpola,lq,ND,LDATOM,IOPOL,GG,lgrad,sgrad,fgrad,alphf, 1 lnmr,lspin,NMRNS,SPINNS,nmrb,spinb,LQUAD,QS,QB) implicit none INTEGER*4 MX,MXB,IOX PARAMETER (MX=4000,MXB=140,IOX=1500) character*80 NAME(IOX) INTEGER*4 okind,IND(MX),IS(MX,MXB),IB(MX,MXB),INDBIG(IOX,MX), 4IBONDS(4,MX),NBOB(MX),IOL(IOX),ISUM,I,IA,IAS,IEND,IERR,II, 1IIND,IIOII,IIX,IL,INDI,INDIB,INDIS,INDJ,INDS,INDJS,IO,NATS_max, 2IOIJ,IP,IPASS,IPO,ISAT,ISTART,IWG,IX,IXX,J,JA,JJ,IERC,IPOOLD, 3JS,JSTART,JX,K,KK,N1,N2,NAT,NATIO,NATS,NF,NO,NS,INDJB,JJX, 4ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,MX),nats_o(IOX), 5ND(*),ITU,IIT,ido,iao,ic,iop,noo,nu,bodyn(IOX),IOPOL, 1bodylist(IOX,MXB),endn(2,2,IOX),endlist(2,2,IOX,MXB),i1,i2,iap, 1neighbor(IOX,2),nt,jend,nn,kx,ID,IG,ie,natname real*8 FFB(3*NAT,3*NAT),RB(3,MX),RS(3,NATS),x,y,z,att(3,3), 1T3(3),TB(3),AM(MX,3,3),B(3,3),TIJ(3),OTIJ(3),PB(MX,3,3), 2AB(MX,3,3),PS(MX,3,3),AS(MX,3,3),WF(IOX),ALPHA(3*MX,3,3), 3ALPHAS(3*MX,3,3),AA(3*MX,3,3,3),AAS(3*MX,3,3,3),G(3*MX,3,3), 4GS(3*MX,3,3),DISTM,SUM,FNEW,PNEW,ANEW,ALPHANEW,GNEW,AANEW, 5DIST,FOLD,POLD,AOLD,AD1,AD2,AD0,CTF,CUTOFF,ata(3,3),ta(3,3), 6RM(IOX,3,3),atad(3,3),sfac,FFS(*),A,XS,XB,TOL,RMS,tt(3,3),rt(3), 7d1,d2,dd,ttt(3,3,3),eps,aname,pname,aaname,gname,alphaname, 1opol(IOX,3,3),coo(IOX,3),tc(3,3),sf,rt2,alphf, 2bohr,gpol,apol,apolname,gpolname,aopol(IOX,3,3,3),alder(3,3), 3gopol(IOX,3,3),alst(3,3,3),glst(3,3,3),tem(3,3),al(3,3), 4aast(3,3,3,3),bigpol(3,3),pi,biggpol(3,3),bigapol(3,3,3), 1tai(3,3),taj(3,3),atai(3,3),ataj(3,3),tv(3),su,zsum(IOX),at(3,3), 1al_l(3,3),ts(3,3),atta(3,3),rname,polname,pol,qname, 1ul(3,3),uj(3,3),ult(3,3),ujt(3,3),ultaj(3,3),ujtad(3,3), 1rl(3),vl(3,3),vlt(3,3),vltaj(3,3),ujts(3,3),rj(3), 1ujtsal(3,3),adt(3,3),adta(3,3),gradname,GG(*),GOLD,xj,yj,zj, 1xk,yk,zk,QS(MX,3,3,3),QB(MX,3,3,3),QNEW,QOLD,AD3 LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LNAMES,LHALTONERROR,lpol,li1,li2,lpolg,lpola,lq,LDATOM,sgrad, 1lgrad,fgrad,LQUAD COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS common/cname/aname(IOX,MX,3,3),polname(IOX,3,3),pol(3,3), 1pname(IOX,MX,3,3),aaname(IOX,3*MX,3,3,3),gname(IOX,3*MX,3,3), 2alphaname(IOX,3*MX,3,3),rname(IOX,3,MX),qname(IOX,MX,3,3,3), 3gpol(3,3),apol(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3), 2gradname(IOX,3*MX),natname(IOX),NAME integer*4 KA,NW,IPOE,JAS real*8 ggp(3) real*8,allocatable::gij(:,:,:),sjk(:,:),sijk(:,:),gsv(:) c NMR variables: logical lnmr,lspin real*8 NMRNS(*),SPINNS(*),nmrb(*),spinb(NAT,NAT) c bohr=0.529177D0 IF(NO.EQ.0)GOTO 19 IF(LDIA)WRITE(3,*)' Diagonal force field will be changed' IF(LOFF)WRITE(3,*)' Off-diagonals will be changed' IF(LABS)WRITE(3,*)' Eletric derivatives will be transformed' IF(LVCD)WRITE(3,*)' Magnetic derivatives will be transformed' IF(LQUAD)WRITE(3,*)' Quadrupole derivatives will be transformed' IF(LROA)WRITE(3,*)' ROA parameters will be transformed' IF(LRAM)WRITE(3,*)' Raman parameters will be transformed' IF(LSTRICT)WRITE(3,*)' Only transformed atoms in the beds' IF(IWG.EQ.1)WRITE(3,*)' Multiple overlaps will be weighted' IF(LINV)WRITE(3,*)' Smaller fragment as optical antipode' IF(LSELECT)WRITE(3,*)' FC for non-transformed atoms set to 0 !' IF(LNAMES)WRITE(3,*)' Named small fragments' IF(LHALTONERROR)WRITE(3,*)' Will stop on error' IF(.not.LHALTONERROR)WRITE(3,*)' Will not stop on error' if(okind.eq.0)write(3,*)'selection on centers of masses' if(okind.eq.1)write(3,*)'selection on RMS distances' if(okind.eq.2)write(3,*)'selection on molecular RMS distances' if(IWG.eq.0)write(3,*)'IWG=0 best overlap only' if(IWG.eq.1)write(3,*)'IWG=1 weighing of overlaps' if(lgrad)write(3,*)'transfer of gradient requested' if(lnmr)write(3,*)'transfer of NMR shift requested' if(lspin)write(3,*)'transfer of NMR J-constants requested' if(sgrad)write(3,*)'sum of all gradients' if(fgrad)write(3,*)'gradient with force field' if(LDATOM)write(3,*)'In mutual polarizability correction,', 1'use distances from differentiated atom' WRITE(3,3999)CUTOFF 3999 format('Distance cutoff ',f12.2,' A') ctf=CUTOFF**2 IOIJ=0 IPOOLD=0 c c group polarizabilities: if(lpol)then c c total polarizability do 607 ix=1,3 do 607 jx=1,3 biggpol(ix,jx)=0.0d0 bigpol(ix,jx)=0.0d0 do 607 kx=1,3 607 bigapol(ix,jx,kx)=0.0d0 c loop over overlaps: DO 551 IO=1,NO c make the transformation matrix for whole overlap: call CRMSo(NAT,INDBIG,IOX,MX,IO,LNAMES,natname,RS,RB,IND, 2 rname,NATS,LINV) c c calculate overlap center: do 53 ix=1,3 53 coo(IO,ix)=0.0d0 noo=0 nu=0 zsum(IO)=0.0d0 do 52 iao=1,NAT if(INDBIG(IO,iao).GT.0)then if(lq)then su=dble(ND(iao)) else su=1.0d0 endif zsum(IO)=zsum(IO)+su noo=noo+1 do 54 ix=1,3 54 coo(IO,ix)=coo(IO,ix)+RB(ix,iao)*su c c look at other overlaps, and find out if this atoms is part of them: ic=0 do 552 IOP=1,NO 552 if(IOP.ne.IO.and.INDBIG(IOP,iao).ne.0)ic=ic+1 if(ic.eq.0)nu=nu+1 endif 52 continue if(noo.gt.0)then do 55 ix=1,3 55 coo(IO,ix)=coo(IO,ix)/zsum(IO)/bohr endif write(3,6009)IO,noo,NATS,nu 6009 format('Overlap',i5,':',i5,' atoms transferred of',i5,',', 1 i5,' unique') c c if named fragments, fill alpha, G, A for each overlap: if(lnames)then call ctd(IO,polname,pol,IOX) if(lpolg)call ctd(IO,gpolname,gpol,IOX) if(lpola)call ctdd(IO,apolname,apol,IOX) endif c c rotate polarizability: if(lpol)call rpp(IOX,IO,opol,pol,A) c c rotate G if(lpolg)call rpp(IOX,IO,gopol,gpol,A) c c rotate A if(lpola)call rppp(IOX,IO,aopol,apol,A) c c construct alpha,G,A of all system: do 5519 IX=1,3 do 5519 JX=1,3 c c plain sum: bigpol(IX,JX) =bigpol(IX,JX) +opol(IO,IX,JX) biggpol(IX,JX)=biggpol(IX,JX)+gopol(IO,IX,JX) do 55191 KX=1,3 55191 bigapol(IX,JX,KX)=bigapol(IX,JX,KX)+aopol(IO,IX,JX,KX) c c polarizability part for origin dependence of G: DO 618 IG=1,3 DO 618 ID=1,3 618 biggpol(IX,JX)=biggpol(IX,JX) 1 -0.5d0*EPS(JX,IG,ID)*coo(IO,IG)*opol(IO,IX,ID) c c polarizability part for origin dependence of A: do 5519 KX=1,3 if(JX.eq.IX)bigapol(IX,JX,KX)=bigapol(IX,JX,KX) 1 -coo(IO,1)*opol(IO,KX,1) 2 -coo(IO,2)*opol(IO,KX,2) 3 -coo(IO,3)*opol(IO,KX,3) 5519 bigapol(IX,JX,KX)=bigapol(IX,JX,KX) 1 +1.5d0*(coo(IO,JX)*opol(IO,KX,IX)+coo(IO,IX)*opol(IO,KX,JX)) 551 continue c correction of polarizability for mutual interaction of groups do 608 IO=1,NO do 608 IOP=1,NO if(IOP.ne.IO)then do 609 ix=1,3 609 rt(ix)=coo(IO,ix)-coo(IOP,ix) rt2=rt(1)**2+rt(2)**2+rt(3)**2 if(rt2.gt.0.0d0)then c c T-tensor call mktt(rt,tt) c c TTT-tensor call mkttt(rt,ttt) c c aTa to a: call ctd(IO,opol,tem,IOX) call ct(ta,tt,tem) call ctd(IOP,opol,tem,IOX) call ct(tc,tem,ta) do 612 ix=1,3 do 612 jx=1,3 612 bigpol(ix,jx)=bigpol(ix,jx)+tc(ix,jx) c ATa-aTA to a: if(lpola)then do 614 i=1,3 do 614 j=1,3 tc(i,j)=0.0d0 do 614 ix=1,3 do 614 jx=1,3 do 614 kx=1,3 614 tc(i,j)=tc(i,j)+ 1 (aopol(IO ,ix,jx,i)*ttt(ix,jx,kx)*opol(IOP,kx,j)- 1 aopol(IOP,ix,jx,j)*ttt(ix,jx,kx)*opol(IO ,kx,i))/3.0d0 call tadd(bigpol,tc) endif c c GTG to alpha: if(lpolg)then pi=4.0d0*atan(1.0d0) sf=(2.0d0*pi/(532.0d0*10.0d0/bohr))**2 do 616 i=1,3 do 616 j=1,3 tc(i,j)=0.0d0 do 616 ix=1,3 do 616 jx=1,3 616 tc(i,j)=tc(i,j)+ 1 gopol(IO ,i,jx)*tt(ix,jx)*gopol(IOP,j,jx)*sf call tadd(bigpol,tc) endif c c contribution of alpha...alpha coupling to G: call ctd(IOP,opol,tem,IOX) call ct(taj,tt,tem) call ctd(IO,opol,tem,IOX) call ct(ata,tem,taj) do 620 ix=1,3 do 620 jx=1,3 tc(ix,jx)=0.0d0 do 620 i=1,3 do 620 j=1,3 620 tc(ix,jx)=tc(ix,jx)-EPS(jx,i,j)/2.0d0*coo(IO,i)*ata(j,ix) call tadd(biggpol,tc) c c contribution of alpha...A coupling to G if(lpola)then do 624 ix=1,3 do 624 jx=1,3 tai(ix,jx)=0.0d0 taj(ix,jx)=0.0d0 do 624 i=1,3 do 624 j=1,3 tai(ix,jx)=tai(ix,jx)+ttt(ix,i,j)*aopol(IO ,i,j,jx) 624 taj(ix,jx)=taj(ix,jx)+ttt(ix,i,j)*aopol(IOP,i,j,jx) do 625 ix=1,3 do 625 jx=1,3 atai(ix,jx)=0.0d0 ataj(ix,jx)=0.0d0 do 625 i=1,3 atai(ix,jx)=atai(ix,jx)+opol(IOP,ix,i)*tai(i,jx) 625 ataj(ix,jx)=ataj(ix,jx)+opol(IO ,ix,i)*taj(i,jx) do 622 ix=1,3 do 622 jx=1,3 tc(ix,jx)=0.0d0 do 622 ie=1,3 do 622 id=1,3 622 tc(ix,jx)=tc(ix,jx)- 1 EPS(jx,ie,id)*coo(IO,ie)*(atai(ix,id)-ataj(id,ix))/6.0d0 call tadd(biggpol,tc) endif c c contribution of Alpha..G to G: call ctd(IOP,opol,tem,IOX) call ct(taj,tt,tem) do 623 ix=1,3 do 623 jx=1,3 623 tc(jx,ix)=-(gopol(IO,jx,1)*taj(1,ix)+ 1 gopol(IO,jx,2)*taj(2,ix)+gopol(IO,jx,3)*taj(3,ix)) call tadd(biggpol,tc) c contribution of alpha..alpha to A if(lpola)then call ctd(IOP,opol,tem,IOX) call ct(taj,tt,tem) call ctd(IO,opol,tem,IOX) tv(1)=coo(IO,1)*tem(1,1)+coo(IO,2)*tem(2,1)+coo(IO,3)*tem(3,1) tv(2)=coo(IO,1)*tem(1,2)+coo(IO,2)*tem(2,2)+coo(IO,3)*tem(3,2) tv(3)=coo(IO,1)*tem(1,3)+coo(IO,2)*tem(2,3)+coo(IO,3)*tem(3,3) do 626 ix=1,3 do 626 jx=1,3 do 626 kx=1,3 su=1.5d0*( 1 (coo(IO,ix)*tem(jx,1)+coo(IO,jx)*tem(ix,1))*taj(1,kx)+ 1 (coo(IO,ix)*tem(jx,2)+coo(IO,jx)*tem(ix,2))*taj(2,kx)+ 1 (coo(IO,ix)*tem(jx,3)+coo(IO,jx)*tem(ix,3))*taj(3,kx)) if(ix.eq.jx) 1 su=su-tv(1)*taj(1,kx)-tv(2)*taj(2,kx)-tv(3)*taj(3,kx) 626 bigapol(ix,jx,kx)=bigapol(ix,jx,kx)+su endif c contribution of alpha...A to A c A_iAab Tijae a_jeg if(lpola)then call ctd(IOP,opol,tem,IOX) call ct(taj,tt,tem) do 627 ix=1,3 do 627 jx=1,3 do 627 kx=1,3 su=aopol(IO,ix,jx,1)*taj(1,kx) 1 + aopol(IO,ix,jx,2)*taj(2,kx)+aopol(IO,ix,jx,3)*taj(3,kx) 627 bigapol(ix,jx,kx)=bigapol(ix,jx,kx)+su endif c endif endif 608 continue call WRITEPOL(bigpol,'POL.TTT') if(lpolg)call WRITEPOL(biggpol,'POL.G.TTT') if(lpola)call WRITEPOLA(bigapol) call w36('Static polarizabilities rotated for each overlap') if(lpolg)call w36('Gs rotated for each overlap') if(lpola)call w36('A for each overlap') if(LWR)then write(3,6301) 6301 format('Isotropic Polarizabilities:') do 555 IO=1,NO 555 write(3,6302)IO,(opol(IO,1,1)+opol(IO,2,2)+opol(IO,3,3))/3.0d0 6302 format(i3,f10.3) write(3,6303)NO 6303 format('Overlap centers:',/,I4) do 556 IO=1,NO 556 write(3,6304)(coo(IO,ix)*bohr,ix=1,3) 6304 format('6 ',3f12.4) endif c c sort atoms in each overlap: do 557 IO=1,NO c c two closest overlaps i1=0 i2=0 d1=9.0d9 d2=9.0d9 do 561 IOP=1,NO if(IOP.ne.IO)then dd=(coo(IO,1)-coo(IOP,1))**2+ 1 (coo(IO,2)-coo(IOP,2))**2+(coo(IO,3)-coo(IOP,3))**2 if(dd.le.d2.and.dd.gt.d1)then d2=dd i2=IOP else if(dd.le.d2.and.dd.le.d1)then d2=d1 i2=i1 d1=dd i1=IOP endif endif endif 561 continue d1=dsqrt(d1)*bohr d2=dsqrt(d2)*bohr c c neighbor= c closests: 111111111111111111111 222222222222 c overlap: XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX c body BBBBBBBBBBBBBBBBBBBBBB c end11 AAA c end12 BBBBBB c end21 AAA c end22 BBBB c body, atoms that are only in IO: bodyn(IO)=0 c end1, end2, atoms that are also part of the closest and second closest c 11-closer to IO, 12- closer to 1, etc.: endn(1,1,IO)=0 endn(1,2,IO)=0 endn(2,1,IO)=0 endn(2,2,IO)=0 neighbor(IO,1)=0 neighbor(IO,2)=0 do 558 ia=1,NAT li1=.false. li2=.false. x=RB(1,ia)/bohr y=RB(2,ia)/bohr z=RB(3,ia)/bohr if(INDBIG(IO,ia).gt.0)then do 559 IOP=1,NO if(IOP.ne.IO.and.(IOP.eq.i1.or.IOP.eq.i2))then do 560 iap=1,NAT if(INDBIG(IOP,iap).gt.0.and.iap.eq.ia)then if(IOP.eq.i1)li1=.true. if(IOP.eq.i2)li2=.true. endif 560 continue endif 559 continue if(.not.li1.and..not.li2)then c (ia is only in IO) if(bodyn(IO).lt.MXB)then bodyn(IO)=bodyn(IO)+1 bodylist(IO,bodyn(IO))=ia else call w36('Warning: bodyn>MXB') endif else c (ia is also in IOP) if(li1.and.li2)then c find if closer to i1 or i2: d1=(x-coo(i1,1))**2+(y-coo(i1,2))**2+(z-coo(i1,3))**2 d2=(x-coo(i2,1))**2+(y-coo(i2,2))**2+(z-coo(i2,3))**2 if(d1.lt.d2)then li2=.false. else li1=.false. endif endif if(li1)then IOP=i1 ii=1 else IOP=i2 ii=2 endif c find if closer to IO or IOP: d1=(x-coo(IO ,1))**2+(y-coo(IO ,2))**2+(z-coo(IO ,3))**2 d2=(x-coo(IOP,1))**2+(y-coo(IOP,2))**2+(z-coo(IOP,3))**2 if(d1.lt.d2)then ic=1 else ic=2 endif if(neighbor(IO,ii).eq.0)neighbor(IO,ii)=IOP if(endn(ii,ic,IO).lt.MXB)then endn(ii,ic,IO)=endn(ii,ic,IO)+1 endlist(ii,ic,IO,endn(ii,ic,IO))=ia else call w36('Warning: end1n>MXB') endif endif endif 558 continue if(LWR)then write(3,*) write(3,*)'Overlap :',IO write(3,*)'Neighbors:',neighbor(IO,1),neighbor(IO,2) write(3,*)'Body atoms:' write(3,30108)(bodylist(IO,ii),ii=1,bodyn(IO)) do 562 IOP=1,2 write(3,3112)IOP,(endn(IOP,ic,IO),ic=1,2) 3112 format('Terminal',i2,':',i4,' and ',i4,' atoms:') do 562 ic=1,2 562 write(3,30108)(endlist(IOP,ic,IO,ii),ii=1,endn(IOP,ic,IO)) 30108 format(10i8) endif 557 continue endif c c if fgrad then decompose gradients on each atom to c projections to others, gi = (j<>i) gi_j * si_j c si_j= r_j - r_i if(fgrad)then if(LNAMES)then IPOE=NO else IPOE=1 endif c gij: c IPOE ... overlap index c NATS_max .. atom in the fragment c 3 ... gradient components c NATS_max ... to these atoms in the fragment allocate(gij(IPOE,NATS_max,NATS_max)) do 651 IPO=1,IPOE NW=NATNAME(IPO)-1 allocate(sjk(NW,NW),sijk(NW,NW),gsv(NW)) do 65 IA=1,NATNAME(IPO) x=rname(IPO,1,IA) y=rname(IPO,2,IA) z=rname(IPO,3,IA) JJ=0 do 66 JA=1,NATNAME(IPO) if(JA.ne.IA)then JJ=JJ+1 xj=rname(IPO,1,JA)-x yj=rname(IPO,2,JA)-y zj=rname(IPO,3,JA)-z call nv7(xj,yj,zj) gsv(JJ)=gradname(IPO,1+3*(IA-1))*xj 1 + gradname(IPO,2+3*(IA-1))*yj 1 + gradname(IPO,3+3*(IA-1))*zj KK=0 do 67 KA=1,NATNAME(IPO) if(KA.ne.IA)then KK=KK+1 xk=rname(IPO,1,KA)-x yk=rname(IPO,2,KA)-y zk=rname(IPO,3,KA)-z call nv7(xk,yk,zk) sjk(JJ,KK)=xj*xk+yj*yk+zj*zk endif 67 continue endif 66 continue do 71 JJ=1,NW 71 sjk(JJ,JJ)=sjk(JJ,JJ)+alphf call INV(sjk,sijk,NW,1.0d-25) JJ=0 do 70 JA=1,NATNAME(IPO) gij(IPO,IA,JA)=0.0d0 if(JA.NE.IA)then JJ=JJ+1 do 68 KK=1,NW 68 gij(IPO,IA,JA)=gij(IPO,IA,JA)+sijk(JJ,KK)*gsv(KK) endif 70 continue 65 continue 651 deallocate(sjk,sijk,gsv) write(3,*)'FGRAD: ',IPOE,' gradients decomposed' write(6,*)'FGRAD: ',IPOE,' gradients decomposed' do 69 KK=1,3*NAT 69 GG(KK)=0.0D0 endif C C Look at all atom pairs i-j: istart=1 iend=NAT if(n1.gt.0)then istart=n1 iend=n2 write(3,*)'first index ',n1,' - ',n2 endif write(6,*)NAT,' atoms:' DO 20 IA= istart,iend write(6,6001)IA 6001 format(i4,$) if(mod(IA,15).eq.0)write(6,*) write(3,*)IA,'. atom of ',NAT DO 20 JA=IA,NAT write(3,*) write(3,3008)IA,JA 3008 format(' Pair ',2I6) IF(.NOT.LOFF.AND..NOT.lspin.AND.IA.NE.JA)GOTO 20 c c look at distance, skip for a large one if(CUTOFF.gt.0.0d0)THEN dist=(RB(1,IA)-RB(1,JA))**2+(RB(2,IA)-RB(2,JA))**2+ 1 (RB(3,IA)-RB(3,JA))**2 if(dist.gt.ctf)then if(LWR)write(3,3998)sqrt(dist) 3998 format('Skipped because of the distance,',f20.2,' A') goto 20 endif endif C C Look at all overlaps io, find that which fits best the i-j pair: c c for okind=2, look at the overlaps only once c cccccccccccccccccccccccccccccccccccccccccccccccccccc IF(okind.ne.2.or.(IA.eq.istart.and.JA.EQ.istart))THEN C C Calculate center of i/j radiusvectors: DO 22 IX=1,3 22 TIJ(IX)=0.5d0*(RB(IX,IA)+RB(IX,JA)) DISTM=1.0d5 NF=0 DO 21 IO=1,NO C C Check, if io contains ij: IF(INDBIG(IO,IA).GT.0.AND.INDBIG(IO,JA).GT.0)THEN NF=NF+1 IOL(NF)=IO C C Calculate center of the overlap and distance from i-j center c i.........TIJ.............j c ----------OTIJ--------------------- DO 23 IX=1,3 23 OTIJ(IX)=0.0D0 NATIO=0 DO 24 II=1,NAT IIOII=INDBIG(IO,II) C If Lstrict, take only transformed atoms IF(LSTRICT.AND.IIOII.LT.0)IIOII=0 IIOII=ABS(IIOII) IF(IIOII.NE.0)THEN NATIO=NATIO+1 DO 25 IX=1,3 25 OTIJ(IX)=OTIJ(IX)+RB(IX,II) ENDIF 24 CONTINUE c c If required, calculate RMS overlap and use this criterium instead c of the distances,remember rotation matrix (RM) if(okind.gt.0)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IO, 1 IA,JA,NBOB,okind,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) do 155 IX=1,3 do 155 JX=1,3 155 RM(NF,IX,JX)=A(IX,JX) ipass_o(NF)=IPASS WF(NF)=DSQRT(RMS) ierc_o(NF)=IERC itu_o(NF)=ITU do 58 I=1,ITU 58 ind_o(NF,I)=IND(I) do 59 I=1,3 do 59 J=1,3 59 RM(NF,I,J)=A(I,J) nats_o(NF)=NATS IF(RMS.LT.DISTM)THEN DISTM=RMS IOIJ=IO ENDIF else c Use distances: DO 26 IX=1,3 26 OTIJ(IX)=OTIJ(IX)/DBLE(NATIO) DIST=(OTIJ(1)-TIJ(1))**2+(OTIJ(2)-TIJ(2))**2+(OTIJ(3)-TIJ(3))**2 IF(DIST.LT.DISTM)THEN DISTM=DIST IOIJ=IO ENDIF WF(NF)=DSQRT(DIST) endif ENDIF 21 CONTINUE c if(NF.EQ.0)then WRITE(3,3003)0,0 else WRITE(3,3003)NF,IOIJ 3003 FORMAT(I5,' overlaps found; best is # ',I4) endif C C Calculate the weighting factors: SUM=0.0d0 DO 43 I=1,NF 43 IF(WF(I).GT.1.0d-9)SUM=SUM+1.0d0/WF(I) DO 42 I=1,NF IF(WF(I).LE.1.0d-9)THEN WF(1)=1.0d0 NF=1 GOTO 44 ENDIF WF(I)=1.0d0/WF(I)/SUM 42 CONTINUE 44 CONTINUE C c cccccccccccccccccccccccccccccccccccccccccccccccccc ENDIF c c if sgrad, construct sum of gradients from all fragments c SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD if(IA.EQ.JA.and.sgrad.and.lgrad)then C Loop over possible overlaps: DO 4500 IP=1,NF IPO=IOL(IP) C C Find the atomic groups to be placed on together for IA-JA: c and the rotation matrix if(okind.eq.0)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1 IA,JA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) else IPASS=ipass_o(IP) IERC=ierc_o(IP) ITU=itu_o(IP) do 5600 I=1,ITU 5600 IND(I)=ind_o(IP,I) do 5700 I=1,3 do 5700 J=1,3 5700 A(I,J)=RM(IP,I,J) NATS=nats_o(IP) endif IF(LWR)WRITE(3,*)IPASS,' passes' if(IERC.ne.0)goto 20 WRITE(3,3004)IPO,1.0d0,ITU,(IND(I),I=1,ITU) WRITE(3,3005)(INDBIG(IPO,IND(I)),I=1,ITU) if(LNAMES)write(3,3399)NAME(IPO) IF(LINV)THEN DO 500 I=1,3 DO 500 J=1,3 500 A(I,J)=-A(I,J) ENDIF IAS=INDBIG(IPO,IA) DO 6100 IX=1,3 GOLD=GG(IX+3*(IA-1)) C Zero out at the start of the summation IF(IP.EQ.1)GG(IX+3*(IA-1))=0.0d0 GNEW=0.0d0 DO 6000 II=1,3 if(LNAMES)then GNEW=GNEW+gradname(IPO,II+3*(IAS-1))*A(IX,II) else GNEW=GNEW+gradname( 1,II+3*(IAS-1))*A(IX,II) endif 6000 continue GG(IX+3*(IA-1))=GNEW+GG(IX+3*(IA-1)) 6100 IF(LWR)WRITE(3,3016)IA,IX,GOLD,GNEW,WF(IP) 3016 FORMAT(' grad ',2I3,F11.6,'/',F11.6,' (',F7.4,')') 4500 CONTINUE endif c SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD SGRAD c C Loop over possible overlaps: ISUM=0 c IPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIP DO 45 IP=1,NF IPO=IOL(IP) IF(IWG.EQ.0)THEN IF(IPO.NE.IOIJ)GOTO 45 WF(IP)=1.0d0 ENDIF ISUM=ISUM+1 C C Find the atomic groups to be placed on together for IA-JA: c and the rotation matrix if(okind.eq.0)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1 IA,JA,NBOB,0,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2 IERC,LINV,LSTRICT,rname,NATS) else IPASS=ipass_o(IP) IERC=ierc_o(IP) ITU=itu_o(IP) do 56 I=1,ITU 56 IND(I)=ind_o(IP,I) do 57 I=1,3 do 57 J=1,3 57 A(I,J)=RM(IP,I,J) NATS=nats_o(IP) endif IF(LWR)WRITE(3,*)IPASS,' passes' if(IERC.ne.0)goto 20 WRITE(3,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU) WRITE(3,3005)(INDBIG(IPO,IND(I)),I=1,ITU) 3004 FORMAT('Overlap ',I2,': weighting factor: ',F7.4,';', 1 I4,' atom bed:',/, 2 ' Big : ',40I3) 3005 FORMAT( ' Small: ',40I3,/) if(LNAMES)write(3,3399)NAME(IPO) 3399 FORMAT(A80) C IF(LINV)THEN DO 50 I=1,3 DO 50 J=1,3 50 A(I,J)=-A(I,J) ENDIF C if(LNAMES)then if(IPOOLD.EQ.0.or.IPOOLD.NE.IPO)then NATS=NATNAME(IPO) IPOOLD=IPO endif ido=IPO ELSE ido=1 ENDIF C Transform the force field: IF((LDIA.AND.IA.EQ.JA).OR.(LOFF.AND.IA.NE.JA))THEN DO 36 IX=1,3 INDIB=3*IA-3+IX JSTART=1 IF(IA.EQ.JA)JSTART=IX DO 36 JX=JSTART,3 INDJB=3*JA-3+JX FOLD=FFB(INDIB,INDJB) C Zero out at the start of the summation IF(ISUM.EQ.1)FFB(INDIB,INDJB)=0.0d0 FNEW=0.0d0 DO 37 IIX=1,3 INDIS=3*INDBIG(IPO,IA)-3+IIX DO 37 JJX=1,3 INDJS=3*INDBIG(IPO,JA)-3+JJX 37 FNEW=FNEW+A(IX,IIX)*A(JX,JJX)* 1 FFS((ido-1)*(3*NATS_max)**2+(INDIS-1)*3*NATS_max+INDJS) IF(LWR)WRITE(3,3001)IA,JA, IX,JX,FOLD,FNEW,IIT,TOL,WF(IP) FFB(INDIB,INDJB)=FNEW*WF(IP)+FFB(INDIB,INDJB) FFB(INDJB,INDIB)=FFB(INDIB,INDJB) 36 CONTINUE ENDIF C C Transform the shifts: IF(lnmr.and.IA.EQ.JA)then FOLD=nmrb(IA) C Zero out at the start of the summation IF(ISUM.EQ.1)nmrb(IA)=0.0d0 IAS=INDBIG(IPO,IA) FNEW=NMRNS((ido-1)*NATS_max+IAS) IF(LWR)WRITE(3,3001)IA,IA, 1,1,FOLD,FNEW,IIT,TOL,WF(IP) nmrb(IA)=FNEW*WF(IP)+nmrb(IA) ENDIF C C Transform the spin-spin constants IF(lspin)then FOLD=spinb(IA,JA) C Zero out at the start of the summation IF(ISUM.EQ.1)spinb(IA,JA)=0.0d0 IAS=INDBIG(IPO,IA) JAS=INDBIG(IPO,JA) FNEW=SPINNS((ido-1)*NATS_max**2+(IAS-1)*NATS_max+JAS) IF(LWR)WRITE(3,3001)IA,JA, 0,0,FOLD,FNEW,IIT,TOL,WF(IP) spinb(IA,JA)=FNEW*WF(IP)+spinb(IA,JA) spinb(JA,IA)=spinb(IA,JA) ENDIF C C C Transform of gradient if(lgrad.and..not.sgrad)then if(fgrad)then c transform gradient together with FF for each pair: c put part of the gradient along the i->j bo IAS=INDBIG(IPO,IA) JAS=INDBIG(IPO,JA) x=rname(IPO,1,JAS)-rname(IPO,1,IAS) y=rname(IPO,2,JAS)-rname(IPO,2,IAS) z=rname(IPO,3,JAS)-rname(IPO,3,IAS) call nv7(x,y,z) ggp(1)=gij(IPO,IAS,JAS)*x ggp(2)=gij(IPO,IAS,JAS)*y ggp(3)=gij(IPO,IAS,JAS)*z c Part of gradient on IA pointing to JA: call trgrp(IA,GG,A,WF,IP,LWR,ggp) ggp(1)=-gij(IPO,JAS,IAS)*x ggp(2)=-gij(IPO,JAS,IAS)*y ggp(3)=-gij(IPO,JAS,IAS)*z c Part of gradient on JA pointing to IA: call trgrp(JA,GG,A,WF,IP,LWR,ggp) else IF(IA.EQ.JA)call trgr(IPO,IA,GG,INDBIG,IOX,MX,A,ISUM,LNAMES, 1 WF,IP,LWR,gradname) endif endif C C Transform the dipole derivatives IF(IA.EQ.JA.AND.(LABS.OR.LVCD))THEN IAS=INDBIG(IPO,IA) DO 40 IX=1,3 DO 40 JX=1,3 POLD=PB(IA,IX,JX) AOLD=AB(IA,IX,JX) C Zero out at the start of the summation IF(ISUM.EQ.1)THEN PB(IA,IX,JX)=0.0d0 AB(IA,IX,JX)=0.0d0 ENDIF PNEW=0.0d0 ANEW=0.0d0 DO 41 II=1,3 AD1=A(IX,II) DO 41 JJ=1,3 AD2=AD1*A(JX,JJ) if(LNAMES)then PNEW=PNEW+pname(IPO,IAS,II,JJ)*AD2 ANEW=ANEW+aname(IPO,IAS,II,JJ)*AD2 else PNEW=PNEW+PS(IAS,II,JJ)*AD2 ANEW=ANEW+AS(IAS,II,JJ)*AD2 endif 41 continue PB(IA,IX,JX)=PNEW*WF(IP)+PB(IA,IX,JX) AB(IA,IX,JX)=ANEW*WF(IP)+AB(IA,IX,JX) IF (.NOT.LVCD)AB(IA,IX,JX)=AOLD IF (.NOT.LVCD)ANEW=AOLD IF (.NOT.LABS)PB(IA,IX,JX)=POLD IF (.NOT.LABS)PNEW=POLD 40 IF(LWR)WRITE(3,3006)IA,IX,JX,POLD,PNEW,AOLD,ANEW,WF(IP) 3006 FORMAT(3I3,F11.6,'/',F11.6,3X,F11.6,'/',F11.6,' (',F7.4,')') if(LQUAD)then C Transform quadrupole derivatives: DO 61 IX=1,3 DO 61 JX=1,3 DO 61 KX=1,3 QOLD=QB(IA,IX,JX,KX) IF(ISUM.EQ.1)QB(IA,IX,JX,KX)=0.0d0 QNEW=0.0d0 DO 62 II=1,3 AD1=A(IX,II) DO 62 JJ=1,3 AD2=AD1*A(JX,JJ) DO 62 KK=1,3 AD3=AD2*A(KX,KK) if(LNAMES)then QNEW=QNEW+qname(IPO,IAS,II,JJ,KK)*AD3 else QNEW=QNEW+QS(IAS,II,JJ,KK)*AD3 endif 62 continue QB(IA,IX,JX,KX)=QNEW*WF(IP)+QB(IA,IX,JX,KX) 61 IF(LWR)WRITE(3,3006)IA,IX,JX,QOLD,QNEW,QOLD,QNEW,WF(IP) endif ENDIF C C Transform the ROA and Raman tensors IF(IA.EQ.JA.AND.(LROA.or.LRAM))THEN if(lwr)then call wrt2(IA,MX,ALPHA,'Old Alpha:') call wrt2(IA,MX,G,'Old G:') call wrt3(IA,MX,AA,'Old A:') endif IAS=INDBIG(IPO,IA) if(LNAMES)then do 599 IXX=1,3 INDS=3*(IAS-1)+IXX do 599 II=1,3 do 599 JJ=1,3 do 5993 KK=1,3 5993 aast(IXX,II,JJ,KK)=aaname(IPO,INDS,II,JJ,KK) alst(IXX,II,JJ)=alphaname(IPO,INDS,II,JJ) 599 glst(IXX,II,JJ)=gname(IPO,INDS,II,JJ) else do 5991 IXX=1,3 INDS=3*(IAS-1)+IXX do 5991 II=1,3 do 5991 JJ=1,3 do 5994 KK=1,3 5994 aast(IXX,II,JJ,KK)=AAS(INDS,II,JJ,KK) alst(IXX,II,JJ)=ALPHAS(INDS,II,JJ) 5991 glst(IXX,II,JJ)=GS(INDS,II,JJ) endif IF(LRAM)THEN DO 46 IX=1,3 IIND=3*(IA-1)+IX DO 46 I=1,3 DO 46 J=1,3 C Zero out at the start of the summation IF(ISUM.EQ.1)ALPHA(IIND,I,J)=0.0d0 ALPHANEW=0.0d0 DO 47 IXX=1,3 AD0=A(IX,IXX) DO 47 II=1,3 AD1=AD0*A(I,II) DO 47 JJ=1,3 47 ALPHANEW=ALPHANEW+alst(IXX,II,JJ)*AD1*A(J,JJ) 46 ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+ALPHANEW*WF(IP) ENDIF IF(LROA)THEN c Rotation of G: DO 461 IX=1,3 IIND=3*(IA-1)+IX DO 461 I=1,3 DO 461 J=1,3 C Zero out at the start of the summation IF(ISUM.EQ.1)G(IIND,I,J)=0.0d0 GNEW=0.0d0 DO 471 IXX=1,3 AD0=A(IX,IXX) DO 471 II=1,3 AD1=AD0*A(I,II) DO 471 JJ=1,3 471 GNEW=GNEW+glst(IXX,II,JJ)*AD1*A(J,JJ) 461 G(IIND,I,J)=G(IIND,I,J)+GNEW*WF(IP) c c Rotation of A: DO 48 IX=1,3 IIND=3*(IA-1)+IX DO 48 I=1,3 DO 48 J=1,3 DO 48 K=1,3 C Zero out at the start of the summation IF(ISUM.EQ.1)AA(IIND,I,J,K)=0.0d0 AANEW=0.0d0 DO 49 IXX=1,3 AD0=A(IX,IXX) INDS=3*(IAS-1)+IXX DO 49 II=1,3 AD1=AD0*A(I,II) DO 49 JJ=1,3 AD2=AD1*A(J,JJ) DO 49 KK=1,3 49 AANEW=AANEW+aast(IXX,II,JJ,KK)*AD2*A(K,KK) 48 AA(IIND,I,J,K)=AA(IIND,I,J,K)+AANEW*WF(IP) ENDIF if(lwr)then call wrt2(IA,MX,ALPHA,'New Alpha:') call wrt2(IA,MX,G,'New G:') call wrt3(IA,MX,AA,'New A:') endif c c ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss c add the static polarizability component if(lpol.and.IOPOL.gt.0)then call ctd(IPO,opol,al_l,IOX) c c loop over other overlaps: do 51 iop=1,NO if(iop.ne.IPO)then if(LNAMES)then nt=NATNAME(iop) else nt=NATS endif c c set scale factor for alpha: if(neighbor(IPO,1).eq.iop.or.neighbor(IPO,2).eq.iop)then if(IOPOL.eq.1)then write(3,3009)IOPOL,IPO,iop 3009 format('IOPOL=',I1,',',i4,' and ',i4,' overlapped - skipped') goto 51 endif if(neighbor(IPO,1).eq.iop)then ix=1 else ix=2 endif c overlapping - only body+ further closer end c and do not correct atoms on IPO termini c is atom IA part of ends2 if IPO?: c ix is the index over the neighbors do 564 ii=1,endn(ix,2,IPO) if(endlist(ix,2,IPO,ii).eq.IA)then write(3,*)IA,' atom on tip, no correction by ',iop,ix goto 51 endif 564 continue c is atom IA part of end1 if ipo?: ic=0 do 563 ix=1,2 do 563 jend=1,2 do 563 ii=1,endn(ix,jend,iop) if(endlist(ix,jend,iop,ii).eq.IA)then if(ix.eq.1)then ic=2 else ic=1 endif goto 5631 endif 563 continue 5631 nn=endn(ic,1,iop)+bodyn(iop) else c non-overlapping or minor overlaps - consider middle part of iop nn=endn(1,1,iop)+endn(2,1,iop)+bodyn(iop) endif sf=dble(nn)/dble(nt) if(LWR)write(3,3111)iop,nn,nt 3111 format(' iop:',i4,' sf:',i4,'/',i4) if(LDATOM)then c distance to differentiated atom: do 157 ix=1,3 rl(ix)=RB(ix,IA)/bohr-RB(ix,IA)/bohr 157 rt(ix)=RB(ix,IA)/bohr-coo(iop,ix) else c distance to overlap center: do 156 ix=1,3 rl(ix)=coo(IPO,ix)-RB(ix,IA)/bohr 156 rt(ix)=coo(IPO,ix)-coo(iop,ix) endif rt2=rt(1)**2+rt(2)**2+rt(3)**2 do 158 ix=1,3 158 rj(ix)=coo(iop,ix)-RB(ix,IA)/bohr if(lq)then sfac=DBLE(ND(IA))/zsum(IPO) else sfac=1.0d0/zsum(IPO) endif if(rt2.gt.0.0d0)then c T and t -tensor: call mktt(rt,tt) call mkttt(rt,ttt) c c alphaj.T: call ctd(iop,opol,al,IOX) call ct(at,al,tt) do 603 IX=1,3 IIND=3*(IA-1)+IX c alphaj.t: call ctd(IX,ttt,ts,3) call ct(att,al,ts) c alphaj. t.alpha_l: call ct(atta,att,al_l) c alpha_l t.alphaj + alphaj.t. alpha_l : call std(atta) c alpha T.alphader: call ctd(IIND,ALPHA,alder,3*MX) call ct(atad,at,alder) c alphader T.alpha + alpha.T alphader : call std(atad) do 604 ig=1,3 do 604 id=1,3 604 tc(ig,id)=atad(ig,id)+atta(ig,id)*sfac call multt(tc,sf*WF(IP)) c alpha=alpha+alder T al + al T alder: call addtoG(IIND,ALPHA,3*MX,tc) c ul_ab = EPS_agd rl_g ader_db call dera(ul,rl,alder) c uj_ab = EPS_agd rj_g a_db call dera(uj,rj,al) c ul.T call ct(ult,ul,tt) c uj.T call ct(ujt,uj,tt) c ul.T.aj call ct(ultaj,ult,al) c uj.T.ader call ct(ujtad,ujt,alder) do 617 ig=1,3 do 617 id=1,3 617 tc(ig,id)=-0.5d0* ujtad(id,ig) c17 tc(ig,id)=-0.5d0*(ultaj(id,ig)+ujtad(id,ig)) c eps rl al_l call dera(vl,rl,al_l) c eps rl al_l t call ct(vlt,vl,ts) c eps rl al_l t aj call ct(vltaj,vlt,al) c eps rl alder t call ct(ujts,uj,ts) c eps rl alder t a_l call ct(ujtsal,ujts,al_l) do 619 ig=1,3 do 619 id=1,3 i1=id+1 if(i1.gt.3)i1=1 i2=i1+1 if(i2.gt.3)i2=1 c if(i1.eq.IX)tc(ig,id)=tc(ig,id)-0.5d0*sfac*al_l(i2,ig) c if(i2.eq.IX)tc(ig,id)=tc(ig,id)+0.5d0*sfac*al_l(i1,ig) 619 tc(ig,id)=tc(ig,id)-0.5d0*( 0.0d0* ujtsal(id,ig))*sfac c19 tc(ig,id)=tc(ig,id)-0.5d0*(vltaj(id,ig)+ujtsal(id,ig))*sfac call multt(tc,sf*WF(IP)) call addtoG(IIND,G,3*MX,tc) c alder.T: call ct(adt,alder,tt) c alder.T .al: call ct(adta,adt,al) c al.t: call ct(att,al,ts) c al.t .al: call ct(atta,att,al) DO 606 IG=1,3 DO 606 ID=1,3 DO 606 IE=1,3 606 AA(IIND,IG,ID,IE)=AA(IIND,IG,ID,IE) 1 +1.5d0*(rj(ID)*adta(IG,IE)+rj(IE)*adta(IG,ID)) 1 +1.5d0*(rj(ID)*atta(IG,IE)+rj(IE)*atta(IG,ID))*sfac DO 605 IG=1,3 DO 605 IE=1,3 605 AA(IIND,IG,IE,IE)=AA(IIND,IG,IE,IE) 1 -(rj(1)*adta(1,IG)+rj(2)*adta(2,IG)+rj(3)*adta(3,IG)) 1 -(rj(1)*atta(1,IG)+rj(2)*atta(2,IG)+rj(3)*atta(3,IG))*sfac 603 continue c endif endif 51 continue c if(lwr)then call wrt2(IA,MX,ALPHA,'Corrected Alpha:') call wrt2(IA,MX,G,'Corrected G:') call wrt3(IA,MX,AA,'Corrected A:') endif endif c ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss ENDIF 45 CONTINUE c IPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIPIP C 20 CONTINUE c IA,JA C C If only interaction among transfered atoms selcted, the other C force constants will be set to zero in final force field: IF(LSELECT)THEN DO 511 IA=1,NAT DO 521 IO=1,NO 521 IF(INDBIG(IO,IA).GT.0)GOTO 511 WRITE(3,*)' Atom ',IA,' taken out of the FF' DO 531 IX=1,3 INDIB=3*IA-3+IX DO 531 JA=1,NAT DO 531 JX=1,3 INDJB=3*JA-3+JX FFB(INDIB,INDJB)=0.0d0 531 FFB(INDJB,INDIB)=0.0d0 511 CONTINUE ENDIF RETURN C 19 WRITE(3,3000) 3000 FORMAT(1X,' Old matrix element New element',/, 1 1X,' atoms xi xj F [a.u.] F',' #it',/, 2 1X,'-------------------------------------------------') DO 1 IA=1,NS DO 3 I=1,3 TB(I) =0.0d0 3 T3(I)=0.0d0 DO 4 ISAT=1,ABS(IND(IA))+1 DO 4 J=1,3 XS(ISAT,J)=RS(J,IS(IA,ISAT)) XB(ISAT,J)=RB(J,IB(IA,ISAT)) TB(J)=TB(J)+XB(ISAT,J) 4 T3(J)=T3(J)+XS(ISAT,J) DO 5 ISAT=1,ABS(IND(IA))+1 DO 5 J=1,3 XS(ISAT,J)=XS(ISAT,J)-T3(J)/DBLE(ABS(IND(IA))+1.0d0) 5 XB(ISAT,J)=XB(ISAT,J)-TB(J)/DBLE(ABS(IND(IA))+1.0d0) IF(IND(IA).LT.0)THEN WRITE(3,*)' Space inversion applied for vicinity of atom ',IA DO 9 ISAT=1,ABS(IND(IA))+1 DO 9 J=1,3 9 XS(ISAT,J)=-XS(ISAT,J) ENDIF ITU=ABS(IND(IA))+1 CALL DOMATRIX(IERR) IF(IERR.NE.0)CALL REPORT('Erroneous stop') DO 10 I=1,3 DO 10 J=1,3 10 AM(IA,I,J)=A(I,J) DO 1 JA=1,IA DO 11 I=1,3 DO 11 J=1,3 11 B(I,J)=AM(JA,I,J) DO 6 I=1,3 INDI=3*IB(IA,1)-3+I JSTART=1 IF(IA.EQ.JA)JSTART=I DO 6 J=JSTART,3 INDJ=3*IB(JA,1)-3+J FOLD=FFB(INDI,INDJ) FFB(INDI,INDJ)=0.0d0 DO 7 II=1,3 IL=3*IS(IA,1)-3+II DO 7 JJ=1,3 JS=3*IS(IA,1)-3+JJ 7 FFB(INDI,INDJ)=FFB(INDI,INDJ)+A(I,II)* 1FFS((IL-1)*3*NATS_max+JS)*B(J,JJ) IF(I.EQ.1.AND.J.EQ.I.AND.JA.EQ.IA)THEN WRITE(3,3001)IA,JA,I,J,FOLD,FFB(INDI,INDJ),IIT ELSE WRITE(3,3001)IA,JA,I,J,FOLD,FFB(INDI,INDJ) ENDIF FFB(INDJ,INDI)=FFB(INDI,INDJ) 6 CONTINUE 3001 FORMAT(4I5,2F11.5,I5,F10.7,F7.4) 1 CONTINUE RETURN END c SUBROUTINE READRX(N0,N,R,ND,IBONDS,NB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),IBONDS(4,N0),NB(N0) READ(2,*) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > MX !') do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(IBONDS(IB,i),IB=1,4) NB(i)=0 DO 1 IB=1,4 1 IF(IBONDS(IB,i).NE.0)NB(i)=IB RETURN END SUBROUTINE READRXS(N0,N,R,ND,NB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),NB(N0) READ(2,*) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > MX !') do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i) 1 NB(i)=0 RETURN END SUBROUTINE WRITERX(N0,N,R,ND,IB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),IB(4,N0) WRITE(2,2000) 2000 FORMAT('FORCE FIELD IMPROVEMENT') write(2,2002)N 2002 format(I6) CHARGE=0.0d0 DO 1 I=1,N 1 WRITE(2,2001)ND(I),R(1,I),R(2,I),R(3,I),(IB(J,I),J=1,4), 10,0,0,CHARGE 2001 FORMAT(I5,3F12.6,7I5,F8.4) RETURN END SUBROUTINE READFF(NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*NA,3*NA) NAT3=3*NA N=NAT3 C Read in the lower triangle of FF, C written in parts as n1,n1 C . . C ln,n1 . ln,ln C . . . C . . . n3,n3 C . . . . C n,n1 . . n,n3 N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 READ(2,17)(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(4X,5D14.6) C DO 31 I=1,N DO 31 J=I,N 31 FCAR(I,J)=FCAR(J,I) RETURN END SUBROUTINE WRITEAFF(MX,NAT,FF,EFC,CUTOFF,R,LWR,j33,a33) IMPLICIT none integer*4 MX,NAT,ia,ja,ii,jj,ix,jx,nn,JSTART,j33(*) REAL*8 FF(3*NAT,3*NAT),EFC,CUTOFF,R(3,MX),ctf,a33(*),ft logical LWR c c zero FF for distant atoms: if(CUTOFF.NE.0.0d0)then ctf=CUTOFF**2 do 1 ia=1,NAT do 1 ja=1,ia if((r(1,ja)-r(1,ia))**2+ 1 (r(2,ja)-r(2,ia))**2+(r(3,ja)-r(3,ia))**2.ge.ctf)then do 2 ix=1,3 ii=3*(ia-1)+ix do 2 jx=1,3 jj=3*(ja-1)+jx FF(jj,ii)=0.0d0 2 FF(ii,jj)=0.0d0 if(LWR)write(3,3998)ia,ja 3998 format(2i5,' skipped because of the distance') endif 1 continue endif c c write non-zero elements open(33,file='AFILE.FC',form='unformatted') do 3 ia=1,NAT do 3 ix=1,3 ii=3*(ia-1)+ix nn=0 do 4 ja=ia,NAT JSTART=1 IF(ia.EQ.ja)JSTART=ix do 4 jx=JSTART,3 jj=3*(ja-1)+jx ft=FF(ii,jj) if(dabs(ft).gt.EFC)then nn=nn+1 j33(nn)=jj if(ii.eq.jj)ft=ft/2.0d0 a33(nn)=ft endif 4 continue 3 write(33)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) close(33) return end SUBROUTINE WRITEFF(NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*NA,3*NA) N=NA*3 N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N 130 WRITE(2,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END SUBROUTINE DOMATRIX(IERR) C This is a Fortran version of Petr Malon's SUBROUTINE Amoeba IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=140) DIMENSION ANG(3),P(4,3),Y(4),PBAR(3),PR(3),PRR(3) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER,RMS IERR=0 C STARTING VALUES FOR THE ITERATION: ANG(1)=2.1d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(1)=0.0d0 ANG(2)=0.1d0 ANG(3)=2.1d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(1)=0.2d0 ANG(2)=2.1d0 ANG(3)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=1.4d0 ANG(2)=0.0d0 ANG(3)=0.1d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=1.0d-9 ITMAX=1500 ITER=0 99999 ILO=1 IF(Y(1).GT.Y(2))THEN IHI=1 INHI=2 ELSE IHI=2 INHI=1 ENDIF DO 5 I=1,4 IF(Y(I).LT.Y(ILO))ILO=I IF(Y(I).GT.Y(IHI))THEN INHI=IHI IHI=I ELSE IF(Y(I).GT.Y(INHI))THEN IF(I.NE.IHI)INHI=I ENDIF ENDIF 5 CONTINUE IF(ABS(Y(IHI))+ABS(Y(ILO)).LT.1.0d-10)THEN RTOL=0.0d0 ELSE RTOL=ABS(Y(IHI)-Y(ILO))/(ABS(Y(IHI))+ABS(Y(ILO)))*2.0d0 ENDIF IF(RTOL.LT.FTOL)then RETURN endif IF(ITER.EQ.ITMAX)THEN WRITE(3,*)' Rotation has not converged !' IERR=1 RETURN ENDIF ITER=ITER+1 DO 55 I=1,3 55 PBAR(I)=0.0d0 DO 6 I=1,4 IF(I.NE.IHI)THEN DO 7 J=1,3 7 PBAR(J)=PBAR(J)+P(I,J) ENDIF 6 CONTINUE DO 8 J=1,3 PBAR(J)=PBAR(J)/3.0d0 8 PR(J)=2.0d0*PBAR(J)-P(IHI,J) YPR=FU(PR) IF(YPR.LE.Y(ILO))THEN DO 9 J=1,3 9 PRR(J)=2.0d0*PR(J)-PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(ILO))THEN DO 10 J=1,3 10 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 11 J=1,3 11 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ELSE IF(YPR.GE.Y(INHI))THEN IF(YPR.LT.Y(IHI))THEN DO 12 J=1,3 12 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF DO 13 J=1,3 13 PRR(J)=0.5d0*P(IHI,J)+0.5d0*PBAR(J) YPRR=FU(PRR) IF(YPRR.LT.Y(IHI))THEN DO 14 J=1,3 14 P(IHI,J)=PRR(J) Y(IHI)=YPRR ELSE DO 15 I=1,4 IF(I.NE.ILO)THEN DO 16 J=1,3 PR(J)=0.5d0*(P(I,J)+P(ILO,J)) 16 P(I,J)=PR(J) Y(I)=FU(PR) ENDIF 15 CONTINUE ENDIF ELSE DO 17 J=1,3 17 P(IHI,J)=PR(J) Y(IHI)=YPR ENDIF ENDIF GOTO 99999 END FUNCTION FU(ANG) implicit none INTEGER*4 MXB,NAT,ITER,I REAL*8 S1,S2,S3,C1,C2,C3,TX,TY,TZ,A,XS,XB,RTOL,RMS,ANG(3),FU PARAMETER (MXB=140) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER,RMS S1=SIN(ANG(1)) S2=SIN(ANG(2)) S3=SIN(ANG(3)) C1=COS(ANG(1)) C2=COS(ANG(2)) C3=COS(ANG(3)) C ang1 = theta C ang2 = phi C ang3 = cappa A(1,1)=C1*C2*C3-S2*S3 A(1,2)=C1*S2*C3+C2*S3 A(1,3)=-S1*C3 A(2,1)=-C1*C2*S3-S2*C3 A(2,2)=-C1*S2*S3+C2*C3 A(2,3)=S1*S3 A(3,1)=S1*C2 A(3,2)=S1*S2 A(3,3)=C1 RMS=0.0d0 DO 1 I=1,NAT TX=A(1,1)*XS(I,1)+A(1,2)*XS(I,2)+A(1,3)*XS(I,3) TY=A(2,1)*XS(I,1)+A(2,2)*XS(I,2)+A(2,3)*XS(I,3) TZ=A(3,1)*XS(I,1)+A(3,2)*XS(I,2)+A(3,3)*XS(I,3) 1 RMS=RMS+(TX-XB(I,1))*(TX-XB(I,1))+ 1 (TY-XB(I,2))*(TY-XB(I,2))+(TZ-XB(I,3))*(TZ-XB(I,3)) FU=RMS RETURN END C SUBROUTINE REPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(3,3000) WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(3,*)RS WRITE(6,*)RS WRITE(3,3001) WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') CLOSE(3) CLOSE(2) STOP END c SUBROUTINE SREPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(3,3000) WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(3,*)RS WRITE(6,*)RS WRITE(3,3001) WRITE(6,3001) 3001 FORMAT(80(1H*),/,/, 1'PROGRAM CONTINUES BECAUSE LHALTONERROR DISABLED') END C SUBROUTINE READTEN(N0,P,A,LAPT,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=4000) DIMENSION P(N0,3,3),A(N0,3,3),AI(MX,3,3),AJ(MX,3,3),PV(MX,3,3), 1R(3,MX),X(3,NAT) LOGICAL LAPT BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR READ (15,*) NOAT NAT=NOAT WRITE(3,*)NAT,' atoms' WRITE(*,*)NAT,' atoms' C C Read dipole derivatives: DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) WRITE(6,*)' Dipole derivatives read - in ' WRITE(3,*)' Dipole derivatives read - in ' C C If no axial tensor, then construct its polar part only: IF(LAPT)THEN DO 2203 L=1,NAT DO 2203 J=1,3 DO 2203 I=1,3 SUM=0.0d0 DO 2204 ID=1,3 DO 2204 IG=1,3 2204 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) AJ(L,J,I)=0.0d0 PV(L,J,I)=P(L,J,I) 2203 AI(L,J,I)=SUM WRITE(6,*)' APT part of AAT constructed' WRITE(3,*)' APT part of AAT constructed' GOTO 1 ENDIF DO 220 L=1,NOAT DO 220 J=1,3 220 READ (15,*) (AI(L,J,I),I=1,3) DO 230 L=1,NOAT DO 230 J=1,3 230 READ (15,*) (AJ(L,J,I),I=1,3) DO 100 L=1,NOAT DO 100 J=1,3 100 READ (15,*) (PV(L,J,I),I=1,3) WRITE(6,*)' VCD parameters read in ...' WRITE(3,*)' VCD parameters read in ...' C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) WRITE(6,*)'DOG used' WRITE(3,*)'DOG used' C C Total axial tensor = electronic + nuclear: DO 4 L=1,NAT DO 4 I=1,3 DO 4 J=1,3 4 A(L,I,J)=AI(L,I,J)+AJ(L,I,J) C C Put axial tensor in the origin on the atom L: DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM-0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(3,*)' AAT shifted to atomic origin' WRITE(6,*)' AAT shifted to atomic origin' RETURN END C SUBROUTINE WRITETEN(N0,P,A,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=4000) DIMENSION P(N0,3,3),A(N0,3,3),R(3,MX),X(3,NAT) BOHR=0.52917705993d0 DO 3 L=1,NAT DO 3 I=1,3 3 R(I,L)=X(I,L)/BOHR DO 2201 L=1,NAT DO 2201 J=1,3 DO 2201 I=1,3 SUM=0.0d0 DO 2202 ID=1,3 DO 2202 IG=1,3 2202 SUM=SUM+0.25d0*EPS(I,IG,ID)*R(IG,L)*P(L,J,ID) 2201 A(L,J,I)=A(L,J,I)+SUM WRITE(3,*)' AAT shifted back to laboratory system' WRITE(6,*)' AAT shifted back to laboratory system' OPEN(15,FILE='FILE.TEN') WRITE(15,1500) NAT,NAT-6,0 1500 FORMAT(3I5) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(15,1501) (P(L,J,I),I=1,3),L 1501 FORMAT(3F14.8,I5) DO 220 L=1,NAT DO 220 J=1,3 220 WRITE(15,1501) (A(L,J,I),I=1,3),L Z0=0.0d0 DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(15,1501) (Z0,I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(15,1501) (P(L,J,I),I=1,3),L CLOSE(15) WRITE(6,*)' FILE.TEN written ...' WRITE(3,*)' FILE.TEN written ...' RETURN END C FUNCTION EPS(I,J,K) REAL*8 EPS INTEGER*4 I,J,K EPS=0.0d0 IF (I.EQ.1.AND.J.EQ.2.AND.K.EQ.3)EPS= 1.0d0 IF (I.EQ.1.AND.J.EQ.3.AND.K.EQ.2)EPS=-1.0d0 IF (I.EQ.2.AND.J.EQ.3.AND.K.EQ.1)EPS= 1.0d0 IF (I.EQ.2.AND.J.EQ.1.AND.K.EQ.3)EPS=-1.0d0 IF (I.EQ.3.AND.J.EQ.1.AND.K.EQ.2)EPS= 1.0d0 IF (I.EQ.3.AND.J.EQ.2.AND.K.EQ.1)EPS=-1.0d0 RETURN END C SUBROUTINE VCDD0(VCD,VEL,DIP,C,NAT) C STOLEN FROM CADPAC, INDICES ARE DIFFERENT !!, VEL WITH NUCLEI IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) INTEGER*4 ALPHA,BETA,GAMMA,DELTA,LAMBDA PARAMETER (MX=4000) DIMENSION VCD(MX,3,3),VEL(MX,3,3),DIP(MX,3,3), 1C(3,NAT) DO 90 BETA =1,3 DO 90 GAMMA=1,3 DO 90 DELTA=1,3 SKEW=0.25D0*EPS(BETA,GAMMA,DELTA) IF (DABS(SKEW).GT.1.0D-10)THEN DO 80 ALPHA=1,3 DO 80 LAMBDA=1,NAT 80 VCD(LAMBDA,ALPHA,BETA)=VCD(LAMBDA,ALPHA,BETA) 1 -SKEW*C(GAMMA,LAMBDA) 2 *(VEL(LAMBDA,ALPHA,DELTA)-DIP(LAMBDA,ALPHA,DELTA)) ENDIF 90 CONTINUE RETURN END C SUBROUTINE READTTT(MX,ALPHA,A,G,NAT) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3) READ(15,*) READ(15,*)NATC if(NAT.NE.NATC)call report('Wrong number of atoms in .TTT file') READ(15,*) READ(15,*) DO 1 I=1,3 READ(15,*) DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 READ(15,*)(ALPHA(IIND,I,J),J=1,2),(ALPHA(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 2 I=1,3 READ(15,*) DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 READ(15,*)(G(IIND,I,J),J=1,2),(G(IIND,I,J),J=1,3) READ(15,*) READ(15,*) DO 3 I=1,3 DO 3 J=1,3 READ(15,*) DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 READ(15,*)(A(IIND,I,J,K),K=1,2),(A(IIND,I,J,K),K=1,3) CLOSE(15) RETURN END C SUBROUTINE WRITETTT(MX,NAT,ALPHA,A,G,s) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION ALPHA(3*MX,3,3),G(3*MX,3,3),A(3*MX,3,3,3) character*(*)s OPEN(2,FILE=s) WRITE(2,2000)NAT 2000 FORMAT(' ROA tensors, cartesian derivatives',/,I4,' atoms',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 WRITE(2,2001)L,IX,(ALPHA(IIND,I,J),J=1,3) 2001 FORMAT(I5,1H ,I1,3F18.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx jy jz') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 WRITE(2,2001)L,IX,(G(IIND,I,J),J=1,3) WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX 3 WRITE(2,2001)L,IX,(A(IIND,I,J,K),K=1,3) CLOSE(2) WRITE(3,*)s//' written ' WRITE(6,*)s//' written ' RETURN END c ============================================================ SUBROUTINE TTTA(ALPHA,A,ICO,coo,iwr) c Tensors, not derivatives C Transforms A to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C coo supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,IA,IB,IC,iwr real*8 ALPHA(3,3),A(3,3,3),coo(3),SIGN,SUM C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C c A .. last index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM= 1coo(1)*ALPHA(IA,2)+coo(2)*ALPHA(IA,2)+coo(3)*ALPHA(IA,3) A(IC,IB,IA)=A(IC,IB,IA)+ 1SIGN*(SUM-1.5d0*(coo(IB)*ALPHA(IA,IC)+coo(IC)*ALPHA(IA,IB))) IF(ICO.EQ.3)A(IC,IB,IA)=0.0d0 3 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(3,*)' A transformed into overlap origin' IF(ICO.EQ.2)WRITE(3,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(3,*)' A set to zero' endif RETURN END c ============================================================ SUBROUTINE TTTG(ALPHA,G,ICO,coo,iwr) c Tensors, not derivatives C Transforms G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C coo supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,IA,IB,IG,ID,iwr real*8 ALPHA(3,3),G(3,3),coo(3),SIGN,SUM,EPS C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*coo(IG)*ALPHA(IA,ID) G(IA,IB)=G(IA,IB)+SUM IF(ICO.EQ.3)G(IA,IB)=0.0d0 2 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(3,*)' G transformed into overlap origin' IF(ICO.EQ.2)WRITE(3,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(3,*)' G set to zero' endif RETURN END c ============================================================ SUBROUTINE TTT1(ALPHA,A,G,ICO,coo,iwr) c Tensors, not derivatives C Transforms A and G to local origins (ICO=1), c common origin (ICO=2), C sets them to zero (ICO=3). C coo supposed to be passed in aus C the static limit G/w is in G C implicit none INTEGER*4 ICO,IA,IB,IG,ID,IC,iwr real*8 ALPHA(3,3),G(3,3),A(3,3,3),coo(3),SIGN,SUM,EPS C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 C DO 2 IA=1,3 DO 2 IB=1,3 SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*0.5d0*EPS(IB,IG,ID)*coo(IG)*ALPHA(IA,ID) G(IA,IB)=G(IA,IB)+SUM IF(ICO.EQ.3)G(IA,IB)=0.0d0 2 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(3,*)' G transformed into overlap origin' IF(ICO.EQ.2)WRITE(3,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(3,*)' G set to zero' endif c A .. last index (IA) is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 SUM=0.0d0 IF(IB.EQ.IC)SUM= 1coo(1)*ALPHA(IA,1)+coo(2)*ALPHA(IA,2)+coo(3)*ALPHA(IA,3) A(IC,IB,IA)=A(IC,IB,IA)+ 1SIGN*(SUM-1.5d0*(coo(IB)*ALPHA(IA,IC)+coo(IC)*ALPHA(IA,IB))) IF(ICO.EQ.3)A(IC,IB,IA)=0.0d0 3 CONTINUE if(iwr.eq.1)then IF(ICO.EQ.1)WRITE(3,*)' A transformed into overlap origin' IF(ICO.EQ.2)WRITE(3,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(3,*)' A set to zero' endif RETURN END c ============================================================ SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R) c Tensors derivatives C Transforms A and G to local origins (ICO=1) or c common origin (ICO=2); C sets A and G to zero for ICO=3. C R supposed to be passed in in A C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) PARAMETER (MX=4000) DIMENSION ALPHA(N,3,3),G(N,3,3),A(N,3,3,3),R(3,NAT),X(3,MX) C C transfer into atomic coordinates: DO 5 I=1,3 DO 5 J=1,NAT 5 X(I,J)=R(I,J)/0.529177D0 C IF(ICO.EQ.1)SIGN=1.0d0 IF(ICO.EQ.2)SIGN=-1.0d0 IF(ICO.EQ.3)SIGN=0.0d0 W=1.0d0 C supposedly the static limit G/w is calculated C DO 2 IA=1,3 DO 2 IB=1,3 DO 2 L=1,NAT DO 2 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 DO 1 IG=1,3 DO 1 ID=1,3 1 SUM=SUM+SIGN*(0.5d0*W)*EPS(IB,IG,ID)*X(IG,L)*ALPHA(IIND,IA,ID) G(IIND,IA,IB)=G(IIND,IA,IB)+SUM IF(ICO.EQ.3)G(IIND,IA,IB)=0.0d0 2 CONTINUE IF(ICO.EQ.1)WRITE(3,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(3,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(3,*)' G set to zero' C c A IA index is electric DO 3 IA=1,3 DO 3 IB=1,3 DO 3 IC=1,3 DO 3 L=1,NAT DO 3 IE=1,3 IIND=3*(L-1)+IE SUM=0.0d0 IF(IB.EQ.IC)THEN DO 4 ID=1,3 4 SUM=SUM+X(ID,L)*ALPHA(IIND,IA,ID)*SIGN ENDIF A(IIND,IA,IB,IC)=A(IIND,IA,IB,IC)+SUM 1-1.5d0*SIGN*(X(IB,L)*ALPHA(IIND,IA,IC)+X(IC,L)*ALPHA(IIND,IA,IB)) IF(ICO.EQ.3)A(IIND,IA,IB,IC)=0.0d0 3 CONTINUE IF(ICO.EQ.1)WRITE(3,*)' A transformed into atomic origins' IF(ICO.EQ.2)WRITE(3,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(3,*)' A set to zero' RETURN END c SUBROUTINE readopt(NS,IB,IS,IND,MXT,MXB,INDBIG,IOX,NO, 1NAT,LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LRAM,LSTRICT,LINV,LSELECT, 2LAPT,LALF,LRDO,LNAMES,NAME,LPOLYMER,LHALTONERROR,CUTOFF, 2n1,n2,okind,LXX3,LXX4,ALIM,L36,L37,lq,LDATOM,IOPOL,lsos,lmut, 1lmpt,IFC,EFC,iweights,lgrad,sgrad,fgrad,alphf,lnmr,lspin,LQUAD, 1SOSNM) implicit none INTEGER*4 NS,IB,IS,IND,MXT,MXB,INDBIG,IOX,NO,IOPOL, 1NAT,n1,n2,okind,I,IA,IAB,IC,IEND,II,IO,ISMALL,IWG,J,JA,NDD,nt, 1IFC,iweights REAL*8 CUTOFF,ALIM,EFC,alphf DIMENSION IB(MXT,MXB),IS(MXT,MXB),IND(MXT),INDBIG(IOX,MXT), 1NDD(20) CHARACTER*7 STR,STRC CHARACTER*80 NAME(IOX) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT,LAPT, 1LALF,LRDO,LPOLYMER,LNAMES,LHALTONERROR,LRAM,LXX3,LXX4,L36, 2L37,loldformat,lq,LDATOM,lsos,lmut,lmpt,lgrad,sgrad,fgrad, 3lspin,lnmr,LQUAD,lnew,SOSNM integer*4,allocatable::bb(:,:) c alphf=1.0d-3 n1=0 n2=0 okind=0 iweights=0 IOPOL=1 IFC=0 ALIM=0.00d0 LDATOM=.false. lsos=.false. SOSNM=.true. lgrad=.false. sgrad=.false. fgrad=.false. lmpt=.false. lnmr=.false. lspin=.false. lmut=.true. LRDO=.FALSE. loldformat=.true. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LQUAD=.FALSE. lnew=.false. LOFF=.TRUE. lq=.false. LDIA=.TRUE. LXX3=.false. LXX4=.false. L36=.false. L37=.false. LHALTONERROR=.TRUE. CUTOFF=-1.0d5 EFC=0.0d0 LROA=.FALSE. LRAM=.FALSE. LSTRICT=.TRUE. LINV=.FALSE. LSELECT=.FALSE. LNAMES=.FALSE. IWG=1 NO=0 9 READ(2,2000)STR 2000 FORMAT(A7) WRITE(3,2000)STR c 20i3 format: IF(STR(1:7).EQ.'OLDFORM')READ(2,*)loldformat c atom interval: IF(STR(1:4).EQ.'N1N2') READ(2,*)n1,n2 c polarizability correction: IF(STR(1:5).EQ.'IOPOL') READ(2,*)IOPOL c limit for anharmonic constants: IF(STR(1:4).EQ.'ALIM') READ(2,*)ALIM c c if okind=0 fragment weighing based on center distances c if okind=1 fragment weighing based on RMS deviations c if okind=2 fragment weighing based on RMS molecular : IF(STR(1:3).EQ.'OKI') READ(2,*)okind c IWG=0 best overlap, IWG=1 weighting: IF(STR(1:3).EQ.'IWG') READ(2,*)IWG c named fragments: IF(STR(1:6).EQ.'LNAMES') READ(2,*)LNAMES c selected atom transfer only: IF(STR.EQ.'LSELECT') READ(2,*)LSELECT c invert fragment (obsolete): IF(STR(1:4).EQ.'LINV') READ(2,*)LINV c IFC=0 FF to FILE.FC c IFC=1 FF to AFILE.FC (based either on CUTOFF or EFC): IF(STR(1:3).EQ.'IFC') READ(2,*)IFC c polarizability distributed according WEIGHTS.TXT: c 0 (def) .. do not use weights c 1.. real c 2.. imaginary c 3.. absolute value c if iweights < 0, read it fro each frequency IF(STR.EQ.'IWEIGHT') READ(2,*)iweights c criterion for AFILE.FC selection: IF(STR(1:3).EQ.'EFC') READ(2,*)EFC c transfor atoms in beds only: IF(STR.EQ.'LSTRICT') READ(2,*)LSTRICT c extensive writing to CCT.OUT: IF(STR(1:3).EQ.'LWR') READ(2,*)LWR c construct atomic axial tensor from atomic polar tensor: IF(STR(1:4).EQ.'LAPT') READ(2,*)LAPT c construct G' and A from electric polarizability: IF(STR(1:4).EQ.'LALF') READ(2,*)LALF c Alpha from SMALL_ALPHA.TTT for small IF(STR(1:4).EQ.'LRDO') READ(2,*)LRDO c Transfer of dynamic polarizabilities: IF(STR(1:4).EQ.'LSOS') READ(2,*)lsos c Dynamic polarizabilities dependent on nm (else cm-1): IF(STR(1:4).EQ.'SOSN') READ(2,*)SOSNM c Use matrix perturbation theory for sos: IF(STR(1:4).EQ.'LMPT') READ(2,*)lmpt c Correction of polarizability for mutual group interaction: IF(STR(1:4).EQ.'LMUT') READ(2,*)lmut c Transfer Raman tensor derivatives: IF(STR(1:4).EQ.'LRAM') READ(2,*)LRAM c Transfer ROA derivatives: IF(STR(1:4).EQ.'LROA') READ(2,*)LROA c Transfer off-diagional force field elements: IF(STR(1:4).EQ.'LOFF') READ(2,*)LOFF c Transfer diagional force field elements: IF(STR(1:4).EQ.'LDIA') READ(2,*)LDIA c Transfer gradient: IF(STR(1:4).EQ.'LGRA') READ(2,*)lgrad c Make sum of gradients from all fragments IF(STR(1:4).EQ.'SGRA') READ(2,*)sgrad c Transfer gradient with FF in pairs IF(STR(1:4).EQ.'FGRA') READ(2,*)fgrad c fgrad method stabilization constants IF(STR(1:4).EQ.'ALPH') READ(2,*)alphf c Transfer NMR shifts: IF(STR(1:4).EQ.'LNMR') READ(2,*)lnmr c Transfer NMR spin-spin (J) constants: IF(STR(1:4).EQ.'LSPI') READ(2,*)lspin c Transfer cubic force field: IF(STR(1:4).EQ.'LXX3') READ(2,*)LXX3 c Transfer quartic force field: IF(STR(1:4).EQ.'LXX4') READ(2,*)LXX4 c Transfer two/three atomic anharmonic FF: IF(STR(1:4).EQ.'L36 ') READ(2,*)L36 c charge center (vs geometry center), for pol interaction: IF(STR(1:4).EQ.'LQ ') READ(2,*)lq c use the L37 (three-atomic) indexing also for L36 (two-atomic): IF(STR(1:4).EQ.'L37 ') READ(2,*)L37 c Transfer absorption derivatives: IF(STR(1:4).EQ.'LABS') READ(2,*)LABS c Transfer vibrational circular dichroism derivatives: IF(STR(1:4).EQ.'LVCD') READ(2,*)LVCD c Transfer quadrupole derivatives: IF(STR(1:4).EQ.'LQUA') READ(2,*)LQUAD c New format of overlap map: IF(STR(1:4).EQ.'LNEW') READ(2,*)lnew c Stop transfer if an error in atom maps occurs: IF(STR(1:7).EQ.'LHALTON')READ(2,*)LHALTONERROR c Cutoff for atomic distances considered: IF(STR(1:6).EQ.'CUTOFF') READ(2,*)CUTOFF c In polarization correction, use distances from c differentiated atom: IF(STR(1:6).EQ.'LDATOM') READ(2,*)LDATOM c outdated definition of the overlap: IF(STR.EQ.'NOPOLYM')THEN LPOLYMER=.false. REWIND 2 READ(2,*)NS DO 1 I=1,NS READ(2,*)IB(I,1),IS(I,1),IND(I) DO 1 J=1,ABS(IND(I)) 1 READ(2,*)IB(I,J+1),IS(I,J+1) C C NS atoms are assigned C IB(I) atom of the big molecule goes with the IS(I) atom of the smaller. C Each of them has IND(I) satelite atoms to define the reference frame. C Reference atoms are given for big as IB and for small as IS strings. C goto 9999 ENDIF IF(STR.EQ.'POLYMER')THEN c automatic definition of overlap: LPOLYMER=.true. READ(2,*)NO WRITE(3,*)NO,' overlaps found, atoms assigned automatically ...' IF(NO.GT.IOX)CALL REPORT(' NO > IOX ! Change the dimensions.') DO 2 I=1,NO READ(2,*)IC IF(IC.NE.I)THEN WRITE(3,*)IC,I WRITE(6,*)IC,I CALL REPORT('Error in overlap reading ') ENDIF c named fragments: IF(LNAMES)READ(2,2900)NAME(I) 2900 FORMAT(a80) c newest format if(lnew)then DO 52 IAB=1,NAT 52 INDBIG(I,IAB)=0 read(2,*)nt allocate(bb(2,nt)) backspace 2 read(2,*)nt,((bb(IO,IAB),IO=1,2),IAB=1,nt) do 53 IAB=1,nt 53 INDBIG(I,bb(1,IAB))=bb(2,IAB) deallocate(bb) else if(loldformat)then READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT) 2001 FORMAT(20I3) DO 5 IAB=1,NAT 5 INDBIG(I,IAB)=0 READ(2,2001)(INDBIG(I,IAB),IAB=1,NAT) else READ(2,*)(INDBIG(I,IAB),IAB=1,NAT) DO 51 IAB=1,NAT 51 INDBIG(I,IAB)=0 READ(2,*)(INDBIG(I,IAB),IAB=1,NAT) endif endif 2 CONTINUE c c Control Output: DO 4 IO=1,NO WRITE(3,2002)IO 2002 FORMAT(' Overlap number ',I5) IA=0 7 IA=IA+1 IEND=IA+19 IF(IEND.GT.NAT)IEND=NAT IF(IA.GT.NAT)GOTO 6 DO 8 II=1,(IEND-IA+1) 8 NDD(II)=II+IA-1 WRITE(3,2001)(NDD(II),II=1,IEND-IA+1) IA=IEND GOTO 7 6 WRITE(3,2001)(INDBIG(IO,IAB),IAB=1,NAT) 4 CONTINUE c c check, that atoms are not defined twice in the overlaps: do io=1,no do ia=1,nat ismall=indbig(io,ia) if(ismall.ne.0)then do ja=ia+1,nat if(ismall.eq.indbig(io,ja))then write(3,*)'overlap ',io write(3,*)'big ',ia,ja write(3,*)'small ',ismall call report('umbiqous overlap') endif enddo endif enddo enddo goto 9999 ENDIF BACKSPACE 2 READ(2,2000)STRC if(STR.eq.STRC)call report('Unknown option '//STRC) goto 9 c c 9999 close(2) RETURN END c SUBROUTINE CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1IA,JA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS) c find covalently bond atoms to IA..JA, oevrlap and make the c rotation matrix c IC=0 ... take atoms two bonds away only if directly c bonded do not suffice c IC=1 ... include all closest atoms by default c IC=2 ... include all transferred atoms implicit none integer*4 IOX,MXB,MX PARAMETER (MXB=140) integer*4 IA,IPASS,II,NAT,INDBIG(IOX,MX),IBA,IBONDS(4,MX), 1NBOB(*),JB,IIA,NATS,natname(*),IX,IND(*),iu,IERC,I,IC,IERR, 2IPO,J,JA real*8 RMS,RS(3,MX),RB(3,MX),TB(3),TS(3),rname(IOX,3,MX) logical LNAMES,LHALTONERROR,LSTRICT,LINV real*8 A,XS,XB,TOL integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS IERC=1 IPASS=0 666 ITU=0 IPASS=IPASS+1 C DO 27 II=1,NAT if(INDBIG(IPO,II).NE.0)then C Take IA and JA right away in the list: IF(II.EQ.IA.OR.II.EQ.JA.or.IC.eq.2)GOTO 277 C Take also atoms connected to IA or JA: DO 29 IBA=1,NBOB(II) 29 IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA)GOTO 277 C For second pass, take also atoms two bonds away: IF(IPASS.EQ.2.or.IC.eq.1)THEN DO 31 IBA=1,NBOB(II) IIA=IBONDS(IBA,II) DO 31 JB=1,NBOB(IIA) 31 IF(IBONDS(JB,IIA).EQ.IA.OR.IBONDS(JB,IIA).EQ.JA)GOTO 277 ENDIF c c This atom is of no use, try another one: GOTO 27 277 IF(INDBIG(IPO,II).LT.0.AND.LSTRICT)GOTO 27 C Include atom II into the transformation set: ITU=ITU+1 IF(ITU.GT.MXB)CALL REPORT('ITU > MXB !') if(LNAMES)NATS=natname(IPO) DO 28 IX=1,3 if(LNAMES)then XS(ITU,IX)=rname(IPO,IX,ABS(INDBIG(IPO,II))) else XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II))) endif 28 XB(ITU,IX)=RB(IX,II) IND(ITU)=II endif 27 CONTINUE c IF(ITU.LT.3)WRITE(3,*)' ITU < 3 !' IF(ITU.LT.3.AND.IPASS.EQ.1.and.IC.eq.0)GOTO 666 if(ITU.LT.3)then write(3,*)'CRMS:' write(3,*)'IA: ',IA,' JA: ',JA write(3,*)'IPO:',IPO write(3,*)'IC:',IC write(3,*)'NAT:',NAT write(3,*)'NATS:',NATS write(3,*)'INDBIG:' write(3,*)(INDBIG(IPO,II),II=1,NAT) DO 271 II=1,NAT write(3,*)'atom ', II if(INDBIG(IPO,II).NE.0)then write(3,*)'part of index' IF(II.EQ.IA.OR.II.EQ.JA.or.IC.eq.2)then write(3,*)'take it ia ja ic' GOTO 2771 endif write(3,*)NBOB(II),' bonds' DO IBA=1,NBOB(II) write(3,*)IBONDS(IBA,II) IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA)then write(3,*)'take it ia ja' goto 2771 endif enddo IF(IPASS.EQ.2.or.IC.eq.1)THEN DO IBA=1,NBOB(II) IIA=IBONDS(IBA,II) DO JB=1,NBOB(IIA) IF(IBONDS(JB,IIA).EQ.IA.OR.IBONDS(JB,IIA).EQ.JA)then write(3,*)'second pass' GOTO 2771 endif ENDDO ENDDO ENDIF ENDIF write(3,*)'This atom is of no use' GOTO 271 2771 IF(INDBIG(IPO,II).LT.0.AND.LSTRICT)then write(3,*)'lstrict' GOTO 271 endif write(3,*)'use this atom' 271 CONTINUE write(3,*)'Only following atoms found in the big molecule:' write(3,*)(IND(iu) ,iu=1,ITU) if(LHALTONERROR)then CALL REPORT(' ITU < 3 !') else CALL SREPORT(' ITU < 3 !') return endif endif c c C Calculate center of the coordinates that will be rotated: DO 33 I=1,3 TB(I) =0.0d0 33 TS(I)=0.0d0 DO 34 II=1,ITU DO 34 J=1,3 TB(J)=TB(J)+XB(II,J) 34 TS(J)=TS(J)+XS(II,J) DO 38 J=1,3 TB(J)=TB(J)/DBLE(ITU) 38 TS(J)=TS(J)/DBLE(ITU) C Subtract mass-center: DO 35 II=1,ITU DO 35 J=1,3 XS(II,J)=XS(II,J)-TS(J) 35 XB(II,J)=XB(II,J)-TB(J) C IF(LINV)THEN DO 39 II=1,ITU DO 39 J=1,3 39 XS(II,J)=-XS(II,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR) IF(IERR.EQ.1.AND.IPASS.EQ.1)THEN write(3,*)'domat error,try other pass' GOTO 666 ENDIF IF(IERR.EQ.1)then if(LHALTONERROR)then CALL REPORT('Program cannot continue') else CALL SREPORT('Program cannot continue') return endif endif IERC=0 return end C SUBROUTINE CRMS34(IPASS,NAT,INDBIG,IOX,MX,IPO, 1IA,JA,KA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,LINV,LSTRICT,rname,NATS) c find covalently bond atoms to IA,JA,KA overlap and make the c rotation matrix c IC=0 ... take atoms two bonds away only if directly c bonded do not suffice c IC=1 ... include all closest atoms by default c IC=2 ... include all transferred atoms implicit none integer*4 IOX,MXB,MX PARAMETER (MXB=140) integer*4 IA,IPASS,II,NAT,INDBIG(IOX,MX),IBA,IBONDS(4,MX), 1NBOB(*),JB,IIA,NATS,natname(*),IX,IND(*),iu,IERC,I,IC,IERR, 2IPO,J,JA,KA real*8 RMS,RS(3,MX),RB(3,MX),TB(3),TS(3),rname(IOX,3,MX) logical LNAMES,LHALTONERROR,LSTRICT,LINV real*8 A,XS,XB,TOL integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS IERC=1 IPASS=0 666 ITU=0 IPASS=IPASS+1 C DO 27 II=1,NAT if(INDBIG(IPO,II).NE.0)then C Take IA and JA right away in the list: IF(II.EQ.IA.OR.II.EQ.JA.or.II.EQ.KA.or.IC.eq.2)GOTO 277 C Take also atoms connected to IA or JA or KA: DO 29 IBA=1,NBOB(II) 29 IF(IBONDS(IBA,II).EQ.IA.OR.IBONDS(IBA,II).EQ.JA 1 .OR.IBONDS(IBA,II).EQ.KA)GOTO 277 C For second pass, take also atoms two bonds away: IF(IPASS.EQ.2.or.IC.eq.1)THEN DO 31 IBA=1,NBOB(II) IIA=IBONDS(IBA,II) DO 31 JB=1,NBOB(IIA) 31 IF(IBONDS(JB,IIA).EQ.IA.OR.IBONDS(JB,IIA).EQ.JA 1 .OR.IBONDS(JB,IIA).EQ.KA)GOTO 277 ENDIF c c This atom is of no use, try another one: GOTO 27 277 IF(INDBIG(IPO,II).LT.0.AND.LSTRICT)GOTO 27 C Include atom II into the transformation set: ITU=ITU+1 IF(ITU.GT.MXB)CALL REPORT('ITU > MXB !') if(LNAMES)NATS=natname(IPO) DO 28 IX=1,3 if(LNAMES)then XS(ITU,IX)=rname(IPO,IX,ABS(INDBIG(IPO,II))) else XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II))) endif 28 XB(ITU,IX)=RB(IX,II) IND(ITU)=II endif 27 CONTINUE c IF(ITU.LT.3)WRITE(3,*)' ITU < 3 !' IF(ITU.LT.3.AND.IPASS.EQ.1.and.IC.eq.0)GOTO 666 if(ITU.LT.3)then write(3,*)'CRMS34:' write(3,*)'Only following atoms found in the big molecule:' write(3,*)(IND(iu) ,iu=1,ITU) if(LHALTONERROR)then CALL REPORT(' ITU < 3 !') else CALL SREPORT(' ITU < 3 !') return endif endif c c C Calculate center of the coordinates that will be rotated: DO 33 I=1,3 TB(I) =0.0d0 33 TS(I)=0.0d0 DO 34 II=1,ITU DO 34 J=1,3 TB(J)=TB(J)+XB(II,J) 34 TS(J)=TS(J)+XS(II,J) DO 38 J=1,3 TB(J)=TB(J)/DBLE(ITU) 38 TS(J)=TS(J)/DBLE(ITU) C Subtract mass-center: DO 35 II=1,ITU DO 35 J=1,3 XS(II,J)=XS(II,J)-TS(J) 35 XB(II,J)=XB(II,J)-TB(J) C IF(LINV)THEN DO 39 II=1,ITU DO 39 J=1,3 39 XS(II,J)=-XS(II,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR) IF(IERR.EQ.1.AND.IPASS.EQ.1)THEN write(3,*)'domat error,try other pass' GOTO 666 ENDIF IF(IERR.EQ.1)then if(LHALTONERROR)then CALL REPORT('Program cannot continue') else CALL SREPORT('Program cannot continue') return endif endif IERC=0 return end c SUBROUTINE vz0(N,cs) integer*4 i,N real*8 cs(*) do 1 i=1,N 1 cs(i)=0.0d0 return end SUBROUTINE READFC(fn,NATS,cs,ncs) integer*4 NATS,ncs,i,j,k,N character*(*) fn real*8 cs(*) open(20,file=fn) N=3*NATS ncs=0 call vz0(N**3,cs) c all constants supposed i=1..N,j=1..N,k=1..N 1 read(20,*,end=200,err=200)i,j,k,cs(i+N*(j-1+N*(k-1))) ncs=ncs+1 goto 1 200 close(20) return end c SUBROUTINE READ36(io,nat,db,ncs) implicit none integer*4 ncs,nat,io,i,j,k,l,ia,ix,ja,jx,ka,kx,la,lx real*8 db(*),x do 1 i=1,81*nat**3 1 db(i)=0.0d0 ncs=0 c supposed constants i<=j<=k<=l, from that we need ia ia ja ka c for all ia and ja<=ka c to make ia<=ja<=ka<=la, but all xyz indices: 2 read(io,*,end=200,err=200)i,j,k,l,x ncs=ncs+1 ia=(i+2)/3 ix=i-3*(ia-1) ja=(j+2)/3 jx=j-3*(ja-1) ka=(k+2)/3 kx=k-3*(ka-1) la=(l+2)/3 lx=l-3*(la-1) c use the L37 (three-atomic indexing) also for two-atom constants: c ============================================================ if(ja.eq.ia)then c II if(ka.eq.ia)then c III if(la.eq.ia)then c IIII* db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ia-4+kx))+3*ia-3+lx)=x db(nat*((9*ia+3*ix+kx-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+lx)=x db(nat*((9*ia+3*kx+ix-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+lx)=x db(nat*((9*ia+3*ix+kx-13)*nat*9+3*(3*ia-4+lx))+3*ia-3+jx)=x db(nat*((9*ia+3*kx+ix-13)*nat*9+3*(3*ia-4+lx))+3*ia-3+jx)=x db(nat*((9*ia+3*kx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+jx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ia-4+kx))+3*ia-3+lx)=x db(nat*((9*ia+3*jx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+lx)=x db(nat*((9*ia+3*kx+jx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+lx)=x db(nat*((9*ia+3*jx+kx-13)*nat*9+3*(3*ia-4+lx))+3*ia-3+ix)=x db(nat*((9*ia+3*kx+jx-13)*nat*9+3*(3*ia-4+lx))+3*ia-3+ix)=x db(nat*((9*ia+3*kx+lx-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+ix)=x db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ia-4+lx))+3*ia-3+kx)=x db(nat*((9*ia+3*ix+lx-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+kx)=x db(nat*((9*ia+3*lx+ix-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+kx)=x db(nat*((9*ia+3*ix+lx-13)*nat*9+3*(3*ia-4+kx))+3*ia-3+jx)=x db(nat*((9*ia+3*lx+ix-13)*nat*9+3*(3*ia-4+kx))+3*ia-3+jx)=x db(nat*((9*ia+3*lx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+jx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ia-4+lx))+3*ia-3+kx)=x db(nat*((9*ia+3*jx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+kx)=x db(nat*((9*ia+3*lx+jx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+kx)=x db(nat*((9*ia+3*jx+lx-13)*nat*9+3*(3*ia-4+kx))+3*ia-3+ix)=x db(nat*((9*ia+3*lx+jx-13)*nat*9+3*(3*ia-4+kx))+3*ia-3+ix)=x db(nat*((9*ia+3*lx+kx-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+ix)=x else c IIIL* = IILI db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ia-4+kx))+3*la-3+lx)=x db(nat*((9*ia+3*ix+kx-13)*nat*9+3*(3*ia-4+jx))+3*la-3+lx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ia-4+kx))+3*la-3+lx)=x db(nat*((9*ia+3*jx+kx-13)*nat*9+3*(3*ia-4+ix))+3*la-3+lx)=x db(nat*((9*ia+3*kx+ix-13)*nat*9+3*(3*ia-4+jx))+3*la-3+lx)=x db(nat*((9*ia+3*kx+jx-13)*nat*9+3*(3*ia-4+ix))+3*la-3+lx)=x db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*la-4+lx))+3*ia-3+kx)=x db(nat*((9*ia+3*ix+kx-13)*nat*9+3*(3*la-4+lx))+3*ia-3+jx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*la-4+lx))+3*ia-3+kx)=x db(nat*((9*ia+3*jx+kx-13)*nat*9+3*(3*la-4+lx))+3*ia-3+ix)=x db(nat*((9*ia+3*kx+ix-13)*nat*9+3*(3*la-4+lx))+3*ia-3+jx)=x db(nat*((9*ia+3*kx+jx-13)*nat*9+3*(3*la-4+lx))+3*ia-3+ix)=x endif else c IIK if(la.eq.ia)then c IIKI* =IIIK call report('forbidden combination IIKI') db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ka-4+kx))+3*ia-3+lx)=x db(nat*((9*ia+3*ix+lx-13)*nat*9+3*(3*ka-4+kx))+3*ia-3+jx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ka-4+kx))+3*ia-3+lx)=x db(nat*((9*ia+3*jx+lx-13)*nat*9+3*(3*ka-4+kx))+3*ia-3+ix)=x db(nat*((9*ia+3*lx+ix-13)*nat*9+3*(3*ka-4+kx))+3*ia-3+jx)=x db(nat*((9*ia+3*lx+jx-13)*nat*9+3*(3*ka-4+kx))+3*ia-3+ix)=x db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ia-4+lx))+3*ka-3+kx)=x db(nat*((9*ia+3*ix+lx-13)*nat*9+3*(3*ia-4+jx))+3*ka-3+kx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ia-4+lx))+3*ka-3+kx)=x db(nat*((9*ia+3*jx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ka-3+kx)=x db(nat*((9*ia+3*lx+ix-13)*nat*9+3*(3*ia-4+jx))+3*ka-3+kx)=x db(nat*((9*ia+3*lx+jx-13)*nat*9+3*(3*ia-4+ix))+3*ka-3+kx)=x else if(la.eq.ka)then c IIKK* =KKII db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ka-4+kx))+3*ka-3+lx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ka-4+kx))+3*ka-3+lx)=x db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ka-4+lx))+3*ka-3+kx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ka-4+lx))+3*ka-3+kx)=x db(nat*((9*ka+3*kx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+jx)=x db(nat*((9*ka+3*lx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ia-3+jx)=x db(nat*((9*ka+3*kx+lx-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+ix)=x db(nat*((9*ka+3*lx+kx-13)*nat*9+3*(3*ia-4+jx))+3*ia-3+ix)=x else c IIKL* =IILK db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*ka-4+kx))+3*la-3+lx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*ka-4+kx))+3*la-3+lx)=x db(nat*((9*ia+3*ix+jx-13)*nat*9+3*(3*la-4+lx))+3*ka-3+kx)=x db(nat*((9*ia+3*jx+ix-13)*nat*9+3*(3*la-4+lx))+3*ka-3+kx)=x endif endif endif c ============================================================ else c IJ if(ka.eq.ia)then c IJI if(la.eq.ia)then c IJII* =IIIJ=IIJI call report('forbidden combination IJII') db(nat*((9*ia+3*ix+kx-13)*nat*9+3*(3*ia-4+lx))+3*ja-3+jx)=x db(nat*((9*ia+3*ix+lx-13)*nat*9+3*(3*ia-4+kx))+3*ja-3+jx)=x db(nat*((9*ia+3*kx+ix-13)*nat*9+3*(3*ia-4+lx))+3*ja-3+jx)=x db(nat*((9*ia+3*kx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+jx)=x db(nat*((9*ia+3*lx+ix-13)*nat*9+3*(3*ia-4+kx))+3*ja-3+jx)=x db(nat*((9*ia+3*lx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+jx)=x db(nat*((9*ia+3*ix+kx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+lx)=x db(nat*((9*ia+3*ix+lx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+kx)=x db(nat*((9*ia+3*kx+ix-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+lx)=x db(nat*((9*ia+3*kx+lx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+ix)=x db(nat*((9*ia+3*lx+ix-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+kx)=x db(nat*((9*ia+3*lx+kx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+ix)=x else if(la.eq.ja)then c IJIJ* =IIJJ=JJII call report('forbidden combination IJIJ') else c IJIL* call report('forbidden combination IJIL') endif endif else if(ka.eq.ja)then c IJJ if(la.eq.ia)then c IJJI* call report('forbidden combination IJJI') else if(la.eq.ja)then c IJJJ* =JJJI=JJIJ db(nat*((9*ja+3*jx+kx-13)*nat*9+3*(3*ja-4+lx))+3*ia-3+ix)=x db(nat*((9*ja+3*jx+lx-13)*nat*9+3*(3*ja-4+kx))+3*ia-3+ix)=x db(nat*((9*ja+3*lx+jx-13)*nat*9+3*(3*ja-4+kx))+3*ia-3+ix)=x db(nat*((9*ja+3*lx+kx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+ix)=x db(nat*((9*ja+3*kx+lx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+ix)=x db(nat*((9*ja+3*kx+jx-13)*nat*9+3*(3*ja-4+lx))+3*ia-3+ix)=x db(nat*((9*ja+3*jx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+lx)=x db(nat*((9*ja+3*jx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+kx)=x db(nat*((9*ja+3*lx+jx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+kx)=x db(nat*((9*ja+3*lx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+jx)=x db(nat*((9*ja+3*kx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+jx)=x db(nat*((9*ja+3*kx+jx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+lx)=x else c IJJL* =JJIL=JJLI db(nat*((9*ja+3*jx+kx-13)*nat*9+3*(3*ia-4+ix))+3*la-3+lx)=x db(nat*((9*ja+3*kx+jx-13)*nat*9+3*(3*ia-4+ix))+3*la-3+lx)=x db(nat*((9*ja+3*jx+kx-13)*nat*9+3*(3*la-4+lx))+3*ia-3+ix)=x db(nat*((9*ja+3*kx+jx-13)*nat*9+3*(3*la-4+lx))+3*ia-3+ix)=x endif endif else c IJK if(la.eq.ia)then c IJKI* call report('forbidden combination IJKI') else if(la.eq.ja)then c IJKJ* call report('forbidden combination IJKJ') else if(la.eq.ka)then c IJKK* =KKIJ=KKJI db(nat*((9*ka+3*lx+kx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+jx)=x db(nat*((9*ka+3*kx+lx-13)*nat*9+3*(3*ia-4+ix))+3*ja-3+jx)=x db(nat*((9*ka+3*lx+kx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+ix)=x db(nat*((9*ka+3*kx+lx-13)*nat*9+3*(3*ja-4+jx))+3*ia-3+ix)=x else c IJKL* write(6,*)ncs,ia,ja,ka,la,ix,jx,kx,lx write(6,*)i,j,k,l,x call report('forbidden combination IJKL') endif endif endif endif endif endif c ============================================================ goto 2 200 continue c 1112 1113 c write(6,*)db(nat*((9+3+1-13)*nat*9+3*(3-4+1))+3-3+2) c write(6,*)db(nat*((9+3+1-13)*nat*9+3*(3-4+1))+3-3+3) c stop return end c SUBROUTINE WRITE36(io,nat,d,ALIM,ncs,L36) implicit none integer*4 ncs,nat,io,i,j,k,l,ia,ix,ja,jx,ka,kx,la,lx,i1,N real*8 d(*),ALIM logical L36 c use the L37 (three-atomic) indexing also for L36 (two-atomic): N=3*nat ncs=0 do 851 i=1,N ia=(i+2)/3 ix=i-3*(ia-1) do 851 j=i,N ja=(j+2)/3 jx=j-3*(ja-1) do 851 k=j,N ka=(k+2)/3 kx=k-3*(ka-1) do 851 l=k,N la=(l+2)/3 lx=l-3*(la-1) if(ia.eq.ja)then i1=nat*((9*ia+3*ix+jx-13)*nat*9+3*(k-1))+l if(L36.and.ka.ne.ia.and.la.ne.ia.and.ka.ne.la)goto 851 else if(ia.eq.ka)then i1=nat*((9*ia+3*ix+kx-13)*nat*9+3*(j-1))+l if(L36.and.ja.ne.ia.and.la.ne.ia.and.ja.ne.la)goto 851 else if(ia.eq.la)then i1=nat*((9*ia+3*ix+lx-13)*nat*9+3*(k-1))+j if(L36.and.ja.ne.ia.and.ka.ne.ia.and.ja.ne.ka)goto 851 else if(ja.eq.ka)then i1=nat*((9*ja+3*jx+kx-13)*nat*9+3*(i-1))+l if(L36.and.ja.ne.ia.and.ja.ne.la.and.ia.ne.la)goto 851 else if(ja.eq.la)then i1=nat*((9*ja+3*jx+lx-13)*nat*9+3*(i-1))+k if(L36.and.ja.ne.ia.and.ja.ne.ka.and.ia.ne.ka)goto 851 else if(ka.eq.la)then i1=nat*((9*ka+3*kx+lx-13)*nat*9+3*(i-1))+j if(L36.and.ka.ne.ia.and.ka.ne.ja.and.ia.ne.ja)goto 851 else goto 851 endif endif endif endif endif endif IF(DABS(d(i1)).GE.ALIM)THEN ncs=ncs+1 WRITE(io,400)i,j,k,l,d(i1) 400 format(4i6,G16.8) ENDIF 851 CONTINUE return end c ============================================================ SUBROUTINE WRITEFC(io,NATS,cs,ALIM,ncs) implicit none integer*4 NATS,io,i,j,k,ncs,N real*8 cs(*),ALIM ncs=0 N=3*NATS do 1 i=1,N do 1 j=1,N do 1 k=1,N if(dabs(cs(i+N*(j-1+N*(k-1)))).ge.ALIM)then ncs=ncs+1 write(io,100)i,j,k,cs(i+N*(j-1+N*(k-1))) 100 format(3i6,G16.8) endif 1 continue close(io) return end c ============================================================ SUBROUTINE readpol(a,fn) character*(*)fn real*8 a(3,3) integer*4 i,j open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(i,j),i=1,j),j=1,3) close(90) do 1 j=1,3 do 1 i=1,j 1 a(j,i)=a(i,j) return end c ============================================================ SUBROUTINE readpolg(a,fn) character*(*)fn real*8 a(3,3) integer*4 i,j open(90,file=fn) read(90,*) read(90,*) read(90,*)((a(j,i),i=1,3),j=1,3) close(90) return end c ============================================================ SUBROUTINE readpola(a,fn) character*(*)fn real*8 a(3,3,3) integer*4 i,j,k open(90,file=fn) read(90,*) read(90,*) do 1 k=1,3 read(90,*)((a(i,j,k),i=1,j),j=1,3) do 1 j=1,3 do 1 i=1,j 1 a(j,i,k)=a(i,j,k) close(90) return end c ============================================================ SUBROUTINE CRMSo(NAT,INDBIG,IOX,MX,IPO,LNAMES,natname,RS,RB,IND, 2rname,NATS,LINV) c make the rotation matrix for all atoms in the overlap IPO implicit none integer*4 IOX,MXB,MX PARAMETER (MXB=140) integer*4 II,NAT,INDBIG(IOX,MX),NATS,natname(*),IX,IND(*),I,IERR, 2IPO,J,ITU,IIT real*8 RMS,RS(3,MX),RB(3,MX),TB(3),TS(3),rname(IOX,3,MX), 1A,XS,XB,TOL logical LNAMES,LINV COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS c if(LNAMES)NATS=natname(IPO) ITU=0 DO 27 II=1,NAT if(INDBIG(IPO,II).NE.0.and.ITU.lt.MXB)then C Include atom II into the transformation set: ITU=ITU+1 IF(ITU.EQ.MXB)CALL w36('Warning: MXB reached in CRMSo') DO 28 IX=1,3 if(LNAMES)then XS(ITU,IX)=rname(IPO,IX,ABS(INDBIG(IPO,II))) else XS(ITU,IX)=RS(IX,ABS(INDBIG(IPO,II))) endif 28 XB(ITU,IX)=RB(IX,II) IND(ITU)=II endif 27 CONTINUE c if(ITU.LT.3) CALL REPORT(' ITU < 3 !') c C Calculate center of the coordinates that will be rotated: DO 33 I=1,3 TB(I) =0.0d0 33 TS(I)=0.0d0 DO 34 II=1,ITU DO 34 J=1,3 TB(J)=TB(J)+XB(II,J) 34 TS(J)=TS(J)+XS(II,J) DO 38 J=1,3 TB(J)=TB(J)/DBLE(ITU) 38 TS(J)=TS(J)/DBLE(ITU) C Subtract mass-centers: DO 35 II=1,ITU DO 35 J=1,3 XS(II,J)=XS(II,J)-TS(J) 35 XB(II,J)=XB(II,J)-TB(J) C IF(LINV)THEN DO 39 II=1,ITU DO 39 J=1,3 39 XS(II,J)=-XS(II,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR) IF(IERR.EQ.1)CALL REPORT('CRMs DoMatrix error') return end c ============================================================ subroutine w36(s) character*(*) s write(6,*)s write(3,*)s return end c ============================================================ subroutine ct(ta,t,o) integer*4 i,j real*8 ta(3,3),t(3,3),o(3,3) do 1 i=1,3 do 1 j=1,3 1 ta(i,j)=t(i,1)*o(1,j)+t(i,2)*o(2,j)+t(i,3)*o(3,j) return end c ============================================================ subroutine ctd(IIND,ALPHA,tem,M) integer*4 IIND,M real*8 ALPHA(M,3,3),tem(3,3) integer*4 i,j do 1 i=1,3 do 1 j=1,3 1 tem(i,j)=ALPHA(IIND,i,j) return end c ============================================================ subroutine ctdd(IIND,A,tem,M) integer*4 IIND,M real*8 A(M,3,3,3),tem(3,3,3) integer*4 i,j,k do 1 i=1,3 do 1 j=1,3 do 1 k=1,3 1 tem(i,j,k)=A(IIND,i,j,k) return end c ============================================================ subroutine addtoG(IIND,G,M,tc) integer*4 IIND,M,I,J real*8 G(M,3,3),tc(3,3) do 1 I=1,3 do 1 J=1,3 1 G(IIND,I,J)=G(IIND,I,J)+tc(I,J) return end c ============================================================ SUBROUTINE WRITEPOL(A,f) real*8 A(3,3) character*(*) f integer*4 IX,IY OPEN(90,FILE=f) if(f.eq.'POL.TTT')then write(90,900)((A(IY,IX),IY=1,IX),IX=1,3) 900 format('Polarizability:',/, 1 ' XX XY YY', 2 ' XZ YZ ZZ',/,6G14.6) endif if(f.eq.'POL.G.TTT')then write(90,901)((A(IX,IY),IY=1,3),IX=1,3) 901 format('G-tensor, first index-electric, second-magnetic:',/, 1 ' XX XY XZ', 2 ' YX YY YZ', 2 ' ZX ZY ZZ',/,9G14.6) endif close(90) WRITE(*,*)f//' written' return end c ============================================================ SUBROUTINE WRITEPOLA(A) real*8 A(3,3,3) integer*4 IX,IY,IZ OPEN(90,FILE='POL.A.TTT') write(90,900) 900 format('A-tensor:',/, 1 ' XX XY YY', 2 ' XZ YZ ZZ') do 1 IZ=1,3 1 write(90,901)((A(IY,IX,IZ),IY=1,IX),IX=1,3) 901 format(6G14.6) close(90) WRITE(*,*)'POL.A.TTT written' return end c ============================================================ subroutine rpp(IOX,IO,gopol,gpol,A) implicit none integer*4 IX,JX,IXP,JXP,IOX,IO real*8 gopol(IOX,3,3),A(3,3),gpol(3,3),t(3,3) do 1 IX=1,3 do 1 JXP=1,3 t(IX,JXP)=0.0d0 do 1 IXP=1,3 1 t(IX,JXP)=t(IX,JXP)+A(IX,IXP)*gpol(IXP,JXP) do 2 IX=1,3 do 2 JX=1,3 gopol(IO,IX,JX)=0.0d0 do 2 JXP=1,3 2 gopol(IO,IX,JX)=gopol(IO,IX,JX)+A(JX,JXP)*t(IX,JXP) return end c ============================================================ subroutine rppp(NO,IO,gopol,gpol,A) implicit none integer*4 IX,JX,JXP,NO,IO,KXP,KX,it real*8 gopol(NO,3,3,3),A(3,3),gpol(3,3,3),t(3,3,3) it=IO do 1 IX=1,3 do 1 JXP=1,3 do 1 KXP=1,3 1 t(IX,JXP,KXP)= A(IX,1)*gpol(1,JXP,KXP)+ 1A(IX,2)*gpol(2,JXP,KXP)+A(IX,3)*gpol(3,JXP,KXP) do 2 IX=1,3 do 2 JX=1,3 do 2 KXP=1,3 2 gopol(it,IX,JX,KXP)=A(JX,1)*t(IX,1,KXP)+ 1A(JX,2)*t(IX,2,KXP)+A(JX,3)*t(IX,3,KXP) call ctdd(it,gopol,t,NO) do 3 IX=1,3 do 3 JX=1,3 do 3 KX=1,3 3 gopol(it,IX,JX,KX)=A(KX,1)*t(IX,JX,1)+ 1A(KX,2)*t(IX,JX,2)+A(KX,3)*t(IX,JX,3) return end c ============================================================ subroutine inipol(p,g,a) implicit none real*8 p(3,3),g(3,3),a(3,3,3) integer*4 i,j,k do 1 i=1,3 do 1 j=1,3 p(i,j)=0.0d0 g(i,j)=0.0d0 do 1 k=1,3 1 a(i,j,k)=0.0d0 return end c ============================================================ subroutine tadd(bigpol,tc) implicit none c B=B+C real*8 bigpol(3,3),tc(3,3) integer*4 ix,jx do 615 ix=1,3 do 615 jx=1,3 615 bigpol(ix,jx)=bigpol(ix,jx)+tc(ix,jx) return end c ============================================================ subroutine mkttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3),rt2,rt5 integer*4 ix,jx,kx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 ttt(ix,jx,kx)=-15.0d0*rt(ix)*rt(jx)*rt(kx)/rt5/rt2 if(ix.eq.jx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(kx)/rt5 if(ix.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(jx)/rt5 613 if(jx.eq.kx)ttt(ix,jx,kx)=ttt(ix,jx,kx)+3.0d0*rt(ix)/rt5 return end c ============================================================ subroutine mktt(rt,tt) implicit none real*8 rt(3),tt(3,3),rt2,rt5 integer*4 ix,jx rt2=rt(1)**2+rt(2)**2+rt(3)**2 rt5=rt2**2*dsqrt(rt2) do 611 ix=1,3 do 611 jx=1,3 611 tt(ix,jx)=3.0d0*rt(ix)*rt(jx)/rt5 do 610 ix=1,3 610 tt(ix,ix)=tt(ix,ix)-rt2/rt5 return end c ============================================================ subroutine dera(era,rt,al) implicit none real*8 era(3,3),rt(3),al(3,3) integer*4 IX DO 159 IX=1,3 era(1,IX)=rt(2)*al(3,IX)-rt(3)*al(2,IX) era(2,IX)=rt(3)*al(1,IX)-rt(1)*al(3,IX) 159 era(3,IX)=rt(1)*al(2,IX)-rt(2)*al(1,IX) return end c ============================================================ subroutine multt(tc,f) implicit none real*8 tc(3,3),f integer*4 I,J do 605 I=1,3 do 605 J=1,3 605 tc(I,J)=tc(I,J)*f return end c ============================================================ subroutine wrt2(IA,MX,ALPHA,s) implicit none integer*4 MX,IIND,IA,IX,I,J real*8 ALPHA(3*MX,3,3) character*(*) s write(3,*)s do 611 IX=1,3 IIND=3*(IA-1)+IX 611 write(3,3007)((ALPHA(IIND,I,J),I=1,3),J=1,3) 3007 format(9f8.3) return end c ============================================================ subroutine wrt3(IA,MX,ALPHA,s) implicit none integer*4 MX,IIND,IA,IX,I,J,K real*8 ALPHA(3*MX,3,3,3) character*(*) s write(3,*)s do 611 IX=1,3 IIND=3*(IA-1)+IX 611 write(3,3007)(((ALPHA(IIND,I,J,K),I=1,3),J=1,3),K=1,3) 3007 format(9f8.3) return end c ============================================================ subroutine std(a) implicit none real*8 a(3,3),t(3,3) integer*4 i,j do 1 i=1,3 do 1 j=1,3 1 t(i,j)=a(i,j)+a(j,i) do 2 i=1,3 do 2 j=1,3 2 a(i,j)=t(i,j) return end c ============================================================ subroutine readsos(a,n,io,np) implicit none integer*4 n,k,io,i,ii,j,l,np real*8 ar(3,3),ai(3,3),a(*),arlv(3,3),ailv(3,3) read(44,*) read(44,*) read(44,*) c dip dip polarizability length-length and length-velocity: do 1 k=1,n read(44,*)ar(1,1),(ar( 1,i),i=1,3),(ai( 1,i),i=1,3), 1 (arlv(1,i),i=1,3),(ailv(1,i),i=1,3) read(44,*) (ar( 2,i),i=1,3),(ai( 2,i),i=1,3), 1 (arlv(2,i),i=1,3),(ailv(2,i),i=1,3) read(44,*) (ar( 3,i),i=1,3),(ai( 3,i),i=1,3), 1 (arlv(3,i),i=1,3),(ailv(3,i),i=1,3) ii=0 do 1 i=1,3 do 1 j=1,3 ii=ii+1 a(np*((io-1)*n+k-1)+ii )=ar(i,j) a(np*((io-1)*n+k-1)+ii+9)=ai(i,j) a(np*((io-1)*n+k-1)+ii+90 )=arlv(i,j) 1 a(np*((io-1)*n+k-1)+ii+99)=ailv(i,j) c 1.. 9 real alpha length-length c 10..18 im alpha length-length c 91..99 real alpha length-velocity c 100.118 im alpha length-velocity read(44,*) read(44,*) do 2 k=1,n c dip mag polarizability: read(44,*)ar(1,1),(ar(1,i),i=1,3),(ai(1,i),i=1,3) read(44,*) (ar(2,i),i=1,3),(ai(2,i),i=1,3) read(44,*) (ar(3,i),i=1,3),(ai(3,i),i=1,3) ii=0 do 2 i=1,3 do 2 j=1,3 ii=ii+1 a(np*((io-1)*n+k-1)+ii+18)=ar(i,j) 2 a(np*((io-1)*n+k-1)+ii+27)=ai(i,j) c 19..27 real G' c 28..36 im G' read(44,*) read(44,*) do 3 k=1,n ii=0 do 3 l=1,3 c dip quad polarizability: if(l.eq.1)then read(44,*)ar(1,1),(ar(1,i),i=1,3),(ai(1,i),i=1,3) else read(44,*) (ar(1,i),i=1,3),(ai(1,i),i=1,3) endif read(44,*) (ar(2,i),i=1,3),(ai(2,i),i=1,3) read(44,*) (ar(3,i),i=1,3),(ai(3,i),i=1,3) do 3 i=1,3 do 3 j=1,3 ii=ii+1 a(np*((io-1)*n+k-1)+ii+36)=ar(i,j) 3 a(np*((io-1)*n+k-1)+ii+63)=ai(i,j) c 37..63 real A c 64..90 im A close(44) return end c ============================================================ subroutine writesos(a,n,wmin,wmax,nel,SOSNM) implicit none integer*4 n,k,i,j,ii,l,nel real*8 a(*),wmin,wmax,dw,w,ar(3,3),ai(3,3),cab,ccd,wau, 1pi,debye,arlv(3,3),ailv(3,3),clight logical SOSNM c clight=137.036d0 dw=(wmax-wmin)/dble(n-1) pi=4.0d0*atan(1.0d0) debye=2.541765d0 cab=debye**2 *108.7d0/pi ccd=debye**2*4.0d0*108.7d0/pi/clight open(44,file='SOS.TTT') open(9,file='spa.prn') open(91,file='sd.prn') write(44,446)n,wmin,wmax 446 format(i6,2f12.2) write(44,443) 443 format(' Dynamic polarizabilities') write(44,444) 444 format(10x,' Alpha ( f ) Alpha ( g) length-length ', 1'and length-velocity : ', 1/,'l(nm) Re ullx Re ully Re ullz Im ullx Im ully', 1' Im ullz Re ulvx Re ulvy Re ulvz Im ulvx Im ulvy', 1' Im ulvz') c w=wmin-dw do 1 k=1,n w=w+dw if(SOSNM)then wau=(1.0d7/w)/219474.63d0 else wau=w/219474.63d0 endif ii=0 do 11 i=1,3 do 11 j=1,3 ii=ii+1 ar( i,j)=a(nel*(k-1)+ii ) ai( i,j)=a(nel*(k-1)+ii+ 9) arlv(i,j)=a(nel*(k-1)+ii+90) 11 ailv(i,j)=a(nel*(k-1)+ii+99) write(9,69)w,ai(1,1)+ai(2,2)+ai(3,3),ar(1,1)+ar(2,2)+ar(3,3) c absorption spectrum only: write(91,69)w,(ai(1,1)+ai(2,2)+ai(3,3))*cab*wau**2 69 format(f7.1,2e14.4) write(44,441)w,(ar( 1,j),j=1,3),(ai( 1,j),j=1,3), 1 (arlv(1,j),j=1,3),(ailv(1,j),j=1,3) write(44,442) (ar( 2,j),j=1,3),(ai( 2,j),j=1,3), 1 (arlv(2,j),j=1,3),(ailv(2,j),j=1,3) 1 write(44,442) (ar( 3,j),j=1,3),(ai( 3,j),j=1,3), 1 (arlv(3,j),j=1,3),(ailv(3,j),j=1,3) 441 format(f8.2,12f12.4) 442 format(8x ,12f12.4) close(9) close(91) write(6,*)' spa.prn, sd.prn written' write(44,445) 445 format(10x,' Gp/w x f Gp/w x g ', 1/,'l(nm) ux uy uz') open(9,file='spg.prn') open(91,file='sr.prn') w=wmin-dw do 2 k=1,n w=w+dw if(SOSNM)then wau=(1.0d7/w)/219474.63d0 else wau=w/219474.63d0 endif ii=0 do 22 i=1,3 do 22 j=1,3 ii=ii+1 ar(i,j)=a(nel*(k-1)+ii+18) 22 ai(i,j)=a(nel*(k-1)+ii+27) write(9,69)w,ai(1,1)+ai(2,2)+ai(3,3),ar(1,1)+ar(2,2)+ar(3,3) write(91,69)w,-(ai(1,1)+ai(2,2)+ai(3,3))*ccd*wau**3 write(44,441)w,(ar(1,j),j=1,3),(ai(1,j),j=1,3) write(44,442) (ar(2,j),j=1,3),(ai(2,j),j=1,3) 2 write(44,442) (ar(3,j),j=1,3),(ai(3,j),j=1,3) close(9) close(91) write(6,*)' spg.prn sr.prn written' write(44,447) 447 format(10x,' A x f A x g ', 1/,'l(nm) ux uy uz') w=wmin-dw do 3 k=1,n w=w+dw ii=0 do 3 i=1,3 do 33 j=1,3 do 33 l=1,3 ii=ii+1 ar(j,l)=a(nel*(k-1)+ii+36) 33 ai(j,l)=a(nel*(k-1)+ii+63) if(i.eq.1)then write(44,441)w,(ar(1,j),j=1,3),(ai(1,j),j=1,3) else write(44,442) (ar(1,j),j=1,3),(ai(1,j),j=1,3) endif write(44,442) (ar(2,j),j=1,3),(ai(2,j),j=1,3) 3 write(44,442) (ar(3,j),j=1,3),(ai(3,j),j=1,3) c 1.. 9 real alpha c 10..18 im alpha c 19..27 real G' c 28..36 im G' c 37..63 real A c 64..90 im A c 91..99 real alpha length-velocity c 100.118 im alpha length-velocity close(44) write(6,*)' SOS.TTT written' return end c ============================================================ subroutine trdpol(np,alsos,coo,io,nel) implicit none integer*4 np,iw,ic,ioa,iog,ioaa,ii,ix,jx,kx,io,nel,ioalv real*8 po(3,3),alsos(*),gpo(3,3),apo(3,3,3),coo(*),polv(3,3) c for each frequency: do 406 iw=1,np c for real and imaginary parts: do 406 ic=1,2 if(ic.eq.1)then ioa=0 ioalv=90 iog=18 ioaa=36 else ioa=9 ioalv=99 iog=27 ioaa=63 endif c rewrite al, G, and A to local variables ii=0 do 407 ix=1,3 do 407 jx=1,3 ii=ii+1 po( ix,jx) =alsos(nel*((io-1)*np+iw-1)+ii+ioa) polv(ix,jx) =alsos(nel*((io-1)*np+iw-1)+ii+ioalv) 407 gpo( ix,jx) =alsos(nel*((io-1)*np+iw-1)+ii+iog) ii=0 do 4071 ix=1,3 do 4071 jx=1,3 do 4071 kx=1,3 ii=ii+1 c apo.. last index electric: 4071 apo(jx,kx,ix)=alsos(nel*((io-1)*np+iw-1)+ii+ioaa) c transform G, to small molecule center: call TTTG(polv,gpo,1,coo,0) c transform A, to small molecule center: call TTTA(po,apo,1,coo,0) c rewrite back: ii=0 do 408 ix=1,3 do 408 jx=1,3 ii=ii+1 alsos(nel*((io-1)*np+iw-1)+ii+ioa)= po(ix,jx) 408 alsos(nel*((io-1)*np+iw-1)+ii+iog)=gpo(ix,jx) ii=0 do 406 ix=1,3 do 406 jx=1,3 do 406 kx=1,3 ii=ii+1 406 alsos(nel*((io-1)*np+iw-1)+ii+ioaa)=apo(jx,kx,ix) return end c ============================================================ subroutine getname(NAME,basicname,NE) implicit none integer*4 NE character*80 NAME,basicname NE=80 do 10 NE=80,1,-1 10 if(NAME(NE:NE).ne.' ')goto 12 12 basicname(1:NE)=NAME(1:NE) return end c ============================================================ subroutine MDP(r,x,Pr,Pi,n,alphar,alphai,gpr,gpi, 1apr,api) implicit none integer*4 n,i,j,k,l,im,ixx,m,kl,ik real*8 x(24*n,24*n),pr(24*n,24*n), 1alphar(3,3),gpr(3,3),apr(3,3,3),Pi(24*n,24*n), 1alphai(3,3),gpi(3,3),api(3,3,3), 3r(n,3),sum,eps,rr(3),delta real*8, allocatable::ai(:,:),ar(:,:),aa(:,:),aai(:,:), 1itemi(:,:),rtemi(:,:),tem(:,:),temi(:,:),ptr(:,:),pti(:,:) c allocate(ar(n*24,n*24),ai(n*24,n*24),aa(n*24,n*24),aai(n*24,n*24), 1rtemi(n*24,n*24),itemi(n*24,n*24),tem(n*24,n*24),temi(n*24,n*24), 1ptr(n*24,n*24),pti(n*24,n*24)) c c make (E-X.P),= ar+i ai =A, real and imaginary parts call mulmg(x,pr,ar,n*24) call negative(ar,24*n) do 991 i=1,n*24 991 ar(i,i)=1.0d0+ar(i,i) call mulmg(x,pi,ai,n*24) call negative(ai,24*n) c c AA=ReA ReA + ImA ImA: call mulmg(ar,ar,tem,n*24) call mulmg(ai,ai,temi,n*24) call summ(tem,temi,aa,n*24) c AA^(-1)=aai call INV(aa,aai,24*n,1.0d-25) c make A^(-1)=(E-X.P)^(-1) = rtemi + i itemi = (ReA-iImA) x aai call mulmg(ar,aai,rtemi,n*24) call mulmg(ai,aai,itemi,n*24) call negative(itemi,24*n) c make total polarizability Pt = P(E-X.P)^(-1)=P temi c =(Re P Re temi - Im P Im temi)+i(Re P Im temi + Im P Re temi) c local origins call mulmg(pr,rtemi,tem,n*24) call mulmg(pi,itemi,temi,n*24) call negative(temi,n*24) call summ(tem,temi,ptr,n*24) call mulmg(pr,itemi,tem,n*24) call mulmg(pi,rtemi,temi,n*24) call summ(tem,temi,pti,n*24) c make total polarizability, common origin: do 900 j=1,3 do 900 k=1,3 alphar(j,k)=0.0d0 alphai(j,k)=0.0d0 gpr(j,k)=0.0d0 gpi(j,k)=0.0d0 do 900 l=1,3 api(j,k,l)=0.0d0 900 apr(j,k,l)=0.0d0 c do 8 im=1,n do 84 ixx=1,3 84 rr(ixx)=r(im,ixx) do 8 j=1,3 do 8 ik=1,n kl=0 do 8 k=1,3 c from electric moment 1..3: c from electric field 1..3: alphar(j,k)=alphar(j,k)+ptr(24*(im-1)+j,24*(ik-1)+k) alphai(j,k)=alphai(j,k)+pti(24*(im-1)+j,24*(ik-1)+k) c from magnetic moment 7-9: c from electric field derivative 4..6: sum=0.0d0 do 85 l=1,3 do 85 m=1,3 85 sum=sum+eps(k,l,m)*rr(l)*ptr(24*(im-1)+3+m,24*(ik-1)+j+3)/2.0d0 gpr(j,k)=gpr(j,k)-ptr(24*(im-1)+6+k,24*(ik-1)+3+j)-sum sum=0.0d0 do 815 l=1,3 do 815 m=1,3 815 sum=sum+eps(k,l,m)*rr(l)*pti(24*(im-1)+3+m,24*(ik-1)+j+3)/2.0d0 gpi(j,k)=gpi(j,k)-pti(24*(im-1)+6+k,24*(ik-1)+3+j)-sum do 8 l=k,3 kl=kl+1 c from quadrupole 13-18 c from electric field 1..3: c a 0 0 G A/3 0 c 0 a -G 0 0 A/3 c 0 -Gt 0 0 0 0 c Gt 0 0 0 0 0 c A/3 0 0 0 0 0 c 0 A/3 0 0 0 0 api(j,k,l)=api(j,k,l)+pti(24*(im-1)+kl+12,24*(ik-1)+j)*3.0d0 1+1.5d0 *(rr(l)*pti(24*(im-1)+k,24*(ik-1)+j) 1 +rr(k)*pti(24*(im-1)+l,24*(ik-1)+j)) 1-delta(k,l)*(rr(1)*pti(24*(im-1)+1,24*(ik-1)+j) 1+ rr(2)*pti(24*(im-1)+2,24*(ik-1)+j) 1+ rr(3)*pti(24*(im-1)+3,24*(ik-1)+j)) api(j,l,k)=api(j,k,l) apr(j,k,l)=apr(j,k,l)+ptr(24*(im-1)+kl+12,24*(ik-1)+j)*3.0d0 1+1.5d0 *(rr(l)*ptr(24*(im-1)+k,24*(ik-1)+j) 1 +rr(k)*ptr(24*(im-1)+l,24*(ik-1)+j)) 1-delta(k,l)*(rr(1)*ptr(24*(im-1)+1,24*(ik-1)+j) 1+ rr(2)*ptr(24*(im-1)+2,24*(ik-1)+j) 1+ rr(3)*ptr(24*(im-1)+3,24*(ik-1)+j)) 8 apr(j,l,k)=apr(j,k,l) return end c ============================================================ subroutine negative(t,n) implicit none integer*4 i,n,j real*8 t(n,n) do 1 i=1,n do 1 j=1,n 1 t(i,j)=-t(i,j) return end c ============================================================ subroutine summ(a,b,c,n) implicit none integer*4 i,j,n real*8 a(n,n),b(n,n),c(n,n) do 1 i=1,n do 1 j=1,n 1 c(i,j)=a(i,j)+b(i,j) return end c ============================================================ subroutine doX(n,r,X) implicit none integer*4 i,n,istart,j,jstart,ix,jx,kx,lx,ii,jj real*8 rij(3),r(n,3),tt(3,3),ttt(3,3,3),tttt(3,3,3,3), 1X(n*24,n*24) do 6 i=1,n istart=24*(i-1)+1 do 6 j=1,n jstart=24*(j-1)+1 if(i.ne.j)then do 61 ix=1,3 61 rij(ix)=r(i,ix)-r(j,ix) c T 0 0 0 -U 0 T=grad(1/rij), U=grad(T), G=grad(U), W=T/c c 0 T 0 0 0 -U c 0 0 W 0 0 0 c 0 0 0 W 0 0 c U 0 0 0 -G 0 c 0 U 0 0 0 -G call mktt(rij,tt) do 62 ix=1,3 do 62 jx=1,3 X(istart+ix-1, jstart+jx-1 )=tt(ix,jx) X(istart+ix-1+3,jstart+jx-1+3)=tt(ix,jx) X(istart+ix-1+6,jstart+jx-1+6)=tt(ix,jx)/137.0d0 62 X(istart+ix-1+9,jstart+jx-1+9)=tt(ix,jx)/137.0d0 call mkttt(rij,ttt) do 65 ix=1,3 ii=0 do 65 jx=1,3 do 65 kx=jx,3 ii=ii+1 X(istart+ix-1,jstart+11+ii)=-ttt(ix,jx,kx) X(istart+11+ii,jstart+ix-1)= ttt(ix,jx,kx) X(istart+ix+2,jstart+17+ii)=-ttt(ix,jx,kx) 65 X(istart+17+ii,jstart+ix+2)= ttt(ix,jx,kx) call mktttt(rij,tttt) ii=0 do 64 ix=1,3 do 64 jx=ix,3 ii=ii+1 jj=0 do 64 kx=1,3 do 64 lx=kx,3 jj=jj+1 X(istart+ii+11,jstart+11+jj)=-tttt(ix,jx,kx,lx) 64 X(istart+ii+17,jstart+17+jj)=-tttt(ix,jx,kx,lx) endif 6 continue return end c ============================================================ function delta(i,j) implicit none real*8 delta integer*4 i,j if(i.eq.j)then delta=1.0d0 else delta=0.0d0 endif return end c ============================================================ subroutine mulmg(a,b,c,n) implicit none integer*4 n,i,j,k real*8 a(n,n),b(n,n),c(n,n) do 1 i=1,n do 1 j=1,n c(i,j)=0.0d0 do 1 k=1,n 1 c(i,j)=c(i,j)+a(i,k)*b(k,j) return end c ============================================================ SUBROUTINE INV(A,AI,N,TOL) IMPLICIT none integer*4 N,oo,ii,jj,iw,io,kk,i2,j2 REAL*8 TOL, A(N,N),AI(N,N),w real*8, allocatable ::E(:,:) allocate(E(N,2*N)) DO 1 ii=1,N DO 1 jj=1,N e(ii,jj)=a(ii,jj) E(II,JJ+N)=0.0D0 1 if (ii.EQ.jj)e(ii,jj+N)=1.0D0 DO 2 ii=1,N-1 iw=ii if (DABS(e(iw,iw)).LT.TOL) then DO 3 io=iw+1,N oo=IO if (io.GT.N)oo=io-N 3 if (DABS(e(oo,iw)).GE.TOL) goto 11 call report('Inverse cannot be done') 11 CONTINUE DO 4 kk=1, 2*N w=e(iw,kk) e(iw,kk)=e(oo,kk) 4 e(oo,kk)=w ENDIF DO 2 jj=ii+1,N DO 6 kk=ii+1, 2*N 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii) 2 e(jj,ii)=0.0D0 DO 7 ii=1, N-1 i2=N-ii+1 DO 7 jj=1,i2-1 j2=i2-jj DO 9 kk=1, N 9 e(j2,kk+N)=e(j2,kk+N)-e(i2,kk+N)*e(j2,i2)/e(i2,i2) 7 e(j2,i2)=0.0D0 DO 10 ii=1,N DO 10 jj=1,N 10 AI(ii,jj)=e(ii,jj+N)/e(ii,ii) deallocate(E) RETURN END c ============================================================ subroutine mktttt(rt,ttt) implicit none real*8 rt(3),ttt(3,3,3,3),r2,r4,r5,r9,delta integer*4 ix,jx,kx,lx r2=rt(1)**2+rt(2)**2+rt(3)**2 r4=r2*r2 r5=r4*dsqrt(r2) r9=r5*r4 do 613 ix=1,3 do 613 jx=1,3 do 613 kx=1,3 do 613 lx=1,3 613 ttt(ix,jx,kx,lx)=3.0d0*(35.0d0*rt(ix)*rt(jx)*rt(kx)*rt(lx) 1-5.0d0*r2*(delta(ix,kx)*rt(jx)*rt(lx) 1 +delta(ix,lx)*rt(jx)*rt(kx) 1 +delta(ix,jx)*rt(kx)*rt(lx) 1 +delta(kx,jx)*rt(lx)*rt(ix) 1 +delta(kx,lx)*rt(ix)*rt(jx) 1 +delta(jx,lx)*rt(ix)*rt(kx)) 1+r4*(delta(ix,lx)*delta(kx,jx) 1 +delta(kx,lx)*delta(ix,jx) 1 +delta(jx,lx)*delta(ix,kx)))/r9 return end c ============================================================ subroutine rweights(NO,O,iweights,weights,NATS_max,NATS,rw,np) implicit none real*8 weights(*),w integer*4 iweights,i,NATS,O,NATS_max,np,nf,NO,ii,iw character*(*) rw if(iweights.gt.0)then c weights constant (for one frequency): nf=1 else c weights changing with frequency nf=np endif open(55,file=rw,status='old') do 2 iw=1,nf ii=(iw-1)*NO*NATS_max+NATS_max*(O-1) c order of weights c 1 real c 2 imaginary c 3 abs c 4 real trimmed c 5 imaginary trimmed c 6 abs trimmed if(iweights.le.0)then read(55,*)w,w backspace 55 else w=0.0d0 endif do 1 i=1,abs(iweights)-1 read(55,*) 1 read(55,*) read(55,*) read(55,*)(weights(ii+i),i=1,NATS) do 3 i=abs(iweights)+1,6 read(55,*) 3 read(55,*) do 21 i=1,NATS 21 weights(ii+i)=weights(ii+i)/100.0d0 write(6,6001)iw,w 6001 format(i6,f12.3) 2 write(6,6002)(weights(ii+i),i=1,NATS) 6002 format(10f7.3) close(55) return end c ============================================================ SUBROUTINE IMPROVESOSw(np,INDBIG,IOX,MX,IO,LNAMES,NATNAME, 1RS,RB,IND,rname,NATS,LINV,NO,NAT,alt,al,lmut,wmin,wmax, 2weights,nel,NATS_max,iweights,CUTOFF,SOSNM) implicit none INTEGER*4 np,iw,ix,jx,kx,IO,IOP,IOX,MX,NATNAME(*), 1INDBIG(IOX,MX),IND(*),NO,iao,NAT,IiO,NATS_max,IiOP, 1ii,jj,i,j,NATS,MXB,iaop,nel,iweights,IOPP parameter (MXB=140) real*8 bohr,bigpolr(3,3),bigpoli(3,3),biggpolr(3,3), 1bigpolvr(3,3),bigpolvi(3,3),CUTOFF,cut2, 1biggpoli(3,3),bigapolr(3,3,3),bigapoli(3,3,3),RS(3,MX), 1RB(3,MX),rname(IOX,3,MX),temr(3,3),temi(3,3),alt(*), 1polr(3,3),poli(3,3),gpolr(3,3),gpoli(3,3),x(3),xp(3), 1apolr(3,3,3),apoli(3,3,3),al(*),rt(3),rt2,tt(3,3), 1ttt(3,3,3),tar(3,3),tai(3,3),tcr1(3,3),tcr2(3,3), 1tci1(3,3),tci2(3,3),tcr(3,3),tci(3,3),pi,sf,taii(3,3), 1tair(3,3),tajr(3,3),taji(3,3),atar1(3,3),atar2(3,3), 1atai1(3,3),atai2(3,3),ataii(3,3),ataji(3,3),atajr(3,3), 1tvi(3),tvr(3),sui,sur,atair(3,3),wmin,wmax,weights(*), 1w1,w2,w,polvr(3,3),polvi(3,3),tx,ty,tz,wo,wau,wnm,dw,CM, 1wcm c alt resultant total big polarizability c al input overlap polarizabilities real*8, allocatable::opolr(:,:,:), 1opoli(:,:,:),gopolr(:,:,:),gopoli(:,:,:), 1opolvi(:,:,:),opolvr(:,:,:), 2aopolr(:,:,:,:),aopoli(:,:,:,:),am(:,:,:) logical LNAMES,LINV,lmut,SOSNM real*8 A,XS,XB,TOL,RMS integer*4 ITU,IIT COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),TOL,ITU,IIT,RMS c allocate(opolr(NO,3,3),opoli(NO,3,3), 1gopolr(NO,3,3),gopoli(NO,3,3), 1opolvr(NO,3,3),opolvi(NO,3,3), 1aopolr(NO,3,3,3),aopoli(NO,3,3,3),am(NO,3,3)) c bohr=0.529177D0 cut2=CUTOFF**2 write(6,8000)wmin,wmax,np,NO,LMUT,NAT 8000 format(/,/,' Dynamic polarizability transfer',/,/, 1 ' wmin:',f10.1,' wmax:',f10.1,' np :',i10,/, 2 ' NO :',i10, ' LMUT:',L10, ' NAT:',i10,/) c c loop over overlaps: c oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo DO 615 IO=1,NO write(6,6010)IO 6010 format(' overlap',i6,', transformation matrix:') c make the transformation matrix for whole overlap: call CRMSo(NAT,INDBIG,IOX,MX,IO,LNAMES,natname,RS,RB,IND, 2rname,NATS,LINV) do 610 ix=1,3 do 610 jx=1,3 610 am(IO,ix,jx)=A(ix,jx) 615 write(6,6019)A 6019 format(3(3f10.3,/)) c oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo c c loop over frequencies CM=219474.63d0 dw=(wmax-wmin)/dble(np-1) wo=wmin-dw do 1 iw=1,np wo=wo+dw if(SOSNM)then wau=(1.0d7/wo)/CM else wau=wo/CM endif wcm=wau*CM wnm=1.0d7/wcm c total polarizability call zp(biggpolr,biggpoli,bigpolr,bigpoli,bigpolvr,bigpolvi, 1bigapolr,bigapoli) c IO, loop over overlaps 551 55155155155155155155155155155155155 DO 551 IO=1,NO do 613 ix=1,3 do 613 jx=1,3 613 A(ix,jx)=am(IO,ix,jx) c c if named fragments, fill alpha, G, A for each overlap: if(lnames)then IiO=IO else IiO=1 endif c transcript polarizabilities: call trpol(polr,poli,polvr,polvi,gpolr,gpoli, 1apolr,apoli,al,nel,np,IiO,iw) c c rotate polarizabilities and write to new strings: call rpp( NO,IO, opolr , polr ,A) call rpp( NO,IO, opolvr, polvr,A) call rpp( NO,IO,gopolr ,gpolr ,A) call rppp(NO,IO,aopolr ,apolr ,A) call rpp( NO,IO, opoli , poli ,A) call rpp( NO,IO, opolvi, polvi,A) call rpp( NO,IO,gopoli ,gpoli ,A) call rppp(NO,IO,aopoli ,apoli ,A) c c construct alpha,G,A of all system: c c loop over atoms in this overlap ppppppppppppppppppppppppppppppp do 4 iao=1,NAT if(INDBIG(IO,iao).GT.0)then if(iweights.lt.0)then if(lnames)then w=weights((iw-1)*NO*NATS_max+NATS_max*(IiO-1)+INDBIG(IO,IAO)) else w=weights((iw-1)*NATS_max+INDBIG(IO,IAO)) endif else w=weights(NATS_max*(IiO-1)+INDBIG(IO,IAO)) endif c plain sum: c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp do 5519 IX=1,3 do 5519 JX=1,3 bigpolr(IX,JX) =bigpolr( IX,JX) +opolr(IO,IX,JX)*w bigpoli(IX,JX) =bigpoli( IX,JX) +opoli(IO,IX,JX)*w bigpolvr(IX,JX) =bigpolvr( IX,JX) +opolvr(IO,IX,JX)*w bigpolvi(IX,JX) =bigpolvi( IX,JX) +opolvi(IO,IX,JX)*w biggpolr(IX,JX)=biggpolr(IX,JX)+gopolr(IO,IX,JX)*w biggpoli(IX,JX)=biggpoli(IX,JX)+gopoli(IO,IX,JX)*w do 5519 KX=1,3 bigapoli(IX,JX,KX)=bigapoli(IX,JX,KX)+aopoli(IO,IX,JX,KX)*w 5519 bigapolr(IX,JX,KX)=bigapolr(IX,JX,KX)+aopolr(IO,IX,JX,KX)*w c pppppppppppppppppppppppppppppppppppppppppppppppppppppppppppp x(1)=RB(1,iao)*w/bohr x(2)=RB(2,iao)*w/bohr x(3)=RB(3,iao)*w/bohr c polarizability part for origin dependence of G call odg(biggpolr,x(1),x(2),x(3),opolr,IO,NO) call odg(biggpoli,x(1),x(2),x(3),opoli,IO,NO) c c polarizability part for origin dependence of A: call oda(bigapolr,x(1),x(2),x(3),opolr,IO,NO) call oda(bigapoli,x(1),x(2),x(3),opoli,IO,NO) endif 4 continue c loop over atoms in this overlap ppppppppppppppppppppppppppppppp 551 continue c IO, loop over overlaps 551 55155155155155155155155155155155155 c correction of polarizability for mutual interaction of groups if(lmut)then do 608 IO=1,NO if(lnames)then IiO=IO else IiO=1 endif c loop over atoms in this overlap 55555555555555555555 do 5 iao=1,NAT if(INDBIG(IO,iao).GT.0)then if(iweights.lt.0)then if(lnames)then w1=weights((iw-1)*NO*NATS_max+NATS_max*(IiO-1)+INDBIG(IO,IAO)) else w1=weights((iw-1)*NATS_max+INDBIG(IO,IAO)) endif else w1=weights(NATS_max*(IiO-1)+INDBIG(IO,IAO)) endif if(w1.gt.0.0d0)then x(1)=RB(1,iao)/bohr x(2)=RB(2,iao)/bohr x(3)=RB(3,iao)/bohr do 6081 IOP=1,NO if(lnames)then IiOP=IOP else IiOP=1 endif if(IOP.ne.IO)then c loop over atoms in this overlap 66666666666666666666 do 6 iaop=1,NAT if(INDBIG(IOP,iaop).GT.0)then if(iweights.lt.0)then if(lnames)then w2=weights((iw-1)*NO*NATS_max+ 1 NATS_max*(IiOP-1)+INDBIG(IOP,IAOP)) else w2=weights((iw-1)*NATS_max+INDBIG(IOP,IAOP)) endif else w2=weights(NATS_max*(IiOP-1)+INDBIG(IOP,IAOP)) endif w=w1*w2 xp(1)=RB(1,iaop)/bohr xp(2)=RB(2,iaop)/bohr xp(3)=RB(3,iaop)/bohr do 609 ix=1,3 609 rt(ix)=x(ix)-xp(ix) rt2=rt(1)**2+rt(2)**2+rt(3)**2 c is this atom pair part of any overlap? do 6082 IOPP=1,NO if(INDBIG(IOPP,iao).gt.0.and.INDBIG(IOPP,iaop).gt.0)then c disable using rt2 variable: rt2=0.0d0 goto 6083 endif 6082 continue 6083 if(rt2.gt.0.0d0)then if((CUTOFF.le.0.0d0.or.rt2.lt.cut2).and.w.ne.0.0d0)then c c T-tensor call mktt(rt,tt) c c TTT-tensor call mkttt(rt,ttt) c c aTa to a: call ctd(IO,opolr,temr,NO) call ctd(IO,opoli,temi,NO) call ct(tar,tt,temr) call ct(tai,tt,temi) call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(tcr1,temr,tar) call ct(tcr2,temi,tai) call ct(tci1,temi,tar) call ct(tci2,temr,tai) do 612 ix=1,3 do 612 jx=1,3 bigpolvi(ix,jx)=bigpolvi(ix,jx)+(tci1(ix,jx)+tci2(ix,jx))*w bigpolvr(ix,jx)=bigpolvr(ix,jx)+(tcr1(ix,jx)-tcr2(ix,jx))*w bigpoli(ix,jx)=bigpoli(ix,jx)+(tci1(ix,jx)+tci2(ix,jx))*w 612 bigpolr(ix,jx)=bigpolr(ix,jx)+(tcr1(ix,jx)-tcr2(ix,jx))*w c ATa-aTA to a: do 614 i=1,3 do 614 j=1,3 tcr(i,j)=0.0d0 tci(i,j)=0.0d0 do 614 ix=1,3 do 614 jx=1,3 do 614 kx=1,3 tcr(i,j)=tcr(i,j)+ 1 (aopolr(IO ,ix,jx,i)*ttt(ix,jx,kx)*opolr(IOP,kx,j)- 1 aopolr(IOP,ix,jx,j)*ttt(ix,jx,kx)*opolr(IO ,kx,i))/3.0d0+ 1 (-aopoli(IO ,ix,jx,i)*ttt(ix,jx,kx)*opoli(IOP,kx,j)+ 1 aopoli(IOP,ix,jx,j)*ttt(ix,jx,kx)*opoli(IO ,kx,i))/3.0d0 614 tci(i,j)=tci(i,j)+ 1 (aopolr(IO ,ix,jx,i)*ttt(ix,jx,kx)*opoli(IOP,kx,j)- 1 aopolr(IOP,ix,jx,j)*ttt(ix,jx,kx)*opoli(IO ,kx,i))/3.0d0+ 1 (aopoli(IO ,ix,jx,i)*ttt(ix,jx,kx)*opolr(IOP,kx,j)- 1 aopoli(IOP,ix,jx,j)*ttt(ix,jx,kx)*opolr(IO ,kx,i))/3.0d0 call scaleww(tcr,tci,w) call tadd(bigpolr,tcr) call tadd(bigpoli,tci) call tadd(bigpolvr,tcr) call tadd(bigpolvi,tci) c c GTG to alpha: pi=4.0d0*atan(1.0d0) sf=(2.0d0*pi/(wnm*10.0d0/bohr))**2 do 616 i=1,3 do 616 j=1,3 tcr(i,j)=0.0d0 tci(i,j)=0.0d0 do 616 ix=1,3 do 616 jx=1,3 tci(i,j)=tci(i,j) 1 +gopoli(IO ,i,jx)*tt(ix,jx)*gopolr(IOP,j,jx)*sf 1 +gopolr(IO ,i,jx)*tt(ix,jx)*gopoli(IOP,j,jx)*sf 616 tcr(i,j)=tcr(i,j) 1 +gopolr(IO ,i,jx)*tt(ix,jx)*gopolr(IOP,j,jx)*sf 1 -gopoli(IO ,i,jx)*tt(ix,jx)*gopoli(IOP,j,jx)*sf call scaleww(tcr,tci,w) call tadd(bigpolr,tcr) call tadd(bigpoli,tci) call tadd(bigpolvr,tcr) call tadd(bigpolvi,tci) c c contribution of alpha...alpha coupling to G: call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(tajr,tt,temr) call ct(taji,tt,temi) call ctd(IO,opolr,temr,NO) call ctd(IO,opoli,temi,NO) call ct(atar1,temr,tajr) call ct(atar2,temi,taji) call ct(atai1,temr,taji) call ct(atai2,temi,tajr) do 620 ix=1,3 call xvp(tx,ty,tz,x(1),x(2),x(3), 1 atar1(1,ix)-atar2(1,ix),atar1(2,ix)-atar2(2,ix), 1 atar1(3,ix)-atar2(3,ix)) tcr(ix,1)=-0.5d0*tx tcr(ix,2)=-0.5d0*ty tcr(ix,3)=-0.5d0*tz call xvp(tx,ty,tz,x(1),x(2),x(3), 1 atai1(1,ix)+atai2(1,ix),atai1(2,ix)+atai2(2,ix), 1 atai1(3,ix)+atai2(3,ix)) tci(ix,1)=-0.5d0*tx tci(ix,2)=-0.5d0*ty 620 tci(ix,3)=-0.5d0*tz call scaleww(tcr,tci,w) call tadd(biggpolr,tcr) call tadd(biggpoli,tci) c c contribution of alpha...A coupling to G do 624 ix=1,3 do 624 jx=1,3 taii(ix,jx)=0.0d0 taji(ix,jx)=0.0d0 tair(ix,jx)=0.0d0 tajr(ix,jx)=0.0d0 do 624 i=1,3 do 624 j=1,3 taii(ix,jx)=taii(ix,jx)+ttt(ix,i,j)*aopoli(IO ,i,j,jx) taji(ix,jx)=taji(ix,jx)+ttt(ix,i,j)*aopoli(IOP,i,j,jx) tair(ix,jx)=tair(ix,jx)+ttt(ix,i,j)*aopolr(IO ,i,j,jx) 624 tajr(ix,jx)=tajr(ix,jx)+ttt(ix,i,j)*aopolr(IOP,i,j,jx) do 625 ix=1,3 do 625 jx=1,3 ataii(ix,jx)=0.0d0 ataji(ix,jx)=0.0d0 atair(ix,jx)=0.0d0 atajr(ix,jx)=0.0d0 do 625 i=1,3 ataii(ix,jx)=ataii(ix,jx)+opolr(IOP,ix,i)*taii(i,jx) 1 +opoli(IOP,ix,i)*tair(i,jx) ataji(ix,jx)=ataji(ix,jx)+opolr(IO ,ix,i)*taji(i,jx) 1 +opoli(IO ,ix,i)*tajr(i,jx) atair(ix,jx)=atair(ix,jx)+opolr(IOP,ix,i)*tair(i,jx) 1 -opoli(IOP,ix,i)*taii(i,jx) 625 atajr(ix,jx)=atajr(ix,jx)+opolr(IO ,ix,i)*tajr(i,jx) 1 -opoli(IO ,ix,i)*taji(i,jx) do 622 ix=1,3 call xvp(tcr(ix,1),tcr(ix,2),tcr(ix,3),-x(1)/6.0d0, 1 -x(2)/6.0d0,-x(3)/6.0d0,atair(ix,1)-atajr(1,ix), 1 atair(ix,2)-atajr(2,ix),atair(ix,3)-atajr(3,ix)) 622 call xvp(tci(ix,1),tci(ix,2),tci(ix,3),-x(1)/6.0d0, 1 -x(2)/6.0d0,-x(3)/6.0d0,ataii(ix,1)-ataji(1,ix), 1 ataii(ix,2)-ataji(2,ix),ataii(ix,3)-ataji(3,ix)) call scaleww(tcr,tci,w) call tadd(biggpoli,tci) call tadd(biggpolr,tcr) c c contribution of Alpha..G to G: call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(taji,tt,temi) call ct(tajr,tt,temr) do 623 ix=1,3 do 623 jx=1,3 tci(jx,ix)=-( 1 gopolr(IO,jx,1)*taji(1,ix)+gopoli(IO,jx,1)*tajr(1,ix)+ 1 gopolr(IO,jx,2)*taji(2,ix)+gopolr(IO,jx,3)*taji(3,ix)+ 1 gopoli(IO,jx,2)*tajr(2,ix)+gopoli(IO,jx,3)*tajr(3,ix)) 623 tcr(jx,ix)=-( 1 gopolr(IO,jx,1)*tajr(1,ix)-gopoli(IO,jx,1)*taji(1,ix)+ 1 gopolr(IO,jx,2)*tajr(2,ix)+gopolr(IO,jx,3)*tajr(3,ix)- 1 gopoli(IO,jx,2)*taji(2,ix)-gopoli(IO,jx,3)*taji(3,ix)) call scaleww(tcr,tci,w) call tadd(biggpolr,tcr) call tadd(biggpoli,tci) c contribution of alpha..alpha to A call ctd(IOP,opolr,temr,NO) call ctd(IOP,opoli,temi,NO) call ct(taji,tt,temi) call ct(tajr,tt,temr) call ctd(IO,opolr,temr,NO) call ctd(IO,opoli,temi,NO) tvi(1)=x(1)*temi(1,1)+x(2)*temi(2,1)+x(3)*temi(3,1) tvi(2)=x(1)*temi(1,2)+x(2)*temi(2,2)+x(3)*temi(3,2) tvi(3)=x(1)*temi(1,3)+x(2)*temi(2,3)+x(3)*temi(3,3) tvr(1)=x(1)*temr(1,1)+x(2)*temr(2,1)+x(3)*temr(3,1) tvr(2)=x(1)*temr(1,2)+x(2)*temr(2,2)+x(3)*temr(3,2) tvr(3)=x(1)*temr(1,3)+x(2)*temr(2,3)+x(3)*temr(3,3) do 626 ix=1,3 do 626 jx=1,3 do 626 kx=1,3 sui=1.5d0*( 1 (x(ix)*temi(jx,1)+x(jx)*temi(ix,1))*tajr(1,kx)+ 1 (x(ix)*temi(jx,2)+x(jx)*temi(ix,2))*tajr(2,kx)+ 1 (x(ix)*temi(jx,3)+x(jx)*temi(ix,3))*tajr(3,kx)+ 1 (x(ix)*temr(jx,1)+x(jx)*temr(ix,1))*taji(1,kx)+ 1 (x(ix)*temr(jx,2)+x(jx)*temr(ix,2))*taji(2,kx)+ 1 (x(ix)*temr(jx,3)+x(jx)*temr(ix,3))*taji(3,kx)) if(ix.eq.jx)sui=sui 1 -tvr(1)*taji(1,kx)-tvr(2)*taji(2,kx)-tvr(3)*taji(3,kx) 1 -tvi(1)*tajr(1,kx)-tvi(2)*tajr(2,kx)-tvi(3)*tajr(3,kx) sur=1.5d0*( 1 (x(ix)*temr(jx,1)+x(jx)*temr(ix,1))*tajr(1,kx)+ 1 (x(ix)*temr(jx,2)+x(jx)*temr(ix,2))*tajr(2,kx)+ 1 (x(ix)*temr(jx,3)+x(jx)*temr(ix,3))*tajr(3,kx)- 1 (x(ix)*temi(jx,1)+x(jx)*temi(ix,1))*taji(1,kx)- 1 (x(ix)*temi(jx,2)+x(jx)*temi(ix,2))*taji(2,kx)- 1 (x(ix)*temi(jx,3)+x(jx)*temi(ix,3))*taji(3,kx)) if(ix.eq.jx)sur=sur 1 -tvr(1)*tajr(1,kx)-tvr(2)*tajr(2,kx)-tvr(3)*tajr(3,kx) 1 +tvi(1)*taji(1,kx)+tvi(2)*taji(2,kx)+tvi(3)*taji(3,kx) bigapoli(ix,jx,kx)=bigapoli(ix,jx,kx)+sui*w 626 bigapolr(ix,jx,kx)=bigapolr(ix,jx,kx)+sur*w c contribution of alpha...A to A do 627 ix=1,3 do 627 jx=1,3 do 627 kx=1,3 sur=aopolr(IO,ix,jx,1)*tajr(1,kx) 1 +aopolr(IO,ix,jx,2)*tajr(2,kx)+aopolr(IO,ix,jx,3)*tajr(3,kx) 1 -aopoli(IO,ix,jx,1)*taji(1,kx) 1 -aopoli(IO,ix,jx,2)*taji(2,kx)-aopoli(IO,ix,jx,3)*taji(3,kx) sui=aopoli(IO,ix,jx,1)*tajr(1,kx) 1 +aopoli(IO,ix,jx,2)*tajr(2,kx)+aopoli(IO,ix,jx,3)*tajr(3,kx) 1 +aopolr(IO,ix,jx,1)*taji(1,kx) 1 +aopolr(IO,ix,jx,2)*taji(2,kx)+aopolr(IO,ix,jx,3)*taji(3,kx) bigapoli(ix,jx,kx)=bigapoli(ix,jx,kx)+sui*w 627 bigapolr(ix,jx,kx)=bigapolr(ix,jx,kx)+sur*w c endif c if(CUTOFF.eq.0.0d0.or.rt2.lt.cut2)then endif c if(rt2.gt.0.0d0)then endif 6 continue c iaop endif c IO<>IOP 6081 continue c IO endif c w1>0 endif 5 continue c iao 608 continue c IOP endif c record to total polarizability: ii=0 jj=0 do 3 ix=1,3 do 3 jx=1,3 ii=ii+1 alt(nel*(iw-1)+ii )=bigpolr( ix,jx) alt(nel*(iw-1)+ii+ 9)=bigpoli( ix,jx) alt(nel*(iw-1)+ii+90)=bigpolvr(ix,jx) alt(nel*(iw-1)+ii+99)=bigpolvi(ix,jx) alt(nel*(iw-1)+ii+18)=biggpolr(ix,jx) alt(nel*(iw-1)+ii+27)=biggpoli(ix,jx) do 3 kx=1,3 jj=jj+1 alt(nel*(iw-1)+jj+36)=bigapolr(jx,kx,ix) 3 alt(nel*(iw-1)+jj+63)=bigapoli(jx,kx,ix) 1 continue c iw-frequency RETURN END c ============================================================ subroutine scaleww(r,i,w) implicit none real*8 r(3,3),i(3,3),w integer*4 u,v do 1 u=1,3 do 1 v=1,3 r(u,v)=r(u,v)*w 1 i(u,v)=i(u,v)*w return end c ============================================================ subroutine zp(biggpolr,biggpoli,bigpolr,bigpoli,bigpolvr,bigpolvi, 1bigapolr,bigapoli) implicit none integer*4 ix,jx,kx real*8 biggpolr(3,3),biggpoli(3,3),bigpolr(3,3), 1bigpoli(3,3),bigpolvr(3,3),bigpolvi(3,3),bigapolr(3,3,3), 2bigapoli(3,3,3) do 607 ix=1,3 do 607 jx=1,3 biggpolr(ix,jx)=0.0d0 biggpoli(ix,jx)=0.0d0 bigpolr(ix,jx)=0.0d0 bigpoli(ix,jx)=0.0d0 bigpolvr(ix,jx)=0.0d0 bigpolvi(ix,jx)=0.0d0 do 607 kx=1,3 bigapolr(ix,jx,kx)=0.0d0 607 bigapoli(ix,jx,kx)=0.0d0 return end c ============================================================ subroutine trpol(polr,poli,polvr,polvi,gpolr,gpoli, 1apolr,apoli,al,nel,np,IiO,iw) implicit none real*8 al(*),polr(3,3),poli(3,3),polvr(3,3),polvi(3,3), 1gpolr(3,3),gpoli(3,3),apolr(3,3,3),apoli(3,3,3) integer*4 ii,jj,ix,jx,kx,nel,np,IiO,iw ii=0 jj=0 do 2 ix=1,3 do 2 jx=1,3 ii=ii+1 polr(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii ) poli(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii+ 9) polvr(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii+90) polvi(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii+99) gpolr(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii+18) gpoli(ix,jx) =al(nel*((IiO-1)*np+iw-1)+ii+27) do 2 kx=1,3 jj=jj+1 apolr(jx,kx,ix)=al(nel*((IiO-1)*np+iw-1)+jj+36) 2 apoli(jx,kx,ix)=al(nel*((IiO-1)*np+iw-1)+jj+63) return end c ============================================================ subroutine odg(G,x1,x2,x3,a,IO,NO) implicit none integer*4 IO,IX,NO real*8 x1,x2,x3,G(3,3),a(NO,3,3) do 1 IX=1,3 G(IX,1)=G(IX,1)-0.5d0*(x2*a(IO,IX,3)-x3*a(IO,IX,2)) G(IX,2)=G(IX,2)-0.5d0*(x3*a(IO,IX,1)-x1*a(IO,IX,3)) 1 G(IX,3)=G(IX,3)-0.5d0*(x1*a(IO,IX,2)-x2*a(IO,IX,1)) return end c ============================================================ subroutine oda(a,x,y,z,al,IO,NO) implicit none integer*4 IO,NO,KX real*8 a(3,3,3),al(NO,3,3),x,y,z,sum,a1,a2,a3 do 1 KX=1,3 a1=al(IO,KX,1) a2=al(IO,KX,2) a3=al(IO,KX,3) sum=x*a1+y*a2+z*a3 a(1,1,KX)=a(1,1,KX)+1.5d0*(x*a1+x*a1)-sum a(1,2,KX)=a(1,2,KX)+1.5d0*(y*a1+x*a2) a(1,3,KX)=a(1,3,KX)+1.5d0*(z*a1+x*a3) a(2,1,KX)=a(2,1,KX)+1.5d0*(x*a2+y*a1) a(2,2,KX)=a(2,2,KX)+1.5d0*(y*a2+y*a2)-sum a(2,3,KX)=a(2,3,KX)+1.5d0*(z*a2+y*a3) a(3,1,KX)=a(3,1,KX)+1.5d0*(x*a3+z*a1) a(3,2,KX)=a(3,2,KX)+1.5d0*(y*a3+z*a2) 1 a(3,3,KX)=a(3,3,KX)+1.5d0*(z*a3+z*a3)-sum return end c ============================================================ subroutine xvp(x1,x2,x3,y1,y2,y3,z1,z2,z3) implicit none real*8 x1,x2,x3,y1,y2,y3,z1,z2,z3 x1=y2*z3-y3*z2 x2=y3*z1-y1*z3 x3=y1*z2-y2*z1 return end c ============================================================ subroutine rgr(s,NAT,IOX,MX,ii,g) implicit none character*(*) s integer*4 NAT,IOX,MX,i,ix,ii real*8 g(IOX,3*MX) open(90,file=s) do 1 i=1,NAT 1 read(90,*)(g(ii,ix+3*(i-1)),ix=1,3) close(90) write(6,*)s//' read' return end c ============================================================ subroutine svgr(g,n) implicit none real*8 g(*) integer*4 n,i,ix open(70,file='FILE.GR') do 1 i=1,n 1 write(70,700)(g(ix+3*(i-1)),ix=1,3) 700 format(3f15.9) close(70) write(6,*) write(6,*)'Gradient written to FILE.GR' write(3,*)'Gradient written to FILE.GR' return end c ============================================================ subroutine trgr(IPO,IA,GG,INDBIG,IOX,MX,A,ISUM,LNAMES,WF,IP,LWR, 1gradname) implicit none integer*4 IPO,IA,IOX,MX,INDBIG(IOX,MX),IX,IAS,ISUM,II,IP real*8 GG(*),A(3,3),GOLD,GNEW,WF(*),gradname(IOX,3*MX) logical LNAMES,LWR IAS=INDBIG(IPO,IA) DO 62 IX=1,3 GOLD=GG(IX+3*(IA-1)) C Zero out at the start of the summation IF(ISUM.EQ.1)GG(IX+3*(IA-1))=0.0d0 GNEW=0.0d0 DO 63 II=1,3 if(LNAMES)then GNEW=GNEW+gradname(IPO,II+3*(IAS-1))*A(IX,II) else GNEW=GNEW+gradname( 1,II+3*(IAS-1))*A(IX,II) endif 63 continue GG(IX+3*(IA-1))=GNEW*WF(IP)+GG(IX+3*(IA-1)) 62 IF(LWR)WRITE(3,3016)IA,IX,GOLD,GNEW,WF(IP) 3016 FORMAT(' grad ',2I3,F11.6,'/',F11.6,' (',F7.4,')') return end c ============================================================ subroutine trgrp(IA,GG,A,WF,IP,LWR,ggp) implicit none integer*4 IA,IX,IP real*8 GG(*),A(3,3),GOLD,GNEW,WF(*),ggp(*) logical LWR DO 62 IX=1,3 GOLD=GG(IX+3*(IA-1)) GNEW=ggp(1)*A(IX,1)+ggp(2)*A(IX,2)+ggp(3)*A(IX,3) GG(IX+3*(IA-1))=GNEW*WF(IP)+GG(IX+3*(IA-1)) 62 IF(LWR)WRITE(3,3016)IA,IX,GOLD,GNEW,WF(IP) 3016 FORMAT(' grad ',2I3,F11.6,'/',F11.6,' (',F7.4,')') return end c ============================================================ subroutine nv7(x,y,z) real*8 x,y,z,f f=x**2+y**2+z**2 if(f.ne.0.0d0)f=1.d0/dsqrt(f) x=x*f y=y*f z=z*f return end c ============================================================ subroutine rdnmr(nmr,s) implicit none integer*4 nat,i real*8 nmr(*) character*(*) s open(9,file=s) read(9,*)nat read(9,*) read(9,*) do 1 i=1,nat 1 read(9,*)nmr(i),nmr(i),nmr(i) read(9,*) close(9) write(6,*)s//' read' return end c ============================================================ subroutine rdspin(s,N,ss) implicit none character*(*) ss integer*4 N,N1,N3,LN,I,J real*8 s(N,N) open(7,file=ss) read(7,*) N1=1 111 N3=N1+4 IF(N3.GT.N)N3=N read(7,*) DO 130 LN=N1,N 130 READ(7,*)s(LN,N1),(s(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 111 DO 31 I=1,N DO 31 J=I,N 31 s(I,J)=s(J,I) close(7) write(6,*)ss//' read' return end c ============================================================ subroutine wrnmr(nat,nmr,iz) implicit none integer*4 nat,i,iz(*) real*8 nmr(*) open(9,file='FILE.NMR') write(9,900)nat 900 format(i5,' NMR isotropic shieldings',/, 1 ' atom sigma') write(9,901) 901 format(60(1h-)) do 1 i=1,nat 1 write(9,902)i,iz(i),nmr(i) 902 format(i5,i3,f12.4) write(9,901) close(9) write(6,*)'FILE.NMR written' return end c ============================================================ subroutine wspin(s,N) implicit none integer*4 N,N1,N3,LN,J real*8 s(N,N) open(88,file='FILE.SPI') write(88,880) 880 format(' NMR spin-spin coupling constants, Hz') N1=1 111 N3=N1+4 IF(N3.GT.N)N3=N write(88,882)(J,J=N1,N3) 882 format(3x,5i14) DO 130 LN=N1,N 130 write(88,881)LN,(s(LN,J),J=N1,MIN(LN,N3)) 881 format(i7,5d14.6) N1=N1+5 IF(N3.LT.N)GOTO 111 close(88) write(6,*)'FILE.SPI written' return end c ============================================================d subroutine asnq(QB,qq,MX,NAT) implicit none integer*4 MX,NAT,ia,ix real*8 QB(MX,3,3,3),qq(3*NAT,6) DO 1 ia=1,NAT DO 1 ix=1,3 QB(ia,ix,1,1)=qq(ix+3*(ia-1),1) QB(ia,ix,2,2)=qq(ix+3*(ia-1),2) QB(ia,ix,3,3)=qq(ix+3*(ia-1),3) QB(ia,ix,1,2)=qq(ix+3*(ia-1),4) QB(ia,ix,2,1)=qq(ix+3*(ia-1),4) QB(ia,ix,1,3)=qq(ix+3*(ia-1),5) QB(ia,ix,3,1)=qq(ix+3*(ia-1),5) QB(ia,ix,2,3)=qq(ix+3*(ia-1),6) 1 QB(ia,ix,3,2)=qq(ix+3*(ia-1),6) return end c ============================================================d subroutine asnqn(qsname,IOX,MX,i,qq,NATS) implicit none integer*4 IOX,MX,i,ia,ix,NATS real*8 qsname(IOX,MX,3,3,3),qq(3*NATS,6) do 32 ia=1,NATS do 32 ix=1,3 qsname(i,ia,ix,1,1)=qq(ix+3*(ia-1),1) qsname(i,ia,ix,2,2)=qq(ix+3*(ia-1),2) qsname(i,ia,ix,3,3)=qq(ix+3*(ia-1),3) qsname(i,ia,ix,1,2)=qq(ix+3*(ia-1),4) qsname(i,ia,ix,2,1)=qq(ix+3*(ia-1),4) qsname(i,ia,ix,1,3)=qq(ix+3*(ia-1),5) qsname(i,ia,ix,3,1)=qq(ix+3*(ia-1),5) qsname(i,ia,ix,2,3)=qq(ix+3*(ia-1),6) 32 qsname(i,ia,ix,3,2)=qq(ix+3*(ia-1),6) return end c ============================================================d subroutine asqq(q,H,nat,ic,ix,ia) implicit none integer*4 ic,nat,a(6),b(6),ix,ia,i real*8 H(3,3),q(3*nat,6) data a/1,2,3,1,1,2/ data b/1,2,3,2,3,3/ if(ic.eq.1)then do 1 i=1,6 H(a(i),b(i))=q(ix+3*(ia-1),i) 1 H(b(i),a(i))=q(ix+3*(ia-1),i) else do 2 i=1,6 2 q(ix+3*(ia-1),i)=H(a(i),b(i)) endif return end c ============================================================d subroutine reqd(fn,q,nat,P,MX,X) implicit none character*(*) fn integer*4 i,ia,ix,nat,MX real*8 q(3*nat,6),P(MX,3,3),X(3,NAT),q0(6) real*8,allocatable::d(:) character*1 xyz(3) data xyz/'x','y','z'/ open(9,file=fn,status='old') do 1 i=1,9 1 read(9,*) do 2 ia=1,nat do 2 ix=1,3 2 read(9,93)(q(ix+3*(ia-1),i),i=1,6) 93 format(9x,3e13.4) close(9) write(6,*)fn//' read' c put to atomic origins: call pqao(X,NAT,1.0d0,P,MX,q) c write the local tensors: do 6 i=1,6 6 q0(i)=0.0d0 allocate(d(18*nat)) do 7 ia=1,nat do 7 ix=1,3 do 7 i=1,6 7 d(ix+3*(ia-1)+3*nat*(i-1))=q(ix+3*(ia-1),i) call wrqd(fn//'.loc',q0,d,nat,xyz) return end c ============================================================d subroutine pqao(X,NAT,si,P,MX,q) c si = 1 .. put quadrupole derivatives to atomic origins c si = -1 .. put quadrupole derivatives to common origin c q = ri ri c q' = q + u R + R u implicit none integer*4 NAT,ia,ix,a,b,MX real*8 X(3,NAT),BOHR,si,P(MX,3,3),H(3,3),HH(3,3), 1q(3*nat,6) real*8,allocatable::R(:,:) BOHR=0.52917705993d0 allocate(R(3,NAT)) DO 3 ia=1,NAT DO 3 ix=1,3 3 R(ix,ia)=X(ix,ia)/BOHR do 4 ia=1,nat do 4 ix=1,3 call asqq(q,H,nat,1,ix,ia) do 5 a=1,3 do 5 b=1,3 5 HH(a,b)=H(a,b)-si*(R(a,ia)*P(ia,ix,b)+P(ia,ix,a)*R(b,ia)) 4 call asqq(q,HH,nat,0,ix,ia) return end c ============================================================d subroutine wrqd(fn,q,d,nat,xyz) implicit none character*(*) fn character*1 xyz(*) integer*4 i,ia,ix,nat real*8 q(*),d(*) open(9,file=fn) write(9,90)nat,'XX',q(1),'YY',q(2),'ZZ',q(3), 1 'XY',q(4),'XZ',q(5),'YZ',q(6) 90 format(i6,' atoms',/, 1 ' Equilibrium quadrupole, sum(i)qi ri ri, atomic units:',/, 12(3(3x,a2,2x,f19.4),/)) write(9,91) 91 format(' Quadrupole derivatives, atomic units: ',/, 1 ' Atom Coord XX YY ZZ ',/, 1 ' XY XZ YZ ',/, 1 ' ------------------------------------------------') do 2 ia=1,nat do 2 ix=1,3 if(ix.eq.1)then write(9,92)ia,xyz(ix),(d(ix+3*(ia-1)+3*nat*(i-1)),i=1,6) 92 format(i5,3x,a1,3e13.4,/,9x,3e13.4) else write(9,93)xyz(ix),(d(ix+3*(ia-1)+3*nat*(i-1)),i=1,6) 93 format(8x,a1,3e13.4,/,9x,3e13.4) endif 2 continue close(9) write(6,*)fn//' written' return end c ============================================================d subroutine wrqq(MX,PB,QB,R,NAT) c writes quadrupole derivatives ion disk implicit none integer*4 MX,NAT,ia,ix,i real*8 QB(MX,3,3,3),R(3,NAT),q0(6),PB(MX,3,3) real*8 ,allocatable::q(:,:),d(:) character*1 xyz(3) data xyz/'x','y','z'/ allocate(q(3*nat,6)) do 2 ia=1,NAT do 2 ix=1,3 q(ix+3*(ia-1),1)=QB(ia,ix,1,1) q(ix+3*(ia-1),2)=QB(ia,ix,2,2) q(ix+3*(ia-1),3)=QB(ia,ix,3,3) q(ix+3*(ia-1),4)=QB(ia,ix,1,2) q(ix+3*(ia-1),5)=QB(ia,ix,1,3) 2 q(ix+3*(ia-1),6)=QB(ia,ix,2,3) call pqao(R,NAT,-1.0d0,PB,MX,q) do 1 i=1,6 1 q0(i)=0.0d0 allocate(d(18*nat)) do 7 ia=1,nat do 7 ix=1,3 do 7 i=1,6 7 d(ix+3*(ia-1)+3*nat*(i-1))=q(ix+3*(ia-1),i) call wrqd('FILE.QEN',q0,d,nat,xyz) return end