C 11 APR 08 - NM,SAN - NEW MODULE CALCULATING G3(MP2,CCSD(T))
C
C*MODULE G3  *DECK G3MP2X
      SUBROUTINE G3MP2X
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     THIS ROUTINE CALCULATES G3(MP2,CCSD) AND G3(MP2,CCSD)//B3LYP
C     MODEL CHEMISTRY.  WRITTEN BY SEAN NEDD AND NIKITA MATSUNAGA AT
C     LONG ISLAND UNIVERSITY IN 2004-2006.
C
      LOGICAL GOPARR,DSKWRK,MASWRK,MPTEST,SG1
      LOGICAL CNVR,G3MP2B3
      CHARACTER*8  CXBASIS
C
      PARAMETER (MXATM=2000, MXRT=100, MXGRID=10)
C
C     ----- ZPE PARAMETERS FOR G3(MP2,CCSD) AND G3(MP2,CCSD)//B3LYP
C
      PARAMETER (SCLZPE1=0.8929D+00,SCLZPE2=0.96D+00)
C
C     ----- G3(MP2,CCSD) PARAMETERS -----
C           CURTISS, ET AL., CPL 314 (1999) 101-107
C
      PARAMETER (HLCA=9.170D-03,HLCB=4.455D-03,
     *           HLCC=9.155D-03,HLCD=1.947D-03)
C
C     ----- PARAMETERS FOR G3(MP2,CCSD)//B3LYP ------
C
      PARAMETER (BHLCA=9.676D-03,BHLCB=4.907D-03,
     *           BHLCC=9.749D-03,BHLCD=2.134D-03)
C
C     ----- ORIGINAL G3(MP2) PARAMETERS -----
C     NOTABEAN: HLC IN CURTISS ET AL. PAPER (JCP, 110 (1999) 4703)
C               HAS ASSIGNED A AND B FOR ATOMS, AND C AND D FOR
C               MOLECULES.  THIS IS WRONG.
C
C     PARAMETER (HLCC=9.345D-03,HLCD=2.021D-03,
C    *           HLCA=9.279D-03,HLCB=4.471D-03)
C
      PARAMETER (TOKCAL=627.51D+00, BOLTZ=1.380662D-23,
     *           AVOGAD=6.022045D+23,
     *           TEMP=2.9815D+02,TOKJ=4.184D+03)
      PARAMETER (ONE5=1.5D+00,TWO5=2.5D+00)
C
      DIMENSION NCOR(103)
C
      COMMON /CCENGY/ ENRG,EREF,EMP,ECCN,ETOT(6),ECORR(6),
     *                DIAGS(3),AMPMX(5,2),IAMPMX(5,4,2),XO1,XO2,
     *                DIFMAX,DIFFENG,ITER,CNVR
      COMMON /CXTHRM/ CXTHERM(11),CXZPE,METHCX,ICXBAS,CXBASIS
      COMMON /DFGRID/ DFTTHR,DFTGTHR,SWOFF,SW0,BSLRD(137),NDFTFG,
     *                NRAD,NTHE,NPHI,NRAD0,NTHE0,NPHI0,
     *                NANGPT(MXGRID),NANGPT0(MXGRID),SG1
      COMMON /DFTPAR/ DFTTYP(20),EXENA,EXENB,EXENC,IDFT34,NAUXFUN,
     *                                                    NAUXSHL
      COMMON /ENRGMP/ EMP2,EMP3,EMP4
      COMMON /ENRGYS/ ENUCR,EELCT,ETOTL,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 /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      DATA NCOR /2*0,
     *             2*1,                  6*1,
     *             2*5,                  6*5,
     *             2*9,          10*9,   6*14,
     *             2*18,         10*18,  6*23,
     *             2*27,  14*27, 10*34,  6*39,
     *             2*43,  14*43,    50/
      DATA CCSDT,RNONE/8HCCSD(T) ,8HNONE    /
C
      ZPESCL=0.0D+00
C
C     ----- PRINT OUT G3(MP2) BANNER -----
C
      IF(MASWRK) THEN
         WRITE(IW,9000)
         IF(MUL.NE.1) THEN
            WRITE(IW,9005)
            CALL ABRT
         END IF
      END IF
      CALL G3INP
      CALL TIMIT(0)
C
C     ----- OPTIMIZATION AND HESSIAN WITH SCF FIRST -----
C
      ICXBAS=1
      MPLEVL=0
      METHCX=1
      IREST =0
      G3MP2B3=.FALSE.
C
      IF(MASWRK) WRITE(IW,9010)
      CALL BASISS
C
      IF(DFTYPE.NE.RNONE) THEN
         G3MP2B3 = .TRUE.
         CALL DFTINP(RNONE)
      END IF
C
      IF(NAT.EQ.1) THEN
         ICXBAS = 1
         MPLEVL = 2
         METHCX = 2
         MPTEST = .FALSE.
         IREST  = 0
         IF(MASWRK) WRITE(IW,9020)
         CALL MP2INP(MPTEST)
        CALL ENERGX
         ENMP2 = EMP2
        GOTO 5000
      ELSE
         CALL SIGX(.FALSE.)
      END IF
C
      CALL TIMIT(1)
      CALL TIMIT(0)
C
C     ----- DONE WITH DFT, SO RESET -----
C
      IF(DFTYPE.NE.RNONE) THEN
         DFTYPE = RNONE
         NDFTFG = 0
         NAUXSHL= 0
         NAUXFUN= 0
         IREST  = 0
         NEVALS = 0
         CALL VCLR(DFTTYP,1,20)
         METHCX = 2
         GOTO 5000
      END IF
C
C     ----- OPTIMIZE WITH MP2/6-31G(D) -----
C
      ICXBAS = 1
      MPLEVL = 2
      METHCX = 2
      MPTEST = .FALSE.
      IREST  = 0
      IF(MASWRK) WRITE(IW,9020)
      CALL MP2INP(MPTEST)
      CALL SIGX(.FALSE.)
      ENMP2= EMP2
      CALL TIMIT(1)
      CALL TIMIT(0)
C
C     ----- SINGLE POINT ENERGY WITH CCSD(T)/6-31G(D) -----
C
 5000 CONTINUE
      IF(MASWRK) WRITE(IW,9030)
      MPLEVL = 0
      CCTYP  = CCSDT
      IREST  = 0
      NEVALS = 0
      CALL CCINP
      CALL ENERGX
      ECCSD = ESCF + ENRG
      ENMP2 = ESCF + EMP
      CALL TIMIT(0)
C
C     ----- SINGLE POINT ENERGY WITH MP2/G3LARGE ------
C
      MPLEVL = 2
      ICXBAS = 2
      CALL BASISS
      IF(MASWRK) WRITE(IW,9500) CXBASIS
C
C     ----- RESETTING BASIS SET -----
C
      CALL RESETB
C
      IREST  = 0
      NEVALS = 0
      CCTYP  = RNONE
      ICXBAS = 2
C
C     ----- SINGLE POINT ENERGY WITH BASIS G3MP2LARGE
C
      IF(MASWRK) WRITE(IW,9040)
      CALL MP2INP(MPTEST)
      CALL ENERGX
      ENMP2L = EMP2
      CALL TIMIT(1)
      CALL TIMIT(0)
C
C     ----- SUMMING UP CONTRIBUTIONS -----
C
      NCORT = 0
      DO I=1,NAT
         ITHAT = INT(ZAN(I))
         NCORT = NCORT + NCOR(ITHAT)
      END DO
      NAVAL = NA - NCORT
      NBVAL = NB - NCORT
C
C     ----- HIGH LEVEL CORRECTIONS -----
C
      IF(NAT.EQ.1) THEN
         FUDGE = -HLCC*NBVAL - HLCD*(NAVAL-NBVAL)
         ICHARG = ICH
         IZATM = INT(ZAN(1))
         CALL G3SPNORB(IZATM,ICHARG,SOVAL)
      ELSE
         FUDGE = -HLCA*NBVAL - HLCB*(NAVAL-NBVAL)
      END IF
