      SUBROUTINE DERNVO(COORD,DXYZ)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION COORD(3,*), DXYZ(*)
***********************************************************************
*
*    IMPLEMENTATION OF ANALYTICAL FORMULATION FOR OPEN SHELL OR CI,
*                      VARIABLES FINITE DIFFERENCE METHODS,
*                      STATISTICAL ESTIMATE OF THE ERRORS,
*                   BY D. LIOTARD
*                      LABORATOIRE DE CHIMIE STRUCTURALE
*                      UNIVERSITE DE PAU ET DES PAYS DE L'ADOUR
*                      AVENUE DE L'UNIVERSITE, 64000, PAU (FRANCE)
*
*
*   MODIFIED BY JJPS TO CONFORM TO MOPAC CONVENTIONS
*   (NOTE BY JJPS:  PROF. LIOTARD'S TECHNIQUE WORKS.  IF THIS
*   IMPLEMENTATION DOES NOT WORK, THE REASON IS A FAULT INTRODUCED
*   BY JJPS, AND DOES NOT REFLECT ON PROF. LIOTARD'S ABILITY)
*
*
*    AS THE WAVE FUNCTION IS NOT VARIATIONALLY OPTIMIZED, I.E.
*    HALF-ELECTRON OR CI, THE DERIVATIVES OF THE 1 AND 2-ELECTRON
*    INTEGRALS IN A.O. BASIS ARE EVALUATED IN CARTESIAN COORDINATES
*    BY A 1 OR 2 POINTS FINITE DIFFERENCE FORMULA AND STORED.
*    THUS ONE GETS THE NON-RELAXED (I.E. FROZEN ELECTRONIC CLOUD)
*    CONTRIBUTION TO THE FOCK EIGENVALUES AND 2-ELECTRON INTEGRALS IN
*    AN M.O. BASIS.  THE NON-RELAXED GRADIENT COMES FROM THE
*    NON-RELAXED C.I. MATRIX DERIVATIVE (SUBROUTINE DERI1).
*    THE DERIVATIVES OF THE M.O. COEFFICIENTS ARE THEN WORKED OUT
*    ITERATIVELY (OK FOR BOTH CLOSED SHELLS AND HALF-ELECTRON CASES)
*    AND STORED. THUS ONE GETS THE ELECTRONIC RELAXATION CONTRIBUTION TO
*    THE FOCK EIGENVALUES AND 2-ELECTRON INTEGRALS IN M.O. BASIS.
*    FINALLY THE RELAXATION CONTRIBUTION TO THE C.I. MATRIX DERIVATIVE
*    GIVES THE RELAXATION CONTRIBUTION TO THE GRADIENT (ROUTINE DERI2).
*
*
*        COORD  HOLDS THE CARTESIAN COORDINATES.
*    INPUT
*        DXYZ   NOT DEFINED.
*    EXIT
*        DXYZ   DERIVATIVES OF ENERGY W.R.T CARTESIAN COORDINATES,
*               IN KCAL/MOL/ANGSTROM (3 * NUMAT OF THESE)
*
***********************************************************************
C
C    NW2 and NW3 should be set to *something*, probably sizes of WORK2 
C    and WORK3 arrays? Since the code works with NW2 and NW3 set to zero
C    by implicit initialization, we can as well do it explicitly.... -P.S.
C
      PARAMETER (NW2=0,NW3=0)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM)
     1               ,NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA
     2               ,NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /GRADNT/ GRAD(MAXPAR),GNORM
      COMMON /CIBITS/ NMOS,LAB,NELEC,NBO(3)
      COMMON /NUMCAL/ NUMCAL
     1       /KEYWRD/ KEYWRD
     2       /VECTOR/ C(MORB2),EIGS(MAXORB),CBETA(MORB2)
     3,EIGB(MAXORB)
      COMMON /FOKMAT/ FDUMY(MPACK), SCALAR(MPACK)
      COMMON /NVOMAT/ DIAG(MPACK/2)
      COMMON /WORK1 / FMOOFF(NPULAY*4), FMOON(NPULAY*4),
     1WORK2(9*NPULAY), WORK3(4*NPULAY)
      DIMENSION FBWO(5*MAXPAR)
      CHARACTER KEYWRD*241, BLANK*60
      DIMENSION DXYZR(MAXPAR), EIGBB(6*MAXPAR)
      LOGICAL DEBUG, DCAR, LARGE, RELAXD, FORCE
      DATA ICALCN /0/
C
C     SELECT THE REQUIRED OPTION AND READ KEYWORDS
C     --------------------------------------------
C
      IF(ICALCN.NE.NUMCAL) THEN
         THROLD=0.08D0
         DEBUG = (INDEX(KEYWRD,'DERNVO') .NE. 0)
         LARGE = (INDEX(KEYWRD,'LARGE') .NE. 0)
         FORCE = (INDEX(KEYWRD,'FORC') .NE. 0)
         DCAR  = (INDEX(KEYWRD,'FORC') + INDEX(KEYWRD,'PREC') .NE. 0)
         IF(DCAR)THROLD=0.004D0
         DO 11 I=1,NVAX
   11    DXYZR(I)=0.D0
