C 11 APR 08 - SAN - SAVE THERMOCHEMISTRY INFORMATION FOR G3(MP2)
C 24 MAR 07 - MWS - PAD FRGINF COMMON BLOCK
C 27 DEC 06 - BN  - THERMO: ADD POSSIBLE ANHARMONIC ZPE COMPUTATION
C 19 SEP 05 - BNJ - USE VSCF ANHARMONIC FREQS IN HARMONIC THERMOCHEM
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - DGF - PRJFCM: USE TFSQU INSTEAD OF TWO MULTIPLIES
C 14 MAR 05 - HL  - CHANGE FMASS IN FRGMSS COMMON
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK FRGINF
C 23 JUL 04 - GDF - THERMO: ALWAYS INCLUDE UNSCALED THERMOCHEMISTRY
C 19 MAY 04 - DGF - CHANGES TO ADD THE FMO METHOD
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  7 AUG 02 - MWS - PRJGRD: SHORT-CIRCUIT THE ROUTINE IF EFP'S USED
C  7 AUG 02 - JHJ - FGMTRX: COMPUTED REDUCED MASS
C 17 APR 02 - MWS - PRJGRD: CORRECTIONS FOR EFP PARTICLES
C 26 MAR 02 - MWS - ADD A ROUTINE FOR PROJECTING T+R FROM GRADIENTS.
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 19 NOV 00 - HL  - THERMO,FGMTRX: CHANGES FOR PARTIAL HESSIAN
C 26 OCT 00 - MWS - INTRODUCE MXAO PARAMETER
C  1 MAY 00 - MWS - FGMTRX: ADD RAMAN INTENSITY EVALUATION
C 16 FEB 00 - GK  - FIX ELECTRONIC ENTROPY
C 21 DEC 99 - MWS - MAKE SYMMOL COMMON CONSISTENT
C 25 OCT 99 - MWS - FGMTRX: DON'T READ DAF GRADIENT IF FRAG ONLY
C  6 JAN 98 - JAB - FGMTRX,DECVIB: ADD NSYMC TO /ZMAT/ COMMON
C                   DECVIB: REPLACE 0.89 SCALE FACTOR WITH SCLFAC
C                   FIX BUG IN SCALING OF INTRINSIC FORCE CONSTANTS
C  1 OCT 97 - MWS - FGMTRX: PRINT FRAGMENT NORMAL MODES FOR HUMANS
C  8 AUG 97 - MWS - THERMO: CORRECT THERMOCHEMISTRY WITH FRAGMENTS
C 16 JUL 97 - GNM - FGMTRX,THERMO: FRAGMENT ONLY CHANGES
C  9 JAN 96 - MWS - EXTEND FRGMSS COMMON
C 30 OCT 95 - MWS - FGMTRX: PRINT WARNING IF NUM. DISP. DIDN'T CONV.
C  1 FEB 95 - MWS - FGMTRX: CORRECT PROJECTION TO BE ONLY NON-STAT.
C 29 DEC 94 - MWS - ALLOW SKIP OF PROJECTION OF ROTATIONS
C 15 DEC 94 - JAB - DECVIB: CORRECT FORCE CONSTANT SCALING
C 17 NOV 94 - WC  - FIX VIBRATIONAL ANALYSIS OF EFFECTIVE FRAGMENTS
C 10 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 27 MAY 94 - PND - FGMTRX: EFFECTIVE FRAGMENT CHANGES
C 10 MAR 94 - FRJ - PRJFC: MOVE GRADIENT NORMALIZATION
C  9 DEC 93 - MWS - DECVIB: CHANGE DAF RECORD NUMBERS
C 19 OCT 93 - MWS - PRJFC,PRJFCM: CLEAN UP SOME, ALLOW FRAG DIM.
C  6 AUG 93 - MWS - ADD ATORQ TO /FGRAD/
C  1 APR 93 - PND - FGMTRX,RAMS: CHANGES FOR EFFECTIVE FRAGMENTS
C 29 MAR 93 - TLW - CHANGES FOR RECTANGULAR B AND B INVERSE
C 24 MAR 93 - MWS - SKIP DECOMPOSITION DURING CHECK RUNS
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C 24 FEB 92 - MWS - THERMO: CHANGE TO USE OF SYMMOL COMMON
C 10 JAN 92 - TLW - PRJFCM: CALL ABRT BEFORE STOP
C  8 JAN 92 - TLW - MAKE WRITES PARALLEL; ADD COMMON PAR
C 29 JUN 91 - JAB - THERMO: REMOVE IFMAX COMPUTATION FOR LINEAR CASE
C 29 APR 91 - MWS - THERMO: ZPE SHOULD INCLUDE THE SCALE FACTOR
C  9 APR 91 - MR  - FGMTRX: NPRINT=6 GIVES NARROW NORMAL MODE OUTPUT
C 11 FEB 91 - MWS - CHANGES IN JAN'S THERMO OPTIONS DUE TO MY STUPIDITY
C 28 NOV 90 - JHJ - THERMO:ADDED TEMP AND SCLFAC OPTION
C  3 NOV 90 - MWS - DUMMY ATOMS HAVE ZERO MASS, NEW OPTION IN RAMS.
C 17 AUG 90 - TRC - RAMS: EXTEND MASS TABLE UP TO LAWRENCIUM
C  7 MAR 90 - MWS - PUNCH NORMAL MODES SUITABLE FOR MOLPLT PROGRAM
C 15 DEC 89 - MWS - FORMAL ARG TO THERMO IS NFIRST, NOT NFIRTS
C  9 DEC 89 - MWS - FGMTRX: CORRECT NUCLEAR CHARGE FOR ECP CASE
C 22 SEP 89 - MWS - CHANGES ANGSTROMS-BOHR CONVERSION FACTOR
C 18 AUG 89 - JAB - CALCULATE INTRINSIC FORCE CONSTANTS IN DECVIB
C 02 MAY 89 - JAB - MOVE DETECTION OF TRANS/ROT MODES TO FGMTRX
C 01 MAY 89 - JAB - PRINT OUT NORMAL MODES IN INTERNALS IN DECVIB
C 10 MAR 89 - TLW - CORRECT THE MASSES OF HE,NA,P,RH,AG IN RAMS
C 21 FEB 89 - STE - DELETE FORMAT 9088 FROM DECVIB
C 18 FEB 89 - MWS - PROJECTION ROUTINE GETS GRADIENT FROM DAF FILE 3.
C  8 FEB 89 - MWS - FIX BUG WITH RECIPROCAL MASSES
C  2 FEB 89 - MWS - NEW ROUTINE INRTIA CREATED, PASS MASSES BY /MASSES/
C 17 JAN 89 - MWS - OLD MODULE -FORCE- SPLIT INTO -VIBANL- AND -HESS-,
C                   ALL HISTORY COMMENTS FROM OLD MODULE KEPT HERE
C 26 OCT 88 - MWS - MAKE NROW=1 IN FORCX WHEN NO Z-MATRIX GIVEN
C 22 AUG 88 - JAB - ADD VCLR CALL IN DECVIB TO ZERO EM MATRIX
C 12 AUG 88 - JAB - WRITE HESSIAN IN INTERNALS TO DAF EVEN IF PURIFY=.F.
C 10 AUG 88 - MWS - MXSH,MXGSH,MXGTOT FROM 120,10,440 TO 1000,30,5000
C 07 AUG 88 - JAB - ADD VIBRATIONAL DECOMPOSITION ROUTINE (DECVIB)
C 25 JUL 88 - MWS - ADD THERMODYNAMICS ROUTINE
C 30 MAY 88 - MWS - USE PARAMETER TO DIMENSION COMMONS
C 13 MAY 88 - MWS - WRITE GRADIENT IN FCM IF STATPT AND $GRAD OMITTED
C  8 MAY 88 - MWS - UP AO-S FROM 255 TO 2047
C  2 MAY 88 - MWS - CALL GUESMO DURING RESTART JOBS
C 12 FEB 88 - MWS - FIX BUG WITH LAB IN FGMTRX
C 16 JAN 88 - STE - FIX PUVIB CALL IN VALFOR
C 15 JAN 88 - MWS - ADD DIPOLES TO $VIB GROUP, ADD PROJECTION OF
C                   F.C.M. ACCORDING TO SHIRO KOSEKI'S ADAPTATION OF
C                   (RENAMED) ROUTINES PRJFC AND PRJFCM FROM 'POLYRATE'.
C 19 NOV 87 - STE - CLEAN UP SOME AND CLARIFY STRUCTURE; USE EXETYP
C 20 OCT 87 - JAM - ADD ROUTINES TFDDI, TFDDMP
C 15 OCT 87 - JAM - CLEAN UP, ADD ROUTINE SETTRA, ADD ROUTINES
C                   DDMIN AND PUDDM TO READ AND PUNCH DIPOLE DERIVS
C  6 OCT 87 - JAM - ADD ROUTINES SETDDM, SYMDDM TO COMPUTE DIPOLE
C                   DERIVATIVES AND TRANSFORM TO NORMAL COORDINATES
C 27 JUN 87 - MWS - PRINT/PUNCH THE VIB 0 GRADIENT
C 13 MAY 87 - MWS - ADD OPTION TO AVOID PURIFICATION
C 11 MAY 87 - STE - FORCX: RESTORE DEFINITION OF VIBWRD
C 11 MAY 87 - MWS - OPTIONALLY PRINT HESSIAN IN INTERNAL COORDS,
C                   PUNCH CORRECT HESSIAN MATRIX (NO MORE ASYMMETRY)
C  5 MAY 87 - STE - FORCX,RDVIB: CLARIFY VARIABLE USAGE, ESP. NTSAVE
C 30 APR 87 - STE - FGMTRX,FORCX: CHANGE IA,SCR FOR USE BY GLDIAG
C 10 APR 87 - MWS - MOVE ROUTINE PUVIB TO IOLIB
C 12 FEB 87 - MWS - RESTART FROM CARDS, ADDING PUVIB,RDVIB. REMOVE
C                   MSREAD,MSWRIT DAF RESTART, GENERAL CLEANUP
C  5 NOV 86 - STE - TIDY UP, USE PARAMETERS, ALLOW FCM WITH N>99
C                   PUHESS CALLS PUSQL NOW; REMOVE NAV REFERENCES
C  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
C 25 JUN 86 - LAM - INITIALIZE VARIABLES IN FORCIN
C 20 JUN 86 - MWS - FIX MASS WARNING MESSAGE IN RAMS
C 17 JUN 86 - MWS - CHANGE CONVERSION FACTOR TO CM**-1 IN FGMTRX,
C                   ADD ROUTINE FORCIN, PURIFY HESSIAN IN
C                   ROUTINE SYMFCM, COMMENT OUT CARTESIAN
C                   SPACE DIAGONALIZATION IN FGMTRX.
C                   JERRY'S VIB. ANALYSIS NOT INCLUDED.
C 31 OCT 85 - STE - GENERIC SQRT,ABS; FORCX: BETTER RESTART
C 13 AUG 85 - MWS - USE MOST ABUNDANT ISOTOPE FOR MASSES
C 22 AUG 84 - STE - REMOVE UNUSED MSREAD/WRITE PARM FROM FORCX
C 20 AUG 84 - MWS - ABORT FORCE CALC. IF SCF NO CONV. IN VALFOR
C  4 MAY 84 - STE - ERROR CONTROL IN RAMS, DICT.RECORD 2 IN VALFOR
C 18 JAN 84 - STE - /SCFOPT/ IN VALFOR MISSING MCONV
C 19 NOV 83 - STE - MAKE 4 CHARACTER HOLLERITH VARIABLES TYPE REAL
C 17 JUN 83 - MWS - PLACE KEN DYKEMA'S SAYVETZ CONDITIONS ONLINE
C  1 JUN 83 - MWS - OPTION TO RETURN MASSES IN RAMS
C 13 MAR 83 - MWS - MODULE GRDDRV SPLIT TO OPTX AND FORCX,
C                   ALLOW INPUT OF ISOTOPES, AND INPUT OF THE
C                   FORCE CONSTANT MATRIX, FIX BUG WITH NORMAL
C                   MODES, MAKE 600K OF STATIC ARRAYS DYNAMIC
C
C*MODULE VIBANL  *DECK FGMTRX
      SUBROUTINE FGMTRX(VEC,FCM,E,SCR,IA,RM,SVTZT,
     *                  SVTZR,SVTZTT,SVTZRT,CC,COM,ZMS,BUF1,BUF2,
     *                  NC1,NC2,NPART,DDM,GOTDDM,ADM,GOTADM,
     *                  PROJCT,STATPT,SCFOK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOTDDM,GOTADM,PROJCT,STATPT,DECOMP,LINEAR,PRTSCN,GOTFRQ,
     *        GOPARR,MASWRK,DSKWRK,EFLDL,PRJROT,PRJGRD,SCFOK,FGONLY,
     *        RTYPRT
C
C        THE FOLLOWING ARE HOLLERITH (1-4 CHARACTER) VARIABLES
      REAL IBLANK,LETI,CLAB(3),LAB(9)
C
      DIMENSION VEC(NC1,NC1),FCM(NC2),E(NC1),SCR(NC1,8),IA(NC1),
     *          RM(*),SVTZT(3,NC1),SVTZR(3,NC1),
     *          SVTZTT(NC1),SVTZRT(NC1),
     *          CC(3,NPART),COM(3,NPART),ZMS(NPART),
     *          BUF1(NC1,NC1),BUF2(NC1,NC1)
C         NOTE THAT NEXT TWO NEED NOT BE DIMENSIONED
C         UNLESS -GOTDM- AND/OR -GOTADM- ARE TRUE
      DIMENSION DDM(3,NC1),ADM(6,NC1)
      DIMENSION CMASS(3),STEP(6)
C
      CHARACTER*8 FRGNME
C
      PARAMETER (MXATM=2000, MXAO=8192, MXPT=2000, MXFRG=1050,
     *           MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /ANHVIB/ ZMASST,NIMAG,NLAST,NSKIP
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /EFLDC / EVEC(3),EFLDL
      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 /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 /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 /MASSES/ ZMASS(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
      COMMON /THERMD/ FREQ(3*MXATM),TEMP(10),SCLFAC,NTEMP,PRTSCN,GOTFRQ
      COMMON /VBDCMP/ DECOMP,NROW,NDEG
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00,
     *           THREE=3.0D+00, FOUR=4.0D+00, SIX=6.0D+00,
     *           SEVEN=7.0D+00, TFACT=2.642461D+07)
C
C        TFACT IS E*E*NA/4*PI*PI*C*C*A0*A0*A0,
C            E =4.803242E-10 SQRT(G*CM**3)/SEC
C            NA=6.022045E+23 /G
C            A0=5.291771E-09 CM
C            C =2.997925E+10 CM/SEC
C
      DATA CHECK/8HCHECK   /
      DATA G3MP2 /8HG3MP2   /
      DATA CLAB /4H   X,4H   Y,4H   Z/
      DATA LETI,IBLANK/2H I,2H  /
      DATA HSSIAN/8HHESSIAN /, VSCF/8HVSCF    /
C
C     ----- WILSON -FG- MATRIX VIBRATIONAL ANALYSIS -----
C     THIS IS CARRIED OUT IN CARTESIAN COORDINATES, BY MASS-WEIGHTING
C     THE FORCE CONSTANT MATRIX.  THIS MATRIX IS THEN DIAGONALIZED,
C     GIVING THE NORMAL MODES AND FREQUENCIES.  SEE THE REFERENCE
C     W.D.GWINN   J.CHEM.PHYS.  55, 477-481 (1971).
C
      FGONLY = NUM.EQ.0  .AND.  NFRG.GT.0
      NCOORD=NAT*3
      NCF=NC1-NCOORD
      NIMAG = 0
      IF(GOTFRQ) GO TO 200
      IF (MASWRK) THEN
         WRITE(IW,9000)
         IF(.NOT.STATPT) WRITE(IW,9010)
         IF(.NOT.SCFOK)  WRITE(IW,9020)
      END IF
C
C     GRADIENT IS READ INTO ARRAY -SVTZT-.
C     FOR EF GRAD USE SVTZTT TEMPORALLY, EFGRDM WILL MOVE THEM TO SVTXT.
C     HESSIAN IS READ INTO ARRAY -VEC-.
C
      IF(NAT.GT.0) CALL DAREAD(IDAF,IODA,SVTZT,NCOORD,3,0)
      INDX=NCOORD+1
      DO 90 IFRG=1,NFRG
        CALL DCOPY(3,DEFT(1,IFRG),1,SVTZTT(INDX),1)
        CALL DCOPY(3,TORQ(1,IFRG),1,SVTZTT(INDX+3),1)
        INDX=INDX+6
   90 CONTINUE
      CALL DAREAD(IDAF,IODA,VEC  ,NC1*NC1,4,0)
C
C     ZMASS=MASSES, RM=TRIPLES OF INVERSE SQUARE ROOTS OF MASSES.
C
      IF (MASWRK .AND. .NOT. FGONLY) THEN
         WRITE(IW,9030)
         DO 100 IAT = 1,NAT
            WRITE(IW,9040) IAT,ANAM(IAT),BNAM(IAT),ZMASS(IAT)
  100    CONTINUE
      END IF
      I0=0
      DO 105 I=1,NAT
         IF(ZMASS(I).GT.ZERO) THEN
            DMY = ONE/SQRT(ZMASS(I))
         ELSE
            DMY = ONE
         END IF
         RM(I0+1) = DMY
         RM(I0+2) = DMY
         RM(I0+3) = DMY
         I0=I0+3
  105 CONTINUE
      CALL EFRMSS(NFRG,RM(NCOORD+1))
C
C     ----- GENERATE MASS WEIGHTED GRADIENT AND HESSIAN ----
C
      I0 = 1
      DO 109 I=1,NAT
        DO 108 J=1,3
           SVTZT(J,I) = RM(I0) * SVTZT(J,I)
  108   CONTINUE
        I0 = I0+3
  109 CONTINUE
      CALL EFGRDM(NCOORD,NCF,NFRG,SVTZTT,RM,SVTZT,0)
C
      CALL HESMAS(NCOORD,NCF,NFRG,VEC,RM,FCM,0)
C
C     ----- COMPUTE THE CENTER OF MASS -----
C     ----- GET CENTER OF MASS IN MASS-WEIGHTED CARTESIAN COORDS -----
C
  200 CONTINUE
      CALL CENMAS(NAT,NFRG,C,COM,ZMASST,CMASS,ZMASS)
C
      IF(GOTFRQ) GO TO 550
C
C     ----- PROJECT THE FORCE CONSTANT MATRIX -----
C
      IF (PROJCT) THEN
         PRJROT = .NOT.EFLDL
         PRJGRD = .NOT.STATPT
         CALL PRJFC(PRJGRD,PRJROT,ZMASST,VEC,SVTZT,COM,RM,FRGMI,
     *              BUF1,BUF2,NAT,NCOORD,NC1)
C
C        PUT FCM INTO SYMMETRIC STORAGE FOR DIAGONALIZATION CALL
C        IF NOT PROJCT, FCM DONE BY HESMAS EARLIER
C
         IJ = 0
         DO 320 I = 1,NC1
            DO 310 J = 1,I
               IJ = IJ+1
               FCM(IJ) = VEC(I,J)
  310       CONTINUE
  320    CONTINUE
      END IF
C
C     ----- GET NORMAL MODES AND FREQUENCIES -----
C
      IERR = 0
      CALL GLDIAG(NC1,NC1,NC1,FCM,SCR,E,VEC,IERR,IA)
      IF (IERR .NE. 0) CALL ABRT
      CALL STFASE(VEC,NC1,NC1,NC1)
C
C     ----- TRANSLATIONAL AND ROTATIONAL SAYVETZ CONDITIONS -----
C
      CALL VCLR(SVTZT,1,NC1*3)
      CALL VCLR(SVTZR,1,NC1*3)
      DO 460 I=1,NC1
         DO 420 J=1,NAT+NFRG
            JJ=MAX(3*(J-1),6*(J-1)-3*NAT)
            AMASS=ONE/RM((J-1)*3+1)
            DO 410 K=1,3
               K1=MOD(K+1,4)+(K+1)/4
               K2=MOD(K+2,4)+(K+2)/4
               SVTZT(K,I)=SVTZT(K,I)+AMASS*VEC(JJ+K,I)
               SVTZR(K,I)=SVTZR(K,I)+COM(K1,J)*VEC(JJ+K2,I)
     *                              -COM(K2,J)*VEC(JJ+K1,I)
  410       CONTINUE
  420    CONTINUE
C
         INDX=NCOORD+3
         JNDX=NCOORD+12*NFRG
         DO 450 J=1,NFRG
            DO 440 K=1,3
               DUM=ZERO
               DO 430 L=1,3
                  DUM=DUM+RM(JNDX+L)*VEC(INDX+L,I)
  430          CONTINUE
               SVTZR(K,I)=SVTZR(K,I)+DUM
               JNDX=JNDX+3
  440       CONTINUE
            INDX=INDX+6
  450    CONTINUE
C
         SVTZTT(I)=SQRT(SVTZT(1,I)**2+SVTZT(2,I)**2+SVTZT(3,I)**2)
         SVTZRT(I)=SQRT(SVTZR(1,I)**2+SVTZR(2,I)**2+SVTZR(3,I)**2)
  460 CONTINUE
C
C     ----- CONVERT NORMAL MODE DISPLACEMENTS -----
C     (FROM MASS WEIGHTED CARTESIAN TO CARTESIAN SPACE)
C
      DO 520 I = 1,NC1
         DO 510 J = 1,NCOORD
            VEC(J,I) = VEC(J,I)*RM(J)
  510    CONTINUE
         CALL EFGRDM(NCOORD,NCF,NFRG,VEC(1,I),RM,BUF1(1,1),0)
         CALL DCOPY(NCF,BUF1(NCOORD+1,1),1,VEC(NCOORD+1,I),1)
  520 CONTINUE
C
C             COMPUTE REDUCED MASS
C
      DO J = 1,NC1
         DD = DDOT(NC1,VEC(1,J),1,VEC(1,J),1)
         SCR(J,4) = ONE/DD
      END DO
C
C     ----- COUNT NUMBER OF NEGATIVE EIGENVALUES -----
C
      DO 530 I = 1,NC1
         IF(E(I).LT.ZERO) NIMAG=NIMAG+1
  530 CONTINUE
C
C     DECIDE WHICH MODES ARE NOT TRUE VIBRATIONS
C
  550 CONTINUE
      NSKIP = 6
      IF(LINEAR  .AND.  NFRG.EQ.0) NSKIP = 5
C
C        WE DO NOT SKIP TRANSLATION/ROTATION "MODES" IF THE
C        SET OF TRUE FREQUENCIES WAS READ IN.
C
      IF(GOTFRQ) THEN
         NFREQ = 3*NAT+6*NFRG - NSKIP
         DO 555 I=1,NFREQ
            IF(FREQ(I).LT.ZERO) NIMAG=NIMAG+1
 555     CONTINUE
         NFIRST = NIMAG+1
         NLAST = NIMAG
         GO TO 900
      END IF
C
      DO 560 IAT=1,NAT
         NUCZ = INT(ZAN(IAT)) + IZCORE(IAT)
         IF(NUCZ.EQ.0) NSKIP=NSKIP+3
  560 CONTINUE
      IF(PROJCT .AND. (.NOT. STATPT)) NSKIP = NSKIP + 1
C
      NLAST = NSKIP
      DO 570 I=1,NIMAG
         N2OF3 = 0
         MAYBE = I+NSKIP
         IF(MAYBE.GT.NC1) GO TO 580
         IF(ABS(E(I)).GT.ABS(E(MAYBE))) N2OF3=N2OF3+1
         IF(PROJCT) THEN
            IF(N2OF3.EQ.1) NLAST = NLAST + 1
         ELSE
            IF(0.01D+00+SVTZRT(I).LT.SVTZRT(MAYBE)) N2OF3=N2OF3+1
            IF(0.01D+00+SVTZTT(I).LT.SVTZTT(MAYBE)) N2OF3=N2OF3+1
            IF(N2OF3.GE.2) NLAST = NLAST + 1
            IF (MASWRK) THEN
               IF(N2OF3.EQ.2) WRITE(IW,9050) I,MAYBE
               IF(N2OF3.EQ.1) WRITE(IW,9050) MAYBE,I
            END IF
            IF(N2OF3.LE.1) GO TO 580
         END IF
  570 CONTINUE
C
  580 CONTINUE
      NFIRST = NLAST - NSKIP + 1
      NIMAG = NFIRST-1
      IF (MASWRK) WRITE(IW,9060) NFIRST,NLAST
C
C     ----- PRINT MESSAGE FOR PARTIAL HESSIAN ANALYSIS -----
C
      NFRZ = 0
      NCOORD=3*NAT
      DO 590 I=1, NCOORD
        IF(IFREEZ(I) .NE. 0) NFRZ = NFRZ + 1
  590 CONTINUE
C
      NFRZ = NFRZ + 6*NFRG
      RTYPRT = RUNTYP.EQ.HSSIAN .OR. RUNTYP.EQ.VSCF
      IF(RTYPRT  .AND.  NFRZ.GE.6.  AND.  MASWRK)THEN
          IF(NFRZ.GE.9)THEN
             WRITE(IW,9063)NLAST+1, NLAST+NFRZ-6, NLAST+NFRZ-5,
     *       NLAST+NFRZ-3, NLAST+1, NLAST+NFRZ-3
                  ELSE
             WRITE(IW,9064) NLAST+1, NLAST+NFRZ-3,
     *       NLAST+1, NLAST+NFRZ-3
          END IF
      END IF
C
C
C     ----- IF REQUESTED, DO VIBRATIONAL ANALYSIS
C
      IF (DECOMP) THEN
         NR3 = NROW*NROW
         CALL VALFM(LOADFM)
         I10 = 1 + LOADFM
         I15 = I10 + NR3
         I20 = I15 + NC1*NVAR
         I30 = I20 + NC1*NVAR
         I60 = I30 + NVAR*NC1
         I70 = I60 + NROW
         I80 = I70 + NVAR*NVAR*NC1
         I90 = I80 + NZMAT
         LAST = I90
         NEED = LAST - I10
         CALL GETFM(NEED)
C
C --- MEMORY ALLOCATION
C
C        I10  FINT (HESSIAN IN INTERNAL COORDS)
C        I15  EM   (TOTAL ENERGY DISTRIBUTION MATRIX)
C        I20  D    (NORMAL MODES IN INTERNAL COORDS)
C        I30  B
C        I60  E1   (VECTOR OF INTRINSIC FREQUENCIES)
C        I70  P    (3-DIM'L MATRIX OF VIBRATIONAL DENSITIES)
C        I80  IZMAT
C ---
C
         IF(EXETYP.NE.CHECK)
     *      CALL DECVIB(VEC,E,X(I10),X(I15),X(I20),X(I30),X(I60),
     *                  X(I70),X(I80),NC1,TFACT,SCLFAC,NFIRST,NLAST)
         CALL RETFM(NEED)
      END IF
C
C     ----- CONVERT FREQUENCIES TO WAVENUMBERS -----
C
      DO 610 I = 1,NC1
         FREQ(I) = SQRT(ABS(TFACT*E(I)))
  610 CONTINUE
C
C     ----- PRINT WARNING FOR SKIPPING FROZEN FREQUENCIES -----
C
      IF(RTYPRT  .AND.  NFRZ.GE.6.  AND.  MASWRK)THEN
            IF(FREQ(NLAST+NFRZ-3).GT.12.AND.FREQ(NLAST+NFRZ-2).LT.30)
     *      WRITE(IW,9065)NLAST+NFRZ-3,NLAST+NFRZ-2
      END IF
C
C     ----- COMPUTE IR INTENSITIES -----
C     PROJECT THE DIPOLE DERIVATIVE TENSOR ONTO EACH NORMAL MODE,
C     AND TAKE THE SQUARE OF THE NORM OF THIS 3 COMPONENT VECTOR
C
      IF (GOTDDM) THEN
         CALL DAREAD (IDAF,IODA,DDM,3*NC1,34,0)
         DO 650 J = 1,NC1
            DDX = DDOT(NC1,VEC(1,J),1,DDM(1,1),3)
            DDY = DDOT(NC1,VEC(1,J),1,DDM(2,1),3)
            DDZ = DDOT(NC1,VEC(1,J),1,DDM(3,1),3)
            SCR(J,1) = DDX*DDX + DDY*DDY + DDZ*DDZ
  650    CONTINUE
      END IF
C
C     ----- COMPUTE RAMAN INTENSITIES -----
C     THE FORMULAE USED BELOW CORRESPOND TO THE EXPERIMENTAL CONDITIONS
C     OF MEASURING THE SCATTERED LIGHT ALONG THE Z AXIS FROM A LINEARLY
C     POLARIZED (ALONG Y AXIS) EXCITATION LASER INCIDENT ON THE SAMPLE
C     FROM THE X DIRECTION.  THE FORMULAE CAN BE FOUND BY G.ZERBI, IN
C        "VIBRATIONAL INTENSITIES IN IR AND RAMAN SPECTROSCOPY",
C        W.B.PERSON AND G.ZERBI, EDS., ELSEVIER NY 1982, PAGE 23.
C     THE SCATTERING INTENSITY WILL BE IN ANGSTROM**4/AMU, AS THE
C     POLARIZABILITY DERIVATIVE IS CONVERTED TO ANGSTROM**2, DOTTED
C     WITH A NORMAL MODE BRINGS UNITS 1/SQRT(AMU), THEN IT IS SQUARED.
C
      TOANGS = 0.52917724924D+00
      UNITS = TOANGS*TOANGS
      IF (GOTADM) THEN
         CALL DAREAD(IDAF,IODA,ADM,6*NC1,249,0)
         CALL DSCAL(6*NC1,UNITS,ADM,1)
         DO 660 J=1,NC1
            AXX = DDOT(NC1,VEC(1,J),1,ADM(1,1),6)
            AXY = DDOT(NC1,VEC(1,J),1,ADM(2,1),6)
            AYY = DDOT(NC1,VEC(1,J),1,ADM(3,1),6)
            AXZ = DDOT(NC1,VEC(1,J),1,ADM(4,1),6)
            AYZ = DDOT(NC1,VEC(1,J),1,ADM(5,1),6)
            AZZ = DDOT(NC1,VEC(1,J),1,ADM(6,1),6)
C  ISOTROPIC AND ANISOTROPIC POLARIZABILITY DERIVATIVE FOR THIS MODE
            APRIM2 = ((AXX+AYY+AZZ)/THREE)**2
            GPRIM2 = ((AXX-AYY)**2 + (AYY-AZZ)**2 + (AZZ-AXX)**2
     *               + SIX*(AXY*AXY + AXZ*AXZ + AYZ*AYZ))/TWO
C  SCATTERING INTENSITY AND DEPOLARIZATION RATIO
            SCR(J,2) = 45.0D+00*APRIM2 + SEVEN*GPRIM2
            DENOM = 45.0D+00*APRIM2 + FOUR*GPRIM2
            SCR(J,3) = ZERO
            IF(DENOM.GT.ZERO) SCR(J,3) = THREE*GPRIM2/DENOM
  660    CONTINUE
      END IF
C
C           PRINT OUT NORMAL COORDINATE INFORMATION
C
      IF(.NOT.MASWRK) GO TO 900
C
      WRITE(IW,9070)
      IF(GOTADM) WRITE(IW,9080)
      IF(.NOT.STATPT) WRITE(IW,9010)
      IF(.NOT.SCFOK)  WRITE(IW,9020)
      MAXCOL = 0
      INCR = 5
      IF(NPRINT.EQ.6) INCR=9
C
      IF(NFRG.GT.0) THEN
         NFGAT=0
         KMASS=0
         DO 705 I=1,NFRG
            DO 700 J=1,NMPTS(I)
               IF(FMASS(J+KMASS).GT.ZERO) NFGAT=NFGAT+1
  700       CONTINUE
            KMASS=KMASS+NMPTS(I)
  705    CONTINUE
         CALL VALFM(LOADFM)
         LFNM = 1    + LOADFM
         LAST = LFNM + 3*NFGAT*INCR
         NEED = LAST - LOADFM -1
         CALL GETFM(NEED)
      END IF
C
C          FIRST, PRINT THE FREQUENCY AND INTENSITY...
C
  710 CONTINUE
      MINCOL = MAXCOL+1
      MAXCOL = MAXCOL+INCR
      IF (MAXCOL .GT. NC1) MAXCOL = NC1
      WRITE (IW,9090)
      WRITE (IW,9100) (J,J = MINCOL,MAXCOL)
      DO 720 J=MINCOL,MAXCOL
         JJ = J + 1 - MINCOL
         LAB(JJ) = LETI
         IF(J.GT.NIMAG) LAB(JJ)=IBLANK
  720 CONTINUE
C
      WRITE(IW,9110) (FREQ(J),LAB(J+1-MINCOL),J = MINCOL,MAXCOL)
      WRITE (IW,9115)    (SCR(J,4),J=MINCOL,MAXCOL)
      IF(GOTDDM) WRITE (IW,9120) (SCR(J,1),J=MINCOL,MAXCOL)
      IF(GOTADM) THEN
          WRITE(IW,9130) (SCR(J,2),J=MINCOL,MAXCOL)
          WRITE(IW,9140) (SCR(J,3),J=MINCOL,MAXCOL)
      END IF
      WRITE(IW,9090)
C
C          THEN THE AB INITIO NORMAL MODE COMPONENTS...
C
      DO 730 IAT = 1,NAT
         I0 = 3*(IAT-1)
         WRITE(IW,9150) IAT,ANAM(IAT),BNAM(IAT),
     *                  CLAB(1),(VEC(I0+1,J),J=MINCOL,MAXCOL)
         WRITE(IW,9160) CLAB(2),(VEC(I0+2,J),J=MINCOL,MAXCOL)
         WRITE(IW,9160) CLAB(3),(VEC(I0+3,J),J=MINCOL,MAXCOL)
  730 CONTINUE
C
C          THEN EFFECTIVE FRAGMENT NORMAL MODE COMPONENTS...
C
      IF(NFRG.GT.0) THEN
         LOC = LFNM
         DO 770 IMODE = MINCOL,MAXCOL
            L=0
            IDEX=3*NAT
            KMASS=0
            DO 760 I=1,NFRG
               DO 740 J=1,6
                  IDEX=IDEX+1
                  STEP(J)=VEC(IDEX,IMODE)
  740          CONTINUE
               DO 750 J=1,NMPTS(I)
                  L=L+1
                  IF (FMASS(J+KMASS).GT.ZERO) THEN
                     CX=EFC(1,L)-EFCENT(1,I)
                     CY=EFC(2,L)-EFCENT(2,I)
                     CZ=EFC(3,L)-EFCENT(3,I)
                     X1=STEP(1)+CZ*STEP(5)-CY*STEP(6)
                     X2=STEP(2)-CZ*STEP(4)+CX*STEP(6)
                     X3=STEP(3)+CY*STEP(4)-CX*STEP(5)
                     X(LOC  ) = X1
                     X(LOC+1) = X2
                     X(LOC+2) = X3
                     LOC = LOC+3
                  END IF
  750          CONTINUE
               KMASS=KMASS+NMPTS(I)
  760       CONTINUE
  770    CONTINUE
         LOC  = LFNM-1
         LTOP = 3*NFGAT*(MAXCOL-MINCOL+1)
         LINC = 3*NFGAT
         IAT  = NAT
         K    = 0
         KMASS= 0
         DO 790 IFRG=1,NFRG
            DO 780 III=1,NMPTS(IFRG)
               K=K+1
               IF(FMASS(III+KMASS).GT.ZERO) THEN
                  IAT=IAT+1
                  WRITE(IW,9170) IAT,IFRG,FRGNME(K),
     *                           CLAB(1),(X(LOC+L),L=1,LTOP,LINC)
                  WRITE(IW,9160) CLAB(2),(X(LOC+L),L=2,LTOP,LINC)
                  WRITE(IW,9160) CLAB(3),(X(LOC+L),L=3,LTOP,LINC)
                  LOC = LOC+3
               END IF
  780       CONTINUE
            KMASS=KMASS+NMPTS(IFRG)
  790    CONTINUE
      END IF
C
C          ...AND FINALLY, THE SAYVETZ CONDITIONS.
C
      WRITE(IW,9090)
      WRITE(IW,9180) CLAB(1),(SVTZT(1,I),I=MINCOL,MAXCOL)
      WRITE(IW,9160) CLAB(2),(SVTZT(2,I),I=MINCOL,MAXCOL)
      WRITE(IW,9160) CLAB(3),(SVTZT(3,I),I=MINCOL,MAXCOL)
      WRITE(IW,9200) (SVTZTT(I),I=MINCOL,MAXCOL)
      WRITE(IW,9090)
      WRITE(IW,9190) CLAB(1),(SVTZR(1,I),I=MINCOL,MAXCOL)
      WRITE(IW,9160) CLAB(2),(SVTZR(2,I),I=MINCOL,MAXCOL)
      WRITE(IW,9160) CLAB(3),(SVTZR(3,I),I=MINCOL,MAXCOL)
      WRITE(IW,9200) (SVTZRT(I),I=MINCOL,MAXCOL)
      IF (MAXCOL.LT.NC1) GO TO 710
C
      IF(NFRG.GT.0) CALL RETFM(NEED)
      WRITE(IW,9210)
      WRITE(IW,9220)
C
C     ----- PUNCH NORMAL MODES -----
C
      WRITE(IP,8000)
      WRITE(IP,8010)
C
C          PUNCH MASSES OF THE ATOMS...
C
      IF (NFRG.NE.0) THEN
         MTOT=0
         DO 810 I=1,NFRG
  810    MTOT=MTOT+NMPTS(I)
         CALL VALFM(LOADFM)
         I10 = 1 + LOADFM
         I20 = I10 +NAT+MTOT
         NEED = I20 -I10
         CALL GETFM(NEED)
C
         IF(NAT.GT.0) CALL DCOPY(NAT,ZMASS,1,X(I10),1)
         IDEX=I10+NAT
         K=0
         KMASS=0
         DO 840 IFRG=1,NFRG
            DO 830 J=1,NMPTS(IFRG)
               K=K+1
               IF (FMASS(J+KMASS).GT.ZERO) THEN
                  X(IDEX)=FMASS(J+KMASS)
                  IDEX=IDEX+1
               END IF
  830       CONTINUE
            KMASS=KMASS+NMPTS(IFRG)
  840    CONTINUE
         WRITE(IP,8020) (X(IAT),IAT=I10,IDEX-1)
         CALL RETFM(NEED)
      ELSE
         WRITE(IP,8020) (ZMASS(IAT),IAT=1,NAT)
      END IF
C
C          PUNCH COMPONENTS OF THE NORMAL MODES...
C
      DO 880 IMODE=1,NC1
         WRITE(IP,8030) IMODE,FREQ(IMODE)
         IF(NAT.GT.0) WRITE(IP,8040) (VEC(I,IMODE),I=1,NAT*3)
         L=0
         IDEX=NAT*3
         KMASS=0
         DO 870 I=1,NFRG
            DO 850 J=1,6
               IDEX=IDEX+1
               STEP(J)=VEC(IDEX,IMODE)
  850       CONTINUE
            DO 860 J=1,NMPTS(I)
               L=L+1
               IF (FMASS(J+KMASS).GT.ZERO) THEN
                  CX=EFC(1,L)-EFCENT(1,I)
                  CY=EFC(2,L)-EFCENT(2,I)
                  CZ=EFC(3,L)-EFCENT(3,I)
                  X1=STEP(1)+CZ*STEP(5)-CY*STEP(6)
                  X2=STEP(2)-CZ*STEP(4)+CX*STEP(6)
                  X3=STEP(3)+CY*STEP(4)-CX*STEP(5)
                  WRITE(IP,8040) X1,X2,X3
               END IF
  860       CONTINUE
            KMASS=KMASS+NMPTS(I)
  870    CONTINUE
  880 CONTINUE
      WRITE(IP,8050)
C
C     ----- THERMOCHEMISTRY ANALYSIS -----
C
  900 CONTINUE
      ANHZPE = ZERO
      DEMP2  = ZERO
      NSWITCH=0
      NM=0
      CALL THERMO(IW,ANHZPE,DEMP2,CC,COM,ZMS,NPART,NSWITCH,NM)
      IF(.NOT.STATPT .AND. .NOT.GOTFRQ .AND. MASWRK) WRITE(IW,9010)
      IF(.NOT.SCFOK  .AND. .NOT.GOTFRQ .AND. MASWRK) WRITE(IW,9020)
C
      IF(RUNTYP.EQ.G3MP2 .AND. NFIRST.GE.3) THEN
         IF(MASWRK) WRITE (IW,9230)
         CALL ABRT
      ENDIF
C
      RETURN
C                            PUNCH FORMATS
 8000 FORMAT('----- START OF NORMAL MODES FOR -MOLPLT- PROGRAM -----')
 8010 FORMAT('ATOMIC MASSES')
 8020 FORMAT(5F12.5)
 8030 FORMAT('MODE',I5,'   FREQUENCY=',F10.5,' (CM**-1)')
 8040 FORMAT(1P,3E17.9)
 8050 FORMAT('----- END OF NORMAL MODES FOR -MOLPLT- PROGRAM -----')
C                            PRINT FORMATS
 9000 FORMAT(/10X,56(1H-)/
     *  10X,'NORMAL COORDINATE ANALYSIS IN THE HARMONIC APPROXIMATION'/
     *        10X,56(1H-))
 9010 FORMAT(/
     *   5X,'*******************************************************'/
     *   5X,'* THIS IS NOT A STATIONARY POINT ON THE MOLECULAR PES *'/
     *   5X,'*     THE VIBRATIONAL ANALYSIS IS NOT VALID !!!       *'/
     *   5X,'*******************************************************')
 9020 FORMAT(/
     *   5X,'*******************************************************'/
     *   5X,'*     THE VIBRATIONAL ANALYSIS IS NOT VALID !!!       *'/
     *   5X,'* THE SCF COMPUTATION DID NOT CONVERGE AT ONE OR MORE *'/
     *   5X,'* OF THE DISPLACED GEOMETRIES.  YOU SHOULD CONVERGE   *'/
     *   5X,'* THESE IN SEPARATE RUNTYP=GRADIENT JOBS, INSERT THE  *'/
     *   5X,'* CORRECT GRADIENTS INTO THE $VIB GROUP FROM YOUR     *'/
     *   5X,'* XXX.IRC FILE, THEN RESUBMIT THIS JOB AS A RESTART.  *'/
     *   5X,'*******************************************************')
 9030 FORMAT(/10X,'ATOMIC WEIGHTS (AMU)'/)
 9040 FORMAT(I5,5X,A8,A2,F15.5)
 9050 FORMAT(/1X,'* * * WARNING, MODE',I2,' HAS BEEN CHOSEN AS A ',
     *          'VIBRATION'/10X,'WHILE MODE',I2,
     *          ' IS ASSUMED TO BE A TRANSLATION/ROTATION.'/
     *       1X,'PLEASE VERIFY THE PROGRAM''S DECISION MANUALLY!'/)
 9060 FORMAT(/1X,'MODES',I2,' TO',I2,' ARE TAKEN AS ROTATIONS',
     *       ' AND TRANSLATIONS.')
 9063 FORMAT(/1X,'MODES',I3,' TO',I3,' ARE INTERNAL ',
     *          'VIBRATIONS ',
     *          'OF FROZEN ATOMS AND FRAGMENTS.',
     *          /1X,'MODES',I3,' TO',I3,' ARE RELATIVE ',
     *          'VIBRATIONS ',
     *          'BETWEEN FROZEN AND UNFROZEN ATOMS.',
     *          /1X,'MODES',I3,' TO',I3,' DO NOT CONTRIBUTE TO ',
     *          'VIBRATIONAL PARTITION FUNCTION',
     *          ' AND ENERGY.',/)
 9064 FORMAT(/1X,'MODES',I3,' TO',I3,' ARE RELATIVE ',
     *          'VIBRATIONS ',
     *          'BETWEEN FROZEN AND UNFROZEN ATOMS.',
     *          /1X,'MODES',I3,' TO',I3,' DO NOT CONTRIBUTE TO ',
     *          'VIBRATIONAL ',
     *          'ENERGIES AND PARTITION FUNCTION.',/)
 9065 FORMAT(//1X,'* * * WARNING ! * * *',//
     *          1X,'MODE',I3,' IS TAKEN AS A RELATIVE VIBRATION ',
     *          'BETWEEN FROZEN AND UNFROZEN ATOMS ',
     *          /1X,'WHILE MODE',I3,' IS TAKEN AS AN INTERNAL ',
     *          'VIBRATION OF UNFROZEN ATOMS.',//
     *          1X,'PLEASE VERIFY THE PROGRAM''S DECISION MANUALLY !'/)
C
 9070 FORMAT(/5X,'FREQUENCIES IN CM**-1, IR INTENSITIES IN DEBYE**2/',
     *           'AMU-ANGSTROM**2,'/5X,'REDUCED MASSES IN AMU.')
 9080 FORMAT(5X,'RAMAN INTENSITIES IN ANGSTROM**4/AMU,',
     *          ' DEPOLARIZATIONS ARE DIMENSIONLESS')
 9090 FORMAT(1X)
 9100 FORMAT(20X,9(4X,I3,5X))
 9110 FORMAT(1X,'      FREQUENCY:',3X,9(F10.2,A2))
 9115 FORMAT(1X,'   REDUCED MASS:',3X,9(F10.5,2X))
 9120 FORMAT(1X,'   IR INTENSITY:',3X,9(F10.5,2X))
 9130 FORMAT(1X,'RAMAN INTENSITY:',3X,9(F10.3,2X))
 9140 FORMAT(1X,' DEPOLARIZATION:',3X,9(F10.3,2X))
 9150 FORMAT(I3,3X,A8,A2,A4,9F12.8)
 9160 FORMAT(16X,A4,9F12.8)
 9170 FORMAT(I3,I3,2X,A8,A4,9F12.8)
 9180 FORMAT(16H TRANS. SAYVETZ ,A4,9F12.8)
 9190 FORMAT(16H   ROT. SAYVETZ ,A4,9F12.8)
 9200 FORMAT(15X,5HTOTAL,9F12.8)
 9210 FORMAT(/1X,'REFERENCE ON SAYVETZ CONDITIONS -',
     *         ' A. SAYVETZ, J.CHEM.PHYS., 7, 383-389(1939).')
 9220 FORMAT(/1X,'NOTE - THE MODES J,K ARE ORTHONORMALIZED',
     *           ' ACCORDING TO'/
     *        1X,'SUM ON I   M(I) * (X(I,J)*X(I,K) + Y(I,J)*',
     *           'Y(I,K) + Z(I,J)*Z(I,K)) = DELTA(J,K)')
 9230 FORMAT(1X,'G3MP2: ONLY FIRST TWO MODE SETS ARE CONSIDERED!')
      END
C*MODULE VIBANL  *DECK CENMAS
      SUBROUTINE CENMAS(NAT,NFRG,C,COM,ZMASST,CMASS,ZMASS)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION C(3,NAT),CMASS(3),COM(3,*),ZMASS(*)
C
      PARAMETER (MXFRG=1050, MXFGPT=12000)
C
      COMMON /FGRAD / DEF(3,MXFGPT),DEFT(3,MXFRG),TORQ(3,MXFRG),
     *                EFCENT(3,MXFRG),FRGMAS(MXFRG),FRGMI(6,MXFRG),
     *                ATORQ(3,MXFRG)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     CALCULATE TOTAL MASS AND CENTER OF MASS
C
      ZMASST=ZERO
      CALL VCLR(CMASS,1,3)
C
      DO 200 I=1,NAT
        AMASS=ZMASS(I)
        ZMASST=ZMASST+AMASS
        DO 100 J=1,3
          CMASS(J)=CMASS(J)+AMASS*C(J,I)
  100   CONTINUE
  200 CONTINUE
C
      DO 400 I=1,NFRG
        AMASS=FRGMAS(I)
        ZMASST=ZMASST+AMASS
        DO 300 J=1,3
          CMASS(J)=CMASS(J)+AMASS*EFCENT(J,I)
  300   CONTINUE
  400 CONTINUE
C
      DO 500 I=1,3
        CMASS(I)=CMASS(I)/ZMASST
  500 CONTINUE
C
C     GET CENTER OF MASS COORDINATES IN MASS-WEIGHTED CARTESIAN
C
      DO 700 I=1,NAT
        IF (ZMASS(I).GT.ZERO) THEN
          AMASS = SQRT(ZMASS(I))
        ELSE
          AMASS = ONE
        END IF
        DO 600 J=1,3
          COM(J,I) = AMASS * (C(J,I)-CMASS(J))
  600   CONTINUE
  700 CONTINUE
C
      II=NAT
      DO 900 I=1,NFRG
        II=II+1
        AMASS = SQRT(FRGMAS(I))
        DO 800 J=1,3
          COM(J,II) = AMASS * (EFCENT(J,I)-CMASS(J))
  800   CONTINUE
  900 CONTINUE
C
      RETURN
      END
C*MODULE VIBANL  *DECK DECVIB
      SUBROUTINE DECVIB(VEC,E,FINT,EM,D,B,E1,P,IZMAT,NC1,TFACT,
     *                  SCLFAC,NFIRST,NLAST)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION VEC(NC1,NC1),E(NC1),FINT(NROW,NROW),EM(NVAR,NC1),
     *          D(NVAR,NC1),B(NROW,*),E1(NROW),
     *          P(NVAR,NVAR,NC1),IZMAT(*)
      DIMENSION LAB(99),ICOOR(7)
C
      LOGICAL LINEAR,DECOMP,GOPARR,MASWRK,DSKWRK
C
      PARAMETER (MXATM=2000)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /VBDCMP/ DECOMP,NROW,NDEG
      COMMON /ZMAT  / NZMAT,NZVAR,NVAR,NSYMC,LINEAR
C
      PARAMETER (ZERO=0.0D+00,ONE=1.0D+00)
      PARAMETER (CONV=1.556914D+01,TM3=1.0D-03)
C
C      ----- "CONV" CONVERTS UNITS OF HARTREE/BOHR**2 TO
C      ----- MILLIDYNES/ANGSTROM.
C      ----- CONV = (4.3598E-08 ERG/HARTREE) *
C      -----       [(1 BOHR/0.529177 ANGSTROM)**2] *
C      -----        (1E+08 ANGSTROM/CM)
C
C
      DATA IBLANK/4H    /,  LETI/4H   I/
      DATA ICOOR/4HSTR.,4HBEND,4HTORS,4HOPBE,4HLIBE,4H5DIH,4HGTRS/
C
C     ----- CALCULATE NORMAL MODES IN INTERNAL COORDINATES (D) -----
C
      CALL DAREAD(IDAF,IODA,B,NROW*NC1,49,0)
      CALL DAREAD(IDAF,IODA,FINT,NROW*NROW,43,0)
      DO 110 I=1,NVAR
         DO 120 J=1,NC1
            SUM = ZERO
            DO 130 K=1,NC1
               SUM = SUM + B(I,K)*VEC(K,J)
  130       CONTINUE
            D(I,J) = SUM
  120    CONTINUE
  110 CONTINUE
C
C     ---- CALCULATE VIBRATIONAL DENSITY MATRICES
C
      CALL VCLR(P,1,NVAR*NVAR*NC1)
      DO 160 I=1,NC1
         IF ((I .GE. NFIRST) .AND. (I .LE. NLAST)) GO TO 160
         DO 150 N=1,NVAR
            DO 140 M=1,NVAR
               P(M,N,I) = D(M,I)*FINT(M,N)*D(N,I)/E(I)
  140       CONTINUE
  150    CONTINUE
  160 CONTINUE
C
C     ---- CHECK TO SEE IF ALL ELEMENTS SUM TO 1.0
C
      CALL VCLR(EM,1,NVAR*NC1)
      IERR = 0
      DO 190 I=1,NC1
         IF ((I .GE. NFIRST) .AND. (I .LE. NLAST)) GO TO 190
         SUM = ZERO
         DO 180 N=1,NVAR
            SUM1 = ZERO
            DO 170 M=1,NVAR
               SUM  = SUM  + P(M,N,I)
               SUM1 = SUM1 + P(M,N,I)
  170       CONTINUE
            EM(N,I) = SUM1
  180    CONTINUE
         IF (ABS(SUM-ONE).GT.TM3 .AND. MASWRK) THEN
            WRITE(IW,9000) I,SUM
            IERR = IERR + 1
         END IF
  190 CONTINUE
C
      IF (IERR.GT.0 .AND. MASWRK) WRITE(IW,9002)
C
C     ---- CALCULATE INTRINSIC VIBRATIONAL FREQUENCIES
C
      DO 258 J=1,NVAR
         SUM2 = ZERO
         DO 256 I=1,NC1
            SUM1 = ZERO
            DO 254 K=1,NVAR
               SUM1 = SUM1 + FINT(J,K)*D(K,I)
  254       CONTINUE
            SUM2 = SUM2 + D(J,I)*SUM1
  256    CONTINUE
         E1(J) = SUM2
  258 CONTINUE
C
C     ---- PRINT INTRINSIC FREQUENCIES
C
      IF (MASWRK) WRITE (IW,9010)  SCLFAC,SCLFAC*SCLFAC
      CALL DAREAD(IDAF,IODA,IZMAT,NZMAT,40,1)
      IF (MASWRK) THEN
      IF (NPRINT .NE. 6) THEN
      WRITE (IW,9014)
      WRITE (IW,9028)
      IADD=1
      NDEGEN=0
      DO 400 I = 1,NVAR
         FREQ=SQRT(ABS(E1(I)*TFACT))
         SFREQ = SCLFAC*FREQ
         IMAG=IBLANK
         IF(E1(I).LT.ZERO) IMAG=LETI
         ITYPE=IZMAT(IADD)
         GO TO (310,320,330,340,350,360,370), ITYPE
  310    CONTINUE
         WRITE(IW,9120) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+3
         GO TO 400
  320    CONTINUE
         WRITE(IW,9130) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+4
         GO TO 400
  330    CONTINUE
         WRITE(IW,9140) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+5
         GO TO 400
  340    CONTINUE
         WRITE(IW,9140) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+5
         GO TO 400
  350    CONTINUE
         WRITE(IW,9130) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),FREQ,IMAG,SFREQ,IMAG
         NDEGEN=NDEGEN+1
         IF(NDEGEN.EQ.2) IADD=IADD+4
         IF(NDEGEN.EQ.2) NDEGEN=0
         GO TO 400
  360    CONTINUE
         WRITE(IW,9150) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),
     *                  IZMAT(IADD+5),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+6
         GO TO 400
  370    CONTINUE
         WRITE(IW,9160) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *                  IZMAT(IADD+6),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+7
  400 CONTINUE
      WRITE (IW,9016)
      WRITE (IW,9028)
      IADD=1
      NDEGEN=0
      DO 450 I = 1,NVAR
         FREQ=SQRT(ABS(E1(I)*TFACT))
         SFREQ = SCLFAC*FREQ
         IMAG=IBLANK
         IF(E1(I).LT.ZERO) IMAG=LETI
         ITYPE=IZMAT(IADD)
         GO TO (432,434,436,438,440,442,444), ITYPE
  432    CONTINUE
C
C      ----- IF INTERNAL COORDINATE IS A STRETCH, COMPUTE
C      ----- INTRINSIC FORCE CONSTANTS.
C
         EM1 = ZMASS(IZMAT(IADD+1))
         EM2 = ZMASS(IZMAT(IADD+2))
         REDMSS = EM1*EM2/(EM1+EM2)
         FCHB = REDMSS*FREQ*FREQ/TFACT
         IF (IMAG .EQ. LETI) FCHB = -FCHB
         FCMA = FCHB*CONV
         SFCHB = SCLFAC*SCLFAC*FCHB
         SFCMA = SCLFAC*SCLFAC*FCMA
         WRITE(IW,9121) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  FCHB,SFCHB,FCMA,SFCMA
         IADD=IADD+3
         GO TO 450
  434    CONTINUE
         IADD=IADD+4
         GO TO 450
  436    CONTINUE
         IADD=IADD+5
         GO TO 450
  438    CONTINUE
         IADD=IADD+5
         GO TO 450
  440    CONTINUE
         NDEGEN=NDEGEN+1
         IF(NDEGEN.EQ.2) IADD=IADD+4
         IF(NDEGEN.EQ.2) NDEGEN=0
         GO TO 450
  442    CONTINUE
         IADD=IADD+6
         GO TO 450
  444    CONTINUE
         IADD=IADD+7
  450 CONTINUE
      ELSE
      WRITE (IW,9012)
      WRITE (IW,9028)
      IADD=1
      NDEGEN=0
      DO 480 I = 1,NVAR
         FREQ=SQRT(ABS(E1(I)*TFACT))
         SFREQ = SCLFAC*FREQ
         IMAG=IBLANK
         IF(E1(I).LT.ZERO) IMAG=LETI
         ITYPE=IZMAT(IADD)
         GO TO (462,464,466,468,470,472,474), ITYPE
  462    CONTINUE
C
C      ----- IF INTERNAL COORDINATE IS A STRETCH, COMPUTE
C      ----- INTRINSIC FORCE CONSTANTS.
C
         EM1 = ZMASS(IZMAT(IADD+1))
         EM2 = ZMASS(IZMAT(IADD+2))
         REDMSS = EM1*EM2/(EM1+EM2)
         FCHB = REDMSS*FREQ*FREQ/TFACT
         IF (IMAG .EQ. LETI) FCHB = -FCHB
         FCMA = FCHB*CONV
         SFCHB = SCLFAC*SCLFAC*FCHB
         SFCMA = SCLFAC*SCLFAC*FCMA
         WRITE(IW,9120) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  FREQ,IMAG,SFREQ,IMAG,FCHB,SFCHB,FCMA,SFCMA
         IADD=IADD+3
         GO TO 480
  464    CONTINUE
         WRITE(IW,9130) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+4
         GO TO 480
  466    CONTINUE
         WRITE(IW,9140) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+5
         GO TO 480
  468    CONTINUE
         WRITE(IW,9140) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+5
         GO TO 480
  470    CONTINUE
         WRITE(IW,9130) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),FREQ,IMAG,SFREQ,IMAG
         NDEGEN=NDEGEN+1
         IF(NDEGEN.EQ.2) IADD=IADD+4
         IF(NDEGEN.EQ.2) NDEGEN=0
         GO TO 480
  472    CONTINUE
         WRITE(IW,9150) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),
     *                  IZMAT(IADD+5),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+6
         GO TO 480
  474    CONTINUE
         WRITE(IW,9160) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *                  IZMAT(IADD+6),FREQ,IMAG,SFREQ,IMAG
         IADD=IADD+7
  480 CONTINUE
      END IF
      END IF
