PROGRAM CCTA IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=900,IOX=25,MXB=40) DIMENSION R(3,MX),ND(MX),NDS(MX),RS(3,MX), 2IBONDSB(4,MX),IBONDSS(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), 5AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3), 6ALPHAT(3*MX,3,3),ixb(MXB),ixs(MXB),XB(MXB,3),XS(MXB,3) LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LINV, 1LAPT,LALF,LRDO,LNAMES logical lxfile CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname common/cname/asname(IOX,MX,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, 3nboname(IOX,MX),ibondsname(IOX,4,MX),NDNAME(IOX,MX) common/var1/FF(3*MX*3*MX),FFS(3*MX*3*MX),pmax(MX*MX), 1imax(MX*MX),jmax(MX*MX),iomax(MX*MX),pmaxnew(MX*MX) c OPEN(3,FILE='CCTA.OUT') WRITE(3,3000) WRITE(6,3000) 3000 FORMAT(/, 1' Transfer of Molecular Tensors in Cartesian Coordinates',/,/, 8' FULLY AUTOMATIC VERSION',/,/, 8' By Petr Bour, Prague 1994-2002',/,/, 9' INPUT:',/, 9' Target: BIG.X coordinates',/, 1' BIG.FC force field (optional)',/, 3' AP.PRO weights for FF (optional)',/, 1' BIG.TEN APT and AAT tensors (optional)',/, 3' BIG.TTT ROA tensors (optional)',/, 3' Source: SMALL.X coordinates',/, 2' SMALL.FC force field',/, 2' SMALL.TEN APT and AAT',/, 4' SMALL.TTT ROA tensors',/, 4' (or named .x,.fc, .ten and .ttt)',/, 4' Options: CCTA.INP',/,/, 5' OUTPUT: FILE.X same as BIG.X',/, 7' FILE.FC improved BIG.FC',/, 3' AP.PRO actual weights for FF',/, 7' FILE.TEN improved BIG.TEN',/, 7' FILE.TTT improved BIG.TTT',/, 8' CCTA.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' WRITE(6,*)NAT, ' atoms in the big molecule' c c Open file with option and overlaps and read it in: OPEN(2,FILE='CCTA.INP',STATUS='OLD') CALL readopt(IOX,NO, 1LWR,LABS,LVCD,LOFF,LDIA,LROA,LINV,LAPT, 2LALF,LRDO,LNAMES,NAME,CUTOFF,n1,n2,ioo) CLOSE(2) c c if tensors constructed from named fragments, load them all: if(LNAMES)then 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 READRX(MX,NATS,RS,NDS,IBONDSS,NBO) NATNAME(I)=NATS do 13 ia=1,nat NDNAME(I,ia)=NDS(ia) NBONAME(I,ia)=NBO(ia) do 131 ib=1,NBO(ia) 131 ibondsname(I,ib,ia)=IBONDSS(ib,ia) do 13 ix=1,3 13 rname(I,ix,ia)=RS(ix,ia) WRITE(6,*)NAT, ' atoms in the fragment' 11 WRITE(3,*)NAT, ' atoms in the fragment' else inquire(file='SMALL.X',exist=lxfile) if(lxfile)then OPEN(2,FILE='SMALL.X',STATUS='OLD') CALL READRX(MX,NATS,RS,NDS,IBONDSS,NBO) CLOSE(2) else call report('SMALL.X not found') endif WRITE(3,*)NATS, ' atoms in the small molecule' WRITE(6,*)NATS, ' atoms in the small molecule' endif c names c INQUIRE(FILE='BIG.FC',EXIST=OPT) IF(.NOT.OPT)THEN call w36('File BIG.FC not found, constants set to zero') DO 1 I=1,3*NAT*3*NAT 1 FF(I)=0.0d0 ELSE OPEN(2,FILE='BIG.FC',STATUS='OLD') CALL READFF(3*MX*3*MX,3*NAT,FF) CLOSE(2) call w36('Force field of the big molecule found') ENDIF c INQUIRE(FILE='AP.PRO',EXIST=OPT) IF(.NOT.OPT)THEN WRITE(3,3003) WRITE(6,3003) 3003 FORMAT(' File AP.PRO not found, probabilities initialized') DO 5 J=1,NAT*NAT 5 pmax(J)=1.0d9 ELSE OPEN(2,FILE='AP.PRO',STATUS='OLD') CALL READFF(MX*MX,NAT,pmax) CLOSE(2) WRITE(3,3008) WRITE(6,3008) 3008 FORMAT(' Probabilities found') ENDIF DO 111 J=1,NAT*NAT pmaxnew(j)=1.0d9 imax(J)=0 iomax(J)=0 111 imax(J)=0 c c Special section for named fragments: if(lnames)then 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 call w36('Force field of small not found') IF(LDIA.OR.LOFF)CALL REPORT('Cannot do force field') ELSE c c just check here, if it is possible to load it: OPEN(2,FILE=fname,STATUS='OLD') CALL READFF(MX*3*3*MX,3*NATS,FFS) CLOSE(2) call w36('Force field of the small molecule found') 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)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) ICO=1 IF(LALF)THEN call w36('A and G obtained from alpha as DO') ICO=3 ENDIF CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,ICO,RS) 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 ENDIF C 14 continue else INQUIRE(FILE='SMALL.FC',EXIST=OPT) IF(.NOT.OPT)THEN call w36('Force field of small not found') IF(LDIA.OR.LOFF)CALL REPORT('Cannot do force field') ELSE OPEN(2,FILE='SMALL.FC',STATUS='OLD') CALL READFF(9*MX*MX,3*NATS,FFS) CLOSE(2) call w36('Force field of the small molecule found') 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)THEN OPEN(15,FILE='SMALL.TTT',STATUS='OLD') CALL READTTT(MX,ALPHAS,AAS,GS,NATS) CLOSE(15) ICO=1 IF(LALF)THEN call w36('A and G obtained from alpha as DO for small') ICO=3 ENDIF CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,ICO,RS) C Now small tensors are in atomic origins C IF(LRDO)THEN call w36('Alpha from ALPHA.TTT for small') OPEN(15,FILE='ALPHA.TTT',STATUS='OLD') CALL READTTT(MX,ALPHAT,AA,G,NAT) DO 4 I=1,3*NAT DO 4 IX=1,3 DO 4 JX=1,3 4 ALPHAS(I,IX,JX)=ALPHAT(I,IX,JX) CLOSE(15) ENDIF endif ENDIF c lnames c IF(LABS.OR.LVCD)THEN INQUIRE(FILE='BIG.TEN',EXIST=OPT) IF(.NOT.OPT)THEN call w36('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) call w36('Tensors BIG.TEN found') CLOSE(15) ENDIF ENDIF C INQUIRE(FILE='BIG.TTT',EXIST=OPT) IF(.NOT.OPT)THEN call w36('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 ICO=1 IF(LALF)THEN call w36('A and G obtained from alpha as DO for big') ICO=3 ENDIF CALL TTTT(3*MX,ALPHA,AA,G,NAT,ICO,R) C call w36('Looking for overlaps ...') CALL IMPROVE(R,RS,NAT,NATS,ND,NDS, 1NO,IBONDSB,NBOB,LWR,LABS,LVCD,LOFF,LDIA, 2LROA,LINV, 3LNAMES,CUTOFF,n1,n2,IBONDSS,NBO,ixb,ixs,XS,XB,ioo) call w36('Transferring the tensors ...') CALL IMPROVES(R,RS,NAT,NATS,ND,NDS, 1IBONDSB,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA, 2LROA,ALPHA,AA,G,ALPHAS,AAS,GS,LINV, 3LNAMES,CUTOFF,n1,n2,IBONDSS,NBO,ixb,ixs,XS,XB,ioo) c OPEN(2,FILE='FILE.FC') CALL WRITEFF(3*MX*3*MX,3*NAT,FF) CLOSE(2) call w36('File FILE.FC written ') c OPEN(2,FILE='AP.PRO') CALL WRITEFF(MX*MX,NAT,pmax) CLOSE(2) call w36('New AP.PRO written') c IF(LABS.OR.LVCD)CALL WRITETEN(MX,PB,AB,R,NAT) IF(LROA)THEN ICO=2 CALL TTTT(3*MX,ALPHA,AA,G,NAT,ICO,R) CALL WRITETTT(MX,NAT,ALPHA,AA,G) call w36('File FILE.TTT written ') ENDIF OPEN(2,FILE='FILE.X') CALL WRITERX(MX,NAT,R,ND,IBONDSB) call w36('File FILE.X written ') CLOSE(2) call report(' Program finished OK.') STOP END c c only find best atoms imax,jmax: SUBROUTINE IMPROVE(RB,RS,NAT,NATS,ND,NDS, 1NO,IBONDS,NBOB,LWR,LABS,LVCD,LOFF,LDIA,LROA, 2LINV,LNAMES, 3CUTOFF,n1,n2,IBONDSS,NBO,ixb,ixs,XS,XB,ioo) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=900,MXB=40,IOX=25) DIMENSION RB(3,NAT),RS(3,NATS),NBO(MX), 1IBONDSS(4,MX),ixb(MXB),ixs(MXB), 4IBONDS(4,MX),NBOB(MX),ND(MX),NDS(MX),XS(MXB,3),XB(MXB,3) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LINV, 1LNAMES COMMON/MATVAR/A(3,3),XST(MXB,3),XBT(MXB,3),TOL,ITU,IIT,FFU character*80 NAME(IOX) common/cname/aname(IOX,MX,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, 3nboname(IOX,MX),ibondsname(IOX,4,MX),NDNAME(IOX,MX) common/var1/FFB(3*MX*3*MX),FFS(3*MX*3*MX),pmax(MX*MX), 1imax(MX*MX),jmax(MX*MX),iomax(MX*MX),pmaxnew(MX*MX) common/opts/ratlim c 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(LINV)WRITE(3,*)' Smaller fragment as optical antipode' IF(LNAMES)WRITE(3,*)' Named small fragments' IF(CUTOFF.gt.0.0d0)WRITE(3,3999)CUTOFF 3999 format('Distance cutoff ',f12.2,' A') C c BIG LOOP OVER FRAGMENTS: if(LNAMES)then nloop=NO c c pre-order bond atoms according to the distance do 5356 IO=1,NO do 5357 IB=1,NATNAME(IO) bx=rname(IO,1,IB) by=rname(IO,2,IB) bz=rname(IO,3,IB) do 538 iab=1,NBONAME(IO,IB) iaa=IBONDSNAME(IO,iab,IB) db=(bx-rname(IO,1,iaa))**2+(by-rname(IO,2,iaa))**2+ 1 (bz-rname(IO,3,iaa))**2 do 538 iac=iab+1,NBONAME(IO,IB) ibb=IBONDSNAME(IO,iac,IB) dc=(bx-rname(IO,1,ibb))**2+(by-rname(IO,2,ibb))**2+ 1 (bz-rname(IO,3,ibb))**2 if(dc.lt.db.or. 1 (dc.eq.db.and.NBONAME(IO,ibb).lt.NBONAME(IO,iaa)))then IBONDSNAME(IO,iab,IB)=ibb IBONDSNAME(IO,iac,IB)=iaa iaa=ibb db=dc endif 538 continue 5357 continue 5356 continue else nloop=1 c c pre-order bond atoms according to the distance do 5355 IB=1,NATS bx=RS(1,IB) by=RS(2,IB) bz=RS(3,IB) do 535 iab=1,NBO(IB) iaa=IBONDSS(iab,IB) db=(bx-RS(1,iaa))**2+(by-RS(2,iaa))**2+(bz-RS(3,iaa))**2 do 535 iac=iab+1,NBO(IB) ibb=IBONDSS(iac,IB) dc=(bx-RS(1,ibb))**2+(by-RS(2,ibb))**2+(bz-RS(3,ibb))**2 if(dc.lt.db.or.(dc.eq.db.and.NBO(ibb).lt.NBO(iaa)))then IBONDSS(iab,IB)=ibb IBONDSS(iac,IB)=iaa iaa=ibb db=dc endif 535 continue 5355 continue endif c c pre-order bond atoms according to the distance do 345 IA=1,NAT ax=RB(1,IA) ay=RB(2,IA) az=RB(3,IA) do 34 iab=1,NBOB(IA) iaa=IBONDS(iab,IA) db=(ax-RB(1,iaa))**2+(ay-RB(2,iaa))**2+(az-RB(3,iaa))**2 do 34 iac=iab+1,NBOB(IA) ibb=IBONDS(iac,IA) dc=(ax-RB(1,ibb))**2+(ay-RB(2,ibb))**2+(az-RB(3,ibb))**2 if(dc.lt.db.or.(dc.eq.db.and.NBOB(ibb).lt.NBOB(iaa)))then IBONDS(iab,IA)=ibb IBONDS(iac,IA)=iaa iaa=ibb db=dc endif 34 continue 345 continue c c LOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOP do IO=1,nloop if(LWR)write(3,*)'Fragment loop ',IO c c load small fc, so that it is available in memory if(LNAMES)then NATS=NATNAME(IO) do 533 ia=1,NATS NDS(ia)=NDNAME(IO,ia) NBO(ia)=nboname(IO,ia) do 534 ib=1,NBO(ia) 534 IBONDSS(ib,ia)=ibondsname(IO,ib,ia) do 533 ix=1,3 533 RS(ix,ia)=rname(IO,ix,ia) 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 c DO 20 IA= istart,iend write(6,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) aix=RB(1,IA) aiy=RB(2,IA) aiz=RB(3,IA) kindia=ND(IA) nboia=NBOB(IA) c call startbedi(itu1,ixb,IA,RB,nboia,IBONDS,xb) c c if ioo>0, add a second layer of bonded atoms if(ioo.gt.0)call addsecond(2,itu1,ixb,RB,IBONDS,xb,NBOB,IA) c c Loops over the atoms in the fragment IO DO 32 IB=1,NATS c exclude different kinds: if(NDS(IB).ne.kindia)goto 32 c exclude atoms with different number of bonds: if(NBO(IB).ne.nboia)goto 32 bix=RS(1,IB) biy=RS(2,IB) biz=RS(3,IB) c c exclude different types do 536 iab=1,nboia 536 if(ND(IBONDS(iab,IA)).ne.NDS(IBONDSS(iab,IB)))goto 32 c c now IA,IB are same kind, are bounded to same kind and number of c atoms c c start bde on IB call startbedi(itu2i,ixs,IB,RS,nboia,IBONDSS,xs) c c if ioo>0, add a second layer of bonded atoms if(ioo.gt.0)call addsecond(2,itu2i,ixs,RS,IBONDSS,xs,NBO,IB) c if(itu2i.ne.itu1)goto 32 c DO 201 JA=1,IA IF(.NOT.LOFF.AND.IA.NE.JA)GOTO 201 ajx=RB(1,JA) ajy=RB(2,JA) ajz=RB(3,JA) dista=sqrt((aix-ajx)**2+(aiy-ajy)**2+(aiz-ajz)**2) if(CUTOFF.gt.0.0d0)THEN if(dista.gt.CUTOFF) goto 201 endif c kindja=ND(JA) nboja=NBOB(JA) c c Add atoms on JA: call addbedja(itu,itu1,JA,IA,ixb,xb,RB,IBONDS,nboja) c if(ioo.gt.0)call addsecond(itu1+2,itu,ixb,RB,IBONDS,xb,NBOB,JA) c c try to find closest atom to enable transfer: Tx=0.5d0*(aix+ajx) Ty=0.5d0*(aiy+ajy) Tz=0.5d0*(aiz+ajz) call incla(itu,NAT,ixb,Tx,Ty,Tz,RB,ND,NDS,0,xb,ixs) c c find if IA,JA are covalently bonded, set cutoff fac=0.5d0 do 350 icc=1,nboia 350 if(IBONDS(icc,IA).eq.JA)goto 351 fac=1.0d0 351 rrlim=ratlim*fac c C Calculate geometric center: TBX=0.0d0 TBY=0.0d0 TBZ=0.0d0 DO 3411 ISAT=1,ITU TBX=TBX+XB(ISAT,1) TBY=TBY+XB(ISAT,2) 3411 TBZ=TBZ+XB(ISAT,3) fc=DBLE(ITU) TBX=TBX/fc TBY=TBY/fc TBZ=TBZ/fc c C Subtract the center: DO 3511 ISAT=1,ITU XBT(ISAT,1)=XB(ISAT,1)-TBX XBT(ISAT,2)=XB(ISAT,2)-TBY 3511 XBT(ISAT,3)=XB(ISAT,3)-TBZ c DO 33 JB=1,NATS if(NDS(JB).ne.kindja)goto 33 if(NBO(JB).ne.nboja)goto 33 bjx=RS(1,JB) bjy=RS(2,JB) bjz=RS(3,JB) distb=sqrt((bix-bjx)**2+(biy-bjy)**2+(biz-bjz)**2) if(distb+dista.lt.1.0d-3)then ratio=0.0d0 else ratio=abs(distb-dista)/(distb+dista) endif if(ratio.gt.rrlim)goto 33 c c exclude different types do 37 iab=1,nboja 37 if(ND(IBONDS(iab,JA)).ne.NDS(IBONDSS(iab,JB)))goto 33 c c transcript the chosen atoms: call addbedja(itu2,itu2i,JB,IB,ixs,xs,RS,IBONDSS,nboja) c if(ioo.gt.0)call addsecond(itu2i+2,itu2,ixs,RS,IBONDSS,xs,NBO,JB) c c include additional atoms Txs=0.5d0*(bix+bjx) Tys=0.5d0*(biy+bjy) Tzs=0.5d0*(biz+bjz) call incla(itu2,NATS,ixs,Txs,Tys,Tzs,RS,ND,NDS,1,xs,ixb) c if(itu2.ne.itu)goto 33 c c now (IA,IB),(JA,JB) are same kind, are bounded to same kind and number of c atoms, and are reasonably close, so try to overlap: C C Calculate geometric center: TSX=0.0d0 TSY=0.0d0 TSZ=0.0d0 DO 3422 ISAT=1,ITU TSX=TSX+XS(ISAT,1) TSY=TSY+XS(ISAT,2) 3422 TSZ=TSZ+XS(ISAT,3) TSX=TSX/fc TSY=TSY/fc TSZ=TSZ/fc c C Subtract the center: DO 35 ISAT=1,ITU XST(ISAT,1)=XS(ISAT,1)-TSX XST(ISAT,2)=XS(ISAT,2)-TSY 35 XST(ISAT,3)=XS(ISAT,3)-TSZ C IF(LINV)THEN DO 39 ISAT=1,ITU DO 39 J=1,3 39 XST(ISAT,J)=-XST(ISAT,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR) c c in case of failure give up IA..JA pair IF(IERR.EQ.1)goto 33 c pp=sqrt(FFU/ITU) c c Best trafo try........................................... if(pp.lt.pmaxnew(IA+(JA-1)*NAT))THEN pmaxnew(IA+(JA-1)*NAT)=pp pmaxnew(JA+(IA-1)*NAT)=pp imax(IA+(JA-1)*NAT)=IB jmax(IA+(JA-1)*NAT)=JB if(LNAMES)iomax(IA+(JA-1)*NAT)=IO endif 33 continue c JB c 201 continue c JA c 32 continue c IB c 20 continue c IA C c LOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOP enddo c return END c c actually do the transformation: SUBROUTINE IMPROVES(RB,RS,NAT,NATS,ND,NDS, 1IBONDS,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,LROA, 2ALPHA,AA,G,ALPHAS,AAS,GS,LINV,LNAMES, 3CUTOFF,n1,n2,IBONDSS,NBO,ixb,ixs,XS,XB,ioo) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=900,MXB=40,IOX=25) DIMENSION RB(3,NAT),RS(3,NATS),NBO(MX), 1ND(MX),NDS(MX), 3IBONDSS(4,MX),ixb(MXB),ixs(MXB),XS(MXB,3),XB(MXB,3), 4IBONDS(4,MX),NBOB(MX),PB(MX,3,3),AB(MX,3,3),PS(MX,3,3), 5AS(MX,3,3),ALPHA(3*MX,3,3),ALPHAS(3*MX,3,3), 6AA(3*MX,3,3,3),AAS(3*MX,3,3,3),G(3*MX,3,3),GS(3*MX,3,3) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LINV, 1LNAMES COMMON/MATVAR/A(3,3),XST(MXB,3),XBT(MXB,3),TOL,ITU,IIT,FFU character*80 NAME(IOX) common/cname/aname(IOX,MX,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, 3nboname(IOX,MX),ibondsname(IOX,4,MX),NDNAME(IOX,MX) common/var1/FFB(3*MX*3*MX),FFS(3*MX*3*MX),pmax(MX*MX), 1imax(MX*MX),jmax(MX*MX),iomax(MX*MX),pmaxnew(MX*MX) 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 c DO 20 IA= istart,iend write(6,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) aix=RB(1,IA) aiy=RB(2,IA) aiz=RB(3,IA) nboia=NBOB(IA) c c assign atomic bed as in IMPROVE(): call startbedi(itu1,ixb,IA,RB,nboia,IBONDS,xb) c c if ioo>0, add a second layer of bonded atoms if(ioo.gt.0)call addsecond(2,itu1,ixb,RB,IBONDS,xb,NBOB,IA) c DO 201 JA=1,IA IF(.NOT.LOFF.AND.IA.NE.JA)GOTO 201 ajx=RB(1,JA) ajy=RB(2,JA) ajz=RB(3,JA) dista=sqrt((aix-ajx)**2+(aiy-ajy)**2+(aiz-ajz)**2) if(CUTOFF.gt.0.0d0)THEN if(dista.gt.CUTOFF) goto 201 endif c pp=pmaxnew(IA+(JA-1)*NAT) if(lwr)write(3,*)'Probability new old:', 1 pp,pmax(IA+(JA-1)*NAT) if(1.00001d0*pp.gt.pmax(IA+(JA-1)*NAT))then if(lwr)write(3,*)'skipped' goto 201 endif pmax(IA+(JA-1)*NAT)=pp pmax(JA+(IA-1)*NAT)=pp c nboja=NBOB(JA) c c assign atomic bed as before in IMPROVE() call addbedja(itu,itu1,JA,IA,ixb,xb,RB,IBONDS,nboja) c if(ioo.gt.0)call addsecond(itu1+2,itu,ixb,RB,IBONDS,xb,NBOB,JA) c c try to find and add closest atom to enable transfer: Tx=0.5d0*(aix+ajx) Ty=0.5d0*(aiy+ajy) Tz=0.5d0*(aiz+ajz) call incla(itu,NAT,ixb,Tx,Ty,Tz,RB,ND,NDS,0,xb,ixs) c c transfer chosen pair IB=imax(IA+(JA-1)*NAT) JB=jmax(IA+(JA-1)*NAT) if(LNAMES)then IO=iomax(IA+(JA-1)*NAT) else IO=1 endif c if(LNAMES)then NATS=NATNAME(IO) NE=80 DO 1500 NE=80,1,-1 1500 if(NAME(IO)(NE:NE).ne.' ')goto 1600 1600 OPEN(2,FILE=NAME(IO)(1:NE)//'.fc') CALL READFF(3*3*MX*MX,3*NATS,FFS) CLOSE(2) do 533 ii=1,NATS NDS(ii)=NDNAME(IO,ii) NBO(ii)=nboname(IO,ii) do 534 ib=1,NBO(ii) 534 IBONDSS(ib,ii)=ibondsname(IO,ib,ii) do 533 ix=1,3 533 RS(ix,ii)=rname(IO,ix,ii) endif c bix=RS(1,IB) biy=RS(2,IB) biz=RS(3,IB) bjx=RS(1,JB) bjy=RS(2,JB) bjz=RS(3,JB) c c define bed on ib call startbedi(itu2,ixs,IB,RS,nboia,IBONDSS,xs) if(ioo.gt.0)call addsecond(2,itu2,ixs,RS,IBONDSS,xs,NBO,IB) c ituold=itu2 call addbedja(itu2,ituold,JB,IB,ixs,xs,RS,IBONDSS,nboja) if(ioo.gt.0)call addsecond(ituold+2,itu2,ixs,RS,IBONDSS,xs,NBO,JB) c c include additional atoms Txs=0.5d0*(bix+bjx) Tys=0.5d0*(biy+bjy) Tzs=0.5d0*(biz+bjz) call incla(itu2,NATS,ixs,Txs,Tys,Tzs,RS,ND,NDS,1,xs,ixb) c if(itu2.ne.itu)call report('itu2<>itu, programming error') c c now (IA,IB),(JA,JB) are same kind, are bounded to same kind and number of c atoms, and are reasonably close, so try to overlap: C if(LWR)then WRITE(3,3004)IA,JA,IO,pmax(IA+(JA-1)*NAT),ITU,(ixb(I),I=1,ITU) WRITE(3,3005)(ixs(i),i=1,itu2) if(LNAMES)write(3,3399)NAME(IO) 3399 FORMAT(A80) 3004 FORMAT('Pair',2i5,', fragment:',I3,', weight = ',F7.4,',', 1 I4,'-atom bed:',/, 2 ' Big : ',20I3) 3005 FORMAT( ' Small: ',20I3,/) endif c C Calculate geometric centes: TBX=0.0d0 TBY=0.0d0 TBZ=0.0d0 TSX=0.0d0 TSY=0.0d0 TSZ=0.0d0 DO 34 ISAT=1,ITU TBX=TBX+XB(ISAT,1) TBY=TBY+XB(ISAT,2) TBZ=TBZ+XB(ISAT,3) TSX=TSX+XS(ISAT,1) TSY=TSY+XS(ISAT,2) 34 TSZ=TSZ+XS(ISAT,3) fac=DBLE(ITU) TBX=TBX/fac TBY=TBY/fac TBZ=TBZ/fac TSX=TSX/fac TSY=TSY/fac TSZ=TSZ/fac c C Subtract the centers: DO 35 ISAT=1,ITU XST(ISAT,1)=XS(ISAT,1)-TSX XST(ISAT,2)=XS(ISAT,2)-TSY XST(ISAT,3)=XS(ISAT,3)-TSZ XBT(ISAT,1)=XB(ISAT,1)-TBX XBT(ISAT,2)=XB(ISAT,2)-TBY 35 XBT(ISAT,3)=XB(ISAT,3)-TBZ C IF(LINV)THEN DO 39 ISAT=1,ITU DO 39 J=1,3 39 XST(ISAT,J)=-XST(ISAT,J) ENDIF C C Calculate the transformation matrix: CALL DOMATRIX(IERR) c IF(LINV)THEN DO 50 I=1,3 DO 50 J=1,3 50 A(I,J)=-A(I,J) ENDIF c if(LWR)write(3,3009)IIT,TOL,FFU 3009 format(i6,' it., err = ',2f10.7) C Transform the force field: IF((LDIA.AND.IA.EQ.JA).OR.(LOFF.AND.IA.NE.JA))THEN if(LWR)write(3,3010) 3010 format('F i j old new:') DO 36 IX=1,3 INDIB=3*(IA-1)+IX ind1=3*NAT*(INDIB-1) JEND=3 IF(IA.EQ.JA)JEND=IX DO 36 JX=1,JEND INDJB=3*(JA-1)+JX indf=ind1+INDJB FNEW=0.0d0 DO 37 IIX=1,3 INDIS=3*(IB-1)+IIX ind1s=3*NATS*(INDIS-1) DO 37 JJX=1,3 INDJS=3*(JB-1)+JJX 37 FNEW=FNEW+A(IX,IIX)*FFS(ind1s+INDJS)*A(JX,JJX) FOLD=FFB(indf) IF(LWR)WRITE(3,3001)IX,JX,FOLD,FNEW 3001 format(1x,2i2,2f11.6) FFB(indf)=FNEW 36 FFB(INDIB+3*NAT*(INDJB-1))=FNEW ENDIF C C Transform the polarization tensors IF(IA.EQ.JA.AND.(LABS.OR.LVCD))THEN if(LWR)write(3,3011) 3011 format(5x,'APT',8x,'old /',8x,'new AAT old / new:') DO 40 IX=1,3 DO 40 JX=1,3 POLD=PB(IA,IX,JX) AOLD=AB(IA,IX,JX) 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(IO,IB,II,JJ)*AD2 ANEW=ANEW+aname(IO,IB,II,JJ)*AD2 else PNEW=PNEW+PS(IB,II,JJ)*AD2 ANEW=ANEW+AS(IB,II,JJ)*AD2 endif 41 continue PB(IA,IX,JX)=PNEW AB(IA,IX,JX)=ANEW 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 3006 FORMAT(i4,2I2,F11.6,' /',F11.6,F12.6,' /',F11.6) ENDIF c C Transform the ROA tensors IF(IA.EQ.JA.AND.LROA)THEN IAS=IB DO 46 IX=1,3 IIND=3*(IA-1)+IX DO 46 I=1,3 DO 46 J=1,3 ALPHAOLD=ALPHA(IIND,I,J) GOLD=G(IIND,I,J) ALPHANEW=0.0d0 GNEW=0.0d0 DO 47 IXX=1,3 AD0=A(IX,IXX) INDS=3*(IAS-1)+IXX DO 47 II=1,3 AD1=AD0*A(I,II) DO 47 JJ=1,3 AD2=AD1*A(J,JJ) if(LNAMES)then ALPHANEW=ALPHANEW+alphaname(IO,INDS,II,JJ)*AD2 GNEW=GNEW+gname(IO,INDS,II,JJ)*AD2 else ALPHANEW=ALPHANEW+ALPHAS(INDS,II,JJ)*AD2 GNEW=GNEW+GS(INDS,II,JJ)*AD2 endif 47 continue ALPHA(IIND,I,J)=ALPHANEW G(IIND,I,J)=GNEW 46 IF(LWR)WRITE(3,3006)IIND,I,J,ALPHAOLD,ALPHANEW, 1 GOLD,GNEW 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 AAOLD=AA(IIND,I,J,K) 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 if(LNAMES)then AANEW=AANEW+aaname(IO,INDS,II,JJ,KK)*AD2*A(K,KK) else AANEW=AANEW+AAS(INDS,II,JJ,KK)*AD2*A(K,KK) endif 49 continue AA(IIND,I,J,K)=AANEW 48 IF(LWR)WRITE(3,3007)IIND,I,J,K,AAOLD,AANEW 3007 FORMAT('A',4I3,F12.6,'/',F12.6) ENDIF c ROA c 201 continue c JA c 20 continue c IA 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) c READ(2,*) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > MX !') do 1 i=1,N READ(2,*)ND(i),(R(j,i),j=1,3),(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 c 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,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N0) 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 LNIND=N*(LN-1) 130 READ(2,17)(FCAR(LNIND+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 IIND=(I-1)*N DO 31 J=I,N 31 FCAR(IIND+J)=FCAR((J-1)*N+I) RETURN END c SUBROUTINE WRITEFF(N0,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N0) c N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N DO 130 LN=N1,N LNIND=(LN-1)*N 130 WRITE(2,17)LN,(FCAR(LNIND+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,FFU c 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=0.0000001d0 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 FFU=ABS(Y(IHI))+ABS(Y(ILO)) IF(FFU.LT.1.d-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)RETURN 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 INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MXB=40) COMMON/MATVAR/A(3,3),XS(MXB,3),XB(MXB,3),RTOL,NAT,ITER,FFU DIMENSION ANG(3) 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 R=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) XB1=XB(I,1) XB2=XB(I,2) XB3=XB(I,3) 1 R=R+(TX-XB1)**2+(TY-XB2)**2+(TZ-XB3)**2 FU=R 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*)) call w36(RS) WRITE(3,3001) WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') CLOSE(3) CLOSE(2) STOP END c SUBROUTINE READTEN(N0,P,A,LAPT,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=900) 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) call w36('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 call w36('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) call w36('VCD parameters read in') C Distributed-origin gauge: 1 CALL VCDD0(AI,PV,P,R,NAT) call w36('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 call w36('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=900) 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 call w36('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) call w36('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=900) 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,*)NAT 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) 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) OPEN(2,FILE='FILE.TTT') 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,3F13.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) RETURN END SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R) C Transforms A and G to local origins (ICO=1) or 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=900) 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)call w36('G transformed into atomic origins') IF(ICO.EQ.2)call w36('G transformed into common origin') IF(ICO.EQ.3)call w36('G set to zero') C 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)call w36('A transformed into atomic origins') IF(ICO.EQ.2)call w36('A transformed into common origin') IF(ICO.EQ.3)call w36('A set to zero') RETURN END c SUBROUTINE readopt(IOX,NO, 1LWR,LABS,LVCD,LOFF,LDIA,LROA,LINV, 2LAPT,LALF,LRDO,LNAMES,NAME,CUTOFF, 2n1,n2,ioo) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*7 STR CHARACTER*80 NAME(IOX) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LINV,LAPT, 1LALF,LRDO,LNAMES common/opts/ratlim c ioo=0 ratlim=0.25d0 n1=0 n2=0 LRDO=.FALSE. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LOFF=.TRUE. LDIA=.TRUE. CUTOFF=-1.0d5 LROA=.FALSE. LINV=.FALSE. LNAMES=.FALSE. 9 READ(2,2000)STR 2000 FORMAT(A7) WRITE(3,2000)STR IF(STR(1:3).EQ.'IOO')THEN READ(2,*)ioo GOTO 9 ENDIF IF(STR(1:4).EQ.'N1N2')THEN READ(2,*)n1,n2 GOTO 9 ENDIF IF(STR(1:6).EQ.'LNAMES')THEN READ(2,*)LNAMES if(LNAMES)then READ(2,*)NO IF(NO.GT.IOX)CALL REPORT(' NO > IOX !') DO 2 I=1,NO READ(2,*)IC IF(IC.NE.I)CALL REPORT('Error in reading ') 2 READ(2,2900)NAME(I) 2900 FORMAT(a80) c c Control Output: DO 4 IO=1,NO WRITE(3,2002)IO 2002 FORMAT(' Overlap number ',I5) 4 WRITE(3,2900)NAME(I) endif GOTO 9 ENDIF IF(STR(1:4).EQ.'LINV')THEN READ(2,*)LINV GOTO 9 ENDIF IF(STR(1:4).EQ.'LAPT')THEN READ(2,*)LAPT GOTO 9 ENDIF IF(STR(1:3).EQ.'LWR')THEN READ(2,*)LWR GOTO 9 ENDIF IF(STR(1:4).EQ.'LALF')THEN READ(2,*)LALF GOTO 9 ENDIF IF(STR(1:4).EQ.'LRDO')THEN READ(2,*)LRDO GOTO 9 ENDIF IF(STR(1:4).EQ.'LROA')THEN READ(2,*)LROA GOTO 9 ENDIF IF(STR(1:4).EQ.'LOFF')THEN READ(2,*)LOFF GOTO 9 ENDIF IF(STR(1:4).EQ.'LDIA')THEN READ(2,*)LDIA GOTO 9 ENDIF IF(STR(1:4).EQ.'LABS')THEN READ(2,*)LABS GOTO 9 ENDIF IF(STR(1:4).EQ.'LVCD')THEN READ(2,*)LVCD GOTO 9 ENDIF IF(STR(1:6).EQ.'CUTOFF')THEN READ(2,*)CUTOFF GOTO 9 ENDIF IF(STR(1:3).EQ.'END')THEN close(2) RETURN ENDIF call w36(STR) call report('Unknown option') END c subroutine w36(s) character*(*)s n=len(s) write(6,600)(s(i:i),i=1,n) write(3,600)(s(i:i),i=1,n) 600 format(1x,79A1) return end c subroutine addbedja(itu,itu1,JA,IA,ixb,xb,RB,IBONDS,nboja) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (MX=900,MXB=40) dimension ixb(MXB),xb(MXB,3),RB(3,MX),IBONDS(4,MX) c itu=itu1 if(JA.ne.IA)then c c always add JA, even if defined, to avoid ambiguity itu=itu+1 if(itu.gt.MXB)call report('itu>MXB!') ixb(itu)=JA do 343 ix=1,3 343 xb(itu,ix)=RB(ix,JA) c do 344 ii=1,nboja jjb=IBONDS(ii,JA) do 345 jj=1,itu1 345 if(ixb(jj).eq.jjb)goto 344 itu=itu+1 if(itu.gt.MXB)call report('itu>MXB!') ixb(itu)=jjb do 346 ix=1,3 346 xb(itu,ix)=RB(ix,jjb) 344 continue endif return end c subroutine incla(itu2,NATS,ixs,Txs,Tys,Tzs,RS,ND,NDS,ic,xs,ixb) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (MX=900,MXB=40) dimension ixs(MXB),xs(MXB,3),RS(3,MX),ND(MX),NDS(MX),ixb(MXB) c 600 if(itu2.lt.3)then dmin=1.0d9 imin=0 do 447 ii=1,NATS do 448 jj=1,itu2 448 if(ii.eq.ixs(jj))goto 447 dd=(Txs-RS(1,ii))**2+(Tys-RS(2,ii))**2+(Tzs-RS(3,ii))**2 if(dd.lt.dmin.and.(ic.eq.0.or.ND(ixb(itu2+1)).eq.NDS(ii)))then dmin=dd imin=ii endif 447 continue if(imin.gt.0)then itu2=itu2+1 ixs(itu2)=imin do 449 ix=1,3 449 xs(itu2,ix)=RS(ix,imin) goto 600 endif endif return end c subroutine startbedi(itu1,ixb,IA,RB,nboia,IBONDS,xb) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (MX=900,MXB=40) dimension ixb(MXB),xb(MXB,3),RB(3,MX),IBONDS(4,MX) c itu1=1 ixb(1)=IA do 342 ix=1,3 342 xb(1,ix)=RB(ix,IA) do 341 iac=1,nboia itu1=itu1+1 if(itu1.gt.MXB)then write(3,*)IA,nboia,MXB do ii=1,itu1 write(3,*)ii,ixb(ii) enddo call report('itu1>MXB!') endif iy=IBONDS(iac,IA) ixb(itu1)=iy do 341 ix=1,3 341 xb(itu1,ix)=RB(ix,iy) return end c subroutine addsecond(istart,itu1,ixb,RB,IBONDS,xb,NB,i0) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) parameter (MX=900,MXB=40) dimension ixb(MXB),xb(MXB,3),RB(3,MX),IBONDS(4,MX),NB(MX) c itu0=itu1 do 1 ii=istart,itu0 ia=ixb(ii) do 2 jj=1,NB(ia) ib=IBONDS(jj,ia) do 3 ic=1,itu1 3 if(ixb(ic).eq.ib)goto 2 itu1=itu1+1 if(itu1.gt.MXB)then write(3,*)MXB,i0 do i5=1,itu1 write(3,*)i5,ixb(i5) enddo call report('itu1>MXB!') endif do 5 ix=1,3 5 xb(itu1,ix)=RB(ix,ib) ixb(itu1)=ib 2 continue 1 continue c c order according to the distance to i0 xa=RB(1,i0) ya=RB(2,i0) za=RB(3,i0) do 6 i=itu0+1,itu1 di=(xa-xb(i,1))**2+(ya-xb(i,2))**2+(za-xb(i,3))**2 do 7 j=i+1,itu1 dj=(xa-xb(j,1))**2+(ya-xb(j,2))**2+(za-xb(j,3))**2 if(dj.lt.di)then it=ixb(i) ixb(i)=ixb(j) ixb(j)=it do 8 ix=1,3 at=xb(i,ix) xb(i,ix)=xb(j,ix) 8 xb(j,ix)=at endif 7 continue 6 continue c return end