C 13 MAR 06 - NM  - ADJUST ARGS TO EGIN CALL
C 14 NOV 05 - JAB - CHANGE WAVEFUNCTION TESTING AND FORBID NUM. GRAD
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  7 MAR 05 - IA  - FIX COMMON BLOCK FRGINF
C 13 FEB 05 - MWS - PAD COMMON BLOCK FRGINF
C  7 SEP 04 - MWS - RAMAN: DELETE CALLS TO GAMGEN
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C 12 DEC 02 - DGF - RAMANX: COMPUTE MEMORY NEED BEFORE FG MATRX
C 17 APR 02 - MWS - SYNCH UP FRGINF COMMON
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C 25 JUN 01 - MWS - ALTER COMMON BLOCK SCFOPT AND WFNOPT
C 11 OCT 00 - MWS - PAD COMMON BLOCK PCMPAR
C 28 JUL 00 - MWS - RELOAD ERF TABLE IF USING POPLE INTS
C  1 MAY 00 - MWS - NEW MODULE TO COMPUTE RAMAN INTENSITY
C
C*MODULE RAMAN   *DECK RAMANX
      SUBROUTINE RAMANX
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL ISEPS,USEPS,MMONLY,QMMM,GOPARR,DSKWRK,MASWRK,
     *        GOTEG,GOTEH,GOTDDM,GOTADM,STATPT,SCFOK,
     *        DECOMP,PRTSCN,GOTFRQ,PROJCT,NOGRAD
C
      CHARACTER*18 BOOBOO
C
      DIMENSION EVALS(19),ALPHA(6)