C
C     ----- CONVERT FREQUENCIES TO WAVENUMBERS -----
C
      NIMAG = 0
      DO 490 I = 1,NC1
         IF(E(I).LT.ZERO) NIMAG=NIMAG+1
         E1(I) = SQRT(ABS(TFACT*E(I)))
  490 CONTINUE
C
C     ----- PRINT TED MATRIX (EM)  -----
C
      IF (MASWRK) THEN
      WRITE (IW,9028)
      WRITE (IW,9018)
      WRITE (IW,9028)
      INCR = 4
      IF(NPRINT.EQ.6) INCR=9
      MAX = 0
  500 CONTINUE
      MIN = MAX+1
      MAX = MAX+INCR
      IF (MAX .GT. NC1) MAX = NC1
      WRITE (IW,9028)
      WRITE (IW,9048) (J,J = MIN,MAX)
      WRITE (IW,9068) (E1(I),I=MIN,MAX)
      IF(NIMAG.LT.MIN) GO TO 520
         DO 510 J=MIN,MAX
            LAB(J) = LETI
            IF(J.GT.NIMAG) LAB(J)=IBLANK
  510    CONTINUE
         WRITE(IW,9078) (LAB(J),J=MIN,MAX)
  520 CONTINUE
      WRITE (IW,9028)
      IADD=1
      NDEGEN=0
      DO 700 I = 1,NVAR
         ITYPE=IZMAT(IADD)
         GO TO (610,620,630,640,650,660,670), ITYPE
  610    CONTINUE
         WRITE(IW,9220) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  (EM(I,J),J=MIN,MAX)
         IADD=IADD+3
         GO TO 700
  620    CONTINUE
         WRITE(IW,9230) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),(EM(I,J),J=MIN,MAX)
         IADD=IADD+4
         GO TO 700
  630    CONTINUE
         WRITE(IW,9240) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *        IZMAT(IADD+3),IZMAT(IADD+4),(EM(I,J),J=MIN,MAX)
         IADD=IADD+5
         GO TO 700
  640    CONTINUE
         WRITE(IW,9240) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *        IZMAT(IADD+3),IZMAT(IADD+4),(EM(I,J),J=MIN,MAX)
         IADD=IADD+5
         GO TO 700
  650    CONTINUE
         WRITE(IW,9230) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),(EM(I,J),J=MIN,MAX)
         NDEGEN=NDEGEN+1
         IF(NDEGEN.EQ.2) IADD=IADD+4
         IF(NDEGEN.EQ.2) NDEGEN=0
         GO TO 700
  660    CONTINUE
         WRITE(IW,9250) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),
     *                  IZMAT(IADD+5),(EM(I,J),J=MIN,MAX)
         IADD=IADD+6
         GO TO 700
  670    CONTINUE
         WRITE(IW,9260) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *                  IZMAT(IADD+6),(EM(I,J),J=MIN,MAX)
         IADD=IADD+7
  700 CONTINUE
      IF (MAX .LT. NC1) GO TO 500
