C 14 NOV 05 - DGF - PAD COMMON BLOCK ENRGYS
C 19 SEP 05 - IA  - SYNCHRONIZE FRGINF COMMON
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - DGF - PLTORB: PUNCH MORE FOR FMO RUNS
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 NSHEL, FRGINF
C  5 FEB 05 - MWS - PAD COMMON MCINP
C  7 SEP 04 - MWS - PAD COMMON BLOCK INTFIL
C 22 DEC 03 - MWS - AIMPAC: FIX FOR CIS RUNS
C  9 DEC 03 - MWS - SYNCH COMMON BLOCK RUNOPT
C  3 SEP 03 - FPA - AIMPAC: AVOID HUGE EIGENVALS, ADD F PERMUTATION
C 26 MAR 03 - MWS - AIMPAC: FIX GVB ORBITAL COUNTER
C 28 JAN 03 - JI  - ADD HOOK FOR ORMAS CI STEP
C 12 DEC 02 - CMA - SYNCH UP MP2PAR COMMON
C  7 AUG 02 - MWS - AIMPAC: WRITE TRAILER BEHIND AIMPAC OUTPUT
C 22 MAY 02 - MWS - NEW CALLING ARG FOR UHFNOS
C 17 APR 02 - MWS - SYNCH UP EFMULT AND FRGINF COMMON
C 26 MAR 02 - MWS - AIMPAC: CORRECT CITYP AND NORMALIZATION FOR F,G
C 24 JAN 02 - BTA - AIMPAC: REGENERATE MP2 NOS, IF APPROPRIATE
C 16 NOV 01 - JAB - PLTORB: ADD LINE FOR PLOTTING PLANE DISPLACEMENT
C  6 SEP 01 - MWS - ADD DUMMY ARGUMENTS TO NAMEIO CALL
C  1 AUG 01 - JI  - TWEAKS FOR GENERAL CI
C 25 JUN 01 - MWS - ALTER COMMON BLOCK WFNOPT
C 20 FEB 01 - FPA - AIMPAC: FORMAT CHANGES, AND THE INFAMOUS RV->RU
C 11 OCT 00 - MWS - UPDATE THE DETWFN COMMON
C 21 DEC 99 - MWS - MAKE GUGWFN COMMON CONSISTENT
C 13 FEB 99 - MWS - AIMPAC: FIX CALL FOR UHF NATURAL ORBITALS
C 12 NOV 98 - MWS - KEEP COMMON MCINP CONSISTENT
C 12 SEP 98 - TLW - AOGRID: CORRECT F,G NORMALIZATIONS
C  6 MAY 98 - MWS - AIMPAC: COUNT OCCUPIED ORBITALS FOR CI/MCSCF
C 28 SEP 97 - MWS - PLT3D: DROP PARALLEL INITIALIZATIONS
C  2 SEP 97 - MWS - MOLPLT: FIX FOR FRAGMENT ONLY RUNS
C 27 MAR 97 - BT  - PLT3D: CORRECT OPEN CALL ARGUMENTS
C 17 SEP 96 - MWS - MOLPLT: CORRECT COMMON DECLARATION
C 18 SEP 96 - MWS - MOLMEM,MOLPLT: CHANGES FOR EFFECTIVE FRAGMENTS
C 13 JUN 96 - VAG - ADD VARIABLE FOR CI TYPE TO SCFOPT COMMON
C 18 MAY 96 - MWS - PLTMEM: REMOVE SOME UNINTERESTING OUTPUT.
C 27 MAR 96 - SK  - PLTORB: 3D MOS TO INTERFACE WITH MOLGRAPH
C  5 MAR 96 - MWS - CHANGE SHELL SYMMETRY PACKING COMMONS
C 26 JUL 95 - VAG - PLTORB: ALLOW FOR F AND G FUNCTIONS
C  6 APR 95 - MWS - ZEALOT: USE FORMAT STATEMENTS
C 12 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 17 DEC 93 - MWS - MOLPLT,PLTORB: CHANGE GTBOND'S CALL, COLOR ATOMS
C 21 SEP 93 - MWS - ZEALX: USE REAL FORMAT STATEMENTS
C 16 JUL 93 - MWS - INCREASE MAXIMUM CI ROOTS TO 100
C 10 NOV 92 - MK  - ADD ZEALOT ROUTINES
C 24 MAY 92 - MWS - FIX AIMPAC AND RPAC NUCLEAR CHARGES IF ECP
C  2 APR 92 - MWS,TLW - COMMON ENRGYS MADE PURE FLOATING POINT
C 12 MAR 92 - MWS - REDIMENSION TO 500 ATOMS
C  9 MAR 92 - JHJ - PLTORB FILES MAY SAY SLATER
C 10 JAN 92 - TLW - PARALLEL I/O CHANGES
C 14 NOV 91 - TB  - RPAC15: DAY SHOULD BE C*9 VARIABLE
C  7 AUG 90 - TLW - ADD CF AND CG TO COMMON NSHEL
C 25 JUL 90 - MWS - NEW MODULE OF INTERFACE ROUTINES CREATED
C
C*MODULE PARLEY  *DECK AIMMEM
      SUBROUTINE AIMMEM
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000)
C
      COMMON /FMCOM / X(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
C
      DATA CHECK/8HCHECK   /
C
      NPRIMS = 0
      DO 130 I = 1,NSHELL
         NPRIMS = NPRIMS + KNG(I) * (KMAX(I)-KMIN(I)+1)
  130 CONTINUE
C
      L0 = NQMT
      L1 = NUM
      L2 = (NUM*NUM+NUM)/2
      L3 = NUM*NUM
C
      CALL VALFM(LOADFM)
      LACO   = LOADFM + 1
      LAIC   = LACO   + NPRIMS
      LEXPON = LAIC   + NPRIMS
      LICENT = LEXPON + NPRIMS
      LITYPE = LICENT + NPRIMS
      LOE    = LITYPE + NPRIMS
      LOCCNO = LOE    + L1
      LVEC   = LOCCNO + L1
      LIWRK  = LVEC   + L3
      LSCR   = LIWRK  + L1
      LDA    = LSCR   + L1*8
      LDB    = LDA    + L2
      LSPHF  = LDB    + L2
      LSMHF  = LSPHF  + L3
      LAST   = LSMHF  + L3
      NEED = LAST - LOADFM - 1
      CALL GETFM(NEED)
      IF(EXETYP.EQ.CHECK) GO TO 200
C
      CALL AIMPAC(X(LACO),X(LAIC),X(LEXPON),X(LICENT),X(LITYPE),
     *            X(LOE),X(LOCCNO),X(LVEC),X(LIWRK),X(LSCR),X(LDA),
     *            X(LDB),X(LSPHF),X(LSMHF),NPRIMS,L0,L1,L2,L3)
C
  200 CONTINUE
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE PARLEY  *DECK AIMPAC
      SUBROUTINE AIMPAC(ACO,AIC,EXPON,ICENT,ITYPE,OE,OCCNO,
     *                  VEC,IWRK,SCR,DA,DB,SPHALF,SMHALF,NPRIMS,
     *                  L0,L1,L2,L3)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION METHOD
C
      CHARACTER*4 ATMLAB(106)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,
     *        CANONC,FCORE,FORS,NOCI,EKT,LINSER
C
      DIMENSION ACO(NPRIMS),AIC(NPRIMS),EXPON(NPRIMS),ICENT(NPRIMS),
     *          ITYPE(NPRIMS),OE(L1),OCCNO(L1),VEC(L1,L1),IWRK(L1),
     *          SCR(L1,8),DA(L2),DB(L2),SPHALF(L3),SMHALF(L3)
C
      PARAMETER (MXGTOT=20000, MXSH=5000, MXATM=2000, MXRT=100,
     *           MXNORO=250, MXAO=8192)
C
      COMMON /DESOCI/ OSPIN(MXRT),NEXTSO,MOSETSO,MAXPSO,KSO,NSOCI
      COMMON /DETWFN/ WSTATE(MXRT),SPINS(MXRT),CRIT,PRTTOL,SDET,SZDET,
     *                GRPDET,STSYM,GLIST,
     *                NFLGDM(MXRT),IWTS(MXRT),NCORSV,NCOR,NACTDT,NORBDT,
     *                NADET,NBDET,KDET,KSTDET,IROOT,IPURES,MAXW1,NITDET,
     *                MAXP,NCIDET,IGPDET,KSTSYM,NFTGCI
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /GUGWFN/ NFZC,NMCC,NDOC,NAOS,NBOS,NALP,NVAL,NEXT,NFZV,
     *                IFORS,IEXCIT,ICICI,NOIRR
      COMMON /GVBWFN/ CICOEF(2,12),F(25),ALPHA(325),BETA(325),NO(10),
     *                NCO,NSETO,NOPEN,NPAIR,NORBGVB,NCONF(MXAO),NHAM
      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 /MCINP / METHOD,CISTEP,ACURCY,ENGTOL,DAMP,MICIT,NWORD,
     *                NORBMC,NOROT(2,MXNORO),MOFRZ(15),NPFLG(10),
     *                NOFO,MCFMO,CANONC,FCORE,FORS,NOCI,EKT,LINSER
      COMMON /MP2PAR/ OSPT,TOL,METHODMP,NWDMP2,MPPROP, NACORMP,NBCORMP,
     *                NOAMP,NOBMP,NORBMP,NBFMP,NOMITMP,MOCPHFMP,MAXITC
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, TWO=2.0D+00)
C
      DATA RHF,UHF,ROHF/8HRHF     ,8HUHF     ,8HROHF    /
      DATA GVB,RMC/8HGVB     ,8HMCSCF   /
      DATA GUGA,ALDET,GENCI,FSOCI,ORMAS,CIS
     *         /8HGUGA    ,8HALDET   ,8HGENCI   ,8HFSOCI   ,
     *          8HORMAS   ,8HCIS     /
      DATA RMP2,UMP2/8HMP2      ,8HUMP2     /
      DATA RNONE/8HNONE    /
C
      DATA ATMLAB/'  H ','  HE','  LI','  BE','  B ','  C ',
     *            '  N ','  O ','  F ','  NE','  NA','  MG',
     *            '  AL','  SI','  P ','  S ','  CL','  AR',
     *            '  K ','  CA','  SC','  TI','  V ','  CR',
     *            '  MN','  FE','  CO','  NI','  CU','  ZN',
     *            '  GA','  GE','  AS','  SE','  BR','  KR',
     *            '  RB','  SR','  Y ','  ZR','  NB','  MO',
     *            '  TC','  RU','  RH','  PD','  AG','  CD',
     *            '  IN','  SN','  SB','  TE','  I ','  XE',
     *            '  CS','  BA','  LA','  CE','  PR','  ND',
     *            '  PM','  SM','  EU','  GD','  TB','  DY',
     *            '  HO','  ER','  TM','  YB','  LU','  HF',
     *            '  TA','  W ','  RE','  OS','  IR','  PT',
     *            '  AU','  HG','  TL','  PB','  BI','  PO',
     *            '  AT','  RN','  FR','  RA','  AC','  TH',
     *            '  PA','  U ','  NP','  PU','  AM','  CM',
     *            '  BK','  CF','  ES','  FM','  MD','  NO',
     *            '  LR','    ','    ','    '/
C
C     ----- PUNCH A "WAVEFUNCTION FILE" FOR THE AIMPAC PROGRAM -----
C
      IA = IP
      IF (MASWRK) WRITE(IW,9000) IA
C
C        FILL -AIC- (CONTRACTION COEFFS), -ITYPE-, -ICENT, -EXPON-
C
      L = 0
      DO 130 I = 1,NSHELL
         JST = KSTART(I)
         JEN = JST + KNG(I) - 1
         MINK = KMIN(I)
         MAXK = KMAX(I)
         DO 120 K = MINK,MAXK
            DO 110 J = JST,JEN
               L = L + 1
               ICENT(L) = KATOM(I)
               ITYPE(L) = K
               EXPON(L) = EX(J)
               IF (K.EQ. 1              ) AIC(L) = CS(J)
               IF (K.GE. 2 .AND. K.LE. 4) AIC(L) = CP(J)
               IF (K.GE. 5 .AND. K.LE.10) AIC(L) = CD(J)
               IF (K.GE.11 .AND. K.LE.20) AIC(L) = CF(J)
               IF (K.GE.21 .AND. K.LE.35) AIC(L) = CG(J)
  110       CONTINUE
  120    CONTINUE
  130 CONTINUE
C
C     OBTAIN THE NATURAL ORBITALS OF THE WAVEFUNCTION,
C     AND THE CORRESPONDING OCCUPATION NUMBERS.
C     FOR RHF WE ALSO REQUIRE THE ORBITAL ENERGIES.
C
      CALL VCLR(OE,1,L1)
      CALL VCLR(OCCNO,1,L1)
C
      CALL DERCHK(NDER)
C
      WFNTYP = SCFTYP
      IF(DFTYPE.NE.RNONE) WFNTYP=DFTYPE