C
      IF(G3MP2B3) THEN
         IF(NAT.EQ.1) THEN
            FUDGE = -BHLCC*NBVAL - BHLCD*(NAVAL-NBVAL)
            ICH   = ICHARG
            IZATM = INT(ZAN(1))
            CALL G3SPNORB(IZATM,ICHARG,SOVAL)
         ELSE
            FUDGE = -BHLCA*NBVAL - BHLCB*(NAVAL-NBVAL)
         END IF
      END IF
C
C     ----- THERMAL ENERGIES -----
C
      IF(NAT.EQ.1) THEN
         ETH = (ONE5)*AVOGAD*BOLTZ*TEMP/TOKJ/TOKCAL
         HTH = TWO5*BOLTZ*TEMP*AVOGAD/TOKJ/TOKCAL
      ELSE
         ETH    = CXTHERM(9)/TOKCAL
         HTH    = CXTHERM(10)/TOKCAL
         GTH    = CXTHERM(11)/TOKCAL
         IF(G3MP2B3) THEN
            ZPESCL = CXZPE*SCLZPE2
         ELSE
            ZPESCL = CXZPE*SCLZPE1
         END IF
      END IF
C
      DELCORR = ECCSD  - ENMP2
      DELBAS  = ENMP2L - ENMP2
C
C-----ZPE INCLUDED IN THERMAL ENERGIES, SO MUST BE REMOVED -----
C
      IF(NAT.EQ.1) THEN
C     ----- FOR ATOMS -----
         E0G3MP2= ENMP2 + DELCORR + DELBAS       + FUDGE + SOVAL
         EG3MP2 = ENMP2 + DELCORR + DELBAS + ETH + FUDGE + SOVAL
         HG3MP2 = ENMP2 + DELCORR + DELBAS + HTH + FUDGE + SOVAL
         GG3MP2 = ENMP2 + DELCORR + DELBAS + GTH + FUDGE + SOVAL
      ELSE
C     ----- FOR MOLECULES -----
         E0G3MP2= ENMP2 + DELCORR + DELBAS               + ZPESCL +FUDGE
         EG3MP2 = ENMP2 + DELCORR + DELBAS + ETH - CXZPE + ZPESCL +FUDGE
         HG3MP2 = ENMP2 + DELCORR + DELBAS + HTH - CXZPE + ZPESCL +FUDGE
         GG3MP2 = ENMP2 + DELCORR + DELBAS + GTH - CXZPE + ZPESCL +FUDGE
      END IF
C
      IF(MASWRK) THEN
         IF(G3MP2B3) THEN
            WRITE(IW,9700)
            WRITE(IW,9605)
            IF(NAT.EQ.1) THEN
               WRITE(IW,9710) BHLCC,BHLCD
            ELSE
               WRITE(IW,9710) BHLCA,BHLCB
            END IF
            WRITE(IW,9715) ENMP2,ECCSD,ENMP2L,DELBAS,ZPESCL,SCLZPE2,
     *                FUDGE,GTH,ETH,HTH,E0G3MP2,EG3MP2,HG3MP2,GG3MP2
         ELSE
            WRITE(IW,9600)
            WRITE(IW,9605)
            IF(NAT.EQ.1) THEN
               WRITE(IW,9610) HLCC,HLCD
            ELSE
               WRITE(IW,9610) HLCA,HLCB
            END IF
            WRITE(IW,9615) ENMP2,ECCSD,ENMP2L,DELBAS,ZPESCL,SCLZPE1,
     *                FUDGE,GTH,ETH,HTH,E0G3MP2,EG3MP2,HG3MP2,GG3MP2
         END IF
      END IF
C
C CALCULATION OF HEATS OF FORMATION
C
C     LET:
      E0G = E0G3MP2
      HG3 = HG3MP2
C
      IF(NAT.GT.1) THEN
        CALL HEATS(E0G,HG3,HFOR0,HFOR298)
        IF(MASWRK) WRITE(IW,9900) HFOR0,HFOR298
      END IF
      CALL TIMIT(1)
      RETURN
C
 9000 FORMAT(/5X,50(1H=),/5X,10(1H-),5X,'G3(MP2) CALCULATIONS',
     *       5X,10(1H-),/5X,50(1H=))
 9005 FORMAT(/5X,'***** ERROR:  G3(MP2) NOT IMPLEMENTED FOR OPEN-',
     *      'SHELLS ******'//,5X,'PLEASE WAIT FOR UCCSD(T)',/
     *      'WILL BE AVAILABLE FROM MICHIGAN',/)
 9010 FORMAT(/5X,50(1H-),/5X,5(1H-),5X,'G3(MP2) HF OPTIMIZATION',
     *      ' AND HESSIAN',/5X,50(1H-))
 9020 FORMAT(/5X,50(1H-),/5X,5(1H-),5X,'G3(MP2) MP2 OPTIMIZATION',
     *      /5X,50(1H-))
 9030 FORMAT(/5X,50(1H-),/5X,5(1H-),5X,'G3(MP2) CCSD(T)/6-31G(D)',
     *      /5X,50(1H-))
 9040 FORMAT(/5X,50(1H-),/5X,5(1H-),5X,'G3(MP2) MP2/G3MP2LARGE',
     *      /5X,50(1H-))
 9500 FORMAT(/1X,30(1H-),/1X,'SWAPPING BASIS SET TO ',A8,/1X,30(1H-))
 9600 FORMAT(/4X,64(1H-),/8X,'DETAILS OF G3(MP2) CALCULATIONS')
 9700 FORMAT(/4X,64(1H-),/8X,'DETAILS OF G3(MP2)//B3LYP CALCULATIONS')
 9605 FORMAT(4X,64(1H-),
     *       /4X,'E(G3(MP2))@0K   = E(CCSD(T))   + DE(BASIS)   + ',
     *       'ZPE(SCALED) + HLC',
     *       /4X,'E(G3(MP2))@298K = E(G3(MP2))@0K + (E(THERMAL) - ZPE)',
     *       /4X,'H(G3(MP2))      = E(G3(MP2))@298K + KT',
     *       /4X,'G(G3(MP2))      = H(G3(MP2)) - TDS',
     *      //8X,'E(CCSD(T))      = CCSD(T)/6-31G(D)//MP2/6-31G(D) ',
     *           'FROZEN CORE',
     *       /8X,'DE(BASIS)       = E(MP2/G3MP2LARGE) - ',
     *           'E(MP2/6-31G(D)) FROZEN CORE',
     *       /8X,'GEOMETRY MP2/6-31G(D) CORE CORRELATION INCLUDED')
 9710 FORMAT(8X,'ZPE AND OTHER THERMAL CALCULATIONS AT',
     *       ' B3LYP/6-31G(D)',
     *       /8X,'HLCS USED ARE:  A = ',F8.6,'     B = ',F8.6,
     *           ' HARTREE')
 9610 FORMAT(8X,'ZPE AND OTHER THERMAL CALCULATIONS AT HF/6-31G(D)',
     *       /8X,'HLCS USED ARE:  A = ',F8.6,'     B = ',F8.6,
     *           ' HARTREE')