C
C     ----- PRINT NORMAL MODES TRANSFORMED TO INTERNAL COORDINATES -----
C
      WRITE (IW,9028)
      WRITE (IW,9098)
      WRITE (IW,9028)
      INCR = 4
      IF (NPRINT.EQ.6) INCR = 9
      MAX = 0
  710 CONTINUE
      MIN = MAX+1
      MAX = MAX+INCR
      IF (MAX .GT. NC1) MAX = NC1
      WRITE (IW,9028)
      WRITE (IW,9048) (J,J = MIN,MAX)
      WRITE (IW,9068) (E1(I),I=MIN,MAX)
      IF(NIMAG.LT.MIN) GO TO 730
         DO 720 J=MIN,MAX
            LAB(J) = LETI
            IF(J.GT.NIMAG) LAB(J)=IBLANK
  720    CONTINUE
         WRITE(IW,9078) (LAB(J),J=MIN,MAX)
  730 CONTINUE
      WRITE (IW,9028)
      IADD=1
      NDEGEN=0
      DO 900 I = 1,NVAR
         ITYPE=IZMAT(IADD)
         GO TO (810,820,830,840,850,860,870), ITYPE
  810    CONTINUE
         WRITE(IW,9220) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  (D(I,J),J=MIN,MAX)
         IADD=IADD+3
         GO TO 900
  820    CONTINUE
         WRITE(IW,9230) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),(D(I,J),J=MIN,MAX)
         IADD=IADD+4
         GO TO 900
  830    CONTINUE
         WRITE(IW,9240) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *        IZMAT(IADD+3),IZMAT(IADD+4),(D(I,J),J=MIN,MAX)
         IADD=IADD+5
         GO TO 900
  840    CONTINUE
         WRITE(IW,9240) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *        IZMAT(IADD+3),IZMAT(IADD+4),(D(I,J),J=MIN,MAX)
         IADD=IADD+5
         GO TO 900
  850    CONTINUE
         WRITE(IW,9230) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),(D(I,J),J=MIN,MAX)
         NDEGEN=NDEGEN+1
         IF(NDEGEN.EQ.2) IADD=IADD+4
         IF(NDEGEN.EQ.2) NDEGEN=0
         GO TO 900
  860    CONTINUE
         WRITE(IW,9250) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),
     *                  IZMAT(IADD+5),(D(I,J),J=MIN,MAX)
         IADD=IADD+6
         GO TO 900
  870    CONTINUE
         WRITE(IW,9260) ICOOR(ITYPE),IZMAT(IADD+1),IZMAT(IADD+2),
     *                  IZMAT(IADD+3),IZMAT(IADD+4),IZMAT(IADD+5),
     *                  IZMAT(IADD+6),(D(I,J),J=MIN,MAX)
         IADD=IADD+7
  900 CONTINUE
      IF (MAX .LT. NC1) GO TO 710
      END IF
      RETURN