C
C            RHF, OR POSSIBLY MP2 OR CIS AFTER RHF
C
      IF(SCFTYP.EQ.RHF) THEN
         IF(CITYP.EQ.CIS  .OR.
     *     (MPLEVL.EQ.2  .AND.
     *      ((NDER.GT.0  .OR.  MPPROP.EQ.1) .AND. NGLEVL.EQ.0))) THEN
            CALL UHFNOS(OCCNO,IWRK,SCR,DA,DB,VEC,SPHALF,SMHALF,
     *                  L0,L1,L2,1)
            WFNTYP = RMP2
            MAXOCC = L0
         ELSE
            CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
            CALL DAREAD(IDAF,IODA,OE ,L1,17,0)
            DO 210 I=1,NA
               OCCNO(I) = TWO
  210       CONTINUE
            MAXOCC = NA
         END IF
      END IF
C
C            UHF, OR POSSIBLY MP2 AFTER UHF
C
      IF(SCFTYP.EQ.UHF) THEN
         CALL UHFNOS(OCCNO,IWRK,SCR,DA,DB,VEC,SPHALF,SMHALF,L0,L1,L2,1)
         IF(MPLEVL.EQ.2  .AND.
     *      ((NDER.GT.0  .OR.  MPPROP.EQ.1) .AND. NGLEVL.EQ.0)) THEN
            WFNTYP = UMP2
            MAXOCC = L0
         ELSE
            MAXOCC = NUM
         END IF
      END IF
C
C            ROHF
C
      IF(SCFTYP.EQ.ROHF) THEN
         CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
         DO 220 I=1,NA
            OCCNO(I) = TWO
            IF(I.GT.NB) OCCNO(I) = ONE
  220    CONTINUE
         MAXOCC = NA
      END IF
C
C            GVB
C
      IF(SCFTYP.EQ.GVB) THEN
         CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
         MAXOCC = NORBGVB
         DO 230 I=1,MAXOCC
            OCCNO(I) = TWO * F(NCONF(I))
  230    CONTINUE
      END IF
C
C            MCSCF, OR CI AFTER ANY TYPE OF SCF
C     NOTE THAT CIS NATURAL ORBITALS ARE TAKEN CARE OF ABOVE, LIKE MP2
C
      IF(SCFTYP.EQ.RMC  .OR.  (CITYP.NE.RNONE .AND. CITYP.NE.CIS)) THEN
         CALL DAREAD(IDAF,IODA,VEC,L3,19,0)
         CALL DAREAD(IDAF,IODA,OCCNO,L1,21,0)
      END IF
      IF(SCFTYP.EQ.RMC) THEN
         IF(CISTEP.EQ.GUGA)  MAXOCC = NMCC+NDOC+NALP+NAOS+NBOS+NVAL+NEXT
         IF(CISTEP.EQ.ALDET) MAXOCC = NCORSV + NACTDT
         IF(CISTEP.EQ.ORMAS) MAXOCC = NCORSV + NACTDT
         IF(CISTEP.EQ.GENCI) MAXOCC = NCORSV + NACTDT
      END IF
      IF(CITYP.NE.RNONE) THEN
         WFNTYP = CITYP
         IF(CITYP.EQ.GUGA)   MAXOCC = NFZC+NDOC+NALP+NAOS+NBOS+NVAL+NEXT
         IF(CITYP.EQ.ALDET)  MAXOCC = NCORSV + NACTDT
         IF(CITYP.EQ.ORMAS)  MAXOCC = NCORSV + NACTDT
         IF(CITYP.EQ.GENCI)  MAXOCC = NCORSV + NACTDT
         IF(CITYP.EQ.FSOCI)  MAXOCC = NCORSV + NACTDT + NEXTSO
         IF(CITYP.EQ.CIS)    MAXOCC = NQMT
      END IF
C
C     ----- WRITE EVERYTHING OUT TO WFNFIL -----
C
      IF (MASWRK) THEN
      WRITE(IA,8000) 'TOP'
      WRITE(IA,8010) (TITLE(I),I=1,10)
C
      WRITE(IA,8020) MAXOCC,NPRIMS,NAT
C
      DO 600 I=1,NAT
         NUCZ = INT(ZAN(I)) + IZCORE(I)
         ZNUC = NUCZ
         WRITE(IA,8030) ATMLAB(NUCZ),I,I,C(1,I),C(2,I),C(3,I),ZNUC
  600 CONTINUE
C
      WRITE(IA,8040) (ICENT(I),I=1,NPRIMS)
      WRITE(IA,8050) (ITYPE(I),I=1,NPRIMS)
      WRITE(IA,8060) (EXPON(I),I=1,NPRIMS)
      END IF
C
C     EXPAND THE NATURAL ORBITALS IN THE AO BASIS
C     INTO THE PRIMITIVE REPRESENTATION
C
      SQRT3 = SQRT(3.0D+00)
      SQRT5 = SQRT(5.0D+00)
      SQRT7 = SQRT(7.0D+00)
      DO 640 II = 1,MAXOCC
         L = 0
         M = 0
         DO 630 I = 1,NSHELL
            KST = KSTART(I)
            KEN = KST + KNG(I) - 1
            JMIN = KMIN(I)
            JMAX = KMAX(I)
            ITYP = KTYPE(I)
            DO 620 J = JMIN,JMAX
               M = M + 1
               DO 610 K = KST,KEN
                  L = L + 1
                  ACO(L) = VEC(M,II) * AIC(L)
                  IF(ITYP.EQ.3) THEN
                     IF(J.GE. 8) ACO(L) = ACO(L) * SQRT3
                  END IF
                  IF(ITYP.EQ.4) THEN
                     IF(J.GE.14) ACO(L) = ACO(L) * SQRT5
                     IF(J.EQ.17) THEN
                        DUMMY    = ACO(L-1)
                        ACO(L-1) = ACO(L)
                        ACO(L)   = DUMMY
                     END IF
                     IF(J.EQ.20) ACO(L) = ACO(L) * SQRT3
                  END IF
                  IF(ITYP.EQ.5) THEN
                     IF(J.GE.24) ACO(L) = ACO(L) * SQRT7
                     IF(J.GE.30) ACO(L) = ACO(L) * SQRT5/SQRT3
                     IF(J.GE.33) ACO(L) = ACO(L) * SQRT3
                  END IF
  610          CONTINUE
  620       CONTINUE
  630    CONTINUE
C
C           MAKE SURE EIGENVALUE DOESN'T PRINT AS STARS.
C
         IF (MASWRK) THEN
            IF(ABS(OE(II)).LT.999.0D+00) THEN
               ORBENG = OE(II)
            ELSE
               ORBENG = SIGN(999.0D+00,OE(II))
            ENDIF
            WRITE(IA,8070) II,OCCNO(II),ORBENG
            WRITE(IA,8080) (ACO(L),L=1,NPRIMS)
         END IF
  640 CONTINUE
C
      VTOT = ENUCR+VEN+VEE
      VIR = ZERO
      IF(EKIN.NE.ZERO) VIR = ABS(VTOT/EKIN)
      IF (MASWRK) THEN
         WRITE(IA,8090)
         WRITE(IA,8100) WFNTYP,ETOT,VIR
         WRITE(IA,8000) 'END'
         CALL FLSHBF(IA)
      END IF
      RETURN
C
 8000 FORMAT(5(1H-),1X,A3,' OF INPUT FILE FOR BADER''S AIMPAC PROGRAM ',
     *       5(1H-))
 8010 FORMAT(10A8)
 8020 FORMAT ('GAUSSIAN',10X,I5,' MOL ORBITALS',1X,I6,' PRIMITIVES',
     *        4X,I5,' NUCLEI')
 8030 FORMAT(A4,I4,4X,'(CENTRE',I3,')',1X,3F12.8,'  CHARGE =',F5.1)
 8040 FORMAT('CENTRE ASSIGNMENTS',2X,20I3)
 8050 FORMAT('TYPE ASSIGNMENTS',4X,20I3)
 8060 FORMAT('EXPONENTS',1X,1P,5E14.7)
 8070 FORMAT('MO',I3,21X,'OCC NO = ',F12.8,' ORB. ENERGY =',F13.8 )
 8080 FORMAT(1P,5E16.8)
 8090 FORMAT('END DATA')
 8100 FORMAT(A8,' ENERGY =',F20.10,'   VIRIAL(-V/T)  =',F13.8)
C
 9000 FORMAT(/1X,'AN AIMPAC INPUT FILE IS BEING WRITTEN TO FILE',I4/)
      END
C*MODULE PARLEY  *DECK MOLMEM
      SUBROUTINE MOLMEM
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      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
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
C        ----- ALLOCATE MEMORY FOR MOLPLT INTERFACE ROUTINE -----
C        NOTE THAT THERE IS A CHARACTER VARIABLE IN -MOLPLT- OF
C        LENGTH MXATM, SO THE MAXIMUM NO. OF ATOM TYPES = MXATM.
C
      MAXBND = 5*NAT + 20*NFRG
      MAXTYP = MXATM
C
      CALL VALFM(LOADFM)
      LIPAIR = LOADFM + 1
      LKINDS = LIPAIR + MAXBND
      LKOLOR = LKINDS + MAXTYP
      LAST   = LKOLOR + MAXTYP
      NEED   = LAST - LOADFM - 1
C
      CALL GETFM(NEED)
      CALL MOLPLT(X(LIPAIR),MAXBND,MAXTYP,X(LKINDS),X(LKOLOR))
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE PARLEY  *DECK MOLPLT
      SUBROUTINE MOLPLT(IPAIRS,MAXBND,MAXTYP,KINDS,KOLOR)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION IPAIRS(2,MAXBND),KINDS(MAXTYP),KOLOR(MAXTYP)
      DIMENSION KOLORS(105)