C
 9615 FORMAT(/4X,64(1H-),/19X,'SUMMARY OF G3(MP2) CALCULATIONS',
     *   /4X,64(1H-),
     *   /4X,'MP2/6-31G(D)    =',F13.6,3X,'CCSD(T)/6-31G(D) =',F13.6,
     *   /4X,'MP2/G3MP2LARGE  =',F13.6,3X,'BASIS CONTRIBUT  =',F13.6,
     *   /4X,'ZPE(HF/6-31G(D))=',F13.6,3X,'ZPE SCALE FACTOR =',F13.6,
     *   /4X,'HLC             =',F13.6,3X,'FREE ENERGY      =',F13.6,
     *   /4X,'THERMAL ENERGY  =',F13.6,3X,'THERMAL ENTHALPY =',F13.6,
     *   /4X,64(1H-),
     *   /4X,'E(G3(MP2)) @ 0K =',F13.6,3X,'E(G3(MP2)) @298K =',F13.6,
     *   /4X,'H(G3(MP2))      =',F13.6,3X,'G(G3(MP2))       =',F13.6,
     *   /4X,64(1H-))
 9715 FORMAT(/4X,64(1H-),/19X,'SUMMARY OF G3(MP2) CALCULATIONS',
     *   /4X,64(1H-),
     *   /4X,'MP2/6-31G(D)    =',F13.6,3X,'CCSD(T)/6-31G(D) =',F13.6,
     *   /4X,'MP2/G3MP2LARGE  =',F13.6,3X,'BASIS CONTRIBUT. =',F13.6,
     *   /4X,'ZPE B3LYP/6-31GD=',F13.6,3X,'ZPE SCALE FACTOR =',F13.6,
     *   /4X,'HLC             =',F13.6,3X,'FREE ENERGY      =',F13.6,
     *   /4X,'THERMAL ENERGY  =',F13.6,3X,'THERMAL ENTHALPY =',F13.6,
     *   /4X,64(1H-),
     *   /4X,'E(G3(MP2)) @ 0K =',F13.6,3X,'E(G3(MP2)) @298K =',F13.6,
     *   /4X,'H(G3(MP2))      =',F13.6,3X,'G(G3(MP2))       =',F13.6,
     *   /4X,64(1H-))
C9800 FORMAT(/4X,64(1H-),/10X,'RESULTS ARE SUSPECT',
C    * /10X,'THERE IS NO OPEN SHELL CCSD(T) CODED CURRENTLY IN GAMESS',
C    *   /10X,'THE MULTIPLICITY IS NOT A SINGLET !!',/)
C9850 FORMAT(/4X,64(1H-),
C    *   /10X,'HEAT OF FORMATION     (0K): ',F10.2,' KCAL/MOL',
C    *   /10X,'HEAT OF FORMATION   (298K): ',F10.2,' KCAL/MOL',
C    *   /10X,'HEAT OF FORMATION (H/298K): ',F10.2,' KCAL/MOL',
C    *   /4X,64(1H-))
 9900 FORMAT(/4X,64(1H-),
     *   /10X,'HEAT OF FORMATION   (0K): ',F10.2,' KCAL/MOL',
     *   /10X,'HEAT OF FORMATION (298K): ',F10.2,' KCAL/MOL',
     *   /4X,64(1H-),
     *   /10X,'HEATS OF FORMATIONS BASED ON NIST DATABASE FROM ',
     */10X,'COMPUTATIONAL CHEMISTRY COMPARISON AND BENCHMARK DATABASE',
     *   /4X,64(1H-),/)
C
      END
C
C*MODULE G3       *DECK CHGBAS
      SUBROUTINE CHGBAS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
      LOGICAL MIDI
C
      PARAMETER (MXSH=5000,MXGSH=30,MXGTOT=20000,MXAO=8192,MXATM=2000,
     *           MXRT=100,ONE=1.0D+00,
     *           ZERO=0.0D+00, TWO=2.0D+00,
     *           TM10=1.0D-10)
C
      CHARACTER*8  CXBASIS
C
      DOUBLE PRECISION LABEL(27)
C
      DIMENSION INTYP(MXSH)
C
      DIMENSION CSINP(MXGTOT),CPINP(MXGTOT),CDINP(MXGTOT),CFINP(MXGTOT)
      DIMENSION NS(MXATM),KS(MXATM),
     *          SCFAC(4),NBFS(27),MINF(27),MAXF(27),NANGM(27)
C
      COMMON /BASSPH/ QMTTOL,ISPHER
      COMMON /CXTHRM/ CXTHERM(11),CXZPE,METHCX,ICXBAS,CXBASIS
      COMMON /ENRGYS/ ENUCR,EELCT,ETOT,SZ,SZZ,ECORE,ESCF,EERD,E1,E2,
     *                VEN,VEE,EPOT,EKIN,ESTATE(MXRT),STATN,EDFT(2)
      COMMON /GEOMOP/ GEOM
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /INTOPT/ ISCHWZ,IECP,NECP,IEXFLD
      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 /RUNLAB/ TITLE(10),A(MXATM),B(MXATM),BFLAB(MXAO)
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /TRANSF/ XSMAL,YSMAL,ZSMAL,XNEW,YNEW,ZNEW,XP,YP,ZP
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
C
      DATA RHF/8HRHF     /
C
C           THIS ROUTINE IS PARTIALLY SET TO READ M AND N SHELLS,
C           BUT THESE ARE NOT IMPLEMETED IN THE INTEGRAL PACKAGE.
C
      DATA LABEL/8HS       ,8HP       ,8HD       ,8HF       ,
     *           8HG       ,8HH       ,8HI       ,8HL       ,
     *           8H1S      ,8H2S      ,8H2P      ,8H2SP     ,
     *           8H3S      ,8H3P      ,8H3D      ,8H3SP     ,8H3SPD    ,
     *           8H4S      ,8H4P      ,8H4D      ,8H4SP     ,8H4SPD    ,
     *           8H5S      ,8H5P      ,8H5D      ,8H5SP     ,8H5SPD    /
C
      DATA NBFS/ 1, 3, 6,10,15,21,28,    4,
     *           1, 1, 3, 4,             1, 3, 6, 4,10,
     *           1, 3, 6, 4,10,          1, 3, 6, 4,10/
      DATA MINF/ 1, 2, 5,11,21,36,57,    1,
     *           1, 1, 2, 1,             1, 2, 5, 1, 1,
     *           1, 2, 5, 1, 1,          1, 2, 5, 1, 1/
      DATA MAXF/ 1, 4,10,20,35,56,84,    4,
     *           1, 1, 4, 4,             1, 4,10, 4,10,
     *           1, 4,10, 4,10,          1, 4,10, 4,10/
      DATA NANGM/1, 2, 3, 4, 5, 6, 7,    2,
     *           1, 1, 2, 2,             1, 2, 3, 2, 3,
     *           1, 2, 3, 2, 3,          1, 2, 3, 2, 3/
C
      DATA DAF/8HDAF     /
C
C     ----- USE SPHERICAL HARMONICS FOR MP2/G3MP2LARGE -----
C
      ISPHER = 1
C
      IERR1=0
      IERR2=0
C
      IGAUSS = 0
      ITYP = 0
C
C     ----- ZERO OUT BASIS SET ARRAYS -----
C
      DO 100 I = 1,MXGTOT
         EX(I) = ZERO
         CS(I) = ZERO
         CP(I) = ZERO
         CD(I) = ZERO
         CF(I) = ZERO
         CG(I) = ZERO
         CH(I) = ZERO
         CI(I) = ZERO
  100 CONTINUE
C
      NFG = 0
      NSHELL = 0
      NGAUSS = 0
      LOC = 0
      ISH = 0
C
      SCFAC(1) = ZERO
      SCFAC(2) = ZERO
      SCFAC(3) = ZERO
      SCFAC(4) = ZERO
C
C     ------------------------------------------------------
C     ----- LOOP OVER ALL ATOMS AND GET BASIS SET -----
C     ------------------------------------------------------
C
      DO 860 ICNT=1,NAT
      NS(ICNT) = 0
      KS(ICNT) = NSHELL+1
      NUCZ=INT(ZAN(ICNT))
C
      IF(IGAUSS .GT. MXGSH) THEN
         IF (MASWRK) WRITE (IW,9010) MXGSH
         CALL ABRT
      END IF