C        ACTUAL SIZES FOR C.I. GRADIENT CALCULATION.
         NBO(1)=NCLOSE
         NBO(2)=NOPEN-NCLOSE
         NBO(3)=NORBS-NOPEN
         MINEAR=NBO(2)*NBO(1)+NBO(3)*NOPEN
         NINEAR=(NMOS*(NMOS+1))/2+1
         ICALCN = NUMCAL
      ENDIF
C        SCALING ROW FACTORS TO SPEED CV OF RELAXATION PROCEDURE.
C#      CALL TIMER('BEFORE DERI0')
      CALL DERI0 (EIGS,NORBS,SCALAR,DIAG,FRACT,NBO)
C#      CALL TIMER('AFTER DERI0')
      NVAX=3*NUMAT
C
C   BECAUSE DERI2 IS CPU INTENSIVE, AND THE CONTRIBUTION TO THE
C   DERVIATIVE DUE TO RELAXATION OF THE ELECTRON CLOUD IS RELATIVELY
C   INSENSITIVE TO CHANGES IN GEOMETRY, WHERE POSSIBLE ONLY CALCULATE
C   THE DERIVATIVE EVERY 2 CALLS TO DERNVO
C
      SUM=0.D0
      IF(DCAR)THEN
         DO 10 I=1,NVAX
   10    DXYZR(I)=0.D0
         RELAXD=.FALSE.
      ENDIF
      DO 20 I=1,NVAX
   20 SUM=SUM+ABS(DXYZR(I))
      RELAXD=(SUM.GT.1.D-7)
C
C  IF DXYZR CONTAINS DATA, USE IT AND FLUSH AFTER USE.
C
      ILAST=0
   30 IFIRST=ILAST+1
      J=2
      IF(MIN(NW2,NW3)/MAX(MINEAR,NINEAR).LT.10)J=1
      ILAST=MIN(NVAX,ILAST+J)
      J=1-MINEAR
      K=1-NINEAR
      DO 40 I=IFIRST,ILAST
         K=K+NINEAR
         J=J+MINEAR
C
C        NON-RELAXED CONTRIBUTION (FROZEN ELECTRONIC CLOUD) IN DXYZ
C        AND NON-RELAXED FOCK MATRICES IN FMOOFF AND FMOON.
C   CONTENTS OF F-MO-OFF: OPEN-CLOSED, VIRTUAL-CLOSED, AND VIRTUAL-OPEN
C   CONTENTS OF F-MO-ON:  CLOSED-CLOSED, OPEN-OPEN AND VIRTUAL-VIRTUAL
C   OVER M.O. INDICES
C
C#      CALL TIMER('BEFORE DERI1')
         CALL DERI1(C,NORBS,COORD,I,CBETA,DXYZ(I),FMOOFF(J),MINEAR
     1              ,FMOON(K),WORK2,WORK2(6*MPACK),WORK3)
C#      CALL TIMER('AFTER DERI1')
   40 CONTINUE
      IF(DEBUG)THEN
         IF(IFIRST.EQ.1.AND.LARGE)THEN
            WRITE(6,*)' CONTENTS OF FMOOFF '
            WRITE(6,*)' OPEN-CLOSED'
            WRITE(6,'(7X,I3,5I12)')(J,J=NCLOSE+1,NOPEN)
            DO 50 I=1,NCLOSE
   50       WRITE(6,'(I3,6F12.6)')I,(FMOOFF(J),J=(I-1)*NBO(2)+1,I*NBO(2)
     1)
C
C
            WRITE(6,*)' VIRTUAL-CLOSED'
            K=NCLOSE*NBO(2)
            WRITE(6,'(7X,I3,5I12)')(J,J=NOPEN+1,MIN(NOPEN+6,NORBS))
            DO 60 I=1,NCLOSE
   60       WRITE(6,'(I3,6F12.6)')I,
     1   (FMOOFF(J+K),J=(I-1)*NBO(3)+1,MIN(6+(I-1)*NBO(3),I*NBO(3)))
            K=NCLOSE*NBO(2)+NBO(3)*NCLOSE
C
C
            WRITE(6,*)' VIRTUAL-OPEN'
            WRITE(6,'(7X,I3,4I12)')(J,J=NCLOSE+1,NOPEN)
            DO 70 I=1,MIN(6,NBO(3))
   70       WRITE(6,'(I3,6F12.6)')I+NOPEN,
     1   (FMOOFF(J+K),J=(I-1)*NBO(2)+1,MIN((I-1)*NBO(2)+6,I*NBO(2)))
            WRITE(6,*)' CONTENTS OF FMOON (ACTIVE-SPACE -- ACTIVE SPACE)
     1'
            K=(NMOS*(NMOS-1))/2
            LL=1
            BLANK=' '
            DO 80 I=1,NMOS
               L=LL+NMOS-I-1
               WRITE(6,'(A,5F12.6)')BLANK(:12*I),(FMOON(J),J=LL,L),FMOON
     1(K+I)
   80       LL=L+1
         ENDIF
      ENDIF