C
      PARAMETER (MXATM=2000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      COMMON /ISEPS / ISEPS, USEPS
      COMMON /FMCOM / XXX(1)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PCMPAR/ IPCM,NFT26,NFT27,IRPPCM,IEF,IP_F,NFMOPCM
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /THERMD/ FREQ(3*MXATM),TEMP(10),SCLFAC,NTEMP,PRTSCN,GOTFRQ
      COMMON /TINOPT/ MMONLY,QMMM
      COMMON /VBDCMP/ DECOMP,NROW,NDEG
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
      COMMON /ZRFPAR/ GZRF,FIND(3),GNUCF,EBORN,DIELEC,IZRF,ICALC
C
      PARAMETER (TOLG=0.0005D+00)
C
      PARAMETER (NNAM=1)
      DIMENSION QNAM(NNAM),KQNAM(NNAM)
      DATA GNAM/8HRAMAN   /
      DATA QNAM/8HEFIELD  /
      DATA KQNAM/3/
C
      DATA RHF,GVB,RMC    /8HRHF     ,8HGVB     ,8HMCSCF   /
      DATA CIS,GUGA/8HCIS     ,8HGUGA    /
      DATA NONE,RNONE/4HNONE,8HNONE    /
C
C         ----- MAIN DRIVER FOR RAMAN SPECTRUM PREDICTION -----
C     THIS PROGRAM WAS WRITTEN BY MIKE SCHMIDT, IN APRIL 2000, AT
C     TOKYO METROPOLITAN UNIVERSITY, HACHIOJI-SHI, TOKYO-TO 192-0397
C
      EFIELD = 0.002D+00
      JRET=0
      CALL NAMEIO(IR,JRET,GNAM,NNAM,QNAM,KQNAM,
     *            EFIELD,
     *            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) THEN
         IF(MASWRK) WRITE(IW,*) 'PLEASE FIX $RAMAN INPUT ERROR'
         CALL ABRT
      END IF
C
C     THE METHOD REQUIRES THE ABILITY TO COMPUTE ANALYTIC GRADIENTS.
C     NO ONE HAS THOUGHT ABOUT THE FINITE DIFFERENCE NUCLEAR GRADIENTS,
C     WHAT DOES IT MEAN TO DISPLACE THE NUCLEI BY FINITE STEPS, IN
C     REGARDS TO THE DIPOLE INTEGRALS USED TO INDUCE THE FINITE FIELDS?
C
      NOGRAD = .FALSE.
      NOGRAD = NGLEVL.NE.0
     *         .OR.NOGRAD
      NOGRAD = MPLEVL.GT.0  .AND.  (SCFTYP.EQ.GVB .OR. SCFTYP.EQ.RMC)
     *         .OR.NOGRAD
      NOGRAD = CITYP.NE.RNONE .AND.
     *                       .NOT. ((SCFTYP.EQ.RHF.AND.CITYP.EQ.GUGA)
     *                        .OR.  (SCFTYP.EQ.RHF.AND.CITYP.EQ.CIS))
     *         .OR.NOGRAD
      NOGRAD = CCTYP.NE.RNONE
     *         .OR.NOGRAD
      IF(NOGRAD) THEN
         IF(MASWRK) WRITE(IW,9000) SCFTYP,MPLEVL,CITYP,CCTYP,NGLEVL
         CALL ABRT
      END IF
C
C        MOPAC DOESN'T HAVE DIPOLE INTEGRALS
C
      IF(MPCTYP.NE.NONE) THEN
         IF(MASWRK) WRITE(IW,9010)
         CALL ABRT
      END IF
C
C        WE DON'T WANT TO THINK ABOUT THE PROGRAMMING CONSEQUENCES OF
C        SOLVENT MODELING, SO JUST MAKE ANY MEDIA EFFECT MODEL BOMB.
C
C        IT IS POSSIBLE THAT SOME OF THESE COULD BE MADE TO WORK IF
C        THEY WERE THOUGHT THROUGH, AS MOST HAVE GRADIENTS PROGRAMMED.
C        FOR EXAMPLE, EFFECTIVE FRAGMENTS WOULD NEED TO HAVE EXTRA
C        STORAGE PROVIDED IN GRADIENT AND HESSIAN FOR THEIR COMPONENTS.
C        SOMETHING TO PONDER: WHAT DO APPLIED FIELDS DO TO DIELECTRIC
C        CONTINUUM MODELS?
C
      IF(QMMM)      BOOBOO='TINKER MM'
      IF(ISEPS)     BOOBOO='COSMO SOLVATION'
      IF(IPCM.NE.0) BOOBOO='PCM SOLVATION'
      IF(IZRF.NE.0) BOOBOO='SCRF SOLVATION'
      IF(NFRG.GT.0) BOOBOO='EFFECTIVE FRAGMENT'
      IF(QMMM  .OR.  ISEPS  .OR.  IPCM.NE.0  .OR.  IZRF.NE.0
     *         .OR.  NFRG.GT.0) THEN
         IF(MASWRK) WRITE(IW,9020) BOOBOO
         CALL ABRT
      END IF
C
C        NUCLEAR HESSIAN MUST HAVE BEEN COMPUTED IN AN EARLIER RUN,
C        AND PROVIDED IN THE INPUT STREAM TO THIS JOB.  IF NOT, BOMB.
C
      CALL VALFM(LOADFM)
      LHESS = LOADFM + 1
      LAST  = LHESS  + 3*NAT*3*NAT
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      CALL FCMIN(XXX(LHESS),3*NAT,GOTEH)
      IF(GOTEH) THEN
         CALL FCMOUT(XXX(LHESS),3*NAT)
      ELSE
         IF(MASWRK) WRITE(IW,9030)
         CALL ABRT
      END IF
      CALL RETFM(NEED)
C
C        ----- POSSIBLE $GRAD, $DIPDR, $ALPDR INPUT GROUPS -----
C        EACH INPUT ROUTINE ZEROS IF NOT FOUND, AND WRITES EITHER
C        THAT ZERO VECTOR OR ELSE THE INPUT VALUES TO DAF RECORDS.
C
      CALL VALFM(LOADFM)
      LGRAD = LOADFM + 1
      LDIPD = LGRAD  + 3*NAT*19
      LALPD = LDIPD  + 3*NAT*3
      LAST  = LALPD  + 3*NAT*6
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      CALL EGIN (XXX(LGRAD),3*NAT,GOTEG,' $GRAD  ')
      CALL DDMIN(XXX(LDIPD),3*NAT,GOTDDM)
      CALL ADMIN(ALPHA,XXX(LALPD),3*NAT,GOTADM)
      SCFOK=.TRUE.
      IF(GOTDDM  .AND.  GOTADM) GO TO 200
C
C        COMPUTE ALPHA, ALPHA DERIVATIVE, AND DIPOLE DERIVATIVE MATRICES
C
      CALL RAMAN(EVALS,XXX(LGRAD),XXX(LDIPD),ALPHA,XXX(LALPD),
     *           NAT,SCFOK,GOTEG,GOTDDM,EFIELD)
      CALL DAWRIT(IDAF,IODA,XXX(LGRAD),  3*NAT,  3,0)
      CALL DAWRIT(IDAF,IODA,XXX(LDIPD),3*3*NAT, 34,0)
      CALL DAWRIT(IDAF,IODA,XXX(LALPD),6*3*NAT,249,0)
C
  200 CONTINUE
      GMAX = XXX(LGRAD-1+IDAMAX(3*NAT,XXX(LGRAD),1))
      GRMS = SQRT(DDOT(3*NAT,XXX(LGRAD),1,XXX(LGRAD),1)/(3*NAT))
      STATPT = GMAX.LT.TOLG  .AND.  GRMS.LT.(TOLG/3.0D+00)
C
      IF(MASWRK) WRITE(IW,9040)
      CALL EGOUT (XXX(LGRAD),NAT)
      CALL DDMOUT(XXX(LDIPD),3*NAT)
      CALL DDMPUN(XXX(LDIPD),3*NAT)
      CALL POLOUT(ALPHA)
      CALL ADMOUT(XXX(LALPD),3*NAT)
      CALL ADMPUN(EVALS,ALPHA,XXX(LALPD),3*NAT)
      CALL RETFM(NEED)
C
C        CARRY OUT NORMAL COORDINATE ANALYSIS TO PRINT IR/RAMAN
C
      NFCM = (3*NAT*3*NAT+3*NAT)/2
      CALL VALFM(LOADFM)
      LVEC  = LOADFM + 1
      LFCM  = LVEC   + 3*NAT*3*NAT
      LE    = LFCM   + NFCM
      LSCR  = LE     + 3*NAT
      LIA   = LSCR   + 3*NAT*8
      LRM   = LIA    + 3*NAT
      LSVT  = LRM    + 3*NAT
      LSVR  = LSVT   + 3*NAT*3
      LSVTT = LSVR   + 3*NAT*3
      LSVRT = LSVTT  + 3*NAT
      LCC   = LSVRT  + 3*NAT
      LCOM  = LCC    + 3*NAT
      LZMS  = LCOM   + 3*NAT
      LBUF1 = LZMS   + NAT
      LBUF2 = LBUF1  + 3*NAT*3*NAT
      LDIPD = LBUF2  + 3*NAT*3*NAT
      LALPD = LDIPD  + 3*NAT*3
      LAST  = LALPD  + 3*NAT*6
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
C
      PROJCT=.FALSE.
      DECOMP=.FALSE.
      NTEMP=1
      TEMP(1)=298.15D+00
      SCLFAC=1.0D+00
      PRTSCN=.FALSE.
      GOTFRQ=.FALSE.
      CALL FGMTRX(XXX(LVEC),XXX(LFCM),XXX(LE),XXX(LSCR),XXX(LIA),
     *            XXX(LRM),XXX(LSVT),XXX(LSVR),XXX(LSVTT),XXX(LSVRT),
     *            XXX(LCC),XXX(LCOM),XXX(LZMS),XXX(LBUF1),XXX(LBUF2),
     *            3*NAT,NFCM,NAT,XXX(LDIPD),.TRUE.,XXX(LALPD),.TRUE.,
     *            PROJCT,STATPT,SCFOK)
      CALL RETFM(NEED)
C
      IF(MASWRK) WRITE(IW,9050)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/1X,'SCFTYP=',A8,' WITH EITHER MPLEVL=',I2,' OR CITYP=',A8,
     *           ' OR CCTYP=',A8/
     *       1X,'HAS NO ANALYTIC GRADIENT. (YOUR NUM.GRADIENT FLAG=',
     *          I4,')'/
     *       1X,'THEREFORE, THE POLARIZABILITY DERIVATIVE CANNOT BE',
     *          ' COMPUTED.')
 9010 FORMAT(1X,'MOPAC HAS NO DIPOLE INTEGRALS, SO ELECTRIC FIELDS',
     *          ' CANNOT BE APPLIED.'/
     *       1X,'THEREFORE, THE POLARIZABILITY DERIVATIVE CANNOT BE',
     *          ' COMPUTED.')
 9020 FORMAT(/1X,'RAMAN SPECTRUM IS PROGRAMMED ONLY FOR PURE AB INITIO',
     *          ' RUNS.'/
     *       1X,'PLEASE REMOVE ',A,' MODEL FROM YOUR INPUT.')
 9030 FORMAT(/1X,'RAMAN SPECTRUM REQUIRES A PREVIOUSLY COMPUTED $HESS'/
     *        1X,'NUCLEAR FORCE CONSTANT MATRIX AS INPUT TO THIS RUN.'/
     *        1X,'PLEASE OBTAIN THIS MATRIX BY A RUNTYP=HESSIAN JOB.')
 9040 FORMAT(/10X,27(1H-)/10X,'ZERO FIELD NUCLEAR GRADIENT'/
     *        10X,27(1H-))
 9050 FORMAT(/1X,'..... DONE WITH NORMAL COORDINATE ANALYSIS .....')
      END