C
 9000 FORMAT(1X,'*** WARNING IN DECVIB: SUM OVER VIB. DENSITY MATRIX ',
     *       I3,' NOT EQUAL TO 1.0, BUT IS ',F12.8)
 9002 FORMAT(/1X,'PERHAPS TRY SETTING PURIFY=.T. IN $FORCE....')
 9010 FORMAT(//1X,79(1H-)/9X,'INTRINSIC VIBRATIONAL FREQUENCIES',
     *          ' AND FORCE CONSTANTS'/1X,'FREQUENCIES AND',
     *          ' FORCE CONSTANTS IN PARENTHESES SCALED BY ',
     *          F5.3,' AND ',F5.3,', RESP.'
     *          /1X,79(1H-))
 9012 FORMAT( 1X,'INTERNAL COORDINATE',8X,'INTRINSIC FREQUENCIES',
     *       17X,'INTRINSIC FORCE CONSTANTS',
     *       15X,'INTRINSIC FORCE CONSTANTS'/
     *       33X,'(CM**-1)',28X,'(HARTREES/BOHR**2)',
     *       28X,'(MDYN/ANG)'/131(1H-))
 9014 FORMAT( 1X,'INTERNAL COORDINATE',8X,'INTRINSIC FREQUENCIES'/
     *       33X,'(CM**-1)'/1X,50(1H-))
 9016 FORMAT(/1X,68(1H-)/1X,'INTERNAL COORDINATE',16X,'INTRINSIC FORCE',
     *           ' CONSTANTS'/
     *       28X,'(HARTREES/BOHR**2)',8X,'(MDYN/ANG)'/1X,68(1H-))
 9018 FORMAT(//6X,36(1H-)/6X,'TOTAL ENERGY DISTRIBUTION MATRIX M',
     *         5X,'(VIBR. CONTRIBUTIONS SUM TO 100 %)'/
     *         6X,36(1H-)//' PULAY AND TOROK, ACTA CHIM. ACAD. ',
     *         'SCI. HUNG. 47, 273-297 (1966)')
 9028 FORMAT(1X)
 9048 FORMAT(24X,9(4X,I3,5X))
 9068 FORMAT(25X,9G12.5)
 9078 FORMAT(27X,9(1X,A4,7X))
 9098 FORMAT(//10X,36(1H-)/10X,'NORMAL MODES IN INTERNAL',
     *       ' COORDINATES',/10X,36(1H-))
 9120 FORMAT(1X,A4,2I3,15X,F6.1,A4,3X,'(',F6.1,A4,')',
     *       16X,F8.4,3X,'(',F8.4,')',21X,F7.3,3X,'(',F7.3,')')
 9121 FORMAT(1X,A4,2I3,
     *       16X,F8.4,2X,'(',F8.4,')',4X,F7.3,2X,'(',F7.3,')')
 9130 FORMAT(1X,A4,3I3,12X,F6.1,A4,3X,'(',F6.1,A4,')')
 9140 FORMAT(1X,A4,4I3, 9X,F6.1,A4,3X,'(',F6.1,A4,')')
 9150 FORMAT(1X,A4,5I3, 6X,F6.1,A4,3X,'(',F6.1,A4,')')
 9160 FORMAT(1X,A4,6I3, 3X,F6.1,A4,3X,'(',F6.1,A4,')')
 9220 FORMAT(1X,A4,2I3,12X,9F12.8)
 9230 FORMAT(1X,A4,3I3, 9X,9F12.8)
 9240 FORMAT(1X,A4,4I3, 6X,9F12.8)
 9250 FORMAT(1X,A4,5I3, 3X,9F12.8)
 9260 FORMAT(1X,A4,6I3,    9F12.8)
      END