C        COMPUTE THE ELECTRONIC RELAXATION CONTRIBUTION.
C
C   DERNVO PROVIDES THE FOLLOWING SCRATCH AREAS TO DERI2: EIGB, WORK2,
C          WORK3, FBWO, CBETA.  THESE ARE DIMENSIONED ON ENTRY TO DERI2
C          WHICH IS WHY THEY ARE NOT DECLARED THERE.  THEY ARE NOT USED
C          AT ALL IN DERNVO.
C
C#      CALL TIMER('BEFORE DERI2')
      IF(.NOT.RELAXD)
     1   CALL DERI2 (C,EIGS,NORBS,MINEAR,FMOOFF
     2              ,FMOON,EIGBB, NINEAR,ILAST-IFIRST+1
     3              ,CBETA,WORK2,NW2,DXYZR(IFIRST)
     4              ,WORK3,NW3,FBWO,THROLD)
C#      CALL TIMER('AFTER DERI2')
      IF (ILAST.LT.NVAX) GO TO 30
      IF(DEBUG)THEN
         SUMX=0.D0
         SUMY=0.D0
         SUMZ=0.D0
         DO 90 I=1,NUMAT
            SUMX=SUMX+DXYZ(I*3-2)
            SUMY=SUMY+DXYZ(I*3-1)
   90    SUMZ=SUMZ+DXYZ(I*3)
         WRITE(6,*)' CARTESIAN DERIVATIVES DUE TO FROZEN CORE'
         WRITE(6,'('' ATOM    X           Y           Z'')')
         DO 100 I=1,NUMAT
  100    WRITE(6,'(I4,3F12.7)')I,DXYZ(I*3-2),DXYZ(I*3-1),DXYZ(I*3)
         WRITE(6,'(/10X,''RESIDUAL ERROR'')')
         WRITE(6,'(4X,3F12.7)')SUMX,SUMY,SUMZ
         WRITE(6,*)
         SUMX=0.D0
         SUMY=0.D0
         SUMZ=0.D0
         DO 110 I=1,NUMAT
            SUMX=SUMX+DXYZR(I*3-2)
            SUMY=SUMY+DXYZR(I*3-1)
  110    SUMZ=SUMZ+DXYZR(I*3)
         WRITE(6,*)' CARTESIAN DERIVATIVES DUE TO RELAXING CORE'
         WRITE(6,'('' ATOM    X           Y           Z'')')
         DO 120 I=1,NUMAT
  120    WRITE(6,'(I4,3F12.7)')I,DXYZR(I*3-2),DXYZR(I*3-1),DXYZR(I*3)
         WRITE(6,'(/10X,''RESIDUAL ERROR'')')
         WRITE(6,'(4X,3F12.7)')SUMX,SUMY,SUMZ
         WRITE(6,*)
      ENDIF
      DO 130 I=1,NVAX
  130 DXYZ(I)=DXYZ(I)+DXYZR(I)
      IF(RELAXD)THEN
         DO 140 I=1,NVAX
  140    DXYZR(I)=0.D0
      ENDIF
      SUMX=0.D0
      SUMY=0.D0
      SUMZ=0.D0
      DO 150 I=1,NUMAT
         SUMX=SUMX+DXYZ(I*3-2)
         SUMY=SUMY+DXYZ(I*3-1)
  150 SUMZ=SUMZ+DXYZ(I*3)
      SUM=MAX(1.D-10,ABS(SUMX)+ABS(SUMY)+ABS(SUMZ))
C
C  HERE IS A ROUGH BUT SIMPLE METHOD FOR DEFINING THROLD FOR DERI2
C  IT MAY NEED MORE WORK
C
      IF(.NOT. FORCE .AND. GNORM .GT. 0.001D0)
     1 THROLD=THROLD*SQRT(GNORM/(SUM*100.D0))
      THROLD=MIN(2.D0,MAX(0.002D0,THROLD))
      IF(DEBUG)THEN
         WRITE(6,*)'CARTESIAN DERIVATIVES FROM ANALYTICAL C.I. CALCULATI
     1ON'
         WRITE(6,'('' ATOM    X           Y           Z'')')
         DO 160 I=1,NUMAT
  160    WRITE(6,'(I4,3F12.7)')I,DXYZ(I*3-2),DXYZ(I*3-1),DXYZ(I*3)
         WRITE(6,'(/10X,''RESIDUAL ERROR'')')
         WRITE(6,'(4X,3F12.7)')SUMX,SUMY,SUMZ
         WRITE(6,*)
      ENDIF
      RETURN
      END