C
      PARAMETER (MXATM=2000, MXAO=8192, MXFRG=1050, MXPT=2000,
     *           MXFGPT=12000,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
C
      CHARACTER*4 ATMSYM(106),SKINDS(MXATM)
      CHARACTER*8 FRGNME
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /EFMULT/ EFC(3,MXFGPT),EFCHG(2,MXFGPT),EFATRM(MXFGPT),
     *                EFBTRM(MXFGPT),EFATRM2(MXFGPT),EFBRTM2(MXFGPT),
     *                EFDIP(3,MXFGPT),EFQAD(6,MXFGPT),
     *                EFOCT(10,MXFGPT),FRGNME(MXFGPT)
      COMMON /FRGINF/ NMPTS(MXFRG),NMTTPT,IEFC,IEFD,IEFQ,IEFO,
     *                NPPTS(MXFRG),NPTTPT,IEFP,
     *                NRPTS(MXFRG),NRTTPT,IREP,ICHGP,NFRG,
     *                NDPPTS(MXDPPT),NDPTTPT,IEFDP
      COMMON /FRGMSS/ FPMASS(MXPT),FMASS(MXFGPT),
     *                FPNUC(MXPT),FGNUC(MXFGPT)
      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 /RUNLAB/ TITLE(10),ANAM(MXATM),BNAM(MXATM),BFLAB(MXAO)
C
      PARAMETER (ZERO=0.0D+00, TOANGS=0.52917724924D+00)
C
      DATA KOLORS/5,14,          2*11,10,1,4,2,15,14,
     *    2*11,10,6,9,7,3,14,    2*11,10*8,4*10,9,14,
     *    2*11,10*8,4*10,12,14,  2*11,14*13,10*8,5*10,14,
     *    2*11,14*13,3*8/
C
C     ----- PUNCH AN INPUT FILE FOR THE -MOLPLT- PROGRAM -----
C
      IF (MASWRK) THEN
         WRITE(IW,9000)
         WRITE(IP,8000)
      END IF
C
      CALL SETLAB(3,ATMSYM)
      CALL GTBOND(IPAIRS,MAXBND,NBONDS,0)
C
      NKINDS=0
      IF(NAT.EQ.0) GO TO 130
C
      NKINDS=1
      NUCZ = INT(ZAN(1)) + IZCORE(1)
      KINDS(1)  = NUCZ
      SKINDS(1) = ATMSYM(NUCZ)
      KOLOR(1)  = KOLORS(NUCZ)
      DO 120 IAT=2,NAT
         NUCZ = INT(ZAN(IAT)) + IZCORE(IAT)
         MATCH=0
         DO 110 I=1,NKINDS
            IF(NUCZ.EQ.KINDS(I)) MATCH=MATCH+1
  110    CONTINUE
         IF(MATCH.EQ.0) THEN
            NKINDS = NKINDS+1
            KINDS(NKINDS)  = NUCZ
            SKINDS(NKINDS) = ATMSYM(NUCZ)
            KOLOR(NKINDS)  = KOLORS(NUCZ)
         END IF
  120 CONTINUE
C
  130 CONTINUE
      NATMS = NAT
C
      IF(NFRG.GT.0) THEN
         IAT=0
         KMASS=0
         DO 160 IFRG=1,NFRG
            DO 150 III=1,NMPTS(IFRG)
               IAT=IAT+1
               IF(FMASS(III+KMASS).GT.ZERO) THEN
                  NUCZ = INT(FGNUC(III)+0.001D+00)
                  NATMS = NATMS+1
                  MATCH=0
                  DO 140 I=1,NKINDS
                     IF(NUCZ.EQ.KINDS(I)) MATCH=MATCH+1
  140             CONTINUE
                  IF(MATCH.EQ.0) THEN
                     NKINDS = NKINDS+1
                     IF(NKINDS.GT.MAXTYP) THEN
                        IF(MASWRK) WRITE(IW,*)
     *                     'TOO MANY ATOM TYPES IN MOLPLT'
                        CALL ABRT
                     END IF
                     KINDS(NKINDS)  = NUCZ
                     SKINDS(NKINDS) = ATMSYM(NUCZ)
                     KOLOR(NKINDS)  = KOLORS(NUCZ)
                  END IF
               END IF
  150          CONTINUE
               KMASS=KMASS+NMPTS(IFRG)
  160       CONTINUE
         END IF
C
C     ----- OPTIONS CARD -----
C
      IF (MASWRK) WRITE(IP,8010) NATMS,NKINDS,NBONDS,TITLE
C
C     ----- PUNCH ATOMIC SYMBOL, KOLOR, SIZE -----
C     THE BALL SIZES ARE DETERMINED BY PLAYING WITH THE RADIAL
C     EXPECTATION VALUES OF THE HIGHEST AO IN C,SI,GE,SN,PB.
C
      CARBON = 0.3D+00
      DO 210 I=1,NKINDS
         NUCZ = KINDS(I)
                        SIZE=0.75D+00* CARBON
         IF(NUCZ.GT. 2) SIZE=          CARBON
         IF(NUCZ.GT.10) SIZE=1.6D+00 * CARBON
         IF(NUCZ.GT.18) SIZE=1.7D+00 * CARBON
         IF(NUCZ.GT.36) SIZE=1.9D+00 * CARBON
         IF(NUCZ.GT.54) SIZE=2.0D+00 * CARBON
         IF(NUCZ.GT.86) SIZE=2.1D+00 * CARBON
         IF (MASWRK) WRITE(IP,8020) SKINDS(I),KOLOR(I),SIZE
  210 CONTINUE
C
C     ----- PUNCH ATOMIC COORDINATES IN ANGSTROMS -----
C
      IF (MASWRK) THEN
         DO 310 IAT=1,NAT
            NUCZ = INT(ZAN(IAT)) + IZCORE(IAT)
            X = TOANGS * C(1,IAT)
            Y = TOANGS * C(2,IAT)
            Z = TOANGS * C(3,IAT)
            WRITE(IP,8030) ATMSYM(NUCZ),X,Y,Z
  310    CONTINUE
C
         IF(NFRG.GT.0) THEN
            IAT=0
            KMASS=0
            DO 330 IFRG=1,NFRG
               DO 320 III=1,NMPTS(IFRG)
                  IAT=IAT+1
                  IF(FMASS(III+KMASS).GT.ZERO) THEN
                     NUCZ = INT(FGNUC(III)+0.001D+00)
                     X = TOANGS * EFC(1,IAT)
                     Y = TOANGS * EFC(2,IAT)
                     Z = TOANGS * EFC(3,IAT)
                     WRITE(IP,8030) ATMSYM(NUCZ),X,Y,Z
                  END IF
  320          CONTINUE
               KMASS=KMASS+NMPTS(IFRG)
  330       CONTINUE
         END IF
      END IF
C
C     ----- PUNCH BONDED ATOM LIST -----
C
      CALL GTBOND(IPAIRS,MAXBND,NBONDS,IP)
      IF (MASWRK) WRITE(IP,8040)
      RETURN
C
 8000 FORMAT('-------- START OF -MOLPLT- INPUT FILE ----------')
 8010 FORMAT('NATOMS=',I4,'   NKINDS=',I4,'    NBONDS=',I4/10A8)
 8020 FORMAT(A4,I2,F5.2)
 8030 FORMAT(A4,3F12.6)
 8040 FORMAT('-------- END OF -MOLPLT- INPUT FILE ----------')
 9000 FORMAT(1X,'AN INPUT FILE FOR -MOLPLT- HAS BEEN PUNCHED.')
      END
C*MODULE PARLEY  *DECK RPACX
      SUBROUTINE RPACX
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      PARAMETER (MXATM=2000, MXSH=5000)
C
      LOGICAL PKSAVE,PK,PANDK,BLOCK,GOPARR,DSKWRK,MASWRK
C
      COMMON /FMCOM / XX(1)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /PKFIL / PK,PANDK,BLOCK
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
C
C
C     ----- INTERFACE TO THE RPAC PROGRAM SYSTEM -----
C     RPAC IS A NMR AND ELECTRONIC EXCITATION PROPERTIES PROGRAM.
C     CONTACT PROFESSOR TOM BOUMAN AT SOUTHERN ILLINOIS UNIVERSITY
C     FOR MORE INFORMATION ABOUT RPAC.
C     NOTE THAT THE RPAC PROGRAM UNDERSTANDS ONLY RHF WAVEFUNCTIONS.
C
      IF (MASWRK) WRITE(IW,9000)
      L1 = NUM
      L3 = L1*L1
C
C         MAKE SURE INTEGRALS ARE IN THE CORRECT FORMAT,
C         WHICH IS NO SYMMETRY AND NOT IN P OR PK FORMAT.
C
      IF(PK  .OR.  NT.GT.1) THEN
         IF (MASWRK) WRITE(IW,9010)
         IS = 8
         IPK= 9
         NPRINT=7
         NTSAVE = NT
         PKSAVE = PK
         NT     = 1
         PK     = .FALSE.
         CALL JANDK
         NT   = NTSAVE
         PK   = PKSAVE
      END IF
C
C         WRITE BINARY INTERFACE FILE
C
      CALL VALFM(LOADFM)
      LVEC   = LOADFM + 1
      LEIG   = LVEC   + L3
      LAST   = LEIG   + L1
      NEED = LAST-LOADFM
      CALL GETFM(NEED)
      CALL RPAC15(XX(LVEC),XX(LEIG),L1,L3)
      CALL RETFM(NEED)
      RETURN
C
 9000 FORMAT(/1X,'PREPARING BINARY FILES FOR -RPAC- INTERFACE')
 9010 FORMAT(/1X,'RECOMPUTING INTEGRALS IN CORRECT FORMAT FOR RPAC...')
      END
C*MODULE PARLEY  *DECK RPAC15
      SUBROUTINE RPAC15(VEC,EIG,L1,L3)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER*8 SCFPRG,STARS,DAYTIM
      CHARACTER*9 DAY
C
      LOGICAL PACK2E,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000, MXAO=8192,
     *           MXRT=100)
C
      DIMENSION VEC(L3),EIG(L1)
      DIMENSION NS(MXATM),NUC(MXATM)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTFIL/ NINTMX,NHEX,NTUPL,PACK2E,INTTYP,IGRDTYP
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
C
      PARAMETER (ZERO=0.0D+00)
C
C     ----- WRITE INTERFACE FILE FOR RPAC -----
C     THIS BINARY FILE MUST BE SAVED TO BE READ BY RPAC.
C     ROUTINE WRITTEN BY TOM BOUMAN IN MAY 1989.
C
      NFT15 = 15
      CALL SEQOPN(NFT15,'WORK15','UNKNOWN',.FALSE.,'UNFORMATTED')
      CALL SEQREW(NFT15)
C
      NOCC  = NE/2
      NVIRT = NUM-NOCC
C
C   WRITE CARD 1 - HEADER RECORD
C
      STARS  = '**SCF2**'
      SCFPRG = 'GAMESS  '
      DAY    = ' MM/DD/YY'
      DAYTIM = 'HH:MM:SS'
      IF (MASWRK) THEN
      WRITE (NFT15) STARS,TITLE,SCFPRG,DAY,DAYTIM
C
C   WRITE CARD 2 - MOLECULE PARAMETERS
C
      WRITE(NFT15) NAT,ICH,MUL,NUM,NOCC,NVIRT,NA,NB,NE
C
C   WRITE CARD 3 - ATOMIC COORDINATES AND CHARGES
C
      DO 5 I=1,NAT
         NUC(I)=INT(ZAN(I)) + IZCORE(I)
    5 CONTINUE
      WRITE(NFT15) (C(1,J),C(2,J),C(3,J),NUC(J),J=1,NAT)
C
C   WRITE CARD 4 - MO COEFFICIENTS
C
      CALL DAREAD(IDAF,IODA,VEC,L3,15,0)
      WRITE (NFT15) VEC
C
C   WRITE CARD 5 - ORBITAL ENERGIES
C
      CALL DAREAD(IDAF,IODA,EIG(1),L1,17,0)
      WRITE (NFT15) (EIG(I),I=1,L1)
C
C   SETUP FOR BASIS OUTPUT (CARD GROUP 6)
C
      DO  10 I = 1,NAT
         NS(I)=0
   10 CONTINUE
      DO  20 J = 1,NSHELL
         IAT=KATOM(J)
         NS(IAT)=NS(IAT)+1
20    CONTINUE
C
C   FLAGS FOR CARTESIAN (=1) OR PURE L (=0) D AND F FUNCTIONS
C     (THESE ARE HARDWIRED IN GAMESS)
C
      IF6D=1
      IF10F=1
      WRITE (NFT15) IF6D, IF10F
C
      ISHELL=0
      ITYPE=0
      DO 100 IAT=1,NAT
C
C        FOR EACH ATOM, WRITE NO. OF SHELLS
C
         WRITE (NFT15) NS(IAT)
         DO 200 ISH=1,NS(IAT)
            ISHELL=ISHELL+1
            IF (KTYPE(ISHELL).EQ.1) ITYPE=1
            IF (KTYPE(ISHELL).EQ.2) THEN
               IF (KMIN(ISHELL).EQ.2) ITYPE=2
               IF (KMIN(ISHELL).EQ.1) ITYPE=3
            END IF
            IF (KTYPE(ISHELL).GE.3) ITYPE=KTYPE(ISHELL)+1
            NG=KNG(ISHELL)
            K1=KSTART(ISHELL)
            K2=K1+NG-1
C
C..... FOR EACH SHELL, WRITE NO. OF GAUSSIANS, TYPE OF SHELL.
C      ALLOWED SHELL TYPES, CORRESPONDING TO ITYPE=1 THROUGH 5,
C      ARE S, P, SP, D, OR F.
C
            WRITE (NFT15) NG,ITYPE
            DO 300 K=K1,K2
               EXT=EX(K)
               GO TO (310,320,330,340),ITYPE
  310          CST=CS(K)
               CPT=ZERO
               GO TO 390
  320          CST=CP(K)
               CPT=ZERO
               GO TO 390
  330          CST=CS(K)
               CPT=CP(K)
               GO TO 390
  340          CST=CD(K)
               CPT=ZERO
C
C..... FOR EACH GAUSSIAN IN SHELL, WRITE EXPONENT AND COEFFICIENT(S)
C
  390          WRITE (NFT15) EXT,CST,CPT
  300       CONTINUE
  200    CONTINUE
  100 CONTINUE
C
C   CARD 7 - RHF ENERGY
C
      WRITE (NFT15) ETOT
C
C   CARD 8 - ADD LENGTH OF 2-ELECTRON INTEGRAL RECORDS
C     THE JK FILE HAS ONE BUFFER OF FLOATING POINT VALUES,
C     ONE BUFFER OF INTEGER LABELS, AND A SINGLE INTEGER VALUE.
C     -IBFDIM- MEASURES FILE SIZE IN INTEGERS.
C
      IBFDIM = NINTMX*(NWDVAR+1) + 1
      WRITE (NFT15) IBFDIM
C
      CALL SEQREW(NFT15)
      WRITE(IW,9000) NFT15
      END IF
      RETURN
C
 9000 FORMAT(/1X,'RHF RESULTS FOR -RPAC- ARE SAVED ON FILE ',I2/
     *        20H   RECORD   CONTENTS/3X,6(1H-),3X,7(1H-)/
     *        6X,1H1,5X,18HTITLE AND JOB DATA/
     *        6X,1H2,5X,19HMOLECULE PARAMETERS/
     *        6X,1H3,5X,18HATOMIC COORDINATES/
     *        6X,1H4,5X,15HMO COEFFICIENTS/
     *        6X,1H5,5X,16HORBITAL ENERGIES/
     *        6X,1H6,5X,18HGAUSSIAN BASIS SET/
     *        6X,1H7,5X,10HSCF ENERGY/
     *        6X,1H8,5X,31HBUFFER LENGTH FOR INTEGRAL FILE//
     *        1X,'YOU MUST SAVE -WORK15- AND -AOINTS- FOR RPAC.')
      END
C*MODULE PARLEY  *DECK ZEALOT
      SUBROUTINE ZEALOT(CODE)
C
C   READ INTERFACING INFO FROM WORK FILE
C   SORT AND APPEND TO PUNCH FILE
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      CHARACTER CODE*4
      CHARACTER PREFIX*4,LINE*1
C
      DIMENSION LINE(80)
C
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
C
      CALL SEQREW(NFT16)
      WRITE(IW,*) 'PUNCHING INPUT FILE FOR ',CODE
      WRITE(IP,8010) CODE
C
  100 CONTINUE
         READ(NFT16,FMT='(A4,80A1)') PREFIX,LINE
         IF(PREFIX.EQ.'ENDD') GO TO 200
         IF(PREFIX.EQ.CODE) THEN
            CALL ZEALTR(LINE,NCHARS)
            WRITE(IP,FMT='(80A1)') (LINE(I),I=1,NCHARS)
         END IF
      GO TO 100
C
  200 CONTINUE
      WRITE(IP,8020) CODE
      RETURN
C
 8010 FORMAT('-------- START OF -',A4,'- INPUT FILE ----------')
 8020 FORMAT('-------- END OF -',A4,'- INPUT FILE ----------')
      END