C*MODULE VIBANL  *DECK EFGRDM
      SUBROUTINE EFGRDM(NCOORD,NCF,NFRG,G,RTRMS,GM,IOPT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION G(*),RTRMS(*),GM(*)
C
      PARAMETER (ZERO=0.0D+00)
C
      IF (NFRG.EQ.0) RETURN
C
C     IOPT=0: TO MASS-WEIGHTED
C     IOPT=1: TO CARTESSIAN
C
      NCF2=NCF/2
      K0=NCOORD
      K1=NCOORD
      DO 900 I=1,NFRG
      DO 500 J=1,3
      K0=K0+1
      K1=K1+1
      IF (IOPT.EQ.0) THEN
      GM(K1)=G(K1)*RTRMS(K0)
      ELSE
      GM(K1)=G(K1)/RTRMS(K0)
      END IF
  500 CONTINUE
      L=K1
      DO 800 J=1,3
      K1=K1+1
      JDX=NCOORD+NCF2+(I-1)*9+IOPT*9*NFRG+(J-1)*3
      DUM=ZERO
      DO 600 M=1,3
      DUM=DUM+RTRMS(JDX+M)*G(L+M)
  600 CONTINUE
      GM(K1)=DUM
  800 CONTINUE
  900 CONTINUE
C
      RETURN
      END
C*MODULE VIBANL  *DECK HESMAS
      SUBROUTINE HESMAS(NCOORD,NCF,NFRG,HESS,RTRMS,A,IOPT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DIMENSION HESS(NCOORD+NCF,NCOORD+NCF),RTRMS(*),A(*)
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
C     IOPT=0: TO MASS-WEIGHTED
C     IOPT=1: TO CARTESION
C     LET A=AB INITIO ATOMS, T=FRAGMENT TRANSLATION, R=FRAGMENT ROTATION
C
      NCF2=NCF/2
      NC2=NCOORD+NCF2
C                          A/A BLOCK
      IJ = 0
      IF (IOPT.EQ.0) THEN
      DO 130 I=1,NCOORD
        RTRMSI = RTRMS(I)
        DO 120 J=1,I
          IJ = IJ + 1
          A(IJ) = RTRMSI * HESS(I,J) * RTRMS(J)
  120   CONTINUE
  130 CONTINUE
      ELSE
      DO 230 I=1,NCOORD
        RTRMSI = ONE/RTRMS(I)
        DO 220 J=1,I
          IJ = IJ + 1
          A(IJ) = RTRMSI * HESS(I,J) / RTRMS(J)
  220   CONTINUE
  230 CONTINUE
      END IF
C
      IF (NFRG.EQ.0) GO TO 1000
C
      DO 500 I=1,NFRG
C                                     A/T BLOCK
        INDXI=NCOORD+(I-1)*6
        INDXJ=NCOORD+(I-1)*3
        DO 370 J=1,3
          IDXH=INDXI+J
          IDXR=INDXJ+J
          DO 310 K=1,NCOORD
            IJ=IJ+1
            IF (IOPT.EQ.0) THEN
              A(IJ) = RTRMS(K) * HESS(K,IDXH) * RTRMS(IDXR)
            ELSE
              A(IJ) = HESS(K,IDXH) / (RTRMS(K) * RTRMS(IDXR))
            END IF
  310     CONTINUE
C                                     T-OTHERS/T BLOCK
          DO 350 K=1,I-1
            INDXK=NCOORD+(K-1)*6
            JNDXK=NCOORD+(K-1)*6+3
            INDXL=NCOORD+(K-1)*3
            DO 320 L=1,3
              IJ=IJ+1
              IF (IOPT.EQ.0) THEN
                A(IJ) = RTRMS(INDXL+L) * HESS(INDXK+L,IDXH) *
     *                  RTRMS(IDXR)
              ELSE
                A(IJ) = HESS(INDXK+L,IDXH) / (RTRMS(INDXL+L) *
     *                  RTRMS(IDXR))
              END IF
  320       CONTINUE
C
C                                     R-OTHERS/T BLOCK
            DO 340 L=1,3
              JNDXL=NC2+(K-1)*9+IOPT*9*NFRG+(L-1)*3
              DUM=ZERO
              IJ=IJ+1
              DO 330 M=1,3
                DUM=DUM+RTRMS(JNDXL+M)*HESS(JNDXK+M,IDXH)
  330         CONTINUE
              IF (IOPT.EQ.0) THEN
                A(IJ) = DUM * RTRMS(IDXR)
              ELSE
                A(IJ) = DUM / RTRMS(IDXR)
              END IF
  340       CONTINUE
  350     CONTINUE
C
C                                     T-SELF/T BLOCK
          DO 360 K=1,J
            IJ=IJ+1
            IF (IOPT.EQ.0) THEN
              A(IJ) = RTRMS(INDXJ+K) * HESS(INDXI+K,IDXH) *
     *                RTRMS(IDXR)
            ELSE
              A(IJ) = HESS(INDXI+K,IDXH) / (RTRMS(INDXJ+K) *
     *                RTRMS(IDXR))
            END IF
  360     CONTINUE
  370   CONTINUE
C
        JNDXI=NCOORD+(I-1)*6+3
        DO 490 J=1,3
          JNDXJ=NC2+(I-1)*9+IOPT*9*NFRG+(J-1)*3
C                                     A/R BLOCK
          DO 420 K=1,NCOORD
            IJ=IJ+1
            DUM=ZERO
            DO 410 L=1,3
              DUM=DUM+RTRMS(JNDXJ+L)*HESS(JNDXI+L,K)
  410       CONTINUE
            IF (IOPT.EQ.0) THEN
              A(IJ) = RTRMS(K) * DUM
            ELSE
              A(IJ) = DUM / RTRMS(K)
            END IF
  420     CONTINUE
C
C                                     T/R BLOCK
          DO 480 K=1,I
            INDXK=NCOORD+(K-1)*6
            INDXL=NCOORD+(K-1)*3
C
            DO 440 L=1,3
              KDXH=INDXK+L
              IJ=IJ+1
              DUM=ZERO
              DO 430 M=1,3
                DUM=DUM+RTRMS(JNDXJ+M)*HESS(JNDXI+M,KDXH)
  430         CONTINUE
              IF (IOPT.EQ.0) THEN
                A(IJ) = RTRMS(INDXL+L) * DUM
              ELSE
                A(IJ) = DUM / RTRMS(INDXL+L)
              END IF
  440       CONTINUE
C
C                                     R/R BLOCK
            LEND=3
            IF (K.EQ.I) LEND=J
            KNDXI=NCOORD+(K-1)*6+3
            DO 470 L=1,LEND
              KNDXJ=NC2+(K-1)*9+IOPT*9*NFRG+(L-1)*3
              IJ=IJ+1
              DUM=ZERO
              DO 460 M=1,3
                DUM1=ZERO
                DO 450 N=1,3
                  DUM1=DUM1+RTRMS(KNDXJ+N)*HESS(KNDXI+N,JNDXI+M)
  450           CONTINUE
                DUM=DUM+DUM1*RTRMS(JNDXJ+M)
  460         CONTINUE
              A(IJ)=DUM
  470       CONTINUE
  480     CONTINUE
  490   CONTINUE
  500 CONTINUE
C
 1000 IJ=0
      DO 1001 I=1,NCOORD+NCF
      DO 1001 J=1,I
      IJ=IJ+1
      HESS(I,J)=A(IJ)
 1001 HESS(J,I)=A(IJ)
      RETURN
      END
C*MODULE VIBANL  *DECK INRTIA
      SUBROUTINE INRTIA(C,COM,ZMASS,VMOI,NPART)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION C(3,NPART),COM(3,NPART),ZMASS(NPART),VMOI(3)
      DIMENSION CMASS(3),TMOI(6),WRK(3,8),IWRK(3),TROT(3,3)
C
      PARAMETER (ZERO=0.0D+00)
C
C     COMPUTE PRODUCT OF INERTIA TENSOR, AND RETURN ITS EIGENVALUES
C     IN VMOI.  NOTE THAT WHAT IS RETURNED ARE THE PRODUCTS OF
C     INERTIA, "MOMENT OF INERTIA" MEANS IXX=YY+ZZ (AND SO ON),
C     WHERE THE PRODUCTS OF INERTIA XX,YY,ZZ ARE THE 3 EIGENVALUES.
C     ON ENTRY, C=NUCLEAR COORDS, ZMASS=NUCLEAR MASSES
C     ON EXIT, COM=CENTER OF MASS COORDS (NOT ROTATED TO PRINCIPAL
C     AXES, WHICH WOULD BE EASY), VMOI=PRODUCTS OF INERTIA.
C
C     ----- COMPUTE THE CENTER OF MASS -----
C
      TOTWT=ZERO
      CMASS(1)=ZERO
      CMASS(2)=ZERO
      CMASS(3)=ZERO
C
      DO 120 I=1,NPART
         WEIGHT=ZMASS(I)
         TOTWT=TOTWT+WEIGHT
         DO 110 J=1,3
            CMASS(J)=CMASS(J)+WEIGHT*C(J,I)
  110    CONTINUE
  120 CONTINUE
C
      DO 140 I=1,3
         CMASS(I)=CMASS(I)/TOTWT
  140 CONTINUE
C
C     ----- TRANSLATE TO CENTER OF MASS -----
C
      DO 180 I=1,NPART
         DO 170 J=1,3
            COM(J,I) = C(J,I)-CMASS(J)
  170    CONTINUE
  180 CONTINUE
C
C     ----- COMPUTE THE MOMENT OF INERTIA TENSOR -----
C
      XX=ZERO
      YY=ZERO
      ZZ=ZERO
      XY=ZERO
      XZ=ZERO
      YZ=ZERO
      DO 200 I=1,NPART
         WEIGHT=ZMASS(I)
         XC=COM(1,I)
         YC=COM(2,I)
         ZC=COM(3,I)
         XX=XX + WEIGHT*(YC*YC+ZC*ZC)
         YY=YY + WEIGHT*(XC*XC+ZC*ZC)
         ZZ=ZZ + WEIGHT*(XC*XC+YC*YC)
         XY=XY - WEIGHT*XC*YC
         XZ=XZ - WEIGHT*XC*ZC
         YZ=YZ - WEIGHT*YC*ZC
  200 CONTINUE
C
      TMOI(1) = XX
      TMOI(2) = XY
      TMOI(3) = YY
      TMOI(4) = XZ
      TMOI(5) = YZ
      TMOI(6) = ZZ
      IERR = 0
      CALL GLDIAG(3,3,3,TMOI,WRK,VMOI,TROT,IERR,IWRK)
      IF(IERR.NE.0) CALL ABRT
C
      RETURN
      END
C*MODULE VIBANL  *DECK PRJFC
      SUBROUTINE PRJFC(PRJGRD,PRJROT,TOTM,F,DX,COORD,RM,FRGMI,
     *                 BUF1,BUF2,NATM,NCC,NCCF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL PRJGRD,PRJROT
C
      DIMENSION COORD(3,*),RM(*),FRGMI(*),DX(3,*),
     *          F(NCCF,NCCF),BUF1(NCCF,NCCF),BUF2(NCCF,NCCF)
C
      PARAMETER (HALF=0.5D+00, ONE=1.0D+00)
C
C     ----- PREPARE FOR PROJECTION OF FORCE CONSTANT MATRIX -----
C
C     PROJECTION OF THE FORCE CONSTANT MATRIX AT AN EQUILIBRIUM
C     GEOMETRY IS DONE WITH AN IDENTICALLY ZERO GRADIENT VECTOR
C     SO THAT THE TRANSLATIONAL AND VIBRATIONAL CONTAMINANTS
C     ONLY ARE ELIMINATED.
C
C     PROJECTION OF THE F.C.M. AT A POINT WITH A NON-ZERO GRADIENT
C     ALSO ELIMINATES THESE CONTAMINANTS, BUT ALSO PROJECTS THE
C     MATRIX SO THAT ONE OF ITS NORMAL MODES LIES PARALLEL TO
C     THE MASS-WEIGHTED GRADIENT.  THIS MODE WILL HAVE A ZERO
C     FREQUENCY, AND THE OTHER 3N-7 MODES WILL BE ORTHOGONAL TO
C     THIS MODE.  THIS TYPE OF PROJECTION IS USEFUL IN DYNAMICS.
C
C     FOR A DESCRIPTION OF THE METHOD, SEE W.H.MILLER, N.C.HANDY,
C     J.E.ADAMS, IN J.CHEM.PHYS. 72, 99-112(1980).
C
C     PRJGRD    = IF .TRUE., THE PROJECTION WILL ZERO THE NORMAL
C                 MODE PARALLEL TO THE GRADIENT VECTOR.  USEFUL
C                 IN TRANSITION STATE THEORY.
C     PRJROT    = IF .TRUE., THE PROJECTION WILL ZERO THE NORMAL
C                 MODES CORRESPONDING TO ROTATIONS
C     TOTM      = TOTAL MASS OF SYSTEM (INCLUDING FRAGMENT)
C     F         = ON ENTRY, MASS-WEIGHTED FORCE CONSTANT MATRIX.
C                 ON EXIT, PROJECTED TO REMOVE T,R AND POSSIBLY
C                 THE GRADIENT'S DEGREE OF FREEDOM
C     DX        = ON ENTRY, MASS-WEIGHTED GRADIENT VECTOR.
C                 ON EXIT, THIS WILL HAVE BEEN DESTROYED
C                 THIS IS NOT USED IF PRJGRD IS FALSE.
C     COORD     = CARTESIAN COORDINATE SHIFTED TO CENTER OF MASS.
C     RM        = INVERSE OF THE SQUARE ROOT MASS.
C     FRGMI     = MOMENTS OF INITIA OF FRAGMENTS. (LOW TRIANGLES)
C     BUF1,BUF2 = WORK BUFFERS (NCCF*NCCF)
C     NATM      = NUMBER OF ATOMS.
C     NCC       = NO OF CARTESIAN COORDINATES (3*NATM)
C     NCCF      = NO OF COORDINATES PLUS 6 PER EFFECTIVE FRAGMENT
C
C     ----- NORMALIZE THE GRADIENT -----
C
      IF(PRJGRD) THEN
         GNORM = DDOT(NCCF,DX,1,DX,1)
         GNORM = ONE/SQRT(GNORM)
         CALL DSCAL(NCCF,GNORM,DX,1)
      END IF
C
C     ----- CARRY OUT THE PROJECTION -----
C
      CALL PRJFCM(PRJGRD,PRJROT,TOTM,F,COORD,DX,RM,FRGMI,
     *            BUF1,BUF2,NATM,NCC,NCCF)
C
C     ----- ENFORCE SYMMETRY UPON F.C.M. -----
C
      DO 320 I=1,NCCF
         DO 310 J=1,I
            AVRG   = HALF*(F(I,J)+F(J,I))
            F(J,I) = AVRG
            F(I,J) = AVRG
  310    CONTINUE
  320 CONTINUE
      RETURN
      END
C*MODULE VIBANL  *DECK PRJFCM
      SUBROUTINE PRJFCM(PRJGRD,PRJROT,TOTM,F,X,DX,RM,FRGMI,
     *                  P,BUF2,NATM,NCC,NCCF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,MASWRK,DSKWRK,PRJGRD,PRJROT
C
      DIMENSION X(*),DX(*),RM(*),FRGMI(6,*),F(NCCF,NCCF),
     *          P(NCCF,NCCF),BUF2(NCCF,NCCF)
      DIMENSION TENS(3,3,3),ROT(3,3),SCR(3,3),ISCR(6),DETERM(2)
C
      COMMON /IOFILE/ IR,IW,IPP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, CUT8=1.0D-08)
C
C         -TENS- IS "THE USUAL TOTALLY ASYMMETRIC CARTESIAN TENSOR"
C
      DATA TENS/ 0.0D+00,  0.0D+00,  0.0D+00,
     X           0.0D+00,  0.0D+00, -1.0D+00,
     X           0.0D+00,  1.0D+00,  0.0D+00,
     Y           0.0D+00,  0.0D+00,  1.0D+00,
     Y           0.0D+00,  0.0D+00,  0.0D+00,
     Y          -1.0D+00,  0.0D+00,  0.0D+00,
     Z           0.0D+00, -1.0D+00,  0.0D+00,
     Z           1.0D+00,  0.0D+00,  0.0D+00,
     Z           0.0D+00,  0.0D+00,  0.0D+00  /
C
C  ----- CALCULATE PROJECTED FORCE CONSTANT MATRIX -----
C  COPIED FROM POLYRATE BY S. KOSEKI, DEC. 2, 1985.
C
C   X : MASS-WEIGHTED COORDINATE
C   DX: NORMALIZED GRADIENT VECTOR
C   F : MASS-WEIGHTED FORCE CONSTANT MATRIX
C   RM: INVERSE OF SQUARE ROOT OF MASS
C   P, BUF2: WORK ARRAYS
C
      NFRG=(NCCF-NCC)/6
C
C
C       --- COMPUTE INERTIA TENSOR ---
C
      CALL VCLR(ROT,1,9)
C
      IF(.NOT.PRJROT) GO TO 200
C
      DO 100 I=1,NATM+NFRG
         L=3*(I-1)+1
         ROT(1,1)=ROT(1,1)+X(L+1)**2+X(L+2)**2
         ROT(1,2)=ROT(1,2)-X(L)*X(L+1)
         ROT(1,3)=ROT(1,3)-X(L)*X(L+2)
         ROT(2,2)=ROT(2,2)+X(L)**2+X(L+2)**2
         ROT(2,3)=ROT(2,3)-X(L+1)*X(L+2)
         ROT(3,3)=ROT(3,3)+X(L)**2+X(L+1)**2
 100  CONTINUE
      DO 101 I=1,NFRG
         ROT(1,1)=ROT(1,1)+FRGMI(1,I)
         ROT(1,2)=ROT(1,2)+FRGMI(2,I)
         ROT(2,2)=ROT(2,2)+FRGMI(3,I)
         ROT(1,3)=ROT(1,3)+FRGMI(4,I)
         ROT(2,3)=ROT(2,3)+FRGMI(5,I)
         ROT(3,3)=ROT(3,3)+FRGMI(6,I)
 101  CONTINUE
      ROT(2,1)=ROT(1,2)
      ROT(3,1)=ROT(1,3)
      ROT(3,2)=ROT(2,3)
C
      CHK=ROT(1,1)*ROT(2,2)*ROT(3,3)
      IF(ABS(CHK).GT.CUT8) GO TO 180
C
C       --- COMPUTE INVERSE IN CASE OF SOME ZERO MOMENTS ---
C
      IF(ABS(ROT(1,1)).GT.CUT8) GO TO 130
      IF(ABS(ROT(2,2)).GT.CUT8) GO TO 120
      IF(ABS(ROT(3,3)).GT.CUT8) GO TO 110
      IF(MASWRK) WRITE(IW,7010) ROT(1,1),ROT(2,2),ROT(3,3)
      RETURN
C
C             X,Y=0 BUT Z.NE.0
  110 ROT(3,3)=ONE/ROT(3,3)
      GO TO 170
C
C             Y.NE.0
  120 IF(ABS(ROT(3,3)).GT.CUT8) GO TO 160
C             X,Z=0 BUT Y.NE.0
      ROT(2,2)=ONE/ROT(2,2)
      GO TO 170
C
C             X.NE.0
  130 IF(ABS(ROT(2,2)).GT.CUT8) GO TO 140
      IF(ABS(ROT(3,3)).GT.CUT8) GO TO 150
C
C             Y,Z=0 BUT X.NE.0
      ROT(1,1)=ONE/ROT(1,1)
      GO TO 170
C
C             X,Y.NE.0 BUT Z=0
  140 DET=ROT(1,1)*ROT(2,2)-ROT(1,2)*ROT(2,1)
      TRP=ROT(1,1)
      ROT(1,1)=ROT(2,2)/DET
      ROT(2,2)=TRP/DET
      ROT(1,2)=-ROT(1,2)/DET
      ROT(2,1)=-ROT(2,1)/DET
      GO TO 170
C
C             X,Z.NE.0 BUT Y=0
  150 DET=ROT(1,1)*ROT(3,3)-ROT(1,3)*ROT(3,1)
      TRP=ROT(1,1)
      ROT(1,1)=ROT(3,3)/DET
      ROT(3,3)=TRP/DET
      ROT(1,3)=-ROT(1,3)/DET
      ROT(3,1)=-ROT(3,1)/DET
      GO TO 170
C
C             Y,Z.NE.0 BUT X=0
  160 DET=ROT(3,3)*ROT(2,2)-ROT(3,2)*ROT(2,3)
      TRP=ROT(3,3)
      ROT(3,3)=ROT(2,2)/DET
      ROT(2,2)=TRP/DET
      ROT(3,2)=-ROT(3,2)/DET
      ROT(2,3)=-ROT(2,3)/DET
C
  170 CONTINUE
      GO TO 200
C
C       --- COMPUTE INVERSE FOR ALL MOMENTS NONZERO CASE ---
C
  180 CONTINUE
      INFO=0
      CALL DGEFA(ROT,3,3,ISCR,INFO)
C
      IF(INFO.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'INERTIA MATRIX IS SINGULAR'
         IF(MASWRK) WRITE(IW,*) 'HESSIAN WILL NOT BE PROJECTED!'
         RETURN
      END IF
C
      CALL DGEDI(ROT,3,3,ISCR,DETERM,SCR,1)
C
C        --- COMPUTE P MATRIX ---
C        THE FORMULA USED IS EQN 4.11 OF MILLER, HANDY, ADAMS.
C
  200 CONTINUE
C
C        ATOM AND FRAG CENTER
C
      DO 280 IP=1,NATM+NFRG
         INDX=3*(IP-1)
         KNDX=MAX(3*(IP-1),6*(IP-1)-3*NATM)
         DO 270 JP=1,IP
            JNDX=3*(JP-1)
            LNDX=MAX(3*(JP-1),6*(JP-1)-3*NATM)
            DO 260 IC=1,3
               II=INDX+IC
               KK=KNDX+IC
               JEND=3
               IF(JP.EQ.IP) JEND=IC
               DO 250 JC=1,JEND
                  JJ=JNDX+JC
                  LL=LNDX+JC
                  SUM=ZERO
                  IF(PRJROT) THEN
                     DO 240 IA=1,3
                        DO 230 IB=1,3
                           IF(TENS(IA,IB,IC).EQ.ZERO) GO TO 230
                           DO 220 JA=1,3
                              DO 210 JB=1,3
                                 IF(TENS(JA,JB,JC).EQ.ZERO) GO TO 210
                                 SUM=SUM+TENS(IA,IB,IC)*TENS(JA,JB,JC)
     *                               *ROT(IA,JA)*X(INDX+IB)*X(JNDX+JB)
  210                         CONTINUE
  220                      CONTINUE
  230                   CONTINUE
  240                CONTINUE
                  END IF
                  P(KK,LL)=SUM
                  IF(PRJGRD) P(KK,LL)=P(KK,LL) + DX(KK)*DX(LL)
                  IF(IC.EQ.JC) P(KK,LL)=P(KK,LL)
     *                         + ONE/(RM(II)*RM(JJ)*TOTM)
  250          CONTINUE
  260       CONTINUE
  270    CONTINUE
  280 CONTINUE
C
C        ATOM AND FRAG CENTER WITH FRAG ROT
C
      INDX=3*NATM+12*NFRG
      KNDX=3*NATM+3
      DO 380 IP=1,NFRG
         DO 370 IC=1,3
            KK=KNDX+IC
            DO 360 JP=1,NATM+NFRG
               JNDX=3*(JP-1)
               LNDX=MAX(3*(JP-1),6*(JP-1)-3*NATM)
               DO 350 JC=1,3
                  LL=LNDX+JC
                  SUM=ZERO
                  IF(PRJROT) THEN
                     DO 340 JA=1,3
                        DO 330 JB=1,3
                           IF(TENS(JA,JB,JC).EQ.ZERO) GO TO 330
                           DO 320 IA=1,3
                              SUM=SUM+TENS(JA,JB,JC)
     *                            *ROT(IA,JA)*X(JNDX+JB)*RM(INDX+IA)
  320                      CONTINUE
  330                   CONTINUE
  340                CONTINUE
                  END IF
                  IF (KK.GE.LL) THEN
                     P(KK,LL)=SUM
                     IF(PRJGRD) P(KK,LL)=P(KK,LL)+DX(KK)*DX(LL)
                  ELSE
                     P(LL,KK)=SUM
                     IF(PRJGRD) P(LL,KK)=P(LL,KK)+DX(KK)*DX(LL)
                  END IF
  350          CONTINUE
  360       CONTINUE
            INDX=INDX+3
  370    CONTINUE
         KNDX=KNDX+6
  380 CONTINUE
C
C        FRAG ROT
C
      INDX0=3*NATM+12*NFRG
      KNDX0=3*NATM+3
      INDX=INDX0
      KNDX=KNDX0
      DO 480 IP=1,NFRG
         DO 470 IC=1,3
            KK=KNDX+IC
            JNDX=INDX0
            LNDX=KNDX0
            DO 460 JP=1,IP
               JCEND=3
               IF (JP.EQ.IP) JCEND=IC
               DO 450 JC=1,JCEND
                  LL=LNDX+JC
                  SUM=ZERO
                  IF(PRJROT) THEN
                     DO 440 IA=1,3
                        DO 430 JA=1,3
                           SUM=SUM+ROT(IA,JA)*RM(INDX+IA)*RM(JNDX+JA)
  430                   CONTINUE
  440                CONTINUE
                  END IF
                  P(KK,LL)=SUM
                  IF(PRJGRD) P(KK,LL)=P(KK,LL)+DX(KK)*DX(LL)
                  JNDX=JNDX+3
  450          CONTINUE
               LNDX=LNDX+6
  460       CONTINUE
            INDX=INDX+3
  470    CONTINUE
         KNDX=KNDX+6
  480 CONTINUE
C
C        --- COMPUTE (I-P) ---
C
      DO 520 I=1,NCCF
         DO 510 J=1,I
            P(I,J)=-P(I,J)
            IF(I.EQ.J) P(I,J) = ONE + P(I,J)
            IF(ABS(P(I,J)).LT.CUT8) P(I,J)=ZERO
            P(J,I)=P(I,J)
  510    CONTINUE
  520 CONTINUE
C
C        --- PROJECT BY COMPUTING F = (I-P)*(F*(I-P)) ---
C
      CALL TFSQU(BUF2,F,P,DX,NCCF,NCCF)
      CALL DCOPY(NCCF*NCCF,BUF2,1,F,1)
C
      RETURN
C
 7010 FORMAT(/1X,'PRJFCM: EVERY DIAGONAL ELEMENT IS ZERO?',3F20.10)
      END
C*MODULE VIBANL  *DECK PRJGRD
      SUBROUTINE PRJGRD(GRAD,C,X,GTOT,RM,P,ZMASS,NATM,NATM3,NCCF)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXFRG=1050, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      LOGICAL GOPARR,MASWRK,DSKWRK
C
      DIMENSION GRAD(NATM3),X(NCCF),GTOT(NCCF),RM(*),P(NCCF,NCCF),
     *          C(3,*),ZMASS(*)
      DIMENSION TENS(3,3,3),ROT(3,3),SCR(3,3),ISCR(6),DETERM(2),TMP(3)
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 /IOFILE/ IR,IW,IPP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, CUT8=1.0D-08)
C
C         -TENS- IS "THE USUAL TOTALLY ASYMMETRIC CARTESIAN TENSOR"
C
      DATA TENS/ 0.0D+00,  0.0D+00,  0.0D+00,
     X           0.0D+00,  0.0D+00, -1.0D+00,
     X           0.0D+00,  1.0D+00,  0.0D+00,
     Y           0.0D+00,  0.0D+00,  1.0D+00,
     Y           0.0D+00,  0.0D+00,  0.0D+00,
     Y          -1.0D+00,  0.0D+00,  0.0D+00,
     Z           0.0D+00, -1.0D+00,  0.0D+00,
     Z           1.0D+00,  0.0D+00,  0.0D+00,
     Z           0.0D+00,  0.0D+00,  0.0D+00  /
C
C    ----- PROJECT ROTATION/TRANSLATION FROM A GRADIENT VECTOR -----
C            COPIED FROM POLYRATE BY S. KOSEKI, DEC. 2, 1985.
C      MODIFIED TO ACT UPON THE GRADIENT VECTOR BY MWS, MARCH 2002
C
C   THE ONLY NECESSARY INITIALIZED INPUTS ARE -GRAD- AND -NATM3- AND
C   -NCCF-, AND IN COMMON -DEFT- AND -TORQ-.  ON EXIT, THE AB INITIO
C   PART OF THE GRADIENT -GRAD- AND THE FRAGMENT COMPONENTS -DEFT-
C   AND -TORQ- WILL HAVE BEEN PROJECTED.
C   C   : CARTESIAN COORDINATES OF ALL ATOMS
C   THE FOLLOWING ARE SCRATCH STORAGE:
C   X   : PROVIDE STORAGE FOR THE MASS-WEIGHTED COORDINATES, AND
C   GTOT: PROVIDE STORAGE FOR COMBINED A.I. AND EFP GRADIENT VECTOR
C   RM  : PROVIDE STORAGE FOR THE INVERSE SQUARE ROOTS OF MASSES
C   P   : PROVIDE STORAGE FOR THE PROJECTOR MATRIX, INCLUDING
C         EFP COMPONENTS, NAMELY NCCF= 3*NAT + 6*NFRG.
C
C        THE STATUS OF THE EFP CODE, AS OF JULY 2, 2002 IS
C          A) SEE NOTES IN MIKE'S DFT FOLDER ABOUT THE NUMERICS
C          B) COMPLETE CONFIDENCE THAT $DATA PART IS CORRECT, BUT
C             ADDING EFPS IS IFFY.  THE CODE OF APRIL OMITTED THE
C             MASS WEIGHTING OF TORQUES, AND COULD OPTIMIZE METHANOL
C             PLUS ONE EFP AND WATER PLUS ONE EFP FAIRLY WELL.
C          C) MASS WEIGHTING OF THE TORQUE WAS ADDED IN JULY AFTER
C             STUDY OF THE -HESMAS- ROUTINE THAT MASS WEIGHTS THE
C             HESSIAN.  THE TEST INPUTS OF WATER PLUS ONE OR PLUS
C             TWO EFPS WERE UNSATISFACTORY  IN GEOMETRY OPTIMIZATION.
C             (APRIL'S INPUT HAD GOTTEN LOST IN THE MEANTIME, SO ITS
C             A DIFFERENT WATER=EFP).  THE BEST RUNS WERE WITH NO
C             PROJECTION AT ALL, SO THE NEXT IF STATEMENT JUST QUITS
C             IF ANY EFP IS PRESENT.  THE CASE OF TWO EFPS IS MUCH
C             HARDER THAN ONE EFP, BY THE WAY.
C          D) PRINTING THE MATRIX I-P BELOW GIVES THE SAME NUMBERS
C             AS IT DOES IF PRINTING IN -PRJFCM- FOR ONE OR TWO EFPS,
C             SO CALLING ARGUMENTS TO THIS ROUTINE ARE CORRECT.
C          E) THE FORMULA FOR P IS UNKNOWN, EXCEPT MAYBE TO PAUL, SO
C             ONE MIGHT WONDER IF THE TWO ROUTINES MAKE THE CORRECT P.
C          F) THE MASS-WEIGHTING IS PROBABLY CORRECT, AS IT SEEMS
C             TO USE M**(-1/2) AND M**(1/2) TO DO THIS (THE COMMENTS
C             IN -EFRMSS- ARE CLEARLY WRONG, AS THE CODE FILLS RM
C             WITH SETS OF THE POWERS OF THESE MOMENTS).  THE CONTENTS
C             OF THE -FRGMI- ARRAYS WERE NOT TESTED FOR CORRECTNESS.
C
      IF(NFRG.GT.0) RETURN