C
C     --- FOR COMPOSITE METHODS WE NEED ALL BASIS SETS HERE ---
C
C     ----- STO-NG -----
C
      IF(CXBASIS .NE. 'CBASSTO ') GO TO 310
      CALL STONG(NUCZ,IGAUSS,CSINP,CPINP,CDINP,IERR1,IERR2,
     *           INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS)
      IF(IERR1.NE.0) GO TO 1560
      IF(IERR2.NE.0) GO TO 1600
      GO TO 200
C
C     ----- N-21G -----
C
  310 CONTINUE
      IF(CXBASIS .NE. 'CBASN21 ') GO TO 320
      CALL N21G(NUCZ,IGAUSS,CSINP,CPINP,CDINP,SCFAC,IERR1,IERR2,
     *          INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS)
      IF(IERR1.NE.0) GO TO 1560
      IF(IERR2.NE.0) GO TO 1600
      GO TO 200
C
C     ----- N-31G BASIS SET -----
C
  320 CONTINUE
      IF(CXBASIS .NE. 'CBASN31 ') GO TO 330
      CALL N31G(NUCZ,IGAUSS,CSINP,CPINP,CDINP,SCFAC,IERR1,IERR2,
     *          INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS)
      IF(IERR1.NE.0) GO TO 1560
      IF(IERR2.NE.0) GO TO 1600
      GO TO 200
C
C     ----- N-311G BASIS SET -----
C
  330 CONTINUE
      IF(CXBASIS .NE. 'CBASN311') GO TO 340
      CALL N311G(NUCZ,IGAUSS,CSINP,CPINP,IERR1,IERR2,
     *           INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS)
      IF(IERR1.NE.0) GO TO 1560
      IF(IERR2.NE.0) GO TO 1600
      GO TO 200
C
C     ----- HUZINAGA GROUP MINI/MIDI BASIS SETS -----
C
  340 CONTINUE
      IF(CXBASIS.NE.'CBASMIN '  .AND.  CXBASIS.NE.'CBASMID ') GO TO 350
      MIDI=.FALSE.
      IF(CXBASIS.EQ.'CBASMID ') MIDI=.TRUE.
      CALL BASHUZ(NUCZ,MIDI,CSINP,CPINP,CDINP,CFINP,SCFAC,IERR1,IERR2,
     *            INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,
     *            EX,CS,CP,CD,CF,KSTART,KATOM,KTYPE,KNG,KLOC,
     *            KMIN,KMAX,NSHELL,MXGTOT,MXSH)
      IF(IERR1 .NE. 0) GO TO 1560
      IF(IERR2 .NE. 0) GO TO 1600
      GO TO 200
C
C     ----- "DOUBLE ZETA" BASIS SETS -----
C
  350 CONTINUE
      IF(CXBASIS.NE.'CBASDH  '  .AND.  CXBASIS.NE.'CBASDZV '
     *      .AND.CXBASIS.NE.'CBASBC  ') GO TO 360
      CALL DZVBAS(NUCZ,CSINP,CPINP,CDINP,SCFAC,IERR1,IERR2,
     *            INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS)
      IF(IERR1 .NE. 0) GO TO 1560
      IF(IERR2 .NE. 0) GO TO 1600
      GO TO 200
C
C     ----- "TRIPLE ZETA" BASIS SETS -----
C
  360 CONTINUE
      IF(CXBASIS.NE.'CBASMC  ' .AND. CXBASIS.NE.'CBASTZV ') GO TO 365
      CALL TZVBAS(NUCZ,CSINP,CPINP,CDINP,SCFAC,IERR1,IERR2,
     *            INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS,.FALSE.)
      IF(IERR1 .NE. 0) GO TO 1560
      IF(IERR2 .NE. 0) GO TO 1600
      GO TO 200
C
C     ----- G3 LARGE BASIS SET -----
C
  365 CONTINUE
      IF(CXBASIS.NE.'CBASG3L ') GO TO 370
      CALL NG3L(NUCZ,ICNT,CSINP,CPINP,CDINP,CFINP,
     *          IERR1,IERR2,INTYP,NANGM,NBFS,MINF,MAXF,LOC,NGAUSS,NS)
      IF(IERR1.NE.0) GO TO 1560
      IF(IERR2.NE.0) GO TO 1600
      GO TO 200
C
  370 CONTINUE
      IF (MASWRK) WRITE(IW,*) 'CHGBAS: IT IS NOT POSSIBLE TO GET HERE!'
      CALL ABRT
C
C     ----- END OF READING BASIS -----
C
  200 CONTINUE
C
  860 CONTINUE
C
C     ----- FORM TRANSFORMATION TABLES FOR ATOMS AND SHELLS.-----
C
      DO 1050 IAT = 1,NAT
         NS1 = KS(IAT)-1
         NS2 = NS(IAT)
         X = C(1,IAT)
         Y = C(2,IAT)
         Z = C(3,IAT)
         XS = X
         YS = Y
         ZS = Z
         CALL LOCAL(X,Y,Z,XS,YS,ZS)
         XSMAL = XS
         YSMAL = YS
         ZSMAL = ZS
         DO 1040 IT = 1,NT
            NN = 9*(IT-1)
            CALL TRANS(NN)
            CALL ROT
            ICTR = -2**20
            DO 1010 I = 1,NAT
               TEST = (XP-C(1,I))**2+(YP-C(2,I))**2+(ZP-C(3,I))**2
               IF(TEST .GT. TM10) GO TO 1010
               ICTR = I
               GO TO 1020
 1010       CONTINUE
 1020       MAPCTR(IAT,IT) = ICTR
            NS3 = KS(ICTR)-1
            DO 1030 ISH = 1,NS2
               MAPSHL(NS1+ISH,IT) = NS3+ISH
 1030       CONTINUE
 1040    CONTINUE
 1050 CONTINUE
C
C     ----- READ COORDINATES FROM DAF IF REQUESTED -----
C     THESE WILL BE USED IN PLACE OF THOSE JUST READ.
C
      GEOM = DAF
      IF(GEOM.EQ.DAF) CALL DAREAD(IDAF,IODA,C,3*NAT,1,0)
C
C     ----- PRINT ATOMIC COORDINATES AND SHELL DESCRIPTIONS -----
C
      IF (MASWRK) THEN
        IF(NFG.EQ.0) THEN
         WRITE (IW,9050)
         DO 1100 IAT = 1,NAT
            WRITE (IW,9060) A(IAT),B(IAT),ZAN(IAT),
     *                      C(1,IAT),C(2,IAT),C(3,IAT)
 1100    CONTINUE
        ELSE
         WRITE(IW,9055)
         DO 1150 IAT = 1,NAT
            WRITE(IW,9065) A(IAT),B(IAT),ZAN(IAT)
 1150    CONTINUE
        END IF
      END IF
C
C     ----- PRINT OUT THE BASIS FUNCTIONS USED -----
C
      IF (MASWRK) THEN
      IF(NFG.EQ.0) THEN
        WRITE (IW,9080)
      ELSE
        WRITE (IW,9085)
      END IF
      IF(NORMP .NE. 1) WRITE (IW,9090)
      IF(NORMF .NE. 1) WRITE (IW,9100)
      WRITE (IW,9110)
      MLDNDA=0
      DO 1420 IAT = 1,NAT
         DO 1260 IT = 1,NT
            IF(MAPCTR(IAT,IT) .GT. IAT) GO TO 1420
 1260    CONTINUE
         WRITE (IW,9120) A(IAT),B(IAT)
         NS1 = KS(IAT)
         NS2 = NS1+NS(IAT)-1
         MLDNDA=MLDNDA+1
C
         DO 1400 ISH = NS1,NS2
            WRITE (IW,9130)
            I1 = KSTART(ISH)
            I2 = I1+KNG(ISH)-1
            ITYP = INTYP(ISH)
C
            DO 1380 IG = I1,I2