C*MODULE PARLEY  *DECK ZEALTR
      SUBROUTINE ZEALTR(LINE,NCHARS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER LINE*1
      DIMENSION LINE(80)
C
C     LOCATE POSITION OF LAST NON-BLANK CHARACTER IN LINE
C     AND RETURN IT IN NCHARS
C
      DO 100 I = 80,1,-1
        IF(LINE(I).NE.' ') THEN
          NCHARS = I
          GO TO 200
        END IF
  100 CONTINUE
C
  200 CONTINUE
      RETURN
      END
C*MODULE PARLEY  *DECK ZEALX
      SUBROUTINE ZEALX
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      CHARACTER*8 LINE08
      LOGICAL AIMPAC,RPAC,PLTORB,MOLPLT
      COMMON /CIFILS/ NFT11,NFT12,NFT13,NFT14,NFT15,NFT16,IDAF20,NEMEMX
      COMMON /INTRFC/ FRIEND,AIMPAC,RPAC,PLTORB,MOLPLT
      DATA IFALLS/4HALL /
C
C     --- FINISH PUNCHING OTHER PROGRAM INPUT FILES ---
C
      WRITE(NFT16,9000)
      WRITE(NFT16,9010)
      WRITE(NFT16,9020)
      WRITE(NFT16,9030)
      WRITE(NFT16,9040)
 9000 FORMAT('GAME#  (1) ENERGY CALCULATION'/'GAMESCFTYPE   RHF'/
     *       'GAMERUNTYPE   SCF'/'GAMELEVEL     2.0'/
     *       'GAMEMAXCYC    60'/'GAMEENTER     1')
 9010 FORMAT('GAME#  (2) PROPERTIES'/'GAMERUNTYPE   ANALYSE'/
     .       'GAMEDMA'/'GAME#  DISTRIBUTED MULTIPOLE ANALYSIS')
 9020 FORMAT('GAME#  (3) GEOMETRY OPTIMIZATION'/
     *       'GAMERUNTYPE   OPTIMIZE'/'GAMEXTOL      0.00005'/
     *       'GAMEVECTORS   1'/'GAMEENTER     1')
 9030 FORMAT('GAME#  (4) VIBRATIONAL FREQUENCIES AT OPTIMIZED GEOMETRY'/
     *       'GAMERUNTYPE   HESSIAN'/'GAMEVECTORS   1'/
     *       'GAMEENTER     1'/'GAMEEOF')
 9040 FORMAT('GAME#   ---   CLEAR DISK SPACE ---'/
     *       'GAMELS -L SCR/XXX.*'/'GAMERM    SCR/XXX.*')
      WRITE(NFT16,FMT='(''ENDD'')')
      CALL FLSHBF(NFT16)
C
      WRITE(UNIT=LINE08,FMT='(A8)') FRIEND
      READ (UNIT=LINE08,FMT='(A4)') IFTEST
      IF (IFTEST.EQ.IFALLS) THEN
         CALL ZEALOT('MELD')
         CALL ZEALOT('HOND')
         CALL ZEALOT('GAUS')
         CALL ZEALOT('GAME')
      ELSE
         CALL ZEALOT(LINE08(1:4))
      END IF
      RETURN
      END
C*MODULE PARLEY  *DECK PLTMEM
      SUBROUTINE PLTMEM
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000, MXAO=8192, MXGSH=30, MXFRG=1050,
     *           MXDFG=5, MXDPPT=MXFRG*MXDFG*12)
      PARAMETER (ZERO=0.0D+00)
C
      LOGICAL   GOPARR,DSKWRK,MASWRK,GRID3D,TOTDEN,XCHECK
C
      DIMENSION BOUND(2,3), NPLT(MXAO), MESH(3)
C
      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
      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 /MACHIN/ NWDVAR,MAXFM,MAXSM,LIMFM,LIMSM
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (NNAM=7)
      DIMENSION QNAM(NNAM), KQNAM(NNAM)
      DATA PLTINP/8HMOLGRF   /
      DATA QNAM/8HGRID3D  ,8HMESH    ,8HBOUND   ,8HNPLOTS  ,8HNPLT    ,
     *          8HCHECK   ,8HTOTAL   /
      DATA KQNAM/0,31,63,1,-1,0,0/
C
      KQNAM(5) = MXAO*10 +1
C
      MAXBND = 5*NAT + 20*NFRG
C
C     $MOLGRF GRID3D=.TRUE.  TOTAL=.TRUE.  CHECK=.TRUE.
C             BOUND(1)=XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
C             MESH(1)=21,21,21   NPLOTS=1  NPLT(1)=10  $END
C
      GRID3D=.FALSE.
      TOTDEN=.FALSE.
      XCHECK=.FALSE.
      MESH(1) = 21
      MESH(2) = 21
      MESH(3) = 21
      BOUND(1,1)=ZERO
      BOUND(2,1)=ZERO
      BOUND(1,2)=ZERO
      BOUND(2,2)=ZERO
      BOUND(1,3)=ZERO
      BOUND(2,3)=ZERO
      NPLOTS=0
      JRET = 0
      CALL NAMEIO(IR,JRET,PLTINP,NNAM,QNAM,KQNAM,
     *            GRID3D,MESH,BOUND,NPLOTS,NPLT,XCHECK,TOTDEN,
     *            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.GT.1) THEN
         IF (MASWRK) WRITE(IW,9001) PLTINP
         CALL ABRT
      END IF
 9001 FORMAT(1X,'**** ERROR IN ',A8,' GROUP...HALTING.')
C
      IF(GRID3D) THEN
      DO 100 J=1,3
         IF(MESH(J).GT.61) THEN
            IF (MASWRK) WRITE(IW,9010) J,MESH(J)
         MESH(J) = 61
         END IF
  100    CONTINUE
 9010 FORMAT(1X,'NUMBER OF GRIDS HAS BEEN SET TO 61.  ',
     *       2I4,' IS TOO LARGE.')
      IF(ABS(BOUND(1,1)-BOUND(2,1)).EQ.ZERO .OR.
     *   ABS(BOUND(1,2)-BOUND(2,2)).EQ.ZERO .OR.
     *   ABS(BOUND(1,3)-BOUND(2,3)).EQ.ZERO ) THEN
         DO 200 J=1,3
            DO 220 IAT=1,NAT
               IF(C(J,IAT).LT.BOUND(1,J)) BOUND(1,J) = C(J,IAT)
               IF(C(J,IAT).GT.BOUND(2,J)) BOUND(2,J) = C(J,IAT)
  220          CONTINUE
            BOUND(1,J) = BOUND(1,J) -3.0D+00
            BOUND(2,J) = BOUND(2,J) +3.0D+00
  200       CONTINUE
         END IF
      IF((NPLOTS.LT.1) .OR. (NPLOTS.GT.10)) THEN
         IF (MASWRK) WRITE(IW,9020) NPLOTS
         NPLOTS  = 1
         NPLT(1) = (NE+1)/2
         END IF
 9020 FORMAT(1X,'NUMBER OF ORBITALS IS NOT REASONABLE (.LT.11).',
     *          '  NPLOTS=',I4)
      DO 300 N=1,NPLOTS
         IF((NPLT(N).LT.1) .OR. (NPLT(N).GT.NUM)) THEN
            IF (MASWRK) WRITE(IW,9030) (NPLT(K),K=1,NPLOTS)
            CALL ABRT
            END IF
  300    CONTINUE
         END IF
 9030 FORMAT(1X,'NOT APPROPRIATE ORBITAL NUMBERS.',
     *      /1X,' NPLT=',10I5)
C
      CALL VALFM(LOADFM)
      CALL GOTFM(NGOT)
      LIPAIR = LOADFM + 1
      LCOEFS = LIPAIR + MAXBND
      LCOEFP = LCOEFS + MXGSH
      LCOEFD = LCOEFP + MXGSH
      LCOEFF = LCOEFD + MXGSH
      LCOEFG = LCOEFF + MXGSH
      LCOEF  = LCOEFG + MXGSH
      LZETA  = LCOEF  + MXGSH*MXGSH
      LTYPE  = LZETA  + MXGSH*MXGSH
      LGAUSS = LTYPE  +(MXGSH*NAT-1)/NWDVAR +1
      LSHELL = LGAUSS +(MXGSH*NAT-1)/NWDVAR +1
      LATOM  = LSHELL +(NAT-1)/NWDVAR +1
      LVEC   = LATOM  + NAT*3
      LNOON  = LVEC   + NUM*NUM
      LGRID  = LNOON  + NUM
      LWORK  = LGRID  +(MESH(1)*MESH(2)*MESH(3)-1) +1
      LAST   = LWORK  +(MESH(1)*MESH(2)*MESH(3)*15-1) +1
      IF(.NOT.GRID3D) THEN
         LAST  = LCOEF
         LZETA = LAST
         LTYPE = LAST
         LGAUSS= LAST
         LSHELL= LAST
         LATOM = LAST
         LVEC  = LAST
         LNOON = LAST
         LGRID = LAST
         LWORK = LAST
         END IF
      NEED   = LAST - LOADFM - 1
C
      IF(NEED.GT.NGOT) THEN
         IF(MASWRK) WRITE(IW,9106) NEED,NGOT
         CALL ABRT
         END IF
      CALL GETFM(NEED)
 9106 FORMAT(/1X,'------------------------------',
     +       /1X,'SORRY, PLEASE INCREASE MEMORY.',
     *       /1X,'------------------------------'/
     *        1X,'NEED,NGOT=',2I20)
C
      CALL PLTORB(X(LIPAIR),MAXBND,X(LCOEFS),X(LCOEFP),X(LCOEFD),
     *            X(LCOEFF),X(LCOEFG),X(LCOEF),X(LZETA),
     *            X(LTYPE),X(LGAUSS),X(LSHELL),
     *            X(LATOM),X(LATOM+NAT),X(LATOM+NAT*2),
     *            X(LVEC),X(LNOON),X(LGRID),X(LWORK),
     *            BOUND,MXGSH,NUM,MESH,NPLOTS,NPLT,
     *            GRID3D,TOTDEN,XCHECK)
      CALL RETFM(NEED)
      RETURN
      END
C*MODULE PARLEY  *DECK PLTORB
      SUBROUTINE PLTORB(IPAIRS,MAXBND,COEFS,COEFP,COEFD,COEFF,COEFG,
     *                  COEF,ZETA,ITYPE,NGAUS,NSHEL,
     *                  XC,YC,ZC,VECS,XNOON,GRID,WORK,
     *                  BOUND,MXGSH,NAOS,MESH,NPLOTS,NPLT,
     *                  GRID3D,TOTDEN,XCHECK)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXSH=5000, MXGTOT=20000, MXATM=2000)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,GRID3D,TOTDEN,XCHECK
C
      DIMENSION IPAIRS(2,MAXBND)
      DIMENSION COEFS(*),COEFP(*),COEFD(*),COEFF(*),COEFG(*)
      DIMENSION COEF(MXGSH,MXGSH),ZETA(MXGSH,*)
      DIMENSION ITYPE(MXGSH,*),NGAUS(MXGSH,*),NSHEL(*)
      DIMENSION XC(*), YC(*), ZC(*), VECS(*), XNOON(*)
      DIMENSION BOUND(*),NPLT(*), MESH(3)
      DIMENSION GRID(*),WORK(*)
C
      CHARACTER*1 TYPE
      CHARACTER*4 ATMSYM(MXATM)
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /EXPONT/ ZS(107),ZP(107),ZD(107),ZF(107),ZG(107)
      COMMON /FMOINF/ NFG,NLAYER,NATFMO,NBDFG,NAOTYP,NBODY
      COMMON /FMORUN/ ESPSCF,IDAFMO,ICURLAY,ICURFG,JCURFG,KCURFG,NCURSH,
     *                NGAU,ICURPOP,IFMOSTP,MONCOR,NAT1E,NEEDR,MODRST,
     *                NORBPROJ
      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 /NSHEL / EX(MXGTOT),CS(MXGTOT),CP(MXGTOT),CD(MXGTOT),
     *                CF(MXGTOT),CG(MXGTOT),CH(MXGTOT),CI(MXGTOT),
     *                KSTART(MXSH),KATOM(MXSH),KTYPE(MXSH),KNG(MXSH),
     *                KLOC(MXSH),KMIN(MXSH),KMAX(MXSH),NSHELL
      COMMON /OUTPUT/ NPRINT,ITOL,ICUT,NORMF,NORMP,NOPK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      PARAMETER (ONE=1.0D+00, HALF=0.5D+00, PT75 = 0.75D+00,
     *           PT187 = 1.875D+00,
     *           PT6562 = 6.5625D+00, TOANGS=0.52917724924D+00)
C
      DATA NONE/4HNONE/
C
C     ----- PUNCH AN INPUT FILE FOR THE -PLTORB- PROGRAM -----
C
      IF (MASWRK) THEN
         WRITE(IW,9000)
         IF(NFG.NE.0) THEN
            WRITE(IP,8005) ICURFG,JCURFG,KCURFG
         ELSE
            WRITE(IP,8000)
         ENDIF
      END IF
      PI   = ACOS(-ONE)
      PI32 = PI * SQRT(PI)
C
C     ----- PUNCH THE OPTION CARDS -----
C
      CALL GTBOND(IPAIRS,MAXBND,NBONDS,0)
      IF(MPCTYP.EQ.NONE) THEN
         IF (MASWRK) WRITE(IP,8010) NAT,NBONDS,NUM
      ELSE
         IF (MASWRK) WRITE(IP,8015) NAT,NBONDS,NUM
      END IF
