PROGRAM CCTNBIG c version for giant molecules and force field, c the final FF matrix is compressed according to flim cutoff IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) integer*4 okind PARAMETER (MX=20000,IOX=2000,MXS=500) c MX maximum number of atoms c MXB maximum number of overlaped atoms for each constant c IOX max # of overlaps DIMENSION R(3,MX),IND(MX),ND(MX), 1NDS(MX),RS(3,MX),FFS(3*MXS,3*MXS), 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), 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),QS(3*MXS,6),QB(3*MX,6) LOGICAL OPT,LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LINV,LSELECT, 1LAPT,LALF,LRDO,LNAMES,LPOLYMER,LHALTONERROR,LRAM,LQQ logical lxfile CHARACTER*80 NAME(IOX),xname,basicname,fname,nname,tname DIMENSION NDNAME(IOX,MXS) common/cname/asname(IOX,MXS,3,3), 1psname(IOX,MXS,3,3),aaname(IOX,3*MXS,3,3,3),gname(IOX,3*MXS,3,3), 2alphaname(IOX,3*MXS,3,3),rname(IOX,3,MXS),natname(IOX), 3qname(IOX,3*MXS,6),NAME c 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-2001',/,/, 9 ' INPUT FILES:',/, 9 ' Target: BIG.X(YZ) coordinates',/, 1 ' BIG.TEN APT and AAT tensors',/, 3 ' BIG.TTT ROA tensors',/, 3 ' Source: SMALL.X(YZ) coordinates',/, 2 ' SMALL. FC force field',/, 2 ' SMALL.TEN APTs and AATs',/, 2 ' SMALLQQ.TEN quad der.',/, 4 ' SMALL.TTT ROA tensors',/, 4 ' (or named .x,.fc, .ten and .ttt)',/, 4 ' Options: CTT.INP atom assignment',/,/, 5 ' OUTPUT:',/, 5 ' FILE.X(YZ) same as BIG.XYZ',/, 7 ' AFILE. FC improved BIG.FC,condensed',/, 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: c .X and .XYZ formates are allowed: 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 OPEN(2,FILE='BIG.XYZ',STATUS='OLD') CALL READR(MX,NAT,R,ND,IBONDSB,NBOB) 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(MX,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, 3flim,LQQ) 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 READRXS(MX,NATS,RS,NDS,NBO) NATNAME(I)=NATS do 13 ia=1,nat NDNAME(I,ia)=NDS(ia) do 13 ix=1,3 13 rname(I,ix,ia)=RS(ix,ia) 11 WRITE(3,*)NAT, ' atoms in the fragment' else if(lxfile)then OPEN(2,FILE='SMALL.X',STATUS='OLD') CALL READRXS(MX,NATS,RS,NDS,NBO) else OPEN(2,FILE='SMALL.XYZ',STATUS='OLD') CALL READRS(MX,NATS,RS,NDS,NBO) endif CLOSE(2) WRITE(3,*)NATS, ' atoms in the small molecule' 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 WRITE(3,3003) WRITE(6,3003) 3003 FORMAT('Initial force constants set to zero') c c Special section for named fragments: c nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn 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 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 c c just check here, if it is possible to load it: OPEN(2,FILE=fname,STATUS='OLD') CALL READFF(MXS,NATS,FFS) CLOSE(2) 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) IF(LQQ)THEN nname=basicname(1:NE)//'qq.ten' write(6,*)nname(1:NE+6) do 191 ia=1,NATS do 191 ix=1,3 191 RS(ix,ia)=rname(i,ix,ia) call READQ(3*MXS,QS,nname) call TRQD(3*MXS,QS,RS,NATS,1,PS,MX) do 181 ia=1,3*NATS do 181 ix=1,6 181 qname(i,ia,ix)=QS(ia,ix) ENDIF ENDIF C IF(LROA.or.LRAM)THEN tname=basicname(1:NE)//'.ttt' write(6,*)tname(1:NE+4) OPEN(15,FILE=tname,STATUS='OLD') CALL READTTT(MX,ALPHAS,AAS,GS,NATS) CLOSE(15) IF(LALF)THEN WRITE(3,*)' A and G obtained from alpha as DO' CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,3,RS) ELSE CALL TTTT(3*MX,ALPHAS,AAS,GS,NATS,1,RS) ENDIF 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 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 OPEN(2,FILE='SMALL.FC',STATUS='OLD') CALL READFF(MXS,NATS,FFS) CLOSE(2) WRITE(3,3002)NATS WRITE(6,3002)NATS ENDIF IF(LABS.OR.LVCD)THEN OPEN(15,FILE='SMALL.TEN',STATUS='OLD') CALL READTEN(MX,PS,AS,LAPT,RS,NATS) CLOSE(15) IF(LQQ)THEN call READQ(3*MXS,QS,'SMALLQQ.TEN') call TRQD(3*MXS,QS,RS,NATS,1,PS,MX) call WRITEQ(3*MXS,QS,NATS,'LOCALQQ.TEN') ENDIF 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 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(6,*)' Transferring the tensors ...' CALL IMPROVE(FFS,R,RS,NAT,NATS,IND, 1INDBIG,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,flim,LQQ,QS,QB) c IF(LABS.OR.LVCD)CALL WRITETEN(MX,PB,AB,R,NAT) IF(LQQ)THEN call TRQD(3*MX,QB,R,NAT,-1,PB,MX) call WRITEQ(3*MX,QB,NAT,'FILEQQ.TEN') ENDIF IF(LROA.or.LRAM)THEN CALL TTTT(3*MX,ALPHA,AA,G,NAT,2,R) CALL WRITETTT(MX,NAT,ALPHA,AA,G) WRITE(3,*)' File FILE.TTT written ' WRITE(6,*)' File FILE.TTT written ' ENDIF if(lxfile)then OPEN(2,FILE='FILE.X') CALL WRITERX(MX,NAT,R,ND,IBONDSB) WRITE(3,*)' File FILE.X written ' WRITE(6,*)' File FILE.X written ' else OPEN(2,FILE='FILE.XYZ') CALL WRITER(MX,NAT,R,ND,IBONDSB) WRITE(3,*)' File FILE.XYZ written ' WRITE(6,*)' File FILE.XYZ written ' endif CLOSE(2) CLOSE(3) WRITE(6,*)' Program finished OK.' STOP END SUBROUTINE IMPROVE(FFS,RB,RS,NAT,NATS,IND,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,flim,LQQ,QS,QB) implicit none INTEGER*4 MX,MXB,IOX,MXS,nn PARAMETER (MX=20000,MXB=40,IOX=2000,MXS=500) INTEGER*4 okind,IND(MX),INDBIG(IOX,MX),j33(MX*3), 4IBONDS(4,MX),NBOB(MX),IOL(IOX),ISUM,I,IA,IAS,IEND,II, 1IIND,IIOII,IIX,INDIS,INDS,INDJS,IO,i6,i66,KX, 2IOIJ,IP,IPASS,IPO,ISTART,IWG,IX,IXX,J,JA,JJ,IERC, 3JSTART,JX,K,KK,N1,N2,NAT,NATIO,NATS,NE,NF,NO,INDJB,JJX, 4ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,MX),nats_o(IOX) real*8 FFS(3*MXS,3*MXS),RB(3,MX),RS(3,NATS), 1TIJ(3),OTIJ(3),PB(MX,3,3),a33(MX*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,POLD,AOLD,AD1,AD2,ALPHAOLD,GOLD,AD0,AAOLD,CTF,CUTOFF, 6RM(IOX,3,3),FFB(3,3*MX),ft,QS(3*MXS,6),QB(3*MX,6) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LRAM,LSTRICT,LINV,LSELECT, 1LNAMES,LHALTONERROR,LQQ real*8 A,XS,XB,TOL,RMS,flim 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,qname,QNEW integer*4 natname character*80 NAME(IOX) common/cname/aname(IOX,MXS,3,3), 1pname(IOX,MXS,3,3),aaname(IOX,3*MXS,3,3,3),gname(IOX,3*MXS,3,3), 2alphaname(IOX,3*MXS,3,3),rname(IOX,3,MXS),natname(IOX), 3qname(IOX,3*MXS,6),NAME 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(LQQ) WRITE(3,*)' Quadrupole 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)then write(3,*)'best overlap on the basis of center of mass distances' else write(3,*)'best overlap on the basis of RMS distances' endif IF(CUTOFF.gt.0.0d0)WRITE(3,3999)CUTOFF 3999 format('Distance cutoff ',f12.2,' A') ctf=CUTOFF**2 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 open(33,file='AFILE.FC',form='unformatted') DO 20 IA= istart,iend write(6,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) DO 361 IX=1,3 DO 361 JA=1,3 DO 361 JX=1,3 361 FFB(IX,3*JA-3+JX)=0.0d0 DO 21 JA=IA,NAT if(LWR)write(3,*) if(LWR)write(3,*)'Pair ',IA,JA IF(.NOT.LOFF.AND.IA.NE.JA)GOTO 21 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 21 endif endif C C Calculate center of i/j radiusvectors: DO 22 IX=1,3 22 TIJ(IX)=0.5d0*(RB(IX,IA)+RB(IX,JA)) C C Look at all overlaps io, find that which fits best the i-j pair: DISTM=1.0d5 NF=0 DO 211 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 231 IX=1,3 231 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.eq.1)then call CRMS(IPASS,NAT,INDBIG,IOX,MX,IO, 1 IA,JA,NBOB,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 c c Use distances: else 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 211 CONTINUE c IF(NF.EQ.0)THEN GOTO 21 ELSE 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 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 21 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 c retrieve FF of fragment IPO NE=80 DO 1500 NE=80,1,-1 1500 if(NAME(IPO)(NE:NE).ne.' ')goto 1600 1600 OPEN(2,FILE=NAME(IPO)(1:NE)//'.fc') CALL READFF(MXS,NATS,FFS) CLOSE(2) ENDIF DO 36 IX=1,3 JSTART=1 IF(IA.EQ.JA)JSTART=IX DO 36 JX=JSTART,3 INDJB=3*JA-3+JX if(ISUM.eq.1)FFB(IX,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)*FFS(INDIS,INDJS)*A(JX,JJX) IF(LWR)WRITE(3,3001)IA,JA,IX,JX,0.0d0,FNEW,IIT,TOL,WF(IP) 3001 format(4i5,2f11.5,i5,F10.7,f7.4) if(IA.EQ.2.and.IX.EQ.1)write(6,*)'kk',JA,JX,IP,FNEW,WF(IP), 1 FFB(IX,INDJB) 36 FFB(IX,INDJB)=FNEW*WF(IP)+FFB(IX,INDJB) 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 quadrupole derivatives IF(IA.EQ.JA.AND.LQQ)THEN IAS=INDBIG(IPO,IA) DO 60 IX=1,3 DO 60 JX=1,3 DO 60 KX=JX,3 if(JX.eq.1.and.KX.eq.1)i6=1 if(JX.eq.2.and.KX.eq.2)i6=2 if(JX.eq.3.and.KX.eq.3)i6=3 if(JX.eq.1.and.KX.eq.2)i6=4 if(JX.eq.1.and.KX.eq.3)i6=5 if(JX.eq.2.and.KX.eq.3)i6=6 IF(ISUM.EQ.1)QB(3*(IA-1)+IX,i6)=0.0d0 QNEW=0.0d0 DO 61 II=1,3 AD1=A(IX,II) DO 61 JJ=1,3 AD2=AD1*A(JX,JJ) DO 61 KK=1,3 if(JJ.eq.1.and.KK.eq.1)i66=1 if(JJ.eq.2.and.KK.eq.2)i66=2 if(JJ.eq.3.and.KK.eq.3)i66=3 if(JJ.eq.1.and.KK.eq.2)i66=4 if(JJ.eq.2.and.KK.eq.1)i66=4 if(JJ.eq.1.and.KK.eq.3)i66=5 if(JJ.eq.3.and.KK.eq.1)i66=5 if(JJ.eq.2.and.KK.eq.3)i66=6 if(JJ.eq.3.and.KK.eq.2)i66=6 if(LNAMES)then QNEW=QNEW+qname(IPO,3*(IAS-1)+II,i66)*AD2*A(KX,KK) else QNEW=QNEW+QS(3*(IAS-1)+II,i66)*AD2*A(KX,KK) endif 61 continue QB(3*(IA-1)+IX,i6)=QNEW*WF(IP)+QB(3*(IA-1)+IX,i6) 60 IF(LWR)WRITE(3,3006)IA,IX,JX,0.0d0,QNEW,0.0d0,0.0d0,WF(IP) ENDIF C C Transform the ROA and Raman tensors IF(IA.EQ.JA.AND.(LROA.or.LRAM))THEN IAS=INDBIG(IPO,IA) 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) C Zero out at the start of the summation IF(ISUM.EQ.1)THEN ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 ENDIF 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(IPO,INDS,II,JJ)*AD2 GNEW=GNEW+gname(IPO,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)=ALPHA(IIND,I,J)+ALPHANEW*WF(IP) G(IIND,I,J)=G(IIND,I,J)+GNEW*WF(IP) IF (.NOT.LRAM)ALPHA(IIND,I,J)=ALPHAOLD IF (.NOT.LRAM)ALPHANEW =ALPHAOLD IF (.NOT.LROA) G(IIND,I,J)=GOLD IF (.NOT.LROA) GNEW=GOLD 46 IF(LWR)WRITE(3,3006)IIND,I,J,ALPHAOLD,ALPHANEW, 1 GOLD,GNEW,WF(IP) 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) 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 if(LNAMES)then AANEW=AANEW+aaname(IPO,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*WF(IP)+AA(IIND,I,J,K) IF (.NOT.LROA)AA(IIND,I,J,K)=AAOLD IF (.NOT.LROA)AANEW=AAOLD 48 IF(LWR)WRITE(3,3007)IIND,I,J,K,AAOLD,AANEW 3007 FORMAT('A',4I3,F12.6,'/',F12.6) ENDIF 45 CONTINUE ENDIF C 21 CONTINUE c JA c c write non-zero elements do 20 IX=1,3 nn=0 do 23 JA=IA,NAT JSTART=1 IF(IA.EQ.JA)JSTART=IX do 23 JX=JSTART,3 ft=FFB(IX,3*JA-3+JX) if(dabs(ft).gt.flim)then nn=nn+1 j33(nn)=3*JA-3+JX if(3*IA-3+IX.eq.3*JA-3+JX)ft=ft/2.0d0 a33(nn)=ft endif 23 continue 20 write(33)nn,(j33(ii),ii=1,nn),(a33(ii),ii=1,nn) c IA c close(33) write(6,*)'FF writen into AFILE.FC' C RETURN C END SUBROUTINE READR(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) N=0 1 N=N+1 READ(2,*)ND(N),R(1,N),R(2,N),R(3,N),(IBONDS(IB,N),IB=1,4) NB(N)=0 DO 3 IB=1,4 3 IF(IBONDS(IB,N).NE.0)NB(N)=IB IF(ND(N).EQ.9999)GOTO 2 GOTO 1 2 N=N-1 IF(N.GT.N0)CALL REPORT(' N > MX !') RETURN END c SUBROUTINE READRS(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) N=0 1 N=N+1 READ(2,*)ND(N),R(1,N),R(2,N),R(3,N) NB(N)=0 IF(ND(N).EQ.9999)GOTO 2 GOTO 1 2 N=N-1 IF(N.GT.N0)CALL REPORT(' N > MX !') 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 WRITER(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') 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),CHARGE WRITE(2,2001)9999,0.0,0.0,0.0,0,0,0,0,0.0 2001 FORMAT(I5,3F12.6,4I5,F8.4) 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 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 READQ(N0,Q,fn) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N0,6) character*(*) fn C C Read quadrupole derivatives open(15,file=fn) read(15,*)NAT read(15,*) DO 10 L=1,NAT*3 10 READ (15,*)idum,idum, (Q(L,I),I=1,6) close(15) c xx yy zz xy xz yz in atomic units WRITE(6,*)' Quadrupole derivatives read - in ' WRITE(6,*)NAT, ' atoms' return end SUBROUTINE WRITEQ(N0,Q,NAT,fn) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N0,6) character*(*) fn open(15,file=fn) write(15,*)NAT write(15,*) DO 10 L=1,NAT DO 10 J=1,3 10 write(15,1515)L,J,(Q(3*(L-1)+J,I),I=1,6) 1515 format(i5,i2,6f12.4) close(15) c xx yy zz xy xz yz in atomic units WRITE(6,*)' Quadrupole derivatives written ' return end C SUBROUTINE TRQD(N0,Q,X,NAT,ic,P,MX) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Q(N0,6),X(3,*),P(MX,3,3) BOHR=0.52917705993d0 if(ic.eq.0)then write(6,*)'no change in quadrupole derivatives' return else if(ic.eq.-1)then sign= 0.5d0/BOHR write(6,*)'quadrupole derivatives shifted back' else sign=-0.5d0/BOHR write(6,*)'quadrupole derivatives put on atoms' endif endif DO 10 L=1,NAT DO 10 J=1,3 LL=3*(L-1)+J Q(LL,1)=Q(LL,1)+sign*(X(1,L)*P(L,J,1)+P(L,J,1)*X(1,L)) Q(LL,2)=Q(LL,2)+sign*(X(2,L)*P(L,J,2)+P(L,J,2)*X(2,L)) Q(LL,3)=Q(LL,3)+sign*(X(3,L)*P(L,J,3)+P(L,J,3)*X(3,L)) Q(LL,4)=Q(LL,4)+sign*(X(1,L)*P(L,J,2)+P(L,J,1)*X(2,L)) Q(LL,5)=Q(LL,5)+sign*(X(1,L)*P(L,J,3)+P(L,J,1)*X(3,L)) 10 Q(LL,6)=Q(LL,6)+sign*(X(2,L)*P(L,J,3)+P(L,J,2)*X(3,L)) c xx yy zz xy xz yz in atomic units return end C SUBROUTINE READTEN(N0,P,A,LAPT,X,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=20000) 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=20000) 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(3I7) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(15,1501) (P(L,J,I),I=1,3),L 1501 FORMAT(3F17.8,I9) 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=20000) 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) 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 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=20000) 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 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(MXT,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,flim,LQQ) implicit none INTEGER*4 MXT,INDBIG,IOX,NO, 1NAT,n1,n2,okind,I,IA,IAB,IC,IEND,II,IO,ISMALL,IWG,JA,NDD REAL*8 CUTOFF,flim DIMENSION 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,LQQ c n1=0 n2=0 okind=0 flim=0.02d-6 LRDO=.FALSE. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LOFF=.TRUE. LQQ=.TRUE. LDIA=.TRUE. LHALTONERROR=.TRUE. CUTOFF=-1.0d5 LROA=.FALSE. LRAM=.FALSE. LSTRICT=.TRUE. LINV=.FALSE. LSELECT=.FALSE. LNAMES=.FALSE. IWG=1 9 READ(2,2000)STR 2000 FORMAT(A7) WRITE(3,2000)STR IF(STR(1:4).EQ.'N1N2') READ(2,*)n1,n2 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.'FLIM') READ(2,*)flim IF(STR(1:4).EQ.'LDIA') READ(2,*)LDIA IF(STR(1:4).EQ.'LABS') READ(2,*)LABS IF(STR(1:3).EQ.'LQQ') READ(2,*)LQQ 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 BACKSPACE 2 READ(2,2000)STRC if(STR.ne.STRC)goto 9 NO=0 c c older laborius 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) 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) 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,2004)(NDD(II),II=1,IEND-IA+1) 2004 FORMAT(20I3) 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 c 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 by default 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)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