C
C       --- SET UP RECIPROCAL MASSES ---
C
      IF(NFG.EQ.0) CALL RAMS(RM,1)
C     RAMS SHOULD ONLY BE CALLED FOR EFP?
      IF(NFRG.GT.0) CALL EFRMSS(NFRG,RM(NATM3+1))
C
C       GET TOTAL MASS -TOTM- AND MASS-WEIGHTED COORDINATES -X-
C
      CALL CENMAS(NATM,NFRG,C,X,TOTM,SCR,ZMASS)
C
C       --- COMPUTE INERTIA TENSOR ---
C
      CALL VCLR(ROT,1,9)
      DO 100 I=1,NATM+NFRG
         L=3*(I-1)+1
         ROT(1,1)=ROT(1,1)+X(L+1)**2+X(L+2)**2
         ROT(1,2)=ROT(1,2)-X(L)*X(L+1)
         ROT(1,3)=ROT(1,3)-X(L)*X(L+2)
         ROT(2,2)=ROT(2,2)+X(L)**2+X(L+2)**2
         ROT(2,3)=ROT(2,3)-X(L+1)*X(L+2)
         ROT(3,3)=ROT(3,3)+X(L)**2+X(L+1)**2
 100  CONTINUE
      DO 101 I=1,NFRG
         ROT(1,1)=ROT(1,1)+FRGMI(1,I)
         ROT(1,2)=ROT(1,2)+FRGMI(2,I)
         ROT(2,2)=ROT(2,2)+FRGMI(3,I)
         ROT(1,3)=ROT(1,3)+FRGMI(4,I)
         ROT(2,3)=ROT(2,3)+FRGMI(5,I)
         ROT(3,3)=ROT(3,3)+FRGMI(6,I)
 101  CONTINUE
      ROT(2,1)=ROT(1,2)
      ROT(3,1)=ROT(1,3)
      ROT(3,2)=ROT(2,3)
C
      CHK=ROT(1,1)*ROT(2,2)*ROT(3,3)
      IF(ABS(CHK).GT.CUT8) GO TO 180
C
C       --- COMPUTE INVERSE IN CASE OF SOME ZERO MOMENTS ---
C
      IF(ABS(ROT(1,1)).GT.CUT8) GO TO 130
      IF(ABS(ROT(2,2)).GT.CUT8) GO TO 120
      IF(ABS(ROT(3,3)).GT.CUT8) GO TO 110
      IF(MASWRK) WRITE(IW,7010) ROT(1,1),ROT(2,2),ROT(3,3)
      RETURN
C
C             X,Y=0 BUT Z.NE.0
  110 ROT(3,3)=ONE/ROT(3,3)
      GO TO 170
C
C             Y.NE.0
  120 IF(ABS(ROT(3,3)).GT.CUT8) GO TO 160
C             X,Z=0 BUT Y.NE.0
      ROT(2,2)=ONE/ROT(2,2)
      GO TO 170
C
C             X.NE.0
  130 IF(ABS(ROT(2,2)).GT.CUT8) GO TO 140
      IF(ABS(ROT(3,3)).GT.CUT8) GO TO 150
C
C             Y,Z=0 BUT X.NE.0
      ROT(1,1)=ONE/ROT(1,1)
      GO TO 170
C
C             X,Y.NE.0 BUT Z=0
  140 DET=ROT(1,1)*ROT(2,2)-ROT(1,2)*ROT(2,1)
      TRP=ROT(1,1)
      ROT(1,1)=ROT(2,2)/DET
      ROT(2,2)=TRP/DET
      ROT(1,2)=-ROT(1,2)/DET
      ROT(2,1)=-ROT(2,1)/DET
      GO TO 170
C
C             X,Z.NE.0 BUT Y=0
  150 DET=ROT(1,1)*ROT(3,3)-ROT(1,3)*ROT(3,1)
      TRP=ROT(1,1)
      ROT(1,1)=ROT(3,3)/DET
      ROT(3,3)=TRP/DET
      ROT(1,3)=-ROT(1,3)/DET
      ROT(3,1)=-ROT(3,1)/DET
      GO TO 170
C
C             Y,Z.NE.0 BUT X=0
  160 DET=ROT(3,3)*ROT(2,2)-ROT(3,2)*ROT(2,3)
      TRP=ROT(3,3)
      ROT(3,3)=ROT(2,2)/DET
      ROT(2,2)=TRP/DET
      ROT(3,2)=-ROT(3,2)/DET
      ROT(2,3)=-ROT(2,3)/DET
C
  170 CONTINUE
      GO TO 200
C
C       --- COMPUTE INVERSE FOR THE GENERAL CASE ALL MOMENTS NONZERO ---
C
  180 CONTINUE
      INFO=0
      CALL DGEFA(ROT,3,3,ISCR,INFO)
C
      IF(INFO.NE.0) THEN
         IF(MASWRK) WRITE(IW,*) 'INERTIA MATRIX IS SINGULAR'
         IF(MASWRK) WRITE(IW,*) 'GRADIENT WILL NOT BE PROJECTED!'
         RETURN
      END IF
C
      CALL DGEDI(ROT,3,3,ISCR,DETERM,SCR,1)
C
C        --- COMPUTE P MATRIX ---
C        THE FORMULA USED IS EQN 4.11 OF MILLER, HANDY, ADAMS.
C
  200 CONTINUE
C
C        ATOM AND FRAG CENTER
C
      DO 280 IP=1,NATM+NFRG
         INDX=3*(IP-1)
         KNDX=MAX(3*(IP-1),6*(IP-1)-3*NATM)
         DO 270 JP=1,IP
            JNDX=3*(JP-1)
            LNDX=MAX(3*(JP-1),6*(JP-1)-3*NATM)
            DO 260 IC=1,3
               II=INDX+IC
               KK=KNDX+IC
               JEND=3
               IF(JP.EQ.IP) JEND=IC
               DO 250 JC=1,JEND
                  JJ=JNDX+JC
                  LL=LNDX+JC
                  SUM=ZERO
                  DO 240 IA=1,3
                     DO 230 IB=1,3
                        IF(TENS(IA,IB,IC).EQ.ZERO) GO TO 230
                        DO 220 JA=1,3
                           DO 210 JB=1,3
                              IF(TENS(JA,JB,JC).EQ.ZERO) GO TO 210
                              SUM=SUM+TENS(IA,IB,IC)*TENS(JA,JB,JC)
     *                               *ROT(IA,JA)*X(INDX+IB)*X(JNDX+JB)
  210                      CONTINUE
  220                   CONTINUE
  230                CONTINUE
  240             CONTINUE
                  P(KK,LL)=SUM
                  IF(IC.EQ.JC) P(KK,LL)=P(KK,LL)
     *                         + ONE/(RM(II)*RM(JJ)*TOTM)
  250          CONTINUE
  260       CONTINUE
  270    CONTINUE
  280 CONTINUE