C
C     ----- PUNCH ATOMIC COORDINATES IN ANGSTROMS -----
C
      CALL SETLAB(2,ATMSYM)
      IF (MASWRK) THEN
      DO 140 IAT=1,NAT
         X = TOANGS * C(1,IAT)
         Y = TOANGS * C(2,IAT)
         Z = TOANGS * C(3,IAT)
         WRITE(IP,8020) ATMSYM(IAT),X,Y,Z
C
         IF(GRID3D) THEN
            XC(IAT) = C(1,IAT)
            YC(IAT) = C(2,IAT)
            ZC(IAT) = C(3,IAT)
            END IF
  140 CONTINUE
C
C     ----- WE CANNOT GUESS WHAT PLANE SHOULD BE PLOTTED -----
C
      WRITE(IP,8030)
      END IF
C
C     ----- PUNCH LIST OF BONDED ATOMS -----
C
      CALL GTBOND(IPAIRS,MAXBND,NBONDS,IP)
C
C     ----- WE CANNOT GUESS WHAT MO-S SHOULD BE PLOTTED -----
C
      IF (MASWRK) WRITE(IP,8060)
C
C     ----- STORE INFORMATION OF BASIS SET FOR GENERATING 3D GRID DATA:
C     IOMO   = TEMP STORAGE OF THE BASIS SET, ON PROGRAM EXIT THE
C              AMPLITUDE GRIDS FOR THE MOS FOR PROGRAM DEND IF
      IF(GRID3D) THEN
         IOMO   = 30
         CALL SEQOPN(  IOMO,'MOGRID','UNKNOWN',.FALSE.,'UNFORMATTED')
         REWIND IOMO
         END IF
C
C     ----- PUNCH THE ATOMIC BASIS SET -----
C
      NSH = 0
      DO 280 IAT=1,NAT
         NUCZ = INT(ZAN(IAT)) + IZCORE(IAT)
C
C             COUNT HOW MANY SHELLS -NSHAT- ARE ON THIS ATOM.
C             PLTORB REGARDS L SHELLS AS SEPARATE S AND P SHELLS
C
         NSHAT = 0
         NLSH  = 0
         DO 220 N=1,NSHELL
            IF(IAT.EQ.KATOM(N)) NSHAT=NSHAT+1
            IF(IAT.EQ.KATOM(N)  .AND.  KMIN(N).EQ.1
     *                          .AND.  KMAX(N).EQ.4) NLSH=NLSH+1
  220    CONTINUE
         NLAYERS = NSHAT+NLSH
         IF (MASWRK) WRITE(IP,8070) ATMSYM(IAT),NLAYERS
         IF(GRID3D) NSHEL(IAT) = NLAYERS
C
         NLSH  = 0
         DO 260 J=1,NSHAT
            NSH = NSH+1
            NLSH=NLSH+1
            NGAUSS = KNG(NSH)
            NCONT  = 1
            IMIN = KMIN(NSH)
            IMAX = KMAX(NSH)
            TYPE = ' '
            IF(IMIN.EQ.1  .AND.  IMAX.EQ. 1) TYPE='S'
            IF(IMIN.EQ.2  .AND.  IMAX.EQ. 4) TYPE='P'
            IF(IMIN.EQ.5  .AND.  IMAX.EQ.10) TYPE='D'
            IF(IMIN.EQ.11 .AND.  IMAX.EQ.20) TYPE='F'
            IF(IMIN.EQ.21 .AND.  IMAX.EQ.35) TYPE='G'
            IF(IMIN.EQ.1  .AND.  IMAX.EQ. 4) TYPE='L'
            IF(TYPE.EQ.' ') THEN
               IF(MASWRK) WRITE(IW,*) 'PLTORB KNOWS S,P,D,F,G SHELLS'
               CALL ABRT
            END IF
C
            NG = KSTART(NSH)-1
C
C               PLTORB NEEDS THE CONTR.COEF. OF NORMALIZED GAUSS.
C
            IF(NORMP .NE. 1) THEN
               DO 240 IG = 1,NGAUSS
                  EE = EX(NG+IG)+EX(NG+IG)
                  FACS = PI32/(EE*SQRT(EE))
                  FACP = HALF*FACS/EE
                  FACD = PT75*FACS/(EE*EE)
                  FACF = PT187*FACS/(EE*EE*EE)
                  FACG = PT6562*FACS/(EE*EE*EE*EE)
                  COEFS(IG) = CS(NG+IG)*SQRT(FACS)
                  COEFP(IG) = CP(NG+IG)*SQRT(FACP)
                  COEFD(IG) = CD(NG+IG)*SQRT(FACD)
                  COEFF(IG) = CF(NG+IG)*SQRT(FACF)
                  COEFG(IG) = CG(NG+IG)*SQRT(FACG)
  240          CONTINUE
            END IF
C
            IF (MASWRK .AND.  MPCTYP.EQ.NONE) THEN
C
            IF(TYPE.NE.'L') THEN
               WRITE(IP,8080) TYPE,NGAUSS,NCONT
               WRITE(IP,8090) (EX(NG+N),N=1,NGAUSS)
               IF(GRID3D) THEN
                  NGAUS(NLSH,IAT) = NGAUSS
                  DO 300 N=1,NGAUSS
                     ZETA(N,NLSH) = EX(NG+N)
  300                CONTINUE
                  END IF
               IF(TYPE.EQ.'S') THEN
                  WRITE(IP,8090) (COEFS(N),N=1,NGAUSS)
                  IF(GRID3D) THEN
                  ITYPE(NLSH,IAT) = 1
                  DO 301 N=1,NGAUSS
                     COEF(N,NLSH) = COEFS(N)
  301                CONTINUE
                  END IF
                  END IF
               IF(TYPE.EQ.'P') THEN
                  WRITE(IP,8090) (COEFP(N),N=1,NGAUSS)
                  IF(GRID3D) THEN
                  ITYPE(NLSH,IAT) = 2
                  DO 302 N=1,NGAUSS
                     COEF(N,NLSH) = COEFP(N)
  302                CONTINUE
                  END IF
                  END IF
               IF(TYPE.EQ.'D') THEN
                  WRITE(IP,8090) (COEFD(N),N=1,NGAUSS)
                  IF(GRID3D) THEN
                  ITYPE(NLSH,IAT) = 3
                  DO 303 N=1,NGAUSS
                     COEF(N,NLSH) = COEFD(N)
  303                CONTINUE
                  END IF
                  END IF
               IF(TYPE.EQ.'F') THEN
                  WRITE(IP,8090) (COEFF(N),N=1,NGAUSS)
                  IF(GRID3D) THEN
                  ITYPE(NLSH,IAT) = 4
                  DO 304 N=1,NGAUSS
                     COEF(N,NLSH) = COEFF(N)
  304                CONTINUE
                  END IF
                  END IF
               IF(TYPE.EQ.'G') THEN
                  WRITE(IP,8090) (COEFG(N),N=1,NGAUSS)
                  IF(GRID3D) THEN
                  ITYPE(NLSH,IAT) = 5
                  DO 305 N=1,NGAUSS
                     COEF(N,NLSH) = COEFG(N)
  305                CONTINUE
                  END IF
                  END IF
            ELSE
               WRITE(IP,8080) 'S',NGAUSS,NCONT
               WRITE(IP,8090) (EX(NG+N),N=1,NGAUSS)
               WRITE(IP,8090) (COEFS(N),N=1,NGAUSS)
               IF(GRID3D) THEN
                  NGAUS(NLSH,IAT) = NGAUSS
                  ITYPE(NLSH,IAT) = 1
                  DO 311 N=1,NGAUSS
                     ZETA(N,NLSH) = EX(NG+N)
                     COEF(N,NLSH) = COEFS(N)
  311                CONTINUE
                  END IF
C
               NLSH = NLSH +1
               WRITE(IP,8080) 'P',NGAUSS,NCONT
               WRITE(IP,8090) (EX(NG+N),N=1,NGAUSS)
               WRITE(IP,8090) (COEFP(N),N=1,NGAUSS)
               IF(GRID3D) THEN
                  NGAUS(NLSH,IAT) = NGAUSS
                  ITYPE(NLSH,IAT) = 2
                  DO 312 N=1,NGAUSS
                     ZETA(N,NLSH) = EX(NG+N)
                     COEF(N,NLSH) = COEFP(N)
  312                CONTINUE
                  END IF
            END IF
C
C     ----- SLATER FUNCTIONS ARE USED -----
C
            ELSE IF (MASWRK .AND.  MPCTYP.NE.NONE) THEN
            IF(TYPE.EQ.'S') THEN
               WRITE(IP,8080) TYPE,1,NCONT
               WRITE(IP,8090) ZS(NUCZ)
               WRITE(IP,8090) ONE
               IF(GRID3D) THEN
                  NGAUS(NLSH,IAT) = 1
                  ITYPE(NLSH,IAT) = 1
                  ZETA(1,NLSH) = ZS(NUCZ)
                  COEF(1,NLSH) = ONE
                  END IF
            ELSE
               WRITE(IP,8080) 'S',1,NCONT
               WRITE(IP,8090)  ZS(NUCZ)
               WRITE(IP,8090)  ONE
               IF(GRID3D) THEN
                  NGAUS(NLSH,IAT) = 1
                  ITYPE(NLSH,IAT) = 1
                  ZETA(1,NLSH) = ZS(NUCZ)
                  COEF(1,NLSH) = ONE
                  END IF
C
               NLSH = NLSH +1
               WRITE(IP,8080) 'P',1,NCONT
               WRITE(IP,8090)  ZP(NUCZ)
               WRITE(IP,8090)  ONE
               IF(GRID3D) THEN
                  NGAUS(NLSH,IAT) = 1
                  ITYPE(NLSH,IAT) = 2
                  ZETA(1,NLSH) = ZP(NUCZ)
                  COEF(1,NLSH) = ONE
                  END IF
            END IF
            END IF
  260    CONTINUE
C
         IF(GRID3D) WRITE(IOMO) ((ZETA(N1,N2),N1=1,MXGSH),N2=1,MXGSH),
     *                          ((COEF(N3,N4),N3=1,MXGSH),N4=1,MXGSH)
  280 CONTINUE
      IF(XCHECK) WRITE(IW,9280) IOMO
 9280 FORMAT(1X,'WRITTEN ZETA AND COEFFICIENTS ONTO ',I2,'-TH FILE.')
C
C     ----- WE CANNOT GUESS WHAT THE ORBITAL TITLE CARDS ARE -----
C
      IF (MASWRK) THEN
         WRITE(IP,8100)
         IF(NFG.NE.0) THEN
            WRITE(IP,8115) ICURFG,JCURFG,KCURFG
         ELSE
            WRITE(IP,8110)
         ENDIF
      END IF
C
C.SK. GENERATE 3D GRID DATA:
      IF(GRID3D)
     *CALL PLT3D(COEF,ZETA,ITYPE,NGAUS,NSHEL,
     *                 XC,YC,ZC,VECS,XNOON,GRID,WORK,ATMSYM,
     *                 BOUND,MXGSH,NAOS,MESH(1),MESH(2),MESH(3),
     *                 NPLOTS,NPLT,IOMO,TOTDEN,XCHECK)
C
      RETURN
C
 8000 FORMAT('------ START OF -PLTORB- INPUT FILE -------')
 8005 FORMAT('------ START OF FMO -PLTORB- ',3I6,' -------')
 8010 FORMAT('NATOMS=',I4,'   NBONDS=',I4,'   NAOS=',I4,
     *       '   NMOS=??   NPLOTS=??'/
     *       'PLANE=????????  KOLOR=1  ANGSTROMS FORMAT (5X,5F15.10)')
 8015 FORMAT('NATOMS=',I4,'   NBONDS=',I4,'   NAOS=',I4,
     *       '   SLATER   NMOS=??   NPLOTS=??'/
     *       'PLANE=????????  KOLOR=1  ANGSTROMS FORMAT (5X,5F15.10)')
 8020 FORMAT(A4,1X,3F20.10)
 8030 FORMAT('?? DESIRED PLOTTING PLANE GOES HERE...'/
     *       '?? PLOTTING PLANE OFFSET GOES HERE...'/
     *       '?? PLOTTING PLANE BOUNDARIES GO HERE...')
 8060 FORMAT('PLOTMOS ??')
 8070 FORMAT(A4,I5)
 8080 FORMAT(A1,2X,2I4)
 8090 FORMAT(1P,5E14.6,:,' >')
 8100 FORMAT('?? ORBITAL TITLE CARDS GO HERE...')
 8110 FORMAT('------ END OF THE -PLTORB- INPUT FILE ------')
 8115 FORMAT('------ END OF FMO -PLTORB- ',3I6,    ' ------')
C
 9000 FORMAT(/1X,'AN INPUT FILE FOR -PLTORB- HAS BEEN PUNCHED.')
      END
C*MODULE PARLEY  *DECK PLT3D
      SUBROUTINE PLT3D(COEF,ZETA,ITYPE,NGAUS,NSHEL,
     *                 XC,YC,ZC,VECS,XNOON,GRID,WORK,ATMSYM,
     *                 BOUND,MXGSH,NAOS,MESH1,MESH2,MESH3,
     *                 NPLOTS,NPLT,IOMO,TOTDEN,CHECK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXATM=2000)
