SUBROUTINE DERIV(GEO,GRAD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INCLUDE "SIZES" DIMENSION GRAD(*), GEO(3,*) COMMON / EULER/ TVEC(3,3), ID COMMON /GEOVAR/ NVAR,LOC(2,MAXPAR), IDUMY, DUMMY(MAXPAR) COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM), 1 NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA, 2 NCLOSE,NOPEN,NDUMY,FRACT COMMON /GEOKST/ NATOMS,LABELS(NUMATM), 1NA(NUMATM),NB(NUMATM),NC(NUMATM) COMMON /GRAVEC/ COSINE COMMON /GEOSYM/ NDEP, IDUMYS(MAXPAR,3) COMMON /PATH / LATOM,LPARAM,REACT(200) COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U COMMON /XYZGRA/ DXYZ(3,NUMATM*27) COMMON /ENUCLR/ ENUCLR COMMON /NUMCAL/ NUMCAL COMMON /DENSTY/ P(MPACK), PA(MPACK), PB(MPACK) COMMON /WMATRX/ WJ(N2ELEC), WK(N2ELEC) COMMON /HMATRX/ H(MPACK) COMMON /ATHEAT/ ATHEAT C*********************************************************************** C C DERIV CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE C INTERNAL COORDINATES. THIS IS DONE BY FINITE DIFFERENCES. C C THE MAIN ARRAYS IN DERIV ARE: C LOC INTEGER ARRAY, LOC(1,I) CONTAINS THE ADDRESS OF THE ATOM C INTERNAL COORDINATE LOC(2,I) IS TO BE USED IN THE C DERIVATIVE CALCULATION. C GEO ARRAY \GEO\ HOLDS THE INTERNAL COORDINATES. C GRAD ON EXIT, CONTAINS THE DERIVATIVES C C*********************************************************************** COMMON /KEYWRD / KEYWRD COMMON /ERRFN / ERRFN(MAXPAR) CHARACTER*80 KEYWRD DIMENSION CHANGE(3), COORD(3,NUMATM), COLD(3,NUMATM*27) 1, XDERIV(3), XPARAM(MAXPAR), XJUC(3), W(N2ELEC) REAL WJ, WK LOGICAL DEBUG, TIMES, HALFE, FAST, SCF1, CI, PRECIS EQUIVALENCE (W,WJ) DATA ICALCN /0/ IF(ICALCN.NE.NUMCAL) THEN I = INDEX(KEYWRD,'PRESS') PRESS=0.D0 IF(I.NE.0) PRESS=READA(KEYWRD,I)*1476.8992D0 IDLO=NATOMS+1 IF(LABELS(NATOMS) .EQ. 107) THEN IDLO=NATOMS IF(LABELS(NATOMS-1) .EQ. 107)THEN IDLO=NATOMS-1 IF(LABELS(NATOMS-2) .EQ. 107)THEN IDLO=NATOMS-2 ENDIF ENDIF ENDIF GRLIM=0.01D0 DEBUG = (INDEX(KEYWRD,'DERIV') .NE. 0) PRECIS= (INDEX(KEYWRD,'PRECIS') .NE. 0) TIMES = (INDEX(KEYWRD,'TIME') .NE. 0) CI = (INDEX(KEYWRD,'C.I.') .NE. 0) SCF1 = (INDEX(KEYWRD,'1SCF') .NE. 0) ICALCN=NUMCAL IF(INDEX(KEYWRD,'RESTART') .EQ. 0) THEN DO 10 I=1,NVAR 10 ERRFN(I)=0.D0 ENDIF GRLIM=0.01D0 IF(PRECIS)GRLIM=0.0001D0 IF(INDEX(KEYWRD,'FULSCF') .GT.0) GRLIM=1.D9 HALFE = (NOPEN.GT.NCLOSE .OR. CI) IDELTA=-7 * * IDELTA IS A MACHINE-PRECISION DEPENDANT INTEGER * IF(HALFE.AND.PRECIS) IDELTA=-3 IF(HALFE.AND..NOT.PRECIS) IDELTA=-3 FAST=.TRUE. CHANGE(1)= 10.D0**IDELTA CHANGE(2)= 10.D0**IDELTA CHANGE(3)= 10.D0**IDELTA C C CHANGE(I) IS THE STEP SIZE USED IN CALCULATING THE DERIVATIVES. C FOR "CARTESIAN" DERIVATIVES, CALCULATED USING DCART,AN C INFINITESIMAL STEP, HERE 0.000001, IS ACCEPTABLE. IN THE C HALF-ELECTRON METHOD A QUITE LARGE STEP IS NEEDED AS FULL SCF C CALCULATIONS ARE NEEDED, AND THE DIFFERENCE BETWEEN THE TOTAL C ENERGIES IS USED. THE STEP CANNOT BE VERY LARGE, AS THE SECOND C DERIVITIVE IN FLEPO IS CALCULATED FROM THE DIFFERENCES OF TWO C FIRST DERIVATIVES. CHANGE(1) IS FOR CHANGE IN BOND LENGTH, C (2) FOR ANGLE, AND (3) FOR DIHEDRAL. C XDERIV(1)= 0.5D0/CHANGE(1) XDERIV(2)= 0.5D0/CHANGE(2) XDERIV(3)= 0.5D0/CHANGE(3) ENDIF GNORM=0.D0 IF(NVAR.EQ.0) RETURN IF(DEBUG)THEN WRITE(6,'('' GEO AT START OF DERIV'')') WRITE(6,'(F19.5,2F12.5)')((GEO(J,I),J=1,3),I=1,NATOMS) ENDIF DO 20 I=1,NVAR XPARAM(I)=GEO(LOC(2,I),LOC(1,I)) 20 GNORM=GNORM+GRAD(I)**2 GNORM=SQRT(GNORM) FAST=(GNORM .GT. GRLIM .AND. .NOT. SCF1 .OR. .NOT. HALFE) TIME1=SECOND() IF(NDEP.NE.0) CALL SYMTRY CALL GMETRY(GEO,COORD) IF( .NOT. FAST ) THEN IF(DEBUG)WRITE(6,'('' DOING FULL SCF''''S IN DERIV'')') CALL HCORE(COORD,H,W, WJ, WK, ENUCLR) IF(NORBS*NELECS.GT.0)THEN CALL ITER(H, W, WJ, WK, AA,.TRUE.,.FALSE.) ELSE AA=0.D0 ENDIF LINEAR=(NORBS*(NORBS+1))/2 DO 30 I=1,LINEAR 30 P(I)=PA(I)*2.D0 AA=(AA+ENUCLR) ENDIF CALL DCART(COORD,DXYZ) IF(NDEP.NE.0) CALL SYMTRY CALL GMETRY(GEO,COORD) IJ=0 DO 70 II=1,NUMAT DO 60 IL=L1L,L1U DO 60 JL=L2L,L2U DO 60 KL=L3L,L3U DO 40 LL=1,3 40 XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*JL+TVEC 1(LL,3)*KL IJ=IJ+1 DO 50 KK=1,3 COLD(KK,IJ)=XJUC(KK) 50 CONTINUE 60 CONTINUE 70 CONTINUE SUM11=1.D-9 SUM22=1.D-9 SUM12=1.D-9 DO 150 I=1,NVAR K=LOC(1,I) L=LOC(2,I) XSTORE=XPARAM(I) DO 80 J=1,NVAR 80 GEO(LOC(2,J),LOC(1,J))=XPARAM(J) GEO(L,K)=XSTORE-CHANGE(L) IF(NDEP.NE.0) CALL SYMTRY CALL GMETRY(GEO,COORD) C# CALL GEOUT C C USE LOOKUP TABLE OF CARTESIAN DERIVATIVES TO WORK OUT INTERNAL C COORDINATE DERIVATIVE. C TOTL=0.D0 IJ=0 DO 130 II=1,NUMAT IF(ID.EQ.0) THEN DO 90 LL=1,3 90 TOTL=TOTL+DXYZ(LL,II)*(COORD(LL,II)-COLD(LL,II)) ELSE DO 120 IL=L1L,L1U DO 120 JL=L2L,L2U DO 120 KL=L3L,L3U DO 100 LL=1,3 100 XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*J 1L+TVEC(LL,3)*KL IJ=IJ+1 DO 110 KK=1,3 TOTL=TOTL+DXYZ(KK,IJ)*(XJUC(KK)-COLD(KK,IJ)) 110 CONTINUE 120 CONTINUE ENDIF 130 CONTINUE TOTL=TOTL*XDERIV(L) C C IF NEEDED, CALCULATE "EXACT" DERIVITIVES. C IF( .NOT. FAST ) THEN CALL HCORE(COORD,H,W, WJ, WK,ENUCLR) IF(NORBS*NELECS.GT.0)THEN CALL ITER(H,W, WJ, WK,EE,.TRUE.,.FALSE.) ELSE EE=0.D0 ENDIF DO 140 II=1,LINEAR 140 P(II)=PA(II)*2.D0 EE=(EE+ENUCLR) TOTL1=(AA-EE)*23.061D0*XDERIV(L)*2.D0 C# WRITE(6,*)AA-EE ERRFN(I)=TOTL1-TOTL ENDIF GEO(L,K)=XSTORE SUM11=SUM11+GRAD(I)**2 SUM22=SUM22+TOTL**2 SUM12=SUM12+TOTL*GRAD(I) GRAD(I)=TOTL+ERRFN(I) 150 CONTINUE IF(DEBUG) THEN WRITE(6,'('' GRADIENTS'')') WRITE(6,'(10F8.3)')(GRAD(I),I=1,NVAR) WRITE(6,'('' ERROR FUNCTION'')') WRITE(6,'(10F8.3)')(ERRFN(I),I=1,NVAR) ENDIF COSINE=SUM12/SQRT(SUM11*SUM22) IF(DEBUG) 1WRITE(6,'('' COSINE OF SEARCH DIRECTION ='',F30.6)')COSINE IF( .NOT. FAST ) COSINE=1.D0 IF(TIMES) 1WRITE(6,'('' TIME FOR DERIVATIVES'',F12.6)')SECOND()-TIME1 RETURN END