C*MODULE RAMAN   *DECK RAMAN
      SUBROUTINE RAMAN(EVALS,GRADS,DIPDER,ALPHA,ALPDER,
     *                 NAT,SCFOK,GOTEG,GOTDDM,EFIELD)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000, MXSH=5000)
      LOGICAL SCFOK,GOTEG,GOTDDM,EFLDL,GOPARR,DSKWRK,MASWRK,
     *        MCCI,SOME
      DIMENSION EVALS(19),GRADS(3*NAT,19),DIPDER(3,3*NAT),
     *          ALPHA(6),ALPDER(6,3*NAT)
      DIMENSION DF(2)
      CHARACTER*1 FIELD(3)
C
      COMMON /EFLDC / EVEC(3),EFLDL
      COMMON /FMCOM / XXX(1)
      COMMON /FUNCT / ENERGY,EG(3,MXATM)
      COMMON /INFOA / NATM,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SCFOPT/ CONVHF,MAXIT,MCONV,NPUNCH,NPREO(4)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
      PARAMETER (SMALL=1.0D-08)
      DATA CHECK/8HCHECK   /, GUGA/8HGUGA    /, RMC/8HMCSCF   /
      DATA FIELD/'X','Y','Z'/
C
C     --- DRIVER TO COMPUTE THE POLARIZABILITY DERIVATIVE MATRIX ---
C     TECHNIQUE IS FINITE DIFFERENCING, UNDER VARIOUS APPLIED FIELDS.
C     SEE A.KOMORNICKI, J.W.MCIVER J.CHEM.PHYS. 70, 2014-2016(1979)
C     D.R.FREDKIN, A.KOMORNICKI, S.R.WHITE, K.R.WILSON
C        J.CHEM.PHYS. 78, 7077-7092(1983)
C     G.B.BACSKAY, S.SAEBO, P.R.TAYLOR CHEM.PHYS. 90, 215-224(1984)
C
      SOME = NPRINT.NE.-5  .AND.  MASWRK
      IF(SOME) WRITE(IW,9000) EFIELD
      MCCI = SCFTYP.EQ.RMC  .OR.  CITYP.EQ.GUGA
      SCFOK=.TRUE.
      L1 = NUM
      L2 = (L1*L1+L1)/2
      L3 = L1*L1