C
      DIMENSION COEF(MXGSH,MXGSH),ZETA(MXGSH,*)
      DIMENSION ITYPE(MXGSH,*),NGAUS(MXGSH,*),NSHEL(*)
      DIMENSION XC(*), YC(*), ZC(*), VECS(NAOS,*), XNOON(*)
      DIMENSION BOUND(2,3),NPLT(*)
      CHARACTER*4  ATMSYM(*)
      DIMENSION GRID(MESH1,MESH2,MESH3),WORK(MESH1,MESH2,MESH3,15)
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 /PRPOPT/ ILOCAL
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      LOGICAL  TOTDEN,CHECK, MCCI
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00)
C
      DATA GUGA     /8HGUGA    /
      DATA AMCSCF /8HMCSCF   /
C
      MCCI = CITYP.EQ.GUGA  .OR.  SCFTYP.EQ.AMCSCF
C
C     THIS PROGRAM DRAWS CONTOUR PLOTS OF MOLECULAR ORBITALS
C     EXPANDED IN A BASIS OF CARTESIAN S,P,D,F OR G GAUSSIANS,
C     OR A BASIS OF S OR P SLATER TYPE ORBITALS.
C
C     THIS PROGRAM DERIVES FROM THE PROGRAM 'CONTOUR5', FROM
C     THE GROUP OF PROFESSOR KLAUS RUEDENBERG AT IOWA STATE
C     UNIVERSITY IN AMES, IOWA.  MARY GILBERT, DAVE FELLER, AND
C     MIKE SCHMIDT OF ISU ARE RESPONSIBLE FOR AN EARLY (1982)
C     VERSION OF THIS PROGRAM, ALSO NAMED 'PLTORB'.
C
C     THIS VERSION IS DUE TO MIKE SCHMIDT, DEPARTMENT OF CHEMISTRY,
C     NORTH DAKOTA STATE UNIVERSITY, FARGO, N.DAK. 58105.
C     THIS VERSION FEATURES (ALMOST) EASY, KEYWORD DRIVEN INPUT,
C     THE ABILITY TO INPUT SHELLS OF BASIS FUNCTIONS, THE ABILITY
C     TO PLOT MO'S EXPANDED IN A 5D BASIS, AND THE READING OF ALL
C     INPUT BEFORE ANY SERIOUS COMPUTATION BEGINS.  MARCH 1987.
C
C           DEFINE FILES USED, OPEN CARD INPUT AND PRINT OUTPUT
C
C     IOPRIM = TEMP STORAGE OF GAUSSIAN PRIMITIVE AMPLITUDE GRIDS
      IOPRIM = 28
      CALL SEQOPN(IOPRIM,'PRGRID','UNKNOWN',.FALSE.,'UNFORMATTED')
C
C     IOCO   = TEMP STORAGE OF CONTRACTED ORBITAL AMPLITUDE GRIDS
      IOCO   = 29
      CALL SEQOPN(  IOCO,'COGRID','UNKNOWN',.FALSE.,'UNFORMATTED')
C
C     IOMO   = TEMP STORAGE OF THE BASIS SET, ON PROGRAM EXIT THE
C              AMPLITUDE GRIDS FOR THE MOS FOR PROGRAM DEND IF
C              (ALREADY OPENED IN PLTORB.)
C
      WRITE(IW,9000)
      IF(TOTDEN) THEN
         WRITE(IW,9040)
         WRITE(IW,9041)
         CALL ABRT
      ELSE
         WRITE(IW,9020)
      END IF
      WRITE(IW,9060) (NPLT(I),I=1,NPLOTS)
C
 9000 FORMAT(/1X,'----------------------------------------',
     *       /1X,'A CARTESIAN GAUSSIAN MO PLOTTING PROGRAM',
     *       /1X,'----------------------------------------')
 9020 FORMAT(/1X,'GENERATE MO GRIDS DATA FOR 3D CG PROGRAM.')
 9040 FORMAT(/1X,'GENERATE TOTAL-DENSITY GRIDS DATA',
     *          ' FOR 3D CG PROGRAM.')
 9041 FORMAT(1X,'UNFORTUNATELY, -MOGRID- HAS INCORRECT PROGRAMMING FOR'/
     *       1X,'THE TOTAL ELECTRON DENSITY, SO RUN IS HALTED.')
 9060 FORMAT( 1X,'ORBITAL NUMBERS ARE ',10I5)
C
C     THE NEXT VAX/VMS SPECIFIC CODE DETERMINES HARD OR SOFT COPY
C     SOFT COPY (TERMINAL) OUTPUT NEEDS TO ERASE THE SCREEN IN
C     BETWEEN PLOTS.  HARD COPY SHOULD OPEN AND CLOSE THE PLOT
C     OUTPUT FILE, SO EACH PLOT IS IN A SEPARATE VERSION NUMBER.
C     GRPHTYPE SHOULD BE 'SOFT' OR 'HARD' ONLY.
C
C     READ THE MO COEFFICIENTS FROM DISK.
      NAOS2 = NAOS*NAOS
      IF(MCCI) THEN
C        POST-HF: 15= CANONICAL ORBITALS.
C                 19= NATURAL ORBITALS AND 21=NOONS.
         CALL DAREAD(IDAF,IODA,VECS,NAOS2,19,0)
         CALL DAREAD(IDAF,IODA,XNOON,NAOS,21,0)
         WRITE(IW,9100)
C
         ELSE
C        HF: 15=ALPHA ORBITALS, 19=BETA ORBITALS.
C            71=ALPHA LMO,      72=BETA LMO.
         IF(ILOCAL.EQ.0) THEN
C           READ CANONICAL ORBITALS:
            CALL DAREAD(IDAF,IODA,VECS,NAOS2,15,0)
            WRITE(IW,9110)
            ELSE
C           READ LOCALIZED ORBITALS:
            CALL DAREAD(IDAF,IODA,VECS,NAOS2,71,0)
            WRITE(IW,9112)
            END IF
         DO 100 I=1,NAOS
            XNOON(I) = ZERO
  100       CONTINUE
         DO 110 IA=1,NA
            XNOON(IA) = ONE
  110       CONTINUE
         DO 120 IB=1,NB
            XNOON(IB) = XNOON(IB) +ONE
  120       CONTINUE
         END IF
      WRITE(IW,9114) (NPLT(I),XNOON(NPLT(I)),I=1,NPLOTS)
C
 9100 FORMAT(/1X,'NATURAL ORBITALS AND NOON ARE READ FROM THE DISK.')
 9110 FORMAT(/1X,'HARTREE-FOCK ORBITALS ARE READ FROM THE DISK.')
 9112 FORMAT(/1X,'LOCALIZED ORBITALS ARE READ FROM THE DISK.')
 9114 FORMAT(3X,5(I7,'=',F7.5),/3X,5(I7,'=',F7.5))
C
C     PRINT OUT THE INPUT DATA
C
      CALL PLTPRT(IW,IOMO,NAT,NAOS,MESH1,MESH2,MESH3,NPLOTS,NPLT,
     *               XC,YC,ZC,ATMSYM,
     *               NSHEL,ITYPE,NGAUS,ZETA,COEF,BOUND,
     *               MXGSH,VECS,XNOON,IP,TOTDEN,CHECK)
C
C     COMPUTE BASIS SET AMPLITUDE GRIDS (3D MESH GRIDS):
      CALL AOGRID(IW,IOPRIM,IOCO,IOMO,MESH1,MESH2,MESH3,NAT,MXGSH,
     *            XC,YC,ZC,BOUND,
     *            NSHEL,NGAUS,ZETA,COEF,ITYPE,GRID,WORK,CHECK,ATMSYM)
      REWIND IOPRIM
C
C     LOOP OVER PLOTTING EACH MOLECULAR ORBITAL
C
      REWIND IOMO
      IPLOT = 1
  410 CONTINUE
         IF(IPLOT.GT.NPLOTS) GO TO 500
C
C        COMPUTE GRID OVER THIS MOLECULAR ORBITAL (3D MESH GRIDS):
         IORB = NPLT(IPLOT)
         CALL MOGRID(IW,IOCO,IOMO,MESH1,MESH2,MESH3,NAOS,IORB,
     *               NPLOTS,NPLT,XNOON,BOUND,
     *               GRID,WORK,VECS,TOTDEN,CHECK,GMIN,GMAX)
C
         IF(TOTDEN) THEN
            WRITE(IW,9402) GMIN,GMAX
            ELSE
            WRITE(IW,9400) IPLOT,IORB,XNOON(IORB),GMIN,GMAX
            END IF
         IPLOT = IPLOT+1
      IF(.NOT.TOTDEN) GO TO 410
 9400 FORMAT(/1X,'PLOT ',I2,' (ORBITAL NO.',I4,' NOON=',F10.5,')',
     *       /1X,'ORBITAL AMPLITUDE MINIMUM,MAXIMUM=',2F12.6)
 9402 FORMAT(/1X,'DENSITY AMPLITUDE MINIMUM,MAXIMUM=',2F12.6)
C
  500 CONTINUE
      REWIND IOCO
      REWIND IOMO
      WRITE(IW,950)
      RETURN
  950 FORMAT(/1X,'EXITING PROGRAM PLTORB NORMALLY')
      END
C*MODULE PARLEY  *DECK PLTPRT
      SUBROUTINE PLTPRT(IOUT,IOMO,NATOMS,NAOS,MESH1,MESH2,MESH3,
     *                  NPLOTS,NPLT,XC,YC,ZC,ATOM,
     *                  NSHEL,ITYPE,NGAUS,ZETA,COEF,BOUND,
     *                  MXSHEL,VECS,XNOON,IP,TOTDEN,CHECK)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL  TOTDEN,CHECK
C
      PARAMETER (MXATM=2000)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
C
      DIMENSION XC(NATOMS),YC(NATOMS),ZC(NATOMS), NPLT(*)
      DIMENSION NSHEL(NATOMS),ITYPE(MXSHEL,*),NGAUS(MXSHEL,*),
     *          ZETA(MXSHEL,*),COEF(MXSHEL,*)
      DIMENSION BOUND(2,3), VECS(NAOS,NAOS), XNOON(*)
      DIMENSION SYMB(5)
      CHARACTER*4  ATOM(*)
C
      PARAMETER (ZERO=0.0D+00, BIG=1.0D+10)
C
      DATA SYMB/8HS       ,8HP       ,8HD       ,8HF       ,8HG       /
C
C     ----- PRINT OUT THE PROGRAM'S INTERPRETATION OF THE INPUT -----
C
      WRITE(IOUT,9000) NATOMS,NAOS,MESH1,MESH2,MESH3,NPLOTS,
     *                 TOTDEN,CHECK
      WRITE(IOUT,9002) (NPLT(I),I=1,NPLOTS)
      IF(TOTDEN) THEN
         WRITE(IP,9004) NATOMS,NAOS,MESH1,MESH2,MESH3,NPLOTS
         ELSE
         WRITE(IP,9005) NATOMS,NAOS,MESH1,MESH2,MESH3,NPLOTS
         END IF
      WRITE(IP  ,9006) (NPLT(I),I=1,NPLOTS)
 9000 FORMAT(/1X,'3D MOLGRF OPTION SETTINGS ARE',
     *       /5X,'NATOMS=',I8,5X,'NAOS  =',I8,
     *       /5X,'GRIDS =',3I8,
     *       /5X,'NPLOTS=',I8,
     *       /5X,'TOTAL =',L8,5X,'CHECK =',L8)
 9002 FORMAT( 5X,'ORBITAL NO.',10I4)
 9004 FORMAT(' $MOLGRF  GRID3D=.TRUE.   TOTAL=.TRUE.  NATOMS=',I3,
     *   /10X,'NAOS=',I4,3X,'MESH(1)=',I2,',',I2,',',I2,3X,'NPLOTS=',I2)
 9005 FORMAT(' $MOLGRF  GRID3D=.TRUE.   NATOMS=',I3,'  NAOS=',I4,
     *      /'          MESH(1)=',I2,',',I2,',',I2,'  NPLOTS=',I2)
 9006 FORMAT('          NPLT(1)=',10(I4,','))
C
      WRITE(IOUT,9010)
 9010 FORMAT(/1X,'THE ATOMIC POSITIONS ARE',
     *       /5X,'ATOM',6X,'X',11X,'Y',11X,'Z (BOHRS)')
      DO 100 K=1,NATOMS
         WRITE(IOUT,9020) ATOM(K),XC(K),YC(K),ZC(K)
  100 CONTINUE
 9020 FORMAT(5X,A4,3F12.7)