C
C        ATOM AND FRAG CENTER WITH FRAG ROT
C
      INDX=3*NATM+12*NFRG
      KNDX=3*NATM+3
      DO 380 IP=1,NFRG
         DO 370 IC=1,3
            KK=KNDX+IC
            DO 360 JP=1,NATM+NFRG
               JNDX=3*(JP-1)
               LNDX=MAX(3*(JP-1),6*(JP-1)-3*NATM)
               DO 350 JC=1,3
                  LL=LNDX+JC
                  SUM=ZERO
                  DO 340 JA=1,3
                     DO 330 JB=1,3
                        IF(TENS(JA,JB,JC).EQ.ZERO) GO TO 330
                        DO 320 IA=1,3
                           SUM=SUM+TENS(JA,JB,JC)
     *                            *ROT(IA,JA)*X(JNDX+JB)*RM(INDX+IA)
  320                   CONTINUE
  330                CONTINUE
  340             CONTINUE
                  IF (KK.GE.LL) THEN
                     P(KK,LL)=SUM
                  ELSE
                     P(LL,KK)=SUM
                  END IF
  350          CONTINUE
  360       CONTINUE
            INDX=INDX+3
  370    CONTINUE
         KNDX=KNDX+6
  380 CONTINUE
C
C        FRAG ROT
C
      INDX0=3*NATM+12*NFRG
      KNDX0=3*NATM+3
      INDX=INDX0
      KNDX=KNDX0
      DO 480 IP=1,NFRG
         DO 470 IC=1,3
            KK=KNDX+IC
            JNDX=INDX0
            LNDX=KNDX0
            DO 460 JP=1,IP
               JCEND=3
               IF (JP.EQ.IP) JCEND=IC
               DO 450 JC=1,JCEND
                  LL=LNDX+JC
                  SUM=ZERO
                  DO 440 IA=1,3
                     DO 430 JA=1,3
                        SUM=SUM+ROT(IA,JA)*RM(INDX+IA)*RM(JNDX+JA)
  430                CONTINUE
  440             CONTINUE
                  P(KK,LL)=SUM
                  JNDX=JNDX+3
  450          CONTINUE
               LNDX=LNDX+6
  460       CONTINUE
            INDX=INDX+3
  470    CONTINUE
         KNDX=KNDX+6
  480 CONTINUE
C
C        --- COMPUTE (I-P) ---
C
      DO 520 I=1,NCCF
         DO 510 J=1,I
            P(I,J)=-P(I,J)
            IF(I.EQ.J) P(I,J) = ONE + P(I,J)
            IF(ABS(P(I,J)).LT.CUT8) P(I,J)=ZERO
            P(J,I)=P(I,J)
  510    CONTINUE
  520 CONTINUE
C
C        --- PROJECT BY COMPUTING G = G*(I-P) ---
C     THE PROJECTION IS MADE UPON THE MASS-WEIGHTED GRADIENT,
C     WHICH IS CONVERTED BACK AFTER THE PROJECTION IS FINISHED.
C
      DO I=1,NATM3
         GTOT(I) = RM(I) * GRAD(I)
      ENDDO
      I=NATM3
      DO IFRG=1,NFRG
         GTOT(I+1)= RM(I+1) * DEFT(1,IFRG)
         GTOT(I+2)= RM(I+2) * DEFT(2,IFRG)
         GTOT(I+3)= RM(I+3) * DEFT(3,IFRG)
         GTOT(I+4)=           TORQ(1,IFRG)
         GTOT(I+5)=           TORQ(2,IFRG)
         GTOT(I+6)=           TORQ(3,IFRG)
         I=I+6
      ENDDO
C                  MASS WEIGHTING FOR THE TORQUE PART
      IBASE=3*NATM
      JBASE=3*NATM+3*NFRG
      DO 610 I=1,NFRG
         JNDXI=IBASE+(I-1)*6+3
         DO J=1,3
            JNDXJ=JBASE+(I-1)*9+(J-1)*3
            DUM=ZERO
            DO L=1,3
               DUM=DUM+RM(JNDXJ+L)*GTOT(JNDXI+L)
            ENDDO
            TMP(J) = DUM
         ENDDO
         DO J=1,3
            GTOT(JNDXI+J) = TMP(J)
         ENDDO
  610 CONTINUE
C
C       LO, THE PROJECTION IS SO SIMPLE, COMPARED TO EVERYTHING ELSE.
C
      DO J=1,NCCF
         SUM=ZERO
         DO K=1,NCCF
            SUM=SUM+GTOT(K)*P(K,J)
         ENDDO
         X(J)=SUM
      ENDDO
C
C        COPY THE GRADIENT BACK, UN-MASS-WEIGHTING AS WE COPY
C
      DO I=1,NATM3
         GRAD(I) = X(I)/RM(I)
C        WRITE(6,*) 'UNMASSING',GRAD(I),X(I),RM(I)
      ENDDO
C
      IBASE=3*NATM
      JBASE=3*NATM+12*NFRG
      DO 670 I=1,NFRG
         JNDXI=IBASE+(I-1)*6+3
         DO J=1,3
            JNDXJ=JBASE+(I-1)*9+(J-1)*3
            DUM=ZERO
            DO L=1,3
               DUM=DUM+RM(JNDXJ+L)*X(JNDXI+L)
            ENDDO
            TMP(J) = DUM
         ENDDO
         DO J=1,3
            X(JNDXI+J) = TMP(J)
         ENDDO
  670 CONTINUE
C
      I=NATM3
      DO IFRG=1,NFRG
         DEFT(1,IFRG)=X(I+1)/RM(I+1)
         DEFT(2,IFRG)=X(I+2)/RM(I+2)
         DEFT(3,IFRG)=X(I+3)/RM(I+3)
         TORQ(1,IFRG)=X(I+4)
         TORQ(2,IFRG)=X(I+5)
         TORQ(3,IFRG)=X(I+6)
         I=I+6
      ENDDO
      RETURN
C
 7010 FORMAT(/1X,'PRJGRD: EVERY DIAGONAL ELEMENT IS ZERO?',3F20.10)
      END
C*MODULE VIBANL  *DECK RAMS
      SUBROUTINE RAMS(RMASS,IOPT)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      DOUBLE PRECISION MASS
C
      LOGICAL FIRST,GOPARR,MASWRK,DSKWRK
C
      PARAMETER (MXATM=2000, NNAM=1)
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      DIMENSION RMASS(*)
      DIMENSION AMS(106),AMASS(MXATM),QNAM(NNAM),KQNAM(NNAM)
C
      SAVE AMASS,FIRST
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 /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      DATA FIRST/.TRUE./
      DATA MASS/8HMASS    /, QNAM/8HAMASS   /
C
C              MASS OF MOST ABUNDANT ISOTOPE, UP TO LR
C              PADDED WITH TWO ZEROS FOR DUMMY ATOMS.
C
      DATA (AMS(I),I=1,54)  /
     *   1.007825D+00,4.0026D+00,7.01600D+00,9.01218D+00,11.00931D+00,
     *   12.0D+00,14.00307D+00,15.99491D+00,18.99840D+00,19.99244D+00,
     *   22.9898D+00,23.98504D+00,26.98153D+00,27.97693D+00,
     *   30.97376D+00,31.97207D+00,34.96885D+00,39.948D+00,
     *   38.96371D+00,39.96259D+00,44.95592D+00,47.90D+00,50.9440D+00,
     *   51.9405D+00,54.9381D+00,55.9349D+00,58.9332D+00,57.9353D+00,
     *   62.9298D+00,63.9291D+00,68.9257D+00,73.9219D+00,74.9216D+00,
     *   79.9165D+00,78.9183D+00,83.9115D+00,
     *   84.9117D+00,87.9056D+00,89.9054D+00,89.9043D+00,92.9060D+00,
     *   97.9055D+00,97.0D+00,101.9037D+00,102.9048D+00,105.9032D+00,
     *   106.9041D+00,113.9036D+00,114.9041D+00,119.9022D+00,
     *   120.9038D+00,129.9067D+00,126.9044D+00,131.9042D+00/
      DATA (AMS(I),I=55,106)  /
     *   132.9054D+00,137.9052D+00,138.9063D+00,139.9054D+00,
     *   140.9076D+00,141.9077D+00,144.9127D+00,151.9197D+00,
     *   152.9212D+00,157.9241D+00,158.9253D+00,163.9292D+00,
     *   164.9303D+00,165.9303D+00,168.9342D+00,173.9389D+00,
     *   174.9408D+00,179.9465D+00,180.9480D+00,183.9509D+00,
     *   186.9557D+00,191.9615D+00,192.9629D+00,194.9648D+00,
     *   196.9665D+00,201.9706D+00,
     *   204.9744D+00,207.9766D+00,208.9804D+00,208.9824D+00,
     *   209.9871D+00,222.0176D+00,
     *   223.0197D+00,226.0254D+00,
     *   227.0278D+00,232.0381D+00,231.0359D+00,238.0508D+00,
     *   237.0482D+00,244.0642D+00,243.0614D+00,247.0703D+00,
     *   247.0703D+00,251.0796D+00,252.0829D+00,257.0751D+00,
     *   258.0986D+00,259.1009D+00,260.1053D+00,261.1087D+00,
     *   2*0.0D+00/
C
C     ----- SET UP ATOMIC MASSES -----
C     IOPT.EQ.0 RETURNS MASS OF EACH ATOM
C     IOPT.EQ.1 RETURNS TRIPLETS OF INVERSE SQUARE ROOT MASSES
C     IOPT.EQ.2 RETURNS PERIODIC TABLE OF THE MASSES
C     THE USER INPUT CAN OVERRIDE THE BUILT IN TABLE,
C     WITH THE DEFAULT BEING THE MOST ABUNDANT ISOTOPE.
C
      IF(IOPT.EQ.0) GO TO 100
      IF(IOPT.EQ.1) GO TO 100
      IF(IOPT.EQ.2) GO TO 500
      IF (MASWRK) WRITE(IW,*) 'ILLEGAL USE OF -RAMS- ROUTINE'
      CALL ABRT
      STOP
C
  100 CONTINUE
      IF(.NOT.FIRST) GO TO 300
C
C        --- READ $MASS GROUP, ONCE ONLY IN THE RUN ---
C        CAUTION - WE ARE READING THE INPUT BEFORE ANY ECP
C        ELECTRONS ARE REMOVED.  ZAN IS THE TRUE CHARGE!
C
      DO 210 IAT = 1,NAT
         ZNUC = ZAN(IAT)
         NUCZ = INT(ZNUC)
         IF(NUCZ.GE.1 .AND. NUCZ.LE.106) THEN
            AMASS(IAT) = AMS(NUCZ)
         ELSE
            AMASS(IAT) = ZERO
         END IF
  210 CONTINUE
C
C       ----- READ NAMELIST $MASS -----
C
      JRET=0
      KQNAM(1) = 10*MXATM + 3
      CALL NAMEIO(IR,JRET,MASS,NNAM,QNAM,KQNAM,
     *            AMASS,
     *            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,0,0,  0,0,0,0,0,  0,0,0,0,0,  0,0,0,0,0)
      IF (JRET.EQ.2) CALL ABRT
      FIRST=.FALSE.
C
  300 CONTINUE
      IF(IOPT.EQ.1) GO TO 400
C
C         RETURN THE NORMAL MASSES
C
      DO 310 I=1,NAT
         RMASS(I) = AMASS(I)
  310 CONTINUE
      RETURN
C
C         TRIPLETS OF THE INVERSE SQUARE ROOTS OF THE
C         MASSES ARE USEFUL FOR MASS WEIGHTED CARTESIANS.
C
  400 CONTINUE
      K = 0
      DO 480 IAT = 1,NAT
         IF(AMASS(IAT).GT.ZERO) THEN
            DUM = ONE/SQRT(AMASS(IAT))
         ELSE
            DUM = ONE
         END IF
         DO 470 J = 1,3
            K = K+1
            RMASS(K) = DUM
  470    CONTINUE
  480 CONTINUE
      RETURN
C
C        RETURN THE BUILT IN TABLE ONLY
C
  500 CONTINUE
      DO 510 I = 1,106
         RMASS(I)=AMS(I)
  510 CONTINUE
      RETURN
      END