C
C        OBTAIN ENERGY AND GRADIENT WITHOUT APPLIED FIELD
C        GOOD ACCURACY IN WAVEFUNCTIONS GIVES MORE ACCURATE DERIVATIVES
C
      IF(CONVHF.GT.1.0D-06) CONVHF=1.0D-06
      EFLDL = .FALSE.
      EVEC(1) = ZERO
      EVEC(2) = ZERO
      EVEC(3) = ZERO
      ENERGY = ZERO
      CALL VCLR(EG,1,3*NAT)
C
C         OBTAIN FIELD FREE WAVEFUNCTION
C
      IF(SOME) WRITE(IW,9010)
      CALL ENERGX
C
      IF(ENERGY.EQ.ZERO .AND. EXETYP.NE.CHECK) THEN
         IF (SOME) WRITE(IW,9020)
         CALL ABRT
      END IF
C
C         SAVE FIELD FREE ORBITALS AND PRISTINE H MATRIX
C
      CALL VALFM(LOADFM)
      LWRK = LOADFM+1
      LAST = LWRK + MAX(L2,L3)
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      CALL FFSVMO(0,XXX(LWRK))
      IF(MCCI) CALL FFSVCI(0)
      CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,11,0)
      CALL DAWRIT(IDAF,IODA,XXX(LWRK),L2,94,0)
      CALL RETFM(NEED)
