C 11 APR 08 - JAB - HSSNUMX,RDVIB2: ALLOW NUMERICAL GRADIENT RESTARTS
C  4 MAR 08 - SPW - HESSX,VALFOR: CHANGES FOR NEO HESS WITH NUCOPT=.F.
C  4 MAR 08 - TJD - HSSANA: ENABLE SA-MCSCF GRADIENT PATHWAY
C  7 DEC 07 - SS  - TDDFT HESSIANS MUST BE SEMI-NUMERIC
C 24 MAR 07 - HM  - MCP HESSIAN CAN NOW BE SEMI-NUMERIC
C 22 DEC 06 - DGF - SUNCHRONISE PSILVL,MP2PAR
C  6 NOV 06 - MWS - ADJUST WAVEFUNCTION COMMON BLOCK
C  7 SEP 06 - TJD - HESSX: ALLOW ORMAS ANALYTIC HESSIAN
C  7 SEP 06 - MWS - HSSNUM,NUMGRDX: SAVE ENTIRE /ENRGYS/ COMMON
C 13 MAR 06 - TJD,NM - HSSANA: MCSCF FORCES GOPARR TRUE, EGIN/EGPUN ARGS
C  2 FEB 06 - MWS - NUMGRDX: SAVE VIB 0 ENERGY TERMS
C 17 JAN 06 - MWS - UPDATE ARGS TO TRFMCX
C 14 NOV 05 - DGF - PAD COMMON BLOCKS ENRGYS AND SCFOPT
C 19 SEP 05 - MWS - CHANGES FOR CPHF=AODDI, PRINT FULLY NUM. HESS COORDS
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 27 JUN 05 - MWS - HESSX: ALLOW ZAPT SEMINUMERICAL HESSIANS
C  1 JUN 05 - MWS - HESSX: MP2 SOLV., VAL. BOND HAS TO BE FULLY NUMERIC
C 14 MAR 05 - HL  - CHANGE FMASS IN FRGMSS COMMON
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - ENSURE THAT SVP MODEL USES FULLY NUMERIC HESSIAN
C  5 FEB 05 - MWS - PAD COMMON MCINP
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  7 APR 04 - RMO - NUMGRDX: ALLOW HIGHER ACCURACY DIFFERENTIATION
C 17 MAR 04 - MWS - HSSANA: CORRECT LOGIC OF AO INTEGRAL REGENERATION
C 22 DEC 03 - MWS - HESSX: TRAP NUMERIC GRADIENT ATTEMPT
C 12 DEC 03 - RMO - NUMGRDX,HSSFUL: ADDED FOR NUM. GRADIENTS + HESSIANS
C  9 DEC 03 - TJD - HSSANA: HOOKS FOR MCSCF ANALYTIC HESSIAN RUNS
C  9 DEC 03 - OQ  - HESSX: SAVE NORMAL COORD INFO FOR EXTENDED TDHF
C  3 SEP 03 - MWS - SET IHESSM BACK TO 0 WHEN HESSIAN IS FINISHED
C  3 JUL 03 - JMM - SUPPRESS PRINTING FOR MONTE CARLO JOBS
C 12 DEC 02 - MWS - HSSNUM: THERE ARE NO T INTEGRALS FOR SEMI-EMPIRICAL
C  7 AUG 02 - MWS - HESNUM: SAVE VIB 0 H,S,T FOR AIMPAC WITH MP2
C  7 AUG 02 - HL  - HESNUM: BUGFIX FOR PARTIAL HESSIANS
C 20 JUN 02 - MWS - HSSNUM: RESTORE VIB0 WAVEFUNC AT END, FIX DFT/EFP
C                   HSSANA: SKIP CPHF=AO INTEGRAL REGENERATION FOR C1
C 22 MAY 02 - MWS - PROVIDE ADDITIONAL TRFMCX/ZMATIN ARGUMENTS
C 17 APR 02 - MWS - SYNCH UP EFMULT AND FRGINF COMMON
C 26 MAR 02 - KRG - USE ABRT CALL
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 25 JUN 01 - MWS - VALFOR: DO VIB 0 DIPOLE EVEN IF NPRINT IS OFF
C 13 JUN 01 - DGF - HESSX: ALLOW FOR GRID DFT
C 20 FEB 01 - JAB - ADD PURADM TO PURIFY $ALPDR, REPLACE TFDDM WITH
C                   GENERIC TFTNSR, HESSX: CHECK FOR PRESENCE OF $ALPDR
C 29 DEC 00 - MWS - HSSANA: REMOVE PARALLEL AO CPHF ABORT
C 19 NOV 00 - HL  - HESSX,SYMFCM,FCMOUT,HSSNUM: PARTIAL HESSIAN CHANGES
C 11 OCT 00 - MWS - PAD COMMON BLOCK PCMPAR
C 28 JUL 00 - MWS - TOTAL CLEANUP OF NUMERICAL HESSIAN COMPUTATION CODE
C  1 MAY 00 - MWS - DDMOUT: ADD UNITS TO OUTPUT, HESSX: MOD FGMTRX ARGS
C 10 APR 00 - MWS - REMOVE STATIC MEMORY FROM COSMO COMMONS
C 25 MAR 00 - KKB,LNB,MWS - NUMERIC HESSIAN FORCED FOR COSMO, RESC, DFT
C 21 DEC 99 - MWS - REMOVE C1 DOWNSHIFT FROM TRANSF. BUT NOT AO INTS
C 25 OCT 99 - MWS - HSSNUM: SAVE VIB0 GRAD ONLY IF ATOMS EXIST
C 29 AUG 99 - CHC - VALFOR: QMMM MODIFICATION
C  6 JUN 99 - MWS - HSSNUM: SKIP MP2 AT VIB0 IF REGENERATING MOS
C  9 APR 99 - MWS - VALFOR: DO DIPOLES AFTER GRADIENT
C 13 MAR 99 - MWS - HSSNUM: SKIP SAVING/RESTORING Q MATRIX IF FRAG ONLY
C 13 FEB 99 - MWS - HSSNUM: RESTORE -Q- MATRIX AND ZERO -IVIB- AT END
C  1 DEC 98 - BMB - CHANGES TO ALLOW ANALYTIC HESSIANS WITH ECP'S
C 12 NOV 98 - GDF - HSSNUM,SYMDDM,SYMFCM: REMOVE ISOOUT USAGE
C 26 OCT 98 - MAF - HSSANA: ALLOW FOR USE OF SPHERICAL HARMONICS
C 12 SEP 98 - MWS - EGIN: EOF/ERR= BRANCHES ADDED.
C  6 MAY 98 - MWS - CHANGE -UPCASE- TO -UPRCAS-
C 27 FEB 98 - MWS - HSSANA: BOMB CONTINUUM SOLVENT RUNS
C  6 JAN 98 - JAB - HESSX:  ADD NSYMC TO COMMON /ZMAT/, DISABLE CALL
C                   TO PUZMAT WHEN SYM. COORDS IN USE.  PURFCM: CONVERT
C                   UNITS OF INTERNAL COORD. HESSIAN TO MDYN/ANG
C 20 NOV 97 - MWS - HSSNUM: FIX MCSCF RESTARTING FROM $VIB
C  1 OCT 97 - MWS - FCMOUT: USE 80 COLUMN OUTPUT, HESSX: SCFMI=NUMERIC
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C  8 AUG 97 - MWS - HESSX: PASS ADDITIONAL STORAGE TO FGMTRX FOR FRAGS
C 16 JUL 97 - GNM - EGMAC,HESSX,HSSNUM: CHANGES FOR FRAGONLY OPTION
C 21 APR 97 - MWS - HESSX: MP2 HESSIAN PERMITS ONLY METHOD=NUMERIC
C 17 JAN 97 - JA,PND - HESSX,FCMPUN: ALL NODES ENTER PUZMAT,ENUC
C  8 JAN 97 - MWS - HSSNUM: FIX EFP CHECK RUNS
C 18 DEC 96 - TLW - HSSANA: CHANGES FOR AO BASIS RESPONSES
C  9 DEC 96 - JHJ - DDINT: 2 NEW ARGUMENTS IN DIPDER CALL
C  3 DEC 96 - BMB - FCMIN: MOVE FNDGRP CALL BELOW DIMENSION SETUP
C 17 OCT 96 - SPW - VALFOR: ADD CALL TO CIGRAD FOR CI RUNS
C 10 SEP 96 - MWS - HESSX: CHECK INPUT SPELLING OF METHOD KEYWORDS
C  1 JUL 96 - MWS - EGIN: BROADCAST EF DATA TO SLAVES
C 13 JUN 96 - VAG - CHANGES TO INTRODUCE CITYP VARIABLE
C 25 MAY 96 - MWS - HSSNUM: DON'T USE SUBGROUPS IN MP2 AND MC HESSIANS
C 12 APR 96 - MWS - HSSNUM: FIX CHECK RUNS
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C  9 FEB 96 - WC  - HSSNUM: FIX FRAGMENT HSSEND RUNS AND $VIB RESTARTS
C 29 JAN 96 - WC  - HSSNUM: CHANGES IN COMPUTATION OF FRAGMENT HESSIAN,
C                   FCMIN: CHANGES IN READING OF HESSIANS WITH FRAGMENTS
C  3 JAN 96 - MWS - HSSNUM: CALL NEW ONE ELECTRON INTEGRAL DRIVER
C  6 DEC 95 - MWS - HESSX: PRINT UNITS USED FOR DDM
C 20 NOV 95 - MWS - HESSX,HSSANA: CHANGE PRINTOUT AMOUNTS
C 16 NOV 95 - PND - HSSNUM: FIX FRAG RESTART, VALFOR: PRINT FRAG COORDS
C 30 OCT 95 - MWS - HESSX,HSSNUM: CORRECTIONS FOR MP2 NUMERICAL HESSIAN
C 28 SEP 95 - MWS - HSSANA: EVALUATE CORRECT DIPOLE INTEGRALS
C 26 JUL 95 - MWS - HSSNUM: STOP GRACEFULLY IF TIME IS GETTING SHORT
C 31 MAY 95 - MWS - HESSX: SET IHESSM TO 1/2 FOR NUMERIC/ANALYTIC RUNS
C 24 MAY 95 - MWS - HSSANA: ALLOW USE OF DISTRIBUTED AO INTEGRALS
C 21 APR 95 - SPW - ADDED POLOUT FOR POLARIZABILITY OUTPUT
C 27 MAR 95 - SPW - HSSANA: CALL NEW DDINT FOR DIPOLE DERIVATIVE
C  4 JAN 95 - MWS - EGMAX: FIX FOR LINEAR CASES
C 29 DEC 94 - TLW - HESSX: ADD CHECK FOR EFIELD WHEN ANALYTIC HESSIANS
C 17 NOV 94 - WC  - FIX EFFECTIVE FRAGMENT NUM. HESS. RESTARTS
C 12 NOV 94 - MWS - REMOVE MOST FTNCHEK WARNINGS
C 27 OCT 94 - MWS - HESSX: DELETE UNUSED COMMON /DAIOLN/
C  6 OCT 94 - PND - HSSNUM: FLOATING POINT ARGUMENTS TO EFDSPL
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 19 JUL 94 - MWS - HESSX: BOMB JOBS WITH F OR G FUNCTIONS
C  1 JUN 94 - MWS - HSSANA: USE NEW INTEGRAL TRANSFORMATION
C 25 APR 94 - MWS - FFSVCI: ALL NODES MUST DO DISK OPERATIONS
C 21 APR 94 - MR  - RDVIB: DEAL WITH PARAGON FORMAT BUG
C  4 APR 94 - MWS - VALFOR: DELETE CALL TO EFPCGD
C 14 MAR 94 - MWS - HSSNUM: USE ONLY C1 FOR MCSCF DISPLACEMENTS
C 25 JAN 94 - MWS - HESSX,PURFCM: PUNCH E,GRAD,HESSIAN IF PRTIFC
C  4 JAN 94 - MWS - HSSANA: ABORT HESSIAN IF ENERGY DOES NOT CONVERGE
C  9 DEC 93 - MWS - PURDDM: CHANGE DAF RECORD NUMBERS
C 22 SEP 93 - MWS - FIX /OPTGRD/ TO MATCH STATPT CHANGES
C  6 AUG 93 - MWS - ADD ATORQ TO /FGRAD/
C  6 MAY 93 - PND - EGIN,EGPUN,HESSX: FIXED PUNCHING/READING GRADIENT
C                   WITH FRAGMENTS, FIXED HESS=CALC
C  2 APR 93 - TLW - HESSX: PURIFY GRADIENT
C  1 APR 93 - PND - EFFECTIVE FRAGMENT MODIFICATIONS.
C 29 MAR 93 - TLW - HESSX,PURDDM,PURFCM,TFDDM: RECTANGULAR B AND B INV
C 29 JAN 93 - JAB - EGPUN,HESSX: WORKAROUND FOR AIX'S IDAMAX BUG
C 28 JAN 93 - MR  - HSSNUM: SKIP VIB0 GEOMETRY DURING HSSEND RUNS
C 11 DEC 92 - MWS - HSSANA: FIX EVALUATION/SKIP WAVEFUNCTION TEST
C 22 NOV 92 - MWS - HESSX: INITIALIZE MESSAGE TYPE IDENTIFIERS
C 23 OCT 92 - MWS - PERMIT CONTROL OVER NUMERICAL HESSIAN PRINT/PUNCH
C 16 SEP 92 - TLW - HSSANA,HESSX: PARALLEL ANALYTIC HESSIANS
C 28 AUG 92 - MWS - FFSVCI: SAVE MORE THAN ONE STATE'S VECTORS
C 22 JUN 92 - MWS - HSSNUM: FIX NUMERICAL RESTARTS
C 15 JUN 92 - MWS - FORCE DISK BASED SCF FOR ANY ANALYTIC HESSIAN RUN
C  5 APR 92 - MWS - USE MSGSET FOR MESSAGE MASKS
C 27 MAR 92 - TLW - RDVIB: PARALLELIZED
C 17 MAR 92 - MWS - RDVIB: MIXED CASE END OF GROUP ALLOWED
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  9 MAR 92 - JHJ - HESSX,VALFOR:IF MPC WFN,NUMERIC HESS,NO DIPOLE INT.
C 24 FEB 92 - MWS - CHANGE TO USE OF COMMON SYMMOL
C 20 FEB 92 - TLW - DDMIN,EGIN,FCMIN,HSSNUM: USE FNDGRP AND ADD BRDCASTS
C 11 JAN 92 - TLW - MAKE READS PARALLEL
C 10 JAN 92 - TLW - CHANGE REWINDS TO CALL SEQREW
C 10 JAN 92 - MWS,TLW - CHANGE OPENCF,OPENDF,OPENPK,OPNIRC TO SEQOPN
C  6 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
C  2 DEC 91 - JHJ - VALFOR:ADD CALL TO EFPCGD.
C 28 NOV 91 - STE - EGIN: USE MORE CHARACTER VARIABLES
C 24 OCT 91 - JHJ - HESSX:NUMERIC HESSIAN FOR EFP/ZRF.
C 14 AUG 91 - TLW - HSSNUM:,SYMDDM:, AND SYMFCM: ADDED FTR AND GTR TO
C                   COMMON SYMSPD
C  6 AUG 91 - MWS - FIX RECYCLING OF VIB 0 VECTOR WHEN READING $VIB
C 26 JUN 91 - MWS - SAVE VIB 0 CI VECTOR FOR NUMERICAL RUNS.
C  9 APR 91 - MR  - HSSNUM: RETAIN SOME SYMMETRY IN NUMERICAL HESSIAN,
C                   RESTARTS USE CONVERGED SCF MOS, NOT HUCKEL GUESS.
C  1 MAR 91 - MWS - CALLING ARG TO HESSX TO FORCE VIBRATIONAL ANALYSIS
C 20 FEB 91 - JAM - FIX CALL TO ZMATIN BY ADDING ARGUMENT
C  7 FEB 91 - MWS - RECOVER JAN'S THERMOCHEMISTRY INPUT OPTIONS
C  5 FEB 91 - JAB - ADD IHREP, FINLHS OPTIONS TO $STATPT.
C  5 OCT 90 - MWS - HSSANA: TRANSFORM BEFORE DOING 2ND DER. INTS
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C 22 JUL 90 - MWS - ALLOW ANALYTIC COMPUTATION OF OS-TCSCF HESSIANS
C 20 FEB 90 - MR  - UPDATE FORMAT 9040 IN HSSNUM
C  9 DEC 89 - MWS - USE CORRECT NUCLEAR CHARGE FOR ECP CASES,
C                   EXETYP=CHECK WITH $VIB USED CORRECTED
C 29 NOV 89 - MWS - OPEN SHELL SCF USES ANALYTIC HESSIAN BY DEFAULT
C 22 SEP 89 - MWS - CHANGE ANGSTROMS-BOHR CONVERSION FACTOR
C 21 SEP 89 - STE - REGULARIZE /INTFIL/
C  3 MAY 89 - MWS - HSSANA: PASS SCFTYP TO CPHFX, HSSNUM: SAVE E0
C 28 APR 89 - JAB - CALL PURFCM IF DECOMP=.T.
C 28 FEB 89 - MWS - PRINT INTERNAL COORD GRADIENT IF PRTIFC CHOSEN
C 26 FEB 89 - STE - DELETE FORMATS 9050,9060 IN VALFOR
C 19 FEB 89 - MWS - DEFAULT PURIFY=F, ALWAYS PUNCH $HESS,$GRAD,$DDM
C 13 FEB 89 - MWS - SAVE 3N-6 HESS ONLY IN PURFCM, UNLESS DECOMP RUN
C 30 JAN 89 - MWS - REMOVE ZMASS STORAGE FROM FGMTRX CALL
C 20 JAN 89 - MWS - NEW MODULE TO DRIVE HESSIAN COMPUTATION
C
C*MODULE HESS    *DECK DDINT
      SUBROUTINE DDINT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL OUT,DBG
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA DEBUG/8HDEBUG   /, DBUGME/8HDIPDER  /, GRD1/8HGRD1    /
      DATA DIPDM/8HDIPDM   /
C
      DBG = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME)  .AND.  MASWRK
      OUT = (EXETYP.EQ.DEBUG  .OR.  EXETYP.EQ.DBUGME  .OR.
     *       EXETYP.EQ.GRD1   .OR.  NPRINT.EQ.-3)
      L1 = NUM
      L2 =(L1*L1+L1)/2
C
      CALL VALFM(LOADFM)
      LDA = LOADFM+1
      LDB = LDA + L2
      LDDV= LDB + L2
      LAST= LDDV+ NAT*9
      NEED = LAST-LOADFM
      CALL GETFM(NEED)
C
C     ---- GET TOTAL DENSITY ----
C
      CALL DENDD1(X(LDA),X(LDB),L2)
C
C     ---- CALCULATE DIPOLE INTEGRAL DERIVATIVES ----
C
      DUMMY = 0.0D+00
      CALL DIPDER(X(LDA),DUMMY,DUMMY,X(LDDV),L2,DUMMY,0,0,0,DIPDM)
      CALL DAWRIT(IDAF,IODA,X(LDDV),9*NAT,34,0)
      IF(OUT.OR.DBG) THEN
          WRITE(IW,9040)
          CALL DDMOUT(X(LDDV),3*NAT)
      END IF
C
      CALL RETFM(NEED)
C
      RETURN
 9040 FORMAT(/10X,42(1H-)/
     *        10X,'DIPOLE INTEGRAL DERIVATIVE CONTRIBUTION TO'/
     *        10X,42(1H-))
      END
