deritr.f

      SUBROUTINE DERITR(ERRFN,GEO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION GEO(3,NUMATM), ERRFN(MAXPAR)
************************************************************************
*
*    DERITR CALCULATES THE DERIVATIVES OF THE ENERGY WITH RESPECT TO THE
*          INTERNAL COORDINATES. THIS IS DONE BY FINITE DIFFERENCES
*          USING FULL SCF CALCULATIONS.
*
*          THIS IS VERY TIME-CONSUMING, AND SHOULD ONLY BE USED WHEN
*          NO OTHER DERIVATIVE CALCULATION WILL DO.
*
*    THE MAIN ARRAYS IN DERIV ARE:
*        LOC    INTEGER ARRAY, LOC(1,I) CONTAINS THE ADDRESS OF THE ATOM
*               INTERNAL COORDINATE LOC(2,I) IS TO BE USED IN THE
*               DERIVATIVE CALCULATION.
*        GEO    ARRAY \GEO\ HOLDS THE INTERNAL COORDINATES.
*
************************************************************************
      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 /GEOSYM/ NDEP, IDUMYS(MAXPAR,3)
      COMMON /UCELL / L1L,L2L,L3L,L1U,L2U,L3U
      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 /KEYWRD / KEYWRD
      CHARACTER*241 KEYWRD
      DIMENSION CHANGE(3), COORD(3,NUMATM)
     1,         XDERIV(3), XPARAM(MAXPAR), XJUC(3), W(N2ELEC)
      DOUBLE PRECISION WJ, WK
      SAVE IDELTA, XDERIV
      LOGICAL DEBUG
      EQUIVALENCE (W,WJ)
      DATA ICALCN /0/
      IF(ICALCN.NE.NUMCAL) THEN
         DEBUG = (INDEX(KEYWRD,'DERITR') .NE. 0)
         ICALCN=NUMCAL
*
*   IDELTA IS A MACHINE-PRECISION DEPENDANT INTEGER
*
         IDELTA=-3
         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    BECAUSE FULL SCF CALCULATIONS ARE BEING DONE QUITE LARGE STEPS
C    ARE NEEDED.  ON THE OTHER HAND, THE STEP CANNOT BE VERY LARGE,
C    AS THE SECOND DERIVITIVE IN FLEPO IS CALCULATED FROM THE
C    DIFFERENCES OF TWO FIRST DERIVATIVES. CHANGE(1) IS FOR CHANGE IN
C    BOND LENGTH, (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
      DO 10 I=1,NVAR
   10 XPARAM(I)=GEO(LOC(2,I),LOC(1,I))
      IF(NDEP.NE.0) CALL SYMTRY
      CALL GMETRY(GEO,COORD)
C
C  ESTABLISH THE ENERGY AT THE CURRENT POINT
C
      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
C
C  RESTORE THE DENSITY MATRIX (WHY?)
C
      DO 20 I=1,LINEAR
   20 P(I)=PA(I)*2.D0
      AA=(AA+ENUCLR)
      IJ=0
      DO 60 II=1,NUMAT
         DO 50 IL=L1L,L1U
            DO 50 JL=L2L,L2U
               DO 50 KL=L3L,L3U
                  DO 30 LL=1,3
   30             XJUC(LL)=COORD(LL,II)+TVEC(LL,1)*IL+TVEC(LL,2)*JL+TVEC
     1(LL,3)*KL
                  IJ=IJ+1
   50    CONTINUE
   60 CONTINUE
      DO 90 I=1,NVAR
         K=LOC(1,I)
         L=LOC(2,I)
         XSTORE=XPARAM(I)
         DO 70 J=1,NVAR
   70    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
C   IF NEEDED, CALCULATE "EXACT" DERIVITIVES.
C
         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 80 II=1,LINEAR
   80    P(II)=PA(II)*2.D0
         EE=(EE+ENUCLR)
         ERRFN(I)=(AA-EE)*23.061D0*XDERIV(L)*2.D0
   90 CONTINUE
      IF(DEBUG)THEN
         WRITE(6,'('' ERROR FUNCTION'')')
         WRITE(6,'(10F8.3)')(ERRFN(I),I=1,NVAR)
      ENDIF
      RETURN
      END