C
C         OBTAIN FIELD FREE GRADIENT, IF IT WAS NOT READ IN
C         WE MUST EXECUTE THE DISTRIBUTED MP2 GRADIENT CODE TO FREE
C         A DISTRIBUTED MATRIX, NO MATTER WHAT.
C
      IF(MPLEVL.EQ.2  .AND.  GOPARR) GOTEG=.FALSE.
      IF(GOTEG) THEN
         IF(SOME) WRITE(IW,9025)
         CALL DCOPY(3*NAT,GRADS,1,EG,1)
      ELSE
         IF (CITYP.EQ.GUGA) THEN
            CALL CIGRAD
         ELSE
            CALL HFGRAD
         END IF
      END IF
      CALL EGOUT(EG,NAT)
      CALL DAWRIT(IDAF,IODA,EG,3*NAT,3,0)
      ISAVE=1
      EVALS(ISAVE) = ENERGY
      CALL DCOPY(3*NAT,EG,1,GRADS(1,ISAVE),1)
      IF(SOME) WRITE(IW,9027)
      CALL TIMIT(1)
C
C         LOOP OVER APPLIED ELECTRIC FIELDS
C
      DF(1) = +EFIELD
      DF(2) = -EFIELD
C
      NPRINT = -5
      NPUNCH = 0
      EFLDL=.TRUE.
      NTSAVE = NT
      CALL SYMOFF
      IF(NTSAVE.GT.1  .AND.  .NOT.MCCI) THEN
         IF(SOME) WRITE(IW,9030)
         CALL JANDK
      END IF
C
C         THE PROPERTY EVALUATION FOR THE FIELD-FREE WAVEFUNCTION
C         HAS LEFT BEHIND CENTER OF MASS ORIGIN DIPOLE INTEGRALS.
C         BELOW WE NEED THEM WITH RESPECT TO THE COORDINATE ORIGIN.
C
      CALL DIPINT(ZERO,ZERO,ZERO,.FALSE.)
C
C             FIELDS FOR X, Y, Z DERIVATIVES - NOTE DOUBLE STEPSIZE
C
      DO 390 IFIELD=1,3
         DO 380 I=1,2
            CALL VALFM(LOADFM)
            LWRK = LOADFM + 1
            LH   = LWRK   + MAX(L2,L3)
            LAST = LH     + L2
            NEED = LAST - LOADFM - 1
            CALL GETFM(NEED)
            CALL FFSVMO(1,XXX(LWRK))
            IF(MCCI) CALL FFSVCI(1)
            CALL ORTHDN
C
            IF(SOME) WRITE(IW,9040) FIELD(IFIELD),TWO*DF(I)
            EVEC(IFIELD) = TWO*DF(I)
C
            CALL DAREAD(IDAF,IODA,XXX(LH),L2,94,0)
            IF(EVEC(1).NE.ZERO) THEN
               CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,95,0)
               CALL DAXPY(L2,EVEC(1),XXX(LWRK),1,XXX(LH),1)
            END IF
            IF(EVEC(2).NE.ZERO) THEN
               CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,96,0)
               CALL DAXPY(L2,EVEC(2),XXX(LWRK),1,XXX(LH),1)
            END IF
            IF(EVEC(3).NE.ZERO) THEN
               CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,97,0)
               CALL DAXPY(L2,EVEC(3),XXX(LWRK),1,XXX(LH),1)
            END IF
            CALL DAWRIT(IDAF,IODA,XXX(LH),L2,11,0)
            CALL RETFM(NEED)
C
            CALL WFN
            IF(ENERGY.EQ.ZERO) SCFOK=.FALSE.
            IF (CITYP.EQ.GUGA) THEN
               CALL CIGRAD
            ELSE
               CALL HFGRAD
            END IF
C
            EVEC(IFIELD) = ZERO
            ISAVE=ISAVE+1
            EVALS(ISAVE) = ENERGY
            CALL DCOPY(3*NAT,EG,1,GRADS(1,ISAVE),1)
            IF(SOME) WRITE(IW,9027)
            CALL TIMIT(1)
  380    CONTINUE
  390 CONTINUE
C
C             CROSS FIELDS FOR XY, XZ, YZ PARTIALS
C
      DO 490 IFIELD=1,3
         DO 480 JFIELD=1,IFIELD
            IF(JFIELD.EQ.IFIELD) GO TO 480
            DO 470 I=1,2
               DO 460 J=1,2
                  CALL VALFM(LOADFM)
                  LWRK = LOADFM + 1
                  LH   = LWRK   + MAX(L2,L3)
                  LAST = LH     + L2
                  NEED = LAST - LOADFM - 1
                  CALL GETFM(NEED)
                  CALL FFSVMO(1,XXX(LWRK))
                  IF(MCCI) CALL FFSVCI(1)
                  CALL ORTHDN
C
                  IF(SOME) WRITE(IW,9050) FIELD(IFIELD),DF(I),
     *                                    FIELD(JFIELD),DF(J)
                  EVEC(IFIELD) = DF(I)
                  EVEC(JFIELD) = DF(J)
