C 30 APR 07 - SPW,TDI - NEW MODULE FOR NEO HESSIANS
C
C*MODULE NEOHSS  *DECK HESSX_NEO
      SUBROUTINE HESSX_NEO(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     --- PROGRAM WRITTEN BY TZVETELIN IRODANOV AND SIMON P. WEBB ---
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(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(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
         IF(NEORUN.EQ.1) THEN
            CALL HSSNUM_NEO(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)
         ELSE
            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)
         END IF
         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
      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)
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 NEOHSS  *DECK HSSNUM_NEO
      SUBROUTINE HSSNUM_NEO(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 /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      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 NEO 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.
CSPW       --- THIS KIND OF STUFF LEFT IN ASSUMING NEO WILL GO EFP SOMETIME
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(IVIB.eq.0) call vibzro(.true.)
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_NEO(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_NEO(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     --- SKIP ANY QM NUCLEAR CENTER DISPLACEMENTS ---
C
      DO 305 I = 1,NQMNUC
          IAT = IQMNUC(I)
          SKIP(IAT) = .TRUE.
  305 CONTINUE
      DO 306 I = 1,NTAUXB
          IAT = IAUXNB(I)
          SKIP(IAT) = .TRUE.
  306 CONTINUE
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_NEO(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
             call vibzro(.false.)
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_NEO(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
C        --- SET GRADIENTS FOR SKIPPED QM NUCLEI GRADS TO ZERO ---
         DO 405 I = 1,NQMNUC
             IAT = IQMNUC(I)
             DO J=1,3
                EG((IAT-1)*3 + J) = 0.0D+00
             END DO
  405    CONTINUE
         DO 406 I = 1,NTAUXB
             IAT = IAUXNB(I)
             DO J=1,3
                EG((IAT-1)*3 + J) = 0.0D+00
             END DO
  406    CONTINUE
C
         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
c        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_NEO(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 NEOHSS  *DECK VALFOR_NEO
      SUBROUTINE VALFOR_NEO(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)
C
      IF(IVIB.EQ.0) THEN
C        --- FOR VIB0 JUST CALCULATE NEO ENERGY HERE
         CALL ENERGX
      ELSE
C        --- FOR CLASSICAL NUCLEUS DISPLACEMENT OPTIMIZE ALL QM-NUC CENETRS ---
         CALL NEOSPT
      END IF
C
      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
      CALL UNFREZ(.TRUE.)
      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 NEOHSS  *DECK VIBZRO
      SUBROUTINE VIBZRO(VIBWRT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL VIBWRT
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)
C
C     WRITE(*,*) 'IN VIBZRO NAT=',NAT
      IF (VIBWRT) THEN
         CALL DAWRIT(IDAF,IODA,C,NAT*3,462,0)
      ELSE
         CALL DAREAD(IDAF,IODA,C,NAT*3,462,0)
      END IF
C
C     CALL RETFM(NEED)
C
      RETURN
      END
C--C*MODULE NEOHSS  *DECK FFSVMO_NEO
C--      SUBROUTINE FFSVMO_NEO(IVIB,WRK,WRK2)
C--C
C--      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C--C
C--      DIMENSION WRK(*),WRK2(*)
C--C
C--      PARAMETER (MXATM=2000)
C--C
C--      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
C--     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C--      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C--      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
C--     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
C--     *                NNA,NNB,NTAUXB
C--      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
C--     *                MPLEVL,MPCTYP
C--C
C--      DATA UHF /8HUHF     /
C--C
C--C     --- THIS MAY BE HOOKED UP AT SOME POINT - SPW ---
C--      L2 = (NUM*NUM+NUM)/2
C--      L3 = NUM*NUM
C--      N2 = (NUMNB*NUMNB+NUMNB)/2
C--      N3 = NUMNB*NUMNB
C--C
C--C     ----- WRITE VIB 0 MO'S AND DENSITY ON SPECIAL DAF FILES -----
C--C
C--      IF (IVIB .GT. 0) GO TO 100
C--C--TDI ELECTRONS
C--      CALL DAREAD(IDAF,IODA,WRK,L3,15,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L3,30,0)
C--      CALL DAREAD(IDAF,IODA,WRK,L2,16,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L2,32,0)
C--C--TDI QM NUCLEI
C--      CALL DAREAD(IDAF,IODA,WRK2,N3,444,0)
C--      CALL DAWRIT(IDAF,IODA,WRK2,N3,406,0)
C--      CALL DAREAD(IDAF,IODA,WRK2,N2,445,0)
C--      CALL DAWRIT(IDAF,IODA,WRK2,N2,508,0)
C--      IF (SCFTYP .NE. UHF) RETURN
C--      CALL DAREAD(IDAF,IODA,WRK,L3,19,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L3,31,0)
C--      CALL DAREAD(IDAF,IODA,WRK,L2,20,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L2,33,0)
C--      RETURN
C--C
C--C     ----- RECOVER VIB 0 MO'S AND DENSITY -----
C--C
C--  100 CONTINUE
C--C--TDI ELECTRONS
C--      CALL DAREAD(IDAF,IODA,WRK,L3,30,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L3,15,0)
C--      CALL DAREAD(IDAF,IODA,WRK,L2,32,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L2,16,0)
C--C--TDI QM NUCLEI
C--      CALL DAREAD(IDAF,IODA,WRK2,N3,406,0)
C--      CALL DAWRIT(IDAF,IODA,WRK2,N3,444,0)
C--      CALL DAREAD(IDAF,IODA,WRK2,N2,508,0)
C--      CALL DAWRIT(IDAF,IODA,WRK2,N2,445,0)
C--      IF (SCFTYP .NE. UHF) RETURN
C--      CALL DAREAD(IDAF,IODA,WRK,L3,31,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L3,19,0)
C--      CALL DAREAD(IDAF,IODA,WRK,L2,33,0)
C--      CALL DAWRIT(IDAF,IODA,WRK,L2,20,0)
C--      RETURN
C--      END
C*MODULE NEOHSS  *DECK SYMFCM_NEO
      SUBROUTINE SYMFCM_NEO(FCM,SKIP2,EG0,NCOORD,NCF,NAT,VIBSIZ,NVIB)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL SKIP2(NAT,NAT)
      LOGICAL SETZRO
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 /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      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,463,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
      I = 0
      DO IAT = 1,NAT
         SETZRO = .FALSE.
         DO IQN = 1,NQMNUC
            IF (IAT.EQ.IQMNUC(IQN)) THEN
               SETZRO = .TRUE.
            END IF
         END DO
         DO IAUX = 1,NTAUXB
            IF (IAT.EQ.IAUXNB(IAUX)) THEN
               SETZRO = .TRUE.
            END IF
         END DO
         DO J=1,3
            I=I+1
            IF(SETZRO) THEN
               DO K=1,NCOORD
                  FCM(I,K) = ZERO
                  FCM(K,I) = ZERO
               END DO
            END IF
         END DO
      END DO
C
C     ----- WRITE CARTESIAN SPACE HESSIAN TO THE DAF -----
C
      CALL DAWRIT(IDAF,IODA,FCM,NC3,4,0)
      RETURN
      END
C*MODULE NEOHSS  *DECK SETFCM_NEO
      SUBROUTINE SETFCM_NEO(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,463,0)
      RETURN
C
C                   STORE COLUMN 'NV' IN THE FCM
C
  200 CONTINUE
      CALL DAREAD(IDAF,IODA,FCM,M*M,463,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,463,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,463,0)
      RETURN
      END
