PROGRAM BMATRIX IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4(I-N) real*8 u(3,3) COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,IS allocatable::X(:,:),B(:,:),UM(:,:),B0(:,:),LV(:),IZ(:) character*3 s3 logical lp WRITE(6,500) 500 format(' (SYMMETRIZED) B-MATRIX GENERATION',/, 1 ' Input: FILE.UMA ... coordinate definition',/, 1 ' FILE.X ... geometry',/, 1 ' FILE.FC ... Cartesian force field',/, 1 ' Output: B.MAT ... B-matrix',/, 1 ' INTY.FC ... internal force field',/, 1 ' INTY.SCR ... INTY.FC binary') icart=0 KEY=-1 lp=.false. if(iargc().gt.0)then call getarg(1,s3) read(s3,*)KEY endif IER=0 ISCAN=0 IS=0 tol=1.0d-8 OPEN(15,FILE='FILE.X') READ(15,*) READ(15,*)NOAT NA=NOAT*3 allocate(X(3,NOAT),B0(NA,NA),LV(NA),IZ(NA)) do 21 I=1,NOAT 21 read(15,*)IZ(I),(X(J,I),J=1,3) close(15) OPEN(15,FILE='FILE.UMA') READ(15,*) READ(15,*) READ(15,*)NOATt,NOB0 if(NOATt.ne.NOAT)call report('Different Nat in UMA and X') WRITE(6,6001)NOAT,NOB0 6001 format(I6,' atoms,',I6,' internal coordinates') allocate(B(NOB0,NA)) B=0.0D0 DO 20 I=1,NOAT 20 READ(15,*) NOB=0 277 NOB=NOB+1 READ(15,*)ICODE,N1,N2,N3,N4,N5,N6,FACTOR write(6,615)ICODE,N1,N2,N3,N4,N5,N6 615 format(7i6) if((N1.lt.0.or.N2.lt.0).and..not.lp)then lp=.true. write(6,*)'Crystal cell detected' open(9,file='CELL.TXT',status='old') read(9,*)(u(1,ix),ix=1,3) read(9,*)(u(2,ix),ix=1,3) read(9,*)(u(3,ix),ix=1,3) close(9) write(6,*)'CELL.TXT read' endif IF(DABS(FACTOR).LT.1.0D-6)FACTOR=1.0D0 IC=IABS(ICODE) IF(IC.LT.1.OR.IC.GT.8)call report('ILLEGAL CODE') GO TO (1,2,3,4,5,6,7,8),IC 1 CALL BOST(NOB0,NA,X,B,u) GO TO 60 2 CALL BEND(NOB0,NA,X,B,u) GO TO 60 3 CALL OPLA(NOB0,NA,X,B,u) GO TO 60 4 CALL TORS(NOB0,NA,X,B,u) GO TO 60 5 IS=IS+1 6 CALL LIBE(NOB0,NA,X,B) GO TO 60 7 CALL INTER(NOB0,IZ,NA,X,B) GO TO 60 8 CALL XYZ(NOB0,NOB,NA,X,B) icart=1 60 IF(IER.NE.0)GO TO 190 IF(DABS(FACTOR-1.0D0).LT.1.0D-8)GO TO 105 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 105 IF(ICODE.GT.0)GO TO 2775 2775 if(NOB.lt.NOB0)goto 277 READ(15,*)IUM IF(IUM.NE.0)THEN NUM=IUM if(IUM.EQ.1)NUM=NA-6 WRITE(6,600)NUM 600 format(' USER-DEFINED U MATRIX,',/, 1 i4,' non-reduntant coordinates expected') allocate(UM(NUM,NOB0)) UM=0.0D0 DO 2772 I=1,NUM 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 c symmetry adapted Si = U(i,j) Ij, Ij .. internal, j=1..NOB0 DO 1556 I=1,NUM DO 1556 J=1,NA B0(I,J)=0.0d0 DO 1556 K=1,NOB0 1556 B0(I,J)=B0(I,J)+UM(I,K)*B(K,J) c B0(NUM,NA)=UM(NUM,NOB0)*B(NOB0,NA) c Si = B0 *X, X .. Cartesian ELSE NUM=NOB0 B0=B ENDIF CLOSE(15) DO 155 I=1,NUM DO 155 J=1,NA 155 IF(DABS(B0(I,J)).LT.tol) B0(I,J) = 0.0D0 WRITE(*,*)' OUTPUT FILE: B.MAT' OPEN(9,FILE='B.MAT') WRITE(9,901)NUM,NA 901 format(I6,' Internal coordinates',/, 1I6,' Cartesian coordinates',/, 1' Transformation matrix, non-zero elements:') do 13 I=1,NUM WRITE(9,9011)I 9011 format(I6,' internal:') nn=0 do 14 J=1,NA 14 if(DABS(B0(I,J)).gt.tol)nn=nn+1 write(9,9012)nn 9012 format(I6) nn=0 do 13 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 13 continue CLOSE(9) write(6,*)KEY,NOB0,NA if(KEY.eq.-1.and.((NOB0.eq.NA-6.and.icart.eq.0). 1 or.(NOB0.eq.NA.and.icart.eq.1)))then WRITE(*,*) 1 ' COMBINE THE B-MATRIX AND CARTESIAN FORCE FIELD (1/0) ?' READ(*,*)KEY IF(KEY.EQ.1)CALL CTI endif WRITE(*,*) WRITE(*,*)' >> PROGRAM TERMINATED <<' 190 IF(IER.EQ.1)WRITE(*,1130) ISCAN=ISCAN+1 IER=0 1130 FORMAT(' ILLEGAL SPECIFICATION OF I, J, K, L, IX OR JX') END SUBROUTINE BOST(NOB0,NA,X,B,u) 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) dimension X(3,NA/3),B(NOB0,NA),u(3,3),xi(3),xj(3) COMMON/BMAT/IC,IR,JR,K,L,IX,JX,NOAT,NOB,IER,IS COMMON/SCHACH/RIJ(3),RJK(3),RJL(3),EIJ(3),EJK(3),EKL(3) C c periodic cell: call setx(IR,NA,xi,X,u,ia1) call setx(JR,NA,xj,X,u,ja1) DO 10 M=1,3 10 RIJ(M)=xj(M)-xi(M) DIJ=DSQRT(RIJ(1)**2+RIJ(2)**2+RIJ(3)**2) II=3*(ia1-1) JJ=3*(ja1-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 END SUBROUTINE BEND(NOB0,NA,X,B,u) 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) dimension X(3,NA/3),B(NOB0,NA),u(3,3),xi(3),xj(3),xk(3) COMMON/BMAT/IC,IR,JR,KR,L,IX,JX,NOAT,NOB,IER,IS COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EKL(3) C call setx(IR,NA,xi,X,u,ia1) call setx(JR,NA,xj,X,u,ja1) call setx(KR,NA,xk,X,u,ka1) 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=xj(M) T=xi(M)-TP RJI(M)=T DJISQ=DJISQ+T*T T=xk(M)-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*(ia1-1) JJ=3*(ja1-1) KK=3*(ka1-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 60 IER=-1 WRITE(*,1000) 1000 FORMAT(' I-J-K IS COLINEAR - USE LINEAR BEND') RETURN END SUBROUTINE OPLA(NOB0,NA,X,B,u) 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) dimension X(3,NA/3),B(NOB0,NA),u(3,3),xi(3),xj(3),xk(3),xl(3) COMMON/BMAT/IC,IR,JR,KR,LR,IX,JX,NOAT,NOB,IER,IS COMMON/SCHACH/RJI(3),RJK(3),RJL(3),EJI(3),EJK(3),EJL(3) DIMENSION C1(3) C call setx(IR,NA,xi,X,u,ia1) call setx(JR,NA,xj,X,u,ja1) call setx(KR,NA,xk,X,u,ka1) call setx(LR,NA,xl,X,u,la1) 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=xj(M) T=xi(M)-TP RJI(M)=T DJISQ=DJISQ+T*T T=xk(M)-TP RJK(M)=T DJKSQ=DJKSQ+T*T T=xl(M)-xj(M) 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*(ia1-1) JJ=3*(ja1-1) KK=3*(ka1-1) LL=3*(la1-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 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 INTER(NOB0,IZ,NA,X,B) 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 IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,NOB0,IS,I,NA, 1IZ(*),MENDELEV,ix,jx,ia,IDIF,IERR,II,is1,is2,KA,KK,KX PARAMETER (MENDELEV=89) real*8 X,B,th1(3,3),th2(3,3),c1(3),c2(3),m1,m2,xi,xj,c0(3), 1u1(3,3),u2(3,3),d1(3),d2(3),PI,STEP,FF,acc, 1v1(3),v2(3),e1(3),e2(3),D,T,A1,A2,B1,B2,D0,T0,A10,A20,B10,B20 dimension X(3,NA/3),B(NOB0,NA) COMMON/BMAT/IC,N1,N2,N3,N4,N5,N6,NOAT,NOB,IER,IS integer*4,allocatable:: list1(:),list2(:) 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/ 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+dble(AM(IZ(list1(i)))) m2=0.0d0 do 13 i=1,N2 13 m2=m2+dble(AM(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: do 3 ix=1,3 c1(ix)=0.0d0 do 4 i=1,N1 4 c1(ix)=c1(ix)+X(ix,list1(i))*dble(AM(IZ(list1(i)))) c1(ix)=c1(ix)/m1 c2(ix)=0.0d0 do 5 i=1,N2 5 c2(ix)=c2(ix)+X(ix,list2(i))*dble(AM(IZ(list2(i)))) 3 c2(ix)=c2(ix)/m2 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(is1.eq.1)B(NOB,ii)= c0(ix)*dble(AM(IZ(ia)))/m1 11 if(is2.eq.1)B(NOB,ii)=-c0(ix)*dble(AM(IZ(ia)))/m2 endif 8 continue endif c calculate moments of inertia: do 1 ix=1,3 do 1 jx=1,3 th1(ix,jx)=0.0d0 th2(ix,jx)=0.0d0 do 2 i=1,N1 ia=list1(i) xi=x(ix,ia)-c1(ix) xj=x(jx,ia)-c1(jx) 2 th1(ix,jx)=th1(ix,jx)+xi*xj*dble(AM(IZ(ia))) do 1 i=1,N2 ia=list2(i) xi=x(ix,ia)-c2(ix) xj=x(jx,ia)-c2(jx) 1 th2(ix,jx)=th2(ix,jx)+xi*xj*dble(AM(IZ(ia))) 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)write(6,607)d1,d2 607 format(' Inertias:',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 write(6,609)D,A1,A2,B1,B2,T if(IDIF.eq.0)then D0=d A10=A1 A20=A2 B10=B1 B20=B2 T0=1 else c B(NOB ,KK)=B(NOB ,KK)+FF*d 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 ENDIF if(IDIF.NE.0)X(KX,KA)=X(KX,KA)-FF*STEP 100 continue 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 609 format(' D A1 A2 B1 B2 T:',F10.4,5F8.2) END SUBROUTINE TORS(NOB0,NA,X,B,u) 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) dimension X(3,NA/3),B(NOB0,NA),u(3,3) COMMON/BMAT/IC,NI,JR,KR,NL,IX,JX,NOAT,NOB,IER,IS 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), 1xil(5,3),xll(5,3),xj(3),xk(3),iil(5),ill(5) C READ(15,*)(IATOM(I),I=1,NI) READ(15,*)(LATOM(L),L=1,NL) do 1 I=1,NI call setx(IATOM(I),NA,xj,X,u,ia1) iil(I)=ia1 do 1 ii=1,3 1 xil(I,ii)=xj(ii) do 2 I=1,NL call setx(LATOM(I),NA,xj,X,u,ia1) ill(I)=ia1 do 2 ii=1,3 2 xll(I,ii)=xj(ii) call setx(JR,NA,xj,X,u,ja1) call setx(KR,NA,xk,X,u,ka1) 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=xk(M)-xj(M) 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*(ja1-1) KK=3*(ka1-1) C C LOOP OVER THE I-TYPE ATOMS C DO 60 N=1,NI DIJSQ=0.D0 DO 40 M=1,3 T=xj(M)-xil(N,M) 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*(iil(N)-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 DLKSQ=0.D0 DO 70 M=1,3 T=xk(M)-xll(N,M) 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*(ill(N)-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(NOB0,NA,X,B) 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) dimension X(3,NA/3),B(NOB0,NA) COMMON/BMAT/IC,I,J,K,L,IX,JX,NOAT,NOB,IER,IS 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 IMPLICIT None INTEGER*4 NOB0,NA,I,J,K real*8 TOL real*8,allocatable::FINT(:,:),B(:,:),T(:,:),F(:,:),BBt(:,:), 1BBtI(:,:),Bti(:,:) WRITE(*,*)' INTRISIC FORCE FIELD GENERATION' NOB0=1 NA=1 allocate(B(NOB0,NA)) call rb(0,NOB0,NA,B) deallocate(B) allocate(FINT(NOB0,NOB0),B(NOB0,NA),F(NA,NA),BBt(NOB0,NOB0), 1BBtI(NOB0,NOB0),T(NOB0,NOB0),Bti(NOB0,NA)) call rb(1,NOB0,NA,B) call READFC(NA,F) BBt=0.0d0 do 1 I=1,NOB0 do 1 J=1,NOB0 do 1 K=1,NA 1 BBt(I,J)=BBt(I,J)+B(I,K)*B(J,K) TOL= 1.0D-13 call sm('BBt:',BBt,NOB0) CALL INV(BBt,BBtI,NOB0,TOL) call sm('BBti:',BBtI,NOB0) T=0.0d0 DO 8 I=1,NOB0 DO 8 J=1,NOB0 DO 8 K=1,NOB0 8 T(I,J)=T(I,J)+BBtI(I,K)*BBt(K,J) DO 42 I=1,NOB0 DO 42 J=1,NOB0 IF ((I.EQ.J).AND.(DABS(T(I,J)-1).GT.0.0001d0))THEN WRITE(*,*)I,J,T(I,J),'-NOT ONE' stop ENDIF IF ((I.NE.J).AND.(DABS(T(I,J)).GT.0.0001d0))THEN WRITE(*,*)I,J,T(I,J),'NOT ZERO' stop ENDIF 42 CONTINUE call sm('BBti x BBt:',T,NOB0) Bti=0.0d0 DO 4 I=1,NOB0 DO 4 J=1,NA DO 4 K=1,NOB0 4 Bti(I,J)=Bti(I,J)+BBtI(I,K)*B(K,J) C THE SYMBOLIC (BT)-1 MATRIX IS IN Bti NOW deallocate(T) allocate(T(NOB0,NA)) T=0.0d0 DO 5 I=1,NOB0 DO 5 J=1,NA DO 5 K=1,NA 5 T(I,J)=T(I,J)+Bti(I,K)*F(K,J) FINT=0.0d0 DO 7 I=1,NOB0 DO 7 J=1,NOB0 DO 7 K=1,NA 7 FINT(I,J)=FINT(I,J)+T(I,K)*Bti(J,K) c I = B X c Fc = Bt Fi B c B Fc = B Bt Fi B c [b Bt]^-1 B Fc [B Bt]-1 = Fi CALL WRITEFF(NOB0,FINT) OPEN(29,FILE='INTY.SCR',FORM='UNFORMATTED') WRITE(29) (NOB0*(NOB0+1))/2,((FINT(I,J),I=1,J),J=1,NOB0) CLOSE(29) write(*,*)'written' RETURN END SUBROUTINE INV(MA,IA,N,TOL) IMPLICIT none INTEGER*4 N,oo,ii,jj,IW,io,kk,i2,j2 real*8 MA(N,N),IA(N,N),TOL,w real*8,allocatable::E(:,:) allocate(E(N,2*N)) DO 1 ii=1,N DO 1 jj=1,N e(ii,jj)=ma(ii,jj) E(II,JJ+N)=0.0D0 1 if (ii.EQ.jj)e(ii,jj+N)=1.0D0 DO 2 ii=1,N-1 IW=II if (DABS(e(iw,iw)).LT.TOL) then DO 3 io=iw+1,N oo=IO if (io.GT.N) oo=io-N if (DABS(e(oo,iw)).GE.TOL) goto 11 3 CONTINUE WRITE(*,*)' INVERSE MATRIX CAN NOT BE FOUND ' STOP 11 CONTINUE DO 4 kk=1, 2*N w=e(iw,kk) e(iw,kk)=e(oo,kk) 4 e(oo,kk)=w ENDIF DO 5 jj=ii+1,N DO 6 kk=ii+1, 2*N 6 e(jj,kk)=e(jj,kk)-e(ii,kk)*e(jj,ii)/e(ii,ii) 5 e(jj,ii)=0.0D0 2 CONTINUE DO 7 ii=1, N-1 i2=N-ii+1 DO 8 jj=1,i2-1 j2=i2-jj DO 9 kk=1, N 9 e(j2,kk+N)=e(j2,kk+N)- 1e(i2,kk+N)*e(j2,i2)/e(i2,i2) e(j2,i2)=0.0D0 8 CONTINUE 7 CONTINUE DO 10 ii=1,N DO 12 jj=1,N 12 iA(ii,jj)=e(ii,jj+N)/e(ii,ii) 10 CONTINUE RETURN END SUBROUTINE WRITEFF(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION FCAR(N,N) OPEN(2,FILE='INTY.FC') WRITE(2,*)' Force constants in internal coordinates:' 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) close(2) 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 subroutine report(s) character*(*) s write(6,*)s stop 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 setx(IR,NA,xi,X,u,ia1) implicit none integer*4 IR,NA,ia,ia1,N real*8 xi(3),X(3,NA/3),a1,a2,a3,u(3,3) if(IR.lt.0)then c number of atom in the 3x3x3 cell: ia=abs(IR) c ia1: corresponding atom in the 000 cell call set123(a1,a2,a3,ia,NA/3,ia1) do 1 N=1,3 1 xi(N)=X(N,ia1)+u(1,N)*a1+u(2,N)*a2+u(3,N)*a3 else ia1=IR do 2 N=1,3 2 xi(N)=X(N,IR) endif return end subroutine set123(a1,a2,a3,a,N,ia1) implicit none integer*4 i1,i2,i3,a,N,ia1,ncube real*8 a1,a2,a3 c number of the cell (0..26): ncube=(a-1)/N c i1,i2,i3: elementary shifts (-1,0 or 1): if(ncube.le.13)then i1=(ncube-10)/9 i2=(ncube-4-9*(i1+1))/3 i3=ncube-9*i1-3*i2-14 else i1=(ncube-9)/9 i2=(ncube-3-9*(i1+1))/3 i3=ncube-9*i1-3*i2-13 endif ia1=a-ncube*N a1=dble(i1) a2=dble(i2) a3=dble(i3) return end SUBROUTINE XYZ(NOB0,NOB,NA,X,B) C C translations/rotations C IMPLICIT none INTEGER*4 NOB0,NOB,NA,i real*8 X(3,NA/3),B(NOB0,NA),p,xx,yy,zz,rr p=1.0d0/dble(NA/3) C c TX: do 1 i=1,NA/3 1 B(NOB,1+3*(i-1))=p c TY: NOB=NOB+1 do 2 i=1,NA/3 2 B(NOB,2+3*(i-1))=p c TZ: NOB=NOB+1 do 3 i=1,NA/3 3 B(NOB,3+3*(i-1))=p C c RX: NOB=NOB+1 do 4 i=1,NA/3 yy=X(2,i) zz=X(3,i) rr=yy**2+zz**2 if(rr.lt.1.0d-8)goto 4 B(NOB,2+3*(i-1))= p*zz/rr B(NOB,3+3*(i-1))=-p*yy/rr 4 continue c RY: NOB=NOB+1 do 5 i=1,NA/3 xx=X(1,i) zz=X(3,i) rr=xx**2+zz**2 if(rr.lt.1.0d-8)goto 5 B(NOB,3+3*(i-1))= p*xx/rr B(NOB,1+3*(i-1))=-p*zz/rr 5 continue c RZ: NOB=NOB+1 do 6 i=1,NA/3 xx=X(1,i) yy=X(2,i) rr=xx**2+yy**2 if(rr.lt.1.0d-8)goto 6 B(NOB,1+3*(i-1))= p*yy/rr B(NOB,2+3*(i-1))=-p*xx/rr 6 continue c RETURN END SUBROUTINE READFC(NA,FCAR) implicit none integer*4 NA,N,NTT,NTE,I,J,NT real*8 FCAR(NA,NA) OPEN(20,FILE='FILE.FC',STATUS='OLD') N=NA 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 return end subroutine rb(ic,NOB0,NA,B) implicit none integer*4 ic,NA,i,NOB0,nn,ii,ia,ix real*8 b(NOB0,NA) open(9,file='B.MAT') read(9,*)nob0 read(9,*)na if(ic.eq.0)goto 99 b=0.0d0 read(9,*) do 13 i=1,nob0 read(9,*) read(9,*)nn do 13 ii=1,nn 13 read(9,*)ia,ix,b(i,ix+3*(ia-1)) 99 close(9) return end subroutine sm(s,B,N) implicit none integer*4 N,i,j real*8 B(N,N) character*(*) s write(6,'(a)')s write(6,600) 600 format(6x,$) do 1 i=1,N 1 if(i.lt.4.or.N-i.lt.4)write(6,601)i 601 format(i12,$) write(6,*) do 2 i=1,N if(i.lt.4.or.N-i.lt.4)then write(6,602)i 602 format(i6,$) do 3 j=1,N 3 if(j.lt.4.or.N-j.lt.4)write(6,603)B(i,j) 603 format(e12.4,$) write(6,*) endif 2 continue return end