C
                  CALL DAREAD(IDAF,IODA,XXX(LH),L2,94,0)
                  IF(EVEC(1).NE.ZERO) THEN
                     CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,95,0)
                     CALL DAXPY(L2,EVEC(1),XXX(LWRK),1,XXX(LH),1)
                  END IF
                  IF(EVEC(2).NE.ZERO) THEN
                     CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,96,0)
                     CALL DAXPY(L2,EVEC(2),XXX(LWRK),1,XXX(LH),1)
                  END IF
                  IF(EVEC(3).NE.ZERO) THEN
                     CALL DAREAD(IDAF,IODA,XXX(LWRK),L2,97,0)
                     CALL DAXPY(L2,EVEC(3),XXX(LWRK),1,XXX(LH),1)
                  END IF
                  CALL DAWRIT(IDAF,IODA,XXX(LH),L2,11,0)
                  CALL RETFM(NEED)
C
                  CALL WFN
                  IF (CITYP.EQ.GUGA) THEN
                     CALL CIGRAD
                  ELSE
                     CALL HFGRAD
                  END IF
C
                  EVEC(IFIELD) = ZERO
                  EVEC(JFIELD) = ZERO
                  ISAVE=ISAVE+1
                  EVALS(ISAVE) = ENERGY
                  CALL DCOPY(3*NAT,EG,1,GRADS(1,ISAVE),1)
                  IF(SOME) WRITE(IW,9027)
                  CALL TIMIT(1)
  460          CONTINUE
  470       CONTINUE
  480    CONTINUE
  490 CONTINUE
C
      EFLDL = .FALSE.
      CALL SYMON
C
C     CENTRAL DIFFERENCING FORMULAE ARE USED FOR ALL QUANTITIES BELOW
C
C     ----- DIPOLE DERIVATIVE MATRIX -----
C     QUANTITY IS ENERGY DERIVATIVE WRT FIELD AND NUCLEAR COORDINATE,
C     OBTAIN BY NUMERICAL FIELD DERIVATIVE OF ANALYTIC NUCLEAR GRADIENT.
C     THIS TENSOR MUST BE IN UNITS OF DEBYE/ANGSTROM ELSEWHERE.
C     IF ONE WAS READ IN, IT MIGHT EVEN BE ANALYTIC, AND CERTAINLY ITS
C     NO WORSE TO USE NUMERICAL NUCLEAR DERIVATIVE OF ANALYTIC DIPOLES
C     INSTEAD OF NUMERICAL FIELD DERIVATIVE OF ANALYTIC NUCLEAR GRADIENT
C
      IF(GOTDDM) THEN
         IF(SOME) WRITE(IW,9060)
      ELSE
         DO 530 IA=1,NAT*3
            X = GRADS(IA,2) - GRADS(IA,3)
            Y = GRADS(IA,4) - GRADS(IA,5)
            Z = GRADS(IA,6) - GRADS(IA,7)
            IF(ABS(X).LT.SMALL) X=ZERO
            IF(ABS(Y).LT.SMALL) Y=ZERO
            IF(ABS(Z).LT.SMALL) Z=ZERO
            DIPDER(1,IA) = X
            DIPDER(2,IA) = Y
            DIPDER(3,IA) = Z
  530    CONTINUE
         DENOM = TWO*TWO*EFIELD
         CALL DSCAL(3*3*NAT,-ONE/DENOM,DIPDER,1)
         TOANGS = 0.52917724924D+00
         DEBYE = 2.541766D+00
         UNITS = DEBYE/TOANGS
         CALL DSCAL(3*3*NAT,UNITS,DIPDER,1)
      END IF