C*MODULE VIBANL  *DECK THERMO
      SUBROUTINE THERMO(IW,ANHZPE,DEMP2,CC,COM,ZMS,NPART,NSWITCH,NM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION KT,JPCAL
C
      LOGICAL LINEAR,PRTSCN,GOTFRQ,GOPARR,MASWRK,DSKWRK,ABEL,RTYPRT
C
      CHARACTER*8 FRGNME,CXBASIS
C
      PARAMETER (MXATM=2000, MXFRG=1050, MXPT=2000, MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      DIMENSION CC(3,NPART),COM(3,NPART),ZMS(NPART),DEMP2(NM)
C         BOTH ASVIB AND TFREQ SHOULD HAVE BEEN DYNAMIC STORAGE
      DIMENSION VMOI(3),ASVIB(3*MXATM),TFREQ(3*MXATM)
C
      COMMON /ANHVIB/ ZMASST,NIMAG,NLAST,NSKIP
      COMMON /CXTHRM/ CXTHERM(11),CXZPE,METHCX,ICXBAS,CXBASIS
      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 /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 /FROZNC/ IFREEZ(3*MXATM),NFRZ
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /THERMD/ FREQ(3*MXATM),TEMP(10),SCLFAC,NTEMP,PRTSCN,GOTFRQ
C
      PARAMETER (PLANCK=6.626176D-34, BOLTZ=1.380662D-23,
     *           PASCAL=1.01325D+05,
     *           AVOGAD=6.022045D+23, TOANGS=0.52917724924D+00,
     *           JPCAL=4.184D+00, CLIGHT=2.99792458D+08,
     *           TOCM=2.194746D+05, TOKCAL=627.51D+00,
     *           TWO=2.0D+00, ONEPT5=1.5D+00, TWOPT5=2.5D+00,
     *           ONE=1.0D+00, ZERO=0.0D+00,   TM3=1.0D-03,
     *           EIGHT=8.0D+00)
C
      DATA HSSIAN/8HHESSIAN /, VSCF/8HVSCF    /
C
C     ----- COMPUTE THERMODYNAMIC PROPERTIES -----
C     WRITTEN BY MIKE SCHMIDT, JULY 1988.
C     ON ENTRY, ZMASS CONTAINS ATOMIC MASSES IN AMU (TOTAL=ZMASST),
C               FREQ CONTAINS VIBRATIONAL FREQENCIES IN CM**-1,
C               NIMAG IS NUMBER OF NEGATIVE FORCE CONSTANTS,
C               ZMASST IS TOTAL MASS OF SYSTEM.
C     MOST INTERMEDIATE CALCULATIONS ARE DONE IN MKS UNITS.
C
C     THE FORMULAE USED HERE ARE TAKEN FROM PAGES 115 AND 123-125
C     OF 'STATISTICAL MECHANICS' BY NORMAN DAVIDSON, MCGRAW-HILL,
C     NEW YORK, 1962.
C
      PI = ACOS(-ONE)
      P = PASCAL
      R = BOLTZ * AVOGAD
      BOHR = 5.29177249D-11
      CATOM = 12.011D+00
C
C     SWITCH FROM HARMONIC TO ANHARMONIC APPROXIMATION THERMOCHEMISTRY
C
      IF(NSWITCH.EQ.0) THEN
        IFMIN = NLAST + 1
        IFMAX = 3*NAT+6*NFRG
        IF(GOTFRQ) IFMAX=IFMAX-NSKIP
      ELSE
        IFMIN=1
        IFMAX=NM
      END IF
      DO I=IFMIN,IFMAX
       IF(NSWITCH.EQ.0) THEN
        TFREQ(I)=FREQ(I)
       ELSE
        TFREQ(I)=DEMP2(I)
       END IF
      END DO
                        NPASS=1
      IF(SCLFAC.NE.ONE) NPASS=2
      SCLSAV = SCLFAC
C
      DO 800 ITEMP=1,NTEMP
      T = TEMP(ITEMP)
      RT = R * T
      KT = BOLTZ * T
C
      IPASS=0
  100 CONTINUE
      IPASS=IPASS+1
      IF(IPASS.LT.NPASS) THEN
         SCLFAC = ONE
      ELSE
         SCLFAC = SCLSAV
      END IF
C
      IF (MASWRK) THEN
       IF(NSWITCH.EQ.0) THEN
         WRITE(IW,9000) T,P
       ELSE
         WRITE(IW,9001) T,P
       END IF
         WRITE(IW,9005) SCLFAC
      END IF
C
C     ----- PRINT MESSAGE FOR PARTIAL HESSIAN ANALYSIS -----
C
      RTYPRT = RUNTYP.EQ.HSSIAN .OR. RUNTYP.EQ.VSCF
      IF(RTYPRT  .AND.  IFREEZ(1).NE.0  .AND.  MASWRK) WRITE(IW,9002)
C
C     ----- ELECTRONIC -----
C     SINCE THE EXCITATION TO EXCITED STATES IS UNKNOWN, AND WE ALSO
C     KNOW NOTHING ABOUT SPIN-ORBIT SPLITTINGS, WE ASSUME THAT ONLY
C     ONE ELECTRONIC STATE CONTRIBUTES TO THE PARTITION FUNCTION.
C     IN ADDITION, SINCE WE DON'T KNOW WHAT THE SPATIAL DEGENERACY OF
C     THIS STATE IS, WE TAKE THE ELECTRONIC DEGENERATCY TO BE JUST
C     THE SPIN MULTIPLICITY.
C
      QELEC = MUL
      QELOG = LOG(QELEC)
      EELEC = ZERO
      HELEC = ZERO
      CVELEC = ZERO
      CPELEC = ZERO
      SELEC = R*QELOG
      GELEC = HELEC - T * SELEC
C
C     ----- TRANSLATION -----
C
      WEIGHT = TM3 * ZMASST / AVOGAD
      TTEMP = TWO*PI*(WEIGHT/PLANCK)*(KT/PLANCK)
      QTRAN = TTEMP**ONEPT5 * KT/P
      QTLOG = LOG(QTRAN)
      ETRAN = ONEPT5 * RT
      HTRAN = TWOPT5 * RT
      CVTRAN = ONEPT5 * R
      CPTRAN = TWOPT5 * R
      STRAN = R * (QTLOG + TWOPT5)
      GTRAN = HTRAN - T * STRAN
C
C     ----- SET UP FOR POSSIBLE MONATOMIC -----
C
      IF(NAT.EQ.1  .AND.  NFRG.EQ.0) THEN
         QROT = ONE
         QRLOG = ZERO
         EROT = ZERO
         HROT = ZERO
         CVROT = ZERO
         CPROT = ZERO
         SROT = ZERO
         GROT = ZERO
         QVIB = ONE
         QVLOG = ZERO
         EVIB = ZERO
         HVIB = ZERO
         CVVIB = ZERO
         CPVIB = ZERO
         SVIB = ZERO
         GVIB = ZERO
         GO TO 700
      END IF
C
C     ----- ROTATION -----
C
C     THE ROTATIONAL SYMMETRY NUMBER SHOULD BE SET CORRECTLY
C     BY THE FOLLOWING CODE.
C     NOTE THAT GROUPS 'CINFV', 'DINFH', 'I', AND 'IH' NEVER GET HERE.
C     -IGROUP- IS A POINTER INTO THE FOLLOWING TABLE,
C     DATA GRP /C1   ,CS   ,CI   ,CN   ,S2N  ,CNH  ,
C    *          CNV  ,DN   ,DNH  ,DND  ,CINFV,DINFH,T
C    *          TH   ,TD   ,O    ,OH   ,I    ,IH   /
C
      SIGMA=1.0D+00
      IF(IGROUP.GE.4) SIGMA=NAXIS
      IF(IGROUP.GE.8) SIGMA=NAXIS+NAXIS
      IF(IGROUP.GE.13) SIGMA=12.0D+00
      IF(IGROUP.GE.16) SIGMA=24.0D+00
C
C        --- LOAD -CC- AND -ZMS- WITH COORDINATES AND MASSES ---
C
      DO 110 I=1,NAT
         CC(1,I) = C(1,I)
         CC(2,I) = C(2,I)
         CC(3,I) = C(3,I)
         ZMS(I) = ZMASS(I)
  110 CONTINUE
C
      K = NAT
      IAT = 0
      KMASS=0
      DO 130 IFRG=1,NFRG
         DO 120 III=1,NMPTS(IFRG)
            IAT = IAT+1
            IF(FMASS(III+KMASS).GT.ZERO) THEN
               K = K + 1
               CC(1,K) = EFC(1,IAT)
               CC(2,K) = EFC(2,IAT)
               CC(3,K) = EFC(3,IAT)
               ZMS(K) = FMASS(III+KMASS)
            END IF
  120    CONTINUE
         KMASS=KMASS+NMPTS(IFRG)
  130 CONTINUE
C
      CALL INRTIA(CC,COM,ZMS,VMOI,NPART)
      LINEAR = VMOI(1).LT.TM3
C
      IF(LINEAR) THEN
         SIGMA=ONE
         DO 210 K=2,NAT
            IF(ABS(COM(1,1)+COM(1,K)).LT.TM3  .AND.
     *         ABS(COM(2,1)+COM(2,K)).LT.TM3  .AND.
     *         ABS(COM(3,1)+COM(3,K)).LT.TM3  .AND.
     *         ZMS(1).EQ.ZMS(K)) SIGMA=TWO
  210    CONTINUE
      END IF
C
      IF (MASWRK) WRITE(IW,9010) (VMOI(I),I=1,3),SIGMA
C
      FACT1 = (CATOM*BOHR*BOHR)/(12.0D+03*AVOGAD)
      FACT2 = PLANCK/(8.0D+09*PI*PI)
      ACONST = ZERO
      BCONST = ZERO
      CCONST = ZERO
      IF(VMOI(1).GT.0.001D+00) ACONST = FACT2/(FACT1*VMOI(1))
      IF(VMOI(2).GT.0.001D+00) BCONST = FACT2/(FACT1*VMOI(2))
      IF(VMOI(3).GT.0.001D+00) CCONST = FACT2/(FACT1*VMOI(3))
      IF (MASWRK) WRITE(IW,9011) ACONST, BCONST, CCONST
C
      FACT = TOANGS*TOANGS/AVOGAD
      VMOI(1) = FACT * VMOI(1)
      VMOI(2) = FACT * VMOI(2)
      VMOI(3) = FACT * VMOI(3)
C
      FACT = EIGHT*PI*PI*(KT/PLANCK)*(1.0D-23/PLANCK)
      IF(LINEAR) THEN
         QROT = FACT*VMOI(3)/SIGMA
         QRLOG = LOG(QROT)
         EROT = RT
         CVROT = R
         CPROT = R
         SROT = R * (QRLOG + ONE)
      ELSE
         QROT = SQRT(PI*(FACT*VMOI(1))
     *                 *(FACT*VMOI(2))
     *                 *(FACT*VMOI(3)))/SIGMA
         QRLOG = LOG(QROT)
         EROT  = ONEPT5 * RT
         CVROT = ONEPT5 * R
         CPROT = ONEPT5 * R
         SROT = R * (QRLOG + ONEPT5)
      END IF
      HROT = EROT
      GROT = HROT - T * SROT
C
C     ----- VIBRATION -----
C
      IF((NIMAG.GT.0) .AND. MASWRK) WRITE(IW,9040) NIMAG
C
C     ----- SKIP FREQUECIES FOR PARTIAL HESSIAN ANALYSIS -----
C
      NFRZ = 0
      NCOORD=3*NAT
      DO 323 I=1, NCOORD
        IF(IFREEZ(I) .NE. 0) NFRZ = NFRZ + 1
  323 CONTINUE
      NFZERO = NFRZ - 3
      IF(NFRZ .EQ. 0) NFZERO = 0
      IFMIN=IFMIN+NFZERO
C
      IF(IFREEZ(1).NE.0  .AND.  RTYPRT) IFMIN=IFMIN+6*NFRG
C
      IF(NSWITCH.EQ.0) THEN
         ZPECM = ZERO
         DO I=IFMIN,IFMAX
            ZPECM = ZPECM + TFREQ(I)
         END DO
         ZPECM  = SCLFAC*ZPECM/TWO
      ELSE
         ZPECM = ANHZPE
      END IF
      ZPE    = ZPECM/TOCM
      ZPEKCL = ZPE*TOKCAL
      ZPEKJ  = ZPEKCL*JPCAL
      IF(NSWITCH.EQ.0) THEN
       IF (MASWRK) WRITE(IW,9050) SCLFAC,ZPE,ZPECM,ZPEKCL,ZPEKJ
      ELSE
       IF (MASWRK) WRITE(IW,9051) SCLFAC,ZPE,ZPECM,ZPEKCL,ZPEKJ
      END IF
C
      QVIB = ONE
      EVIB = ZERO
      CVVIB = ZERO
      SVIB = ZERO
      DO 330 I=IFMIN,IFMAX
         IF(TFREQ(I).EQ.ZERO) THEN
            IF (MASWRK) WRITE(IW,9200) I
            CALL ABRT
         END IF
         W = TFREQ(I)*SCLFAC
         U = (PLANCK*W*CLIGHT*1.0D+02)/KT
         EXPU=1.0D+35
         IF(U.LT.80.0D+00) EXPU = EXP(U)
         EXPMU = EXP(-U)
         QVIB = QVIB/(ONE-EXPMU)
         EVIB = EVIB + RT*U/(EXPU-ONE)
         CVVIB = CVVIB + R*U*U*(EXPU/(EXPU-ONE))/(EXPU-ONE)
         ASVIB(I)= R*(U/(EXPU-ONE) - LOG(ONE-EXPMU))
         SVIB = SVIB + ASVIB(I)
  330 CONTINUE
      QVLOG = LOG(QVIB)
      CPVIB = CVVIB
      ETHERM = EVIB
      EVIB = EVIB + ZPEKJ*1.0D+03
      HVIB = EVIB
      GVIB = HVIB - T * SVIB
C
C     ----- TOTALS -----
C
  700 CONTINUE
      QTOT = QELEC * QTRAN * QROT  * QVIB
      QLOG = QELOG + QTLOG + QRLOG + QVLOG
      IF (MASWRK) THEN
         IF(NFRG.GT.0) WRITE(IW,9065)
         WRITE(IW,9070)
         WRITE(IW,9080) 'ELEC. ',QELEC,QELOG
         WRITE(IW,9080) 'TRANS.',QTRAN,QTLOG
         WRITE(IW,9080) 'ROT.  ',QROT,QRLOG
         WRITE(IW,9080) 'VIB.  ',QVIB,QVLOG
         WRITE(IW,9080) 'TOT.  ',QTOT,QLOG
      END IF
      EELEC = TM3 * EELEC
      HELEC = TM3 * HELEC
      GELEC = TM3 * GELEC
      ETRAN = TM3 * ETRAN
      HTRAN = TM3 * HTRAN
      GTRAN = TM3 * GTRAN
      EROT = TM3 * EROT
      HROT = TM3 * HROT
      GROT = TM3 * GROT
      EVIB = TM3 * EVIB
      HVIB = TM3 * HVIB
      GVIB = TM3 * GVIB
      ETOT = EELEC + ETRAN + EROT + EVIB
      HTOT = HELEC + HTRAN + HROT + HVIB
      GTOT = GELEC + GTRAN + GROT + GVIB
      STOT = SELEC + STRAN + SROT + SVIB
      CVTOT = CVELEC + CVTRAN + CVROT + CVVIB
      CPTOT = CPELEC + CPTRAN + CPROT + CPVIB
      IF (MASWRK) THEN
      WRITE(IW,9090)
      WRITE(IW,9100) '    KJ/MOL','    KJ/MOL','    KJ/MOL',
     *               '   J/MOL-K','   J/MOL-K','   J/MOL-K'
      WRITE(IW,9110) 'ELEC. ',EELEC,HELEC,GELEC,CVELEC,CPELEC,SELEC
      WRITE(IW,9110) 'TRANS.',ETRAN,HTRAN,GTRAN,CVTRAN,CPTRAN,STRAN
      WRITE(IW,9110) 'ROT.  ',EROT,HROT,GROT,CVROT,CPROT,SROT
      WRITE(IW,9110) 'VIB.  ',EVIB,HVIB,GVIB,CVVIB,CPVIB,SVIB
      WRITE(IW,9110) 'TOTAL ',ETOT,HTOT,GTOT,CVTOT,CPTOT,STOT
      WRITE(IW,9115) ETHERM
      END IF
      FACT = ONE/JPCAL
      EELEC = FACT*EELEC
      HELEC = FACT*HELEC
      GELEC = FACT*GELEC
      CVELEC = FACT*CVELEC
      CPELEC = FACT*CPELEC
      SELEC = FACT*SELEC
      ETRAN = FACT*ETRAN
      HTRAN = FACT*HTRAN
      GTRAN = FACT*GTRAN
      CVTRAN = FACT*CVTRAN
      CPTRAN = FACT*CPTRAN
      STRAN = FACT*STRAN
      EROT = FACT*EROT
      HROT = FACT*HROT
      GROT = FACT*GROT
      CVROT = FACT*CVROT
      CPROT = FACT*CPROT
      SROT = FACT*SROT
      ETHERM = FACT*ETHERM
      EVIB = FACT*EVIB
      HVIB = FACT*HVIB
      GVIB = FACT*GVIB
      CVVIB = FACT*CVVIB
      CPVIB = FACT*CPVIB
      SVIB = FACT*SVIB
      ETOT = FACT*ETOT
      HTOT = FACT*HTOT
      GTOT = FACT*GTOT
      CVTOT = FACT*CVTOT
      CPTOT = FACT*CPTOT
      STOT = FACT*STOT
      IF (MASWRK) THEN
      WRITE(IW,9090)
      WRITE(IW,9100) '  KCAL/MOL','  KCAL/MOL','  KCAL/MOL',
     *               ' CAL/MOL-K',' CAL/MOL-K',' CAL/MOL-K'
      WRITE(IW,9110) 'ELEC. ',EELEC,HELEC,GELEC,CVELEC,CPELEC,SELEC
      WRITE(IW,9110) 'TRANS.',ETRAN,HTRAN,GTRAN,CVTRAN,CPTRAN,STRAN
      WRITE(IW,9110) 'ROT.  ',EROT,HROT,GROT,CVROT,CPROT,SROT
      WRITE(IW,9110) 'VIB.  ',EVIB,HVIB,GVIB,CVVIB,CPVIB,SVIB
      WRITE(IW,9110) 'TOTAL ',ETOT,HTOT,GTOT,CVTOT,CPTOT,STOT
      WRITE(IW,9115) ETHERM
C
C     ---- SAVE SOME OF ABOVE FOR COMPOSITE THERMOCHEMISTRY RUNS ----
C
      CXTHERM( 1) = EELEC
      CXTHERM( 2) = HELEC
      CXTHERM( 3) = ETRAN
      CXTHERM( 4) = HTRAN
      CXTHERM( 5) = EROT
      CXTHERM( 6) = HROT
      CXTHERM( 7) = EVIB
      CXTHERM( 8) = HVIB
      CXTHERM( 9) = ETOT
      CXTHERM(10) = HTOT
      CXTHERM(11) = GTOT
      CXZPE = ZPE
C
C        IF REQUESTED, WRITE OUT VIBRATIONAL ENTROPY CONTRIBUTION
C        OF EACH FREQUENCY IN KCAL/MOL-K
C
      IF(PRTSCN) THEN
         WRITE(IW,9055)
         WRITE(IW,9060) (TFREQ(I),ASVIB(I)*FACT,
     *           ASVIB(I)*FACT*100/SVIB,I=IFMIN,IFMAX)
         WRITE(IW,9061)
      END IF
      END IF
      IF(IPASS.LT.NPASS) GO TO 100
  800 CONTINUE
      RETURN
C
 9000 FORMAT(/5X,31(1H-)/5X,'THERMOCHEMISTRY AT T=',F8.2,' K'/
     *    5X,31(1H-)//
     *    1X,'USING IDEAL GAS, RIGID ROTOR, HARMONIC NORMAL MODE',
     *       ' APPROXIMATIONS.'/
     *    1X,'P=',1P,E13.5,' PASCAL.')
 9001 FORMAT(/5X,31(1H-)/5X,'THERMOCHEMISTRY AT T=',F8.2,' K'/
     *    5X,31(1H-)//
     *    1X,'USING IDEAL GAS, RIGID ROTOR, ANHARMONIC',
     *       ' APPROXIMATIONS.'/
     *    1X,'P=',1P,E13.5,' PASCAL.')
 9002 FORMAT(//1X,'PARTIAL HESSIAN ANALYSIS HAS BEEN DONE. ',
     *'FROZEN ATOMS AND FRAGMENTS',/
     *1X,'MAKE NO CONTRIBUTION TO ZERO POINT ENERGY, ',
     *'VIBRATIONAL PARTITION FUNCTION'/
     *1X,'AND THERMOCHEMICAL FUNCTIONS. BUT THEY STILL CONTRIBUTE ',
     *'TO TRANSLATIONAL'/
     *1X,'AND ROTATIONAL ENERGIES AND PARTITION',
     *' FUNCTIONS AS IF THEY WERE UNFROZEN.',//)
 9005 FORMAT(1X,'ALL FREQUENCIES ARE SCALED BY',2X,F8.5)
 9010 FORMAT(1X,'THE MOMENTS OF INERTIA ARE (IN AMU*BOHR**2)'/
     *       1X,3F12.5/
     *       1X,'THE ROTATIONAL SYMMETRY NUMBER IS',F5.1)
 9011 FORMAT(1X,'THE ROTATIONAL CONSTANTS ARE (IN GHZ)',/
     *       1X,3F12.5)
 9040 FORMAT(1X,I3,' IMAGINARY FREQUENCY VIBRATION(S) IS(ARE) EXCLUDED'/
     *       5X,'FROM THE FOLLOWING THERMOCHEMICAL ANALYSIS.')
 9050 FORMAT(1X,'THE HARMONIC ZERO POINT ENERGY IS (SCALED BY',
     *          F8.3,')'/
     *       1X,F15.6,' HARTREE/MOLECULE',
     *       1X,F15.6,' CM**-1/MOLECULE '/
     *       1X,F15.6,' KCAL/MOL        ',
     *       1X,F15.6,' KJ/MOL')
 9051 FORMAT(1X,'THE ANHARMONIC ZERO POINT ENERGY IS (SCALED BY',
     *          F8.3,')'/
     *       1X,F15.6,' HARTREE/MOLECULE',
     *       1X,F15.6,' CM**-1/MOLECULE '/
     *       1X,F15.6,' KCAL/MOL        ',
     *       1X,F15.6,' KJ/MOL')
 9055 FORMAT(11X,33(1H-)/11X,'VIBRATIONAL-ENTROPY DECOMPOSITION',/
     *    11X,33(1H-)//
     *    10X,'FREQUENCY',3X,'ENTROPY',3X,'%-CONTRIBUTION'/
     *    10X,9(1H-),3X,7(1H-),3X,14(1H-))
 9060 FORMAT(11X,F7.2,3X,F7.3,10X,F6.2)
 9061 FORMAT(11X,33(1H-))
 9065 FORMAT(/1X,'WARNING: THE HESSIAN MATRIX FOR AN EFFECTIVE',
     *           ' FRAGMENT RUN'/
     *        1X,'DOES NOT INCLUDE THE INTERNAL VIBRATIONS OF THE',
     *           ' FRAGMENTS.'/
     *        1X,'THIS MEANS THE ZERO POINT CONTRIBUTIONS OMIT THESE',
     *           ' FREQUENCIES,'/
     *        1X,'AND IN TURN THE THERMOCHEMICAL QUANTITIES OMIT THESE',
     *           ' TOO.')
 9070 FORMAT(/1X,14X,'Q',15X,'LN Q')
 9080 FORMAT(1X,A6,1P,E15.5,0P,F15.6)
 9090 FORMAT(/14X,'E',9X,'H',9X,'G',9X,'CV',8X,'CP',8X,'S')
 9100 FORMAT(1X,6X,6A10)
 9110 FORMAT(1X,A6,6F10.3)
 9115 FORMAT(1X,'E(T)-E(0) = H(T)-H(0) =',F10.3)
 9200 FORMAT(1X,'FREQUENCY NUMBER',I4,' IS UNEXPECTEDLY ZERO!')
      END