C
C     PLOTTING BOUNDARY:
C     THE ORIGIN IS SET AS A GRID POINT.
      DX = (BOUND(2,1)-BOUND(1,1))/(MESH1-1)
      DY = (BOUND(2,2)-BOUND(1,2))/(MESH2-1)
      DZ = (BOUND(2,3)-BOUND(1,3))/(MESH3-1)
      DMIN = BIG
      XMIN = ZERO
      YMIN = ZERO
      ZMIN = ZERO
      MINX = 0
      MINY = 0
      MINZ = 0
      ZZ = BOUND(1,3) -DZ
      DO 120 K=1,MESH3
         ZZ = ZZ + DZ
         YY = BOUND(1,2) -DY
         DO 120 J=1,MESH2
            YY = YY + DY
            XX = BOUND(1,1) -DX
            DO 120 I=1,MESH1
               XX = XX + DX
               DIST = XX*XX + YY*YY + ZZ*ZZ
               IF(DIST.LT.DMIN) THEN
                  DMIN = DIST
                  XMIN = XX
                  YMIN = YY
                  ZMIN = ZZ
                  MINX = I
                  MINY = J
                  MINZ = K
                  END IF
  120          CONTINUE
      DO 122 I=1,2
         BOUND(I,1) = BOUND(I,1) -XMIN
         BOUND(I,2) = BOUND(I,2) -YMIN
         BOUND(I,3) = BOUND(I,3) -ZMIN
  122    CONTINUE
      WRITE(IOUT,9030) BOUND(1,1),BOUND(1,2),BOUND(1,3),
     *                 BOUND(2,1),BOUND(2,2),BOUND(2,3)
      WRITE(IOUT,9032) MINX,MINY,MINZ
 9030 FORMAT(/1X,'PLOTTING AREA:     X         Y         Z',
     *       /1X,'          MIN=(',3F10.3,' )',
     *       /1X,'          MAX=(',3F10.3,' )')
 9032 FORMAT(/1X,'       ORIGIN= ',3(I6,4X))
C
C     PUNCH OUT SOME DATA FOR MOLGRAPH...
      WRITE(IP  ,9008) BOUND(1,1),BOUND(2,1),BOUND(1,2),BOUND(2,2),
     *                 BOUND(1,3),BOUND(2,3)
 9008 FORMAT('          BOUND(1)=',6(F7.3,','),'  $END',
     *      /'----- COORDINATES FOR MOLGRAPH ',39('-'))
      DO 140 K=1,NATOMS
         XNUCZ = ZAN(K) + IZCORE(K)
         WRITE(IP  ,9141) ATOM(K),XNUCZ,XC(K),YC(K),ZC(K)
  140 CONTINUE
      WRITE(IP  ,9142)
 9141 FORMAT(A4,F5.1,3F20.10)
 9142 FORMAT(70('-'))
C
      IF(.NOT.CHECK) RETURN
      WRITE(IOUT,9080)
 9080 FORMAT(/1X,'MOLECULAR BASIS SET (ONLY THE FIRST EIGHT',
     *     ' CONTRACTIONS FIT ON THE PAGE)',
     *       /5X,'ATOM',1X,'NSH',1X,'TYPE',8X,'ZETA',8X,'CONTR.COEFS')
      REWIND IOMO
      DO 230 K=1,NATOMS
         NSH=NSHEL(K)
         READ(IOMO) ((ZETA(N1,N2),N1=1,MXSHEL),N2=1,MXSHEL),
     *              ((COEF(N3,N4),N3=1,MXSHEL),N4=1,MXSHEL)
         DO 220 I=1,NSH
            NG=NGAUS(I,K)
            DO 210 KI=1,NG
               IF(KI.EQ.1) THEN
                  WRITE(IOUT,9090) ATOM(K),I,SYMB(ITYPE(I,K)),
     *                           ZETA(KI,I),COEF(KI,I)
               ELSE
                  WRITE(IOUT,9100) ZETA(KI,I),COEF(KI,I)
               END IF
  210       CONTINUE
  220    CONTINUE
  230 CONTINUE
 9090 FORMAT(6X,A4,I3,2X,A4,1X,F15.7,3X,8F11.7)
 9100 FORMAT(              20X,F15.7,3X,8F11.7)
C
      WRITE(IOUT,9110)
 9110 FORMAT(/1X,'THE MO''S TO BE PLOTTED ARE')
      JMAX=0
      DO 330 N=1,NPLOTS,10
         JMIN=JMAX+1
         JMAX=MIN(JMAX+10,NPLOTS)
         WRITE(IOUT,9120) (NPLT(J),J=JMIN,JMAX)
         WRITE(IOUT,9122) (XNOON(J),J=JMIN,JMAX)
         DO 310 I=1,NAOS
            WRITE(IOUT,9130) I,(VECS(I,ABS(NPLT(J))),J=JMIN,JMAX)
  310    CONTINUE
  330 CONTINUE
      RETURN
 9120 FORMAT(         13X,10(I3,9X))
 9122 FORMAT(1X,'NOON',2X,10F12.7)
 9130 FORMAT(1X,I4,    2X,10F12.7)
      END
C*MODULE PARLEY  *DECK AOGRID
      SUBROUTINE AOGRID(IOUT,IOPRIM,IOCO,IOMO,MESH1,MESH2,MESH3,
     *            NATOMS,MXSHEL, XC,YC,ZC,BOUND,
     *            NSHEL,NGAUS,ZETA,COEF,ITYPE,GRID,WORK,PRTGRD,ATOM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
C
      LOGICAL  PRTGRD,SLATER
C
      COMMON /ECP2  / CLP(400),ZLP(400),NLP(400),KFIRST(MXATM,6),
     *                KLAST(MXATM,6),LMAX(MXATM),LPSKIP(MXATM),
     *                IZCORE(MXATM)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /WFNOPT/ SCFTYP,CITYP,DFTYPE,CCTYP,MPLEVL,MPCTYP
C
      DIMENSION XC(NATOMS),YC(NATOMS),ZC(NATOMS)
      DIMENSION NSHEL(NATOMS),NGAUS(MXSHEL,*),ITYPE(MXSHEL,*)
      DIMENSION ZETA(MXSHEL,*),COEF(MXSHEL,*)
      DIMENSION BOUND(2,3)
      DIMENSION GRID(MESH1,MESH2,MESH3),WORK(MESH1,MESH2,MESH3,15)
C
      DIMENSION IQN(15,3,5),IDEGEN(5),IATPQ(104)
C
      CHARACTER*4 ATOM(NATOMS)
      CHARACTER*1 TYPE(5)
      PARAMETER (ZERO=0.0D+00, THREE=3.0D+00,
     *           FIVE=5.0D+00, SEVEN=7.0D+00)
      DATA IQN/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,
     *         1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     *         0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
     *         0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
     *         2,0,0,1,1,0,0,0,0,0,0,0,0,0,0,
     *         0,2,0,1,0,1,0,0,0,0,0,0,0,0,0,
     *         0,0,2,0,1,1,0,0,0,0,0,0,0,0,0,
     *         3,0,0,2,2,1,0,1,0,1,0,0,0,0,0,
     *         0,3,0,1,0,2,2,0,1,1,0,0,0,0,0,
     *         0,0,3,0,1,0,1,2,2,1,0,0,0,0,0,
     *         4,0,0,3,3,1,0,1,0,2,2,0,2,1,1,
     *         0,4,0,1,0,3,3,0,1,2,0,2,1,2,1,
     *         0,0,4,0,1,0,1,3,3,0,2,2,1,1,2/
      DATA IDEGEN/1,3,6,10,15/
      DATA TYPE/'S','P','D','F','G'/
      DATA IATPQ/2*1,8*2,8*3,18*4,18*5,32*6,18*7/
      DATA NONE/4HNONE/
C
C     ----- COMPUTE AMPLITUDE GRID FOR EACH CONTRACTED AO -----
C
C     THIS ROUTINE USES MOST OF THE COMPUTER TIME.
C     FIRST, A GRID OVER EACH GAUSSIAN IN A PARTICULAR SHELL
C     IS COMPUTED.  THEN, THESE ARE SUMMED WITH THE CONTRACTION
C     COEFFICIENTS TO PRODUCE AN AMPLITUDE GRID FOR EACH
C     CONTRACTED AO IN THE MOLECULE.
C
      SLATER = MPCTYP.NE.NONE
C
      SQRT3 = SQRT(THREE)
      SQRT5 = SQRT(FIVE)
      SQRT7 = SQRT(SEVEN)
C
      DX = (BOUND(2,1)-BOUND(1,1))/(MESH1-1)
      DY = (BOUND(2,2)-BOUND(1,2))/(MESH2-1)
      DZ = (BOUND(2,3)-BOUND(1,3))/(MESH3-1)
C
C        LOOP OVER EACH ATOM, AND EACH SHELL ON EACH ATOM
C
C        UNIT IOMO HOLDS THE BASIS SET, ONE RECORD PER ATOM
C        UNIT IOPRIM WILL HOLD THE GRIDS FOR ALL GAUSSIANS IN A SHELL
C        UNIT IOCO WILL RECIEVE THE GRIDS FOR EACH CONTRACTED ORBITAL
C
      REWIND IOCO
      REWIND IOMO
C
C     LOOP OVER EACH ATOM, AND EACH SHELL ON THAT ATOM
      IPQNUM=0
      DO 720 K=1,NATOMS
C
C        MOPAC...
         IF(SLATER) THEN
            NUMATM = INT(ZAN(K)+0.5D+00) + IZCORE(K)
            IPQNUM=IATPQ(NUMATM)
            END IF
C
         READ(IOMO) ((ZETA(N1,N2),N1=1,MXSHEL),N2=1,MXSHEL),
     *              ((COEF(N3,N4),N3=1,MXSHEL),N4=1,MXSHEL)
         XAT = XC(K)
         YAT = YC(K)
         ZAT = ZC(K)
         NSH = NSHEL(K)
C
         DO 710 I=1,NSH
            ITYP = ITYPE(I,K)
            NG   = NGAUS(I,K)
C
C                 LOOP OVER GAUSSIAN PRIMITIVES IN THIS SHELL
C
            REWIND IOPRIM
            DO 390 KI=1,NG
               ZET = ZETA(KI,I)
C
C                 LOOP OVER THE DEGENERACY OF CARTESIAN GAUSSIAN
C
               IDEG=IDEGEN(ITYP)
               DO 240 LMN=1,IDEG
                  LQN=IQN(LMN,1,ITYP)
                  MQN=IQN(LMN,2,ITYP)
                  NQN=IQN(LMN,3,ITYP)
C
                  DO 230 II=1,MESH1
                     XA = BOUND(1,1) +(II-1)*DX -XAT
                     DO 220 JJ=1,MESH2
                        YA = BOUND(1,2) +(JJ-1)*DY -YAT
                        DO 200 KK=1,MESH3
                           ZA = BOUND(1,3) +(KK-1)*DZ -ZAT
C
      GAUS = PSI(ZET,LQN,MQN,NQN,XA,YA,ZA,SLATER,ITYP,IPQNUM)
      IF(ITYP.EQ.3 .AND. LMN.GE.4)   GAUS=SQRT3*GAUS
      IF(ITYP.EQ.4) THEN
         IF(LMN.GE.4 .AND. LMN.LE.9) GAUS=SQRT5*GAUS
         IF(LMN.EQ.10)               GAUS=SQRT3*SQRT5*GAUS
      END IF
      IF(ITYP.EQ.5) THEN
         IF(LMN.GE.4  .AND. LMN.LE.9)  GAUS=SQRT7*GAUS
         IF(LMN.GE.10 .AND. LMN.LE.12) GAUS=(SQRT5*SQRT7/SQRT3)*GAUS
         IF(LMN.GE.13 .AND. LMN.LE.15) GAUS=SQRT5*SQRT7*GAUS
      END IF
                           WORK(II,JJ,KK,LMN) = GAUS
  200                   CONTINUE
  220                CONTINUE
  230             CONTINUE
  240          CONTINUE
C
C               WRITE PRIMITIVE GRIDS TO DISK, WITH OPTIONAL PRINTOUT
C
                DO 360 LMN=1,IDEG
                 CALL SQWRIT(IOPRIM,WORK(1,1,1,LMN),MESH1*MESH2*MESH3)
                   IF(PRTGRD) THEN
                     WRITE(IOUT,910) LMN,KI,ZET,I,TYPE(ITYP),K,ATOM(K)
                     WRITE(IOUT,920) (JJ,JJ=1,MESH2,3)
                     KK=MESH3/2
                     DO 350 II=1,MESH1,3
                WRITE(IOUT,930) II,(WORK(II,JJ,KK,LMN),JJ=1,MESH2,3)
  350                CONTINUE
                   END IF
  360           CONTINUE
  390       CONTINUE
  910 FORMAT(/1X,'GRID FOR SUBSPECIES',I3,' FOR ',I3,
     *       '(',F11.5,') OF SHELL',I3,'(',A1,')',
     *       ', ON ATOM',I3,'(',A4,')')
  920 FORMAT( 5X,21I6)
  930 FORMAT( 1X,I2,2X,21F6.2)
C
C           PERFORM NORMALIZATION OF THE PRIMITIVES,FOR EXAMPLE
C           THE NORMALIZATION FOR XY,XZ,YZ GAUSSIANS IS SQRT(3)
C           LARGER THAN XX,YY,ZZ, AND IS ACCOUNTED FOR ABOVE.
C
            DO 420 KI=1,NG
               ZET = ZETA(KI,I)
               CNST = CNORM(ZET,ITYP,SLATER,IPQNUM)
               COEF(KI,I) = CNST * COEF(KI,I)
  420       CONTINUE
C
C           COMPUTE GRID FOR EACH CONTRACTED ORBITAL
C
            IDEG = IDEGEN(ITYP)
               REWIND IOPRIM
               DO 530 LMN=1,IDEG
                  DO 525 KK=1,MESH3
                     DO 520 JJ=1,MESH2
                        DO 510 II=1,MESH1
                           WORK(II,JJ,KK,LMN) = ZERO
  510                   CONTINUE
  520                CONTINUE
  525             CONTINUE
  530          CONTINUE
C
               DO 590 KI=1,NG
                  COEFX = COEF(KI,I)
                  DO 580 LMN=1,IDEG
                     IF(COEFX.EQ.ZERO) THEN
                           READ(IOPRIM)
                        ELSE
                           READ(IOPRIM) GRID
                           DO 560 II=1,MESH1
                              DO 550 JJ=1,MESH2
                                 DO 555 KK=1,MESH3
      WORK(II,JJ,KK,LMN)=WORK(II,JJ,KK,LMN) + COEFX * GRID(II,JJ,KK)
  555                            CONTINUE
  550                         CONTINUE
  560                      CONTINUE
                        END IF
  580             CONTINUE
  590          CONTINUE
C
C              PRINT GRID OVER CONTRACTED ATOMIC ORBITAL
C
               DO 650 LMN=1,IDEG
                CALL SQWRIT(IOCO,WORK(1,1,1,LMN),MESH1*MESH2*MESH3)
                  IF(PRTGRD) THEN
                     WRITE(IOUT,940) LMN,I,TYPE(ITYP),K,ATOM(K)
                     WRITE(IOUT,920) (JJ,JJ=1,MESH2,3)
                     KK=MESH3/2
                     DO 640 II=1,MESH1,3
               WRITE(IOUT,930) II,(WORK(II,JJ,KK,LMN),JJ=1,MESH2,3)
  640                CONTINUE
                  END IF
  650          CONTINUE
C
  710    CONTINUE
  720 CONTINUE
  940 FORMAT(/1X,'GRID FOR SUBSPECIES',I3,
     *          ' OF SHELL',I3,'(',A1,') ON ATOM',I3,'(',A4,')')
      RETURN
      END
C*MODULE PARLEY  *DECK MOGRID
      SUBROUTINE MOGRID(IOUT,IOCO,IOMO,MESH1,MESH2,MESH3,NAOS,IORB,
     *               NPLOTS,NPLT,XNOON,BOUND,
     *               GRID,WORK,VECS,TOTDEN,PRTGRD,GMIN,GMAX)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL   TOTDEN,PRTGRD
      DIMENSION GRID(MESH1,MESH2,MESH3),WORK(MESH1,MESH2,MESH3)
      DIMENSION XNOON(*),     NPLT(*)
      DIMENSION VECS(NAOS,*), BOUND(2,3)
      PARAMETER (ZERO=0.0D+00, BIG=1.0D+10)
C
C     ----- COMPUTE MOLECULAR ORBITAL AMPLITUDE GRID -----
C
      REWIND IOCO
C
      DO 120 K=1,MESH3
         DO 110 J=1,MESH2
            DO 100 I=1,MESH1
               GRID(I,J,K) = ZERO
  100          CONTINUE
  110       CONTINUE
  120    CONTINUE
C
C     COMPUTE ORBITAL GRID
C
      MO = ABS(IORB)
      CMIN =  BIG
      CMAX = -BIG
      DO 220 I=1,NAOS
         COEF = VECS(I,MO)
C
C        TOTAL DENSITY:
         IF(TOTDEN) THEN
            COEF = ZERO
C                      THE FOLLOWING LOOP ASSUMES THE MO COEFS
C                      OBEY C-DAGGER * C = I, BUT THIS IS NOT
C                      SO, ACTUALLY C-DAGGER * S * C = I.
            DO 230 IX=1,NPLOTS
               MOX = NPLT(IX)
               COEF = COEF + XNOON(MOX)*VECS(I,MOX)*VECS(I,MOX)
  230          CONTINUE
            END IF
C
         IF(COEF.EQ.ZERO) THEN
               READ(IOCO)
               GO TO 220
            END IF
C
         READ(IOCO) WORK
C
         IF(TOTDEN) THEN
         DO 240 L=1,MESH3
            DO 242 K=1,MESH2
               DO 244 J=1,MESH1
                  GRID(J,K,L) = GRID(J,K,L) + COEF* WORK(J,K,L)**2
  244             CONTINUE
  242          CONTINUE
  240       CONTINUE
         ELSE
         DO 204 L=1,MESH3
            DO 202 K=1,MESH2
               DO 200 J=1,MESH1
                  GRID(J,K,L) = GRID(J,K,L) + COEF*WORK(J,K,L)
  200             CONTINUE
  202          CONTINUE
  204       CONTINUE
         END IF
         CMAX = MAX(CMAX,COEF)
         CMIN = MIN(CMIN,COEF)
  220    CONTINUE
C
C     POSSIBLE ORBITAL PHASE CHANGES
C
      IF(.NOT.TOTDEN) THEN
      IFLIP=1
      IF(-CMIN.GT.CMAX) THEN
            WRITE(IOUT,900) MO
            IFLIP = -IFLIP
         END IF
      IF(IORB.LT.0) THEN
            WRITE(IOUT,910) MO
            IFLIP = -IFLIP
         END IF
      IF(IFLIP.GT.0) GO TO 400
      DO 312 K=1,MESH3
         DO 310 J=1,MESH2
            DO 300 I=1,MESH1
               GRID(I,J,K) = -GRID(I,J,K)
  300       CONTINUE
  310    CONTINUE
  312 CONTINUE
      END IF
  900 FORMAT(1X,'AUTOMATIC ORBITAL PHASE CHANGE PERFORMED FOR',
     *          ' ORBITAL',I5)
  910 FORMAT(1X,'REQUESTED PHASE CHANGE MADE FOR ORBITAL',I5)
C
C     OPTIONAL PRINTOUT OF THE GRID
C
  400 CONTINUE
      IF(PRTGRD) THEN
         IF(TOTDEN) THEN
            WRITE(IOUT,922) (NPLT(I),I=1,NPLOTS)
            WRITE(IOUT,924) (XNOON(NPLT(I)),I=1,NPLOTS)
            END IF
         WRITE(IOUT,920) MO,XNOON(MO)
         WRITE(IOUT,930) (JJ,JJ=1,MESH2,3)
         K = MESH3/2
         DO 420 I=1,MESH1,3
            WRITE(IOUT,940) I,(GRID(I,J,K),J=1,MESH2,3)
  420    CONTINUE
      END IF
  920 FORMAT(/1X,'AMPLITUDE GRID FOR MOLECULAR ORBITAL',I4,
     *       /1X,'          NOON =',F10.5)
  922 FORMAT( 1X,5I10)
  924 FORMAT( 1X,5F10.5)
  930 FORMAT( 5X,21I6)
  940 FORMAT( 1X,I2,2X,21F6.2)
C
C     STORE MOLECULAR ORBITAL GRID ON DISK
C
      IF(TOTDEN) IORB=9999
      WORK(1,1,1) = BOUND(2,1)-BOUND(1,1)
      WORK(2,1,1) = BOUND(2,2)-BOUND(1,2)
      WORK(3,1,1) = BOUND(2,3)-BOUND(1,3)
      WRITE(IOMO) WORK(1,1,1),WORK(2,1,1),WORK(3,1,1),IORB
C
      WORK(1,1,1) = XNOON(MO)
      WRITE(IOMO) MESH1,MESH2,MESH3,WORK(1,1,1)
C
      DX = (BOUND(2,1)-BOUND(1,1))/(MESH1-1)
      DY = (BOUND(2,2)-BOUND(1,2))/(MESH2-1)
      DZ = (BOUND(2,3)-BOUND(1,3))/(MESH3-1)
      HALF= 0.5D+00
      LXG = INT(BOUND(1,1)/DX+HALF)-1
      KXG = LXG + MESH1-1
      LYG = INT(BOUND(1,2)/DY+HALF)-1
      KYG = LYG + MESH2-1
      LZG = INT(BOUND(1,3)/DZ+HALF)-1
      KZG = LZG + MESH3-1
      WRITE(IOMO) LXG,KXG,LYG,KYG,LZG,KZG
C
      DO 500 K=1,MESH3
         DO 500 J=1,MESH2
            DO 500 I=1,MESH1
               WRITE(IOMO) GRID(I,J,K)
  500          CONTINUE
C
C     FIND MINIMUM AND MAXIMUM AMPLITUDES
C
      GMIN = BIG
      GMAX =-BIG
      DO 620 K=1,MESH3
         DO 610 J=1,MESH2
            DO 612 I=1,MESH1
               GMAX = MAX(GMAX,GRID(I,J,K))
               GMIN = MIN(GMIN,GRID(I,J,K))
  612       CONTINUE
  610    CONTINUE
  620 CONTINUE
      RETURN
      END
C*MODULE PARLEY  *DECK PSI
      DOUBLE PRECISION FUNCTION PSI(ZETA,L,M,N,XAT,YAT,ZAT,
     *                              SLATER,IPQNUM,ITYP)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL SLATER
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, SMALL=5.0D-05,
     *           FIFTY=50.0D+00)