C
C     ----- ALPHA POLARIZABILITY TENSOR -----
C     QUANTITY IS ENERGY DERIVATIVE WRT FIELD TWICE
C     OBTAIN BY 2ND NUMERICAL FIELD DERIVATIVE OF THE ENERGY.
C     WE DON'T ACTUALLY NEED THIS FOR RAMAN SPECTRUM,
C     BUT ITS AN INTERESTING QUANTITY TO OBTAIN ANYWAY.
C
      XX = EVALS(2) - TWO * EVALS(1) + EVALS(3)
      YY = EVALS(4) - TWO * EVALS(1) + EVALS(5)
      ZZ = EVALS(6) - TWO * EVALS(1) + EVALS(7)
      XY = EVALS( 8) - EVALS( 9) - EVALS(10) + EVALS(11)
      XZ = EVALS(12) - EVALS(13) - EVALS(14) + EVALS(15)
      YZ = EVALS(16) - EVALS(17) - EVALS(18) + EVALS(19)
      IF(ABS(XX).LT.SMALL) XX=ZERO
      IF(ABS(XY).LT.SMALL) XY=ZERO
      IF(ABS(YY).LT.SMALL) YY=ZERO
      IF(ABS(XZ).LT.SMALL) XZ=ZERO
      IF(ABS(YZ).LT.SMALL) YZ=ZERO
      IF(ABS(ZZ).LT.SMALL) ZZ=ZERO
      ALPHA(1) = XX
      ALPHA(2) = XY
      ALPHA(3) = YY
      ALPHA(4) = XZ
      ALPHA(5) = YZ
      ALPHA(6) = ZZ
      DENOM = TWO*TWO*EFIELD*EFIELD
      CALL DSCAL(6,-ONE/DENOM,ALPHA,1)
C
C     ----- ALPHA DERIVATIVE TENSOR -----
C     QUANTITY IS ENERGY DERIVATIVE WRT FIELD TWICE, NUCLEAR COORDS ONCE
C     OBTAIN BY 2ND NUM. FIELD DERIVATIVE OF ANALYTIC NUCLEAR GRADIENT.
C
      DO 580 IA=1,NAT*3
         XX = GRADS(IA,2) - TWO * GRADS(IA,1) + GRADS(IA,3)
         YY = GRADS(IA,4) - TWO * GRADS(IA,1) + GRADS(IA,5)
         ZZ = GRADS(IA,6) - TWO * GRADS(IA,1) + GRADS(IA,7)
         XY = GRADS(IA, 8) - GRADS(IA, 9) - GRADS(IA,10) + GRADS(IA,11)
         XZ = GRADS(IA,12) - GRADS(IA,13) - GRADS(IA,14) + GRADS(IA,15)
         YZ = GRADS(IA,16) - GRADS(IA,17) - GRADS(IA,18) + GRADS(IA,19)
         IF(ABS(XX).LT.SMALL) XX=ZERO
         IF(ABS(XY).LT.SMALL) XY=ZERO
         IF(ABS(YY).LT.SMALL) YY=ZERO
         IF(ABS(XZ).LT.SMALL) XZ=ZERO
         IF(ABS(YZ).LT.SMALL) YZ=ZERO
         IF(ABS(ZZ).LT.SMALL) ZZ=ZERO
         ALPDER(1,IA) = XX
         ALPDER(2,IA) = XY
         ALPDER(3,IA) = YY
         ALPDER(4,IA) = XZ
         ALPDER(5,IA) = YZ
         ALPDER(6,IA) = ZZ
  580 CONTINUE
      DENOM = TWO*TWO*EFIELD*EFIELD
      CALL DSCAL(6*3*NAT,-ONE/DENOM,ALPDER,1)
C
      IF(SOME) WRITE(IW,9070)
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/5X,57(1H-)/5X,'NUMERICAL COMPUTATION OF POLARIZABILITY',
     *          ' DERIVATIVE MATRIX'/5X,57(1H-)//
     *        5X,'APPLIED ELECTRIC FIELD STRENGTH=',F8.4,' A.U.')
 9010 FORMAT(/1X,'OBTAINING GRADIENT IN THE ABSENCE OF APPLIED FIELDS')
 9020 FORMAT(1X,'*** FIELD-FREE SCF CALCULATION DID NOT CONVERGE ***'/
     *       1X,'TO FIX, PLEASE OBTAIN BETTER STARTING ORBITALS, ETC.')
 9025 FORMAT(/1X,'USING INPUT $GRAD AS THE FIELD FREE GRADIENT VECTOR.')
 9027 FORMAT(1X,'... DONE WITH GRADIENT ...')
 9030 FORMAT(/1X,'RECOMPUTING 2E- INTEGRALS WITHOUT SYMMETRY FOR',
     *           ' APPLIED FIELD COMPUTATIONS.')
 9040 FORMAT(/1X,'OBTAINING GRADIENT FOR ',A1,' FIELD AT STRENGTH',
     *          F6.3/)
 9050 FORMAT(/1X,'OBTAINING GRADIENT FOR ',A1,' FIELD AT STRENGTH',
     *          F6.3,' AND ',A1,' FIELD AT',F6.3/)
 9060 FORMAT(/1X,'ACCEPTING INPUT DIPOLE DERIVATIVE TENSOR INSTEAD',
     *           ' OF COMPUTING IT NUMERICALLY')
 9070 FORMAT(/1X,'..... DONE COMPUTING RAMAN DERIVATIVES .....')
      END
