PROGRAM ict implicit none INTEGER*4 MX,NAT,L,I,ip,jp,kp,k,ix,jx,ixp,jxp,ID,IG, 1ia,IIND,j,ja,IOX,IAB,IAX,ii,KX,n1,n2,nats,NE,IINDP PARAMETER (MX=1000,IOX=25) c MX maximum number of atoms c IOX max # of overlaps INTEGER*4 IBOND(4,MX),NB(MX),ND(MX),NDS(MX),NBS(MX), 1iii_bf(4,MX),iii_bten(4,MX),iii_bttt(4,MX), 1iii_sf(4,MX),iii_sten(4,MX),iii_sttt(4,MX), 1iiis(IOX,4,MX) REAL*8 FI(3*MX,3*MX),R(3,MX),FC(3*MX,3*MX),RS(3,MX), 1FSI(3*MX,3*MX),SUM,EPS, 1A(MX,3,3),P(MX,3,3),ALPHA(3*MX,3,3),AA(3*MX,3,3,3),G(3*MX,3,3), 1AI(MX,3,3),PI(MX,3,3),ALPHAI(3*MX,3,3),AAI(3*MX,3,3,3), 1GI(3*MX,3,3),u(3,3,MX), 1ASI(MX,3,3),PSI(MX,3,3),ALPHASI(3*MX,3,3),AASI(3*MX,3,3,3), 1GSI(3*MX,3,3) logical lxfile,lff,lten,lttt INTEGER*4 IND,INDBIG,NO,okind,IWG REAL*8 CUTOFF,BOHR DIMENSION IND(MX),INDBIG(IOX,MX) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LAPT, 1LALF,LRDO,LNAMES,LHALTONERROR CHARACTER*80 NAME(IOX),xname,basicname,nname,tname real*8 aname,pname,aaname,gname,alphaname,rname integer*4 natname,NDNAME(IOX,MX) 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 data BOHR/0.52917705993d0/ c WRITE(6,3000) 3000 FORMAT(/, 1 ' Internal Coordinate Tensor Transfer' ,/,/, 9 ' INPUT FILES: BIG.X coordinates' ,/, 3 ' CCT.OPT options+overlap' , 4 ' definition' ,/, 1 ' Optional: BIGI.FC force field' ,/, 1 ' BIGI.TEN APT and AAT tensors',/, 3 ' BIGI.TTT ROA tensors' ,/,/, 1 ' SMALLI.FC force field' ,/, 1 ' SMALLI.TEN APT and AAT tensors',/, 3 ' SMALLI.TTT ROA tensors' ,/,/, 1 ' OUTPUT FILES:' ,/, 1 ' Internal Tensors: INT.FC force field' ,/, 1 ' INT.TEN APT and AAT tensors',/, 3 ' INT.TTT ROA tensors' ,/, 1 ' Cartesian Tensors: FILE.FC force field' ,/, 1 ' FILE.TEN APT and AAT tensors',/, 3 ' FILE.TTT ROA tensors' ,/, 3 ' FILE.X geometry' ,/) c inquire(FILE='BIG.X',EXIST=lxfile) if(lxfile)then OPEN(2,FILE='BIG.X',STATUS='OLD') CALL READRX(MX,NAT,R,ND,IBOND,NB) CLOSE(2) WRITE(6,*)NAT, ' atoms in FILE.X' else call report('BIG.X does not exist') endif call readopt(MX,INDBIG,IOX,NO, 1NAT,LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LSTRICT, 2LAPT,LALF,LRDO,LNAMES,NAME,LHALTONERROR,CUTOFF, 2okind,NATS,n1,n2) if(.not.LOFF.and..not.LDIA.and..not.LROA.and..not.LVCD. 1and..not.LABS)call report('nothing to do') c INQUIRE(FILE='BIGI.FC',EXIST=lff) IF(lff)THEN OPEN(2,FILE='BIGI.FC') CALL READFF(MX,NAT,FI) call readatms(2,NAT,iii_bf,MX) CLOSE(2) WRITE(6,3001)NAT 3001 FORMAT(' Force field found in BIGI.FC') ELSE do 101 i=1,3*NAT do 101 j=1,3*NAT 101 FI(i,j)=0.0d0 WRITE(6,*)' Force field zeroed-out' ENDIF c INQUIRE(FILE='BIGI.TEN',EXIST=lten) IF(lten)THEN write(6,*)'Tensors BIGI.TEN found' OPEN(15,FILE='BIGI.TEN') CALL READTENI(MX,PI,AI,NAT) call readatms(15,NAT,iii_bten,MX) CLOSE(15) ELSE CALL ZEROTEN(MX,PI,AI,NAT) ENDIF C INQUIRE(FILE='BIGI.TTT',EXIST=lttt) IF(lttt)THEN write(6,*)'Tensors BIGI.TTT found' OPEN(15,FILE='BIGI.TTT') CALL READTTT(MX,ALPHAI,AAI,GI,NAT) call readatms(15,NAT,iii_bttt,MX) CLOSE(15) ELSE CALL ZEROTTT(MX,ALPHAI,AAI,GI,NAT) ENDIF c c if named fragments, load all geometries and local atom defs: 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,NBS) 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(LDIA.or.LOFF)THEN NE=80 DO 1500 NE=80,1,-1 1500 if(NAME(I)(NE:NE).ne.' ')goto 1600 1600 OPEN(2,FILE=NAME(I)(1:NE)//'.fc') CALL READFF(MX,NATS,FSI) call readatms(2,NATS,iii_sf,MX) CLOSE(2) do 81 ii=1,4 do 81 ia=1,NATS 81 iiis(I,ii,ia)=iii_sf(ii,ia) ENDIF IF(LABS.OR.LVCD)THEN nname=basicname(1:NE)//'.ten' write(6,*)nname(1:NE+4) OPEN(15,FILE=nname,STATUS='OLD') CALL READTENI(MX,PSI,ASI,NATS) call readatms(15,NATS,iii_sten,MX) CLOSE(15) do 18 ia=1,NATS do 18 ix=1,3 do 18 jx=1,3 aname(i,ia,ix,jx)=ASI(ia,ix,jx) 18 pname(i,ia,ix,jx)=PSI(ia,ix,jx) do 82 ii=1,4 do 82 ia=1,NATS 82 iiis(I,ii,ia)=iii_sten(ii,ia) ENDIF IF(LROA)THEN tname=basicname(1:NE)//'.ttt' write(6,*)tname(1:NE+4) OPEN(15,FILE=tname,STATUS='OLD') CALL READTTT(MX,ALPHASI,AASI,GSI,NATS) call readatms(15,NATS,iii_sttt,MX) CLOSE(15) do 20 iax=1,3*NATS do 20 ix=1,3 do 20 jx=1,3 alphaname(i,iax,ix,jx)=ALPHASI(iax,ix,jx) gname(i,iax,ix,jx)=GSI(iax,ix,jx) do 20 kx=1,3 20 aaname(i,iax,ix,jx,kx)=AASI(iax,ix,jx,kx) do 83 ii=1,4 do 83 ia=1,NATS 83 iiis(I,ii,ia)=iii_sttt(ii,ia) ENDIF 11 WRITE(6,*)NATS, ' atoms in the fragment' else OPEN(2,FILE='SMALL.X',STATUS='OLD') CALL READRXS(MX,NATS,RS,NDS,NBS) CLOSE(2) WRITE(6,*)NATS, ' atoms in the small molecule' c c Read tensors IF(LROA)then OPEN(15,FILE='SMALLI.TTT',status='old') CALL READTTT(MX,ALPHASI,AASI,GSI,NATS) call readatms(15,NATS,iii_sttt,MX) CLOSE(15) write(6,*)'SMALLI.TTT read' do 84 I=1,NO do 84 ii=1,4 do 84 ia=1,NATS 84 iiis(I,ii,ia)=iii_sttt(ii,ia) endif IF(LABS.or.LVCD)THEN OPEN(15,FILE='SMALLI.TEN',status='old') CALL READTENI(MX,PSI,ASI,NATS) call readatms(15,NATS,iii_sten,MX) CLOSE(15) write(6,*)'SMALLI.TEN read' do 85 I=1,NO do 85 ii=1,4 do 85 ia=1,NATS 85 iiis(I,ii,ia)=iii_sten(ii,ia) ENDIF IF(LDIA.or.LOFF)THEN if(NATS.eq.0)call report('NATS not set') OPEN(2,FILE='SMALLI.FC',STATUS='OLD') CALL READFF(MX,NATS,FSI) call readatms(2,NATS,iii_sf,MX) CLOSE(2) WRITE(6,*)'SMALLI.FC read' do 86 I=1,NO do 86 ii=1,4 do 86 ia=1,NATS 86 iiis(I,ii,ia)=iii_sf(ii,ia) ENDIF 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(6,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(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 CALL IMPROVE(FI,FSI,R,RS,NAT,NATS,IND, 1INDBIG,NO,IBOND,NB,LWR,PI,AI,LABS,LVCD,ASI,PSI,LOFF,LDIA, 2IWG,LROA,ALPHAI,AAI,GI,ALPHASI,AASI,GSI,LSTRICT, 3LNAMES,LHALTONERROR,CUTOFF,n1,n2,okind, 4iiis,iii_bf) OPEN(2,FILE='FILE.X') CALL WRITERX(MX,NAT,R,ND,IBOND) close(2) C IF(ldia.or.loff)THEN c Define transformation matrix for FF OPEN(2,FILE='INT.FC') CALL WRITEFF(MX,NAT,FI) call writeatms(NAT,iii_bf,MX) CLOSE(2) WRITE(6,*)' File INT.FC written ' WRITE(6,*)' Transforming FF to Cartesians' call dtm(iii_bf,NAT,u,MX,R) do 3 ia=1,NAT do 3 ja=1,NAT do 3 ix=1,3 do 3 jx=1,3 FC(3*(ia-1)+ix,3*(ja-1)+jx)=0.0d0 do 3 ixp=1,3 do 3 jxp=1,3 3 FC(3*(ia-1)+ix,3*(ja-1)+jx)=FC(3*(ia-1)+ix,3*(ja-1)+jx) 1 +FI(3*(ia-1)+ixp,3*(ja-1)+jxp)*u(ixp,ix,ia)*u(jxp,jx,ia) OPEN(2,FILE='FILE.FC') CALL WRITEFF(MX,NAT,FC) CLOSE(2) WRITE(6,*)' File FILE.FC written ' ENDIF IF(labs.or.lvcd)then c Define transformation matrix for FF OPEN(2,FILE='INT.TEN') CALL WRITETEN(MX,PI,AI,NAT) call writeatms(NAT,iii_bf,MX) CLOSE(2) WRITE(6,*)' File INT.TEN written ' WRITE(6,*)' Transforming APT,AAT to cartesians' call dtm(iii_bf,NAT,u,MX,R) DO 14 L=1,NAT DO 14 J=1,3 DO 14 I=1,3 P(L,J,I)=0.0d0 A(L,J,I)=0.0d0 DO 14 JP=1,3 DO 14 IP=1,3 A(L,J,I)=A(L,J,I)+AI(L,JP,IP)*u(JP,J,L)*u(IP,I,L) 14 P(L,J,I)=P(L,J,I)+PI(L,JP,IP)*u(JP,J,L)*u(IP,I,L) c c Transform to common origin: 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/BOHR WRITE(6,*)' AAT shifted to common origin' OPEN(2,FILE='FILE.TEN') CALL WRITETEN(MX,P,A,NAT) CLOSE(2) WRITE(6,*)' File FILE.TEN written ' endif IF(lroa)then OPEN(2,FILE='INT.TTT') CALL WRITETTT(MX,NAT,ALPHAI,AAI,GI) call writeatms(NAT,iii_bf,MX) close(2) WRITE(6,*)' File INT.TTT written ' WRITE(6,*)' Transforming ALPHA, GP, A to Cartesians' call dtm(iii_bf,NAT,u,MX,R) DO 15 L=1,NAT DO 15 IX=1,3 IIND=3*(L-1)+IX DO 15 I=1,3 DO 15 J=1,3 ALPHA(IIND,I,J)=0.0d0 G(IIND,I,J)=0.0d0 DO 15 IXP=1,3 IINDP=3*(L-1)+IXP DO 15 IP=1,3 DO 15 JP=1,3 ALPHA(IIND,I,J)=ALPHA(IIND,I,J) 1 +ALPHAI(IINDP,IP,JP)*u(IXP,IX,L)*u(IP,I,L)*u(JP,J,L) 15 G(IIND,I,J)=G(IIND,I,J) 1 +GI(IINDP,IP,JP)*u(IXP,IX,L)*u(IP,I,L)*u(JP,J,L) DO 16 L=1,NAT DO 16 IX=1,3 IIND=3*(L-1)+IX DO 16 I=1,3 DO 16 J=1,3 DO 16 K=1,3 AA(IIND,I,J,K)=0.0d0 DO 16 IXP=1,3 IINDP=3*(L-1)+IXP DO 16 IP=1,3 DO 16 JP=1,3 DO 16 KP=1,3 16 AA(IIND,I,J,K)=AA(IIND,I,J,K) 1 +AAI(IINDP,IP,JP,KP)*u(IXP,IX,L)*u(IP,I,L)*u(JP,J,L)*u(KP,K,L) CALL TTTT(3*MX,ALPHA,AA,G,NAT,2,R) OPEN(2,FILE='FILE.TTT') CALL WRITETTT(MX,NAT,ALPHA,AA,G) close(2) WRITE(6,*)' File FILE.TTT written ' ENDIF call report(' Program finished OK.') END c SUBROUTINE READRX(N0,N,R,ND,IBONDS,NB) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION R(3,N0),ND(N0),IBONDS(4,N0),NB(N0) READ(2,*) read(2,*)N IF(N.GT.N0)CALL REPORT(' N > MX !') do 1 i=1,N READ(2,*)ND(i),R(1,i),R(2,i),R(3,i),(IBONDS(IB,i),IB=1,4) NB(i)=0 DO 1 IB=1,4 1 IF(IBONDS(IB,i).NE.0)NB(i)=IB RETURN END SUBROUTINE 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 C SUBROUTINE REPORT(RS) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3000) STOP END c SUBROUTINE READTENI(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (MX=1000) DIMENSION P(N0,3,3),A(N0,3,3),AI(MX,3,3),AJ(MX,3,3) READ (15,*) NOAT if(NOAT.ne.NAT)call report('NAT<>NOAT') DO 10 L=1,NAT DO 10 J=1,3 10 READ (15,*) (P(L,J,I),I=1,3) 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,*) 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 RETURN END SUBROUTINE ZEROTEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) DO 10 L=1,NAT DO 10 J=1,3 DO 10 I=1,3 A(L,J,I)=0.0d0 10 P(L,J,I)=0.0d0 WRITE(6,*)' AAT, APT zeroed ' RETURN END C SUBROUTINE WRITETEN(N0,P,A,NAT) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION P(N0,3,3),A(N0,3,3) WRITE(2,1500) NAT,NAT-6,0 1500 FORMAT(3I5) DO 10 L=1,NAT DO 10 J=1,3 10 WRITE(2,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(2,1501) (A(L,J,I),I=1,3),L Z0=0.0d0 DO 230 L=1,NAT DO 230 J=1,3 230 WRITE(2,1501) (Z0,I=1,3),L DO 100 L=1,NAT DO 100 J=1,3 100 WRITE(2,1501) (P(L,J,I),I=1,3),L 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 ZEROTTT(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) DO 1 I=1,3 DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX DO 1 J=1,3 G(IIND,I,J)=0.0d0 1 ALPHA(IIND,I,J)=0.0d0 DO 3 I=1,3 DO 3 J=1,3 DO 3 L=1,NAT DO 3 IX=1,3 IIND=3*(L-1)+IX DO 3 K=1,3 3 A(IIND,I,J,K)=0.0d0 write(6,*)' Alpha, Gp, A zeroed' RETURN END 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) 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) 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) RETURN END SUBROUTINE TTTT(N,ALPHA,A,G,NAT,ICO,R) C Transforms A and G to local origins (ICO=1) c 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=1000) 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(6,*)' G transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' G transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' 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(6,*)' A transformed into atomic origins' IF(ICO.EQ.2)WRITE(6,*)' A transformed into common origin' IF(ICO.EQ.3)WRITE(6,*)' A set to zero' RETURN END c subroutine norm(v) implicit none real*8 v(*),n n=dsqrt(v(1)**2+v(2)**2+v(3)**2) v(1)=v(1)/n v(2)=v(2)/n v(3)=v(3)/n return end subroutine writeatms(NAT,iii,MX) integer*4 MX,NAT,iii(4,MX),ia do 4 ia=1,NAT 4 write(2,200)ia,(iii(ii,ia),ii=1,4) 200 format(5i5) return end subroutine readatms(io,NAT,iii,MX) integer*4 MX,NAT,iii(4,MX),ia,io,idumm do 1 ia=1,NAT 1 read(io,*)idumm,(iii(ii,ia),ii=1,4) return end SUBROUTINE readopt(MXT,INDBIG,IOX,NO, 1NAT,LWR,LABS,LVCD,LOFF,LDIA,IWG,LROA,LSTRICT, 2LAPT,LALF,LRDO,LNAMES,NAME,LHALTONERROR,CUTOFF, 2okind,NATS,n1,n2) implicit none INTEGER*4 MXT,INDBIG,IOX,NO,NATS,n1,n2, 1NAT,okind,I,IA,IAB,IC,IEND,II,IO,ISMALL,IWG,JA,NDD REAL*8 CUTOFF DIMENSION INDBIG(IOX,MXT),NDD(20) CHARACTER*7 STR CHARACTER*80 NAME(IOX) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT,LAPT, 1LALF,LRDO,LNAMES,LHALTONERROR c n1=0 n2=0 okind=1 LRDO=.FALSE. LALF=.FALSE. LWR=.FALSE. LAPT=.FALSE. LABS=.FALSE. LVCD=.FALSE. LOFF=.TRUE. LDIA=.TRUE. NATS=0 LHALTONERROR=.TRUE. CUTOFF=-1.0d5 LROA=.FALSE. LSTRICT=.TRUE. LNAMES=.FALSE. IWG=1 open(2,file='CCT.OPT') 9 READ(2,2000,end=99,err=99)STR 2000 FORMAT(A7) WRITE(6,2000)STR IF(STR(1:3).EQ.'OKI')THEN READ(2,*)okind GOTO 9 ENDIF IF(STR(1:4).EQ.'NATS')THEN READ(2,*)NATS GOTO 9 ENDIF IF(STR(1:3).EQ.'IWG')THEN READ(2,*)IWG GOTO 9 ENDIF IF(STR(1:6).EQ.'LNAMES')THEN READ(2,*)LNAMES GOTO 9 ENDIF IF(STR(1:4).EQ.'LAPT')THEN READ(2,*)LAPT GOTO 9 ENDIF IF(STR.EQ.'LSTRICT')THEN READ(2,*)LSTRICT 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:7).EQ.'LHALTON')THEN READ(2,*)LHALTONERROR GOTO 9 ENDIF IF(STR(1:6).EQ.'CUTOFF')THEN READ(2,*)CUTOFF GOTO 9 ENDIF IF(STR(1:6).EQ.'POLYME')THEN READ(2,*)NO WRITE(6,*)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(6,*)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(6,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(6,2004)(NDD(II),II=1,IEND-IA+1) 2004 FORMAT(20I3) IA=IEND GOTO 7 6 WRITE(6,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(6,*)'overlap ',io write(6,*)'big ',ia,ja write(6,*)'small ',ismall call report('umbiqous overlap') endif enddo endif enddo enddo GOTO 9 ENDIF c 99 close(2) RETURN END subroutine dtm(ii,NAT,u,MX,R) implicit none integer*4 MX,NAT,ii(4,MX),ia,i1,i2,j1,j2,ix,jx real*8 u(3,3,MX),R(3,MX),v1(3),v2(3),v3(3) do 1 ia=1,NAT i1=ii(1,ia) i2=ii(2,ia) j1=ii(3,ia) j2=ii(4,ia) if(i1.gt.0.and.i2.gt.0.and.j1.gt.0.and.j2.gt.0)then do 2 ix=1,3 v1(ix)=R(ix,i2)-R(ix,i1) 2 v2(ix)=R(ix,j2)-R(ix,j1) v3(1)=v1(2)*v2(3)-v1(3)*v2(2) v3(2)=v1(3)*v2(1)-v1(1)*v2(3) v3(3)=v1(1)*v2(2)-v1(2)*v2(1) call norm(v1) call norm(v3) v2(1)=v3(2)*v1(3)-v3(3)*v1(2) v2(2)=v3(3)*v1(1)-v3(1)*v1(3) v2(3)=v3(1)*v1(2)-v3(2)*v1(1) do 13 ix=1,3 u(1,ix,ia)=v1(ix) u(2,ix,ia)=v2(ix) 13 u(3,ix,ia)=v3(ix) else do 14 ix=1,3 do 14 jx=1,3 14 u(jx,ix,ia)=0.0d0 endif 1 continue return end c ============================================================ SUBROUTINE IMPROVE(FI,FFS,RB,RS,NAT,NATS,IND,INDBIG, 1NO,IBONDS,NBOB,LWR,PB,AB,LABS,LVCD,AS,PS,LOFF,LDIA,IWG,LROA, 2ALPHA,AA,G,ALPHAS,AAS,GS,LSTRICT,LNAMES, 3LHALTONERROR,CUTOFF,n1,n2,okind, 4iiis,iiib) implicit none INTEGER*4 MX,MXB,IOX PARAMETER (MX=1000,IOX=25,MXB=40) INTEGER*4 okind,IND(MX),INDBIG(IOX,MX), 4IBONDS(4,MX),NBOB(MX),IOL(IOX),ISUM,I,IA,IAS,IEND,II, 1IIND,IIOII,INDIB,INDIS,INDS,INDJS,IO,i4,iaa,ipp, 2IOIJ,IP,IPASS,IPO,ISTART,IWG,IX,J,JA,IERC, 3JSTART,JX,K,N1,N2,NAT,NATIO,NATS,NE,NF,NO,INDJB, 4ipass_o(IOX),ierc_o(IOX),itu_o(IOX),ind_o(IOX,MX),nats_o(IOX), 5iiis(IOX,4,MX),iiib(4,MX) real*8 FI(3*MX,3*MX),FFS(3*MX,3*MX),RB(3,MX),RS(3,NATS), 1TIJ(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,ALPHAOLD,GOLD,AAOLD,CTF,CUTOFF, 6RM(IOX,3,3) LOGICAL LWR,LABS,LVCD,LOFF,LDIA,LROA,LSTRICT, 1LNAMES,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 integer*4 natname 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 c IF(LDIA)WRITE(6,*)' Diagonal force field will be changed' IF(LOFF)WRITE(6,*)' Off-diagonals will be changed' IF(LABS)WRITE(6,*)' Eletric derivatives will be transformed' IF(LVCD)WRITE(6,*)' Magnetic derivatives will be transformed' IF(LROA)WRITE(6,*)' ROA parameters will be transformed' c IF(LSTRICT)WRITE(6,*)' Only transformed atoms in the beds' IF(IWG.EQ.1)WRITE(6,*)' Multiple overlaps will be weighted' IF(LNAMES)WRITE(6,*)' Named small fragments' IF(LHALTONERROR)WRITE(6,*)' Will stop on error' IF(.not.LHALTONERROR)WRITE(6,*)' Will not stop on error' if(okind.eq.0)then write(6,*)'best overlap on the basis of center of mass distances' else write(6,*)'best overlap on the basis of RMS distances' endif IF(CUTOFF.gt.0.0d0)WRITE(6,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(6,*)'first index ',n1,' - ',n2 endif DO 20 IA= istart,iend write(6,*)IA,'. atom of ',NAT open(56,file='progress') write(56,*)IA,'. atom of ',NAT close(56) DO 20 JA=IA,NAT if(LWR)write(6,*) if(LWR)write(6,*)'Pair ',IA,JA 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(6,3998)sqrt(dist) 3998 format('Skipped because of the distance,',f20.2,' A') goto 20 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 21 IO=1,NO C C Check, if overlap io contains ij: c 1)check, if transfered atoms are part IF(INDBIG(IO,IA).GT.0.AND.INDBIG(IO,JA).GT.0)THEN c 2)check, if atoms for coordinate definitions are part ipp=0 do 101 i4=1,4 do 101 iaa=1,NAT if(INDBIG(IO,iaa).eq.iiis(IO,i4,INDBIG(IO,IA)))then ipp=ipp+1 iiib(i4,IA)=iaa endif if(INDBIG(IO,iaa).eq.iiis(IO,i4,INDBIG(IO,JA)))then ipp=ipp+1 iiib(i4,JA)=iaa endif 101 continue IF(ipp.eq.8)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.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,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 ENDIF 21 CONTINUE c IF(NF.EQ.0)THEN GOTO 20 ELSE if(LWR)WRITE(6,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,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(6,*)IPASS,' passes' if(IERC.ne.0)goto 20 if(LWR)WRITE(6,3004)IPO,WF(IP),ITU,(IND(I),I=1,ITU) if(LWR)WRITE(6,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.and.LWR)write(6,3399)NAME(IPO) 3399 FORMAT(A80) C C C Transfer 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(MX,NATS,FFS) CLOSE(2) 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=FI(INDIB,INDJB) C Zero out at the start of the summation IF(ISUM.EQ.1)FI(INDIB,INDJB)=0.0d0 INDIS=3*INDBIG(IPO,IA)-3+IX INDJS=3*INDBIG(IPO,JA)-3+JX FNEW=FFS(INDIS,INDJS) IF(LWR)WRITE(6,3001)IA,JA, IX,JX,FOLD,FNEW,IIT,TOL,WF(IP) 3001 FORMAT(4I5,2F11.5,I5,F10.7,F7.4) FI(INDIB,INDJB)=FNEW*WF(IP)+FI(INDIB,INDJB) FI(INDJB,INDIB)=FI(INDIB,INDJB) 36 CONTINUE ENDIF C C Transform the polarization tensors 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 if(LNAMES)then PNEW=pname(IPO,IAS,IX,JX) ANEW=aname(IPO,IAS,IX,JX) else PNEW=PS(IAS,IX,JX) ANEW=AS(IAS,IX,JX) endif 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(6,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 Transform the ROA tensors IF(IA.EQ.JA.AND.LROA)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 INDS=3*(IAS-1)+IX if(LNAMES)then ALPHANEW=alphaname(IPO,INDS,I,J) GNEW=gname(IPO,INDS,I,J) else ALPHANEW=ALPHAS(INDS,I,J) GNEW=GS(INDS,I,J) endif ALPHA(IIND,I,J)=ALPHA(IIND,I,J)+ALPHANEW*WF(IP) G(IIND,I,J)=G(IIND,I,J)+GNEW*WF(IP) 46 IF(LWR)WRITE(6,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 INDS=3*(IAS-1)+IX if(LNAMES)then AANEW=aaname(IPO,INDS,I,J,K) else AANEW=AAS(INDS,I,J,K) endif AA(IIND,I,J,K)=AANEW*WF(IP)+AA(IIND,I,J,K) 48 IF(LWR)WRITE(6,3007)IIND,I,J,K,AAOLD,AANEW 3007 FORMAT('A',4I3,F12.6,'/',F12.6) ENDIF 45 CONTINUE ENDIF C 20 CONTINUE c IA,JA C RETURN END c ============================================================ subroutine CRMS(IPASS,NAT,INDBIG,IOX,MX,IPO, 1IA,JA,NBOB,IC,IBONDS,LNAMES,natname,RS,RB,IND,LHALTONERROR, 2IERC,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 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(6,*)' ITU < 3 !' IF(ITU.LT.3.AND.IPASS.EQ.1.and.IC.eq.0)GOTO 666 if(ITU.LT.3)then write(6,*)'Only following atoms found in the big molecule:' write(6,*)(IND(iu) ,iu=1,ITU) if(LHALTONERROR)then CALL REPORT(' ITU < 3 !') else write(6,*)' ITU < 3 !' return endif endif 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 C Calculate the transformation matrix: CALL DOMATRIX(IERR) IF(IERR.EQ.1.AND.IPASS.EQ.1)THEN write(6,*)'domat error,try other pass' GOTO 666 ENDIF IF(IERR.EQ.1)then if(LHALTONERROR)then CALL REPORT('Program cannot continue') else write(6,*)'DOMAT error' return endif endif IERC=0 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 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(6,*)' 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 C 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 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