C
C         -ITYP- IS AN INDEX INTO THE TABLE OF -LABEL- VALUES ABOVE.
C         1-7=S,P,D,F,G,H,I AND 8=L ARE THE BASIC TYPES OF MOST SETS.
C         9-27 ARE ARCHAIC UNDOCUMENTED STO SUBSHELLS STRINGS, WHICH
C         CONSIST OF S,P,D,L AND UNIMPLEMENTED -M- SHELLS.
C
               GO TO (1280,1290,1300,1310,
     *                1320,1330,1340,1350,
     *                1280,
     *                1280,1290,     1350,
     *                1280,1290,1300,1350,1360,
     *                1280,1290,1300,1350,1360,
     *                1280,1290,1300,1350,1360), ITYP
C
C -----PRINTOUTS ARE BASED ON COMMON BLOCK ELEMENTS------
C -----RATHER THAN THROUGH ROUTINE CALLS-----
C
 1280          C2=CS(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1290          C2=CP(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1300          C2=CD(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1310          C2=CF(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1320          C2=CG(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1330          C2=CH(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1340          C2=CI(IG)
               WRITE(IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2
               GO TO 1380
C
 1350          C2=CS(IG)
               C4=CP(IG)
               WRITE (IW,9140) ISH,LABEL(ITYP),IG,EX(IG),C2,C4
               GO TO 1380
C
 1360          IF(MASWRK) WRITE(IW,*) 'STO -M- SHELLS DO NOT EXIST'
               CALL ABRT
C
 1380       CONTINUE
 1400    CONTINUE
 1420 CONTINUE
C
      END IF
C
C     KILL JOB IF DIMENSIONS HAVE BEEN EXCEEDED
C
      IF(NAT .LE. MXATM) GO TO 1560
      IF (MASWRK) WRITE (IW,9230) MXATM
      CALL ABRT
C
 1560 IF(NSHELL .LE. MXSH) GO TO 1580
      IF (MASWRK) WRITE (IW,9240) MXSH
      CALL ABRT
C
 1580 IF(LOC .LE. MXAO) GO TO 1600
      IF (MASWRK) WRITE (IW,9250) MXAO
      CALL ABRT
C
 1600 IF(NGAUSS .LE. MXGTOT) GO TO 1620
      IF (MASWRK) WRITE (IW,9260) MXGTOT
      CALL ABRT
C
C     ESTABLISH NUMBER OF ELECTRONS, MULTIPLICITY, ETC.
C     LATER, NQMT MAY BE RESET TO VALUES SMALLER THAN NUM.
C
 1620 CONTINUE
      NUM  = LOC
      NQMT = LOC
      NE = NE-ICH
      NA = (NE+MUL-1)/2
      NB = (NE-MUL+1)/2
      IF (MASWRK) THEN
         WRITE(IW,9270) NSHELL,NUM
         IF(ISPHER.EQ.1) WRITE(IW,9271)
         WRITE(IW,9272) NE,ICH,MUL,NA,NB,NAT
      END IF
C
C     ----- CALCULATE "CORRECT" NRE FIRST -----
C     THIS UTILIZES TRUE NUCLEAR CHARGES, NOT OMITTING
C     ANY PROTONS YET FOR SEMI-EMPIRICAL OR FOR ECP RUNS.
C
      ENUCR  = ENUC(NAT,ZAN,C)
      IF (MASWRK) THEN
         WRITE(IW,9280) ENUCR
         IF(IECP.NE.0) WRITE(IW,9282)
      END IF
         IF(NA+NB .NE. NE) THEN
         IF (MASWRK) WRITE(IW,9290) NE,ICH,MUL
         CALL ABRT
         STOP
      END IF
      IF(SCFTYP.EQ.RHF .AND. MUL.NE.1) THEN
         IF (MASWRK) WRITE(IW,9300)
         CALL ABRT
         STOP
      END IF
C 850 CONTINUE
C
C     WRITE ENERGY, MULTIPLICITY RECORD TO DAF
C
      EELCT= ZERO
      ETOT = ENUCR
      SZ   = (MUL+1)/TWO
      SZZ  = SZ*(SZ+ONE)
      ECORE= ZERO
      ESCF = ZERO
      EERD = ZERO
      E1   = ZERO
      E2   = ZERO
      VEN  = ZERO
      VEE  = ZERO
      EPOT = ZERO
      EKIN = ZERO
      DO 1640 I=1,10
         ESTATE(I)=ZERO
 1640 CONTINUE
      STATN=ZERO
C
      RETURN
C
 9010 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' GAUSSIANS PER SHELL')
 9050 FORMAT(/1X,'ATOM',6X,'ATOMIC',22X,'COORDINATES (BOHR)'/
     *         11X,'CHARGE',9X,'X',19X,'Y',19X,'Z')
 9055 FORMAT(/1X,'ATOM',3X,'ATOMIC CHARGE')
 9060 FORMAT(1X,A8,A2,F5.1,F17.10,2F20.10)
 9065 FORMAT(1X,A8,A2,F5.1)
 9080 FORMAT(/5X,'ATOMIC BASIS SET'/5X,16(1H-))
 9085 FORMAT(/5X,'ATOMIC BASIS SET LIBRARY'/5X,16(1H-))
 9090 FORMAT(1X,'THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN',
     *          ' UNNORMALIZED')
 9100 FORMAT(1X,'THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED',
     *          ' TO UNITY')
 9110 FORMAT(/1X,'SHELL TYPE PRIMITIVE    EXPONENT',
     *       10X,'CONTRACTION COEFFICIENTS')
 9120 FORMAT(/1X,A8,A2)
 9130 FORMAT(1X)
 9140 FORMAT(1X,I3,3X,A2,I4,F20.7,4F20.12)
 9230 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL ATOMS')
 9240 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL SHELLS')
 9250 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL BASIS FUNCTIONS')
 9260 FORMAT(1X,'*** ERROR, NO MORE THAN',I6,' TOTAL GAUSSIANS')
 9270 FORMAT(/1X,'TOTAL NUMBER OF BASIS SET SHELLS             =',I5/
     *        1X,'NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =',I5)
 9271 FORMAT(1X,'NOTE: THIS RUN WILL RESTRICT THE MO VARIATION SPACE',
     *          ' TO SPHERICAL HARMONICS.'/
     *       1X,'THE NUMBER OF ORBITALS KEPT IN THE VARIATIONAL SPACE',
     *          ' WILL BE PRINTED LATER.')
 9272 FORMAT( 1X,'NUMBER OF ELECTRONS                          =',I5/
     *        1X,'CHARGE OF MOLECULE                           =',I5/
     *        1X,'SPIN MULTIPLICITY                            =',I5/
     *        1X,'NUMBER OF OCCUPIED ORBITALS (ALPHA)          =',I5/
     *        1X,'NUMBER OF OCCUPIED ORBITALS (BETA )          =',I5/
     *        1X,'TOTAL NUMBER OF ATOMS                        =',I5)
 9280 FORMAT(1X,'THE NUCLEAR REPULSION ENERGY IS',F20.10)
 9282 FORMAT(1X,'NOTE THIS RUN IS USING CORE POTENTIALS, AND THE',
     *          ' NUMBER OF ELECTRONS,'/
     *       1X,'OCCUPIED ORBITALS, AND NUCLEAR REPULSION',
     *          ' ENERGY WILL BE ADJUSTED BELOW'/
     *       1X,'AFTER REMOVAL OF THE CORE CHARGES.')
 9290 FORMAT(//1X,'*** CHECK YOUR INPUT CHARGE AND MULTIPLICITY ***'/
     *       1X,'THERE ARE',I6,' ELECTRONS, WITH CHARGE ICHARG=',I3/
     *       1X,'BUT YOU SELECTED MULTIPLICITY MULT=',I3)
 9300 FORMAT(//1X,'*** SCFTYP=RHF MUST HAVE MULT=1 ***')
      END
C
C*MODULE G3      *DECK RESETB
      SUBROUTINE RESETB
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK,ABEL
C
      PARAMETER (MXAO=8192, MXATM=2000, MXSH=5000, MXRT=100)
C
      CHARACTER*8  CXBASIS
      CHARACTER*4  CDUMMY
C
      COMMON /BASSPH/ QMATOL,ISPHER
      COMMON /CXTHRM/ CXTHERM(11),CXZPE,METHCX,ICXBAS,CXBASIS
      COMMON /IJPAIR/ IA(MXAO)
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /MASSES/ ZMASS(MXATM)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      COMMON /RELWFN/ RMETHOD,QRQMT,CLIG,CLIG2,QRTOL,IQRORD,MODQR,NESOC,
     *                NRATOM,NUMU,NQMTR,NQRDAF,MORDA,NDARELB
      COMMON /RESTAR/ TIMLIM,IREST,NREC,INTLOC,IST,JST,KST,LST
      COMMON /RUNOPT/ RUNTYP,EXETYP,NEVALS,NGLEVL,NHLEVL
      COMMON /SYMMOL/ GROUP,COMPLEX,IGROUP,NAXIS,ILABMO,ABEL
      COMMON /SYMTRY/ MAPSHL(MXSH,48),MAPCTR(MXATM,48),
     *                T(432),INVT(48),NT
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      COMMON /ZMTALT/ NZMAT2,NZVAR2,NVAR2,NZMTRD,ICOORD
C
      DATA C1/8HC1      /
      DATA RHF/8HRHF     /
      DATA G3MP2/8HG3MP2   /
C
C     ----- BEGINNING OF RUN -----
C     1. READ IN BASIS SET + OPTIONS
C     2. GENERATE INITIAL MO'S
C
C     IREST = 0     NORMAL START + NORMAL RUNNING CONDITION.
C     IREST = 1     2E-INTEGRAL RESTART ( 1E +MO'S SAVED)
C     IREST = 2     SCF RESTART ( 1E + MO'S SAVED; 2E SAVED)
C     IREST = 3     1E-GRADIENT RESTART ( MO'S SAVED; NO GRADIENT SAVED)
C     IREST = 4     2E-GRAD RESTART ( MO'S, 1-E GRAD, ORDERED DM2 SAVED)
C
C     ----- KEEP NECESSARY PARAMETERS AND VARIABLES HERE! -----
C     ----- CLOSING DICTIONARY FILE! -----
C
      NOSYM =0
      GROUP1 =GROUP
      IGROUP1=IGROUP
      NAXIS1 =NAXIS
      NT1    =NT
C9060 FORMAT(1X,A8,A2,F5.1,F17.10,2F20.10)
C
      CLOSE(UNIT=10, STATUS='DELETE')
C
C     ---- NOT A RESTART JOB, THEREFORE IREST = 0
C
      IREST  = 0
      NEVALS = 0
C
C     ----- OPEN THE DIRECT ACCESS FILE -----
C
      CALL OPENDA(0)
C
C     ----- CHANGING THE BASIS TO G3MP2LARGE -----
C
C     ----- KEEPING ORIGINAL SYMMETRY SETUP------
      GROUP  =GROUP1
      IGROUP =IGROUP1
      NAXIS  =NAXIS1
      NT     =NT1
C     ----------------------------------------
C
C     ----- RESETING NE TO ORIGINAL VALUE-----
C
      NE = NE + ICH
      ENUCR = ENUC(NAT,ZAN,C)
      CALL DAWRIT(IDAF,IODA,ENUCR,MXRT+15,2,0)
      CALL DAWRIT(IDAF,IODA,C,3*NAT,1,0)
C
C     ----- INITIALIZING INTEGRALS-------
C
      CALL INTIN
C
C     ----- CHANGING THE BASIS SET-------
C
      CALL CHGBAS
C
C     ----- FOR PRINTING PURPOSE CHANGE LABELING -----
C
      CALL SETLAB(1,CDUMMY)
C
      NUMU=NUM
C
      CALL SPDTR
C
C     ----- FILL ATOMIC MASS TABLE (REQUIRES MOLECULE READ IN) -----
C
      CALL RAMS(ZMASS,0)
C
C     ----- FORCE NO SYMMETRY IF REQUESTED -----
C
      IF(NOSYM.EQ.1 .OR. ICOORD.EQ.4) THEN
         GROUP=C1
         IGROUP=1
         NAXIS=1
         NT=1
      END IF
C
C     ----- KILL IF OPEN SHELL, BUT RHF TYPE -----
C
      IF(NA.NE.NB  .AND.  SCFTYP.EQ.RHF) THEN
         IF(MASWRK) WRITE (IW,9095) NA,NB
         CALL ABRT
         STOP
      END IF
C
      IREST = 0
C     IF(IREST.LT.0) IREST = 0
C
C     ----- SET UP TRIANGULAR INDEX MATRIX -----
C
      DO 900 I = 1,MXAO
        IA(I) = (I*I-I)/2
  900 CONTINUE
C
C     ----- PREPARE FOR ORBITAL SYMMETRY ASSIGNMENTS -----
C     AND RESTORE /SYMSPD/ CARTESIAN SPACE TRANSFORMATIONS, IF NEEDED.
C     THIS IS KEPT HERE, AS SOME RUNS, LIKE PCM, TURN SYMMETRY OFF
C
C     ----- G3MP2LARGE REQUIRES ISPHER=1 -----
C
      IF(RUNTYP.EQ.G3MP2.AND.ICXBAS.EQ.2) THEN
       ISPHER = 1
      END IF
      IF (ICOORD.NE.4) THEN
         CALL SYMORB
C        IF(ISPHER.GT.0) CALL SPDTR
      END IF
C
      IF (MASWRK) WRITE(IW,9090)
      CALL TIMIT(1)
C
      RETURN
C
 9090 FORMAT(/1X,'..... DONE SETTING UP THE RUN .....')
 9095 FORMAT(/1X,'AN ODD NUMBER OF ELECTRONS IS IMPOSSIBLE WITH RHF.'/
     *        1X,'THIS RUN HAS',I4,' ALPHA AND',I4,' BETA ELECTRONS.'/
     *        1X,'PLEASE REVIEW YOUR SCFTYP,MULT,ICHARG INPUT VALUES,'/
     *        1X,'AND CHECK THE NUMBER OF ATOMS GENERATED IN $DATA')
      END
C
C*MODULE G3  *DECK G3SPNORB
      SUBROUTINE G3SPNORB(IZATM,ICHARG,SOVAL)
C
C    ADDITION OF SPIN-ORBIT CORRECTIONS TO ATOM CALCULATIONS
C    OF G3(MP2).  REFERENCE: JCP, 110 (1999) 4703 TABLE 1
C    ORIGINAL UNITS MAU CHANGED TO AU
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00)
C
      IF(IZATM.EQ.1)     SOVAL=  ZERO
      IF(IZATM.EQ.2) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.3) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
         IF(ICHARG.EQ.-1) SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.4) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.5) THEN
         IF(ICHARG.EQ.0)  SOVAL= -5.0D-05
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
         IF(ICHARG.EQ.-1) SOVAL= -3.0D-05
      ENDIF
      IF(IZATM.EQ.6) THEN
         IF(ICHARG.EQ.0)  SOVAL= -1.4D-04
         IF(ICHARG.EQ.1)  SOVAL= -2.0D-04
         IF(ICHARG.EQ.-1) SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.7) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL= -4.3D-04
      ENDIF
      IF(IZATM.EQ.8) THEN
         IF(ICHARG.EQ.0)  SOVAL= -3.6D-04
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
         IF(ICHARG.EQ.-1) SOVAL= -2.6D-04
      ENDIF
      IF(IZATM.EQ.9) THEN
         IF(ICHARG.EQ.0)  SOVAL= -6.1D-04
         IF(ICHARG.EQ.1)  SOVAL= -6.7D-04
         IF(ICHARG.EQ.-1) SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.10) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL= -1.19D-03
      ENDIF
      IF(IZATM.EQ.11) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
         IF(ICHARG.EQ.-1) SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.12) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.13) THEN
         IF(ICHARG.EQ.0)  SOVAL= -3.4D-04
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
         IF(ICHARG.EQ.-1) SOVAL= -2.8D-04
      ENDIF
      IF(IZATM.EQ.14) THEN
         IF(ICHARG.EQ.0)  SOVAL= -6.8D-04
         IF(ICHARG.EQ.1)  SOVAL= -9.3D-04
         IF(ICHARG.EQ.-1) SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.15) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL= -1.43D-03
         IF(ICHARG.EQ.-1) SOVAL= -4.5D-04
      ENDIF
      IF(IZATM.EQ.16) THEN
         IF(ICHARG.EQ.0)  SOVAL= -8.9D-04
         IF(ICHARG.EQ.1)  SOVAL=  ZERO
         IF(ICHARG.EQ.-1) SOVAL= -8.8D-04
      ENDIF
      IF(IZATM.EQ.17) THEN
         IF(ICHARG.EQ.0)  SOVAL= -1.34D-03
         IF(ICHARG.EQ.1)  SOVAL= -1.68D-03
         IF(ICHARG.EQ.-1) SOVAL=  ZERO
      ENDIF
      IF(IZATM.EQ.18) THEN
         IF(ICHARG.EQ.0)  SOVAL=  ZERO
         IF(ICHARG.EQ.1)  SOVAL= -2.18D-03
      ENDIF
      IF(IZATM.GE.19) THEN
         SOVAL = ZERO
         IF(MASWRK) WRITE(IW,9010) IZATM
      END IF
      RETURN
 9010 FORMAT(1X,'WE DO NOT HAVE STORED SPIN-ORBIT COUPLING FOR ATOM',I5/
     *       1X,'ASSIGNING SPIN-ORBIT COUPLING TO BE ZERO.')
      END
C
C*MODULE G3  *DECK G3INP
      SUBROUTINE G3INP
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
C     ----- READ IN $G3MP2 GROUP -----
C        SEVERAL SCENERIOS FOR GEOMETRY, ORBITALS, GRADIENT, AND HESSIAN
C     FOR G3(MP2,CCSD)
C     1) HF/6-31G(D) ORBITALS AT FIRST GEOMETRY KNOWN
C     2) HF/6-31G(D) GRADIENT OF FIRST GEOMETRY KNOWN
C     3) HF/6-31G(D) HESSIAN KNOWN
C     4) MP2/6-31G(D) SECOND GEOMETRIES KNOWN
C     5) HF/6-31G(D) ORBITALS AT MP2/6-31G(D) SECOND GEOM KNOWN
C     6) RHF/G3MP2LARGE ORBITALS KNOWN
C
C     FOR G3(MP2,CCSD)//B3
C     1) B3LYP/6-31G(D) ORBITALS KNOWN
C     2) B3LYP/6-31G(D) HESSIAN KNOWN
C     3) HF/6-31G(D) ORBITALS AT B3LYP/6-31G(D) GEOM KNOWN
C     4) HF/G3MP2LARGE ORBITALS KNOWN
C
      COMMON /WFNOPT/ SCFTYP,VBTYP,DFTYPE,TDDFTYP,CITYP,CCTYP,
     *                MPLEVL,MPCTYP
      DATA RNONE /8HNONE    /
      IF(DFTYPE.EQ.RNONE) THEN
      ENDIF
      RETURN
      END
