PROGRAM CCTN IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind PARAMETER (MX=2691,MXB=40,IOX=199,n30=170) 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 common/mc/FF(3*MX,3*MX),R(3,MX),IND(MX),IB(MX,MXB),ND(MX), 1NDS(MX),RS(3,MX),FT(3*MX,3*MX),IS(MX,MXB),AM(MX,3,3), 2INDBIG(IOX,MX),IBONDSB(4,MX), 3NBO(MX),NBOB(MX),PS(MX,3,3),AS(MX,3,3),PB(MX,3,3),AB(MX,3,3), 4ALPHA(3*MX,3,3),ALPHAS(3*MX,3,3),AA(3*MX,3,3,3),alder(3,3), 5AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3), 5AAT(3*MX,3,3,3),GT(3*MX,3,3) 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 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 cb(:),cs(:),ct(:),FFS(:) 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),natname(IOX),NAME, 3gpo(3,3),apo(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3) c call inipol(po,gpo,apo) bohr=0.529177D0 iffs=-1 lpol=.false. lpolg=.false. lpola=.false. OPEN(3,FILE='CCT.OUT') WRITE(3,3000) 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',/, 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',/) 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) CLOSE(2) c c if tensors constructed from named fragments, load them all: if(LNAMES)then NATS_max=0 DO 11 I=1,NO NE=80 do 10 NE=80,1,-1 10 if(NAME(I)(NE:NE).ne.' ')goto 12 12 basicname(1:NE)=NAME(I)(1: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 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 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(MX,NAT,FF) CLOSE(2) WRITE(3,3001)NAT WRITE(6,3001)NAT 3001 FORMAT(' Force field of the big molecule found,',/,I10,' atoms') ENDIF c c Special section for named fragments: c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn if(lnames)then allocate(FFS(1:NO*(3*NATS_max)**2),stat=iffs) if(iffs.ne.0)call report('FFS cannot be allocated.') do 14 i=1,NO NATS=natname(i) NE=80 do 15 NE=80,1,-1 15 if(NAME(I)(NE:NE).ne.' ')goto 16 16 basicname=NAME(I)(1:NE) write(6,*) write(6,*)' Fragment ',basicname(1:NE) c fname=basicname(1:NE)//'.fc' write(6,*)fname(1:NE+3) INQUIRE(FILE=fname,EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,*)'Force field of small not found' WRITE(6,*)'Force field of small not found' IF(LDIA.OR.LOFF)CALL REPORT('Cannot do force field') ELSE OPEN(2,FILE=fname,STATUS='OLD') CALL READFF(MX,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) WRITE(3,3002)NATS WRITE(6,3002)NATS 3002 FORMAT(' Force field found,',/,I10,' atoms') ENDIF c IF(LABS.OR.LVCD)THEN nname=basicname(1:NE)//'.ten' write(6,*)nname(1:NE+4) OPEN(15,FILE=nname,STATUS='OLD') do 19 ia=1,NATS do 19 ix=1,3 19 RS(ix,ia)=rname(i,ix,ia) 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) 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 20 iax=1,3*NATS do 20 ix=1,3 do 20 jx=1,3 alphaname(i,iax,ix,jx)=ALPHAS(iax,ix,jx) gname(i,iax,ix,jx)=GS(iax,ix,jx) do 20 kx=1,3 20 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) 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) ENDIF C 14 continue c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn else 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),stat=iffs) if(iffs.ne.0)call report('FFS cannot be allocated.') OPEN(2,FILE='SMALL.FC',STATUS='OLD') CALL READFF(MX,NATS,FT) CLOSE(2) WRITE(3,3002)NATS WRITE(6,3002)NATS do 25 ii=1,3*NATS do 25 jj=1,3*NATS 25 FFS((ii-1)*3*NATS+jj)=FT(ii,jj) ENDIF IF(LABS.OR.LVCD)THEN OPEN(15,FILE='SMALL.TEN',STATUS='OLD') CALL READTEN(MX,PS,AS,LAPT,RS,NATS) CLOSE(15) ENDIF 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) 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 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 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 WRITE(6,*)' Transferring the tensors ...' CALL IMPROVE(FF,FFS,R,RS,NAT,NATS,NS,IND,IB,IS, 1AM,INDBIG,NO,IBONDSB,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA, 2IWG,LROA,LRAM,ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LINV,LSELECT, 3LNAMES,LHALTONERROR,CUTOFF,n1,n2,okind,NATS_max,lpol, 1lpolg,lpola,lq,ND,LDATOM) if(iffs.ne.-1)deallocate(FFS) c OPEN(2,FILE='FILE.FC') CALL WRITEFF(MX,NAT,FF) CLOSE(2) WRITE(3,*)' File FILE.FC written ' WRITE(6,*)' File FILE.FC written ' IF(LABS.OR.LVCD)CALL WRITETEN(MX,PB,AB,R,NAT) 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 24 i=1,3*NAT do 24 ix=1,3 do 24 jx=1,3 G(i,ix,jx)=0.0d0 do 24 kx=1,3 24 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,*)' File FILE.X written ' WRITE(6,*)' File FILE.X written ' c 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') ic=0 allocate(cs((3*NATS)**3),stat=ic) if(ic.ne.0)call report('Cijk(small) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('Cijk(names) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('Cijk(t) cannot be allocated.') 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),stat=ic) IF(OPT)then if(ic.ne.0)call report('Cijk(big) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('Diijk(small) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('Diijk(names) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('Diijk(t) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('Diijk(big) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('37iijk(small) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('37iijk(names) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('37iijk(t) cannot be allocated.') 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),stat=ic) if(ic.ne.0)call report('37iijk(b) cannot be allocated.') 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=2691,MXB=40,IOX=199) 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 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),natname(IOX),NAME, 3gpo(3,3),apo(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3) 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=2691,MXB=40,IOX=199) 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 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),natname(IOX),NAME, 3gpo(3,3),apo(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3) 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 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) implicit none INTEGER*4 MX,MXB,IOX PARAMETER (MX=2691,MXB=40,IOX=199) 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(*) real*8 FFB(3*MX,3*MX),RB(3,MX),RS(3,NATS),x,y,z, 1TS(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), 6RM(IOX,3,3),adta(3,3) real*8 FFS(*) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LNAMES,LHALTONERROR,lpol,li1,li2,lpolg,lpola,lq,LDATOM real*8 A,XS,XB,TOL,RMS,tt(3,3),rt(3),d1,d2,dd,ttt(3,3,3),eps integer*4 ITU,IIT,ido,iao,ic,iop,noo,nu,bodyn(IOX), 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 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, 1opol(IOX,3,3),coo(IOX,3),ta(3,3),tc(3,3),era(3,3),sf,rt2,rt5, 2bohr,sf1,gpol,apol,apolname,gpolname,aopol(IOX,3,3,3), 3gopol(IOX,3,3),alst(3,3,3),glst(3,3,3),tem(3,3),tg(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,at(3,3) 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),natname(IOX),NAME, 3gpol(3,3),apol(3,3,3),apolname(IOX,3,3,3),gpolname(IOX,3,3) 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(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(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=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=zsum+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/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(lpolg)call rpp(IOX,IO,opol,pol,A) c c rotate G if(lpolg)call rpp(IOX,IO,gopol,gpol,A) write(6,*)IO,gopol(IO,1,2),gpol(1,2) 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)then bigapol(IX,JX,KX)=bigapol(IX,JX,KX)-coo(IO,1)*opol(IO,KX,1) 1 -coo(IO,2)*opol(IO,KX,2)-coo(IO,3)*opol(IO,KX,3) endif 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 rt5=rt2**2*dsqrt(rt2) c c T-tensor call mktt(rt,tt) c c TTT-tensor 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 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 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(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))*ta(1,kx)+ 1 (coo(IO,ix)*tem(jx,2)+coo(IO,jx)*tem(ix,2))*ta(2,kx)+ 1 (coo(IO,ix)*tem(jx,3)+coo(IO,jx)*tem(ix,3))*ta(3,kx)) if(ix.eq.jx)su=su-tv(1)*ta(1,kx)-tv(2)*ta(2,kx)-tv(3)*ta(3,kx) 626 bigapol(ix,jx,kx)=bigapol(ix,jx,kx)+su endif c contribution of alpha...A to A if(lpola)then 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,*)'Neighbor',neighbor(IO,1),neighbor(IO,2) write(3,*)' body:' 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 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 DO 20 IA= istart,iend write(6,*)IA,'. atom of ',NAT write(3,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) DO 20 JA=IA,NAT write(3,*) write(3,3008)IA,JA 3008 format(' Pair ',2I6) IF(.NOT.LOFF.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 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 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 : ',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 field: IF((LDIA.AND.IA.EQ.JA).OR.(LOFF.AND.IA.NE.JA))THEN if(LNAMES)then c retrieve FF of fragment IPO c c load only if not in memory already if(IPOOLD.EQ.0.or.IPOOLD.NE.IPO)then NATS=NATNAME(IPO) IPOOLD=IPO endif ido=IPO ELSE ido=1 ENDIF 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 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,')') 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)then 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 c overlapping - only body+ further closer end c and do nto correct atoms on IPO termini c is atom IA part of ends2 if IPO?: c ix is the index over the neighbors do 564 ix=1,2 do 564 jend=1,endn(ix,2,IPO) 564 if(jend.eq.IA)goto 51 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 distance to differentiated atom: do 156 ix=1,3 56 rt(ix)=coo(iop,ix)-RB(ix,IA)/bohr else c distance to overlap center: do 156 ix=1,3 156 rt(ix)=coo(iop,ix)-coo(IPO,ix) endif rt2=rt(1)**2+rt(2)**2+rt(3)**2 if(rt2.gt.0.0d0)then rt5=rt2**2*dsqrt(rt2) c T-tensor call mktt(rt,tt) c call ctd(iop,opol,al,IOX) c T.alpha: call ct(ta,tt,al) c alpha.T: call ct(at,al,tt) do 603 IX=1,3 IIND=3*(IA-1)+IX call ctd(IIND,ALPHA,alder,3*MX) c alphader T.alpha: call ct(adta,alder,ta) c alphader T.alpha + alpha.T alphader : call std(adta) do 604 ig=1,3 do 604 id=1,3 604 tc(ig,id)=adta(ig,id) call multt(tc,sf*WF(IP)) c alpha=alpha+alder T al + al T alder: call addtoG(IIND,ALPHA,3*MX,tc) c ERa = EPS_bgd R_g Al_da call dera(era,rt,al) c ERad = EPS_bgd R_g Alder_da call dera(erad,rt,alder) c eps R al T: call ct(tc,era,tt) c eps R al T alder: call ct(tg,tc,alder) c eps R alder T: call ct(tc,erad,tt) c eps R alder T al: call ct(tem,tc,al) c eps R al T alder + eps R alder T al: call tadd(tg,tem) call multt(tg,-sf*WF(IP)*0.5d0) call addtoG(IIND,G,3*MX,tg) DO 604 IG=1,3 DO 604 ID=1,3 DO 604 IE=1,3 604 AA(IIND,IG,ID,IE)=AA(IIND,IG,ID,IE) +1.5d0*(rt(IG)*adta(ID,IE)+rt(ID)*adta(IG,IE)) DO 605 IG=1,3 DO 605 IE=1,3 605 AA(IIND,IG,IG,IE)=AA(IIND,IG,IG,IE) -rt(1)*adta(1,IE)-rt(2)*adta(2,IE)-rt(3)*adta(3,IE) 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 51 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',5H' #it,/, 2 1X,'-------------------------------------------------') DO 1 IA=1,NS DO 3 I=1,3 TB(I) =0.0d0 3 TS(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 TS(J)=TS(J)+XS(ISAT,J) DO 5 ISAT=1,ABS(IND(IA))+1 DO 5 J=1,3 XS(ISAT,J)=XS(ISAT,J)-TS(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) CHARACTER*78 TITDUM READ(2,2000) TITDUM 2000 FORMAT(A78) 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) CHARACTER*78 TITDUM READ(2,2000) TITDUM 2000 FORMAT(A78) 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(N0,NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*N0,3*N0) 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 WRITEFF(MX,NA,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(3*MX,3*MX) 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=40) 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)=0.2d0 ANG(2)=0.0d0 ANG(3)=0.0d0 DO 1 I=1,3 1 P(1,I)=ANG(I) Y(1)=FU(ANG) ANG(2)=ANG(1) ANG(1)=0.0d0 DO 2 I=1,3 2 P(2,I)=ANG(I) Y(2)=FU(ANG) ANG(3)=ANG(2) ANG(2)=0.0d0 DO 3 I=1,3 3 P(3,I)=ANG(I) Y(3)=FU(ANG) ANG(1)=ANG(3) ANG(2)=ANG(3) ANG(3)=0.0d0 DO 4 I=1,3 4 P(4,I)=ANG(I) Y(4)=FU(ANG) FTOL=1.0d-7 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=40) 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=2691) 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=2691) 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=2691) 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,*)LDUM,IXDUM,(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,*)LDUM,IXDUM,(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,*)LDUM,IXDUM,(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 SUBROUTINE TTT1(ALPHA,A,G,ICO,coo) 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 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(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' 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(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' 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=2691) 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) implicit none INTEGER*4 NS,IB,IS,IND,MXT,MXB,INDBIG,IOX,NO, 1NAT,n1,n2,okind,I,IA,IAB,IC,IEND,II,IO,ISMALL,IWG,J,JA,NDD REAL*8 CUTOFF,ALIM 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 c n1=0 n2=0 okind=0 ALIM=0.00d0 LDATOM=.false. LRDO=.FALSE. loldformat=.true. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LOFF=.TRUE. lq=.false. LDIA=.TRUE. LXX3=.false. LXX4=.false. L36=.false. L37=.false. LHALTONERROR=.TRUE. CUTOFF=-1.0d5 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 IF(STR(1:7).EQ.'OLDFORM')READ(2,*)loldformat IF(STR(1:4).EQ.'N1N2') READ(2,*)n1,n2 IF(STR(1:4).EQ.'ALIM') READ(2,*)ALIM IF(STR(1:3).EQ.'OKI') READ(2,*)okind IF(STR(1:3).EQ.'IWG') READ(2,*)IWG IF(STR(1:6).EQ.'LNAMES') READ(2,*)LNAMES IF(STR.EQ.'LSELECT') READ(2,*)LSELECT IF(STR(1:4).EQ.'LINV') READ(2,*)LINV IF(STR(1:4).EQ.'LAPT') READ(2,*)LAPT IF(STR.EQ.'LSTRICT') READ(2,*)LSTRICT IF(STR(1:3).EQ.'LWR') READ(2,*)LWR IF(STR(1:4).EQ.'LALF') READ(2,*)LALF IF(STR(1:4).EQ.'LRDO') READ(2,*)LRDO IF(STR(1:4).EQ.'LRAM') READ(2,*)LRAM IF(STR(1:4).EQ.'LROA') READ(2,*)LROA IF(STR(1:4).EQ.'LOFF') READ(2,*)LOFF IF(STR(1:4).EQ.'LDIA') READ(2,*)LDIA IF(STR(1:4).EQ.'LXX3') READ(2,*)LXX3 IF(STR(1:4).EQ.'LXX4') READ(2,*)LXX4 IF(STR(1:4).EQ.'L36 ') READ(2,*)L36 IF(STR(1:4).EQ.'LQ ') READ(2,*)LQ IF(STR(1:4).EQ.'L37 ') READ(2,*)L37 IF(STR(1:4).EQ.'LABS') READ(2,*)LABS IF(STR(1:4).EQ.'LVCD') READ(2,*)LVCD IF(STR(1:7).EQ.'LHALTON')READ(2,*)LHALTONERROR IF(STR(1:6).EQ.'CUTOFF') READ(2,*)CUTOFF IF(STR(1:6).EQ.'LDATOM') READ(2,*)LDATOM IF(STR.EQ.'NOPOLYM')THEN c outdated definition of the overlap: 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 CALL REPORT('Error in reading ') ENDIF IF(LNAMES)READ(2,2900)NAME(I) 2900 FORMAT(a80) 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 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=40) 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,*)'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=40) 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,*)'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) 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) 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) 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=40) 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) 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) 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) 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 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(IOX,IO,gopol,gpol,A) implicit none integer*4 IX,JX,JXP,IOX,IO,KXP,KX real*8 gopol(IOX,3,3,3),A(3,3),gpol(3,3,3),t(3,3,3) 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(IO,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(IO,gopol,t,IOX) do 3 IX=1,3 do 3 JX=1,3 do 3 KX=1,3 3 gopol(IO,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) real*8 p(3,3),g(3,3),a(3,3,3) 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) c B=B+C real*8 bigpol(3,3),tc(3,3) do 615 ix=1,3 do 615 jx=1,3 615 bigpol(ix,jx)=bigpol(ix,jx)+tc(ix,jx) return end c ============================================================ subroutine mktt(rt,tt) 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) real*8 era(3,3),r(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) real*8 tc(3,3),t 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) 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(9f7.3) return end c ============================================================ subroutine wrt3(IA,MX,ALPHA,s) 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(9f7.3) return end c ============================================================ call std(a) 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