PROGRAM BUMATD c get diagonal intrinsic force constants IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,NAT0), 1B(MX,MX),NOB0,IS,UM(MX,MX) DIMENSION B0(MX,MX),LV(MX),IZ(NAT0) character*7 COO(MX),cty(7) integer*4 d1(MX),d2(MX) data cty/'STRETCH','BEND ','XXXXXXX','TORSION','XXXXXXX', 1 'XXXXXXX','INTERMO'/ integer*4,allocatable::hl(:) tol=1.0d-8 WRITE(*,*)' (SYMMETRIZED) B-MATRIX GENERATION' WRITE(*,*)' FILE.UMA ... INPUT FILE' OPEN(15,FILE='FILE.UMA') READ(15,*) READ(15,*) READ(15,*)NOAT,NOB0 WRITE(6,*)NOAT,' atoms',NOB0,' internal coordinates' IF(NOAT.GT.NAT0)call report(' Too many atoms') if(NOB0.gt.MX)call report('Too many internal coordinates') DO 201 I=1,NOAT 201 READ(15,1020)(X(J,I),J=1,3),IZ(I),IZ(I) 1020 FORMAT(3F12.6,2i4) NA=NOAT*3 B=0.0d0 if(NOB0.ne.NA-6.or.NOB0.ne.NA-5)write(6,61000) 61000 format(/,' Warning: Number of Coordinates <> 3*NAT-6/5.',/) IER=0 IS=0 NOB=0 277 NOB=NOB+1 READ(15,*)ICODE,N1,N2,N3,N4,N5,N6,FACTOR write(6,1050)NOB,ICODE,N1,N2,N3,N4,N5,N6,FACTOR 1050 FORMAT(8I6,F12.6) IF(DABS(FACTOR).LT.1.0D-6)FACTOR=1.0D0 if(ICODE.le.7)COO(NOB)=cty(ICODE) if(ICODE.eq.4)then d1(NOB)=n2 d2(NOB)=n3 else if(ICODE.eq.7)then if(NOB+5.gt.MX)call report('Too many internal coordinates') COO(NOB )='MSTRETC' COO(NOB+1)='MANGLE1' COO(NOB+2)='MANGLE2' COO(NOB+3)='MANGLE3' COO(NOB+4)='MANGLE4' COO(NOB+5)='MTORSIO' do 141 i=NOB,NOB+5 d1(i)=n3 141 d2(i)=n4 else if(ICODE.ge.8.and.ICODE.LE.13)then if(ICODE.eq. 8)COO(NOB)='MSTRETC' if(ICODE.eq. 9)COO(NOB)='ALPHA 1' if(ICODE.eq.10)COO(NOB)='ALPHA 2' if(ICODE.eq.11)COO(NOB)='BETA 1' if(ICODE.eq.12)COO(NOB)='BETA 2' if(ICODE.eq.13)COO(NOB)='TAU ' d1(NOB)=n3 d2(NOB)=n4 else if(ICODE.ge.14.and.ICODE.le.15)then if(ICODE.eq.14)COO(NOB)='HEL-TX ' if(ICODE.eq.15)COO(NOB)='HEL-TY ' if(ICODE.eq.16)COO(NOB)='HEL-TZ ' if(ICODE.eq.17)COO(NOB)='HEL-RX ' if(ICODE.eq.18)COO(NOB)='HEL-RY ' if(ICODE.eq.19)COO(NOB)='HEL-RZ ' if(ICODE.eq.20)COO(NOB)='HEL-STR' if(ICODE.eq.21)COO(NOB)='HEL-BRE' if(ICODE.eq.22)COO(NOB)='HEL-DEF' if(ICODE.eq.23)COO(NOB)='HEL-UWI' d1(NOB)=n1 d2(NOB)=n2 else if(ICODE.ge.24.and.ICODE.le.29)then if(ICODE.eq.24)COO(NOB)='TR - X ' if(ICODE.eq.25)COO(NOB)='TR - Y ' if(ICODE.eq.26)COO(NOB)='TR - Z ' if(ICODE.eq.27)COO(NOB)='ROT - X' if(ICODE.eq.28)COO(NOB)='ROT - Y' if(ICODE.eq.29)COO(NOB)='ROT - Z' d1(NOB)=n1 d2(NOB)=n2 else d1(NOB)=n1 d2(NOB)=n2 endif endif endif endif endif IC=IABS(ICODE) IF(IC.lt.1.or.IC.gt.29)then write(6,*)ICODE call report('Unknown CODE') endif GOTO (1, 2, 3, 4, 5, 6, 7, 8, 9,10,11,12,13,14,15,16,17,18,19,20, 1 21,22,23,24,25,26,27,28,29),IC 1 CALL BOST GO TO 60 2 CALL BEND GO TO 60 3 CALL OPLA GO TO 60 4 CALL TORS GO TO 60 5 IS=IS+1 6 CALL LIBE GO TO 60 c intermolecular, together: 7 CALL INTER(IZ,ICODE) GO TO 60 c intermolecular, separate: 8 CALL INTER(IZ,ICODE) GO TO 60 9 CALL INTER(IZ,ICODE) GO TO 60 10 CALL INTER(IZ,ICODE) GO TO 60 11 CALL INTER(IZ,ICODE) GO TO 60 12 CALL INTER(IZ,ICODE) GO TO 60 13 CALL INTER(IZ,ICODE) GO TO 60 c helix/cigar: 14 allocate(hl(N1)) CALL HELIX(IZ,hl) GO TO 60 15 CALL HELIX(IZ,hl) GO TO 60 16 CALL HELIX(IZ,hl) GO TO 60 17 CALL HELIX(IZ,hl) GO TO 60 18 CALL HELIX(IZ,hl) GO TO 60 19 CALL HELIX(IZ,hl) GO TO 60 20 CALL HELIX(IZ,hl) GO TO 60 21 CALL HELIX(IZ,hl) GO TO 60 22 CALL HELIX(IZ,hl) GO TO 60 23 CALL HELIX(IZ,hl) deallocate(hl) GO TO 60 c monomolecular: 24 CALL MONO(IZ,ICODE) GO TO 60 25 CALL MONO(IZ,ICODE) GO TO 60 26 CALL MONO(IZ,ICODE) GO TO 60 27 CALL MONO(IZ,ICODE) GO TO 60 28 CALL MONO(IZ,ICODE) GO TO 60 29 CALL MONO(IZ,ICODE) GO TO 60 60 IF(IER.NE.0)call report('Input error') IF(DABS(FACTOR-1.0D0).gT.1.0D-8)then IF(IC.EQ.5)GO TO 80 70 ISW=0 I=NOB GO TO 90 80 ISW=1 I=NOB0+IS 90 DO 100 J=1,NA 100 B(I,J)=B(I,J)*FACTOR IF(ISW.EQ.1)GO TO 70 endif if(NOB.lt.NOB0)goto 277 write(6,*)'Coordinate definition read' UM=0.0d0 READ(15,*)IUM IF(IUM.NE.0)THEN DO 2772 I=1,NA-6 READ(15,*)J SUM=0.0D0 DO 2774 IK=1,J READ(15,*)IJ,UM(I,IJ) LV(IK)=IJ 2774 SUM=SUM+UM(I,IJ)*UM(I,IJ) DO 2773 IK=1,J 2773 UM(I,LV(IK))=UM(I,LV(IK))/sqrt(SUM) 2772 CONTINUE WRITE(*,*)' USER-DEFINED U MATRIX' M=NA-6 ELSE DO 2771 I=1,NA 2771 UM(I,I)=1.0D0 M=NOB0 ENDIF CLOSE(15) IF(IUM.NE.0)THEN DO 1556 I=1,NA-6 DO 1556 J=1,NA SUM=0.0D0 DO 1557 K=1,NOB0 1557 SUM=SUM+UM(I,K)*B(K,J) 1556 B0(I,J)=SUM ELSE DO 1558 I=1,NOB0 DO 1558 J=1,NA 1558 B0(I,J)=B(I,J) ENDIF DO 155 I=1,NOB0 DO 155 J=1,NA 155 IF(DABS(B(I,J)).LT.tol) B(I,J) = 0.0D0 WRITE(*,*)' OUTPUT FILE: B.MAT' OPEN(9,FILE='B.MAT') WRITE(9,901)M 901 format(I6,' Internal coordinates') WRITE(9,903)NA 903 format(I6,' Cartesian coordinates') WRITE(9,902) 902 format(' Transformation matrix, non-zero elements:') do 150 I=1,M WRITE(9,9011)I 9011 format(I6,' internal:') nn=0 do 140 J=1,NA 140 if(DABS(B0(I,J)).gt.tol)nn=nn+1 write(9,9012)nn 9012 format(I6) nn=0 do 150 J=1,NA if(DABS(B0(I,J)).gt.tol)then nn=nn+1 ia=(J+2)/3 ix=J-3*(ia-1) WRITE(9,9100)ia,ix,B0(I,J) 9100 format(2I5,F15.8) endif 150 continue CLOSE(9) CALL CTI(d1,d2,COO) END SUBROUTINE BOST C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A BOND STRETCH C AS DEFINED BY WILSON. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RIJ(3),RJK(3),RJL(3),EIJ(3),EJK(3),EKL(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 30 IF(J.LE.0.OR.J.GT.NOAT)GO TO 30 IF(K.NE.0)GO TO 30 IF(L.NE.0)GO TO 30 IF(IX.NE.0)GO TO 30 IF(JX.NE.0)GO TO 30 DIJSQ=0.D0 DO 10 M=1,3 T=X(M,J)-X(M,I) RIJ(M)=T 10 DIJSQ=DIJSQ+T*T DIJ=DSQRT(DIJSQ) II=3*(I-1) JJ=3*(J-1) DO 20 M=1,3 T=RIJ(M) IF(DABS(T).LT.1.D-8)GO TO 20 T=T/DIJ B(NOB,II+M)=-T B(NOB,JJ+M)=T 20 CONTINUE RETURN 30 IER=1 WRITE(*,200) 200 FORMAT('BOST ERROR') RETURN END SUBROUTINE BEND C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS OF A VALENCE C ANGLE BENDING COORDINATE AS DEFINED BY WILSON. C I AND K ARE THE NUMBERS OF THE END ATOMS. C J IS THE NUMBER OF THE CENTRAL ATOM. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EKL(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 50 IF(J.LE.0.OR.J.GT.NOAT)GO TO 50 IF(K.LE.0.OR.K.GT.NOAT)GO TO 50 IF(L.NE.0)GO TO 50 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 50 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 50 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DOTJ=0.D0 DO 30 M=1,3 T=RJI(M)/DJI EJI(M)=T TP=RJK(M)/DJK EJK(M)=TP 30 DOTJ=DOTJ+T*TP IF(DABS(DOTJ).GT.0.99995D0)GO TO 60 SINJ=DSQRT(1.D0-DOTJ*DOTJ) II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) DO 40 M=1,3 SMI=DX*(DOTJ*EJI(M)-EJK(M))/(DJI*SINJ) IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=SMI SMK=DX*(DOTJ*EJK(M)-EJI(M))/(DJK*SINJ) IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=SMK SUM=SMI+SMK 40 IF(DABS(SUM).GE.1.D-8)B(NOB,JJ+M)=-SUM RETURN 50 IER=1 *C*C* ADDED BY KATH SPRING 1991 TO FIND INPUT ERROR * NEXT TWO LINES WRITE(*,200) 200 FORMAT('BEND ERROR') RETURN 60 IER=-1 WRITE(*,1000) RETURN C 1000 FORMAT(' I-J-K IS COLINEAR - USE LINEAR BEND') C END SUBROUTINE OPLA C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR AN OUT OF C PLANE WAGGING COORDINATE AS DEFINED BY DECIUS, MCINTOSH, MICHAELIAN C AND PETERSON. SUBROUTINE CODED BY M PETERSON, UNIV OF TORONTO. C I IS THE END ATOM (ATOM WAGGED WITH RESPECT TO J-K-L PLANE). C J IS THE APEX ATOM (ATOMS I, K AND L ARE ATTACHED TO J). C K AND L ARE THE ANCHOR ATOMS (DEFINE THE J-K-L PLANE). C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EJL(3) DIMENSION C1(3) C IF(I.LE.0.OR.I.GT.NOAT)GO TO 60 IF(J.LE.0.OR.J.GT.NOAT)GO TO 60 IF(K.LE.0.OR.K.GT.NOAT)GO TO 60 IF(L.LE.0.OR.L.GT.NOAT)GO TO 60 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 60 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 60 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DJLSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,L)-X(M,J) RJL(M)=T DJLSQ=DJLSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DJL=DSQRT(DJLSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 COSI=0.D0 COSK=0.D0 COSL=0.D0 DO 30 M=1,3 T=RJI(M)/DJI EJI(M)=T TP=RJK(M)/DJK EJK(M)=TP TPP=RJL(M)/DJL EJL(M)=TPP COSI=COSI+TP*TPP COSK=COSK+T*TPP 30 COSL=COSL+T*TP IF(DABS(COSI).GT.0.99995D0)GO TO 70 SINSIN=1.D0-COSI*COSI SINI=DSQRT(SINSIN) C1(1)=EJK(2)*EJL(3)-EJK(3)*EJL(2) C1(2)=EJK(3)*EJL(1)-EJK(1)*EJL(3) C1(3)=EJK(1)*EJL(2)-EJK(2)*EJL(1) DOT=EJI(1)*C1(1)+EJI(2)*C1(2)+EJI(3)*C1(3) SINT=DOT/SINI IF(DABS(SINT).GT.0.00005D0)WRITE(*,1020) IF(DABS(SINT).GT.0.99995D0)GO TO 80 COST=DSQRT(1.D0-SINT*SINT) TANT=SINT/COST II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) LL=3*(L-1) COSSIN=COST*SINI DO 50 M=1,3 T=C1(M)/COSSIN SMI=(T-TANT*EJI(M))/DJI IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=DX*SMI SMK=T*(COSI*COSK-COSL)/(SINSIN*DJK) IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=DX*SMK SML=T*(COSI*COSL-COSK)/(SINSIN*DJL) IF(DABS(SML).GE.1.D-8)B(NOB,LL+M)=DX*SML SUM=SMI+SMK+SML 50 IF(DABS(SUM).GE.1.D-8)B(NOB,JJ+M)=-DX*SUM RETURN 60 IER=1 *C*C ADDED BY KATH, NEXT TWO LINES, SPRING 1991, LOOKING FOR ERRORS * WRITE(*,200) 200 FORMAT('OPLA ERROR') RETURN 70 IER=-1 WRITE(*,1000) RETURN 80 IER=-1 WRITE(*,1010) RETURN C 1000 FORMAT(' K-J-L IS COLINEAR (NO PLANE DEFINED FOR WAG OF I)') 1010 FORMAT(' I IS PERPENDICULAR TO J-K-L PLANE - USE VALENCE ANGLE 1 BENDS') 1020 FORMAT('+',86X,'*** WARNING: WAG OF A NON-PLANAR SYSTEM ***') C END SUBROUTINE TORS C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR TORSION AS DEFIN C BY R L HILDERBRANDT IN J MOLEC SPEC, 44, 599 (1972). C SUBROUTINE CODED BY M PETERSON, DEPT OF CHEMISTRY, UNIV OF TORONTO. C J AND K DEFINE THE BOND UNDER TORSION. C NI AND NL ARE THE NUMBER OF ATOMS ATTACHED TO J AND K RESPECTIVELY C (NI, NL <= 5). 2 DATA CARDS ARE READ: (1) CONTAINS NI ATOM NUMBERS C FOR THE I-TYPE ATOMS, AND (2) CONTAINS NL ATOM NUMBERS FOR THE L-TY C ATOMS (BOTH CARDS ARE IN 5I4 FORMAT). C C IATOM, LATOM: ATOM NUMBERS FOR THE I- AND L-TYPE ATOMS (SIZE: 5) C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,NI,J,K,NL,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RIJ(3),RJK(3),RLK(3),EIJ(3),EJK(3),ELK(3) DIMENSION CR(3),IATOM(5),LATOM(5),SJ(3),SK(3) C READ(15,*)(IATOM(I),I=1,NI) READ(15,*)(LATOM(L),L=1,NL) IF(NI.LE.0.OR.NI.GT.5)GO TO 110 IF(J.LE.0.OR.J.GT.NOAT)GO TO 110 IF(K.LE.0.OR.K.GT.NOAT)GO TO 110 IF(NL.LE.0.OR.NL.GT.5)GO TO 110 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 110 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 110 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJKSQ=0.D0 DXSQ=0.D0 DO 20 M=1,3 SJ(M)=0.D0 SK(M)=0.D0 T=X(M,K)-X(M,J) RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) 20 DXSQ=DXSQ+T*T DJK=1.D0/DSQRT(DJKSQ) DX=DSQRT(DXSQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DO 30 M=1,3 30 EJK(M)=RJK(M)*DJK JJ=3*(J-1) KK=3*(K-1) C C LOOP OVER THE I-TYPE ATOMS C DO 60 N=1,NI I=IATOM(N) IF(I.LE.0.OR.I.GT.NOAT)GO TO 110 DIJSQ=0.D0 DO 40 M=1,3 T=X(M,J)-X(M,I) RIJ(M)=T 40 DIJSQ=DIJSQ+T*T DIJ=1.D0/DSQRT(DIJSQ) COSJ=0.D0 DO 50 M=1,3 T=RIJ(M)*DIJ EIJ(M)=T 50 COSJ=COSJ-T*EJK(M) IF(DABS(COSJ).GT.0.99995D0)GO TO 120 SIN2J=(1.D0-COSJ*COSJ)*DFLOAT(NI) II=3*(I-1) CR(1)=EIJ(2)*EJK(3)-EIJ(3)*EJK(2) CR(2)=EIJ(3)*EJK(1)-EIJ(1)*EJK(3) CR(3)=EIJ(1)*EJK(2)-EIJ(2)*EJK(1) DO 60 M=1,3 T=CR(M)/SIN2J SMI=T*DIJ IF(DABS(SMI).GE.1.D-8)B(NOB,II+M)=-DX*SMI SMK=T*COSJ*DJK SK(M)=SK(M)+SMK SMJ=SMI-SMK 60 SJ(M)=SJ(M)+SMJ C C LOOP OVER THE L-TYPE ATOMS C DO 90 N=1,NL L=LATOM(N) IF(L.LE.0.OR.L.GT.NOAT)GO TO 110 DLKSQ=0.D0 DO 70 M=1,3 T=X(M,K)-X(M,L) RLK(M)=T 70 DLKSQ=DLKSQ+T*T DLK=1.D0/DSQRT(DLKSQ) COSK=0.D0 DO 80 M=1,3 T=RLK(M)*DLK ELK(M)=T 80 COSK=COSK+EJK(M)*T IF(DABS(COSK).GT.0.99995D0)GO TO 120 SIN2K=(1.D0-COSK*COSK)*DFLOAT(NL) LL=3*(L-1) CR(1)=ELK(3)*EJK(2)-ELK(2)*EJK(3) CR(2)=ELK(1)*EJK(3)-ELK(3)*EJK(1) CR(3)=ELK(2)*EJK(1)-ELK(1)*EJK(2) DO 90 M=1,3 T=CR(M)/SIN2K SML=T*DLK IF(DABS(SML).GE.1.D-8)B(NOB,LL+M)=-DX*SML SMJ=T*COSK*DJK SJ(M)=SJ(M)+SMJ SMK=SML-SMJ 90 SK(M)=SK(M)+SMK DO 100 M=1,3 SMJ=SJ(M) IF(DABS(SMJ).GE.1.D-8)B(NOB,JJ+M)=SMJ*DX SMK=SK(M) 100 IF(DABS(SMK).GE.1.D-8)B(NOB,KK+M)=SMK*DX RETURN 110 IER=1 WRITE(*,200) I,L 200 FORMAT('TORS ERROR AT I=',I2,' AND L=',I2) RETURN 120 IER=-1 WRITE(*,1030) RETURN 1030 FORMAT(' I-J-K OR J-K-L IS COLINEAR (NO TORSION POSSIBLE)') END SUBROUTINE LIBE C C THIS SUBROUTINE COMPUTES THE B MATRIX ELEMENTS FOR A LINEAR BEND C OR FOR A PAIR OF PERPENDICULAR LINEAR BENDS. C I AND K ARE THE END ATOMS. C J IS THE CENTRAL ATOM. C C A GIVES THE CARTESIAN COORDINATES OF A POINT IN SPACE, SUCH C THAT THE VECTOR FROM ATOM J TO POINT A IS PERPENDICULAR TO C THE LINE I-J-K AND SERVES TO ORIENT THE COORDINATES IN SPACE. C IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,X(3,NAT0),B(MX,MX) 1,NOB0,IS,UM(MX,MX) COMMON/SCHACH/RJI(3),RJK(3),EJK(3),UP(3),UN(3),UNIT(3) DIMENSION A(3) C READ(15,1000)A IF(I.LE.0.OR.I.GT.NOAT)GO TO 60 IF(J.LE.0.OR.J.GT.NOAT)GO TO 60 IF(K.LE.0.OR.K.GT.NOAT)GO TO 60 IF(L.NE.0)GO TO 60 IF(IX.LT.0.OR.IX.GT.NOAT)GO TO 60 IF(JX.LT.0.OR.JX.GT.NOAT)GO TO 60 IF(IX.NE.0.AND.JX.NE.0)GO TO 10 IX=1 JX=1 10 DJISQ=0.D0 DJKSQ=0.D0 DXSQ=0.D0 DJASQ=0.D0 DO 20 M=1,3 TP=X(M,J) T=X(M,I)-TP RJI(M)=T DJISQ=DJISQ+T*T T=X(M,K)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=X(M,JX)-X(M,IX) DXSQ=DXSQ+T*T T=A(M)-TP UN(M)=T 20 DJASQ=DJASQ+T*T DJI=DSQRT(DJISQ) DJK=DSQRT(DJKSQ) DX=DSQRT(DXSQ) DJA=DSQRT(DJASQ) IF(DABS(DX).LT.1.0D-10)DX=1.0D0 DOTJ=0.D0 DOTP=0.D0 DO 30 M=1,3 T=RJI(M)/DJI TP=RJK(M)/DJK EJK(M)=TP DOTJ=DOTJ+T*TP TP=UN(M)/DJA UNIT(M)=TP 30 DOTP=DOTP+T*TP TEST=DABS(DOTJ)-1.D0 IF(DABS(TEST).GT.0.00005D0)GO TO 70 IF(DABS(DOTP).GT.0.00005D0)GO TO 80 II=3*(I-1) JJ=3*(J-1) KK=3*(K-1) DO 40 M=1,3 T=UNIT(M) IF(DABS(T).LT.1.D-8)GO TO 40 T=-DX*T SMI=T/DJI B(NOB,II+M)=SMI SMK=T/DJK B(NOB,KK+M)=SMK B(NOB,JJ+M)=-SMI-SMK 40 CONTINUE IF(IC.EQ.6)RETURN NOB=NOB0+IS UP(1)=EJK(2)*UNIT(3)-EJK(3)*UNIT(2) UP(2)=EJK(3)*UNIT(1)-EJK(1)*UNIT(3) UP(3)=EJK(1)*UNIT(2)-EJK(2)*UNIT(1) DO 50 M=1,3 T=UP(M) IF(DABS(T).LT.1.D-8)GO TO 50 T=-DX*T SMI=T/DJI B(NOB,II+M)=SMI SMK=T/DJK B(NOB,KK+M)=SMK B(NOB,JJ+M)=-SMI-SMK 50 CONTINUE RETURN 60 IER=1 RETURN 70 IER=-1 WRITE(*,1020) RETURN 80 IER=-1 WRITE(*,1030) RETURN 1000 FORMAT(3D23.16) 1020 FORMAT(' I-J-K NOT COLINEAR - USE VALENCE ANGLE BEND') 1030 FORMAT(' ATOM A NOT PERPENDICULAR TO I-J-K AT J') END SUBROUTINE CTI(d1,d2,COO) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) PARAMETER (NAT0=1000,MX=3*NAT0,NX=MX) DIMENSION FINT(MX,MX),B(MX,NX),FCAR(MX,MX),nb(MX),lb(MX,NX) allocatable::am(:,:),ami(:,:),cm(:),t(:,:) character*7 COO(MX) integer*4 d1(MX),d2(MX) logical lex inquire(file='FILE.FC',exist=lex) if (.not.lex) 1call report(' FILE.FC not available, program stopped.') WRITE(*,*)' INTRISIC FORCE FIELD GENERATION' WRITE(*,*) WRITE(*,*)' FILE.FC ... CARTESIAN FORCE FIELD' WRITE(*,*) open(9,file='B.MAT') read(9,*)nob0 read(9,*)na read(9,*) do 13 i=1,nob0 read(9,*) read(9,*)nb(i) if(nb(i).gt.NX)call report('too many B-matrix elements') do 13 ii=1,nb(i) read(9,*)ia,ix,b(i,ii) 13 lb(i,ii)=ix+3*(ia-1) close(9) write(6,*)'B.MAT read, ',na,' atomic ',nob0,' internal' NAT3=NA OPEN(20,FILE='FILE.FC',STATUS='OLD') N=NAT3 NTT= -4 NTE= 5 IF (NTE.GT.N) GO TO 120 110 NTT=NTT+5 DO 100 I=NTT,N NT=I IF (NT.GT.NTE) NT=NTE 100 READ(20,27) (FCAR(I,J),J=NTT,NT) 27 FORMAT(4X,5D14.6) NTE=NTE+5 IF (NTE.LE.N) GOTO 110 120 NTT=NTT+5 DO 130 I=NTT,N 130 READ(20,27) (FCAR(I,J),J=NTT,I) CLOSE(20) DO 31 I=1,N DO 31 J=I+1,N 31 FCAR(I,J)=FCAR(J,I) DO 6 I=1,N DO 6 J=1,N C 6 FCAR(I,J)=FCAR(I,J)* 8.23D0 / 0.527D0 6 FCAR(I,J)=FCAR(I,J)*4.3597482D0/(0.529177249D0*0.529177249D0) c 1Hartree=4.3597482e-18 J; 1Bohr=0.529177249A write(6,*)'FILE.FC read' TOL= 1.0D-33 allocate(am(nob0,nob0),ami(nob0,nob0),t(nob0,nob0)) write(6,*)'fields allocated' alpha=1.0d-6 c c do 51 i=1,nob0 do 51 j=1,nob0 sum=0.0d0 do 511 lc=1,nb(i) do 511 kc=1,nb(j) 511 if(lb(i,lc).eq.lb(j,kc))sum=sum+b(i,lc)*b(j,kc) 51 t(i,j)=sum write(6,*)'t made' do 5 i=1,nob0 do 5 j=1,nob0 5 am(i,j)=t(i,j)**2 do 52 i=1,nob0 52 am(i,i)=am(i,i)+alpha write(6,*)'a done' CALL INV(am,ami,nob0,TOL) deallocate(am) allocate(cm(nob0)) do 7 i=1,nob0 cm(i)=0.0d0 do 7 lc=1,nb(i) do 7 ic=1,nb(i) 7 cm(i)=cm(i)+b(i,ic)*b(i,lc)*FCAR(lb(i,ic),lb(i,lc)) write(6,*)'c done' do 4 i=1,nob0 do 4 j=1,nob0 FINT(i,j)=0 if(i.eq.j)then do 41 ip=1,nob0 41 FINT(i,j)=FINT(i,j)+ami(i,ip)*cm(ip) endif 4 continue deallocate(cm,t) deallocate(ami) OPEN(2,FILE='INTY.FC') WRITE(2,*)' Force constants in internal coordinates:' CALL WRITEFF(MX,nob0,FINT) close(2) WRITE(*,*)'INTY.FC : OUTPUT WITH THE INTRISIC FORCE CONSTANTS' OPEN(2,FILE='IFC.TXT') do 8 i=1,nob0 8 write(2,800)i,COO(i),d1(i),d2(i),FINT(i,i) 800 format(i6,' ',a7,i6,i6,f12.6) close(2) RETURN END SUBROUTINE INV(MA,IA,N,TOL) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) REAL*8 MA,IA INTEGER*4 oo,n2 DIMENSION IA(N,N),MA(N,N) real*8, allocatable ::e(:) write(6,*)N*2*N/8,' bytes',N,2*N allocate(e(N*2*N)) n2=2*N write(6,*)'Inversion started' DO 1 ii=1,N DO 1 jj=1,N e(n2*(ii-1)+jj)=ma(ii,jj) e(n2*(ii-1)+jj+n)=0.0D0 1 if (ii.EQ.jj)e(n2*(ii-1)+jj+n)=1.0D0 write(6,*)'E initialized' DO 2 ii=1,n-1 IW=II if (DABS(e(n2*(iw-1)+iw)).LT.TOL) then DO 3 io=iw+1,n oo=IO if (io.GT.n)oo=io-n if (DABS(e(n2*(oo-1)+iw)).GE.TOL) goto 11 3 CONTINUE do 13 i=1,N 13 write(6,90)(MA(i,j),j=1,N) 90 format(12g12.4) WRITE(*,*)' INVERSE MATRIX CAN NOT BE FOUND ' STOP 11 CONTINUE DO 4 kk=1, 2*n w=e(n2*(iw-1)+kk) e(n2*(iw-1)+kk)=e(n2*(oo-1)+kk) 4 e(n2*(oo-1)+kk)=w ENDIF DO 5 jj=ii+1,n af=e(n2*(jj-1)+ii)/e(n2*(ii-1)+ii) DO 6 kk=ii+1, 2*n 6 e(n2*(jj-1)+kk)=e(n2*(jj-1)+kk)-e(n2*(ii-1)+kk)*af 5 e(n2*(jj-1)+ii)=0.0D0 2 CONTINUE write(6,*)'Inversion first loop' DO 7 ii=1,n-1 i2=n-ii+1 DO 8 jj=1,i2-1 j2=i2-jj af=e(n2*(j2-1)+i2)/e(n2*(i2-1)+i2) DO 9 kk=1, n 9 e(n2*(j2-1)+kk+n)=e(n2*(j2-1)+kk+n)-e(n2*(i2-1)+kk+n)*af e(n2*(j2-1)+i2)=0.0D0 8 CONTINUE 7 CONTINUE write(6,*)'Inversion second loop' DO 10 ii=1,n DO 12 jj=1,n 12 iA(ii,jj)=e(n2*(ii-1)+jj+n)/e(n2*(ii-1)+ii) 10 CONTINUE deallocate(E) RETURN END SUBROUTINE WRITEFF(MX,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(MX,MX) N1=1 1 N3=N1+4 IF(N3.GT.N)N3=N write(2,18)(ij,ij=N1,N3) 18 FORMAT(4x,5I14) DO 130 LN=N1,N 130 WRITE(2,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3)) N1=N1+5 IF(N3.LT.N)GOTO 1 17 FORMAT(I4,5D14.6) RETURN END subroutine report(s) character*(*) s write(6,*)s stop end SUBROUTINE HELIX(IZ,hl) c helical coordinates implicit none integer*4 NAT0,MX PARAMETER (NAT0=1000,MX=3*NAT0) integer*4 IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,j,ja, 1NOB0,IS,I,ix,jx,kx,IERR,ia,KK,IZ(*),hl(*) real*8 X,B,UM,A(3,3),c(3),U(3,3),EV(3),amf,XX,YY,ZZ,Xm,Ym,Zm, 1dwi,R2,R2j,wi,wj,Xj,Yj,Zj,YYj,ZZj COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,NAT0),B(MX,MX), 1NOB0,IS,UM(MX,MX) c read helix atoms only for the first helical coordinate: if(IC.EQ.14)READ(15,*)(hl(I),I=1,N1) c center of mass call ccm(c,N1,hl,IZ,X,NAT0) c moment of inertia: call cmi(A,N1,hl,c,IZ,X,NAT0) CALL TRED12(3,A,U,EV,2,IERR) C DIAGONAL TD(I,J) = U(IS,I)*A(IS,JS)*U(JS,j) write(6,6001)EV 6001 format(' Inertia moments: ',3f12.2) if(IC.eq.14.or.IC.EQ.15.or.IC.eq.16)then c translations of mass center jx=IC-13 do 7 ia=1,N1 do 7 ix=1,3 KK=ix+3*(hl(ia)-1) 7 B(NOB,KK)=B(NOB,KK)+U(ix,jx)*amf(IZ(hl(ia))) endif if(IC.eq.17.or.IC.eq.18.or.IC.eq.19)then c rotation around jx ix=IC-16 jx=ix+1 if(jx.gt.3)jx=1 kx=jx+1 if(kx.gt.3)kx=1 do 8 ia=1,N1 c R - radius in inertia moment coordinate system: i=hl(ia) Xm=(X(1,i)-c(1))*amf(IZ(i)) Ym=(X(2,i)-c(2))*amf(IZ(i)) Zm=(X(3,i)-c(3))*amf(IZ(i)) YY=U(1,jx)*Xm+ U(2,jx)*Ym+U(3,jx)*Zm ZZ=U(1,kx)*Xm+ U(2,kx)*Ym+U(3,kx)*Zm R2=YY**2+ZZ**2 if(R2.gt.1.0d-3)then do 81 ix=1,3 KK=ix+3*(i-1) 81 B(NOB,KK)=(YY*U(ix,jx)-ZZ*U(ix,kx))/R2 endif 8 continue endif if(IC.eq.20)then c moment derivative along the first=longest axis: c C = sum(i) xi xi mi= sum(i) Uax ai Ubx bi mi c dC/dgj = 2 Uax aj Ugx mj do 4 ia=1,N1 do 4 ix=1,3 KK=ix+3*(hl(ia)-1) B(NOB,KK)=0.0d0 do 4 jx=1,3 XX=X(jx,hl(ia))-c(jx) 4 B(NOB,KK)=B(NOB,KK)+U(ix,1)*U(jx,1)*XX*amf(IZ(hl(ia))) endif if(IC.eq.21)then c moment derivative along the other directions," breathing" do 5 ia=1,N1 do 5 ix=1,3 KK=ix+3*(hl(ia)-1) B(NOB,KK)=0.0d0 do 5 jx=1,3 Xm=(X(jx,hl(ia))-c(jx))*amf(IZ(hl(ia))) 5 B(NOB,KK)=B(NOB,KK)+(U(ix,2)*U(jx,2)+U(ix,3)*U(jx,3))*Xm endif if(IC.eq.22)then c moment derivative along the other directions," def/bend" do 6 ia=1,N1 do 6 ix=1,3 KK=ix+3*(hl(ia)-1) B(NOB,KK)=0.0d0 do 6 jx=1,3 Xm=(X(jx,hl(ia))-c(jx))*amf(IZ(hl(ia))) 6 B(NOB,KK)=B(NOB,KK)+(U(ix,2)*U(jx,2)-U(ix,3)*U(jx,3))*Xm endif if(IC.eq.23)then c rotation difference around X, "unwinding" do 9 ia=1,N1 c R - radius in inertia moment coordinate system: i=hl(ia) Xm=(X(1,i)-c(1))*amf(IZ(i)) Ym=(X(2,i)-c(2))*amf(IZ(i)) Zm=(X(3,i)-c(3))*amf(IZ(i)) YY=U(1,2)*Xm+ U(2,2)*Ym+U(3,2)*Zm ZZ=U(1,3)*Xm+ U(2,3)*Ym+U(3,3)*Zm R2=YY**2+ZZ**2 if(R2.gt.1.0d-3)then wi=atan2(YY,ZZ) do 91 ix=1,3 KK=ix+3*(i-1) dwi=(YY*U(ix,2)-ZZ*U(ix,3))/R2 B(NOB,KK)=0.0d0 do 91 ja=1,N1 if(ja.ne.ia)then j=hl(ja) Xj=(X(1,j)-c(1))*amf(IZ(j)) Yj=(X(2,j)-c(2))*amf(IZ(j)) Zj=(X(3,j)-c(3))*amf(IZ(j)) YYj=U(1,2)*Xj+ U(2,2)*Yj+U(3,2)*Zj ZZj=U(1,3)*Xj+ U(2,3)*Yj+U(3,3)*Zj R2j=YYj**2+ZZj**2 if(R2j.gt.1.0d-3)then wj=atan2(Yj,Zj) B(NOB,KK)=B(NOB,KK)+(wi-wj)*dwi endif endif 91 continue endif 9 continue endif return end SUBROUTINE MONO(IZ,ICODE) c coordinates defining 6 molecular coordinates c c (list1) c atoms 1....N1 c (molecule 1) c implicit none integer*4 IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,NOB0,IS,I,ICODE, 1IZ(*),ix,ia,NAT0,iy,zz,MX PARAMETER (NAT0=1000,MX=3*NAT0) real*8 X,B,UM,m1,amf,c1(3),z,y COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,NAT0),B(MX,MX), 1NOB0,IS,UM(MX,MX) integer*4,allocatable:: list1(:) allocate(list1(N1)) READ(15,*)(list1(I),I=1,N1) c total mass: m1=0.0d0 do 12 i=1,N1 12 m1=m1+amf(IZ(list1(i))) c c centers of mass: call ccm(c1,n1,list1,IZ,X,NAT0) do 1 i=1,N1 ia=list1(i) c translation: if(ICODE.ge.24.and.ICODE.le.26)then ix=ICODE-23 B(NOB,ix+3*(ia-1))=amf(IZ(ia))/m1 endif c rotation: if(ICODE.ge.27.and.ICODE.le.29)then ix=ICODE-26 iy=ix+1 if(iy.gt.3)iy=1 zz=iy+1 if(zz.gt.3)zz=1 y=X(iy,ia)-c1(iy) z=X(zz,ia)-c1(zz) B(NOB,iy+3*(ia-1))=B(NOB,iy+3*(ia-1))+z/(y**2+z**2) B(NOB,zz+3*(ia-1))=B(NOB,zz+3*(ia-1))-y/(y**2+z**2) endif 1 continue return END SUBROUTINE INTER(IZ,ICODE) c coordinates defining 6 inter-molecular coordinates c c (list1) (list2) c atoms 1....N1 atoms 1..N2 c (molecule 1) (molecule 2) c implicit none integer*4 MX,IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,NOB0,IS,I,ICODE, 1IZ(*),ix,ia,IDIF,IERR,II,is1,is2,KA,KK,KX,NAT0 PARAMETER (NAT0=1000,MX=3*NAT0) real*8 X,B,UM,th1(3,3),th2(3,3),c1(3),c2(3),m1,m2,c0(3), 1u1(3,3),u2(3,3),d1(3),d2(3),PI,STEP,FF,acc,amf, 1v1(3),v2(3),e1(3),e2(3),D,T,A1,A2,B1,B2,D0,T0,A10,A20,B10,B20, 1u01(3,3),u02(3,3),ph COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,X(3,NAT0),B(MX,MX), 1NOB0,IS,UM(MX,MX) integer*4,allocatable:: list1(:),list2(:) D0=0.0d0 A10=0.0d0 A20=0.0d0 B10=0.0d0 B20=0.0d0 T0=0.0d0 PI=4.0d0*atan(1.0d0) allocate(list1(N1),list2(N2)) C READ(15,*)(list1(I),I=1,N1) READ(15,*)(list2(I),I=1,N2) m1=0.0d0 do 12 i=1,N1 12 m1=m1+amf(IZ(list1(i))) m2=0.0d0 do 13 i=1,N2 13 m2=m2+amf(IZ(list2(i))) c c numerical derivatives STEP=0.01D0 do 100 IDIF=0,6*NOAT if(IDIF.gt.0)THEN if(IDIF.GT.3*NOAT)THEN KK=IDIF-3*NOAT FF=-1.0D0 ELSE KK=IDIF FF=1.0D0 ENDIF KA=(KK+2)/3 KX=KK-3*(KA-1) X(KX,KA)=X(KX,KA)+FF*STEP ELSE FF=0.0d0 KA=0 KX=0 ENDIF c calculate centers of mass: call ccm(c1,n1,list1,IZ,X,NAT0) call ccm(c2,n2,list2,IZ,X,NAT0) d=0.0d0 do 6 ix=1,3 c0(ix)=c1(ix)-c2(ix) 6 d=d+c0(ix)**2 d=dsqrt(d) do 7 ix=1,3 7 c0(ix)=c0(ix)/d c c inter-molecular stretching - analytically: if(IDIF.EQ.0)then do 8 ia=1,NOAT c is atom ia part of the first or second molecule?: is1=0 is2=0 do 9 i=1,N1 9 if(list1(i).eq.ia)is1=1 do 10 i=1,N2 10 if(list2(i).eq.ia)is2=1 if(is1.ne.0.or.is2.ne.0)then do 11 ix=1,3 ii=ix+3*(ia-1) if(ICODE.EQ.7.or.ICODE.EQ.8)then if(is1.eq.1)B(NOB,ii)= c0(ix)*amf(IZ(ia))/m1 if(is2.eq.1)B(NOB,ii)=-c0(ix)*amf(IZ(ia))/m2 endif 11 continue endif 8 continue endif c calculate moments of inertia: call cmi(th1,N1,list1,c1,IZ,X,NAT0) call cmi(th2,N2,list2,c2,IZ,X,NAT0) c diagonalize moments of inertia: CALL TRED12(3,th1,u1,d1,2,IERR) CALL TRED12(3,th2,u2,d2,2,IERR) if(IDIF.EQ.0)then c record for equilibrium: u01=u1 u02=u2 else c adjust phase to equilibrium do 103 i=1,3 ph=u1(1,i)*u01(1,i)+u1(2,i)*u01(2,i)+u1(3,i)*u01(3,i) if(ph.lt.0.0d0)then u1(1,i)=-u1(1,i) u1(2,i)=-u1(2,i) u1(3,i)=-u1(3,i) endif ph=u2(1,i)*u02(1,i)+u2(2,i)*u02(2,i)+u2(3,i)*u02(3,i) if(ph.lt.0.0d0)then u2(1,i)=-u2(1,i) u2(2,i)=-u2(2,i) u2(3,i)=-u2(3,i) endif 103 continue endif if(IDIF.EQ.0)write(6,607)IDIF,d1,d2 607 format(' Inertias:',i2,6F10.3) c angles between longest moment axes and the distance vecto A1=acc(c0(1)*u1(1,1)+c0(2)*u1(2,1)+c0(3)*u1(3,1)) A2=acc(c0(1)*u2(1,1)+c0(2)*u2(2,1)+c0(3)*u2(3,1)) c c arbitrary vectors defining rotating axes v1(1)=u1(2,1)*c0(3)-u1(3,1)*c0(2) v1(2)=u1(3,1)*c0(1)-u1(1,1)*c0(3) v1(3)=u1(1,1)*c0(2)-u1(2,1)*c0(1) v2(1)=u2(2,1)*c0(3)-u2(3,1)*c0(2) v2(2)=u2(3,1)*c0(1)-u2(1,1)*c0(3) v2(3)=u2(1,1)*c0(2)-u2(2,1)*c0(1) e1(1)=v1(2)*u1(3,1)-v1(3)*u1(2,1) e1(2)=v1(3)*u1(1,1)-v1(1)*u1(3,1) e1(3)=v1(1)*u1(2,1)-v1(2)*u1(1,1) e2(1)=v2(2)*u2(3,1)-v2(3)*u2(2,1) e2(2)=v2(3)*u2(1,1)-v2(1)*u2(3,1) e2(3)=v2(1)*u2(2,1)-v2(2)*u2(1,1) c angles between the second longest u and the C0U plane B1=acc(e1(1)*u1(1,2)+e1(2)*u1(2,2)+e1(3)*u1(3,2)) B2=acc(e2(1)*u2(1,2)+e2(2)*u2(2,2)+e2(3)*u2(3,2)) c torsion angle T=acc(v1(1)*v2(2)+v1(2)*v2(2)+v1(3)*v2(3)) if(v1(1)*v2(1)+v1(2)*v2(2)+v1(3)*v2(3).lt.0.0d0)T=-T c try to keep conventions: c torsion -180...180 if(T.lt.-180.0d0)T= 360.0d0+T if(T.gt.+180.0d0)T=-360.0d0+T c axis angles: the sharp ones if(A1.gt.90.0d0)A1=180.0d0-A1 if(A2.gt.90.0d0)A2=180.0d0-A2 if(B1.gt.90.0d0)B1=180.0d0-B1 if(B2.gt.90.0d0)B2=180.0d0-B2 if(IDIF.eq.0)then D0=d A10=A1 A20=A2 B10=B1 B20=B2 T0=T else c for the rare occasion +/-180, watch angles: call rare(T,T0,' t') call rare(A1,A10,'a1') call rare(A2,A20,'a2') call rare(B1,B10,'b1') call rare(B2,B20,'b2') c B(NOB ,KK)=B(NOB ,KK)+FF*d if(ICODE.EQ.7)then B(NOB+1,KK)=B(NOB+1,KK)+FF*A1 B(NOB+2,KK)=B(NOB+2,KK)+FF*A2 B(NOB+3,KK)=B(NOB+3,KK)+FF*B1 B(NOB+4,KK)=B(NOB+4,KK)+FF*B2 B(NOB+5,KK)=B(NOB+5,KK)+FF*T else if(ICODE.EQ. 9)B(NOB,KK)=B(NOB,KK)+FF*A1 if(ICODE.EQ.10)B(NOB,KK)=B(NOB,KK)+FF*A2 if(ICODE.EQ.11)B(NOB,KK)=B(NOB,KK)+FF*B1 if(ICODE.EQ.12)B(NOB,KK)=B(NOB,KK)+FF*B2 if(ICODE.EQ.13)B(NOB,KK)=B(NOB,KK)+FF*T endif ENDIF if(IDIF.NE.0)X(KX,KA)=X(KX,KA)-FF*STEP 100 continue if(ICODE.EQ.7)then do 101 i=1,5 do 101 KK=1,3*NOAT 101 B(NOB+i,KK)=B(NOB+i,KK)/2.0d0/STEP*PI/180.0d0 NOB=NOB+5 write(6,609)D0,A10,A20,B10,B20,T0,NOB-6,NOB 609 format(' D A1 A2 B1 B2 T:',F10.4,5F8.2,i4,' -',i4) else if(ICODE.GT.8)then do 102 KK=1,3*NOAT 102 B(NOB,KK)=B(NOB,KK)/2.0d0/STEP*PI/180.0d0 endif write(6,609)D0,A10,A20,B10,B20,T0,NOB endif END subroutine rare(T,T0,s2) real*8 T,T0 character*2 s2 if(dabs(T+360.0d0-T0).lt.dabs(T-T0))then T=T+360.0d0 write(6,*)'Periodicity change '//s2,T endif if(dabs(T-360.0d0-T0).lt.dabs(T-T0))then T=T-360.0d0 write(6,*)'Periodicity change '//s2,T endif return end SUBROUTINE TRED12(N,A,Z,D,IEIG,IERR) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION A(N,N),Z(N,N),D(*) real*8,allocatable::E(:) allocate(E(N)) IF (IEIG .LT. 0) GO TO 110 DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C 110 ICODE = IABS(IEIG) IF (N .EQ. 1) GO TO 320 C :::::::::: FOR I=N STEP -1 UNTIL 2 DO -- :::::::::: DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C :::::::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED) :::::::::: DO 120 K = 1, L 120 SCALE = SCALE + DABS(Z(I,K)) C IF (DABS(SCALE) .GT. 1.0D-10) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L F = Z(I,K) / SCALE Z(I,K) = F H = H + F * F 150 CONTINUE C G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0D0 C DO 240 J = 1, L IF (ICODE .EQ. 2) Z(J,I) = Z(I,J) / (SCALE * H) G = 0.0D0 C :::::::::: FORM ELEMENT OF A*U :::::::::: DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C :::::::::: FORM ELEMENT OF P :::::::::: 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C :::::::::: FORM REDUCED A :::::::::: DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C DO 280 K = 1, L 280 Z(I,K) = SCALE * Z(I,K) C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0D0 E(1) = 0.0D0 C :::::::::: ACCUMULATION OF TRANSFORMATION MATRICES :::::::::: IF (ICODE .NE. 2) GO TO 600 DO 500 I = 1, N L = I - 1 IF (DABS(D(I)) .LT. 1.0D-10) GO TO 380 C DO 360 J = 1, L G = 0.0D0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0D0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 620 CALL TQL12 (N,Z,D,IERR,ICODE,E) RETURN C :::::::::: ALTERNATE FINAL LOOP FOR EIGENVALUES ONLY :::::::::: 600 DO 610 I=1,N 610 D(I) = Z(I,I) GO TO 620 END C SUBROUTINE TQL12(N,Z,D,IERR,ICODE,E) IMPLICIT REAL*8(A-H,O-Z) IMPLICIT INTEGER*4 (I-N) DIMENSION Z(N,N),D(*),E(*) REAL*8 MACHEP EPS = 1.0D0 10 EPS = 0.50D0*EPS TOL1 = EPS + 1.0D0 IF((TOL1.GT.1.0D0).AND.(TOL1-EPS.EQ.1.0D0)) GO TO 10 IF(TOL1-EPS.EQ.1.0D0) EPS = EPS + EPS MACHEP = EPS C C MACHEP=16.0D0**(-13) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C F = 0.0D0 B = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = MACHEP * (DABS(D(L)) + DABS(E(L))) IF (B .LT. H) B = H C :::::::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT :::::::::: DO 110 M = L, N IF (DABS(E(M)) .LE. B) GO TO 120 C :::::::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP :::::::::: 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C :::::::::: FORM SHIFT :::::::::: L1 = L + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) R = DSQRT(P*P+1.0D0) D(L) = E(L) / (P + DSIGN(R,P)) H = G - D(L) C DO 140 I = L1, N 140 D(I) = D(I) - H C F = F + H C :::::::::: QL TRANSFORMATION :::::::::: P = D(M) C = 1.0D0 S = 0.0D0 MML = M - L C :::::::::: FOR I=M-1 STEP -1 UNTIL L DO -- :::::::::: DO 200 II = 1, MML I = M - II G = C * E(I) H = C * P IF (DABS(P) .LT. DABS(E(I))) GO TO 150 C = E(I) / P R = DSQRT(C*C+1.0D0) E(I+1) = S * P * R S = C / R C = 1.0D0 / R GO TO 160 150 C = P / E(I) R = DSQRT(C*C+1.0D0) E(I+1) = S * E(I) * R S = 1.0D0 / R C = C * S 160 P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C :::::::::: FORM VECTOR :::::::::: IF (ICODE .NE. 2) GO TO 200 DO 180 K = 1, N H = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * H Z(K,I) = C * Z(K,I) - S * H 180 CONTINUE C 200 CONTINUE C E(L) = S * P D(L) = C * P IF (DABS(E(L)) .GT. B) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C :::::::::: ORDER EIGENVALUES AND EIGENVECTORS :::::::::: DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .LE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C IF (ICODE .NE. 2) GO TO 300 DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C :::::::::: SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS :::::::::: 1000 IERR = L 1001 RETURN END function acc(p) implicit none real*8 A,acc,p,PI PI=4.0d0*atan(1.0d0) A=0.0D0 IF(p.LE.-1.0D0)A=PI IF(DABS(p).LT.1.0D0)A=ACOS(p) acc=180.0D0*A/PI return end subroutine cmi(t,N,l,c,IZ,X,NAT0) implicit none integer*4 ix,jx,i,N,l(*),IZ(*),ia,NAT0 real*8 t(3,3),c(3),xi,xj,X(3,NAT0),amf do 1 ix=1,3 do 1 jx=1,3 t(ix,jx)=0.0d0 do 1 i=1,N ia=l(i) xi=X(ix,ia)-c(ix) xj=X(jx,ia)-c(jx) 1 t(ix,jx)=t(ix,jx)+xi*xj*amf(IZ(ia)) return end subroutine ccm(c,n,l,z,X,NAT0) implicit none integer*4 n,l(*),z(*),ix,i,NAT0 real*8 c(3),X(3,NAT0),amf,mm mm=0.0d0 do 3 i=1,n 3 mm=mm+amf(z(l(i))) do 1 ix=1,3 c(ix)=0.0d0 do 2 i=1,n 2 c(ix)=c(ix)+X(ix,l(i))*amf(z(l(i))) 1 c(ix)=c(ix)/mm return end function amf(z) implicit none real*8 amf integer*4 z,MENDELEV PARAMETER (MENDELEV=89) real*4 AM(MENDELEV) data AM/1.008,4.003, 2 6.941, 9.012, 10.810,12.011,14.007,15.999,18.998,20.179, 3 22.990,24.305, 26.981,28.086,30.974,32.060,35.453,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ if(z.lt.1.or.z.gt.MENDELEV)then write(6,*)z call report('invalid atomic number') endif amf=dble(AM(z)) return end