C*MODULE RAMAN   *DECK ADMIN
      SUBROUTINE ADMIN(ALPHA,ALPDER,NCOORD,GOTADM)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,GOTADM
C
      DIMENSION ALPHA(6),ALPDER(6,NCOORD)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- READ ALPHA POLARIZABILITY DERIVATIVE MATRIX FROM CARDS -----
C
      GOTADM = .FALSE.
      CALL VCLR(ALPDER,1,6*NCOORD)
C
C           POSITION INPUT FILE TO $ALPDR
C
      CALL SEQREW(IR)
      CALL FNDGRP(IR,' $ALPDR ',IEOF)
      IF (IEOF.EQ.1) GO TO 400
      GOTADM=.TRUE.
C
      IF (MASWRK) THEN
         READ(IR,*)
         READ(IR,*)
         READ(IR,9010) (ALPHA(I),I=1,6)
         READ(IR,*)
         DO 240 J=1,NCOORD
            READ(IR,9010) (ALPDER(I,J),I=1,6)
  240    CONTINUE
      END IF
C
C SEND ALPHA AND ALPDER TO ALL SLAVES
C
      IF (GOPARR) THEN
         CALL DDI_BCAST(450,'F',ALPHA ,6       ,MASTER)
         CALL DDI_BCAST(451,'F',ALPDER,6*NCOORD,MASTER)
      END IF
C
  400 CONTINUE
      CALL DAWRIT(IDAF,IODA,ALPHA ,6       ,251,0)
      CALL DAWRIT(IDAF,IODA,ALPDER,6*NCOORD,249,0)
      RETURN
C
 9010 FORMAT(6E13.6)
      END
C*MODULE RAMAN   *DECK ADMOUT
      SUBROUTINE ADMOUT(ALPDER,NCOORD)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION ALPDER(6,NCOORD)
C
      PARAMETER (MXATM=2000, MXAO=8192)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
C     ----- PRINT POLARIZABILITY DERIVATIVE MATRIX -----
C
      IF (.NOT.MASWRK) RETURN
C
      WRITE(IW,9040)
      DO 480 J=1,NAT
         WRITE(IW,9060) ANAM(J),BNAM(J),(ALPDER(I,3*J-2),I=1,6)
         WRITE(IW,9080)                 (ALPDER(I,3*J-1),I=1,6)
         WRITE(IW,9100)                 (ALPDER(I,3*J  ),I=1,6)
  480 CONTINUE
      RETURN
C
 9040 FORMAT(/10X,38(1H-),5X,7(1H-)/
     *        10X,'ALPHA POLARIZABILITY DERIVATIVE TENSOR',5X,'BOHR**2'/
     *        10X,38(1H-),5X,7(1H-)//
     *       3X,'ATOM',14X,'ALP-XX',4X,'ALP-XY',4X,'ALP-YY',4X,
     *          'ALP-XZ',4X,'ALP-YZ',4X,'ALP-ZZ')
 9060 FORMAT(1X,A8,A2,' D/DX ',6F10.5)
 9080 FORMAT(11X,     ' D/DY ',6F10.5)
 9100 FORMAT(11X,     ' D/DZ ',6F10.5)
      END
C*MODULE RAMAN   *DECK ADMPUN
      SUBROUTINE ADMPUN(ENERGY,ALPHA,ALPDER,NCOORD)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION ALPHA(6),ALPDER(6,NCOORD)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C      ----- PUNCH THE POLARIZABILITY DERIVATIVE MATRIX -----
C
      IF (.NOT.MASWRK) RETURN
C
      WRITE(IP,9001) ENERGY
      WRITE(IP,9002) (ALPHA(I),I=1,6)
      WRITE(IP,9003)
      DO 10 J=1,NCOORD
         WRITE(IP,9002) (ALPDER(I,J),I=1,6)
   10 CONTINUE
      WRITE(IP,9004)
      RETURN
C
 9001 FORMAT(' $ALPDR'/'E=',F20.10/'ALPHA POLARIZABILITY IN A.U.')
 9002 FORMAT(1P,6E13.6)
 9003 FORMAT('ALPHA POLARIZABILITY NUCLEAR DERIVATIVE IN A.U.')
 9004 FORMAT(' $END')
      END