C
C*MODULE G3 *DECK HEATS
      SUBROUTINE HEATS(E0G,HG3,HFOR0,HFOR298)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
C
C     DIMENSION G3H(MXATM),HG3H(MXATM),
      DIMENSION SO(MXATM),ATMFMN(MXATM),HATFMN(MXATM)
C
C     IF (HF.EQ.'FORM') THEN
         CALL DLHFMN(E0G,HG3,SO,ATMFMN,HATFMN,
     *               HFOR0,HFOR298)
C     ELSE
C        CALL DLHRXN( )
C     END IF
C
      RETURN
      END
C
C*MODULE G3 *DECK DLHFMN
      SUBROUTINE DLHFMN(E0G,HG3,SO,ATMFMN,HATFMN,
     *                  HFOR0,HFOR298)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
      DIMENSION ATMFMN(MXATM),SO(MXATM),HATFMN(MXATM)
      DIMENSION G3H(MXATM),HG3H(MXATM)
C
      CALL RDINPF(HG3H,G3H)
      CALL FILATM(ATMFMN,HATFMN,SO)
      CALL GETH(E0G,HG3,HG3H,G3H,ATMFMN,HATFMN,
     *          SO,HFOR0,HFOR298)
      RETURN
      END
C
C*MODULE G3 *DECK RDINPF
      SUBROUTINE RDINPF(HG3H,G3H)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MXATM=2000)
      PARAMETER (TOKCAL=627.51D+00, BOLTZ=1.380662D-23,
     *           AVOGAD=6.022045D+23,
     *           TEMP=2.9815D+02,TOKJ=4.184D+03)
      PARAMETER (TWO5=2.5D+00)