C
C     ----- EVALUATE THE SLATER FUNCTION OR GAUSSIAN PRIMITIVE -----
C
      XA = XAT
      YA = YAT
      ZA = ZAT
      IF(ABS(XA).LE.SMALL) XA=ZERO
      IF(ABS(YA).LE.SMALL) YA=ZERO
      IF(ABS(ZA).LE.SMALL) ZA=ZERO
      RHO=SQRT(XA*XA+YA*YA+ZA*ZA)
C
      IF(RHO.LE.SMALL) RHO=ZERO
      ISPD = L + M + N
      IF(RHO.EQ.ZERO .AND. ISPD.EQ.0) GO TO 200
      IF(RHO.EQ.ZERO .AND. ISPD.GT.0) GO TO 100
C
C        GENERAL, OFF NUCLEUS CASE
C
      IF(SLATER) THEN
         CHECK = ZETA*RHO
      ELSE
         CHECK = ZETA*RHO*RHO
      END IF
C
      IF(CHECK.GT.FIFTY) GO TO 100
      A = ONE
      B = ONE
      C = ONE
      IF(L.EQ.1) A = XA
      IF(M.EQ.1) B = YA
      IF(N.EQ.1) C = ZA
      IF(L.EQ.2) A = XA*XA
      IF(M.EQ.2) B = YA*YA
      IF(N.EQ.2) C = ZA*ZA
      IF(L.EQ.3) A = XA*XA*XA
      IF(M.EQ.3) B = YA*YA*YA
      IF(N.EQ.3) C = ZA*ZA*ZA
      IF(L.EQ.4) A = XA*XA*XA*XA
      IF(M.EQ.4) B = YA*YA*YA*YA
      IF(N.EQ.4) C = ZA*ZA*ZA*ZA
      PSI =  A * B * C * EXP(-CHECK)
      IF(SLATER) PSI = PSI * RHO**(IPQNUM-ITYP)
      RETURN
C
C        CUSP OF ANY OTHER ORBITAL
C
  100 CONTINUE
      PSI = ZERO
      RETURN
C
C        CUSP OF AN S ORBITAL
C
  200 CONTINUE
      PSI = ONE
      RETURN
      END
C*MODULE PARLEY  *DECK CNORM
      DOUBLE PRECISION FUNCTION CNORM(ZETA,ITYP,SLATER,IPQNUM)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION LMNQN(5),FACTR(5),CONST(7),FACT2N(7)
      PARAMETER (PITO32= 5.568327996831707D+00, TWO=2.0D+00,
     *           ONE=1.0D+00, SQRT3 = 1.7320508D+00,
     *           PI=3.141592653589793D+00)
      LOGICAL FIRST,SLATER
      SAVE FIRST,CONST
C
      DATA LMNQN/0,1,2,3,4/
      DATA FACTR/1.0D+00,1.0D+00,  3.0D+00, 15.0D+00, 105.0D+00/
      DATA FIRST/.TRUE./
C
      IF(FIRST) THEN
         FIRST=.FALSE.
         FOURPI = (TWO+TWO)*PI
         FACT2N(1)=TWO
         DO 20 N=1,7
            TWON = N+N
            IF(N.GT.1) FACT2N(N) = FACT2N(N-1) * (TWON-ONE) * TWON
            CONST(N) = SQRT(((TWO)**(TWON+ONE))/(FOURPI * FACT2N(N)))
   20    CONTINUE
      END IF
C
C     ITYP=1,2,3,4,5 MEANS S,P,D,F,G FUNCTION
C
      IF(SLATER) THEN
C
C     ----- NORMALIZATION FACTOR FOR A SLATER -----
C        RADIAL NORMALIZATION= SQRT[(2*ZETA)**(2*N+1)] / (2N)!
C            ANGULAR FUNCTION= SQRT(1/4*PI) FOR S, SQRT(3/4*PI) FOR P
C
         CNORM = CONST(IPQNUM) * ZETA**(IPQNUM+0.5D+00)
         IF(ITYP.EQ.2) CNORM = CNORM * SQRT3
C
C     ----- NORMALIZATION FACTOR FOR A GAUSSIAN -----
C     GENERAL NORMALIZATION FORMULA FOR X**L,Y**M,Z**N
C     CARTESIAN GAUSSIAN IS: (FACT2=DOUBLE FACTORIAL)
C        A=ZETA**(L+M+N+1.5)
C        B=2**(2*L+2*M+2*N+1.5)
C        C=FACT2(L+L-1)*FACT2(M+M-1)*FACT2(N+N-1)*PITO32
C
      ELSE
         LMN = LMNQN(ITYP)
         A=ZETA**(LMN+1.5D+00)
         B=TWO**(LMN+LMN+1.5D+00)
         C=FACTR(ITYP)*PITO32
         CNORM=SQRT(A*B/C)
      END IF
      RETURN
      END
