C 30 APR 07 - SPW - NEW MODULES FOR NEO VIBRATIONAL ANALYSIS
C
C*MODULE NEOVIB  *DECK FGMTRX_NEO
      SUBROUTINE FGMTRX_NEO(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
      LOGICAL NUCOPT,CLSTEP,NEOTS
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
CSPW  COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
CSPW *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
CSPW *                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 /NEOOPT/ NUCOPT,CLSTEP,NEOTS
      COMMON /NUCMOI/ NUNIQN,IUNIQN(20),IUNIQT(20),NQMNUC,IQMNUC(20),
     *                IQNTYP(20),NUMNB,NUCST,NAUXNB,IAUXNB(20),NUMULT,
     *                NNA,NNB,NTAUXB
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /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 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
C     --- MASSES FOR QM NUCLEI NEED TO BE SET TO ZERO ---
C
      DO 95 I = 1,NQMNUC
          IAT = IQMNUC(I)
          ZMASS(IAT) = ZERO
   95 CONTINUE
      DO 96 I = 1,NTAUXB
          ZMASS(IAT) = ZERO
   96 CONTINUE
C
C     --- NOW PRINT OUT 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
         PROJCT = .FALSE.
         WRITE(IW,*) 'FOR NEO BEST NOT TO DO PROJECTION'
         WRITE(IW,*) 'SETTING PROJCT = .FALSE.'
      END IF
      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
CSPW  DO 560 IAT=1,NAT
CSPW     NUCZ = INT(ZAN(IAT)) + IZCORE(IAT)
CSPW     IF(NUCZ.EQ.0) NSKIP=NSKIP+3
C 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
C
      IF(NUCOPT .AND. MASWRK) WRITE(IW,9058)
      IF (MASWRK) WRITE(IW,9060) NFIRST,NLAST
C
C     --- FOR NUCOPT=.T. NEO QM NUCLEI ARE SET FROZEN FOR VIB ANALYSIS ---
C
      IF(NUCOPT) THEN
         II = 0
         DO I=1,NQMNUC
            IAT = IQMNUC(I)
            DO J=1,3
               II = II + 1
               ICRD = (IAT-1)*3 + J
               IFREEZ(II) = ICRD
            END DO
         END DO
C
         DO I=1,NTAUXB
            IAT = IAUXNB(I)
            DO J=1,3
               II = II + 1
               ICRD = (IAT-1)*3 + J
               IFREEZ(II) = ICRD
            END DO
         END DO
      END IF
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
CSPW  --- NEO DIPOLES ARE CORRECT YET SO SKIP INTENSITIES ---
      GOTDDM = .FALSE.
      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)
      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!'/)
 9058 FORMAT(/10X,'NEO QM NUCLEI ARE EXCLUDED FROM THE VIB ANALYSIS',
     *       /10X,'AS BY DEFINITION THEIR MODES HAVE ZERO FREQUENCIES.',
     *       /10X,'GAMESS ALSO TREATS THEM AS FROZEN COODINATES.')
 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)')
      END