C
      DIMENSION HG3H(MXATM),G3H(MXATM)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      HTH = TWO5*BOLTZ*TEMP*AVOGAD/TOKJ/TOKCAL
C
      DO J=1,NAT
        NAZ = INT(ZAN(J))
        IF(NAZ.EQ.1) THEN
          G3H(J) = -0.50177D+00
        ELSEIF(NAZ.EQ.2) THEN
          G3H(J) = -2.90235D+00
        ELSEIF(NAZ.EQ.3) THEN
          G3H(J) = -7.43397D+00
        ELSEIF(NAZ.EQ.4) THEN
          G3H(J) = -14.62907D+00
        ELSEIF(NAZ.EQ.5) THEN
          G3H(J) = -24.60679D+00
        ELSEIF(NAZ.EQ.6) THEN
          G3H(J) = -37.78899D+00
        ELSEIF(NAZ.EQ.7) THEN
          G3H(J) = -54.52478D+00
        ELSEIF(NAZ.EQ.8) THEN
          G3H(J) = -74.98920D+00
        ELSEIF(NAZ.EQ.9) THEN
          G3H(J) = -99.64020D+00
        ELSEIF(NAZ.EQ.10) THEN
          G3H(J) = -128.82775D+00
        ELSEIF(NAZ.EQ.11) THEN
          G3H(J) = -161.84793D+00
        ELSEIF(NAZ.EQ.12) THEN
          G3H(J) = -199.65065D+00
        ELSEIF(NAZ.EQ.13) THEN
          G3H(J) = -241.93666D+00
        ELSEIF(NAZ.EQ.14) THEN
          G3H(J) = -288.93907D+00
        ELSEIF(NAZ.EQ.15) THEN
          G3H(J) = -340.82623D+00
        ELSEIF(NAZ.EQ.16) THEN
          G3H(J) = -397.66322D+00
        ELSEIF(NAZ.EQ.17) THEN
          G3H(J) = -459.68658D+00
        END IF
         HG3H(J) = G3H(J) + HTH
      END DO
C
      RETURN
C
      END
C
C*MODULE G3 *DECK FILATM
      SUBROUTINE FILATM(ATMFMN,HATFMN,SO)
C
C ATOM HEATS OF FORMATION FROM CCCBDB:
C HTTP://SRDATA.NIST.GOV/CCCBDB/
C SPIN ORBIT CORRECTIONS:
C ENTHALPIES: J. CHEM. PHYS. 106 (3), 15 JANUARY, 1997 (P.1063-1079)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000)
      DIMENSION ATMFMN(MXATM),HATFMN(MXATM),SO(MXATM)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      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