C*MODULE HESS    *DECK DDMIN
      SUBROUTINE DDMIN(DDM,NCOORD,GOTDDM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOTDDM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION DDM(3,NCOORD)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- READ DIPOLE DERIVATIVE MATRIX FROM CARDS -----
C
      GOTDDM = .FALSE.
      CALL VCLR(DDM,1,3*NCOORD)
C
C           POSITION INPUT FILE TO $DIPDR
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $DIPDR ',IEOF)
      IF (IEOF.EQ.1) GO TO 400
C
      GOTDDM=.TRUE.
      IF (MASWRK) THEN
         DO 240 J=1,NCOORD
            READ(IR,9010) (DDM(I,J),I=1,3)
  240    CONTINUE
      END IF
C
C SEND DDM TO ALL SLAVES
C
      IF (GOPARR) CALL DDI_BCAST(450,'F',DDM,3*NCOORD,MASTER)
C
  400 CONTINUE
      CALL DAWRIT(IDAF,IODA,DDM,3*NCOORD,34,0)
      RETURN
C
 9010 FORMAT(1X,5F15.8)
      END
C*MODULE HESS    *DECK DDMOUT
      SUBROUTINE DDMOUT(DDM,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION DDM(3,NCOORD)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT DIPOLE DERIVATIVE MATRIX -----
C
      IF (MASWRK) THEN
      WRITE(IW,9040)
      DO 480 J=1,NAT
         WRITE(IW,9060) ANAM(J),BNAM(J),(DDM(I,3*J-2),I=1,3)
         WRITE(IW,9080)                 (DDM(I,3*J-1),I=1,3)
         WRITE(IW,9100)                 (DDM(I,3*J  ),I=1,3)
  480 CONTINUE
      END IF
      RETURN
C
 9040 FORMAT(/10X,24(1H-),10X,16(1H-)/
     *        10X,'DIPOLE DERIVATIVE TENSOR',10X,'(DEBYE/ANGSTROM)'/
     *        10X,24(1H-),10X,16(1H-)//
     *        8X,'ATOM',17X,'MU-X',11X,'MU-Y',11X,'MU-Z')
 9060 FORMAT(5X,A8,A2,' D/DX ',3F15.9)
 9080 FORMAT(15X,     ' D/DY ',3F15.9)
 9100 FORMAT(15X,     ' D/DZ ',3F15.9)
      END
C*MODULE HESS    *DECK DDMPUN
      SUBROUTINE DDMPUN(DDM,NCOORD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION DDM(3,NCOORD)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C      ----- PUNCH THE DIPOLE DERIVATIVE MATRIX -----
C
      IF (MASWRK) THEN
      WRITE(IP,FMT='(7H $DIPDR)')
      DO 10 J=1,NCOORD
         WRITE(IP,9002)  DDM(1,J),DDM(2,J),DDM(3,J)
   10 CONTINUE
      WRITE(IP,FMT='(5H $END)')
      END IF
      RETURN
C
 9002 FORMAT(1X,1P,3E15.8)
      END
C*MODULE HESS    *DECK EGIN
      SUBROUTINE EGIN(EG,NCOORD,GOTEG,GRDWRD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      CHARACTER*8 GRDWRD
C
      LOGICAL GOTEG,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION EG(NCOORD)
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FUNCT / E,EGRAD(3,MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- READ GRADIENT VECTOR FROM CARDS - $GRAD GROUP -----
C     IF NOT PRESENT IN INPUT, A ZERO GRADIENT VECTOR IS GENERATED.
C
      CALL VCLR(EG,1,NCOORD)
      IF(NFRG.NE.0) THEN
         CALL VCLR(DEFT,1,3*NFRG)
         CALL VCLR(TORQ,1,3*NFRG)
      END IF
      GOTEG=.FALSE.
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,GRDWRD,IEOF)
      IF (IEOF.EQ.1) GO TO 300
C
      GOTEG=.TRUE.
      IF (MASWRK) THEN
         READ(IR,9030,END=700,ERR=700) E
         READ(IR,9040,END=700,ERR=700) EG
         IF(NFRG.NE.0) THEN
            DO 150 I=1,NFRG
               READ(IR,9040,END=700,ERR=700) (DEFT(J,I),J=1,3)
               READ(IR,9040,END=700,ERR=700) (TORQ(J,I),J=1,3)
 150        CONTINUE
         END IF
         ISTAT=0
      END IF
C
C        SEND ENERGY AND ENERGY GRADIENT TO ALL SLAVES
C
      IF (GOPARR) THEN
         CALL DDI_BCAST(400,'I',ISTAT,1  ,MASTER)
         IF(ISTAT.NE.0) CALL ABRT
         CALL DDI_BCAST(401,'F',E ,1     ,MASTER)
         CALL DDI_BCAST(402,'F',EG,NCOORD,MASTER)
         IF(NFRG.NE.0) THEN
            CALL DDI_BCAST(403,'F',DEFT,3*NFRG,MASTER)
            CALL DDI_BCAST(404,'F',TORQ,3*NFRG,MASTER)
         END IF
      END IF
C
  300 CONTINUE
      CALL DAWRIT(IDAF,IODA,EG,NCOORD, 3,0)
      RETURN
C
C         MASTER HIT SOME KIND OF ERROR READING $GRAD
C
  700 CONTINUE
      WRITE(IW,9060)
      ISTAT=1
      CALL DDI_BCAST(400,'I',ISTAT,1,MASTER)
      CALL ABRT
      STOP
C
 9030 FORMAT(2X,F20.10)
 9040 FORMAT(15X,3F20.10)
 9060 FORMAT(1X,'**** ERROR READING $GRAD GROUP *****')
      END
C*MODULE HESS    *DECK EGMAX
      SUBROUTINE EGMAX(EG,NCOORD,DEFT,TORQ,NFRG,GMAX,GRMS)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (ZERO = 0.0D+00)
C
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
C
      DIMENSION EG(NCOORD),DEFT(*),TORQ(*)
C
      IF (ICOORD .EQ.4) THEN
         GMAX = ZERO
         GRMS = ZERO
      ELSE
         LOC = IDAMAX(NCOORD,EG,1)
         GMAX = ABS(EG(LOC))
         GRMS = DDOT(NCOORD,EG,1,EG,1)
      END IF
      IF(NFRG.NE.0) THEN
         NCF2  = NFRG*3
         LOC   = IDAMAX(NCF2,DEFT,1)
         GMAX1 = ABS(DEFT(LOC))
         LOC   = IDAMAX(NCF2,TORQ,1)
         GMAX2 = ABS(TORQ(LOC))
         GMAX  = MAX(GMAX,GMAX1,GMAX2)
         GRMS1 = DDOT(NCF2,DEFT,1,DEFT,1)
         GRMS2 = DDOT(NCF2,TORQ,1,TORQ,1)
         GRMS = GRMS + GRMS1 + GRMS2
      END IF
      NCTOT = NCOORD + 6*NFRG
      GRMS  = SQRT(GRMS/NCTOT)
      RETURN
      END
C*MODULE HESS    *DECK EGOUT
      SUBROUTINE EGOUT(EG,NAT)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION EG(3,NAT)
      DIMENSION DNAM(3),EGT(3),ETT(3)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /INFOA / NDUM,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
      DATA DNAM /4HE'X ,4HE'Y ,4HE'Z /
C
C     ----- PRINT ENERGY GRADIENT VECTOR -----
C
      IF(NFRG.EQ.0) GO TO 200
C
      DO 50 J=1,3
         EGT(J)=0.0D+00
         ETT(J)=0.0D+00
   50 CONTINUE
C
      DO 90 I=1,NAT
         ETT(1)=ETT(1)+EG(3,I)*(C(2,I)-EFCENT(2,1))
     *                -EG(2,I)*(C(3,I)-EFCENT(3,1))
         ETT(2)=ETT(2)+EG(1,I)*(C(3,I)-EFCENT(3,1))
     *                -EG(3,I)*(C(1,I)-EFCENT(1,1))
         ETT(3)=ETT(3)+EG(2,I)*(C(1,I)-EFCENT(1,1))
     *                -EG(1,I)*(C(2,I)-EFCENT(2,1))
         DO 80 J=1,3
            EGT(J)=EGT(J)+EG(J,I)
   80    CONTINUE
   90 CONTINUE
C
  200 CONTINUE
      IF (MASWRK) THEN
         WRITE(IW,9000) (DNAM(J),J=1,3)
         DO 220 I=1,NAT
            WRITE(IW,9010) I,ANAM(I),BNAM(I),(EG(J,I),J=1,3)
  220    CONTINUE
         IF(NFRG.GT.0) THEN
            WRITE(IW,9020) (EGT(J),J=1,3)
            WRITE(IW,9023) (ETT(J),J=1,3)
         END IF
      END IF
      RETURN
C
 9000 FORMAT(/6X,'ATOM',17X,A4,14X,A4,14X,A4)
 9010 FORMAT(1X,I4,1X,A8,A2,3F18.9)
 9020 FORMAT(6X,'TOTALS',4X,3F18.9)
 9023 FORMAT(6X,'TORQUE',4X,3F18.9)
      END
C*MODULE HESS    *DECK EGPUN
      SUBROUTINE EGPUN(EG,NAT,GWRD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 GWRD
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION EG(3,NAT)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FUNCT / E,EGRAD(3,MXATM)
      COMMON /INFOA / NATM,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
C
C     ----- PUNCH GRADIENT VECTOR -----
C
      CALL EGMAX(EG,3*NAT,DEFT,TORQ,NFRG,GMAX,GRMS)
C
      IF (MASWRK  .AND.  NPRTGO.GT.0) THEN
         WRITE(IP,9040) GWRD
         WRITE(IP,9010) E,GMAX,GRMS
         DO 100 I=1,NAT
            ZNUC = ZAN(I) + IZCORE(I)
            WRITE(IP,9000) ANAM(I),BNAM(I),ZNUC,(EG(J,I),J=1,3)
  100    CONTINUE
         IF(NFRG.NE.0) THEN
            DO 150 I=1,NFRG
               WRITE(IP,9020) I,(DEFT(J,I),J=1,3)
               WRITE(IP,9030) I,(TORQ(J,I),J=1,3)
 150        CONTINUE
         END IF
         WRITE(IP,FMT='(5H $END)')
      END IF
      RETURN
 9000 FORMAT(A8,A2,F5.0,1P,3E20.10)
 9010 FORMAT('E=',F20.10,'  GMAX=',F12.7,'  GRMS=',F12.7)
 9020 FORMAT('FRAG',I2,' TRAN',4X,1P,3E20.10)
 9030 FORMAT('FRAG',I2,' TORQ',4X,1P,3E20.10)
 9040 FORMAT(A8)
      END
C*MODULE HESS    *DECK FCMIN
      SUBROUTINE FCMIN(FCM,NC1,GOTEH)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXPT=2000, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL GOTEH,GOPARR,DSKWRK,MASWRK,EXPAND,CONSTR
C
      DIMENSION FCM(*)
      DIMENSION TEXT1(10)
      CHARACTER*8 FRGNME
C
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DATA SMALL/1.0D-02/
C
C     ----- READ HESSIAN MATRIX FROM CARDS -----
C
C     SET UP FOR READING AB PART ONLY OR ALL AB COOR IN THE CASE OF
C     FRAGMENT OPT
C
      EXPAND=.FALSE.
      CONSTR=.FALSE.
      NAT3=3*NAT
      NFRG6=6*NFRG
      NCT=NC1
      NC2=NC1
      ISTART=1
      IF (NFRG.NE.0.AND.NC1.EQ.NAT3) EXPAND=.TRUE.
      IF (NFRG.NE.0.AND.NC1.GT.(NAT3+NFRG6)) CONSTR=.TRUE.
      IF (EXPAND.OR.CONSTR) NCT=NAT3+NFRG6
      IF (EXPAND) NC2=NCT
      NCT2=NCT**2
      IF (CONSTR) ISTART=NCT2+1
      IF (CONSTR) CALL EFCM
C
      CALL VCLR(FCM,1,NCT2)
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $HESS  ',IEOF)
      IF (IEOF.EQ.1) THEN
         GOTEH=.FALSE.
         GO TO 800
      ELSE
         GOTEH=.TRUE.
      END IF
C
      IF (MASWRK) THEN
         READ (IR,9000) TEXT1
         WRITE(IW,9010) TEXT1
C
         II = 0
         ICC = 0
         DO 240 I=1,NC1
            IC=0
            DO 220 MINCOL = 1, NC1, 5
               MAXCOL=MINCOL+4
               IF(MAXCOL.GT.NC1) MAXCOL=NC1
               IC=IC+1
               READ(IR,9020) II,ICC,
     *           (FCM(ISTART+(J-1)*NC2),J=MINCOL,MAXCOL)
               IMD100 = MOD(I,100)
               IF (II .NE. IMD100 .AND. ICC .NE. IC) THEN
                  WRITE(IW,9030) I,IC,II,ICC
                  CALL ABRT
               END IF
  220       CONTINUE
            ISTART=ISTART+1
  240    CONTINUE
C
        IF (EXPAND) THEN
          HESSII1=0.01D+00
          HESSII2=0.03D+00
C  FRAG TRANS
          DO 350 I=1,NFRG
            DO 330 J=1,3
              FCM(ISTART*(NC2+1)-NC2)=HESSII1
              ISTART=ISTART+1
  330       CONTINUE
C  FRAG ROT
            DO 340 J=1,3
              FCM(ISTART*(NC2+1)-NC2)=HESSII2
              ISTART=ISTART+1
  340       CONTINUE
  350     CONTINUE
        END IF
C
        IF (CONSTR) THEN
C   AB-AB
          DO 410 I=1,NAT3
            INDX1=NCT2+(I-1)*NC2+1
            INDX2=(I-1)*NCT+1
            CALL DCOPY(NAT3,FCM(INDX1),1,FCM(INDX2),1)
  410     CONTINUE
C   EF-AB AND AB-EF
          DO 460 I=1,NAT3
            INDX4=NCT2+(I-1)*NC2+NAT3
            INDX5=0
            INDX6=NAT3*NCT+I
            KMASS=0
            DO 450 J=1,NFRG
              INDX3=(I-1)*NCT+NAT3+(J-1)*6
C   EF-AB
              DO 430 K=1,NMPTS(J)
                INDX5=INDX5+1
                IF (ABS(FMASS(K+KMASS)).LE.SMALL) GO TO 430
                DO 420 L=1,3
                  M=MOD(L+1,4)+L/3
                  N=MOD(L+2,4)+(L+1)/3
                  FCM(INDX3+L)=FCM(INDX3+L)+FCM(INDX4+L)
                  FCM(INDX3+L+3)=FCM(INDX3+L+3)
     *                     -FCM(INDX4+M)*(EFC(N,INDX5)-EFCENT(N,J))
     *                     +FCM(INDX4+N)*(EFC(M,INDX5)-EFCENT(M,J))
  420           CONTINUE
                INDX4=INDX4+3
  430         CONTINUE
C   AB-EF
              DO 440 K=1,6
                FCM(INDX6)=FCM(INDX3+K)
                INDX6=INDX6+NCT
  440         CONTINUE
              KMASS=KMASS+NMPTS(J)
  450       CONTINUE
  460     CONTINUE
C   EF-EF
          INDX1=NCT2+(NC2+1)*NAT3
          INDX2=(NCT+1)*NAT3
          INDX3=(NCT+1)*NAT3+3*NCT
          JC=0
          KMASS=0
          DO 590 I=1,NFRG
            DO 580 J=1,NMPTS(I)
              JC=JC+1
              IF (ABS(FMASS(J+KMASS)).LE.SMALL) GO TO 580
              DO 570 K=1,3
                K1=MOD(K+1,4)+K/3
                K2=MOD(K+2,4)+(K+1)/3
                INDX4=INDX1+(K-1)*NC2
                INDX5=INDX1+(K1-1)*NC2
                INDX6=INDX1+(K2-1)*NC2
                INDX7=INDX2+(K-1)*NCT
                INDX8=INDX3+(K-1)*NCT
                CK1=EFC(K1,JC)-EFCENT(K1,I)
                CK2=EFC(K2,JC)-EFCENT(K2,I)
                MC=0
                KMASS2=0
                DO 560 L=1,NFRG
                  DO 550 M=1,NMPTS(L)
                    MC=MC+1
                    IF (ABS(FMASS(M+KMASS2)).LE.SMALL) GO TO 550
                    DO 540 N=1,3
                     N1=MOD(N+1,4)+N/3
                     N2=MOD(N+2,4)+(N+1)/3
                     FCM(INDX7+N)=FCM(INDX7+N)+FCM(INDX4+N)
                     FCM(INDX7+N+3)=FCM(INDX7+N+3)
     *                -FCM(INDX4+N1)*(EFC(N2,MC)-EFCENT(N2,L))
     *                +FCM(INDX4+N2)*(EFC(N1,MC)-EFCENT(N1,L))
                     FCM(INDX8+N)=FCM(INDX8+N)
     *                -FCM(INDX5+N)*CK2+FCM(INDX6+N)*CK1
                     FCM(INDX8+N+3)=FCM(INDX8+N+3)
     *                +FCM(INDX5+N1)*CK2*(EFC(N2,MC)-EFCENT(N2,L))
     *                -FCM(INDX5+N2)*CK2*(EFC(N1,MC)-EFCENT(N1,L))
     *                -FCM(INDX6+N1)*CK1*(EFC(N2,MC)-EFCENT(N2,L))
     *                +FCM(INDX6+N2)*CK1*(EFC(N1,MC)-EFCENT(N1,L))
  540               CONTINUE
                    INDX4=INDX4+3
                    INDX5=INDX5+3
                    INDX6=INDX6+3
  550             CONTINUE
                  INDX7=INDX7+6
                  INDX8=INDX8+6
                  KMASS2=KMASS2+NMPTS(L)
  560           CONTINUE
  570         CONTINUE
              INDX1=INDX1+NC2*3
  580       CONTINUE
            INDX2=INDX2+6*NCT
            INDX3=INDX3+6*NCT
            KMASS=KMASS+NMPTS(I)
  590     CONTINUE
C
        END IF
C
      END IF
C
C SEND FCM TO ALL SLAVES
C
      IF (GOPARR) THEN
         CALL DDI_BCAST(425,'F',FCM,NCT2,MASTER)
      END IF
C
  800 CONTINUE
      CALL DAWRIT(IDAF,IODA,FCM,NCT2,4,0)
      IF (CONSTR) CALL FCMPUN(FCM(1),NCT)
      RETURN
C
 9000 FORMAT(10A8)
 9010 FORMAT(/1X,'$HESS GROUP READ FROM CARDS'/1X,10A8)
 9020 FORMAT(I2,I3,5F15.8)
 9030 FORMAT(' *** ERROR READING FORCE CONSTANT MATRIX ELEMENT',2I4/
     *       '           INPUT WAS',2I4/
     *  'DO YOU NEED TO INCLUDE A TITLE CARD BELOW THE $HESS CARD?'/
     *  'A BLANK TITLING CARD IS PERFECTLY OK.')
      END
C*MODULE HESS    *DECK FCMOUT
      SUBROUTINE FCMOUT(FCM,NCOORD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION FCM(NCOORD,NCOORD)
      DIMENSION CLAB(3)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FROZNC/ IFREEZ(3*MXATM),NFRZ
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CLAB /4H   X,4H   Y,4H   Z/
      DATA HSSIAN/8HHESSIAN /
C
      IF(.NOT.MASWRK) RETURN
C
      WRITE(IW,9000)
C
      IF(IFREEZ(1).NE.0.AND.RUNTYP.EQ.HSSIAN.AND.MASWRK)
     *  WRITE (IW, 9005)
C
      INCR= 2
      IF(NPRINT.EQ.6) INCR=4
      DO 140 MINCOL=1,NAT,INCR
         MAXCOL = MINCOL+INCR-1
         IF (MAXCOL .GT. NAT) MAXCOL = NAT
         WRITE(IW,9010) (N,N = MINCOL,MAXCOL)
         WRITE(IW,9020) (ANAM(N),BNAM(N),N = MINCOL,MAXCOL)
         WRITE(IW,9030) ((CLAB(M),M = 1,3),N = MINCOL,MAXCOL)
         J0 = 3*(MINCOL-1)+1
         J1 = 3*MAXCOL
         IFC=1
         DO 120 IAT = MINCOL,NAT+2*NFRG
            I0 = 3*(IAT-1)
            IF(IAT.LE.NAT) THEN
               WRITE(IW,9040) IAT,ANAM(IAT),BNAM(IAT),CLAB(1),
     *                         (FCM(I0+1,J),J = J0,J1)
            ELSE
               IFCT=MOD(IAT-NAT,2)
               IF(IFCT.EQ.1) THEN
                 WRITE(IW,9050) IFC,CLAB(1),
     *                         (FCM(I0+1,J),J = J0,J1)
               ELSE IF(IFCT.EQ.0) THEN
                 WRITE(IW,9060) IFC,CLAB(1),
     *                         (FCM(I0+1,J),J = J0,J1)
               IFC=IFC+1
               END IF
            END IF
            WRITE(IW,9070) CLAB(2),(FCM(I0+2,J),J = J0,J1)
            WRITE(IW,9070) CLAB(3),(FCM(I0+3,J),J = J0,J1)
  120    CONTINUE
  140 CONTINUE
C                    TWO SETS OF COLUMNS (TRANSLATION AND ROTATION)
C                    FOR PRINTING THE FRAGMENT/FRAGMENT BLOCK
      DO 150 IFRG=1,NFRG
         WRITE(IW,9080) IFRG
         WRITE(IW,9030) ((CLAB(M),M = 1,3),N = 1,2)
         J0 = 3*NAT + 6*(IFRG-1) + 1
         J1 = J0 + 5
         DO 130 JFRG = IFRG,NFRG
            I0 = 3*NAT + 6*(JFRG-1)
            WRITE(IW,9050) JFRG,CLAB(1),(FCM(I0+1,J),J = J0,J1)
            WRITE(IW,9070)      CLAB(2),(FCM(I0+2,J),J = J0,J1)
            WRITE(IW,9070)      CLAB(3),(FCM(I0+3,J),J = J0,J1)
            WRITE(IW,9060) JFRG,CLAB(1),(FCM(I0+4,J),J = J0,J1)
            WRITE(IW,9070)      CLAB(2),(FCM(I0+5,J),J = J0,J1)
            WRITE(IW,9070)      CLAB(3),(FCM(I0+6,J),J = J0,J1)
  130    CONTINUE
  150 CONTINUE
      RETURN
C
 9000 FORMAT(/10X,31(1H-)/10X,'CARTESIAN FORCE CONSTANT MATRIX'/
     *        10X,31(1H-))
 9005 FORMAT(//1X,'THERE ARE FROZEN COORDINATES, ',
     *'PARTIAL HESSIAN ANALYSIS WILL BE DONE.'//1X,
     *'FOR THE FROZEN ATOMS AND FRAGMENTS, ',
     *'DIAGONAL MATRIX ELEMENTS ARE '/1X,'SET TO BE 1.0D-08,',
     *' ALL THE OTHER MATRIX ELEMENTS RELATED TO FROZEN '/1X,
     *'ATOMS AND FRAGMENTS ARE SET TO BE ZERO.'//)
 9010 FORMAT(/20X,4(12X,I4,11X))
 9020 FORMAT( 20X,4(11X,A8,A2,6X))
 9030 FORMAT( 20X,4(3X,A4,5X,A4,5X,A4,2X))
 9040 FORMAT(I3,3X,A8,A2,A4,12F9.6)
 9050 FORMAT(1X,'FRAG.',I3,1X,'TRANS.',A4,12F9.6)
 9060 FORMAT(1X,'FRAG.',I3,1X,'ROT.  ',A4,12F9.6)
 9070 FORMAT(16X,A4,12F9.6)
 9080 FORMAT(/20X,'FRAGMENT',I3,3X,'TRANSLATION',13X,'ROTATION ')
      END
C*MODULE HESS    *DECK FCMPUN
      SUBROUTINE FCMPUN(FCM,NCOORD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION FCM(NCOORD,NCOORD)
C
      PARAMETER (MXATM=2000, MXRT=100)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SIMDAT/ NACC,NREJ,IGOMIN,NRPA,IBWM,NACCT,NREJT,NRPAT,
     *                NPRTGO,IDPUNC,IGOFLG
C
C     ----- PUNCH HESSIAN -----
C
      IF (MASWRK  .AND.  NPRTGO.GT.0) THEN
         WRITE(IP,FMT='(6H $HESS)')
         WRITE(IP,8000) E,ENUCR
         CALL PUSQL(FCM,NCOORD,NCOORD,NCOORD)
         WRITE(IP,FMT='(5H $END)')
      END IF
      RETURN
 8000 FORMAT('ENERGY IS',F20.10,' E(NUC) IS',F20.10,5X)
      END
C*MODULE HESS    *DECK FFSVCI
      SUBROUTINE FFSVCI(IVIB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SVDSKW
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /FMCOM / X(1)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        --- PRESERVE/RESTORE THE VIB 0 CI EIGENVECTORS ---
C
      SVDSKW = DSKWRK
      DSKWRK = .FALSE.
C
      NFT17=17
C
      IF(IVIB.EQ.0) THEN
         NFTIN = NFT12
         NFTOUT= NFT17
         CALL SEQOPN(NFTOUT,'CSFSAVE','UNKNOWN',.FALSE.,'UNFORMATTED')
      ELSE
         NFTIN = NFT17
         NFTOUT= NFT12
      END IF
C
      CALL SEQREW(NFTIN)
      IF(MASWRK) READ(NFTIN) NSTATE,NWKS
      CALL DDI_BCAST(1232,'I',NSTATE,1,MASTER)
      CALL DDI_BCAST(1233,'I',NWKS  ,1,MASTER)
C
      CALL VALFM(LOADFM)
      LCI  = LOADFM + 1
      LAST = LCI    + NWKS
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      CALL SEQREW(NFTOUT)
      IF(MASWRK) WRITE(NFTOUT) NSTATE,NWKS
C
      DO I=1,NSTATE
         CALL SQREAD(NFTIN ,X(LCI),NWKS)
         CALL SQWRIT(NFTOUT,X(LCI),NWKS)
      ENDDO
C
      CALL SEQREW(NFTIN)
      CALL SEQREW(NFTOUT)
C
      CALL RETFM(NEED)
      DSKWRK = SVDSKW
      RETURN
      END
C*MODULE HESS    *DECK FFSVMO
      SUBROUTINE FFSVMO(IVIB,WRK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION WRK(*)
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PSILVL/ IPSI,ISKPRP
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA RHF,RMC /8HRHF     ,8HMCSCF   /
      DATA RNONE/8HNONE    /
C
      L2 = (NUM*NUM+NUM)/2
      L3 = NUM*NUM
C
C     ----- WRITE VIB 0 MO'S AND DENSITY ON SPECIAL DAF FILES -----
C
      IF (IVIB .GT. 0) GO TO 100
      CALL DAREAD(IDAF,IODA,WRK,L3,15,0)
      CALL DAWRIT(IDAF,IODA,WRK,L3,30,0)
      CALL DAREAD(IDAF,IODA,WRK,L2,16,0)
      CALL DAWRIT(IDAF,IODA,WRK,L2,32,0)
      IF(SCFTYP.EQ.RHF) RETURN
      IF(SCFTYP.EQ.RMC) RETURN
      IF(CITYP.NE.RNONE  .AND.  IPSI.EQ.1) RETURN
      CALL DAREAD(IDAF,IODA,WRK,L3,19,0)
      CALL DAWRIT(IDAF,IODA,WRK,L3,31,0)
      CALL DAREAD(IDAF,IODA,WRK,L2,20,0)
      CALL DAWRIT(IDAF,IODA,WRK,L2,33,0)
      RETURN
C
C     ----- RECOVER VIB 0 MO'S AND DENSITY -----
C
  100 CONTINUE
      CALL DAREAD(IDAF,IODA,WRK,L3,30,0)
      CALL DAWRIT(IDAF,IODA,WRK,L3,15,0)
      CALL DAREAD(IDAF,IODA,WRK,L2,32,0)
      CALL DAWRIT(IDAF,IODA,WRK,L2,16,0)
      IF(SCFTYP.EQ.RHF) RETURN
      IF(SCFTYP.EQ.RMC) RETURN
      IF(CITYP.NE.RNONE  .AND.  IPSI.EQ.1) RETURN
      CALL DAREAD(IDAF,IODA,WRK,L3,31,0)
      CALL DAWRIT(IDAF,IODA,WRK,L3,19,0)
      CALL DAREAD(IDAF,IODA,WRK,L2,33,0)
      CALL DAWRIT(IDAF,IODA,WRK,L2,20,0)
      RETURN
      END
C*MODULE HESS    *DECK HESSX
      SUBROUTINE HESSX(VIBOVR,GOTWFN)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION METHOD,METHNR
C
      LOGICAL PRTIFC,PURIFY,PROJCT,VIBANL,RDHESS,DECOMP,LINEAR,
     *        GOTEG,GOTEH,GOTDDM,GOTADM,STATPT,SOME,PRTSCN,GOTFRQ,
     *        VIBOVR,GOTWFN,GOPARR,DSKWRK,MASWRK,TEST,
     *        PULCOR,EFLDL,POLAR,HSSEND,SCFOK,REDOVB,OK,FGONLY,
     *        QMMM,MMONLY,CARTSN,CANONC,FCORE,FORS,NOCI,EKT,LINSER,
     *        MPGRD,CIGRD,DIISON
C
      PARAMETER (MXATM=2000, MXAO=8192, MXPT=2000, MXFRG=1050,
     *           MXFGPT=12000, MXNORO=250, MAXNFRAGS=10,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (NNAM=18)
C
      DIMENSION QNAM(NNAM), KQNAM(NNAM), APOL(6), RAMALF(6)
C
      COMMON /CHFINF/ CHFSLV,CPTYPE,NWDCHF,POLAR
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /FROZNC/ IFREEZ(3*MXATM),NFRZ
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /HSSCTL/ IHESSM,IHREP,HSSEND
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEFLD
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MCINP / METHNR,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /MIOPT / DTOL,CNVLOC,DIISTL,ALPHAX,
     *                NF(MAXNFRAGS),MF(MAXNFRAGS),NFRAGS,MSHIFT,ITERMI,
     *                ITLOC,IGUESS,MXDIIS,
     *                IOPT,MIFLAG,DIISON
      COMMON /MP2PAR/ OSPT,CODEMP,SCSPT,TOLMP,METHMP,NWDMP2,MEMPRI,
     *                MPPROP,NACORMP,NBCORMP,NOAMP,NOBMP,NORBMP,NBFMP,
     *                NOMITMP,MPCPHF,MAXITCMP
      COMMON /NEOJOB/ NEORUN,NELERM
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,NFMOPCM
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SVPOPT/ ISVP,NSVP
      COMMON /THERMD/ FREQ(3*MXATM),TEMP(10),SCLFAC,NTEMP,PRTSCN,GOTFRQ
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /VBDCMP/ DECOMP,NROW,NDEG
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
C     COSMO CHANGES
C
      PARAMETER (NPPA=1082)
      LOGICAL ISEPS,USEPS
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /COSFRQ/ FCOORD(3,MXATM),POS0(NPPA),
     *                COSZAN0(NPPA),EDIEL0,FINDEX,ICFREQ,ICFREQ1
C
      PARAMETER (TOLG=0.0005D+00, ZERO=0.0D+00)
C
C     ----- SET UP NAMELIST SIMULATION -----
C
      DATA FORCE/8HFORCE   /
      DATA QNAM/8HNVIB    ,8HVIBSIZ  ,8HPRTIFC  ,8HPURIFY  ,8HPROJCT  ,
     *          8HMETHOD  ,8HVIBANL  ,8HRDHESS  ,8HDECOMP  ,8HSCLFAC  ,
     *          8HTEMP    ,8HFREQ    ,8HPRTSCN  ,8HNPRT    ,8HNPUN    ,
     *          8HPULCOR  ,8HTEST    ,8HREDOVB  /
      DATA KQNAM/1,3,0,0,0,5,0,0,0,3,103,-3,0,1,1,0,0,0/
C
      DATA HSSIAN/8HHESSIAN /
      DATA ANAL,SNUM,FNUM,RNUM
     *      /8HANALYTIC,8HSEMINUM ,8HFULLNUM ,8HNUMERIC /
      DATA RHF,ROHF,UHF/8HRHF     ,8HROHF    ,8HUHF     /
      DATA      GVB,RMC/8HGVB     ,8HMCSCF   /
      DATA GENCI,GUGA,CIS/8HGENCI   ,8HGUGA    ,8HCIS     /
      DATA ZAPT/8HZAPT    /
      DATA NONE,RNONE/4HNONE,8HNONE    /
      DATA CHECK/8HCHECK   /
C
C     ----- MAIN DRIVER TO OBTAIN HESSIAN MATRIX -----
C
C     THE ENERGY, GRADIENT, HESSIAN, AND DIPOLE DERIVATIVE TENSOR
C     ARE OBTAINED BY READING THEM FROM CARDS, OR COMPUTING THEM.
C     THE ALPHA POLARIZABILITY DERIVATIVE TENSOR IS OBTAINED ONLY
C     BY READING IT FROM CARDS.  IT IS COMPUTED ONLY BY RUNTYP=RAMAN.
C
C     THE CALLING ARGUMENTS PERMIT YOU TO OVERRIDE THE VALUE OF
C     VIBANL GIVEN IN $FORCE BELOW, AND INDICATE THE SUCCESFUL
C     CONCLUSION OF A GEOMETRY SEARCH, USEFUL FOR $STATPT 'HSSEND'
C
      SOME = NPRINT.NE.-5 .AND. MASWRK
      NCOORD = 3*NAT
      NCF=6*NFRG
      GOTEG  = .FALSE.
      GOTEH  = .FALSE.
      GOTDDM = .FALSE.
      GOTADM = .FALSE.
C
C     ----- READ $FORCE INPUT GROUP -----
C
C        THE STRATEGY FOR DETERMINING THE DEFAULT DIFFERENTIATION
C        IS TO PERFORM AN ANALYTIC CALCULATION IF POSSIBLE, OR TO
C        CHOOSE SEMI-NUMERICAL DIFFERENTIATION OF ANALYTIC GRADIENTS.
C        IF THE ANALYTIC GRADIENT IS NOT CODED, THE USER IS PERMITTED
C        TO REQUEST FULLY NUMERIC DIFFERENTIATION, BUT BECAUSE THIS
C        IS REASONABLE ONLY FOR SMALL MOLECULES, WE NEVER PICK THAT
C        AS A DEFAULT (THEY MUST TYPE IT).
C
      METHOD=ANAL
C
C     ---- THE FOLLOWING METHODS ONLY ALLOW SEMI-NUMERIC HESSIANS ----
C        BASIS SETS WITH F OR G FUNCTIONS
C        UHF, OR GVB WITH MORE THAN ONE PERFECT PAIR
C        MCSCF IF IT DOES NOT USE FULL ACTIVE SPACE DETERMINANT CODE
C        CI, MP2 (RHF OR UHF), DFT
C        MOPAC WAVEFUNCTIONS
C        SIMOMM MODEL
C        SOLVENT MODELS: EFP, ONSAGER, PCM, OR COSMO
C        SCALAR RELATIVISTIC CORRECTIONS
C        APPLIED EXTERNAL ELECTRIC FIELD
C        SCF-MI
C     ALL OF WHICH HAVE ANALYTIC GRADIENTS PROGRAMMED.
C     IN ADDITION, THE PARTIAL HESSIAN ANALYSIS IS SEMI-NUMERIC
C
      CALL BASCHK(LMAX)
      IF(LMAX.GT.2)                             METHOD=SNUM
      IF(SCFTYP.EQ.UHF)                         METHOD=SNUM
      IF(SCFTYP.EQ.GVB  .AND.  NPAIR.GT.1)      METHOD=SNUM
      IF(SCFTYP.EQ.RMC  .AND.  CISTEP.EQ.GUGA)  METHOD=SNUM
      IF(SCFTYP.EQ.RMC  .AND.  CISTEP.EQ.GENCI) METHOD=SNUM
      IF(MPLEVL.GT.0)                           METHOD=SNUM
      IF(CITYP.NE.RNONE)                        METHOD=SNUM
      IF(DFTYPE.NE.RNONE)                       METHOD=SNUM
      IF(SCFTYP.EQ.RHF  .AND. TDDFTYP.NE.RNONE) METHOD=SNUM
      IF(MPCTYP.NE.NONE)                        METHOD=SNUM
      IF(IECP.EQ.5)                             METHOD=SNUM
      IF(QMMM)                                  METHOD=SNUM
      ISOLV = NFRG + IZRF + IPCM
      IF(ISOLV.GE.1  .OR.  ISEPS)               METHOD=SNUM
      IF(RMETHOD.NE.RNONE)                      METHOD=SNUM
      IF(EFLDL)                                 METHOD=SNUM
      IF(MIFLAG.EQ.1)                           METHOD=SNUM
      IF(RUNTYP.EQ.HSSIAN.AND.IFREEZ(1).NE.0)   METHOD=SNUM
      IF(NEORUN.EQ.1)                           METHOD=SNUM
C
C         IF WE HAVE THE ENERGY, BUT NOT THE ANALYTIC GRADIENTS...
C
      MPGRD = SCFTYP.EQ.RHF  .OR.  SCFTYP.EQ.UHF
     *                       .OR. (SCFTYP.EQ.ROHF  .AND.  OSPT.EQ.ZAPT)
      CIGRD = SCFTYP.EQ.RHF .AND. (CITYP.EQ.GUGA  .OR.  CITYP.EQ.CIS)
      IF(MPLEVL.GT.0  .AND.  .NOT.MPGRD)           METHOD=RNONE
      IF(CITYP.NE.RNONE  .AND.  .NOT.CIGRD)        METHOD=RNONE
      IF(CCTYP.NE.RNONE)                           METHOD=RNONE
      IF(VBTYP.NE.RNONE)                           METHOD=RNONE
      IF(SCFTYP.NE.RHF  .AND. TDDFTYP.NE.RNONE)    METHOD=RNONE
      IF(ISVP.EQ.1)                                METHOD=RNONE
C        NOTE THAT WE ARE LETTING EFP SLIDE HERE,
C        FOR TRADITIONAL REASONS,
C        BUT WE DON'T HAVE MP2/CI+EFP GRADIENTS CORRECTLY CODED.
      IF(MPLEVL.GT.0  .OR.  CITYP.NE.RNONE) THEN
         IF(IZRF.GE.1 .OR.  IPCM.EQ.1  .OR. ISEPS) METHOD=RNONE
      END IF
C
C         REMEMBER THIS DEFAULT FOR ERROR CHECKING LATER
C
      IANAL=0
      ISNUM=0
      IFNUM=0
      IF(METHOD.EQ.ANAL)  IANAL=1
      IF(METHOD.EQ.SNUM)  ISNUM=1
      IF(METHOD.EQ.RNONE) IFNUM=1
C
      FGONLY = NUM.EQ.0  .AND.  NFRG.GT.0
      IF(ISEPS) ICFREQ=1
      RDHESS=.FALSE.
      NVIB = 1
      MPRINT=NPRINT
      VIBSIZ=0.01D+00
      PURIFY=.FALSE.
      PRTIFC=.FALSE.
      VIBANL=.FALSE.
      IF(RUNTYP.EQ.HSSIAN) VIBANL=.TRUE.
      DECOMP=.FALSE.
      PROJCT=.FALSE.
      IF(IFREEZ(1).NE.0.AND.RUNTYP.EQ.HSSIAN) PROJCT=.TRUE.
      SCLFAC = 1.0D+00
      DO 100 I=1,10
         TEMP(I) = ZERO
  100 CONTINUE
      TEMP(1) = 298.15D+00
      KQNAM(12) = NCOORD*10+3
      DO 110 I=1,NCOORD
         FREQ(I) = ZERO
  110 CONTINUE
      PRTSCN = .FALSE.
      NPRT = 0
      NPUN = 0
      PULCOR=.FALSE.
      REDOVB=.TRUE.
C
C     "TEST" IS AN UNDOCUMENTED WAY TO AVOID THE CHECK MADE
C     BELOW AS TO WHEN ANALYTIC HESSIANS CAN BE REQUESTED.
C
      TEST=.FALSE.
C
      JRET=0
      CALL NAMEIO(IR,JRET,FORCE,NNAM,QNAM,KQNAM,
     *            NVIB,VIBSIZ,PRTIFC,PURIFY,PROJCT,METHOD,VIBANL,
     *            RDHESS,DECOMP,SCLFAC,TEMP,FREQ,PRTSCN,NPRT,NPUN,
     *            PULCOR,TEST,REDOVB,
     *            0,  0,0,0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,
     *    0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      NERR=0
      IF(JRET.EQ.2) NERR=NERR+1
C
C       'NUMERIC ' IS TO BE CONSIDERED AN ACCEPTABLE TYPO FOR 'SEMINUM '
      IF(METHOD.EQ.RNUM) METHOD=SNUM
C
      IF(NPRT.LT.0) NPRT=0
      IF(NPRT.GT.1) NPRT=1
      IF(NPUN.LT.0) NPUN=0
      IF(NPUN.GT.2) NPUN=2
C
      IF(SOME) WRITE(IW,9010) METHOD,NVIB,VIBSIZ,
     *                        RDHESS,PURIFY,PRTIFC,
     *                        VIBANL,DECOMP,PROJCT,
     *                        SCLFAC,PRTSCN,NPRT,
     *                        PULCOR,NPUN,REDOVB
C
      NTEMP=0
      DO 200 I=1,10
         IF(TEMP(I).NE.ZERO) NTEMP=NTEMP+1
  200 CONTINUE
      IF(SOME  .AND.  VIBANL) WRITE(IW,9012) NTEMP,(TEMP(I),I=1,NTEMP)
C
      NFREQ=0
      DO 210 I=1,NCOORD
         IF(FREQ(I).NE.ZERO) NFREQ=NFREQ+1
  210 CONTINUE
      IF(NFREQ.GT.0) THEN
         IF (SOME) WRITE(IW,9014) NFREQ,(FREQ(I),I=1,NFREQ)
         GOTFRQ=.TRUE.
      END IF
C
      IF(NVIB.LE.0  .OR.  NVIB.GT.2) NERR=NERR+1
      IF(ABS(VIBSIZ).GT.0.25D+00) NERR=NERR+1
C
      OK=.FALSE.
      IF(METHOD.EQ.ANAL) OK=.TRUE.
      IF(METHOD.EQ.SNUM) OK=.TRUE.
      IF(METHOD.EQ.FNUM) OK=.TRUE.
      IF(.NOT.OK) THEN
         IF(MASWRK) WRITE(IW,*)
     *       '*** ERROR: METHOD = ANALYTIC, SEMINUM, OR FULLNUM ONLY'
         NERR=NERR+1
      END IF
      IF(METHOD.EQ.ANAL) NHLEVL=0
      IF(METHOD.EQ.SNUM) NHLEVL=1
      IF(METHOD.EQ.FNUM) NHLEVL=2
C
      IF(NGLEVL.EQ.1) THEN
         IF(MASWRK) WRITE(IW,9090)
         NERR=NERR+1
      END IF
C
C        CHECK TO SEE IF DIFFERENTIATION METHOD IS REASONABLE
C     THE UNDOCUMENTED KEYWORD -TEST- ALLOWS SKIPPING THESE CHECKS.
C     OF COURSE, -TEST- SHOULD BE CHOSEN ONLY BY SOMEONE WHO IS
C     TRYING TO PROGRAM A MISSING METHOD.
C
      IF(.NOT.TEST) THEN
         IF(METHOD.EQ.ANAL  .AND.  (ISNUM+IFNUM).GT.0) THEN
            IF(MASWRK) WRITE(IW,9100)
            IF(MASWRK) WRITE(IW,9105) METHOD
            NERR=NERR+1
         END IF
         IF(METHOD.EQ.SNUM  .AND.  IFNUM.GT.0) THEN
            IF(MASWRK) WRITE(IW,9100)
            IF(MASWRK) WRITE(IW,9105) METHOD
            NERR=NERR+1
         END IF
      END IF
C
C        PRINT SOME FRIENDLY ADVISE IF A HIGHER METHOD IS FEASIBLE
C
      IF(MASWRK) THEN
         IF(METHOD.EQ.SNUM  .AND.  IANAL.EQ.1) WRITE(IW,9106) ANAL,SNUM
         IF(METHOD.EQ.FNUM  .AND.  IANAL.EQ.1) WRITE(IW,9106) ANAL,FNUM
         IF(METHOD.EQ.FNUM  .AND.  ISNUM.EQ.1) WRITE(IW,9106) SNUM,FNUM
      END IF
C
      IF(NERR.GT.0) THEN
         IF(MASWRK) WRITE(IW,*)
     *      'ERROR IN $FORCE INPUT, PLEASE FIX AND TRY AGAIN.'
         CALL ABRT
      END IF
C
      IF(GOTFRQ) GO TO 800
C
C     ----- READ INTERNALS, GENERATE B MATRIX -----
C     THIS MUST BE DONE AFTER -DECOMP- IS KNOWN, AS THE DECOMPOSITION
C     SOMETIMES USES MORE THAN 3N-6 COORDINATES.  THE Z-MATRIX
C     CODE PERMITS MORE THAN 3N-6 COORDINATES ONLY IF DECOMP=.T.
C
      CARTSN = NZVAR.EQ.0
      IF(NZVAR.GT.0 .AND. RUNTYP.EQ.HSSIAN) THEN
         CALL ZMATIN(CARTSN)
         CALL BANDBI
         CALL PZANDG(F,0)
      END IF
C
C     ----- OBTAIN THE INFORMATION FROM CARD INPUT -----
C
      IF(RDHESS) THEN
         CALL VALFM(LOADFM)
         LFCM   = LOADFM+1
         LDDM   = LFCM   + (NCOORD+NCF)*(NCOORD+NCF)
         LADM   = LDDM   + 3*(NCOORD+NCF)
C
C            NOTE THAT $ALPDR IS NOT SET UP TO WORK WITH ANY
C        --- SOLVATION MODELS, SO ONLY 6*NCOORDS WORDS ARE
C            NEEDED INSTEAD OF 6*(NCOORD+NCF)
C
         LEG    = LADM   + 6*NCOORD
         LAST   = LEG    + NCOORD
         NEED = LAST-LOADFM
         CALL GETFM(NEED)
C
C           LOOK FOR OPTIONAL GRADIENT VECTOR
C           LOOK FOR HESSIAN MATRIX
C           LOOK FOR OPTIONAL DIPOLE DERIVATIVE TENSOR
C           LOOK FOR OPTIONAL ALPHA POLARIZABILITY DERIVATIVE TENSOR
C
         CALL EGIN(X(LEG),NCOORD,GOTEG,' $GRAD  ')
         CALL FCMIN(X(LFCM),NCOORD+NCF,GOTEH)
         CALL DDMIN(X(LDDM),NCOORD+NCF,GOTDDM)
         CALL ADMIN(RAMALF,X(LADM),NCOORD,GOTADM)
         CALL RETFM(NEED)
         IF(GOTEH) THEN
            SCFOK=.TRUE.
            GO TO 700
         ELSE
            IF (MASWRK) WRITE(IW,9030)
            CALL ABRT
         END IF
      END IF
C
C     ----- OTHERWISE, THEY MUST BE COMPUTED -----
C
      IF(METHOD.EQ.ANAL) THEN
         IHESSM=2
         GOTDDM=.TRUE.
         GOTEG =.TRUE.
         CALL VALFM(LOADFM)
         LDDM = LOADFM+1
         LAST = LDDM  + 3*NCOORD
         NEED = LAST-LOADFM
         CALL GETFM(NEED)
         CALL VCLR(X(LDDM),1,3*NCOORD)
         CALL DAWRIT(IDAF,IODA,X(LDDM),3*NCOORD,34,0)
         CALL RETFM(NEED)
         CALL HSSANA(GOTWFN)
         SCFOK=.TRUE.
      END IF
C
      IF(METHOD.EQ.SNUM) THEN
         IHESSM=1
         GOTDDM=.TRUE.
         GOTEG =.TRUE.
         CALL VALFM(LOADFM)
         LFCM   = LOADFM + 1
         LDDM   = LFCM   +   (NCOORD+6*NFRG)**2
         LIST   = LDDM   + 3*(NCOORD+6*NFRG)
         LSKIP  = LIST   + 6*NFRG*NVIB
         LSKIP2 = LSKIP  + NAT+2*NFRG
         LWRK   = LSKIP2 + NAT*NAT
         LAST   = LWRK   + NUM*NUM
         NEED = LAST-LOADFM
         CALL GETFM(NEED)
         NATM = NAT
         CALL HSSNUM(X(LFCM),X(LDDM),X(LSKIP),X(LSKIP2),X(LWRK),
     *               X(LIST),NCOORD+6*NFRG,NAT+2*NFRG,NATM,
     *               NVIB,VIBSIZ,NPRT,NPUN,SCFOK,GOTWFN,REDOVB)
         CALL RETFM(NEED)
      END IF
C
      IF(METHOD.EQ.FNUM) THEN
         IHESSM=0
         IF(NFRG.GT.0) THEN
            IF(MASWRK) WRITE(IW,*)
     *          'FULLY NUMERIC HESSIAN NOT AVAILABLE WITH EFP'
            CALL ABRT
         END IF
         GOTDDM = .FALSE.
         GOTEG  = .FALSE.
         MYNC1  = 3*NAT
C           THE ACTUAL NUMBER OF POINTS DONE CANNOT EXCEED MXSP
         MXSP   = 1 + 24*NAT + 18*NAT*NAT
         CALL VALFM(LOADFM)
         LFCM   = LOADFM + 1         ! START OF FCM
         LDDM   = LFCM   + MYNC1**2  ! START OF DDM
         LSKIP  = LDDM   + 3*MYNC1   ! START OF SKIP ARRAY
         LATMP  = LSKIP  + NAT       ! START OF ATOM PAIR ARRAY
         LWRK   = LATMP  + NAT*NAT   ! STORT OF VIB0 ORBITALS
         LESP   = LWRK   + NUM*NUM   ! START OF SINGLE POINT ENERGIES
         LEG    = LESP   + MXSP      ! START OF GRADIENT
         LAST   = LEG    + NCOORD    ! LAST ADDRESS ON STACK
         NEED   = LAST - LOADFM
         CALL GETFM(NEED)
         CALL HSSFUL(X(LESP),X(LFCM),X(LSKIP),X(LATMP),X(LWRK),
     *               NAT,MYNC1,MXSP,VIBSIZ,GOTEG,SCFOK)
         CALL RETFM(NEED)
      END IF
C
C         THE HESSIAN COMPUTATION IS NOW FINISHED
C
  700 CONTINUE
C
C         HESSIAN MAY BE COMPUTED AS PART OF A GEOMETRY OPTIMIZATION,
C         ETC., AND THOSE RUNS MAY NOW TEST -IHESSM- TO KNOW THEY ARE
C         NO LONGER IN A HESSIAN CALCULATION PHASE...
C
      IHESSM=0
C
C     --- PURIFY THE HESSIAN,DIPOLE DERIVATIVES,ALPHA POLARIZABILITY ---
C
      NPRINT=MPRINT
      IF(NZVAR.GT.0 .AND. (PURIFY.OR.PRTIFC.OR.DECOMP.OR.PULCOR)) THEN
         CALL VALFM(LOADFM)
         LFCM   = LOADFM + 1
         LFCMI  = LFCM  + NCOORD*NCOORD
         LEG    = LFCMI + NROW*NROW
         LDDM   = LEG   + MAX(NCOORD,NROW)
         LDDN   = LDDM  + 3*NCOORD
         LADM   = LDDN  + 3*NVAR
         LADMI  = LADM  + 6*NCOORD
         LBMAT  = LADMI + 6*NVAR
         LBINV  = LBMAT + NVAR*NCOORD
         LWRK   = LBINV + NCOORD*NVAR
         LZMAT  = LWRK  + MAX(NCOORD,NROW)
         LIZMAT = LZMAT + NVAR
         LIWRK  = LIZMAT+ NZMAT
         LAST   = LIWRK + NVAR
         NEED = LAST - LOADFM - 1
         CALL GETFM(NEED)
C
C    PURIFY GRADIENT AND/OR PRINT INTERNAL GRADIENT
C
         IF(PRTIFC.OR.PURIFY) THEN
            CALL DAREAD(IDAF,IODA,X(LEG),NCOORD,3,0)
            CALL TRANG(X(LEG),NVAR,NCOORD)
            IF (PRTIFC) THEN
               CALL PZANDG(X(LEG),1)
               IF (NSYMC .EQ. 0)
     *            CALL PUZMAT(X(LZMAT),X(LIZMAT),X(LIWRK),NVAR,NZMAT)
            END IF
            IF (PURIFY) CALL TRANGB(X(LEG),NVAR,NCOORD)
         END IF
C
         CALL PURFCM(X(LFCM),X(LFCMI),X(LIZMAT),NCOORD,NROW,NVAR,
     *               NZMAT,PURIFY,PRTIFC,DECOMP)
C
         IF (PULCOR) THEN
C---         CALL HSSCOR(X(LFCM),X(LFCMI),X(LBINV),X(LEG),NCOORD,NVAR)
             IF(MASWRK) WRITE(IW,*) 'OPTION -PULCOR- IS INACTIVE'
         END IF
C
         IF(GOTDDM  .AND.  PURIFY)
     *      CALL PURDDM(X(LDDM),X(LDDN),X(LBMAT),X(LBINV),X(LWRK),
     *                  NCOORD,NVAR)
         IF(GOTADM .AND. PURIFY)
     *      CALL PURADM(X(LADM),X(LADMI),X(LBMAT),X(LBINV),X(LWRK),
     *                  NCOORD,NVAR)
         CALL RETFM(NEED)
      END IF
C
C     ----- PRINT/PUNCH GRADIENT, HESSIAN, DIPOLE DERIVATIVES,
C     ----- ALPHA POLARIZABILITY DERIVATIVES.
C
      CALL VALFM(LOADFM)
      LFCM = LOADFM + 1
      LEG  = LFCM   + (NCOORD+6*NFRG)*(NCOORD+6*NFRG)
      LDDM = LEG    + NCOORD
      LADM = LDDM   + 3*(NCOORD+6*NFRG)
      LAST = LADM   + 6*NCOORD
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      IF(NAT.GT.0) CALL DAREAD(IDAF,IODA,X(LEG) ,NCOORD,3,0)
      CALL DAREAD(IDAF,IODA,X(LFCM),(NCOORD+6*NFRG)**2, 4,0)
      IF(GOTDDM) THEN
         CALL DAREAD(IDAF,IODA,X(LDDM),3*(NCOORD+6*NFRG) ,34,0)
      ELSE
         CALL VCLR(X(LDDM),1,3*(NCOORD+6*NFRG))
      END IF
C
      IF(GOTEG .AND. .NOT. FGONLY) THEN
         GRMS = DDOT(NCOORD,X(LEG),1,X(LEG),1)
         GRMS = SQRT(GRMS/NCOORD)
         LOCMX = IDAMAX(NCOORD,X(LEG),1)
         GMAX = X(LEG-1+LOCMX)
         STATPT = GMAX.LT.TOLG  .AND.  GRMS.LT.(TOLG/3.0D+00)
      ELSE
         STATPT = .TRUE.
      END IF
C
      IF(RUNTYP.EQ.HSSIAN  .OR.  GOTWFN) THEN
         IF(GOTEG.AND.MASWRK) WRITE(IW,9040)
         IF(GOTEG)  CALL  EGOUT(X(LEG) ,NAT)
                    CALL FCMOUT(X(LFCM),NCOORD+6*NFRG)
         IF(POLAR  .AND.  EXETYP.NE.CHECK) THEN
            IF(IFREEZ(1).EQ.0) THEN
               CALL DAREAD(IDAF,IODA,APOL,6,251,0)
               CALL POLOUT(APOL)
            END IF
         END IF
         IF(GOTDDM .AND. .NOT.FGONLY) THEN
            IF(IFREEZ(1).EQ.0) THEN
               CALL DDMOUT(X(LDDM),NCOORD+6*NFRG)
            END IF
         END IF
         IF(GOTADM .AND. .NOT.FGONLY) THEN
            CALL POLOUT(RAMALF)
            CALL ADMOUT(X(LADM),NCOORD)
         END IF
      END IF
C
      IF(GOTEG)  CALL  EGPUN(X(LEG) ,NAT,' $GRAD  ')
                 CALL FCMPUN(X(LFCM),NCOORD+6*NFRG)
      IF(GOTDDM) CALL DDMPUN(X(LDDM),NCOORD+6*NFRG)
      IF(GOTADM) CALL ADMPUN(E,RAMALF,X(LADM),NCOORD)
      CALL RETFM(NEED)
C
C     ----- NORMAL COORDINATE VIBRATIONAL ANALYSIS -----
C
  800 CONTINUE
      IF(.NOT.VIBANL  .AND.  .NOT.VIBOVR) RETURN
      NC1 = NCOORD+6*NFRG
      NC2 = (NC1**2+NC1)/2
      NC3 = NC1**2
C
      NPART = NAT
      KMASS=0
      DO 820 IFRG=1,NFRG
         DO 810 III=1,NMPTS(IFRG)
            IF(FMASS(III+KMASS).GT.ZERO) NPART=NPART+1
  810    CONTINUE
         KMASS=KMASS+NMPTS(IFRG)
  820 CONTINUE
C
      CALL VALFM(LOADFM)
      LVEC   = LOADFM + 1
      LFCM   = LVEC   + NC3
      LE     = LFCM   + NC2
      LSCR   = LE     + NC1
      LIA    = LSCR   + NC1*8
      LRM    = LIA    + NC1
      LSVT   = LRM    + (NCOORD+21*NFRG)
      LSVR   = LSVT   + NC1*3
      LSVTT  = LSVR   + NC1*3
      LSVRT  = LSVTT  + NC1
      LCC    = LSVRT  + NC1
      LCOM   = LCC    + 3*NPART
      LZMS   = LCOM   + 3*NPART
      LBUF1  = LZMS   +   NPART
      LBUF2  = LBUF1  + NC3
      LDDM   = LBUF2  + NC3
      LADM   = LDDM   + NC1*3
      LAST   = LADM   + NC1*6
      NEED = LAST-LOADFM
      CALL GETFM(NEED)
      NATM = NAT+2*NFRG
      IF(NEORUN.EQ.0) THEN
         CALL FGMTRX(X(LVEC),X(LFCM),X(LE),X(LSCR),X(LIA),X(LRM),
     *               X(LSVT),X(LSVR),X(LSVTT),X(LSVRT),
     *               X(LCC),X(LCOM),X(LZMS),X(LBUF1),X(LBUF2),
     *               NC1,NC2,NPART,X(LDDM),GOTDDM,X(LADM),GOTADM,
     *               PROJCT,STATPT,SCFOK)
      ELSE
         CALL FGMTRX_NEO(X(LVEC),X(LFCM),X(LE),X(LSCR),X(LIA),X(LRM),
     *               X(LSVT),X(LSVR),X(LSVTT),X(LSVRT),
     *               X(LCC),X(LCOM),X(LZMS),X(LBUF1),X(LBUF2),
     *               NC1,NC2,NPART,X(LDDM),GOTDDM,X(LADM),GOTADM,
     *               PROJCT,STATPT,SCFOK)
      END IF
C
C        SAVE HARMONIC FREQUENCIES AND NORMAL MODES FOR EXTENDED TDHF
C
      IF(IGETOLI(1).NE.0) THEN
        CALL SVGFNM(FREQ,X(LVEC),NC1,LINEAR)
      END IF
C
      CALL RETFM(NEED)
      IF (MASWRK) WRITE(IW,9020)
      CALL TIMIT(1)
      RETURN
C
 9010 FORMAT(/5X,33(1H-)/
     *        5X,'HESSIAN MATRIX CONTROL PARAMETERS'/5X,33(1H-)/
     *        5X,'METHOD=',A8,3X,'NVIB  =',I8,3X,'VIBSIZ=',F8.5/
     *        5X,'RDHESS=',L8,3X,'PURIFY=',L8,3X,'PRTIFC=',L8/
     *        5X,'VIBANL=',L8,3X,'DECOMP=',L8,3X,'PROJCT=',L8/
     *        5X,'SCLFAC=',F8.5,3X,'PRTSCN=',L8,3X,'NPRT  =',I8/
     *        5X,'PULCOR=',L8,3X,'NPUN  =',I8,3X,'REDOVB=',L8)
 9012 FORMAT(5X,'THERMOCHEMISTRY WILL BE PRINTED FOR',I5,
     *          ' TEMPERATURES:'/(3X,5F12.5/))
 9014 FORMAT(5X,'HESSIAN COMPUTATION WILL BE SKIPPED AS'/
     *       5X,'A TOTAL OF',I5,' FREQUENCIES WERE INPUT:'/
     *      (5X,5F12.5/))
 9020 FORMAT(1X,'......END OF NORMAL COORDINATE ANALYSIS......')
 9030 FORMAT(1X,'RDHESS IS TRUE, BUT NO HESSIAN MATRIX WAS FOUND IN',
     *          ' YOUR INPUT FILE.')
 9040 FORMAT(/10X,15(1H-)/10X,'ENERGY GRADIENT'/10X,15(1H-))
 9090 FORMAT(/1X,'PLEASE DO NOT USE NUMGRD=.TRUE. DURING HESSIAN JOBS,'/
     *   1X,'INSTEAD USE METHOD=ANALYTIC, SEMINUM, OR FULLNUM ONLY.'/
     *   1X,'(IF NECESSARY, RUN THE HESSIAN IN A SEPARATE JOB THAN'/
     *   1X,'A GEOMETRY SEARCH (I.E. AVOID HESS=CALC OR HSSEND=.TRUE.)')
 9100 FORMAT(/10X,'**** ERROR IN DIFFERENTIATION METHOD ****'//
     *  1X,'ANALYTIC COMPUTATION OF THE HESSIAN IS IMPLEMENTED'/
     *  5X,'FOR S,P,D BASIS SETS (INCLUDING ECP USAGE) FOR AB INITIO',
     *     ' RHF, ROHF,'/
     *  5X,'GVB (OPEN SHELL OR TCSCF) AND MCSCF (CISTEP=ALDET)',
     *     ' WAVEFUNCTIONS.'//
     *  1X,'SEMI-NUMERIC HESSIANS CAN BE COMPUTED WHEN ANALYTIC',
     *     ' GRADIENTS ARE AVAILABLE:'/
     *  5X,'IF F,G FUNCTIONS ARE IN THE BASIS SET,'/
     *  5X,'IF USING GRID-BASED OR GRID-FREE DFT,'/
     *  5X,'IF UHF, OR IF GVB WITH MORE THAN ONE PERFECT PAIR,'/
     *  5X,'IF MCSCF USES SOME CISTEP OTHER THAN ALDET,'/
     *  5X,'IF MPLEVL=2 FOR RHF OR UHF WAVEFUNCTIONS,'/
     *  5X,'IF EXCITED STATE CIS FOR RHF WAVEFUNCTIONS,'/
     *  5X,'IF GROUND OR EXCITED STATE GUGA CI FOR RHF WAVEFUNCTIONS,'/
     *  5X,'IF SEMI-EMPIRICAL MODELS AM1, PM3, OR MNDO ARE USED,'/
     *  5X,'IF EFP, PCM, ONSAGER, OR COSMO SOLVENT MODELS ARE USED,'/
     *  5X,'IF THERE IS AN APPLIED EXTERNAL ELECTRIC FIELD,'/
     *  5X,'IF SCALAR RELATIVISTIC TRANSFORMATIONS ARE IN USE, OR'/
     *  5X,'IF RUNNING SCF-MI COMPUTATIONS.'//
     *  1X,'FULLY NUMERIC HESSIANS CAN BE COMPUTED WHEN ONLY THE',
     *     ' ENERGY IS AVAILABLE:'/
     *  5X,'WHEN MODEL CORE POTENTIAL INTEGRALS ARE USED,'/
     *  5X,'WHEN MPLEVL=2 FOR SCFTYP=ROHF/OSPT=RMP OR SCFTYP=MCSCF,'/
     *  5X,'WHEN CITYP IS SELECTED AND THE SCFTYP IS NOT RHF,'/
     *  5X,'WHEN THE CITYP SELECTED IS NOT CIS OR GUGA, OR'/
     *  5X,'WHEN ANY CCTYP IS SELECTED'/)
 9105 FORMAT(1X,'YOUR INPUT MISTAKENLY REQUESTED METHOD=',A8)
 9106 FORMAT(/15X,'* * * EFFICIENCY NOTE * * *'/
     *        1X,'THIS CALCULATION CAN BE RUN WITH DIFFERENTIATION',
     *           ' METHOD=',A8/
     *        1X,'THIS MIGHT BE MORE ACCURATE AND/OR FASTER THAN',
     *           ' YOUR CHOICE,'/
     *        1X,'BUT THE REQUESTED CALCULATION METHOD=',A8,
     *           ' IS PROCEEDING.'/)
      END
C*MODULE HESS    *DECK HSSANA
      SUBROUTINE HSSANA(GOTWFN)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL FIRST,SECND,CPHF,BOTH,MFIRST,MSECND,MCPHF,
     *        PK,PANDK,BLOCK,SVGPAR,SVDSKW,MODRV,POLAR,
     *        GOTWFN,DIRSCF,DIRTRF,FDIFF,GOPARR,DSKWRK,MASWRK,
     *        ABEL,ABELPT,REGEN
C
      PARAMETER (MXATM=2000, MXAO=8192, MXSH=5000, MXRT=100)
C
      COMMON /CHFINF/ CHFSLV,CPTYPE,NWDCHF,POLAR
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,S,SZ,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACT,NDTORB,
     *                MNA,MNB,K,KSTDET,IROOT,IPURES,MAXW1,NITER,MAXP,
     *                NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /FUNCT / E,EG(3*MXATM)
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORB,NCONF(MXAO),NHAM
      COMMON /HSSPAR/ FIRST,SECND,CPHF,BOTH,MFIRST,MSECND,MCPHF
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /OPTSCF/ DIRSCF,FDIFF
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IKREP,IEF,IP_F,NFMOPCM
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRFOPT/ CUTTRF,NWDTRF,MPTRAN,ITRFAO,NOSYMT,IPURTF,DIRTRF
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (ZERO=0.0D+00, TOL=1.0D-09)
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF     /
      DATA GVB,RMC/8HGVB     ,8HMCSCF   /
      DATA TYPEMO/8HMO      /
      DATA CHECK/8HCHECK   /
C
C     ----- COMPUTE HESSIAN MATRIX ANALYTICALLY -----
C     (FOR GRADIENT ONLY, SET ONLY FIRST AND MFIRST TRUE.)
C
      FIRST=.TRUE.
      SECND=.TRUE.
      CPHF=.TRUE.
      BOTH=.TRUE.
      MFIRST=.TRUE.
      MSECND=.TRUE.
      MCPHF=.TRUE.
C
C     ----- IF SAMCHF GRADIENTS OR NACME, SHUT OFF HESS -----
C
      IF(SCFTYP.EQ.RMC) THEN
         NSTATS = 0
         DO ISTAT=1,MXRT
           IF(WSTATE(ISTAT).GT.TOL) NSTATS = NSTATS+1
         ENDDO
C
         IF(NSTATS.GT.1) THEN
           SECND = .FALSE.
           BOTH = .FALSE.
           MSECND = .FALSE.
         ENDIF
      END IF
C
C     ----- READ $CPHF INPUT -----
C
      CALL CPINP
      MODRV = CPTYPE.EQ.TYPEMO
C
C           NOT PROGRAMMED FOR CONTINUUM SOLVENT MODELS AT PRESENT
C
      IF(IPCM.NE.0  .OR.  IZRF.NE.0) THEN
         IF(MASWRK) WRITE(IW,9010)
         CALL ABRT
      END IF
C
C     ----- GET THE ENERGY AND THE WAVEFUNCTION -----
C
      IF(.NOT.GOTWFN) THEN
         CALL ENERGX
         IF(E.EQ.ZERO  .AND.  EXETYP.NE.CHECK) THEN
            IF(MASWRK) WRITE(IW,9040)
            CALL ABRT
         END IF
      END IF
C
C         MCSCF HESSIAN IS PURELY A DISTRIBUTED MEMORY PARALLEL
C         INTEGRAL TRANSFORMATION AND RESPONSE EQUATION CODE.
C
      IF(SCFTYP.EQ.RMC) THEN
         GOPARR=.TRUE.
         GO TO 200
      END IF
C
C     ----- PREPARE CORRECT INTEGRAL LIST ---
C     THIS MUST BE -J- FORMAT, WHICH SHOULD HAVE BEEN ANTICIPATED.
C
      IF(PK) THEN
         IF(MASWRK) WRITE(IW,*) 'HSSANA: SOMEHOW THIS IS A -PK- FILE'
         CALL ABRT
      END IF
C
C     PARALLEL TRANSFORMATION MAY HAVE REQUESTED DISTRIBUTED AO STORAGE.
C     NON-ABELIAN POINT GROUPS NEED C1 INTEGRAL FILE FOR TRANSFORMATION.
C     AO-BASED CPHF MUST HAVE A C1 INTEGRAL LIST FOR THE RESPONSE STEP.
C     IN THESE CASES, THE AO INTEGRAL FILE NEEDS TO BE REGENERATED.
C
      ABEL = ABELPT()
      NTSAVE = NT
      SVDSKW = DSKWRK
      SVGPAR = GOPARR
      IF(GOPARR) DSKWRK = .TRUE.
C
      IF(MODRV) THEN
         REGEN= (GOPARR.AND.ITRFAO.EQ.1)  .OR.   (.NOT.ABEL)
      ELSE
         REGEN = NT.GT.1
      END IF
C
      IF(CPHF  .AND.  REGEN  .AND.  .NOT.DIRSCF) THEN
         IF(MODRV) THEN
            IF(ITRFAO.EQ.1) THEN
               IF(MASWRK) WRITE(IW,9060)
               GOPARR = .FALSE.
            END IF
            IF(.NOT.ABEL) THEN
               IF(MASWRK) WRITE(IW,9070)
               NT=1
            END IF
         ELSE
            IF(MASWRK) WRITE(IW,9050)
            NT = 1
         END IF
         CALL JANDK
      END IF
C
      GOPARR = SVGPAR
      NT = NTSAVE
C
C     ----- INTEGRAL TRANSFORMATION -----
C     WE DON'T NEED THE TRANFORMED INTEGRALS UNTIL THE CPHF STEP,
C     BUT THE JOB BOTTLENECK IS DISK STORAGE DURING THE ORDERING
C     AND TRANSFORMATION.  WE'D RATHER FIND THAT WAS LIMITED BEFORE
C     THE MAIN CPU BOTTLENECK, WHICH IS THE HESSIAN INTEGRALS.
C     OF COURSE, AO BASIS RESPONSE EQUATIONS CAN SKIP THE TRANSFORM.
C
  200 CONTINUE
      NOCC = 0
      IF(SCFTYP.EQ.RHF)  NOCC=NA
      IF(SCFTYP.EQ.UHF)  NOCC=NA
      IF(SCFTYP.EQ.ROHF) NOCC=NA
      IF(SCFTYP.EQ.GVB)  NOCC=NORB
      IF(SCFTYP.EQ.RMC)  NOCC=NCORSV+NACT
C
C          THE MCSCF RESPONSE EQUATIONS USE DISTRIBUTED INTEGRALS
C          WHEN RUNNING IN PARALLEL, ALL OTHERS USE DISK INTEGRALS.
C          HOWEVER, RHF CAN RUN FROM AO INTEGRALS IN WHICH CASE
C          WE SKIP THE TRANSFORMATION ENTIRELY.
C
      IF(CPHF) THEN
        NPR = 0
        IF(NPRINT.EQ.-5) NPR=-5
        IF(DIRSCF) DIRTRF=.TRUE.
        IF(SCFTYP.EQ.RMC .AND. GOPARR) THEN
          CALL TRFMCX(NPR,0,NOCC,NQMT,.FALSE.,.FALSE.,
     *                .TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE.,
     *                .FALSE.,.FALSE.,.TRUE.)
        ELSE IF(MODRV) THEN
          CALL TRFMCX(NPR,0,NOCC,NQMT,.FALSE.,.FALSE.,
     *                .FALSE.,.FALSE.,.FALSE.,.FALSE.,.FALSE.,
     *                .FALSE.,.FALSE.,.TRUE.)
        END IF
      END IF
C
C     ----- OPEN DERIVATIVE FOCK (DERIVATIVE LAGRANGIAN) MATRIX FILE
C
      NFT18 = 18
      IF(GOPARR) DSKWRK = .FALSE.
      CALL SEQOPN(NFT18,'FOCKDER','UNKNOWN',.FALSE.,'UNFORMATTED')
C
C     ----- ONE ELECTRON INTEGRAL DERIVATIVE TERMS -----
C
      CALL STVDD
C
C     ----- DIPOLE INTEGRALS AND THEIR DERIVATIVES -----
C
      CALL CALCOM(XC,YC,ZC)
      CALL DIPINT(XC,YC,ZC,.FALSE.)
      CALL DDINT
C
C     ----- TWO ELECTRON INTEGRAL DERIVATIVE TERMS -----
C
      CALL DDERJK
C
C     ----- SOLVE THE COUPLED PERTURBED EQUATIONS -----
C
      IF(GOPARR) DSKWRK=.TRUE.
C
      IF(CPHF) THEN
         IF(SCFTYP.EQ.RMC) THEN
            CALL CPMCX
         ELSE
            CALL CPHFX
         END IF
      END IF
C
      DSKWRK = SVDSKW
      RETURN
C
 9010 FORMAT(//1X,'*** ERROR ***'/
     *       1X,'NO ANALYTIC HESSIAN FOR -SCRF- OR -PCM- AT PRESENT')
 9040 FORMAT(//1X,'ENERGY DID NOT CONVERGE...ABORTING HESSIAN'//)
 9050 FORMAT(/1X,'REGENERATING AO INTEGRAL LIST IN C1 POINT GROUP',
     *           ' FOR AO-DRIVEN CPHF')
 9060 FORMAT(/1X,'GENERATING REPLICATED AO INTEGRAL LISTS FOR THE',
     *           ' TRANSFORMATION STEP.')
 9070 FORMAT(/1X,'REGENERATING AO INTEGRAL LIST IN C1 POINT GROUP',
     *           ' FOR TRANSFORMATION STEP.')
      END
C*MODULE HESS    *DECK HSSCOR
C---      SUBROUTINE HSSCOR(FCM,FCMI,BINV,EGI,NCOORD,NVAR)
C---C
C---      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C---C
C---      DIMENSION FCM(NCOORD,NCOORD),FCMI(NVAR,NVAR)
C---      DIMENSION BINV(NCOORD,NVAR),EGI(NVAR)
C---C
C---C        THIS IS CODE THERESA ADDED TO WORK ON PULAY'S
C---C        EMPIRICAL GRADIENT CORRECTIONS
C---C
C---      RETURN
C---      END
C*MODULE HESS    *DECK HSSNUM
      SUBROUTINE HSSNUM(FCM,DDM,SKIP,SKIP2,WRK,LIST,NC1,NSKIP,NATM,
     *                  NVIB,VIBSIZ,NPRT,NPUN,SCFOK,GOTWFN,REDOVB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000, MXAO=8192,
     *           MXFRG=1050, MXFGPT=12000, MXRT=100,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL SKIP(NSKIP),SKIP2(NATM,NATM),SCFOK,REDOVB,ABEL
      LOGICAL GOPARR,DSKWRK,MASWRK,MCCI,RSTART,EOF,FGONLY,RHFEFP,GOTWFN
C
      DIMENSION FCM(NC1,NC1),DDM(3,NC1),WRK(*),LIST(*)
      DIMENSION D(2),SVENGY(MXRT+15),ENGYIN(MXRT+15)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FGROL / DEFTO(3,MXFRG),TORQO(3,MXFRG)
      COMMON /FMCOM / XX(1)
      COMMON /FRAME / U1,U2,U3,V1,V2,V3,W1,W2,W3,X0,Y0,Z0
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FROZNC/ IFREEZ(3*MXATM),NFRZ
      COMMON /FUNCT / E,EG(3*MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3*MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /NEOJOB/ NEORUN,NELERM
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /SOOPT / NOSO
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TMVALS/ TI,TX,TIM
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /XYZPRP/ DUM(3),DIP(3),DUMM(32)
C
      EQUIVALENCE (ENGYIN(1),ENUCR)
C
      PARAMETER (UNIT=0.52917724924D+00, ONE=1.0D+00,TOL=1.0D-10)
      PARAMETER (SMALL=1.0D-06)
C
      DATA RHF,RMC/8HRHF     ,8HMCSCF   /
      DATA GUGA/8HGUGA    /, RNONE/8HNONE    /, NONE/4HNONE/
      DATA CHECK/8HCHECK   /
      DATA READMO/8HMOREAD  /
      DATA HSSIAN/8HHESSIAN /
C
C     ----- COMPUTE HESSIAN MATRIX NUMERICALLY -----
C
      MCCI = SCFTYP.EQ.RMC  .OR.  CITYP.EQ.GUGA
      EOF = .FALSE.
      D(1) =  VIBSIZ
      D(2) = -VIBSIZ
      IVIB   = 0
      IATOM  = 0
      ICOORD = 0
      NFRG2  = 2*NFRG
      NCF    = 6*NFRG
      NCOORD = NC1-NCF
      SCFOK = .TRUE.
      DEL = VIBSIZ*NVIB*UNIT
      IREAD=0
C
C          THERE ARE SPECIAL CASE CODES IF THERE ARE ONLY FRAGMENTS,
C          OR FOR RHF-LEVEL SOLUTES WITH FRAGMENTS PRESENT.
C
      FGONLY = NUM.EQ.0  .AND.  NFRG.GT.0
      RHFEFP = NUM.GT.0  .AND.  NFRG.GT.0  .AND.
     *         (SCFTYP.EQ.RHF .AND. MPLEVL.EQ.0 .AND.
     *          CITYP.EQ.RNONE  .AND.  DFTYPE.EQ.RNONE)
      IF (MASWRK  .AND.  RHFEFP) WRITE(IW,9010)
C
C     ----- OPEN SUPPLEMENTAL PUNCH FILE FOR RESTART DATA -----
C
      IFCM=4
      CALL SEQOPN(IFCM,'IRCDATA','NEW',.FALSE.,'FORMATTED')
      IF (MASWRK) WRITE(IFCM,8010) TITLE
C
C     ----- CHECK FOR A RESTART RUN (PRESENCE OF $VIB)  -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $VIB   ',IEOF)
      IF(IEOF.EQ.1) THEN
         RSTART=.FALSE.
      ELSE
         RSTART=.TRUE.
         IF(MASWRK) WRITE(IW,9020)
      END IF
      IF(GOTWFN) RSTART=.FALSE.
C
C     ----- INITIAL GEOMETRY (THE VIB 0 POINT) -----
C
C     EVALUATE THE VIB 0 WAVEFUNCTION AND GRADIENT, IF NOT RESTARTING.
C     THE FLAG -GOTWFN- INDICATES HESSIAN AFTER GEOMETRY SEARCH,
C     IN WHICH CASE WE KNOW THE WAVEFUNCTION, ENERGY, AND GRADIENT.
C
      IF(.NOT.RSTART  .AND.  .NOT.GOTWFN) THEN
         IF(RHFEFP) THEN
            CALL EFPENG(IVIB,IATOM,ICOORD,SCFOK)
            CALL STVDER
         ELSE
            IF(FGONLY) THEN
               CALL EFPENG(IVIB,IATOM,ICOORD,SCFOK)
               CALL EFGRAD
            ELSE
               CALL VALFOR(NCOORD,IVIB,IATOM,ICOORD,NPRT,NPUN,SCFOK)
            END IF
         END IF
      END IF
C
C     FOR A RESTART, WE MAY NEED TO REEVALUATE THE VIB 0 WAVEFUNCTION.
C     BECAUSE WE HANDLE THE INITIAL GUESS HERE, WE INCREMENT NEVALS
C     SO THAT CALLS TO ENERGX JUST EXECUTE ORTHDN, NOT FULL GUESS.
C
      IF(RSTART) THEN
         IF(FGONLY) THEN
            IF(MASWRK) WRITE(IW,9030)
            CALL EFPENG(IVIB,IATOM,ICOORD,SCFOK)
         ELSE
            CALL ONEEI
            CALL GUESMO(GUESS)
            NEVALS = NEVALS + 1
C
C   FOR SCF, IF WE ARE GIVEN THE ORBITALS, WE ASSUME THEY ARE THE
C   CONVERGED VIB0 ORBITALS.  FOR MCSCF, WE REPEAT THE COMPUTATION
C   OF THE VIB0 POINT ANYWAY, SINCE WE NEED THE CI EIGENVECTOR TOO.
C   NO POINT IN REEVALUATING MP2 ENERGY IF WE REGENERATE VIB0 ORBS.
C
C   IN THE CASE OF FRAGMENTS, 2E- INTEGRALS NEED BE CALCULATED ONCE,
C   AS THE VIB 0 INTEGRALS ARE USED BY ALL FRAGMENT DISPLACEMENTS.
C
            IF(GUESS.NE.READMO  .OR.  MCCI) THEN
               IF(MASWRK) WRITE(IW,9040)
               MPSAVE = MPLEVL
               MPLEVL = 0
               CALL ENERGX
               MPLEVL = MPSAVE
            ELSE
               IF(NFRG.GT.0) CALL JANDK
            END IF
         END IF
C
C           REPOSITION TO THE $VIB GROUP, AND READ VIB 0 GRADIENT
C
         CALL SEQREW(IR)
         CALL FNDGRP(IR,' $VIB   ',IEOF)
         IF(MASWRK) WRITE(IW,9050)
         CALL RDVIB(IR,IW,IREAD,NCOORD,IVIB,IATOM,ICOORD,E,EG,DIP,EOF)
         IF(EOF) THEN
            IF(MASWRK) WRITE(IW,*) 'PROBLEM EXTRACTING VIB0 GRADIENT'
            CALL ABRT
         END IF
      END IF
C
C     ----- INITIALIZE HESSIAN AND DIPOLE DERIVATIVE MATRICES -----
C          WITH THE VIB 0 GEOMETRY'S GRADIENT AND DIPOLE VECTOR
C
      CALL SETDDM(DDM,DIP,DEL,0,NVIB,NC1,NVIB)
      IF(.NOT.RHFEFP  .OR.  RSTART  .OR.  GOTWFN) THEN
         IF(MASWRK) WRITE(IW,9060)
         CALL EGOUT(EG,NAT)
         CALL SETFCM(FCM,NC1,NCOORD,EG,IVIB)
         CALL PUVIB(IFCM,IW,RSTART,NCOORD,IVIB,IATOM,ICOORD,E,EG,DIP)
      END IF
C
C        SAVE THE VIB 0 WAVEFUNCTION, GRADIENT, AND Q-MATRIX.
C        NOTE THAT AB INITIO ATOM GRADIENT COMPONENTS ARE ON DAF,
C        BUT THE FRAGMENT GRADIENT COMPONENTS SAVED IN COMMON: UGH.
C
      E0 = E
      DO I=1,MXRT+15
         SVENGY(I) = ENGYIN(I)
      ENDDO
C
      NUM2 = (NUM*NUM+NUM)/2
      NUM3 = NUM*NUM
      IF(.NOT.FGONLY) THEN
         CALL FFSVMO(IVIB,WRK)
         IF(MCCI) CALL FFSVCI(IVIB)
         CALL DAREAD(IDAF,IODA,WRK,NUM3, 45,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM3,276,0)
         CALL DAREAD(IDAF,IODA,WRK,NUM2, 11,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM2,277,0)
         CALL DAREAD(IDAF,IODA,WRK,NUM2, 12,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM2,278,0)
         IF(MPCTYP.EQ.NONE) THEN
            CALL DAREAD(IDAF,IODA,WRK,NUM2, 13,0)
            CALL DAWRIT(IDAF,IODA,WRK,NUM2,279,0)
         END IF
      END IF
      IF(.NOT.RHFEFP  .OR.  GOTWFN  .OR.  RSTART) THEN
         IF(.NOT.FGONLY) CALL DAWRIT(IDAF,IODA,EG,NCOORD,29,0)
         DO IFR=1,NFRG
            DO J=1,3
               DEFTO(J,IFR)=DEFT(J,IFR)
               TORQO(J,IFR)=TORQ(J,IFR)
            ENDDO
         ENDDO
      END IF
C
      IF(EXETYP.EQ.CHECK) GO TO 500
C
C     ----- MINIMIZE PRINT/PUNCH OUTPUT AT DISPLACED GEOMETRIES -----
C
      NPRINT = -5
                    NPUNCH = 0
      IF(NPUN.EQ.1) NPUNCH = 1
      IF(NPUN.EQ.2) NPUNCH = 2
C
      IF(TIM.GE.TIMLIM) GO TO 800
      TVIB0 = TIM
C
C     ----- FIND WHICH SYMMETRY UNIQUE ATOMS MUST BE DISPLACED -----
C     FRAGMENT RUNS ARE ALWAYS IN C1 SYMMETRY, ALL ATOMS MUST DISPLACE.
C
      IF(NFRG.EQ.0) THEN
         DO 260 IAT = 1,NAT
            SKIP(IAT) = .TRUE.
            DO IT = 1,NT
               IC = MAPCTR(IAT,IT)
               IF(IC.GT.IAT) GO TO 260
            ENDDO
            SKIP(IAT) = .FALSE.
  260    CONTINUE
      ELSE
         DO IAT = 1,NAT+NFRG2
            SKIP(IAT) = .FALSE.
         ENDDO
      END IF
C
C     ----- COMPUTE FRAGMENT'S PORTION OF HESSIAN MATRIX -----
C
      IF(NFRG.GT.0) THEN
C
         CALL EFPRST(IFCM,FCM,NFRG,NAT,DEL,NVIB,DDM,RSTART,
     *               NVST,NFTODO,IREAD,LIST,REDOVB,GOTWFN)
C
C     THERE IS SPECIAL CODE HERE FOR THE CASE OF RHF WITH FRAGMENTS.
C
         IF(RHFEFP) THEN
            IF(NFTODO.GT.0) THEN
               L2=(NUM*NUM+NUM)/2
               CALL VALFM(LOADFM)
               IGRD=LOADFM+1
               IDEN=IGRD+NFTODO*NC1
               IENG=IDEN+NFTODO*L2
               IDIP=IENG+NFTODO
               LAST=IDIP+NFTODO*3
               NEED=LAST-LOADFM-1
               CALL GETFM(NEED)
C
C              FIRST, EVALUATE THE WAVEFUNCTION AND THE ONE ELECTRON
C              TERM OF THE GRADIENT FOR ALL FRAGMENT DISPLACEMENTS.
C
               CALL EFPVIB(NAT,NUM,DEL,NVIB,NPUN,D,IFCM,FCM,DDM,WRK,
     *                     XX(IDEN),XX(IGRD),XX(IENG),XX(IDIP),
     *                     NVST,NFTODO,LIST,
     *                     SCFOK,FGONLY,RHFEFP,MCCI,RSTART)
C
C              THEN, USING THE STORED 1-PARTICLE DENSITY MATRICES FOR
C              EACH FRAGMENT DISPLACEMENT, EVALUATE THE 2E- PARTS OF
C              THE GRADIENT IN 1 PASS OVER THE 2E-GRADIENT INTEGRALS.
C
               CALL EFPGRD(IFCM,NAT,NVIB,NPUN,FCM,
     *                     XX(IDEN),XX(IGRD),XX(IENG),XX(IDIP),
     *                     NVST,NFTODO,LIST)
               CALL RETFM(NEED)
            END IF
         ELSE
            CALL EFPVIB(NAT,NUM,DEL,NVIB,NPUN,D,IFCM,FCM,DDM,WRK,
     *                  DUMMY,DUMMY,DUMMY,DUMMY,
     *                  NVST,NFTODO,LIST,
     *                  SCFOK,FGONLY,RHFEFP,MCCI,RSTART)
         END IF
C               PREPARE TO READ POSSIBLE ATOMIC DISPLACEMENTS IN $VIB
         IF (IREAD.EQ.1) RSTART=.TRUE.
      END IF
C
C       SAVE THE POINT GROUP
C     FOR SOME RUNS, ALL DISPLACEMENTS MUST BE LOWERED CLEAR TO C1.
C
      NTSAVE = NT
      NOSOSV = NOSO
      IF(MPLEVL.NE.0)    CALL SYMOFF
      IF(CITYP.NE.RNONE) CALL SYMOFF
      IF(SCFTYP.EQ.RMC)  CALL SYMOFF
      IF(ISPHER.GT.0)    CALL SYMOFF
      IF(NEORUN.EQ.1)    CALL SYMOFF
C
C     ----- IDENTIFY FROZEN ATOMS -----
C           IF THERE IS ANY ATOM OR COORDINATE FROZEN,
C           THE FRAGMENTS (IF ANY)  WILL BE FROZEN AUTOMATICALLY.
C
      DO 290 I = 1, NCOORD
        IF(IFREEZ(I).NE.0) THEN
          IAT = NINT(1.0+AINT((REAL(IFREEZ(I))-0.1)/3.0))
          SKIP(IAT) = .TRUE.
        END IF
  290 CONTINUE
C
      IF((IFREEZ(1).NE.0) .AND. (RUNTYP.EQ.HSSIAN)) THEN
         DO 300 IAT = NAT+1, NAT+NFRG2
           SKIP(IAT) = .TRUE.
  300    CONTINUE
      END IF
C
C
C     ----- COMPUTE AB INITIO ATOM'S PORTION OF HESSIAN MATRIX -----
C
      NDISPL=0
      IRST=0
      DO 390 IVIB = 1,NVIB
         DO 380 IATOM = 1,NAT
           IF(SKIP(IATOM)) GO TO 380
           NVA = 3*(IATOM-1)
           DO 370 ICOORD = 1,3
             NV = NVA+ICOORD
             C(NV) = C(NV)+D(IVIB)
C
             IF(NT.EQ.1)   GO TO 360
             IF(NFRG.NE.0) GO TO 360
C
C   RETAIN SOME SYMMETRY FOR SOME DISPLACEMENTS IN SOME CIRCUMSTANCES.
C   THIS IS POSSIBLE ONLY FOR CS, CN, CNV, CNH, DNH, AND FOR THESE,
C   ONLY FOR NON-CI, NON-MP2, NON-MCSCF, NON-SPHERICAL HARMONICS.
C   FOR UNSUPPORTED CASES, SYMMETRY WAS ALREADY TURNED OFF ABOVE.
C
             NT = 1
             IF (
     *          ( IGROUP.EQ.2
     *       .AND. ABS((C(NVA+1)-X0)*(PTR(1,4)-ONE)
     *                +(C(NVA+2)-Y0)*(PTR(2,5)-ONE)
     *                +(C(NVA+3)-Z0)*(PTR(3,6)-ONE))   .LE. TOL
     *       .AND. ABS(PTR(1,4)+PTR(2,5)+PTR(3,6)-ONE) .LE. TOL ) .OR.
     *            ( (IGROUP.EQ.4 .OR. IGROUP.EQ.7)
     *       .AND. ABS((C(NVA+1)-X0)*(PTR(1,4)-ONE)
     *                +(C(NVA+2)-Y0)*(PTR(2,5)-ONE)
     *                +(C(NVA+3)-Z0)*(PTR(3,6)-ONE)) .LE. TOL
     *       .AND. (ABS(PTR(1,4)-ONE) .LE. TOL .OR.
     *              ABS(PTR(2,5)-ONE) .LE. TOL .OR.
     *              ABS(PTR(3,6)-ONE) .LE. TOL ) )
     *                        ) NT = NTSAVE
             IF (
     *          ( (IGROUP.EQ.6 .OR. IGROUP.EQ.9)
     *       .AND. ABS((C(NVA+1)-X0)*(PTR(1,4)-ONE)
     *                +(C(NVA+2)-Y0)*(PTR(2,5)-ONE)
     *                +(C(NVA+3)-Z0)*(PTR(3,6)-ONE)) .LE. TOL
     *       .AND. (ABS(PTR(1,4)-ONE) .LE. TOL .OR.
     *              ABS(PTR(2,5)-ONE) .LE. TOL .OR.
     *              ABS(PTR(3,6)-ONE) .LE. TOL ) )
     *                        ) NT = MIN ( NAXIS , NTSAVE )
C
C        FOR A SUBGROUP OF THE FULL GROUP, WE CAN EXPLOIT AO INTEGRAL
C        SYMMETRY (VIA -NT-), BUT WE CANNOT USE -Q- MATRIX SYMMETRY.
C
             NOSO  = NOSOSV
             IF(NT.LT.NTSAVE) NOSO=1
C
C        IF WE ARE DOING A RESTART FROM A $VIB GROUP, ATTEMPT
C        TO READ THIS ATOM'S GRADIENT AND DIPOLE VECTORS.
C
  360        CONTINUE
             IF(RSTART) THEN
C
C        MCSCF RUNS MAY HAVE READ OTHER INPUT GROUPS, AND SO WE
C        MUST RELOCATE THE CORRECT POSITION WITHIN THE $VIB GROUP
C
                IF(SCFTYP.EQ.RMC  .AND.  IRST.EQ.1) THEN
                   IRST=0
                   IF(MASWRK) WRITE(IW,9070)
                   CALL SEQREW(IR)
                   CALL FNDGRP(IR,' $VIB   ',IEOF)
                   IF(IEOF.EQ.1) CALL ABRT
                   KCOORD=0
                   KVIB=0
                   KKATOM=0
                   CALL RDVIB(IR,IW,IREAD,NCOORD,KVIB,KKATOM,KCOORD,
     *                        E,EG,DIP,EOF)
                   IF(EOF) THEN
                      IF(MASWRK) WRITE(IW,*) 'REPOS TO VIB0 PROBLEM'
                      CALL ABRT
                   END IF
                   KVREAD = 0
                   DO 363 KVIB = 1,NVIB
                      DO 362 KKATOM = 1,NAT
                         IF(SKIP(KKATOM)) GO TO 362
                         DO 361 KCOORD = 1,3
                            IF(KVREAD.EQ.NDISPL) GO TO 365
                            KVREAD = KVREAD+1
                            CALL RDVIB(IR,IW,IREAD,NCOORD,
     *                                 KVIB,KKATOM,KCOORD,
     *                                 E,EG,DIP,EOF)
                            IF(EOF) THEN
                               IF(MASWRK) WRITE(IW,*)
     *                               'VIB REPOSITION PROBLEM'
                               CALL ABRT
                            END IF
  361                    CONTINUE
  362                 CONTINUE
  363              CONTINUE
  365              CONTINUE
                   IF(MASWRK) WRITE(IW,9080)
                END IF
C
C        READ THE CORRECT ENTRY FOR THIS DISPLACEMENT
C
                CALL RDVIB(IR,IW,IREAD,NCOORD,IVIB,IATOM,ICOORD,
     *                     E,EG,DIP,EOF)
             END IF
C
C        IF WE HIT END OF FILE IN $VIB GROUP, OR FOUND THAT THIS
C        ENERGY WAS NOT CONVERGED (ENERGY READ AS ZERO), WE NEED
C        TO PRETEND WE DID NOT MANAGE TO READ THIS DISPLACEMENT.
C
             IF(EOF) RSTART=.FALSE.
             IF(ABS(E).LT.SMALL  .AND.  RSTART  .AND.  REDOVB) THEN
                IF(MASWRK) WRITE(IW,9090)
                IRST=1
                RSTART=.FALSE.
             END IF
C
C        IF WE WERE UNABLE TO READ THE GRADIENT OF THIS VIBRATION,
C        WE MUST NOW SET ABOUT COMPUTING IT...
C
             IF(.NOT.RSTART) THEN
                CALL VALFOR(NCOORD,IVIB,IATOM,ICOORD,NPRT,NPUN,SCFOK)
                CALL FFSVMO(IVIB,WRK)
                IF(MCCI) CALL FFSVCI(IVIB)
             END IF
C
C        MOVE THIS DISPLACED ATOM BACK TO WHERE IT CAME FROM
C
             C(NV) = C(NV)-D(IVIB)
C
C        SAVE GRADIENT AND DIPOLE OF THIS DISPLACEMENT
C
             CALL PUVIB(IFCM,IW,RSTART,NCOORD,IVIB,IATOM,ICOORD,
     *                  E,EG,DIP)
             CALL SETFCM(FCM,NC1,NV,EG,IVIB)
             CALL SETDDM(DDM,DIP,DEL,NV,IVIB,NC1,NVIB)
C
             IF(IRST.EQ.1) THEN
                RSTART=.TRUE.
             END IF
             NDISPL=NDISPL+1
             TAVG = (TIM-TVIB0)/NDISPL
             IF (TIMLIM-TIM  .LT.  TAVG) GO TO 800
  370      CONTINUE
  380    CONTINUE
  390 CONTINUE
C
C     ----- WE ARE DONE WITH ALL ATOM AND FRAGMENT DISPLACEMENTS -----
C
C     ---- RESTORE THE VIB 0 ENERGY AND GRADIENT ----
C     (SET FROZEN ATOMS' GRADIENTS TO ZERO)
C
      E = E0
      DO I=1,MXRT+15
         ENGYIN(I) = SVENGY(I)
      ENDDO
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
C
      IF(NAT.GT.0) THEN
         CALL DAREAD(IDAF,IODA,EG,NCOORD,29,0)
         IF((IFREEZ(1).NE.0).AND.(RUNTYP.EQ.HSSIAN)) THEN
            DO 400 I = 1, NCOORD
               IF(IFREEZ(I).NE.0) EG(IFREEZ(I)) = 0.0D+00
  400       CONTINUE
         END IF
         CALL DAWRIT(IDAF,IODA,EG,NCOORD, 3,0)
      END IF
      DO 76 IFR=1,NFRG
      DO 76 J=1,3
         DEFT(J,IFR)=DEFTO(J,IFR)
         TORQ(J,IFR)=TORQO(J,IFR)
 76   CONTINUE
C
C        RESTORE THE VIB 0 POINT GROUP
C
      NT = NTSAVE
      IF(MPLEVL.NE.0)    CALL SYMON
      IF(CITYP.NE.RNONE) CALL SYMON
      IF(SCFTYP.EQ.RMC)  CALL SYMON
      IF(ISPHER.GT.0)    CALL SYMON
      IF(NEORUN.EQ.1)    CALL SYMON
      NOSO = NOSOSV
C
C        RESTORE THE VIB 0 ORBITALS AND DENSITY MATRIX
C        RESTORE THE VIB 0 Q-MATRIX AND THE H,S,T INTEGRALS
C
      IF(NAT.GT.0) THEN
         IVIB=3
         CALL FFSVMO(IVIB,WRK)
         IF(MCCI) CALL FFSVCI(IVIB)
         CALL DAREAD(IDAF,IODA,WRK,NUM3,276,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM3, 45,0)
         CALL DAREAD(IDAF,IODA,WRK,NUM2,277,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM2, 11,0)
         CALL DAREAD(IDAF,IODA,WRK,NUM2,278,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM2, 12,0)
         IF(MPCTYP.EQ.NONE) THEN
            CALL DAREAD(IDAF,IODA,WRK,NUM2,279,0)
            CALL DAWRIT(IDAF,IODA,WRK,NUM2, 13,0)
         END IF
      END IF
C
C     ----- COMPLETE BOTH NUMERICAL DIFFERENTIATIONS -----
C     ----- SYMMETRIZE THE FORCE CONSTANT MATRIX -----
C     ----- SYMMETRIZE THE DIPOLE DERIVATIVE TENSOR -----
C
  500 CONTINUE
      IF(EXETYP.EQ.CHECK) THEN
         CALL VCLR(FCM,1,NC1*NC1)
         DO 510 I=1,NC1
            FCM(I,I) = I
  510    CONTINUE
         CALL DAWRIT(IDAF,IODA,FCM,NC1*NC1,4,0)
         CALL VCLR(DDM,1,3*NCOORD)
         CALL DAWRIT(IDAF,IODA,DDM,3*NC1,34,0)
      ELSE
         CALL SYMFCM(FCM,SKIP2,EG,NCOORD,NCF,NAT,VIBSIZ,NVIB)
         CALL SYMDDM(DDM,SKIP,NC1,NAT)
      END IF
C
      IF (MASWRK) WRITE(IW,9100)
      CALL TIMIT(1)
      RETURN
C
C     ----- HIT TIME LIMIT -----
C
  800 CONTINUE
      IF (MASWRK) WRITE(IW,9110)
      CALL ABRT
      STOP
C
 8010 FORMAT('ENERGY/GRADIENT/DIPOLE RESTART DATA FOR RUNTYP=HESSIAN'/
     *         10A8)
C
 9010 FORMAT(/1X,'SPECIAL CODE WILL BE USED FOR RHF + FRAGMENT',
     *           ' DIFFERENTIATION.'/
     *        1X,'WAVEFUNCTIONS FOR ALL FRAGMENT DISPLACEMENTS WILL',
     *           ' BE EVALUATED FIRST,'/
     *        1X,'TOGETHER WITH 1E- AND FRAGMENT DERIVATIVE',
     *           ' CONTRIBUTIONS.'/
     *        1X,'THEN ALL 2E- CONTRIBUTIONS TO THE GRADIENT WILL BE',
     *           ' DONE IN ONE GRADIENT PASS.')
 9020 FORMAT(/5X,'A $VIB GROUP WAS FOUND IN YOUR INPUT.'/
     *        5X,'THEREFORE, THIS IS A NUMERICAL HESSIAN RESTART.')
 9030 FORMAT(/1X,'THIS RESTART WILL REEVALUATE THE FRAGMENT ONLY',
     *           ' ENERGY AT VIB 0 GEOMETRY.'/)
 9040 FORMAT(/5X,'EITHER THIS IS AN MCSCF FUNCTION, IN WHICH CASE THE'/
     *     5X,'VIB 0 CI EIGENVECTOR IS UNKNOWN, OR YOU DID NOT'/
     *     5X,'READ IN CONVERGED VIB 0 ORBITALS BY GUESS=MOREAD.'//
     *     5X,'THEREFORE THE VIB 0 WAVEFUNCTION MUST BE REEVALUATED.')
 9050 FORMAT(/5X,'READING $VIB RESTART DATA FROM $VIB CARDS...')
 9060 FORMAT(/1X,'THE VIB 0 GEOMETRY GRADIENT IS')
 9070 FORMAT(/5X,'REPOSITIONING TO THE CORRECT POINT IN $VIB GROUP...')
 9080 FORMAT(5X,'CORRECT POSITION IN $VIB GROUP IS NOW RELOCATED.'/)
 9090 FORMAT(1X,'RECALCULATING THE GRADIENT FOR THIS POINT SINCE THE',
     *          ' ENERGY WAS'/1X,'NOT CONVERGED IN THE PREVIOUS RUN')
 9100 FORMAT(' ...... END OF NUMERICAL HESSIAN CALCULATION ......')
 9110 FORMAT(//1X,'*** INTERNAL TIME LIMIT REACHED ****'/
     *         1X,'THERE IS NOT ENOUGH TIME LEFT TO START THE NEXT',
     *            ' NUMERICAL DISPLACEMENT.'/
     *         1X,'THIS RUN MUST BE RESTARTED WITH A $VIB GROUP.'//)
      END
C*MODULE HESS    *DECK HSSFUL
      SUBROUTINE HSSFUL(ESP,FCM,SKIP,ITRANS,WRK,
     *                  MYNAT,NC1,MXSP,VIBSIZ,GOTEG,SCFOK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,MCCI,GOTEG,SCFOK
C
      DIMENSION ESP(MXSP),FCM(NC1,NC1),ITRANS(MYNAT,MYNAT),WRK(*)
      DIMENSION AVIB(3,2),BVIB(3,2),CVIB(3,2),DVIB(3,2),RVIB(3,2)
      DIMENSION ADIS(3,2),BDIS(3,2),CDIS(3,2),DDIS(3,2),RDIS(3,2)
      DIMENSION V(3,3)
C
      CHARACTER*2 XYZ(3)
C
      LOGICAL SKIP(MYNAT)
      LOGICAL SYM,SYMFLP,AEQB,AEQC,AEQD,BEQC
C
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TMVALS/ TI,TX,TIM
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, TWO=2.0D+00, TOL=1.0D-09)
C
      DATA RMC/8HMCSCF   /
      DATA GUGA/8HGUGA    /
      DATA NONE/4HNONE/
      DATA CHECK/8HCHECK   /
      DATA ENDWRD/8H $END   /
      DATA XYZ/' X',' Y',' Z'/
C
C     ----- COMPUTE HESSIAN MATRIX NUMERICALLY FROM ENERGY VALUES -----
C       WRITTEN NOVEMBER 2003 BY RYAN OLSON AT IOWA STATE UNIVERSITY
C
      IVIB = 0
      MCCI = SCFTYP.EQ.RMC  .OR.  CITYP.EQ.GUGA
      CALL VCLR(ESP,1,MXSP)
      CALL VCLR(FCM,1,NC1*NC1)
C
C     ----- FIND AND COUNT THE SYMMETRY UNIQUE ATOMS -----
C
      NUAT = 0
      DO 50 IAT = 1,NAT
         SKIP(IAT) = .TRUE.
         DO IT = 1,NT
            IC = MAPCTR(IAT,IT)
            IF(IC.GT.IAT) GO TO 50
         ENDDO
         NUAT = NUAT + 1
         SKIP(IAT) = .FALSE.
   50 CONTINUE
C
C     CLEAR THE ITRANS MATRIX
C     -----------------------
      DO IATOM = 1,NAT
        DO JATOM = 1,NAT
          ITRANS(IATOM,JATOM) = 0
        END DO
      END DO
C
C     ----- LOOK FOR A RESTART FILE -----
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $VIB2  ',IEOF)
      IF(IEOF.NE.1 .AND. MASWRK) THEN
        WRITE(IW,*)
     *    'READING RESTART DATA FROM $VIB2 GROUP IN INPUT FILE...'
C               READ VIB0 RESTART ENERGY
        READ(IR,8050) WORD,ISP,EVAL
        CALL UPRCAS(WORD,8)
C
        DO WHILE(WORD.NE.ENDWRD)
           IF(ISP .GT. MXSP) THEN
             WRITE(IW,9040)
             CALL ABRT
           END IF
           IF(ISP.GT.0) ESP(ISP) = EVAL
           READ(IR,8050) WORD,ISP,EVAL
           CALL UPRCAS(WORD,8)
        END DO
      END IF
      IF(GOPARR) CALL DDI_BCAST(472,'F',ESP,MXSP,MASTER)
C
C     ----- GET ENERGY AND WAVEFUNCTION AT INPUT (VIB0) GEOMETRY -----
C
      SCFOK  = .TRUE.
      CALL ENERGX
      IF(E.EQ.ZERO .AND. EXETYP.NE.CHECK) THEN
         IF(MASWRK) WRITE(IW,*)
     *      'NO FULLY NUMERIC HESSIAN, SINCE ENERGY DID NOT CONVERGE.'
         CALL ABRT
      END IF
C
      IF(TIM.GE.TIMLIM) GO TO 800
      TVIB0 = TIM
C
C     ----- OPEN SUPPLEMENTAL PUNCH FILE FOR RESTART DATA -----
C
      IFCM=4
      CALL SEQOPN(IFCM,'IRCDATA','NEW',.FALSE.,'FORMATTED')
      IF (MASWRK) WRITE(IFCM,8000)
      IF (MASWRK) WRITE(IFCM,8010) 0,E
C
C     ----- SAVE THE VIB0 ENERGY, WAVEFUNCTION, ET CETERA -----
C
      E0 = E
      CALL FFSVMO(IVIB,WRK)
      IF(MCCI) CALL FFSVCI(IVIB)
      NUM2 = (NUM*NUM+NUM)/2
      CALL DAREAD(IDAF,IODA,WRK,NUM*NUM, 45,0)
      CALL DAWRIT(IDAF,IODA,WRK,NUM*NUM,249,0)
      CALL DAREAD(IDAF,IODA,WRK,NUM2   , 11,0)
      CALL DAWRIT(IDAF,IODA,WRK,NUM2   ,277,0)
      CALL DAREAD(IDAF,IODA,WRK,NUM2   , 12,0)
      CALL DAWRIT(IDAF,IODA,WRK,NUM2   ,278,0)
      IF(MPCTYP.EQ.NONE) THEN
         CALL DAREAD(IDAF,IODA,WRK,NUM2, 13,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM2,279,0)
      END IF
C
C  FOR THE TWO BODY INTERACTIONS NEEDED FOR THE HESSIAN, WE TAKE THE
C  DIRECT PRODUCT OF SET OF ATOMS WITH ITSELF.  THIS GIVE US ALL PAIR
C  INTERACTIONS.  SINCE THESE INTERACTION ARE SYMMETRIC, WE ONLY USE
C  THE LOWER TRIANGLUAR PORTION WHERE JATOM <= IATOM.  THE FOLLOWING
C  BLOCK OF CODE LOOPS OVER THE SYMMETRY OPERATIONS OF THE GROUP AND
C  APPLIES EACH SYMMETRY OPERATION TO THE PAIR.  TO DECIDE WHICH PAIR
C  INTERACTIONS TO CALCULATE, WE ADOPT THE "MAXIMUM PROCRASTINATION"
C  STRATEGY USED IN GAMESS.  IF ANY TRANSFORMATION MAPS THE IATOM TO
C  A SYMMETRY EQUIVALENT ATOM WITH A HIGHER VALUE OR THE TRANSFORMATION
C  MAPS THE IATOM TO A HIGHER OR EQUAL VALUE AND THE JATOM TO A HIGHER
C  VALUE, THAT BLOCK IS SKIPPED BECAUSE IT WILL BE CALCULATED LATER.
C
C  THUS ONLY THE IATOM,JATOM PAIRS THAT HAVE AN ITRANS VALUE EQUAL TO 1
C  (THE IDENTITY OPERATION) NEED TO BE CALCULATED.
C
      DO IATOM = 1,NAT
        DO JATOM = 1,IATOM
          IMAX = IATOM
          JMAX = JATOM
          ITRANS(IATOM,JATOM) = 1
          DO IT = 1,NT
            IF(MAPCTR(IATOM,IT).GT.IMAX) THEN
              IMAX = MAPCTR(IATOM,IT)
              JMAX = MAPCTR(JATOM,IT)
              ITRANS(IATOM,JATOM) = IT
            END IF
            IF(MAPCTR(IATOM,IT).GE.IMAX
     *        .AND. MAPCTR(JATOM,IT).GT.JMAX) THEN
              IMAX = MAPCTR(IATOM,IT)
              JMAX = MAPCTR(JATOM,IT)
              ITRANS(IATOM,JATOM) = IT
            END IF
          END DO
        END DO
      END DO
C
C     ----- ESTIMATE THE NUMBER OF SINGLE POINTS NEEDED -----
C     THE CODE BELOW MAY ELIMINATE EVEN MORE POINTS THAN THIS ESTIMATION
C
      NTSP = 1 + 6*NUAT
      DO IATOM = 1,NAT
        DO JATOM = 1,IATOM
          IF(ITRANS(IATOM,JATOM).EQ.1) THEN
          DO ICOORD = 1,3
            DO JCOORD = 1,3
              IF(.NOT.(JATOM.EQ.IATOM .AND. JCOORD.LE.ICOORD)) THEN
                NTSP = NTSP + 4
              END IF
            END DO
          END DO
          END IF
        END DO
      END DO
      IF(MASWRK) WRITE(IW,9000) NTSP-1,MXSP-1
C
C     ----- MINIMIZE PRINT/PUNCH OUTPUT AT DISPLACED GEOMETRIES -----
C
      NPRINT = -5
      NPUNCH = 0
C
      NSP = 0
      DENOM = TWO * VIBSIZ
      DENOM = DENOM * DENOM
      HTERM = 1 / DENOM
      DO IATOM = 1,NAT
        DO 100 JATOM = 1,IATOM
          IF(ITRANS(IATOM,JATOM).NE.1) GO TO 100
          DO ICOORD = 1,3
            DO 110 JCOORD = 1,3
C
C        BEFORE DISPLACING "THE" COORDINATES, SAVE'EM
C
              CISAVE = C(ICOORD,IATOM)
              CJSAVE = C(JCOORD,JATOM)
C
              AEQB = .FALSE.
              AEQC = .FALSE.
              AEQD = .FALSE.
              BEQC = .FALSE.
C     CLEAR DISPLACEMENT VECTORS
              CALL VCLR(AVIB,1,6)
              CALL VCLR(BVIB,1,6)
              CALL VCLR(CVIB,1,6)
              CALL VCLR(DVIB,1,6)
C     INITIALIZE DISPLACEMENT VECTORS
              AVIB(ICOORD,1) =  VIBSIZ
              AVIB(JCOORD,2) =  VIBSIZ
              BVIB(ICOORD,1) =  VIBSIZ
              BVIB(JCOORD,2) = -VIBSIZ
              CVIB(ICOORD,1) = -VIBSIZ
              CVIB(JCOORD,2) =  VIBSIZ
              DVIB(ICOORD,1) = -VIBSIZ
              DVIB(JCOORD,2) = -VIBSIZ
C
C                FORM BVIB,CVIB AND DVIB
C
              DO I = 1,3
                ADIS(I,1) = AVIB(I,1) + C(I,IATOM)
                ADIS(I,2) = AVIB(I,2) + C(I,JATOM)
                BDIS(I,1) = BVIB(I,1) + C(I,IATOM)
                BDIS(I,2) = BVIB(I,2) + C(I,JATOM)
                CDIS(I,1) = CVIB(I,1) + C(I,IATOM)
                CDIS(I,2) = CVIB(I,2) + C(I,JATOM)
                DDIS(I,1) = DVIB(I,1) + C(I,IATOM)
                DDIS(I,2) = DVIB(I,2) + C(I,JATOM)
              END DO
C
C     DETERMINE WHETHER THE COORDINATES ARE SYMMETRIC
C     CREATE AVIB WITH THE COORDINATES (ICOORD/JCOORD) REVERSED
C
              CALL VCLR(RVIB,1,6)
              RVIB(JCOORD,1) = VIBSIZ
              RVIB(ICOORD,2) = VIBSIZ
              SYM    = IATOM.EQ.JATOM
              SYMFLP = .FALSE.
C
              DO 115 IT = 1,NT
                IF(SYM) GO TO 115
                IATOM_T = MAPCTR(IATOM,IT)
                JATOM_T = MAPCTR(JATOM,IT)
                IF(.NOT.((IATOM.EQ.IATOM_T .AND. JATOM.EQ.JATOM_T) .OR.
     *             (IATOM.EQ.JATOM_T .AND. JATOM.EQ.IATOM_T))) GO TO 115
C
C                 TRANSFORM RVIB UNDER IT
C
                CALL VROT2(RVIB,IT,IATOM,JATOM,RDIS)
C
C          CHECK IF RDIS IS EQUIVALENT TO ANY OF THE DISPLACEMENTS
C          WHEN THE COORDINATES ARE NOT REVERSED
C          CHECK RDIS == ADIS
C
                CALL VECCMP(RDIS,ADIS,6,TOL,SYM)
                IF(SYM) GO TO 115
C
C          CHECK RDIS == BDIS
C
                CALL VECCMP(RDIS,BDIS,6,TOL,SYM)
                IF(SYM) SYMFLP = .TRUE.
                IF(SYM) GO TO 115
C
C          CHECK RDIS == CDIS
C
                CALL VECCMP(RDIS,CDIS,6,TOL,SYM)
                IF(SYM) SYMFLP = .TRUE.
                IF(SYM) GO TO 115
C
C           CHECK RDIS == DDIS
C
                CALL VECCMP(RDIS,DDIS,6,TOL,SYM)
C
  115         CONTINUE
C
C           SKIP SYMMETRIC COORDINATES
C
              IF(SYM .AND. JCOORD.GT.ICOORD) GO TO 110
C
              DO 120 IT = 1,NT
                IATOM_T = MAPCTR(IATOM,IT)
                JATOM_T = MAPCTR(JATOM,IT)
                IF(.NOT.((IATOM.EQ.IATOM_T .AND. JATOM.EQ.JATOM_T) .OR.
     *             (IATOM.EQ.JATOM_T .AND. JATOM.EQ.IATOM_T))) GO TO 120
C
C             TRANSFORM T (=A) UNDER IT
C
                CALL VROT2(AVIB,IT,IATOM,JATOM,ADIS)
C
C           CHECK IF A == D
C
                IF(.NOT.AEQD) CALL VECCMP(ADIS,DDIS,6,TOL,AEQD)
C
C           IATOM == JATOM && ICOORD == JCOORD ==> NO MORE TESTS
C
                IF(IATOM.EQ.JATOM .AND. ICOORD.EQ.JCOORD) GO TO 120
C
C           CHECK IF A == B
C
                IF(.NOT.AEQB) CALL VECCMP(ADIS,BDIS,6,TOL,AEQB)
C
C           CHECK IF A == C
C
                IF(.NOT.AEQC) CALL VECCMP(ADIS,CDIS,6,TOL,AEQC)
C
  120         CONTINUE
C
C           DEAL WITH THE DIAGONAL ELEMENTS
C
              IF(IATOM.EQ.JATOM.AND.ICOORD.EQ.JCOORD) THEN
                C(ICOORD,IATOM) = CISAVE + 2*VIBSIZ
                IF(MASWRK) THEN
                   IF(ESP(NSP+1).EQ.ZERO) WRITE(IW,9200) NSP+1
                   WRITE(IW,9210) IATOM,XYZ(ICOORD),C(ICOORD,IATOM)
                END IF
                CALL NUMNRG(ESP,NSP,E1,WRK,MCCI,SCFOK,IFCM)
                C(ICOORD,IATOM) = CISAVE
C
                E2 = E1
                IF(.NOT.AEQD) THEN
                   C(ICOORD,IATOM) = CISAVE - 2*VIBSIZ
                   IF(MASWRK) THEN
                      IF(ESP(NSP+1).EQ.ZERO) WRITE(IW,9200) NSP+1
                      WRITE(IW,9210) IATOM,XYZ(ICOORD),C(ICOORD,IATOM)
                   END IF
                   CALL NUMNRG(ESP,NSP,E2,WRK,MCCI,SCFOK,IFCM)
                   C(ICOORD,IATOM) = CISAVE
                END IF
C
C           CALC CONTRIBUTION TO THE FCM IMMEDIATELY
C
                IX = 3*(IATOM-1) + ICOORD
                FCM(IX,IX) = HTERM*(E1 + E2 - 2*E0)
C
                IF(.NOT.GOTEG) EG(ICOORD,IATOM) = (E1 - E2)/(4*VIBSIZ)
C
                TAVG = (TIM-TVIB0)/NSP
                IF (TIMLIM-TIM  .LT.  1.1D+00*TAVG) GO TO 800
C
                GO TO 110
              END IF
C
C            RESET T AS DISPLACEMENT B
C
              DO 130 IT = 1,NT
                IATOM_T = MAPCTR(IATOM,IT)
                JATOM_T = MAPCTR(JATOM,IT)
                IF(.NOT.((IATOM.EQ.IATOM_T .AND. JATOM.EQ.JATOM_T) .OR.
     *             (IATOM.EQ.JATOM_T .AND. JATOM.EQ.IATOM_T))) GO TO 130
C
C           TRANSFORM T (=B) UNDER IT
C
                CALL VROT2(BVIB,IT,IATOM,JATOM,BDIS)
C
C           CHECK T (=B) == C
C
                IF(.NOT.BEQC) CALL VECCMP(BDIS,CDIS,6,TOL,BEQC)
C
  130         CONTINUE
C
C           DEAL WITH THE OFF-DIAGONAL ELEMENTS
C           IF A == B || A == C, THEN THE FCM ELEMENT IS ZERO
C
              IF(AEQB .OR. AEQC) GO TO 110
C
C           CALCULATE AVIB
C
              C(ICOORD,IATOM) = CISAVE + VIBSIZ
              C(JCOORD,JATOM) = CJSAVE + VIBSIZ
              IF(MASWRK) THEN
                 IF(ESP(NSP+1).EQ.ZERO) WRITE(IW,9200) NSP+1
                 WRITE(IW,9220) IATOM,XYZ(ICOORD),C(ICOORD,IATOM),
     *                          JATOM,XYZ(JCOORD),C(JCOORD,JATOM)
              END IF
              CALL NUMNRG(ESP,NSP,E1,WRK,MCCI,SCFOK,IFCM)
              C(ICOORD,IATOM) = CISAVE
              C(JCOORD,JATOM) = CJSAVE
C
C           GET OR CALCULATE BVIB
C
              E2 = E1
              IF(.NOT.(AEQB .OR. (AEQC .AND. BEQC))) THEN
                C(ICOORD,IATOM) = CISAVE + VIBSIZ
                C(JCOORD,JATOM) = CJSAVE - VIBSIZ
                IF(MASWRK) THEN
                   IF(ESP(NSP+1).EQ.ZERO) WRITE(IW,9200) NSP+1
                   WRITE(IW,9220) IATOM,XYZ(ICOORD),C(ICOORD,IATOM),
     *                            JATOM,XYZ(JCOORD),C(JCOORD,JATOM)
                END IF
                CALL NUMNRG(ESP,NSP,E2,WRK,MCCI,SCFOK,IFCM)
                C(ICOORD,IATOM) = CISAVE
                C(JCOORD,JATOM) = CJSAVE
              END IF
C
C          GET OR CALCULATE CVIB
C
              E3 = E1
              IF(BEQC) E3 = E2
              IF(.NOT.(AEQC .OR. BEQC)) THEN
                C(ICOORD,IATOM) = CISAVE - VIBSIZ
                C(JCOORD,JATOM) = CJSAVE + VIBSIZ
                IF(MASWRK) THEN
                   IF(ESP(NSP+1).EQ.ZERO) WRITE(IW,9200) NSP+1
                   WRITE(IW,9220) IATOM,XYZ(ICOORD),C(ICOORD,IATOM),
     *                            JATOM,XYZ(JCOORD),C(JCOORD,JATOM)
                END IF
                CALL NUMNRG(ESP,NSP,E3,WRK,MCCI,SCFOK,IFCM)
                C(ICOORD,IATOM) = CISAVE
                C(JCOORD,JATOM) = CJSAVE
              END IF
C
C          GET OR CALCULATE DVIB
C
              E4 = E1
              IF(.NOT.AEQD) THEN
                C(ICOORD,IATOM) = CISAVE - VIBSIZ
                C(JCOORD,JATOM) = CJSAVE - VIBSIZ
                IF(MASWRK) THEN
                   IF(ESP(NSP+1).EQ.ZERO) WRITE(IW,9200) NSP+1
                   WRITE(IW,9220) IATOM,XYZ(ICOORD),C(ICOORD,IATOM),
     *                            JATOM,XYZ(JCOORD),C(JCOORD,JATOM)
                END IF
                CALL NUMNRG(ESP,NSP,E4,WRK,MCCI,SCFOK,IFCM)
                C(ICOORD,IATOM) = CISAVE
                C(JCOORD,JATOM) = CJSAVE
              END IF
C
C          CALC CONTRIBUTION TO THE FCM IMMEDIATELY
C
              TMP = HTERM*(E1-E2-E3+E4)
              IX = 3*(IATOM-1) + ICOORD
              IY = 3*(JATOM-1) + JCOORD
              FCM(IX,IY) = TMP
              FCM(IY,IX) = TMP
C
              IF(SYM) THEN
                IF(SYMFLP) TMP = -TMP
                IX = 3*(IATOM-1) + JCOORD
                IY = 3*(JATOM-1) + ICOORD
                FCM(IX,IY) = TMP
                FCM(IY,IX) = TMP
              END IF
C
              TAVG = (TIM-TVIB0)/NSP
              IF (TIMLIM-TIM  .LT.  1.1D+00*TAVG) GO TO 800
C
  110       CONTINUE ! END JCOORD LOOP
          END DO     ! END ICOORD LOOP
  100   CONTINUE     ! END JATOM LOOP
      END DO         ! END IATOM LOOP
C
      IF (MASWRK) WRITE(IFCM,8040)
C
C     --- RECOVER THE VIB0 ENERGY, WAVEFUNCTION, ET CETERA
C
      E = E0
      IVIB = 3
      CALL FFSVMO(IVIB,WRK)
      IF(MCCI) CALL FFSVCI(IVIB)
      CALL DAREAD(IDAF,IODA,WRK,NUM*NUM,249,0)
      CALL DAWRIT(IDAF,IODA,WRK,NUM*NUM, 45,0)
      CALL DAREAD(IDAF,IODA,WRK,NUM2   ,277,0)
      CALL DAWRIT(IDAF,IODA,WRK,NUM2   , 11,0)
      CALL DAREAD(IDAF,IODA,WRK,NUM2   ,278,0)
      CALL DAWRIT(IDAF,IODA,WRK,NUM2   , 12,0)
      IF(MPCTYP.EQ.NONE) THEN
         CALL DAREAD(IDAF,IODA,WRK,NUM2,279,0)
         CALL DAWRIT(IDAF,IODA,WRK,NUM2, 13,0)
      END IF
C
C     IF EXETYP=CHECK, GENERATE A DUMMY HESSIAN AND QUIT
C
      IF(EXETYP.EQ.CHECK) THEN
        IF(MASWRK) WRITE(IW,9030) NSP
        CALL VCLR(EG,1,3*NAT)
        GOTEG=.TRUE.
        CALL VCLR(FCM,1,NC1*NC1)
        DO I = 1,3*NAT
           FCM(I,I) = I
        END DO
        CALL DAWRIT(IDAF,IODA,EG, 3*NC1  ,3,0)
        CALL DAWRIT(IDAF,IODA,FCM,NC1*NC1,4,0)
        RETURN
      END IF
C
C     SYMMETRIZE THE SKELETON FCM AND GRADIENT
C
C     SYMMETRIZE THE GRADIENT, IF NEEDED.
C     -----------------------------------
      IF(.NOT.GOTEG) THEN
        DO 500 IATOM = 1,NAT
          IF(.NOT.SKIP(IATOM)) GO TO 500
          DO 510 IT = 1,NT
            JATOM = MAPCTR(IATOM,IT)
            IF(SKIP(JATOM)) GO TO 510
            DO 520 JT = 1,NT
              IF(MAPCTR(JATOM,JT).NE.IATOM) GO TO 520
              N = 3*(JT-1)
              DO I = 1,3
                TMP = ZERO
                DO M = 1,3
                  TMP = TMP + EG(M,JATOM)*PTR(I,N+M)
                END DO
                EG(I,IATOM) = TMP
              END DO
              GO TO 500
  520       CONTINUE
  510     CONTINUE
  500   CONTINUE
      END IF
C
C     ----- SYMMETRIZE THE FORCE CONSTANT MATRIX -----
C
      DO IATOM = 1,NAT
        DO JATOM = 1,IATOM
          IF(ITRANS(IATOM,JATOM).GT.1) THEN
            N = 3*(ITRANS(IATOM,JATOM)-1)
            IX = 3*(MAPCTR(IATOM,ITRANS(IATOM,JATOM))-1)
            IY = 3*(MAPCTR(JATOM,ITRANS(IATOM,JATOM))-1)
            JX = 3*(IATOM-1)
            JY = 3*(JATOM-1)
            DO I = 1,3
              DO J = 1,3
                TEMP = ZERO
                DO K = 1,3
                  TEMP = TEMP + FCM(IX+I,IY+K)*PTR(K,N+J)
                END DO
                V(I,J) = TEMP
              END DO
            END DO
            DO I = 1,3
              DO J = 1,3
                TEMP = ZERO
                DO K = 1,3
                  TEMP = TEMP + PTR(K,N+I)*V(K,J)
                END DO
                FCM(JX+I,JY+J) = TEMP
                FCM(JY+J,JX+I) = TEMP
              END DO
            END DO
          END IF
        END DO
      END DO
C
C     DEBUG. PRINT OUT THE FULL FCM.
C     ------------------------------
C     WRITE(IW,*) "FULL FCM"
C     CALL FCMOUT(FCM,3*NAT)
C
C     SAVE GRADIENT, IF NEEDED.
C     -------------------------
      IF(.NOT.GOTEG) THEN
         CALL DAWRIT(IDAF,IODA,EG,3*NAT,3,0)
         GOTEG = .TRUE.
      END IF
C
C     SAVE FCM.
C     ---------
      CALL DAWRIT(IDAF,IODA,FCM,NC1*NC1,4,0)
      IF (MASWRK) WRITE(IW,9100)
      RETURN
C
C     ----- HIT TIME LIMIT -----
C
  800 CONTINUE
      IF (MASWRK) WRITE(IW,9110)
      CALL ABRT
      STOP
C
 8000 FORMAT('ENERGY RESTART DATA FOR RUNTYP=HESSIAN, METHOD=FULLNUM'/
     *       ' $VIB2  ')
 8010 FORMAT(8X,1X,I5,' E=',F20.10)
 8040 FORMAT(' $END   ')
 8050 FORMAT(A8,1X,I5,3X,F20.10)
C
 9000 FORMAT(/1X,'BEGINNING HESSIAN COMPUTATION BY TWICE-DIFFERENCING',
     *          ' THE ENERGY...'/
     * 1X,'THE NUMBER OF ENERGIES AT DISPLACED GEOMETRIES WILL BE  ',I8,
     *    ' AT MOST,'/
     * 1X,'BUT ADDITIONAL SYMMETRY CONSIDERATIONS TO BE MADE AT EACH'/
     * 1X,'DISPLACED GEOMETRY MAY REDUCE THIS VALUE SUBSTANTIALLY.'/
     * 1X,'A CHECK JOB WILL PRINT THE PRECISE NUMBER OF POINTS NEEDED.'/
     * 1X,'WITHOUT SYMMETRY, THE NUMBER OF ENERGIES NEEDED WOULD BE',I8)
 9030 FORMAT(/1X,'THIS JOB WILL COMPUTE A TOTAL OF',I8,
     *           ' ENERGIES AT DISPLACED GEOMETRIES.')
 9040 FORMAT(1X,'ERROR IN $VIB2 GROUP')
 9100 FORMAT(' ...... END OF NUMERICAL HESSIAN CALCULATION ......')
 9110 FORMAT(//1X,'*** INTERNAL TIME LIMIT REACHED ****'/
     *         1X,'THERE IS NOT ENOUGH TIME LEFT TO START THE NEXT',
     *            ' ENERGY DISPLACEMENT.'/
     *         1X,'THIS RUN MUST BE RESTARTED WITH A $VIB2 GROUP.'//)
 9200 FORMAT(/1X,'************************************************',
     *       /1X,'* CALCULATING ENERGY ',I5,' FOR NUMERIC HESSIAN *',
     *       /1X,'************************************************')
 9210 FORMAT(1X,'MOVING ATOM',I4,' COORDINATE',A2,' TO',F20.10,' A.U.')
 9220 FORMAT(1X,'MOVING ATOM',I4,' COORDINATE',A2,' TO',F20.10,', AND'/
     *       1X,'MOVING ATOM',I4,' COORDINATE',A2,' TO',F20.10,' A.U.')
      END
C*MODULE HESS    *DECK VROT2
      SUBROUTINE VROT2(VIB,ISYMOP,IATOM,JATOM,DIS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER(MXATM=2000,MXSH=5000)
      PARAMETER(ZERO=0.0D+00)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      DIMENSION VIB(3,2),DIS(3,2),TMP(3,2)
C
C       TRANSFORMS A SET TWO CARTESIAN DISPLACEMENTS RELATED TO
C       IATOM AND JATOM BY THE ISYMOP TRANSFORMATION MATRIX
C
C     TRANSFORM THE DISPLACEMENT VECTORS BY ISYMOP
C
      N = 3*(ISYMOP-1)
      DO I = 1,3
        TMP1 = ZERO
        TMP2 = ZERO
        DO J = 1,3
          TMP1 = TMP1 + PTR(I,N+J)*VIB(J,1)
          TMP2 = TMP2 + PTR(I,N+J)*VIB(J,2)
        END DO
        TMP(I,1) = TMP1
        TMP(I,2) = TMP2
      END DO
C
C---  WRITE(6,8000)VIB(1,1),VIB(2,1),VIB(3,1),VIB(1,2),VIB(2,2),VIB(3,2)
C---  WRITE(6,8001)TMP(1,1),TMP(2,1),TMP(3,1),TMP(1,2),TMP(2,2),TMP(3,2)
C
C         DISPLACE THE MOLECULAR COORDINATES BY TVIB
C
      IATOM_T = MAPCTR(IATOM,ISYMOP)
      JATOM_T = MAPCTR(JATOM,ISYMOP)
      J1 = 1
      J2 = 2
      IF(IATOM .NE. IATOM_T) THEN
        J1 = 2
        J2 = 1
      END IF
      DO I = 1,3
        DIS(I,J1) = TMP(I,1) + C(I,IATOM_T)
        DIS(I,J2) = TMP(I,2) + C(I,JATOM_T)
      END DO
C
C---  WRITE(6,8002) C(1,IATOM),C(2,IATOM),C(3,IATOM),
C--- *             C(1,JATOM),C(2,JATOM),C(3,JATOM)
C---  WRITE(6,8003)DIS(1,1),DIS(2,1),DIS(3,1),DIS(1,2),DIS(2,2),DIS(3,2)
      RETURN
C8000 FORMAT(1X,'INPUT VIB ',3F6.3,' | ',3F6.3)
C8001 FORMAT(1X,'SYMOP*VIB ',3F6.3,' | ',3F6.3)
C8002 FORMAT(1X,'ORIG COOR ',3F6.3,' | ',3F6.3)
C8003 FORMAT(1X,'DISP COOR ',3F6.3,' | ',3F6.3)
      END
C*MODULE HESS    *DECK VECCMP
      SUBROUTINE VECCMP(A,B,N,TOL,TEST)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER(ZERO=0.0D+00)
      LOGICAL TEST
      DIMENSION A(N),B(N)
C
C          COMPARES TWO VECTORS (A & B) OF LENGTH N
C          RETURNS TEST = .TRUE. IF THE SUM((A-B)^2) < TOL
C
      TMP2 = ZERO
      DO I = 1,N
        TMP1 = A(I) - B(I)
        TMP2 = TMP2 + TMP1*TMP1
      END DO
C
      TEST = TMP2.LT.TOL
      RETURN
      END
C*MODULE HESS    *DECK NUMNRG
      SUBROUTINE NUMNRG(ESP,NSP,EVAL,WRK,MCCI,SCFOK,IFCM)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ESP(*),WRK(*)
      LOGICAL GOPARR,DSKWRK,MASWRK,MCCI,SCFOK
C
      PARAMETER (MXATM=2000)
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
      NSP  = NSP + 1
      EVAL = ZERO
C
      IF(EXETYP.EQ.CHECK) THEN
        IF(MASWRK) WRITE(IFCM,1000) NSP,ZERO
        RETURN
      END IF
C
      IF(ESP(NSP).EQ.ZERO) THEN
        CALL SYMOFF
        CALL FFSVMO(1,WRK)
        IF(MCCI) CALL FFSVCI(1)
        CALL ENERGX
        EVAL     = E
        ESP(NSP) = E
        IF(E.EQ.ZERO) SCFOK = .FALSE.
        CALL SYMON
      ELSE
        EVAL = ESP(NSP)
        IF(MASWRK) WRITE(6,9100) NSP,EVAL
      END IF
C
      IF(MASWRK) THEN
         WRITE(IFCM,1000) NSP,EVAL
         CALL FLSHBF(IFCM)
      END IF
      RETURN
C
 1000 FORMAT(8X,1X,I5,' E=',F20.10)
 9100 FORMAT(/1X,'************************************************',
     *       /1X,'* ENERGY ',I5,' E=',F20.10,10X,'*',
     *       /1X,'************************************************')
      END
C*MODULE HESS    *DECK POLOUT
      SUBROUTINE POLOUT(APOL)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION APOL(6),APOLANG(6)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (THREE=3.0D+00, UNIT=0.52917724924D+00)
C
C     ----- CONVERT UNITS FROM BOHRS TO ANGSTROMS ----
C
      DO 100 I=1,6
         APOLANG(I)=APOL(I)*UNIT**THREE
  100 CONTINUE
C
C     ---- CALCULATE THE MEAN POLARIZABILITY ----
C
      ALPHAM=(APOLANG(1)+APOLANG(3)+APOLANG(6))/THREE
C
C     ----- PRINT ALPHA POLARIZABILITY MATRIX -----
C
      IF (MASWRK) THEN
         WRITE(IW,9040)
         WRITE(IW,9060)  APOLANG(1)
         WRITE(IW,9080) (APOLANG(I),I=2,3)
         WRITE(IW,9100) (APOLANG(I),I=4,6)
         WRITE(IW,9110) ALPHAM
      END IF
      RETURN
C
 9040 FORMAT(/10X,42(1H-)/
     *        10X,'ALPHA POLARIZABILITY TENSOR (ANGSTROMS**3)'/
     *        10X,42(1H-)//
     *        23X,'X',14X,'Y',14X,'Z')
 9060 FORMAT(10X,'X ',3F15.5)
 9080 FORMAT(10X,'Y ',3F15.5)
 9100 FORMAT(10X,'Z ',3F15.5)
 9110 FORMAT(/10X,     'MEAN ALPHA POLARIZABILITY = ',F15.5)
      END
C*MODULE HESS    *DECK PURADM
      SUBROUTINE PURADM(ADMXYZ,ADMINT,BMAT,BINV,WRK,NCOORD,NVAR)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION ADMXYZ(6,NCOORD),BMAT(NVAR,NCOORD),WRK(*),
     *          ADMINT(6,NVAR),BINV(NCOORD,NVAR)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ----- PURIFY THE ALPHA POLARIZABILITY DERIVATIVE TENSOR -----
C
      NQ1 = NCOORD
      NQ3 = NVAR*NQ1
      CALL DAREAD(IDAF,IODA,ADMXYZ,6*NCOORD,249,0)
C
C     ----- TRANSFORM TO INTERNALS WITH B-INVERSE -----
C
      CALL DAREAD(IDAF,IODA,BINV,NQ3,50,0)
      CALL TFTNSR(ADMXYZ,ADMINT,BINV,WRK,NQ1,NVAR,6)
C
C     ----- AND BACK TO CARTESIANS WITH THE B MATRIX -----
C
      CALL DAREAD(IDAF,IODA,BMAT,NQ3,49,0)
      CALL TFTNSR(ADMINT,ADMXYZ,BMAT,WRK,NVAR,NQ1,6)
C
      CALL DAWRIT(IDAF,IODA,ADMXYZ,6*NCOORD,249,0)
      RETURN
      END
C*MODULE HESS    *DECK PURDDM
      SUBROUTINE PURDDM(DDM,DDN,BMAT,BINV,WRK,NCOORD,NVAR)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION DDM(3,NCOORD),BMAT(NVAR,NCOORD),WRK(*),
     *          DDN(3,NVAR),BINV(NCOORD,NVAR)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
C     ----- PURIFY THE DIPOLE DERIVATIVE TENSOR -----
C
      NQ1 = NCOORD
      NQ3 = NVAR*NQ1
      CALL DAREAD(IDAF,IODA,DDM,3*NCOORD,34,0)
C
C     ----- TRANSFORM TO INTERNALS WITH B-INVERSE -----
C
      CALL DAREAD(IDAF,IODA,BINV,NQ3,50,0)
      CALL TFTNSR(DDM,DDN,BINV,WRK,NQ1,NVAR,3)
C
C     ----- AND BACK TO CARTESIANS WITH THE B MATRIX -----
C
      CALL DAREAD(IDAF,IODA,BMAT,NQ3,49,0)
      CALL TFTNSR(DDN,DDM,BMAT,WRK,NVAR,NQ1,3)
C
      CALL DAWRIT(IDAF,IODA,DDM,3*NCOORD,34,0)
      RETURN
      END
C*MODULE HESS    *DECK PURFCM
      SUBROUTINE PURFCM(FCM,FCMINT,IZMAT,NCOORD,NROW,NVAR,
     *                  NZMAT,PURIFY,PRTIFC,DECOMP)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL PURIFY,PRTIFC,DECOMP,GOPARR,DSKWRK,MASWRK
C
      DIMENSION FCM(NCOORD,NCOORD),FCMINT(NROW,NROW),IZMAT(NZMAT)
C
C     ----- "HT2MDA" IS THE CONVERSION FACTOR FROM HARTREE/BOHR**2
C     ----- TO MILLIDYNE-ANGSTROM.
C     ----- HT2MDA = (4.3598E-08 ERG/HARTREE) *
C     -----         [(1 BOHR/0.529177 ANGSTROM)**2] *
C     -----         (1E+08 ANGSTROM/CM)
C     -----
C     ----- THE CONVERSIONS FROM HARTREE/BOHR-RADIAN TO MDYN/RADIAN
C     ----- AND HARTREE/RADIAN**2 TO MDYN-ANG/RADIAN**2 ARE
C     ----- HT2MDA*ANGSTR AND HT2MDA*ANGSTR**2, RESPECTIVELY.
C
      PARAMETER (HT2MDA=1.556914D+01, ANGSTR=0.52917724924D+00)
      PARAMETER (MXSIZE=30)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- PURIFY THE FORCE CONSTANT MATRIX -----
C     TRANSFORM HESSIAN TO INTERNALS, AND BACK TO CARTESIANS
C
      CALL DAREAD(IDAF,IODA,FCM,NCOORD*NCOORD,4,0)
      CALL TFHF(FCM,FCMINT,NVAR,NCOORD,NROW,NCOORD,50)
      NQ3 = NVAR*NVAR
      IF(DECOMP) NQ3=NROW*NROW
      CALL DAWRIT(IDAF,IODA,FCMINT,NQ3,43,0)
      IF(PRTIFC) THEN
         IF (MASWRK) WRITE(IW,9010)
         IF (MASWRK .AND. NCOORD.GT.MXSIZE)
     +       WRITE(IW,9030) HT2MDA,HT2MDA*ANGSTR*ANGSTR,HT2MDA*ANGSTR
         CALL PRSQ(FCMINT,NVAR,NVAR,NROW)
         CALL PUSQL(FCMINT,NVAR,NVAR,NROW)
C
C        ----- IF THERE ARE LESS THAN 30 COORDINATES, CONVERT TO
C        ----- UNITS OF MILLIDYNES/ANG, MDYN-ANG/RAD, AND MDYN/RAD**2
C        ----- AND PRINT INTERNAL COORDINATE HESSIAN AGAIN.
C
         IF (NCOORD.LE.MXSIZE) THEN
C
C           ----- BEFORE CONVERTING UNITS, FINISH THE PURIFICATION
C           ----- AND STORE ON THE DAF.
C
            IF (PURIFY) THEN
               CALL TFHB(FCMINT,FCM,NVAR,NCOORD,NROW,NCOORD,49)
               CALL DAWRIT(IDAF,IODA,FCM,NCOORD*NCOORD,4,0)
            END IF
C
            IF (MASWRK) WRITE(IW,9020)
C
C           ----- READ IN THE INTERNAL COORDINATE DEFINITIONS
C
            CALL DAREAD(IDAF,IODA,IZMAT,NZMAT,48,1)
C
C           ----- COMPRESS IZMAT SO THAT IT CONTAINS ONLY THE
C           ----- INTERNAL/SYMMETRY COORDINATE TYPES.
C
            IADD = 1
            IVAR = 1
  100       CONTINUE
            ITYPE = IZMAT(IADD)
            IF (ITYPE.LT.0) ISKIP=1
            IF (ITYPE.EQ.1) ISKIP=3
            IF (ITYPE.EQ.2) ISKIP=4
            IF (ITYPE.EQ.3) ISKIP=5
            IF (ITYPE.EQ.4) ISKIP=5
            IF (ITYPE.EQ.5) ISKIP=4
            IF (ITYPE.EQ.6) ISKIP=6
            IF (ITYPE.EQ.7) ISKIP=7
            IADD = IADD + ISKIP
            IZMAT(IVAR) = ITYPE
            IF (ITYPE.EQ.5) THEN
               IVAR = IVAR + 1
               IZMAT(IVAR) = ITYPE
            END IF
            IVAR = IVAR + 1
            IF (IVAR .LE. NVAR) GO TO 100
C                               *********
C
            DO 210 J=1,NVAR
               FACT1 = HT2MDA
               IF (ABS(IZMAT(J)).GT.1) FACT1=FACT1*ANGSTR
               DO 200 I=1,NVAR
                  FACT2 = FACT1
                  IF (ABS(IZMAT(I)).GT.1) FACT2 = FACT2*ANGSTR
                  FCMINT(I,J) = FCMINT(I,J)*FACT2
  200          CONTINUE
  210       CONTINUE
            CALL PRSQ(FCMINT,NVAR,NVAR,NROW)
C
C           ---- RESTORE THE ORIGINAL UNITS
C
            CALL DAREAD(IDAF,IODA,FCMINT,NQ3,43,0)
         END IF
      END IF
      IF(.NOT.PURIFY) RETURN
      CALL TFHB(FCMINT,FCM,NVAR,NCOORD,NROW,NCOORD,49)
      CALL DAWRIT(IDAF,IODA,FCM,NCOORD*NCOORD,4,0)
      IF (MASWRK) WRITE(IW,9000)
      RETURN
C
 9000 FORMAT(1X,'THE HESSIAN MATRIX HAS BEEN PURIFIED')
 9010 FORMAT(/1X,'HESSIAN MATRIX IN INTERNAL COORDINATES'/
     *       1X,'UNITS ARE HARTREE/BOHR**2, HARTREE/RADIAN**2, OR',
     *       1X,'HARTREE/BOHR*RADIAN')
 9020 FORMAT(/1X,'HESSIAN MATRIX IN INTERNAL COORDINATES'/
     *       1X,'UNITS ARE MILLIDYNE/ANGSTROM, MDYN-ANG/RADIAN**2,',
     *       1X,'OR MDYN/RADIAN')
 9030 FORMAT(/1X,'TO CONVERT UNITS, USE THE FOLLOWING CONVERSION ',
     *           'FACTORS'/
     *        5X,'MILLIDYNE/ANGSTROM',11X,'= ',F10.6,
     *           ' * HARTREE/BOHR**2'/
     *        5X,'MILLIDYNE-ANGSTROM/RADIAN**2 = ',F10.6,' * HARTREE/',
     *           'RADIAN**2'/
     *        5X,'MILLIDYNE/RADIAN ',12X,'= ',F10.6,
     *           ' * HARTREE/BOHR*RADIAN')
      END
C*MODULE HESS    *DECK PUZMAT
      SUBROUTINE PUZMAT(ZMAT,IZMAT,IWRK,NVAR,NZMAT)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION ZMAT(NVAR),IZMAT(NZMAT),IWRK(NVAR)
      LOGICAL GOPARR,DSKWRK,MASWRK
      PARAMETER (MXATM=2000)
      COMMON /FUNCT / E,EGRAD(3,MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C        PUNCH INTERNAL COORDINATE INFORMATION FOR PROGRAM -ESTTS-
C          (HESSIAN IS PUNCHED LATER BY -PURFCM-)
C
      CALL DAREAD(IDAF,IODA,ZMAT ,NVAR ,39,0)
      CALL DAREAD(IDAF,IODA,IZMAT,NZMAT,40,1)
C
      IF(MASWRK) WRITE(IP,8000) E
C
C          PUNCH COORDINATE VALUES AND TYPES
C
      II = 1
      NDEG = 0
      DO 100 I=1,NVAR
         ITYP = IZMAT(II)
         IWRK(I) = ITYP
         IF(ITYP.LE.0  .OR.  ITYP.GE.8) THEN
            IF(MASWRK) WRITE(IW,*) 'PUZMAT: ILLEGAL ITYP'
            CALL ABRT
         END IF
         IF(ITYP.EQ.1) NCHK = 3
         IF(ITYP.EQ.2) NCHK = 4
         IF(ITYP.EQ.3) NCHK = 5
         IF(ITYP.EQ.4) NCHK = 5
         IF(ITYP.EQ.5) NCHK = 4
         IF(ITYP.EQ.5) NDEG = NDEG+1
         IF(ITYP.EQ.6) NCHK = 6
         IF(ITYP.EQ.7) NCHK = 7
         IF(ITYP.NE.5  .OR.  NDEG.EQ.2) THEN
            II = II + NCHK
            NDEG = 0
         END IF
  100 CONTINUE
C
      IF(MASWRK) WRITE(IP,8010) ZMAT
      IF(MASWRK) WRITE(IP,8020) IWRK
      RETURN
C
 8000 FORMAT('ENERGY, COORDINATES, TYPES, INTERNAL COORD HESSIAN',
     *       ' (ALL ARE A.U.):'/F20.10)
 8010 FORMAT(5F15.7,:,' >')
 8020 FORMAT(15I5,:,' >')
      END
C*MODULE HESS    *DECK RDVIB
      SUBROUTINE RDVIB(IR,IW,IRST,NCOORD,IVIB,IATOM,ICOORD,E,EG,DIP,EOF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      CHARACTER*8 WORD,ENDWRD
      CHARACTER*80 LINE
C
      PARAMETER (MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      DIMENSION EG(NCOORD),DIP(3),JCON(2)
C
      LOGICAL EOF,GOPARR,DSKWRK,MASWRK
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      SAVE JVIB,JATOM,JCOORD,E0
      DATA ENDWRD/' $END   '/
C
C     ----- READ RESTART DATA FROM UNIT IR -----
C
      IZERO = 0
      IONE  = 1
      ITWO  = 2
      IF (IRST.EQ.0) THEN
         JVIB = 0
         JATOM = 0
         JCOORD = 0
         EOF=.FALSE.
         E0=ZERO
      ELSE
         IF(MASWRK) E=E0
         GO TO 100
      END IF
C
C        JCON(1)=0,1,2 MEANS EOF,$END,VIB FOUND
C
      IF(MASWRK) THEN
        READ(IR,FMT='(A80)',END=900) LINE
        READ(LINE,9000,ERR=800) WORD,JVIB,JATOM,JCOORD,E
        CALL UPRCAS(WORD,8)
        IF(WORD.EQ.ENDWRD) THEN
           EOF=.TRUE.
        ELSE IF (NFRG.NE.0.AND.IVIB.EQ.2.AND.JVIB.EQ.1) THEN
           EOF=.TRUE.
           E0=E
           IRST=1
        END IF
        JCON(1)=2
        IF(EOF) JCON(1)=1
        JCON(2)=IRST
        IF (GOPARR) CALL DDI_BCAST(475,'I',JCON,2,MASTER)
        IF(EOF) RETURN
      ELSE
        IF (GOPARR) THEN
           CALL DDI_BCAST(475,'I',JCON,2,MASTER)
           IRST=JCON(2)
           IF (JCON(1).EQ.IZERO) GO TO 800
           IF (JCON(1).EQ.IONE) THEN
              EOF=.TRUE.
              RETURN
           END IF
        END IF
      END IF
C
C        JCON(1)=0,1,2 MEANS EOF, MIXED UP DATA, GOOD DATA FOUND
C
  100 IRST=0
      IF(MASWRK) THEN
         IF(JVIB.NE.IVIB .OR. JATOM.NE.IATOM .OR. JCOORD.NE.ICOORD) THEN
            WRITE(IW,9040)
     *         IVIB,IATOM,ICOORD,JVIB,JATOM,JCOORD
            IF (GOPARR) CALL DDI_BCAST(476,'I',IONE,1,MASTER)
            CALL ABRT
         END IF
C
         WRITE(IW,9060) JVIB,JATOM,JCOORD,E
         IF(NCOORD.GT.0) READ(IR,9080,END=900) EG
         IF(NFRG.GT.0) THEN
            DO 620 IFR=1,NFRG
               READ(IR,9080,END=900) (DEFT(IFRC,IFR),IFRC=1,3)
               READ(IR,9080,END=900) (TORQ(IFRC,IFR),IFRC=1,3)
 620        CONTINUE
         END IF
         READ(IR,9080,END=900) DIP
         IF (GOPARR) CALL DDI_BCAST(476,'I',ITWO,1,MASTER)
      ELSE
         IF (GOPARR) THEN
            CALL DDI_BCAST(476,'I',JCON,1,MASTER)
            IF (JCON(1).EQ.IZERO) GO TO 900
            IF (JCON(1).EQ.IONE) CALL ABRT
         END IF
      END IF
C
      IF (GOPARR) THEN
                         CALL DDI_BCAST(477,'F',E   ,1     ,MASTER)
                         CALL DDI_BCAST(478,'F',DIP ,3     ,MASTER)
         IF(NCOORD.GT.0) CALL DDI_BCAST(479,'F',EG  ,NCOORD,MASTER)
         IF(NFRG.GT.0)   CALL DDI_BCAST(478,'F',DEFT,3*NFRG,MASTER)
         IF(NFRG.GT.0)   CALL DDI_BCAST(477,'F',TORQ,3*NFRG,MASTER)
      END IF
      RETURN
C
  800 CONTINUE
      IF (MASWRK) THEN
         IF (GOPARR) CALL DDI_BCAST(476,'I',IZERO,1,MASTER)
         WRITE(IW,9030)
      END IF
      CALL ABRT
C
  900 CONTINUE
      IF (MASWRK) THEN
         IF (GOPARR) CALL DDI_BCAST(476,'I',IZERO,1,MASTER)
         WRITE(IW,9020)
      END IF
      CALL ABRT
      RETURN
C
 9000 FORMAT(A8,6X,I4,7X,I4,8X,I4,3X,F20.10)
 9020 FORMAT(1X,'**** ERROR, UNEXPECTED END OF FILE READING $VIB GROUP')
 9030 FORMAT(1X,'**** ERROR, READING $VIB GROUP.')
 9040 FORMAT(1X,'**** ERROR, EXPECTING IVIB,IATOM,ICOORD=',3I10/
     *       20X,'GOT JVIB,JATOM,JCOORD=',3I10/
     *       1X,'INSPECT YOUR $VIB GROUP.')
 9060 FORMAT(1X,'IVIB=',I4,' IATOM=',I4,' ICOORD=',I4,' E=',F20.10)
 9080 FORMAT(5E16.9)
      END
C*MODULE HESS    *DECK RDVIB2
      SUBROUTINE RDVIB2(IOREAD,IOWRIT,IMODE,IVIB,EVIB2)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ZERO=0.0D+00)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      EVIB2 = ZERO
      IMDERR = 0
      IVBERR = 0
      IF (MASWRK) THEN
         READ(IOREAD,END=800,ERR=900,IOSTAT=IBOMBD,FMT=8000)
     *        KMODE,KVIB,EVIB2
C
C     --- VERIFY THAT THE VALUES KMODE AND KVIB READ FROM $VIB2
C     --- MATCH THE EXPECTED DISPLACEMENT DEFINED BY IMODE AND IVIB
C
         IF (KMODE .NE. IMODE) IMDERR=1
         IF (KVIB  .NE. IVIB ) IVBERR=1
         IF (IMDERR+IVBERR .EQ. 0) THEN
            WRITE(IW,9000) KMODE,KVIB,EVIB2
            WRITE(IOWRIT,9010) KMODE,KVIB,EVIB2
            CALL FLSHBF(IOWRIT)
         END IF
      END IF
      IF (GOPARR) CALL DDI_BCAST(481,'F',EVIB2,1,MASTER)
C
  800 CONTINUE
C
C     ******
      RETURN
C     ******
C
C     --- ERROR SECTION
C
  900 CONTINUE
      IF (MASWRK) THEN
         WRITE(IW,9100)
         IF (IMDERR .GT. 0) WRITE(IW,9110) 'IMODE',IMODE,KMODE
         IF (IVBERR .GT. 0) WRITE(IW,9110) 'IVIB ', IVIB,KVIB
         IF (IBOMBD .NE. 0) WRITE(IW,9120) IOREAD,IBOMBD
      END IF
C
C     *********
      CALL ABRT
C     *********
C
 8000 FORMAT(1X,7X,I4,6X,I2,3X,F20.10)
 9000 FORMAT(1X,'FOUND IMODE= ',I4,', IVIB= ',I1,', E= ',F20.10)
 9010 FORMAT(1X,'IMODE= ',I4,' IVIB= ',I1,' E= ',F20.10)
 9100 FORMAT(1X,'*** ERROR(S) IN SUBROUTINE RDVIB2 ****')
 9110 FORMAT(1X,'RDVIB2:  ILLEGAL VALUE OF ',A5,' FOUND IN $VIB2.',
     *          ' EXPECTED ',I5,' BUT FOUND ',I5)
 9120 FORMAT(1X,'RDVIB2:  I/O ERROR WHILE READING UNIT ',I5,
     *          '.  IOSTAT ERROR CODE = ',I5)
      END
C*MODULE HESS    *DECK SETDDM
      SUBROUTINE SETDDM(DDM,DIP,DEL,NP,IVIB,NCOORD,NVIB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (ONE=1.0D+00)
C
      DIMENSION DDM(3,NCOORD),DEQ(3),DIP(3)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      SAVE DEQ
C
      DELI = ONE/DEL
      NCORD3 = 3 * NCOORD
      IF(NP .GT. 0) GO TO 100
C
C        INITIALIZE DIPOLE DERIVATIVE MATRIX
C
         DEQ(1) = DIP(1)
         DEQ(2) = DIP(2)
         DEQ(3) = DIP(3)
         CALL VCLR(DDM,1,NCORD3)
         CALL DAWRIT(IDAF,IODA,DDM,NCORD3,34,0)
         RETURN
C
C     UPDATE DIPOLE DERIVATIVE MATRIX
C
  100 CONTINUE
      CALL DAREAD(IDAF,IODA,DDM,NCORD3,34,0)
      IF(IVIB .EQ. 2) GO TO 125
      IF(NVIB .EQ. 2) GO TO 115
      DO 110 I=1,3
         DDM(I,NP) = (DIP(I) - DEQ(I))*DELI
  110 CONTINUE
      GO TO 135
C
  115 CONTINUE
      DO 120 I=1,3
         DDM(I,NP) = DIP(I)
  120 CONTINUE
      GO TO 135
C
  125 CONTINUE
      DO 130 I=1,3
         DDM(I,NP) = (DDM(I,NP) - DIP(I))*DELI
  130 CONTINUE
C
  135 CONTINUE
      CALL DAWRIT(IDAF,IODA,DDM,NCORD3,34,0)
      RETURN
      END
C*MODULE HESS    *DECK SETFCM
      SUBROUTINE SETFCM(FCM,M,NV,EG,IVIB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION FCM(M,M),EG(*)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
C
C     ----- SET ELEMENTS OF THE FORCE CONSTANT MATRIX -----
C
      NCF=6*NFRG
      MNF=M-NCF
C
      IF(IVIB .GT. 0) GO TO 200
C
C                   INITIALIZE FCM TO ZERO
C
      CALL VCLR(FCM,1,M*M)
      CALL DAWRIT(IDAF,IODA,FCM,M*M,4,0)
      RETURN
C
C                   STORE COLUMN 'NV' IN THE FCM
C
  200 CONTINUE
      CALL DAREAD(IDAF,IODA,FCM,M*M,4,0)
      IF(IVIB.GT.1) GO TO 300
C                   FIRST DIFFERENCING
         IFR=0
         IFRPT=0
      DO 220 I = 1,M
         IF(I.LE.MNF) THEN
          FCM(I,NV) = EG(I)
         ELSE
           IFC=MOD(IFRPT,6)
           IF(IFC.EQ.0)IFR=IFR+1
           IF(IFC.LE.2) THEN
              IX=IFC+1
              FCM(I,NV)=DEFT(IX,IFR)
           ELSE
              IX=IFC-2
              FCM(I,NV)=TORQ(IX,IFR)
           END IF
             IFRPT=IFRPT+1
         END IF
  220 CONTINUE
      CALL DAWRIT(IDAF,IODA,FCM,M*M,4,0)
      RETURN
C     ******
C                   SECOND DIFFERENCING
  300 CONTINUE
         IFR=0
         IFRPT=0
      DO 320 I = 1,M
         IF(I.LE.MNF) THEN
         FCM(I,NV) = FCM(I,NV)-EG(I)
         ELSE
             IFC=MOD(IFRPT,6)
             IF(IFC.EQ.0)IFR=IFR+1
           IF(IFC.LE.2) THEN
              IX=IFC+1
              FCM(I,NV)=FCM(I,NV)-DEFT(IX,IFR)
           ELSE
              IX=IFC-2
              FCM(I,NV)=FCM(I,NV)-TORQ(IX,IFR)
           END IF
             IFRPT=IFRPT+1
         END IF
  320 CONTINUE
      CALL DAWRIT(IDAF,IODA,FCM,M*M,4,0)
      RETURN
      END
C*MODULE HESS    *DECK SYMDDM
      SUBROUTINE SYMDDM(DDM,SKIP,NCOORD,NAT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION DDM(3,NCOORD)
      DIMENSION V(3,3)
C
      LOGICAL SKIP(*),OUT,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      PARAMETER (ZERO=0.0D+00, TM9=1.0D-09)
C
C     ----- SYMMETRIZE THE DIPOLE DERIVATIVE TENSOR -----
C
      OUT = NPRINT.EQ.1 .AND. MASWRK
C
C     ----- OBTAIN RAW DIPOLE DERIVATIVE MATRIX -----
C
      CALL DAREAD(IDAF,IODA,DDM,3*NCOORD,34,0)
      IF(OUT) WRITE(IW,9020)((DDM(I,J),I=1,3),J=1,NCOORD)
C
C     ----- GET TRANFORMATION MATRICES OF COORDINATES. -----
C
      IF(NT.EQ.1) GO TO 400
      CALL DAREAD(IDAF,IODA,PTR,432,7,0)
C
C     ---- SEE WHICH BLOCKS HAVE ALREADY BEEN COMPUTED ----
C
      DO 220 IC = 1,NAT
         SKIP(IC) = .FALSE.
         DO 210 IT = 1,NT
            IF (MAPCTR(IC,IT) .GT. IC) GO TO 220
  210    CONTINUE
         SKIP(IC) = .TRUE.
  220 CONTINUE
C
C     ----- SYMMETRIZE DIPOLE DERIVATIVES -----
C
      DO 360 IC = 1,NAT
         IF(.NOT.SKIP(IC)) GO TO 360
         DO 340 IT = 1,NT
            ICNU = MAPCTR(IC,IT)
            IF(SKIP(ICNU)) GO TO 340
            INDX = 3*(IC-1)
            N = 3*(IT-1)
            DO 280 JT = 1,3
               DO 260 L=1,3
                  DUM = ZERO
                  DO 240 M=1,3
                     DUM = DUM + DDM(JT,INDX+M)*PTR(L,N+M)
  240             CONTINUE
                  V(JT,L) = DUM
  260          CONTINUE
  280       CONTINUE
            INDX = 3*(ICNU-1)
            DO 320 JT = 1,3
               DO 310 L=1,3
                  DUM = ZERO
                  DO 300 M=1,3
                     DUM = DUM + PTR(JT,N+M)*V(M,L)
  300             CONTINUE
                  DDM(JT,INDX+L) = DUM
  310          CONTINUE
  320       CONTINUE
            SKIP(ICNU) = .TRUE.
  340    CONTINUE
  360 CONTINUE
C
  400 CONTINUE
      DO 420 I=1,NAT
         DO 410 J=1,3
            IF(ABS(DDM(J,I)).LT.TM9) DDM(J,I)=ZERO
  410    CONTINUE
  420 CONTINUE
C
C     ----- STORE IT ON THE DAF -----
C
      CALL DAWRIT(IDAF,IODA,DDM,3*NCOORD,34,0)
      RETURN
C
 9020 FORMAT(/' RAW DIPOLE DERIVATIVE MATRIX'/(3E15.6))
      END
C*MODULE HESS    *DECK SYMFCM
      SUBROUTINE SYMFCM(FCM,SKIP2,EG0,NCOORD,NCF,NAT,VIBSIZ,NVIB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SKIP2(NAT,NAT)
C
      PARAMETER (MXFRG=1050, MXSH=5000, MXATM=2000)
C
      DIMENSION FCM(NCOORD+NCF,NCOORD+NCF),EG0(NCOORD)
      DIMENSION O(3,3),V(3)
C
      COMMON /FGROL / DEFTO(3,MXFRG),TORQO(3,MXFRG)
      COMMON /FROZNC/ IFREEZ(3*MXATM),NFRZ
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMSPD/ PTR(3,144),DTR(6,288),FTR(10,480),GTR(15,720)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
      PARAMETER (ZERO=0.0D+00, HALF=0.5D+00, ONE=1.0D+00)
C
      DATA HSSIAN/8HHESSIAN /
C
C     ----- COMPLETE COMPUTATION OF THE FORCE CONSTANT MATRIX -----
C     SYMMETRY DUPLICATE BLOCKS ARE FORMED, MATRIX IS SYMMETRIZED.
C
      NC1 = NCOORD
      NC3 = (NCOORD+NCF)*(NCOORD+NCF)
C
C     ----- COMPLETE THE FINITE DIFFERENCING -----
C
      CALL DAREAD(IDAF,IODA,FCM,NC3,4,0)
      IF (NVIB .EQ. 2) GO TO 120
C
      IF(NAT.GT.0) CALL DAREAD(IDAF,IODA,EG0,NC1,29,0)
      DO 110 J = 1,NCOORD+NCF
         IFR=0
         IFRPT=0
         DO 100 I = 1,NCOORD+NCF
          IF(I.LE.NCOORD) THEN
            FCM(I,J) = FCM(I,J)-EG0(I)
          ELSE
            IFC=MOD(IFRPT,6)
            IF(IFC.EQ.0)IFR=IFR+1
            IF(IFC.LE.2) THEN
              IX=IFC+1
              FCM(I,J)=FCM(I,J)-DEFTO(IX,IFR)
            ELSE
              IX=IFC-2
              FCM(I,J)=FCM(I,J)-TORQO(IX,IFR)
            END IF
            IFRPT=IFRPT+1
          END IF
  100    CONTINUE
  110 CONTINUE
C
  120 CONTINUE
      DUM = ONE/(VIBSIZ*NVIB)
      DO 140 I = 1,NCOORD+NCF
         DO 130 J = 1,NCOORD+NCF
            FCM(I,J) = DUM*FCM(I,J)
  130    CONTINUE
  140 CONTINUE
C
C     ----- GET TRANFORMATION MATRICES OF COORDINATES. -----
C
      IF(NT.EQ.1) GO TO 570
      CALL DAREAD(IDAF,IODA,PTR,432,7,0)
C
C     ----- CHECK WHICH BLOCKS OF THE -FCM- HAVE BEEN COMPUTED -----
C
      DO 220 JAT = 1,NAT
         DO 210 IAT = 1,NAT
           SKIP2(IAT,JAT) = .FALSE.
  210    CONTINUE
  220 CONTINUE
      DO 280 JAT = 1,NAT
         DO 240 IT = 1,NT
            IF (MAPCTR(JAT,IT) .GT. JAT) GO TO 280
  240    CONTINUE
         DO 260 IAT = 1,NAT
            SKIP2(IAT,JAT) = .TRUE.
  260    CONTINUE
  280 CONTINUE
C
C     ----- GENERATE SYMMETRY EQUIVALENT BLOCKS -----
C
      DO 560 IAT = 1,NAT
      DO 540 JAT = 1,NAT
      IF (SKIP2(IAT,JAT)) GO TO 540
      DO 520 IT = 1,NT
      KAT = MAPCTR(IAT,IT)
      LAT = MAPCTR(JAT,IT)
      IF ( .NOT. SKIP2(KAT,LAT)) GO TO 520
      N = 3*(IT-1)
      LOCI = 3*(IAT-1)
      LOCJ = 3*(JAT-1)
      LOCK = 3*(KAT-1)
      LOCL = 3*(LAT-1)
      DO 360 K = 1,3
      DO 360 L = 1,3
  360 O(K,L) = FCM(LOCK+K,LOCL+L)
      DO 420 K = 1,3
      DO 400 L = 1,3
      DUM = ZERO
      DO 380 M = 1,3
  380 DUM = DUM+O(K,M)*PTR(M,N+L)
  400 V(L) = DUM
      DO 420 L = 1,3
  420 O(K,L) = V(L)
      DO 480 L = 1,3
      DO 460 K = 1,3
      DUM = ZERO
      DO 440 M = 1,3
  440 DUM = DUM+PTR(M,N+K)*O(M,L)
  460 V(K) = DUM
      DO 480 K = 1,3
  480 O(K,L) = V(K)
      DO 500 K = 1,3
      DO 500 L = 1,3
  500 FCM(LOCI+K,LOCJ+L) = O(K,L)
      SKIP2(IAT,JAT) = .TRUE.
      GO TO 540
  520 CONTINUE
  540 CONTINUE
  560 CONTINUE
C
C     ----- SYMMETRIZE THE FORCE CONSTANT MATRIX -----
C
  570 CONTINUE
      DO 590 I = 2,NCOORD+NCF
         DO 580 J = 1,I-1
            AVE = (FCM(I,J)+FCM(J,I))*HALF
            FCM(J,I) = AVE
            FCM(I,J) = AVE
  580    CONTINUE
  590 CONTINUE
C
C     ----- ZERO OFF FCM ELEMENTS FOR PARTIAL HESSIAN ANALYSIS -----
C
      DO 600 I = 1, NCOORD
        DO 610 J = 1, NCOORD
          IF(IFREEZ(J).NE.0)FCM(I,IFREEZ(J))=ZERO
          IF(IFREEZ(J).NE.0)FCM(IFREEZ(J),I)=ZERO
  610   CONTINUE
  600 CONTINUE
C
      IF((IFREEZ(1).NE.0) .AND. (RUNTYP.EQ.HSSIAN) ) THEN
         DO 620 J = 1, NCOORD+NCF
            DO 630 I = NCOORD+1, NCOORD+NCF
              FCM(I,J) = ZERO
              FCM(J,I) = ZERO
  630       CONTINUE
  620    CONTINUE
      END IF
C
C     ---- SET DIAGONAL FCM ELEMENTS BE 1.0D-08 FOR PARTIAL HESSIAN ----
C
      DO 650 I = 1, NCOORD
        IF(IFREEZ(I).NE.0)FCM(IFREEZ(I),IFREEZ(I))=1.0D-08
  650 CONTINUE
C
      IF((IFREEZ(1).NE.0) .AND. (RUNTYP.EQ.HSSIAN) ) THEN
        DO 660 I = NCOORD+1, NCOORD+NCF
          FCM(I,I) = 1.0D-08
  660   CONTINUE
      END IF
C
C
C
C     ----- WRITE CARTESIAN SPACE HESSIAN TO THE DAF -----
C
      CALL DAWRIT(IDAF,IODA,FCM,NC3,4,0)
      RETURN
      END
C*MODULE HESS    *DECK TFTNSR
      SUBROUTINE TFTNSR(TNSR1,TNSR2,Q,DUM,N,M,LDIM)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION TNSR1(LDIM,N),TNSR2(LDIM,M),Q(N,M),DUM(M)
C
C     TRANSFORM TNSR2 = TNSR1 * Q
C
      INCX = LDIM
      DO 30 I=1,INCX
         DO 10 K=1,M
            DUM(K) = DDOT(N,TNSR1(I,1),INCX,Q(1,K),1)
   10    CONTINUE
         DO 20 J=1,M
            TNSR2(I,J) = DUM(J)
   20    CONTINUE
   30 CONTINUE
      RETURN
      END
C
C*MODULE HESS    *DECK VALFOR
      SUBROUTINE VALFOR(NCOORD,IVIB,IATOM,ICOORD,NPRT,NPUN,SCFOK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,SCFOK
C
      PARAMETER (MXATM=2000, MXAO=8192, MXPT=2000, MXFRG=1050,
     *           MXFGPT=12000, MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      CHARACTER*6 FRGNAM
      CHARACTER*8 FRGNME,PTNAM
C
      LOGICAL MMONLY,QMMM
C
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBTRM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FMCOM / X(1)
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
      COMMON /FRAGMT/ XCRD(3,MXFRG),YCRD(3,MXFRG),ZCRD(3,MXFRG),
     *                PTNAM(3,MXFRG),FRGNAM(MXFRG)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      COMMON /FUNCT / E,EG(3*MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /TMVALS/ TI,TX,TIM
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, TEENY=1.0D-08)
C
      DATA CHECK /8HCHECK   /
      DATA UHF/8HUHF     /
      DATA GUGA   /8HGUGA    /
      DATA NONE/4HNONE/
C
C     ----- CALCULATE ENERGY + DIPOLE MOMENTS + GRADIENT
C               FOR VIBRATIONAL ANALYSIS -----
C
C         PRINT THE DISPLACEMENT WE ARE WORKING ON NOW
C
      IF (MASWRK) THEN
      WRITE(IW,9008) IATOM,ICOORD,IVIB
         DO 100 I = 1,NAT
            WRITE(IW,9028) I,ANAM(I),BNAM(I),C(1,I),C(2,I),C(3,I)
  100    CONTINUE
C
         IF(NFRG.GT.0) THEN
            WRITE(IW,9011)
            WRITE(IW,9021)
            IAT=0
            KMASS=0
            DO 120 IFRG=1,NFRG
               WRITE(IW,9022) FRGNAM(IFRG)
               DO 110 III=1,NMPTS(IFRG)
                  IAT=IAT+1
                  IF(FMASS(III+KMASS).GT.ZERO) THEN
                     XF = EFC(1,IAT)
                     YF = EFC(2,IAT)
                     ZF = EFC(3,IAT)
                     WRITE(IW,9031) FRGNME(IAT),XF,YF,ZF
                  END IF
  110         CONTINUE
              KMASS=KMASS+NMPTS(IFRG)
  120      CONTINUE
         END IF
      END IF
C
C     ----- MINIMIZE PRINT/PUNCH OUTPUT AT DISPLACED GEOMETRIES -----
C
      IF(IVIB.GT.0) THEN
         NPRINT = -5
                       NPUNCH = 0
         IF(NPUN.EQ.1) NPUNCH = 1
         IF(NPUN.EQ.2) NPUNCH = 2
         IF(NPUN.GT.0  .AND.  MASWRK) WRITE(IP,8010) IVIB,IATOM,ICOORD
      END IF
C
C     ----- CALCULATE ENERGY -----
C     IF THE ENERGY DOES NOT CONVERGE (E=0) FOR THE DISPLACED GEOMETRIES
C     THE CALCULATION IS ALLOWED TO PROCEED ANYWAY.  BUT A WARNING WILL
C     BE PRINTED SO THE USER KNOWS RESTART WITH FIXED $VIB IS NECESSARY.
C
      E = ZERO
      CALL VCLR(EG,1,NCOORD)
      CALL ENERGX
      IF(E.EQ.ZERO) SCFOK=.FALSE.
C
      IF(E.EQ.ZERO   .AND.  IVIB.EQ.0  .AND. EXETYP.NE.CHECK) THEN
         IF (MASWRK) WRITE(IW,9100)
         CALL ABRT
      END IF
C
      IF(IVIB.GT.0  .AND.  NPRT.EQ.1) THEN
         CALL VALFM(LOADFM)
         LVEC = LOADFM + 1
         LEIG = LVEC   + NUM*NUM
         LAST = LEIG   + NUM
         NEED = LAST - LOADFM -1
         CALL GETFM(NEED)
         IF(SCFTYP.EQ.UHF) THEN
            CALL DAREAD(IDAF,IODA,X(LVEC),NUM*NUM,15,0)
            CALL DAREAD(IDAF,IODA,X(LEIG),NUM    ,17,0)
            IF(MASWRK) WRITE(IW,*) ' '
            IF(MASWRK) WRITE(IW,*) 'ALPHA ORBITALS'
            IF(MASWRK) CALL PREV(X(LVEC),X(LEIG),NUM,NUM,NUM)
            CALL DAREAD(IDAF,IODA,X(LVEC),NUM*NUM,19,0)
            CALL DAREAD(IDAF,IODA,X(LEIG),NUM    ,21,0)
            IF(MASWRK) WRITE(IW,*) ' '
            IF(MASWRK) WRITE(IW,*) 'BETA ORBITALS'
            IF(MASWRK) CALL PREV(X(LVEC),X(LEIG),NUM,NUM,NUM)
         ELSE
            CALL DAREAD(IDAF,IODA,X(LVEC),NUM*NUM,15,0)
            CALL DAREAD(IDAF,IODA,X(LEIG),NUM    ,17,0)
            IF(MASWRK) WRITE(IW,*) ' '
            IF(MASWRK) WRITE(IW,*) 'MOLECULAR ORBITALS'
            IF(MASWRK) CALL PREV(X(LVEC),X(LEIG),NUM,NUM,NUM)
         END IF
         CALL RETFM(NEED)
      END IF
C
C     ----- CALCULATE GRADIENT -----
C
      IF (CITYP.EQ.GUGA) THEN
         CALL CIGRAD
      ELSE
         CALL HFGRAD
      END IF
C
C     ---- QMMM OPTIMIZATION ----
C     IF QMMM, DO MM OPT AND FORM HYBRID GRADIENT.
C
      IF(QMMM) CALL TOYS
C
      DO 160 I = 1,NCOORD
         IF(ABS(EG(I)).LT.TEENY) EG(I) = ZERO
  160 CONTINUE
      DO 180 I = 1,NFRG
         DO 170 J = 1,3
            IF(ABS(DEFT(J,I)).LT.TEENY) DEFT(J,I) = ZERO
            IF(ABS(TORQ(J,I)).LT.TEENY) TORQ(J,I) = ZERO
  170    CONTINUE
  180 CONTINUE
C
C     ----- CALCULATE DIPOLES (ENERGX DOESN'T ALWAYS DO VIB 0 DIPOLE)
C     THIS COMES AFTER THE GRADIENT SO THAT IN CASE OF CI OR MP2,
C     WE WILL HAVE THE CORRECT DENSITY MATRIX FOR THE PROPERTIES.
C
      IF(IVIB.GT.0  .OR.  (IVIB.EQ.0  .AND.  NPRINT.EQ.-5)) THEN
         IF(MPCTYP.EQ.NONE) THEN
            CALL ELMOMC
         ELSE
            CALL MPCPRP(.TRUE.,.FALSE.)
         END IF
      END IF
      IF (TIM .GE. TIMLIM) RETURN
C
C     ----- CHECK FOR STATIONARY POINT -----
C
      IF(EXETYP.EQ.CHECK) RETURN
      IF(IVIB.GT.0) RETURN
      CALL EGMAX(EG,3*NAT,DEFT,TORQ,NFRG,GMAX,GRMS)
      IF (MASWRK) THEN
         WRITE(IW,9040) GRMS
         WRITE(IP,9040) GRMS
      END IF
      RETURN
C
 8010 FORMAT('----- IVIB=',I2,'  IATOM=',I5,'  ICOORD=',I2)
C
 9008 FORMAT(1H1,7H   ATOM,I3/ 40X,19(1H-)/8H   COORD,I3,29X,
     *     19HNUCLEAR COORDINATES/ 40X,19(1H-)/8H     VIB,I3//22X,
     *     5H ATOM,13X, 1HX,14X,1HY,14X,1HZ/)
 9011 FORMAT(1X,'COORDINATES OF FRAGMENT MULTIPOLE CENTERS (BOHR)')
 9021 FORMAT(1X,'MULTIPOLE NAME',8X,'X',14X,'Y',14X,'Z'/
     *       1X,60(1H-))
 9022 FORMAT(1X,'FRAGNAME=',A6)
 9028 FORMAT(15X,I3,2X,A8,A2,3F15.6)
 9031 FORMAT(1X,A8,7X,3F15.10)
 9100 FORMAT(' *** NO FORCE FIELD, SCF DOES NOT CONVERGE AT VIB0 POINT')
 9040 FORMAT(1X,'THE VIB 0 GRADIENT RMS =',1P,E15.8)
      END
C*MODULE HESS    *DECK NUMGRDX
      SUBROUTINE NUMGRDX
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXRT=100, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION SVENGY(MXRT+15),ENGYIN(MXRT+15)
C
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /FMCOM / X(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP,LSTMPTS(MXFRG)
      COMMON /FUNCT / E,EG(3,MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /MASSES/ ZMASS(MXATM)
C
      CHARACTER*8  WORD,ENDWRD
      EQUIVALENCE (ENGYIN(1),ENUCR)
C
      LOGICAL FIRST,MCCI,RSTART
      DIMENSION ESP(5),A(5),PTS(5)
      PARAMETER (TOL=1.0D-04,ONE=1.0D+00,ZERO=0.0D+00)
C
      DATA  RMC /8HMCSCF   /
      DATA GUGA /8HGUGA    /
      DATA CHECK/8HCHECK   /
      DATA NONE/4HNONE/
      DATA ENDWRD/8H $END   /
      DATA FIRST/.TRUE./
C
      SAVE FIRST,N0VIB
C
      IZERO = 0
C
C        GENERATE GRADIENT BY NUMERICAL DIFFERENTIATION OF THE ENERGY
C        CODE WRITTEN BY RYAN OLSON, DECEMBER 2003
C        INSPIRATION FROM DMITRI FEDOROV
C
      IF(MASWRK) WRITE(IW,9000)
C
      IF(NFRG.GT.0) THEN
         IF(MASWRK) WRITE(IW,9030)
         CALL ABRT
      END IF
C
      MCCI   = SCFTYP.EQ.RMC  .OR.  CITYP.EQ.GUGA
      VIBSIZ = 0.01D+00
      IVIB   = 0
      NUM2 = (NUM*NUM+NUM)/2
C
C        DYNAMIC MEMORY ALLOCATION
C
      NCCF = 3*NAT
      CALL VALFM(LOADFM)
      LX   = LOADFM+ 1
      LGTOT= LX    + NCCF
      LRM  = LGTOT + NCCF
      LP   = LRM   + NCCF + 21*NFRG
      LGE  = LP    + NCCF*NCCF
      LDR  = LGE   + 3*MXATM
      LT   = LDR   + NCCF
      LWRK = LT    + NCCF*NCCF
      LTMP = LWRK  + MAX(NUM*NUM,NCCF)
      LAST = LTMP  + 6*NAT
      NEEDP= LAST-LOADFM-1
      CALL GETFM(NEEDP)
C
C     SAVE CURRENT ENERGY, ORBITALS, AND OTHER INFORMATION
C     ----------------------------------------------------
      E0 = E
      DO I=1,MXRT+15
         SVENGY(I) = ENGYIN(I)
      ENDDO
      CALL FFSVMO(IVIB,X(LWRK))
      IF(MCCI) CALL FFSVCI(IVIB)
      CALL DAREAD(IDAF,IODA,X(LWRK),NUM*NUM, 45,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK),NUM*NUM,249,0)
      CALL DAREAD(IDAF,IODA,X(LWRK),NUM2   , 11,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK),NUM2   ,277,0)
      CALL DAREAD(IDAF,IODA,X(LWRK),NUM2   , 12,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK),NUM2   ,278,0)
      IF(MPCTYP.EQ.NONE) THEN
         CALL DAREAD(IDAF,IODA,X(LWRK),NUM2, 13,0)
         CALL DAWRIT(IDAF,IODA,X(LWRK),NUM2,279,0)
      END IF
C
C     SAVE THE CURRENT PRINT VALUES
C     -----------------------------
      ISVPRT = NPRINT
      ISVPUN = NPUNCH
      NPRINT = -5
      NPUNCH = 0
C
C     FROM PRJGRD WE ONLY WANT THE P MATRIX (AND THE REDUCED MASSES).
C
      CALL VCLR(X(LGE),1,3*MXATM)
      CALL PRJGRD(X(LGE),C,X(LX),X(LGTOT),X(LRM),X(LP),ZMASS,NAT,3*NAT,
     *            NCCF)
C
C     UN-MASS-WEIGH
C
      DO I=1,NCCF
        CALL DSCAL(NCCF,ONE/X(LRM+I-1),X(LP+(I-1)*NCCF),1)
      END DO
C
C     PROJECT EACH ROW ON THE TOTALLY SYMMETRIC REPRESENTATION
C
      DO I=1,NCCF
        CALL DCOPY(NCCF,X(LP+I-1),NCCF,X(LDR),1)
        CALL SYMDR(X(LDR))
        RABS=SQRT(DDOT(NCCF,X(LDR),1,X(LDR),1))
C           NORMALISE
        IF(RABS.GT.TOL) CALL DSCAL(NCCF,ONE/RABS,X(LDR),1)
        CALL DCOPY(NCCF,X(LDR),1,X(LP+I-1),NCCF)
      END DO
C
C     THROW AWAY TRANSLATIONS, ROTATIONS AND NONSYMMETRIC STRETCHES.
C     STRETCHES ARE STORED AS RAWS, SO TRANSPOSE NOW AND BACK
C     (RAWS BECAUSE COORDINATES ARE TRANSFORMED X'=AX: RAWS OF A;
C     GRADIENT IS D/DX=D/DX' * DX'/DX = D/DX' * A).
C
      CALL TRPOSQ(X(LP),NCCF)
      CALL UNLINDEP(NCCF,NCCF,X(LP),NCCF,X(LT),NCCF,X(LWRK),TOL,NMODE)
      CALL TRPOSQ(X(LP),NCCF)
      IF(MASWRK) WRITE(IW,9040) NMODE,NCCF,2*NMODE
C-DBG-      CALL PRSQ(X(LP),NCCF,NMODE,NCCF)
C
C        SAVE COORDINATES
C
      CALL VCLR(EG,1,3*NAT)
      CALL DCOPY(3*NAT,C,1,X(LGE),1)
C
C     ----- OPEN SUPPLEMENTAL PUNCH FILE FOR RESTART DATA -----
C
      IFCM=4
      IF (MASWRK) CALL SEQOPN(IFCM,'IRCDATA','NEW',.FALSE.,'FORMATTED')
C
C     ----- LOOK FOR A RESTART FILE ON THE FIRST CALL TO THIS ROUTINE
C
      RSTART = .FALSE.
      IFCM=4
      IF (FIRST) THEN
         FIRST = .FALSE.
         N0VIB = 2
         CALL SEQREW(IR)
         CALL FNDGRP(IR,' $VIB2  ',IEOF)
         IF(IEOF.NE.1) THEN
            RSTART=.TRUE.
            IF(MASWRK) THEN
              WRITE(IW,*)
     *       'READING RESTART DATA FROM $VIB2 GROUP IN INPUT FILE...'
C
C       --- AT THIS STAGE, READ ONLY THE VALUE OF NVIB.
C       --- DISPLACEMENT DATA WILL BE READ LATER ON.
C
              READ(IR,8030) WORD,N0VIB
              CALL UPRCAS(WORD,8)
              IF (WORD.NE.ENDWRD) THEN
                 IF ((N0VIB .NE. 2) .AND. (N0VIB .NE. 5)) THEN
                    WRITE(IW,9140)
                    WRITE(IW,9150) N0VIB
                    CALL ABRT
                 END IF
              END IF
            END IF
            IF (GOPARR) CALL DDI_BCAST(480,'I',N0VIB,1,MASTER)
         END IF
      END IF
C
C     SET -NVIB- TO 5 TO PERFORM A MORE ACCURATE CALCULATION, WHICH
C     MAY BE USEFUL IN THE DEVELOPING ANALYTIC GRADIENT PROGRAMS.
C     NVIB=2 USES 2 POINTS (ON EITHER SIDE OF VIB, BUT NOT VIB0)
C     NVIB=5 USES 5 POINTS (TWO ON EITHER SIDE, AND INCLUDES VIB0)
C
C     JUST CHANGE THE CODE IF YOU DON'T MIND DOUBLING THE EXECUTION
C     TIME OF NUMERICAL GRADIENTS, OR USE JERRY'S TRICK BELOW:
C     --- NOTE:  THIS IS AN UNDOCUMENTED WAY TO USE NVIB=5 WITHOUT
C     ---        HAVING TO CHANGE THE CODE; I.E., ADD $VIB2 TO THE
C     ---        INPUT FILE TO TRIGGER A RESTART, AND THEN PUT
C     ---        "NVIB= 5" AS THE FIRST (AND ONLY!) LINE OF $VIB2
C
      NVIB = N0VIB
C
      NDISPL = 2 * NMODE
      IF (NVIB .EQ. 5) NDISPL = 4 * NMODE
      IF(EXETYP.EQ.CHECK) THEN
         IF(MASWRK) WRITE(IW,*)
     *       'NUMBER OF DISPLACEMENTS TO CALCULATE WITH NVIB= ',NVIB,
     *       ' = ',NDISPL
         CALL DAWRIT(IDAF,IODA,EG,3*NAT,3,0)
         CALL RETFM(NEEDP)
         RETURN
      END IF
C
      IOREAD = IR
      IOWRIT = IFCM
C
C     ----- OPEN SUPPLEMENTAL PUNCH FILE FOR RESTART DATA -----
C
      EVIB2 = ZERO
      IF (RSTART) CALL RDVIB2(IOREAD,IOWRIT,IZERO,IZERO,EVIB2)
      IF (EVIB2 .EQ. ZERO) THEN
         RSTART=.FALSE.
         IF (MASWRK) THEN
            WRITE(IFCM,8000)
            WRITE(IFCM,8010) NVIB
            WRITE(IFCM,8020) 0,0,E0
            CALL FLSHBF(IFCM)
         END IF
      END IF
C
C        START LOOP OVER DISPLACEMENTS
C
      IF(NVIB.EQ.2) THEN
         DO IMODE = 1,NMODE
           DO IVIB = 1,2
             DIFF = VIBSIZ
             IF(IVIB.EQ.2) DIFF = -VIBSIZ
C
             EVIB2 = ZERO
             IF (RSTART) CALL RDVIB2(IOREAD,IOWRIT,IMODE,IVIB,EVIB2)
             IF (EVIB2 .EQ. ZERO) THEN
                RSTART=.FALSE.
C
C    DISPLACE ALL COORDINATES ALONG THE CURRENT MODE
C
                CALL DAXPY(3*NAT,DIFF,X(LP+IMODE-1),NCCF,C,1)
C
                IF(MASWRK) WRITE(IW,9100) IMODE,IVIB
                IF(MASWRK) CALL PRATM(IW,1)
C
C     CALCULATE ENERGY AT DISPLACED MODE
C
                CALL FFSVMO(IVIB,X(LWRK))
                IF(MCCI) CALL FFSVCI(IVIB)
                CALL ENERGX
                IF(E.EQ.ZERO) THEN
                   IF(MASWRK) WRITE(IW,9020)
                   CALL ABRT
                END IF
                ESP(IVIB) = E
                IF (MASWRK) THEN
                   WRITE(IFCM,8020) IMODE,IVIB,E
                   CALL FLSHBF(IFCM)
                END IF
C
C      RESTORE ORIGINAL COORDINATES
C
                CALL DCOPY(3*NAT,X(LGE),1,C,1)
             ELSE
                ESP(IVIB) = EVIB2
             END IF
           END DO
           X(LTMP+IMODE-1) = (ESP(1)-ESP(2))/(2.0D+00*VIBSIZ)
         END DO
C
C        --- CURRENT NUMERICAL GRADIENT STEP IS COMPLETE.
C        --- CLOSE/DELETE THE IRCDATA FILE SO A NEW ONE
C        --- CAN BE CREATED FOR THE NEXT GRADIENT.
C
         IF (MASWRK) CALL SEQCLO(IFCM,'DELETE')
      ELSE
         DO IMODE = 1,NMODE
           DO 50 IVIB = 1,5
             PTS(IVIB) = (IVIB-3)*VIBSIZ
             ESP(IVIB) = E0
             IF(IVIB.EQ.3) GO TO 50
             EVIB2 = ZERO
             IF (RSTART) CALL RDVIB2(IOREAD,IOWRIT,IMODE,IVIB,EVIB2)
             IF (EVIB2 .EQ. ZERO) THEN
                RSTART=.FALSE.
C
C    DISPLACE ALL COORDINATES ALONG THE CURRENT MODE
C
                CALL DAXPY(3*NAT,PTS(IVIB),X(LP+IMODE-1),NCCF,C,1)
C
                IF(MASWRK) WRITE(IW,9100) IMODE,IVIB
                IF(MASWRK) CALL PRATM(IW,1)
C
C     CALCULATE ENERGY AT DISPLACED MODE
C
                CALL FFSVMO(IVIB,X(LWRK))
                IF(MCCI) CALL FFSVCI(IVIB)
                CALL ENERGX
                IF(E.EQ.ZERO) THEN
                   IF(MASWRK) WRITE(IW,9020)
                   CALL ABRT
                END IF
                ESP(IVIB) = E
                IF (MASWRK) THEN
                   WRITE(IFCM,8020) IMODE,IVIB,E
                   CALL FLSHBF(IFCM)
                END IF
C
C      RESTORE ORIGINAL COORDINATES
C
                 CALL DCOPY(3*NAT,X(LGE),1,C,1)
              ELSE
                 ESP(IVIB) = EVIB2
              END IF
   50      CONTINUE
           DUMMY = 0.0D+00
           CALL POLFIT(5,5,0,PTS,ESP,A,DUMMY,CHISQR)
           X(LTMP+IMODE-1) = A(2)
         END DO
C
C        --- CURRENT NUMERICAL GRADIENT STEP IS COMPLETE.
C        --- CLOSE/DELETE THE IRCDATA FILE SO A NEW ONE
C        --- CAN BE CREATED FOR THE NEXT GRADIENT.
C
         IF (MASWRK) CALL SEQCLO(IFCM,'DELETE')
      END IF
C
C     CONVERT THE GRADIENT BACK TO CARTESIANS
C
      CALL MRARBR(X(LTMP),1,1,NMODE,X(LP),NCCF,NAT*3,EG,1)
C
C     RECOVER THE E0 ENERGY, WAVEFUNCTION AND ET CETERA
C     -------------------------------------------------
      NPRINT = ISVPRT
      NPUNCH = ISVPUN
      E = E0
      DO I=1,MXRT+15
         ENGYIN(I) = SVENGY(I)
      ENDDO
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
      IVIB = 3
      CALL FFSVMO(IVIB,X(LWRK))
      IF(MCCI) CALL FFSVCI(IVIB)
      CALL DAREAD(IDAF,IODA,X(LWRK),NUM*NUM,249,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK),NUM*NUM, 45,0)
      CALL DAREAD(IDAF,IODA,X(LWRK),NUM2   ,277,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK),NUM2   , 11,0)
      CALL DAREAD(IDAF,IODA,X(LWRK),NUM2   ,278,0)
      CALL DAWRIT(IDAF,IODA,X(LWRK),NUM2   , 12,0)
      IF(MPCTYP.EQ.NONE) THEN
         CALL DAREAD(IDAF,IODA,X(LWRK),NUM2,279,0)
         CALL DAWRIT(IDAF,IODA,X(LWRK),NUM2, 13,0)
      END IF
C
      CALL DAWRIT(IDAF,IODA,EG,3*NAT,3,0)
      CALL RETFM(NEEDP)
      RETURN
C
 8000 FORMAT('ENERGY RESTART DATA FOR RUNTYP=OPTIMIZE, NUMGRD=.TRUE.'/
     *       ' $VIB2  ')
 8010 FORMAT(1X,'NVIB=',I2)
 8020 FORMAT(1X,'IMODE= ',I4,' IVIB= ',I1,' E= ',F20.10)
 8030 FORMAT(A6,I2)
 9000 FORMAT(/1X,'COMPUTING THE NUCLEAR GRADIENT BY NUMERIC',
     *           ' DIFFERENTIATION')
 9020 FORMAT(/1X,'ENERGY FAILED TO CONVERGE, NUMERICAL GRADIENT',
     *           ' CANNOT PROCEED')
 9030 FORMAT(/1X,'NUMERICAL GRADIENT IS NOT PROGRAMMED FOR THE EFP',
     *           ' MODEL')
 9040 FORMAT(1X,'PROJECTION FINDS',I4,
     *          ' TOTALLY SYMMETRIC MODES OUT OF 3N=',I6/
     *       1X,'SO',I6,' SINGLE POINT ENERGIES MUST BE COMPUTED TO',
     *          ' GET A GRADIENT.')
 9100 FORMAT(1X,'WORKING ON IMODE=',I4,' IVIB=',I1,'...')
 9140 FORMAT(1X,'ERROR IN $VIB2 GROUP')
 9150 FORMAT(1X,'ILLEGAL VALUE OF NVIB FOUND (',I5,'); MUST BE 2 OR 5.')
      END
