SUBROUTINE FMAT(FMATRX, NREAL, TSCF, TDER, DELDIP, HEAT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE 'SIZES' DIMENSION FMATRX(*), DELDIP(3,*) *********************************************************************** * * VALUE CALCULATES THE SECOND-ORDER OF THE ENERGY WITH * RESPECT TO THE CARTESIAN COORDINATES I AND J AND PLACES IT * IN FMATRX * * ON INPUT NATOMS = NUMBER OF ATOMS IN THE SYSTEM. * XPARAM = INTERNAL COORDINATES OF MOLECULE STORED LINEARLY * * VARIABLES USED * COORDL = ARRAY OF CARTESIAN COORDINATES, STORED LINEARLY. * I = INDEX OF CARTESIAN COORDINATE. * J = INDEX OF CARTESIAN COORDINATE. * * ON OUTPUT FMATRX = SECOND DERIVATIVE OF THE ENERGY WITH RESPECT TO * CARTESIAN COORDINATES I AND J. *********************************************************************** COMMON /KEYWRD/ KEYWRD COMMON /GEOKST/ NATOMS,LABELS(NUMATM), 1 NA(NUMATM),NB(NUMATM),NC(NUMATM) COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR),IDUMY, DUMY(MAXPAR) COMMON /DENSTY/ P(MPACK),PDUMY(2,MPACK) COMMON /TIMDMP/ TLEFT, TDUMP COMMON /ATMASS/ ATMASS(NUMATM) COMMON /TIME / TIME0 COMMON /CORE / CORE(107) COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 2 NCLOSE,NOPEN,NDUMY,FRACT COMMON /COORD / COORD(3,NUMATM) COMMON /NLLCOM/ EVECS(MAXPAR*MAXPAR),BMAT(MAXPAR,MAXPAR*2) DIMENSION COLD(MAXPAR), GRAD(MAXPAR), 1GROLD(MAXPAR), COORDL(MAXPAR), Q(NUMATM), DEL2(3), G2OLD(MAXPAR) 2, EIGS(MAXPAR), G2RAD(MAXPAR), 3 FCONST(MAXPAR) CHARACTER*241 KEYWRD SAVE SPACE, DOTT, ZERO, NINE, FACT LOGICAL DEBUG, RESTRT, PRNT, RESFIL, ANALYT, PRECIS, BIG, LOG CHARACTER SPACE*1, DOTT*1, ZERO*1, NINE*1 EQUIVALENCE (COORD(1,1),COORDL(1)) DATA SPACE,DOTT,ZERO,NINE /' ','.','0','9'/ DATA FACT/6.95125D-3/ C C FACT IS THE CONVERSION FACTOR FROM KCAL/MOLE TO ERGS C C SET UP CONSTANTS AND FLAGS NA(1)=99 C C SET UP THE VARIABLES IN XPARAM ANDLOC,THESE ARE IN CARTESIAN COORDINA C NUMAT=0 C$DOIT ASIS DO 10 I=1,NATOMS IF(LABELS(I).NE.99.AND.LABELS(I).NE.107) THEN NUMAT=NUMAT+1 LABELS(NUMAT)=LABELS(I) ENDIF 10 CONTINUE NATOMS=NUMAT C C THIS IS A QUICK, IF CLUMSY, WAY TO CALCULATE NUMAT, AND TO REMOVE C THE DUMMY ATOMS FROM THE ARRAY LABELS. C NVAR=NUMAT*3 DO 20 I=1,NUMAT LOC(1,(I-1)*3+1)=I LOC(2,(I-1)*3+1)=1 C LOC(1,(I-1)*3+2)=I LOC(2,(I-1)*3+2)=2 C LOC(1,(I-1)*3+3)=I LOC(2,(I-1)*3+3)=3 20 CONTINUE LIN=(NVAR*(NVAR+1))/2 DO 30 I=1,LIN 30 FMATRX(I)=0.D0 PRNT =(INDEX(KEYWRD,'IRC=') .EQ. 0) LOG =(INDEX(KEYWRD,'NOLOG') .EQ. 0) PRECIS =(INDEX(KEYWRD,'PREC') .NE. 0) ANALYT =(INDEX(KEYWRD,'ANALYT') .NE. 0) RESTRT =(INDEX(KEYWRD,'RESTART') .NE. 0) IF(INDEX(KEYWRD,'NLLSQ') .NE. 0) RESTRT=.FALSE. DEBUG =(INDEX(KEYWRD,'FMAT') .NE. 0) BIG =(INDEX(KEYWRD,'LARGE') .NE. 0 .AND. DEBUG) IF(PRNT)WRITE(6,'(//4X,''FIRST DERIVATIVES WILL BE USED IN THE'' 1,'' CALCULATION OF SECOND DERIVATIVES'')') TLAST=TLEFT RESFIL=.FALSE. IF(RESTRT) THEN DO 40 I=1,NVAR 40 COLD(I)=COORDL(I) ISTART = 0 I=0 CALL FORSAV(TOTIME,DELDIP,ISTART,I,FMATRX, COORD, NVAR,HEAT, 1 EVECS,JSTART,FCONST) KOUNTF=(ISTART*(ISTART+1))/2 ISTART=ISTART+1 JSTART=JSTART+1 TIME2 = SECOND() ELSE KOUNTF=0 TOTIME=0.D0 IF (TSCF.GT.0.D0)TLEFT=TLEFT-TSCF-TDER ISTART=1 ENDIF C CALCULATE FMATRX IF(ISTART.GT.1) THEN ESTIME=(NVAR-ISTART+1)*TOTIME/(ISTART-1.D0) ELSE ESTIME=NVAR*(TSCF+TDER)*2.D0 IF (PRECIS) ESTIME=ESTIME*2.D0 ENDIF IF(TSCF.GT.0) 1WRITE(6,'(/10X,''ESTIMATED TIME TO COMPLETE CALCULATION ='' 2,F9.2,'' SECONDS'')')ESTIME IF(RESTRT) THEN IF(ISTART.LE.NVAR) 1 WRITE(6,'(/10X,''STARTING AGAIN AT LINE'',18X,I4)')ISTART WRITE(6,'(/10X,''TIME USED UP TO RESTART ='',F22.2)')TOTIME ENDIF LU=KOUNTF TIME1 = SECOND() NUMAT=NVAR/3 DO 50 I=1,NVAR 50 EIGS(I)=0.D0 DO 120 I=ISTART,NVAR TIME2 = SECOND() DELTA=1.D0/120.D0 IF(PRECIS)THEN C C DETERMINE A GOOD STEP SIZE C G2OLD(1)=100.D0 COORDL(I)=COORDL(I)+DELTA CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.) COORDL(I)=COORDL(I)-DELTA DELTA=DELTA*10.D0/SQRT(DOT(G2OLD,G2OLD,NVAR)) C C CONSTRAIN DELTA TO A 'REASONABLE' VALUE C DELTA=MIN(0.05D0,MAX(0.005D0,DELTA)) IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT 1A G2OLD(1)=100.D0 COORDL(I)=COORDL(I)+DELTA CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., G2OLD, .TRUE.) IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +1.0*DELTA', 1SQRT(DOT(G2OLD,G2OLD,NVAR)) COORDL(I)=COORDL(I)-DELTA*2.D0 G2RAD(1)=100.D0 CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., G2RAD, .TRUE.) COORDL(I)=COORDL(I)+DELTA IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -1.0*DELTA', 1SQRT(DOT(G2RAD,G2RAD,NVAR)) ELSE IF(DEBUG)WRITE(6,'(A,I3,A,F12.5)')' STEP:',I,' DELTA :',DELT 1A ENDIF COORDL(I)=COORDL(I)+0.5D0*DELTA GROLD(1)=100.D0 CALL COMPFG(COORDL, .TRUE., ESCF, .TRUE., GROLD, .TRUE.) IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM +0.5*DELTA', 1SQRT(DOT(GROLD,GROLD,NVAR)) CALL CHRGE(P,Q) DO 60 II=1,NUMAT 60 Q(II)=CORE(LABELS(II))-Q(II) SUM = DIPOLE(P,Q,COORDL,DELDIP(1,I),0) COORDL(I)=COORDL(I)-DELTA GRAD(1)=100.D0 CALL COMPFG(COORDL, .TRUE., HEATAA, .TRUE., GRAD, .TRUE.) IF(DEBUG)WRITE(6,'(A,F12.5)')' GNORM -0.5*DELTA', 1SQRT(DOT(GRAD,GRAD,NVAR)) CALL CHRGE(P,Q) DO 70 II=1,NUMAT 70 Q(II)=CORE(LABELS(II))-Q(II) SUM = DIPOLE(P,Q,COORDL,DEL2,0) COORDL(I)=COORDL(I)+DELTA*0.5D0 DELDIP(1,I)=(DELDIP(1,I)-DEL2(1))*0.5D0/DELTA DELDIP(2,I)=(DELDIP(2,I)-DEL2(2))*0.5D0/DELTA DELDIP(3,I)=(DELDIP(3,I)-DEL2(3))*0.5D0/DELTA LL=LU+1 LU=LL+I-1 L=0 IF(PRECIS)THEN DO 80 KOUNTF=LL,LU L=L+1 C C G2OLD = X + 1.0*DELTA C GROLD = X + 0.5*DELTA C GRAD = X - 0.5*DELTA C G2RAD = X - 1.0*DELTA C DUMY(L)= (8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L))) 1 /DELTA*FACT/24.D0 EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L))) 1 /DELTA**3*FACT/56.D0 C C CORRECT FOR 4'TH ORDER CONTAMINATION C C# CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0) C# DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L)) FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L) 80 CONTINUE L=L-1 DO 90 K=I,NVAR L=L+1 KK=(K*(K-1))/2+I DUMY(L)=(8.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L))) 1 /DELTA*FACT/24.D0 EIGS(L)=(2.D0*(GROLD(L)-GRAD(L))-(G2OLD(L)-G2RAD(L))) 1 /DELTA**3*FACT/56.D0 C C CORRECT FOR 4'TH ORDER CONTAMINATION C C# CORR=MIN(ABS(DUMY(L)),ABS(EIGS(L))*0.0001D0) C# DUMY(L)=DUMY(L)-SIGN(CORR,DUMY(L)) FMATRX(KK)=FMATRX(KK)+DUMY(L) 90 CONTINUE ELSE DO 100 KOUNTF=LL,LU L=L+1 DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT FMATRX(KOUNTF)=FMATRX(KOUNTF)+DUMY(L) 100 CONTINUE L=L-1 DO 110 K=I,NVAR L=L+1 KK=(K*(K-1))/2+I DUMY(L)=((GROLD(L)-GRAD(L)))*0.25D0/DELTA*FACT FMATRX(KK)=FMATRX(KK)+DUMY(L) 110 CONTINUE ENDIF IF(BIG)THEN WRITE(6,'(A)')' CONTRIBUTIONS TO F-MATRIX' WRITE(6,'(A)')' ELEMENT +1.0*DELTA +0.5*DELTA -0.5*DEL' 1//'TA -1.0*DELTA 2''ND ORDER 4TH ORDER' WRITE(6,'(I7,6F12.6)')(L,G2OLD(L),GROLD(L),GRAD(L),G2RAD(L), 1DUMY(L),EIGS(L),L=1,NVAR) ENDIF TIME3 = SECOND() TSTEP=TIME3-TIME2 TLEFT= MAX(0.1D0,TLEFT-TSTEP) IF(TSTEP.GT.1.D6)TSTEP=TSTEP-1.D6 TOTIME= TOTIME+TSTEP IF(RESFIL)THEN WRITE(6,'('' STEP:'',I4,'' RESTART FILE WRITTEN, INTEGRAL =' 1',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT IF(LOG)WRITE(11,'('' STEP:'',I4,'' RESTART FILE WRITTEN, '', +''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TOTIME,TLEFT RESFIL=.FALSE. ELSE WRITE(6,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, INTEGRAL = 1'',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT IF(LOG) WRITE(11,'('' STEP:'',I4,'' TIME ='',F9.2,'' SECS, '', +''INTEGRAL ='',F10.2,'' TIME LEFT:'',F10.2)')I,TSTEP,TOTIME,TLEFT ENDIF ESTIM = TOTIME/I IF(TLAST-TLEFT.GT.TDUMP)THEN TLAST=TLEFT RESFIL=.TRUE. JSTART=1 II=I CALL FORSAV(TOTIME,DELDIP,II,NVAR,FMATRX, COORD,NVAR,HEAT, 1 EVECS,JSTART,FCONST) ENDIF IF(I.NE.NVAR.AND.TLEFT-10.D0 .LT. ESTIM) THEN WRITE(6,'(//10X,''- - - - - - - TIME UP - - - - - - -'',//)' 1) WRITE(6,'(/10X,'' POINT REACHED ='',I4)')I WRITE(6,'(/10X,'' RESTART USING KEY-WORD "RESTART"'')') WRITE(6,'(10X,''ESTIMATED TIME FOR THE NEXT STEP ='',F8.2, 1'' SECONDS'')')ESTIM JSTART=1 II=I CALL FORSAV(TOTIME,DELDIP,II,NVAR,FMATRX, COORD,NVAR,HEAT, 1 EVECS,JSTART,FCONST) WRITE(6,'(//10X,''FORCE MATRIX WRITTEN TO DISK'')') NREAL=-1 RETURN ENDIF 120 CONTINUE DO 130 I=1,NATOMS IF(ATMASS(I).LT.1.D-20.AND.LABELS(I).LT.99)THEN CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT, 1 EVECS,ILOOP,FCONST) WRITE(6,'(A)')' AT LEAST ONE ATOM HAS A ZERO MASS. A RESTART 1' WRITE(6,'(A)')' FILE HAS BEEN WRITTEN AND THE JOB STOPPED' STOP ENDIF 130 CONTINUE IF(ISTART.LE.NVAR .AND. INDEX(KEYWRD,'ISOT') .NE. 0) 1CALL FORSAV(TOTIME,DELDIP,NVAR,NVAR,FMATRX, COORD,NVAR,HEAT, 2 EVECS,ILOOP,FCONST) RETURN END