C--------------------------------------------------------------
C
C     HEATS OF FORMATION OF ATOM IS GIVEN BY NIST DATABASE
C        IN KCAL/MOL
C
C     DO I=1,NAT
C     NATM = INT(ZAN(I))
C        IF(NATM.EQ.1) THEN
C           IF(IDON.EQ.1) THEN
C              ATMFMN(I) = 50.62D+00
C           ELSE
C              ATMFMN(I) = 51.63D+00
C           END IF
C               SO(I) = 0.0D+00
C        ELSEIF(NATM.EQ.6) THEN
C            IF(IDON.EQ.1) THEN
C               ATMFMN(I) = 169.73D+00
C            ELSE
C               ATMFMN(I) = 169.98D+00
C            END IF
C               SO(I) = -0.14D+00
C        ELSEIF(NATM.EQ.7) THEN
C           ATMFMN(I) = 112.97D+00
C               SO(I) = 0.00D+00
C        ELSEIF(NATM.EQ.8) THEN
C           ATMFMN(I) = 59.553D+00
C               SO(I) = -0.36D+00
C        ELSE IF(NATM.EQ.9) THEN
C           ATMFMN(I) = 18.97D+00
C               SO(I) = -0.61D+00
C        ELSE IF(NATM.EQ.14) THEN
C           ATMFMN(I) = 107.55D+00
C               SO(I) = -0.68D+00
C        ELSE IF(NATM.EQ.15) THEN
C           ATMFMN(I) = 75.619D+00
C               SO(I) =  0.00D+00
C        ELSE IF(NATM.EQ.16) THEN
C           ATMFMN(I) = 66.200D+00
C               SO(I) = -0.89D+00
C        ELSE IF(NATM.EQ.17) THEN
C           ATMFMN(I) = 28.991D+00
C               SO(I) = -1.34D+00
C        ELSE
C           WRITE(6,1000) I
C        END IF
C     END DO
C
C     RETURN
C1000 FORMAT(/'NO SO COUPLING DATA FOR ATOM #',I3,/)
C     END
C
C--------------------------------------------------------------
C--------------------------------------------------------------
C ENTHALPIES: J. CHEM. PHYS. 106 (3), 15 JANUARY, 1997 (P.1063-1079)
C
C     DO I=1,NAT
C     NATM = INT(ZAN(I))
C        IF(NATM.EQ.1) THEN
C           ATMFMN(I) = 51.63D+00
C           HATFMN(I) = 52.64D+00
C               SO(I) = 0.0D+00
C        ELSEIF(NATM.EQ.6) THEN
C           ATMFMN(I) = 169.98D+00
C           HATFMN(I) = 170.23D+00
C               SO(I) = -0.14D+00
C        ELSEIF(NATM.EQ.7) THEN
C           ATMFMN(I) = 112.53D+00
C           HATFMN(I) = 113.57D+00
C               SO(I) = 0.00D+00
C        ELSEIF(NATM.EQ.8) THEN
C           ATMFMN(I) = 58.99D+00
C           HATFMN(I) = 60.03D+00
C               SO(I) = -0.36D+00
C        ELSE IF(NATM.EQ.9) THEN
C           ATMFMN(I) = 18.47D+00
C           HATFMN(I) = 19.52D+00
C               SO(I) = -0.61D+00
C        ELSE IF(NATM.EQ.14) THEN
C           ATMFMN(I) = 106.6D+00
C           HATFMN(I) = 107.36D+00
C               SO(I) = -0.68D+00
C        ELSE IF(NATM.EQ.15) THEN
C           ATMFMN(I) = 75.42D+00
C           HATFMN(I) = 76.70D+00
C               SO(I) =  0.00D+00
C        ELSE IF(NATM.EQ.16) THEN
C           ATMFMN(I) = 65.66D+00
C           HATFMN(I) = 66.71D+00
C               SO(I) = -0.89D+00
C        ELSE IF(NATM.EQ.17) THEN
C           ATMFMN(I) = 28.59D+00
C           HATFMN(I) = 29.69D+00
C               SO(I) = -1.34D+00
C        ELSE
C           WRITE(6,1000) I
C        END IF
C     END DO
C
C     RETURN
C1000 FORMAT(/'NO SO COUPLING DATA FOR ATOM #',I3,/)
C     END
C
C--------------------------------------------------------------
C--------------------------------------------------------------
C     REVISED NIST USING CCCBD VALUES
C
      DO I=1,NAT
      NATM = INT(ZAN(I))
         IF(NATM.EQ.1) THEN
            ATMFMN(I) = 51.63D+00
            HATFMN(I) = 52.103D+00
                SO(I) = 0.0D+00
         ELSEIF(NATM.EQ.6) THEN
            ATMFMN(I) = 169.98D+00
            HATFMN(I) = 171.29D+00
                SO(I) = -0.14D+00
         ELSEIF(NATM.EQ.7) THEN
            ATMFMN(I) = 112.53D+00
            HATFMN(I) = 112.97D+00
                SO(I) = 0.00D+00
         ELSEIF(NATM.EQ.8) THEN
            ATMFMN(I) = 58.99D+00
            HATFMN(I) = 59.56D+00
                SO(I) = -0.36D+00
         ELSE IF(NATM.EQ.9) THEN
            ATMFMN(I) = 18.47D+00
            HATFMN(I) = 18.97D+00
                SO(I) = -0.61D+00
         ELSE IF(NATM.EQ.14) THEN
            ATMFMN(I) = 106.52D+00
            HATFMN(I) = 107.55D+00
                SO(I) = -0.68D+00
         ELSE IF(NATM.EQ.15) THEN
            ATMFMN(I) = 75.45D+00
            HATFMN(I) = 75.65D+00
                SO(I) =  0.00D+00
         ELSE IF(NATM.EQ.16) THEN
            ATMFMN(I) = 65.71D+00
            HATFMN(I) = 66.25D+00
                SO(I) = -0.89D+00
         ELSE IF(NATM.EQ.17) THEN
            ATMFMN(I) = 28.59D+00
            HATFMN(I) = 28.99D+00
                SO(I) = -1.34D+00
         ELSE
            IF(MASWRK) WRITE(6,1000) I
         END IF
      END DO
C
      RETURN
 1000 FORMAT(/'NO SO COUPLING DATA FOR ATOM #',I3,/)
      END
C
C--------------------------------------------------------------
C
C*MODULE G3 *DECK GETH
      SUBROUTINE GETH(E0G,HG3,HG3H,G3H,ATMFMN,HATFMN,
     *                SO,HFOR0,HFOR298)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      PARAMETER (MXATM=2000,TM3=1.0D+03,TOKCAL=627.51D+00)
C
      COMMON /INFOA / NAT,ICH,MUL,NUM,NQMT,NE,NA,NB,
     *                ZAN(MXATM),C(3,MXATM),IAN(MXATM)
C
      DIMENSION HG3H(MXATM),G3H(MXATM),
     *       ATMFMN(MXATM),HATFMN(MXATM),SO(MXATM)
C
      E0ATMS = 0.0D+00
C     EHFOR0 = 0.0D+00
C     EHFOR298 = 0.0D+00
      HATMS  = 0.0D+00
      FATMS  = 0.0D+00
      FATMSH = 0.0D+00
      SOSUM  = 0.0D+00
C
      EMOL0  = E0G
      EMOLH  = HG3
C
      DO I=1,NAT
         E0ATMS = E0ATMS + G3H(I)
         HATMS  = HATMS  + HG3H(I)
         FATMS  = FATMS  + ATMFMN(I)
         FATMSH = FATMSH + HATFMN(I)
         SOSUM  = SOSUM  + SO(I)
      END DO
C
      SOSUM    = (SOSUM*TOKCAL)/TM3
C     EHFOR0   = TOKCAL*(EMOL0 - E0ATMS)
C     EHFOR298 = TOKCAL*(EMOLH - HATMS )
      HFOR0    = TOKCAL*(EMOL0 - E0ATMS) + FATMS
      HFOR298  = TOKCAL*(EMOLH - HATMS ) + FATMSH
C
      RETURN
      END
