(20 Feb 01)

                      *                               *
                      * Section 2 - Input Description *
                      *                               *

This is a copy of the input section of the GAMESS manual, including some specific changes made during the installation at the TU Graz.


$CONTRL group (optional)

          This is a free format group specifying global switches.

          SCFTYP             together with MPLEVL or CITYP specifies
                             the wavefunction.  You may choose from

                 = RHF       Restricted Hartree Fock calculation

                 = UHF       Unrestricted Hartree Fock calculation

                 = ROHF      Restricted open shell Hartree-Fock.
                             (high spin, see GVB for low spin)

                 = GVB       Generalized valence bond wavefunction
                             or OCBSE type ROHF. (needs $SCF input)

                 = MCSCF     Multiconfigurational SCF wavefunction
                             (this requires $DET or $DRT input)

                 = NONE      indicates a single point computation,
                             rereading a converged SCF function.
                             This option requires that you select
                             CITYP=GUGA or ALDET, RUNTYP=ENERGY
                             or TRANSITN, and GUESS=MOREAD.

          MPLEVL =           chooses Moller-Plesset perturbation
                             theory level, after the SCF.
                             See $MP2 and $MCQDPT input groups.
                 = 0         skips the MP computation (default)
                 = 2         performs a second order energy
                             correction.  MP2 is implemented only
                             for RHF, UHF, ROHF, and MCSCF wave
                             functions.  Gradients are available
                             only for RHF, so for the others you
                             may pick from RUNTYP=ENERGY, TRUDGE,
                             SURFACE, or FFIELD only.

          CITYP  =           chooses CI computation after the SCF,
                             for any SCFTYP except UHF.
                 = NONE      skips the CI. (default)
                 = GUGA      runs the Unitary Group CI package,
                             which requires $CIDRT input.
                             Gradients are available only for RHF,
                             so for other SCFTYPs, you may choose
                             only RUNTYP=ENERGY, TRUDGE, SURFACE,
                             FFIELD, TRANSITN.
                 = ALDET     runs the Ames Laboratory determinant
                             full CI package, requiring $CIDET
                             input.  RUNTYP=ENERGY only.

          Obviously, at most one of MPLEVL or CITYP may be chosen.

          RUNTYP             specifies the type of computation, for
                             example at a single geometry point:

                 = ENERGY    Molecular energy. (default)
                 = GRADIENT  Molecular energy plus gradient.
                 = HESSIAN   Molecular energy plus gradient plus
                             second derivatives, including harmonic
                             harmonic vibrational analysis.  See the
                             $FORCE and $CPHF input groups.

                             multiple geometry options:

                 = OPTIMIZE  Optimize the molecular geometry using
                             analytic energy gradients. See $STATPT.
                 = TRUDGE    Non-gradient total energy minimization.
                             See groups $TRUDGE and $TRURST.
                 = SADPOINT  Locate saddle point (transition state).
                             See the $STATPT group.
                 = IRC       Follow intrinsic reaction coordinate.
                             See the $IRC group.
                 = VSCF      Compute anharmonic vibrational
                             corrections (see $VSCF)
                 = DRC       Follow dynamic reaction coordinate.
                             See the $DRC group.
                 = GLOBOP    global optimization of effective fragment
                             positions via Monte Carlo.  See $GLOBOP.
                 = GRADEXTR  Trace gradient extremal.
                             See the $GRADEX group.
                 = SURFACE   Scan linear cross sections of the
                             potential energy surface.  See $SURF.

                             single geometry property options:

                 = PROP      Properties will be calculated.  A $DATA
                             deck and converged $VEC group should be
                             input.  Optionally, orbital localization
                             can be done.  See $ELPOT, etc.
                 = RAMAN     computes Raman intensities, see $RAMAN.
                 = MOROKUMA  Performs monomer energy decomposition.
                             See the $MOROKM group.
                 = TRANSITN  Compute radiative transition moment or
                             spin-orbit coupling.  See $TRANST group.
                 = FFIELD    applies finite electric fields, most
                             commonly to extract polarizabilities.
                             See the $FFCALC group.
                 = TDHF      analytic computation of time dependent
                             polarizabilities.  See the $TDHF group.
                 = MAKEFP    creates an effective fragment potential.

              * * * * * * * * * * * * * * * * * * * * * * * * * *
              Note that RUNTYPs involving the energy gradient,
              GLOBOP, IRC, GRADEXTR, and DRC, cannot be used for
              any CI or MP2 computation, except when SCFTYP=RHF.
              * * * * * * * * * * * * * * * * * * * * * * * * * *

          EXETYP = RUN       Actually do the run. (default)
                 = CHECK     Wavefunction and energy will not be
                             evaluated.  This lets you speedily
                             check input and memory requirements.
                             See the overview section for details.
                 = DEBUG     Massive amounts of output are printed,
                             useful only if you hate trees.
                 = routine   Maximum output is generated by the
                             routine named.  Check the source for
                             the routines this applies to.

          MAXIT  =           Maximum number of SCF iteration cycles.
                             Pertains only to RHF, UHF, ROHF, or
                             GVB runs.  See also MAXIT in $MCSCF.
                             (default = 30)

                           * * * * * * *

          ICHARG =           Molecular charge.  (default=0, neutral)

          MULT   =           Multiplicity of the electronic state
                 = 1         singlet (default)
                 = 2,3,...   doublet, triplet, and so on.

             ICHARG and MULT are used directly for RHF, UHF, ROHF.
             For GVB, these are implicit in the $SCF input, while
             for MCSCF or CI, these are implicit in $DRT/$CIDRT or
             $DET/$CIDET input.  You must still give them correctly.

                           * * * * * * *

          ECP    =           effective core potential control.
                 = NONE      all electron calculation (default).
                 = READ      read the potentials in $ECP group.
                 = SBKJC     use Stevens, Basch, Krauss, Jasien,
                             Cundari potentials for all heavy
                             atoms (Li-Rn are available).
                 = HW        use Hay, Wadt potentials for all the
                             heavy atoms (Na-Xe are available).

                           * * * * * * *

          RELWFN = NONE (default)  See also $RELWFN input group.
                 = NESC normalised elimination of small component,
                        the method of K. Dyall
                 = RESC relativistic elimination of small component,
                        the method of T. Nakajima and K. Hirao.

            * * * the next three control molecular geometry * * *

          COORD  = choice for molecular geometry in $DATA.
                 = UNIQUE    only the symmetry unique atoms will be
                             given, in Cartesian coords (default).
                 = HINT      only the symmetry unique atoms will be
                             given, in Hilderbrandt style internals.
                 = CART      Cartesian coordinates will be input.
                             Please read the warning just below!!!
                 = ZMT       GAUSSIAN style internals will be input.
                 = ZMTMPC    MOPAC style internals will be input.
                 = FRAGONLY  means no part of the system is treated
                             by ab initio means, hence $DATA is not
                             given.  The system is specified by $EFRAG.

            Note that the CART, ZMT, ZMTMPC choices require input of
            all atoms in the molecule.  These three also orient the
            molecule, and then determine which atoms are unique.  The
            reorientation is very likely to change the order of the
            atoms from what you input.  When the point group contains
            a 3-fold or higher rotation axis, the degenerate moments
            of inertia often cause problems choosing correct symmetry
            unique axes, in which case you must use COORD=UNIQUE
            rather than Z-matrices.

            Warning:  The reorientation into principal axes is done
            only for atomic coordinates, and is not applied to the
            axis dependent data in the following groups: $VEC, $HESS,
            $GRAD, $DIPDR, $VIB, nor Cartesian coords of effective
            fragments in $EFRAG.  COORD=UNIQUE avoids reorientation,
            and thus is the safest way to read these.

            Note that the choices CART, ZMT, ZMTMPC require the use
            of a $BASIS group to define the basis set.  The first
            two choices might or might not use $BASIS, as you wish.

          UNITS  = distance units, any angles must be in degrees.
                 = ANGS      Angstroms (default)
                 = BOHR      Bohr atomic units

          NZVAR  = 0  Use Cartesian coordinates (default).
                 = M  If COORD=ZMT or ZMTMPC and a $ZMAT is not given:
                      the internal coordinates will be those defining
                      the molecule in $DATA.  In this case, $DATA must
                      not contain any dummy atoms.  M is usually 3N-6,
                      or 3N-5 for linear.
                 = M  For other COORD choices, or if $ZMAT is given:
                      the internal coordinates will be those defined
                      in $ZMAT.  This allows more sophisticated
                      internal coordinate choices.  M is ordinarily
                      3N-6 (3N-5), unless $ZMAT has linear bends.

            NZVAR refers mainly to the coordinates used by OPTIMIZE
            or SADPOINT runs, but may also print the internal's
            values for other run types.  You can use internals to
            define the molecule, but Cartesians during optimizations!

          LOCAL  =          controls orbital localization.
                 = NONE     Skip localization (default).
                 = BOYS     Do Foster-Boys localization.
                 = RUEDNBRG Do Edmiston-Ruedenberg localization.
                 = POP      Do Pipek-Mezey population localization.
                            See the $LOCAL group.   Localization
                            does not work for SCFTYP=GVB or CITYP.

          ISPHER =      Spherical Harmonics option
                 = -1   Use Cartesian basis functions to construct
                        symmetry-adapted linear combination (SALC)
                        of basis functions.  The SALC space is the
                        linear variation space used.  (default)
                 = 0    Use spherical harmonic functions to create
                        SALC functions, which are then expressed
                        in terms of Cartesian functions.  The
                        contaminants are not dropped, hence this
                        option has EXACTLY the same variational
                        space as ISPHER=-1.  The only benefit to
                        obtain from this is a population analysis
                        in terms of pure s,p,d,f,g functions.
                 = +1   Same as ISPHER=0, but the function space
                        is truncated to eliminate all contaminant
                        Cartesian functions [3S(D), 3P(F), 4S(G),
                        and 3D(G)] before constructing the SALC
                        functions.  The computation corresponds
                        to the use of a spherical harmonic basis.

          QMTTOL = linear dependence threshhold
                   Any functions in the SALC variational space whose
                   eigenvalue of the overlap matrix is below this
                   tolerence is considered to be linearly dependent.
                   Such functions are dropped from the variational
                   space.  What is dropped is not individual basis
                   functions, but rather some linear combination(s)
                   of the entire basis set that represent the linear
                   dependent part of the function space.  The default
                   is a reasonable value for most purposes, 1.0E-6.

                   When many diffuse functions are used, it is common
                   to see the program drop some combinations.  On
                   occasion, in multi-ring molecules, we have raised
                   QMTTOL to 3.0E-6 to obtain SCF convergence, at the
                   cost of some energy.

                 * * * interfaces to other programs * * *

          MOLPLT = flag that produces an input deck for a molecule
                   drawing program distributed with GAMESS.
                   (default is .FALSE.)

          PLTORB = flag that produces an input deck for an orbital
                   plotting program distributed with GAMESS.
                   (default is .FALSE.)

          AIMPAC = flag to create an input deck for Bader's atoms
                   in molecules properties code. (default=.FALSE.)
                   For information about this program, see the URL

          FRIEND = string to prepare input to other quantum
                   programs, choose from
                 = HONDO    for HONDO 8.2
                 = MELDF    for MELDF
                 = GAMESSUK for GAMESS (UK Daresbury version)
                 = GAUSSIAN for Gaussian 9x
                 = ALL      for all of the above

          PLTORB, MOLPLT, and AIMPAC decks are written to file
          PUNCH at the end of the job.  Thus all of these correspond
          to the final geometry encountered during jobs such as

          In contrast, selecting FRIEND turns the job into a
          CHECK run only, no matter how you set EXETYP.  Thus the
          geometry is that encountered in $DATA.  The input is
          added to the PUNCH file, and may require some (usually
          minimal) massaging.

          PLTORB and MOLPLT are written even for EXETYP=CHECK.
          AIMPAC requires at least RUNTYP=PROP.

             The NBO program of Frank Weinhold's group can be
          attached to GAMESS.  The input to control the natural
          bond order analysis is read by the add in code, so is
          not described here.  The NBO program is available by
          anonymous FTP to ftp.osc.edu, in the directory

                 * * * computation control switches * * *

             For the most part, the default is the only sensible
          value, and unless you are sure of what you are doing,
          these probably should not be touched.

          NPRINT =           Print/punch control flag
                             See also EXETYP for debug info.
                             (options -7 to 5 are primarily debug)
                 = -7        Extra printing from Boys localization.
                 = -6        debug for geometry searches
                 = -5        minimal output
                 = -4        print 2e-contribution to gradient.
                 = -3        print 1e-contribution to gradient.
                 = -2        normal printing, no punch file
                 =  1        extra printing for basis,symmetry,ZMAT
                 =  2        extra printing for MO guess routines
                 =  3        print out property and 1e- integrals
                 =  4        print out 2e- integrals
                 =  5        print out SCF data for each cycle.
                             (Fock and density matrices, current MOs
                 =  6        same as 7, but wider 132 columns output.
                             This option isn't perfect.
                 =  7        normal printing and punching (default)
                 =  8        more printout than 7. The extra output
                             is (AO) Mulliken and overlap population
                             analysis, eigenvalues, Lagrangians, ...
                 =  9        everything in 8 plus Lowdin population
                             analysis, final density matrix.

          NOSYM  = 0     the symmetry specified in $DATA is used
                         as much as possible in integrals, SCF,
                         gradients, etc.  (this is the default)
                 = 1     the symmetry specified in the $DATA group
                         is used to build the molecule, then
                         symmetry is not used again.   Some GVB
                         or MCSCF runs (those without a totally
                         symmetric charge density) require you
                         request no symmetry.

          INTTYP = POPLE use fast Pople-Hehre routines for sp integral
                         blocks, and HONDO Rys polynomial code for
                         all other integrals.  (default)
                 = HONDO use HONDO/Rys integrals for all integrals.
                         This option produces slightly more accurate
                         integrals but is also slower.

                   When diffuse functions are used, the inaccuracy in
                   Pople/Hehre sp integrals shows up as inaccurate
                   LCAO coefficients in virtual orbitals.  This means
                   the error in SCF (meaning RHF to MCSCF) energies is
                   expected to be about 5d-8 Hartree, but the error in
                   computations that OCCUPY the virtual orbitals may
                   be much larger.  We have seen an energy error of
                   1d-4 in an MP2 energy when diffuse functions were
                   used.  We recommend that all MP2 or CI jobs with
                   diffuse functions select INTTYP=HONDO.

          NORMF  = 0     normalize the basis functions (default)
                 = 1     no normalization

          NORMP  = 0     input contraction coefficients refer to
                         normalized Gaussian primitives. (default)
                 = 1     the opposite.

          ITOL   =       primitive cutoff factor (default=20)
                 = n     products of primitives whose exponential
                         factor is less than 10**(-n) are skipped.

          ICUT   = n     integrals less than 10.0**(-n) are not
                         saved on disk. (default = 9)

                      * * * restart options * * *

          IREST  =       restart control options
                         (for OPTIMIZE run restarts, see $STATPT)
                         Note that this option is unreliable!
                 = -1    reuse dictionary file from previous run,
                         useful with GEOM=DAF and/or GUESS=MOSAVED.
                         Otherwise, this option is the same as 0.
                 = 0     normal run (default)
                 = 1     2e restart (1-e integrals and MOs saved)
                 = 2     SCF restart (1-,2-e integrls and MOs saved)
                 = 3     1e gradient restart
                 = 4     2e gradient restart

          GEOM   =       select where to obtain molecular geometry
                 = INPUT from $DATA input (default for IREST=0)
                 = DAF   read from DICTNRY file (default otherwise)

              As noted in the first chapter, binary file restart is
          not a well tested option!



$SYSTEM group (optional)

              This group provides global control information for
          your computer's operation.  This is system related input,
          and will not seem particularly chemical to you!

          TIMLIM =  time limit, in minutes.  Set to about 95 percent
                    of the time limit given to the batch job so that
                    GAMESS can stop itself gently.  (default=600.0)

          MWORDS =  the maximum replicated memory which your job can
                    use, on every node.  This is given in units of
                    1,000,000 words (as opposed to 1024*1024 words),
                    where a word is always a 64 bit quantity.  Most
                    systems allocate this memory at run time, but
                    some more primitive systems may have an upper
                    limit chosen at compile time.  (default=1)
                    In case finer control over the memory is needed,
                    this value can be given in units of words by
                    using the keyword MEMORY instead of MWORDS.

          MEMDDI =  the grand total memory needed for the distributed
                    data interface (DDI) storage, given in units of
                    1,000,000 words.  See Chapter 5 of this manual for
                    an extended explanation of running with MEMDDI.

          note: the memory required on each node for a run using
                p processors is therefore MWORDS + MEMDDI/p.

          PARALL =  a flag to cause the distributed data parallel
                    MP2 program to execute the parallel algorithm
                    even if you are running on only one node.
                    The main purpose of this is to allow you to
                    do EXETYP=CHECK runs to learn what the correct
                    value of MEMDDI needs to be.

          KDIAG  =    diagonalization control switch
                 = 0  use a vectorized diagonalization routine
                      if one is available on your machine,
                      else use EVVRSP. (default)
                 = 1  use EVVRSP diagonalization.  This may
                      be more accurate than KDIAG=0.
                 = 2  use GIVEIS diagonalization
                      (not as fast or reliable as EVVRSP)
                 = 3  use JACOBI diagonalization
                      (this is the slowest method)

          COREFL =  a flag to indicate whether or not GAMESS
                    should produce a "core" file for debugging
                    when subroutine ABRT is called to kill
                    a job.  This variable pertains only to
                    UNIX operating systems.  (default=.FALSE.)

          * * * the next three refer to parallel GAMESS * * *

          The next three apply only to parallel runs, and as they
          are more or less obsolete, their use is discourged.

          BALTYP =  Parallel load balence scheme
                    LOOP turns off dynamic load balancing (DLB)
                    NXTVAL uses dynamic load balancing
                    (default = LOOP)

          XDR    =  a flag to indicate whether or not messages
                    should be converted into a generic format
                    known as external data representation.
                    If true, messages can exchange between
                    machines of different vendors, at the cost
                    of performing the data type conversions.
                    (default=.FALSE.)  --inactive at present--

          PTIME  =  a logical flag to print extra timing info
                    during parallel runs.  This is not currently



$BASIS group (optional)

              This group allows certain standard basis sets to be
          easily given.  If this group is omitted, the basis set
          must be given instead in the $DATA group.

          GBASIS =        Name of the Gaussian basis set.
                 = MINI - Huzinaga's 3 gaussian minimal basis set.
                          Available H-Rn.
                 = MIDI - Huzinaga's 21 split valence basis set.
                          Available H-Rn.
                 = STO  - Pople's STO-NG minimal basis set.
                          Available H-Xe, for NGAUSS=2,3,4,5,6.
                 = N21  - Pople's N-21G split valence basis set.
                          Available H-Xe, for NGAUSS=3.
                          Available H-Ar, for NGAUSS=6.
                 = N31  - Pople's N-31G split valence basis set.
                          Available H-Ne,P-Cl for NGAUSS=4.
                          Available H-He,C-F for NGAUSS=5.
                          Available H-Ar, for NGAUSS=6.
                          For Ga-Kr, N31 selects the BC basis.
                 = N311 - Pople's "triple split" N-311G basis set.
                          Available H-Ne, for NGAUSS=6.
                          Selecting N311 implies MC for Na-Ar.
                 = DZV  - "double zeta valence" basis set.
                          a synonym for DH for H,Li,Be-Ne,Al-Cl.
                          (14s,9p,3d)/[5s,3p,1d] for K-Ca.
                          (14s,11p,5d/[6s,4p,1d] for Ga-Kr.
                 = DH   - Dunning/Hay "double zeta" basis set.
                          (3s)/[2s] for H.
                          (9s,4p)/[3s,2p] for Li.
                          (9s,5p)/[3s,2p] for Be-Ne.
                          (11s,7p)/[6s,4p] for Al-Cl.
                 = TZV  - "triple zeta valence" basis set.
                          (5s)/[3s] for H.
                          (10s,3p)/[4s,3p] for Li.
                          (10s,6p)/[5s,3p] for Be-Ne.
                          a synonym for MC for Na-Ar.
                          (14s,9p)/[8s,4p] for K-Ca.
                          (14s,11p,6d)/[10s,8p,3d] for Sc-Zn.
                 = MC   - McLean/Chandler "triple split" basis.
                          (12s,9p)/[6s,5p] for Na-Ar.
                          Selecting MC implies 6-311G for H-Ne.

              additional values for GBASIS are on the next page.

               * * * the next two are ECP bases only * * *

          GBASIS = SBKJC- Stevens/Basch/Krauss/Jasien/Cundari
                          valence basis set, for Li-Rn.  This choice
                          implies an unscaled -31G basis for H-He.
                 = HW   - Hay/Wadt valence basis.
                          This is a -21 split, available Na-Xe,
                          except for the transition metals.
                          This implies a 3-21G basis for H-Ne.

               * * * semiempirical basis sets * * *

                   The elements for which these exist can be found
                   in the 'further information' section of this
                   manual.  If you pick one of these, all other data
                   in this group is ignored.  Semi-empirical runs
                   actually use valence-only STO bases, not GTOs.

          GBASIS = MNDO - selects MNDO model hamiltonian

                 = AM1  - selects AM1 model hamiltonian

                 = PM3  - selects PM3 model hamiltonian

          NGAUSS = the number of Gaussians (N).   This parameter
                   pertains only to GBASIS=STO, N21, N31, or N311.

          NDFUNC = number of heavy atom polarization functions to
                   be used.  These are usually d functions, except
                   for MINI/MIDI.  The term "heavy" means Na on up
                   when GBASIS=STO, HW, or N21, and from Li on up
                   otherwise.  The value may not exceed 3.  The
                   variable POLAR selects the actual exponents to
                   be used, see also SPLIT2 and SPLIT3. (default=0)

          NFFUNC = number of heavy atom f type polarization
                   functions to be used on Li-Cl.  This may only
                   be input as 0 or 1.  (default=0)

          NPFUNC = number of light atom, p type polarization
                   functions to be used on H-He.  This may not
                   exceed 3, see also POLAR.  (default=0)

          DIFFSP = flag to add diffuse sp (L) shell to heavy atoms.
                   Heavy means Li-F, Na-Cl, Ga-Br, In-I, Tl-At.
                   The default is .FALSE.

          DIFFS  = flag to add diffuse s shell to hydrogens.
                   The default is .FALSE.

          Warning: if you use diffuse functions, please read QMTTOL
          and INTTYP in the $CONTRL group for numerical concerns.

          POLAR  = exponent of polarization functions
                 = POPLE     (default for all other cases)
                 = POPN311   (default for GBASIS=N311, MC)
                 = DUNNING   (default for GBASIS=DH, DZV)
                 = HUZINAGA  (default for GBASIS=MINI, MIDI)
                 = HONDO7    (default for GBASIS=TZV)

          SPLIT2 = an array of splitting factors used when NDFUNC
                   or NPFUNC is 2.  Default=2.0,0.5

          SPLIT3 = an array of splitting factors used when NDFUNC
                   or NPFUNC is 3.  Default=4.00,1.00,0.25


          The splitting factors are from the Pople school, and are
          probably too far apart.  See for example the Binning and
          Curtiss paper.  For example, the SPLIT2 value will usually
          cause an INCREASE over the 1d energy at the HF level for

          The actual exponents used for polarization functions, as
          well as for diffuse sp or s shells, are described in the
          'Further References' section of this manual.  This section
          also describes the sp part of the basis set chosen by
          GBASIS fully, with all references cited.

          Note that GAMESS always punches a full $DATA group.  Thus,
          if $BASIS does not quite cover the basis you want, you can
          obtain this full $DATA group from EXETYP=CHECK, and then
          change polarization exponents, add Rydbergs, etc.



$DATA group (required)

$DATAS group (if NESC chosen, gives small component basis)

$DATAL group (if NESC chosen, gives large component basis)

              This group describes the global molecular data such as
          point group symmetry, nuclear coordinates, and possibly
          the basis set.  It consists of a series of free format
          card images.  See $RELWFN for more informaton on large and
          small component basis sets.  The input structure of $DATAS
          and $DATAL is identical to the COORD=UNIQUE $DATA input.


          -1-   TITLE     a single descriptive title card.


          -2-   GROUP, NAXIS

          GROUP is the Schoenflies symbol of the symmetry group,
          you may choose from
              C1, CS, CI, CN, S2N, CNH, CNV, DN, DNH, DND,
              T, TH, TD, O, OH.

          NAXIS is the order of the highest rotation axis, and
          must be given when the name of the group contains an N.
          For example, "Cnv 2" is C2v.  "S2n 3" means S6.

          For linear molecules, choose either CNV or DNH, and enter
          NAXIS as 4.  Enter atoms as DNH with NAXIS=2.  If the
          electronic state of either is degenerate, check the note
          about the effect of symmetry in the electronic state
          in the SCF section of REFS.DOC.


              In order to use GAMESS effectively, you must be able
          to recognize the point group name for your molecule.  This
          presupposes a knowledge of group theory at about the level
          of Cotton's "Group Theory", Chapter 3.

              Armed with only the name of the group, GAMESS is able
          to exploit the molecular symmetry throughout almost all of
          the program, and thus save a great deal of computer time.
          GAMESS does not require that you know very much else about
          group theory, although a deeper knowledge (character
          tables, irreducible representations, term symbols, and so
          on) is useful when dealing with the more sophisticated

          Cards -3- and -4- are quite complicated, and are rarely
          given.  A *SINGLE* blank card may replace both cards -3-
          and -4-, to select the 'master frame', which is defined on
          the next page.   If you choose to enter a blank card, skip
          to the bottom of the next page.

          If the point group is C1 (no symmetry), skip over cards
          -3- and -4- (which means no blank card).


          -3-  X1, Y1, Z1, X2, Y2, Z2

          For C1 group, there is no card -3- or -4-.
          For CI group, give one point, the center of inversion.
          For CS group, any two points in the symmetry plane.
          For axial groups, any two points on the principal axis.
          For tetrahedral groups, any two points on a two-fold axis.
          For octahedral groups, any two points on a four-fold axis.


          -4-  X3, Y3, Z3, DIRECT

          third point, and a directional parameter.
          For CS group, one point of the symmetry plane,
                        noncollinear with points 1 and 2.
          For CI group, there is no card -4-.

          For other groups, a generator sigma-v plane (if any) is
          the (x,z) plane of the local frame (CNV point groups).

          A generator sigma-h plane (if any) is the (x,y) plane of
          the local frame (CNH and dihedral groups).

          A generator C2 axis (if any) is the x-axis of the local
          frame (dihedral groups).

          The perpendicular to the principal axis passing through
          the third point defines a direction called D1.  If
          DIRECT='PARALLEL', the x-axis of the local frame coincides
          with the direction D1.  If DIRECT='NORMAL', the x-axis of
          the local frame is the common perpendicular to D1 and the
          principal axis, passing through the intersection point of
          these two lines.  Thus D1 coincides in this case with the
          negative y axis.


              The 'master frame' is just a standard orientation for
          the molecule.  By default, the 'master frame' assumes that
              1.   z is the principal rotation axis (if any),
              2.   x is a perpendicular two-fold axis (if any),
              3.  xz is the sigma-v plane (if any), and
              4.  xy is the sigma-h plane (if any).
          Use the lowest number rule that applies to your molecule.

                  Some examples of these rules:
          Ammonia (C3v): the unique H lies in the XZ plane (R1,R3).
          Ethane (D3d): the unique H lies in the YZ plane (R1,R2).
          Methane (Td): the H lies in the XYZ direction (R2).  Since
                   there is more than one 3-fold, R1 does not apply.
          HP=O (Cs): the mirror plane is the XY plane (R4).

          In general, it is a poor idea to try to reorient the
          molecule.  Certain sections of the program, such as the
          orbital symmetry assignment, do not know how to deal with
          cases where the 'master frame' has been changed.

          Linear molecules (C4v or D4h) must lie along the z axis,
          so do not try to reorient linear molecules.

          You can use EXETYP=CHECK to quickly find what atoms are
          generated, and in what order.  This is typically necessary
          in order to use the general $ZMAT coordinates.

                               * * * *

          Depending on your choice for COORD in $CONTROL,

              if COORD=UNIQUE, follow card sequence U
              if COORD=HINT,   follow card sequence U
              if COORD=CART,   follow card sequence C
              if COORD=ZMT,    follow card sequence G
              if COORD=ZMTMPC, follow card sequence M

          Card sequence U is the only one which allows you to define
          a completely general basis here in $DATA.

          Recall that UNIT in $CONTRL determines the distance units.

          -5U-   Atom input.  Only the symmetry unique atoms are
          input, GAMESS will generate the symmetry equivalent atoms
          according to the point group selected above.

             if COORD=UNIQUE   NAME, ZNUC, X, Y, Z

          NAME  = 10 character atomic name, used only for printout.
                  Thus you can enter H or Hydrogen, or whatever.
          ZNUC  = nuclear charge.  It is the nuclear charge which
                  actually defines the atom's identity.
          X,Y,Z = Cartesian coordinates.

             if COORD=HINT


          NAME = 10 character atomic name (used only for print out).
          ZNUC = nuclear charge.
          CONX = connection type, choose from
            'LC'   linear conn.               'CCPA' central conn.
            'PCC'  planar central conn.              with polar atom
            'NPCC' non-planar central conn.   'TCT'  terminal conn.
            'PTC'  planar terminal conn.             with torsion
          R    = connection distance.
          ALPHA= first connection angle
          BETA = second connection angle
          SIGN = connection sign, '+' or '-'
          POINT1, POINT2, POINT3 =
               connection points, a serial number of a previously
               input atom, or one of 4 standard points: O,I,J,K
               (origin and unit points on axes of master frame).
               defaults:  POINT1='O', POINT2='I', POINT3='J'

          ref- R.L. Hilderbrandt, J.Chem.Phys. 51, 1654 (1969).
          You cannot understand HINT input without reading this.

          Note that if ZNUC is negative, the internally stored
          basis for ABS(ZNUC) is placed on this center, but the
          calculation uses ZNUC=0 after this.  This is useful
          for basis set superposition error (BSSE) calculations.

          * * * If you gave $BASIS, continue entering cards -5U-
                until all the unique atoms have been specified.
                When you are done, enter a " $END " card.
          * * * If you did not, enter cards -6U-, -7U-, -8U-.

          -6U-  GBASIS, NGAUSS, (SCALF(i),i=1,4)

          GBASIS has exactly the same meaning as in $BASIS.  You may
          choose from MINI, MIDI, STO, N21, N31, N311, DZV, DH, BC,
          TZV, MC, SBKJC, or HW.  In addition, you may choose S, P,
          D, F, G, or L to enter an explicit basis set.  Here, L
          means both an s and p shell with a shared exponent.

          NGAUSS is the number of Gaussians (N) in the Pople style
          basis, or user input general basis.  It has meaning only
          for GBASIS=STO, N21, N31, or N311, and S,P,D,F,G, or L.

          Up to four scale factors may be entered.  If omitted,
          standard values are used.  They are not documented as
          every GBASIS treats these differently.  Read the source
          code if you need to know more.  They are seldom given.

          * * * If GBASIS is not S,P,D,F,G, or L, either add more
                shells by repeating card -6U-, or go on to -8U-.
          * * * If GBASIS=S,P,D,F,G, or L, enter NGAUSS cards -7U-.

          -7U- IG, ZETA, C1, C2

                IG = a counter, IG takes values 1, 2, ..., NGAUSS.
              ZETA = Gaussian exponent of the IG'th primitive.
                C1 = Contraction coefficient for S,P,D,F,G shells,
                     and for the s function of L shells.
                C2 = Contraction coefficient for the p in L shells.

          * * * For more shells on this atom, go back to card -6U-.
          * * * If there are no more shells, go on to card -8U-.

          -8U-    A blank card ends the basis set for this atom.

          Continue entering atoms with -5U- through -8U- until all
          are given, then terminate the group with a " $END " card.

                 --- this is the end of card sequence U ---

          COORD=CART input:


          -5C- Atom input.

          Cartesian coordinates for all atoms must be entered.  They
          may be arbitrarily rotated or translated, but must possess
          the actual point group symmetry.  GAMESS will reorient the
          molecule into the 'master frame', and determine which
          atoms are the unique ones.  Thus, the final order of the
          atoms may be different from what you enter here.

                NAME, ZNUC, X, Y, Z

          NAME  = 10 character atomic name, used only for printout.
                  Thus you can enter H or Hydrogen, or whatever.
          ZNUC  = nuclear charge.  It is the nuclear charge which
                  actually defines the atom's identity.
          X,Y,Z = Cartesian coordinates.


          Continue entering atoms with card -5C- until all are
          given, and then terminate the group with a " $END " card.

                 --- this is the end of card sequence C ---

          COORD=ZMT input:       (GAUSSIAN style internals)


          -5G-      ATOM

          Only the name of the first atom is required.
          See -8G- for a description of this information.

          -6G-      ATOM  i1 BLENGTH

          Only a name and a bond distance is required for atom 2.
          See -8G- for a description of this information.

          -7G-      ATOM  i1 BLENGTH  i2 ALPHA

          Only a name, distance, and angle are required for atom 3.
          See -8G- for a description of this information.

          -8G-      ATOM  i1 BLENGTH  i2 ALPHA  i3 BETA i4

          ATOM    is the chemical symbol of this atom.  It can be
                  followed by numbers, if desired, for example Si3.
                  The chemical symbol implies the nuclear charge.
          i1      defines the connectivity of the following bond.
          BLENGTH is the bond length "this atom-atom i1".
          i2      defines the connectivity of the following angle.
          ALPHA   is the angle "this atom-atom i1-atom i2".
          i3      defines the connectivity of the following angle.
          BETA    is either the dihedral angle "this atom-atom i1-
                  atom i2-atom i3", or perhaps a second bond
                  angle "this atom-atom i1-atom i3".
          i4      defines the nature of BETA,
                  If BETA is a dihedral angle, i4=0 (default).
                  If BETA is a second bond angle, i4=+/-1.
                  (sign specifies one of two possible directions).

           o  Repeat -8G- for atoms 4, 5, ...
           o  The use of ghost atoms is possible, by using X or BQ
              for the chemical symbol.  Ghost atoms preclude the
              option of an automatic generation of $ZMAT.
           o  The connectivity i1, i2, i3 may be given as integers,
              1, 2, 3, 4, 5,...  or as strings which match one of
              the ATOMs.  In this case, numbers must be added to the
              ATOM strings to ensure uniqueness!

           o  In -6G- to -8G-, symbolic strings may be given in
              place of numeric values for BLENGTH, ALPHA, and BETA.
              The same string may be repeated, which is handy in
              enforcing symmetry.  If the string is preceeded by a
              minus sign, the numeric value which will be used is
              the opposite, of course.  Any mixture of numeric data
              and symbols may be given.  If any strings were given
              in -6G- to -8G-, you must provide cards -9G- and
              -10G-, otherwise you may terminate the group now with
              a " $END " card.


          -9G-   A blank line terminates the Z-matrix section.


          -10G-   STRING VALUE

          STRING is a symbolic string used in the Z-matrix.
          VALUE  is the numeric value to substitute for that string.


          Continue entering -10G- until all STRINGs are defined.
          Note that any blank card encountered while reading -10G-
          will be ignored.  GAMESS regards all STRINGs as variables
          (constraints are sometimes applied in $STATPT).  It is not
          necessary to place constraints to preserve point group
          symmetry, as GAMESS will never lower the symmetry from
          that given at -2-.  When you have given all STRINGs a
          VALUE, terminate the group with a " $END " card.

                 --- this is the end of card sequence G ---

                                * * * *

              The documentation for sequence G above and sequence M
          below presumes you are reasonably familiar with the input
          to GAUSSIAN or MOPAC.  It is probably too terse to be
          understood very well if you are unfamiliar with these.  A
          good tutorial on both styles of Z-matrix input can be
          found in Tim Clark's book "A Handbook of Computational
          Chemistry", published by John Wiley & Sons, 1985.

              Both Z-matrix input styles must generate a molecule
          which possesses the symmetry you requested at -2-.  If
          not, your job will be terminated automatically.

          COORD=ZMTMPC input:       (MOPAC style internals)


          -5M-     ATOM

          Only the name of the first atom is required.
          See -8M- for a description of this information.

          -6M-     ATOM BLENGTH

          Only a name and a bond distance is required for atom 2.
          See -8M- for a description of this information.

          -7M-     ATOM BLENGTH j1 ALPHA j2

          Only a bond distance from atom 2, and an angle with repect
          to atom 1 is required for atom 3.  If you prefer to hook
          atom 3 to atom 1, you must give connectivity as in -8M-.
          See -8M- for a description of this information.

          -8M-     ATOM BLENGTH j1 ALPHA j2 BETA j3 i1 i2 i3

          ATOM, BLENGTH, ALPHA, BETA, i1, i2 and i3 are as described
          at -8G-.  However, BLENGTH, ALPHA, and BETA must be given
          as numerical values only.  In addition, BETA is always a
          dihedral angle.   i1, i2, i3 must be integers only.

          The j1, j2 and j3 integers, used in MOPAC to signal
          optimization of parameters, must be supplied but are
          ignored here.  You may give them as 0, for example.

          Continue entering atoms 3, 4, 5, ... with -8M- cards until
          all are given, and then terminate the group by giving a
          " $END " card.

                 --- this is the end of card sequence M ---

                         This is the end of $DATA!

          If you have any doubt about what molecule and basis set
          you are defining, or what order the atoms will be
          generated in, simply execute an EXETYP=CHECK job to find



$ZMAT group (required if NZVAR is nonzero in $CONTRL)

              This group lets you define the internal coordinates in
          which the gradient geometry search is carried out.  These
          need not be the same as the internal coordinates used in
          $DATA.  The coordinates may be simple Z-matrix types,
          delocalized coordinates, or natural internal coordinates.

              You must input a total of M=3N-6 internal coordinates
          (M=3N-5 for linear molecules).  NZVAR in $CONTRL can be
          less than M IF AND ONLY IF you are using linear bends.  It
          is also possible to input more than M coordinates if they
          are used to form exactly M linear combinations for new
          internals.  These may be symmetry coordinates or natural
          internal coordinates.  If NZVAR > M, you must input IJS and
          SIJ below to form M new coordinates.  See DECOMP in $FORCE
          for the only circumstance in which you may enter a larger
          NZVAR without giving SIJ and IJS.

             **** IZMAT defines simple internal coordinates ****

          IZMAT is an array of integers defining each coordinate.
          The general form for each internal coordinate is
                code number,I,J,K,L,M,N

          IZMAT =1 followed by two atom numbers. (I-J bond length)
                =2 followed by three numbers. (I-J-K bond angle)
                =3 followed by four numbers. (dihedral angle)
                   Torsion angle between planes I-J-K and J-K-L.
                =4 followed by four atom numbers. (atom-plane)
                   Out-of-plane angle from bond I-J to plane J-K-L.
                =5 followed by three numbers. (I-J-K linear bend)
                   Counts as 2 coordinates for the degenerate bend,
                   normally J is the center atom.  See $LIBE.
                =6 followed by five atom numbers. (dihedral angle)
                   Dihedral angle between planes I-J-K and K-L-M.
                =7 followed by six atom numbers. (ghost torsion)
                   Let A be the midpoint between atoms I and J, and
                   B be the midpoint between atoms M and N.  This
                   coordinate is the dihedral angle A-K-L-B.  The
                   atoms I,J and/or M,N may be the same atom number.
                   (If I=J AND M=N, this is a conventional torsion).
                   Examples: N2H4, or, with one common pair, H2POH.

          Example - a nonlinear triatomic, atom 2 in the middle:
                $ZMAT IZMAT(1)=1,1,2,  2,1,2,3,  1,2,3  $END
          This sets up two bonds and the angle between them.
          The blanks between each coordinate definition are
          not necessary, but improve readability mightily.

               **** the next define delocalized coordinates ****

          DLC    is a flag to request delocalized coordinates.
                 (default is .FALSE.)

          AUTO   is a flag to generate all redundant coordinates,
                 automatically.  The DLC space will consist of all
                 non-redundant combinations of these which can be
                 found.  The list of redundant coordinates will
                 consist of bonds, angles, and torsions only.
                 (default is .FALSE.)

          NONVDW is an array of atom pairs which are to be joined
                 by a bond, but might be skipped by the routine
                 that automatically includes all distances shorter
                 than the sum of van der Waals radii.  Any angles
                 and torsions associated with the new bond(s) are
                 also automatically included.

          The format for IXZMAT, IRZMAT, IFZMAT is that of IZMAT:

          IXZMAT is an extra array of simple internal coordinates
                 which you want to have added to the list generated
                 by AUTO.  Unlike NONVDW, IXZMAT will add only the
                 coordinate(s) you specify.

          IRZMAT is an array of simple internal coordinates which
                 you would like to remove from the AUTO list of
                 redundant coordinates.  It is sometimes necessary
                 to remove a torsion if other torsions around a bond
                 are being frozen, to obtain a nonsingular G matrix.

          IFZMAT is an array of simple internal coordinates which
                 you would like to freeze.  See also FVALUE below.
                 Note that IFZMAT/FVALUE work only with DLC, see the
                 IFREEZ option in $STATPT to freeze coordinates if
                 you wish to freeze simple or natural coordinates.

          FVALUE is an array of values to which the internal
                 coordinates should be constrained.  It is not
                 necessary to input $DATA such that the initial
                 values match these desired final values, but it is
                 helpful if the initial values are not too far away.

             **** SIJ,IJS define natural internal coordinates ****

          SIJ is a transformation matrix of dimension NZVAR x M,
              used to transform the NZVAR internal coordinates in
              IZMAT into M new internal coordinates.  SIJ is a
              sparse matrix, so only the non-zero elements are
              given, by using the IJS array described below.
              The columns of SIJ will be normalized by GAMESS.
              (Default: SIJ = I, unit matrix)

          IJS is an array of pairs of indices, giving the row and
              column index of the entries in SIJ.

          example - if the above triatomic is water, using
               IJS(1) = 1,1, 3,1,   1,2, 3,2,   2,3
               SIJ(1) = 1.0, 1.0,   1.0,-1.0,   1.0

              gives the matrix S=  1.0   1.0   0.0
                                   0.0   0.0   1.0
                                   1.0  -1.0   0.0

          which defines the symmetric stretch, asymmetric stretch,
          and bend of water.

          references for natural internal coordinates:
            P.Pulay, G.Fogarasi, F.Pang, J.E.Boggs
               J.Am.Chem.Soc. 101, 2550-2560(1979)
            G.Fogarasi, X.Zhou, P.W.Taylor, P.Pulay
               J.Am.Chem.Soc. 114, 8191-8201(1992)
          reference for delocalized coordinates:
            J.Baker, A. Kessi, B.Delley
               J.Chem.Phys. 105, 192-212(1996)



$LIBE group (required if linear bends are used in $ZMAT)

          A degenerate linear bend occurs in two orthogonal planes,
          which are specified with the help of a point A.  The first
          bend occurs in a plane containing the atoms I,J,K and the
          user input point A.  The second bend is in the plane
          perpendicular to this, and containing I,J,K.  One such
          point must be given for each pair of bends used.

          APTS(1)= x1,y1,z1,x2,y2,z2,...  for linear bends 1,2,...

          Note that each linear bend serves as two coordinates, so
          that if you enter 2 linear bends (HCCH, for example), the
          correct value of NZVAR is M-2, where M=3N-6 or 3N-5, as



$SCF group relevant if SCFTYP = RHF, UHF, or ROHF, required if SCFTYP = GVB)

              This group of parameters provides additional control
          over the RHF, UHF, ROHF, or GVB SCF steps.  It must be
          used for GVB open shell or perfect pairing wavefunctions.

          DIRSCF = a flag to activate a direct SCF calculation,
                   which is implemented for all the Hartree-Fock
                   type wavefunctions:  RHF, ROHF, UHF, and GVB.
                   This keyword also selects direct MP2 computation.
                   .FALSE. stores integrals on disk storage for a 
                   conventional SCF calculation.
                   The default value at the TU Graz is .TRUE. .

          FDIFF  = a flag to compute only the change in the Fock
                   matrices since the previous iteration, rather
                   than recomputing all two electron contributions.
                   This saves much CPU time in the later iterations.
                   This pertains only to direct SCF, and has a
                   default of .TRUE.  This option is implemented
                   only for the RHF, ROHF, UHF cases.

                   Cases with many diffuse functions in the basis
                   set sometimes oscillate at the end, rather than
                   converging.  Turning this parameter off will
                   normally give convergence.

          ---- The next flags affect convergence rates.

          EXTRAP = controls Pople extrapolation of the Fock matrix.
          DAMP   = controls Davidson damping of the Fock matrix.
          SHIFT  = controls level shifting of the Fock matrix.
          RSTRCT = controls restriction of orbital interchanges.
          DIIS   = controls Pulay's DIIS interpolation.
          SOSCF  = controls second order SCF orbital optimization.
                   (default=.TRUE. for RHF, Abelian group ROHF, GVB)
                   (default=.FALSE. for UHF, non-Abelian group ROHF)
          DEM    = controls direct energy minimization, which is
                   implemented only for RHF.  (default=.FALSE.)

          defaults for     EXTRAP  DAMP  SHIFT RSTRCT  DIIS  SOSCF
          ab initio:         T      F      F      F      T    T/F
          semiempirical:     T      F      F      F      F     F

               The above parameters are implemented for all SCF
          wavefunction types, except that DIIS will work for GVB
          only for those cases with NPAIR=0 or NPAIR=1.  If both
          DIIS and SOSCF are chosen, SOSCF is stronger than DIIS,
          and so DIIS will not be used.

               Once either DIIS or SOSCF are initiated, any other
          accelerator in effect is put in abeyance.

          ---- These parameters fine tune the various convergers.

          NCONV  = n  SCF density convergence criteria.
                   Convergence is reached when the density change
                   between two consecutive SCF cycles is less than
                   10.0**(-n) in absolute value.  One more cycle
                   is executed after reaching convergence.   Less
                   accuracy in NCONV gives questionable gradients.
                   (default is n=5, except CI or MP2 gradients n=6)

          SOGTOL = second order gradient tolerance.  SOSCF will be
                   initiated when the orbital gradient falls below
                   this threshold.  (default=0.25 au)

          ETHRSH = energy error threshold for initiating DIIS.  The
                   DIIS error is the largest element of e=FDS-SDF.
                   Increasing ETHRSH forces DIIS on sooner.
                   (default = 0.5 Hartree)

          MAXDII = Maximum size of the DIIS linear equations, so
                   that at most MAXDII-1 Fock matrices are used
                   in the interpolation.  (default=10)

          DEMCUT = Direct energy minimization will not be done
                   once the density matrix change falls below
                   this threshold.  (Default=0.5)

          DMPCUT = Damping factor lower bound cutoff.  The damping
                   damping factor will not be allowed to drop
                   below this value. (default=0.0)
             note: The damping factor need not be zero to achieve
                   valid convergence (see Hsu, Davidson, and
                   Pitzer, J.Chem.Phys., 65, 609 (1976), see
                   especially the section on convergence control),
                   but it should not be astronomical either.

                  ----- miscellaneous options -----

          UHFNOS = flag controlling generation of the natural
                   orbitals of a UHF function. (default=.FALSE.)

          MVOQ   = 0  Skip MVO generation (default)
                 = n  Form modified virtual orbitals, using a cation
                      with n electrons removed.   Implemented for
                      RHF, ROHF, and GVB.   If necessary to reach a
                      closed shell cation, the program might remove
                      n+1 electrons.  Typically, n will be about 6.

          NPUNCH = SCF punch option
                 =  0  do not punch out the final orbitals
                 =  1  punch out the occupied orbitals
                 =  2  punch out occupied and virtual orbitals
                       The default is NPUNCH = 2.

                  ----- options for virial scaling -----

          VTSCAL =   A flag to request that the virial theorem be
                     satisfied.  An analysis of the total energy
                     as an exact sum of orbital kinetic energies
                     is printed.  The default is .FALSE.

             This option is implemented for RHF, UHF, and ROHF,
             for RUNTYP=ENERGY, OPTIMIZE, or SADPOINT.  Related
             input is as follows:

          SCALF  =   initial exponent scale factor when VTSCAL is
                     in use, useful when restarting.  The default
                     is 1.0.

          MAXVT  =   maximum number of iterations (at a single
                     geometry) to satisfy the energy virial theorem.
                     The default is 20.

          VTCONV =   convergence criterion for the VT, which is
                     satisfied when 2 +  + R x dE/dR is less
                     than VTCONV.  The default is 1.0D-6 Hartree.

          For more information on this option, which is most
          economically employed during a geometry search, see
          M.Lehd and F.Jensen, J.Comput.Chem. 12, 1089-1096(1991).

              The next parameters define the GVB wavefunction.  Note
          that ALPHA and BETA also have meaning for ROHF.  See also
          MULT in the $CONTRL group.  The GVB wavefunction assumes
          orbitals are in the order core, open, pairs.

          NCO    =   The number of closed shell orbitals.  The
                     default almost certainly should be changed!

          NSETO  =   The number of sets of open shells in the
                     function.  Maximum of 10. (default=0)

          NO     =   An array giving the degeneracy of each open
                     shell set.  Give NSETO values.

          NPAIR  =   The number of geminal pairs in the -GVB-
                     function.  Maximum of 12.  The default
                     corresponds to open shell SCF (default=0).

          CICOEF =   An array of ordered pairs of CI coefficients
                     for the -GVB- pairs.  For example, a two pair
                     case for water, say, might be
                     CICOEF(1)=0.95,-0.05,0.95,-0.05.  If not
                     normalized, as in the default, they will be.
                     This parameter is useful in restarting a GVB
                     run, with the current CI coefficients.
                     (default = 0.90,-0.20,0.90,-0.20,...)

          COUPLE =   A switch controlling the input of F, ALPHA,
                     and BETA.  The default is to use internally
                     stored values for these variables.   Note
                     ALPHA and BETA can be given for -ROHF-, as
                     well as -GVB-.  (Default=.FALSE.)

          F      =   An vector of fractional occupations.

          ALPHA  =   An array of A coupling coefficients given in
                     lower triangular order.

          BETA   =   An array of B coupling coefficients given in
                     lower triangular order.

              Note:  The default for F, ALPHA, and BETA depends on
          the state chosen.  Defaults for the most commonly occuring
          cases are internally stored.



$SCFMI group (optional, relevant if SCFTYP=RHF)

              The SCF-MI method is a modification of the Roothaan
          equations that avoids basis set superposition error (BSSE)
          in intermolecular interaction calculations, by expanding
          each monomer's orbitals using only its own basis set.
          Thus, the resulting orbitals are not orthogonal.  The
          presence of a $SCFMI group in the input triggers the use
          of this option.

              The implementation is limited to two monomers, treated
          at the RHF level.  The energy, gradient, and therefore
          numerical hessian are available.  The SCF step may be run
          in direct SCF mode.  The first 4 parameters must be given.
          All atoms of monomer A must be given in $DATA before the
          atoms of monomer B.

          NA        = number of doubly occupied MOs on fragment A.
          NB        = number of doubly occupied MOs on fragment B.
          MA        = number of basis functions on fragment A.
          MB        = number of basis functions on fragment B.

          ITER      = maximum number of SCF-MI cycles, overriding
                      the usual MAXIT value.  (default is 50).

          DTOL      = SCF-MI density convergence criteria.
                      (default is 1.0d-10)

          ALPHA     = possible level shift parameter.
                      (default is 0.0, meaning shifting is not used)

          IOPT      =   prints additional debug information.
                    = 0 standard outout (default)
                    = 1 print for each SCF-MI cycle MOs, overlap
                        between the MOs, CPU times.
                    = 2 print some extra informations in secular
                        systems solution.

          MSHIFT    = debugging option that permits to shift all
                      the memory pointer of the SCF-MI section
                      of code of the quantity MSHIFT (default is 0).


             "Modification of Roothan Equations to Exclude BSSE
                 from Molecular Interaction Calculations"
              E. Gianinetti, M. Raimondi, E. Tornaghi
              Int. J. Quantum Chem. 60, 157 (1996)

              A. Famulari, E. Gianinetti, M. Raimondi, and M. Sironi
              Int. J. Quantum Chem. (1997), submitted.



$DFT group (relevant if SCFTYP=RHF,UHF,ROHF)

               This group permits the use of various empirical one-
          electron operators instead of the correct many electron
          Hamiltonian.  The implementation is based on the use of
          the resolution of the identity to simplify integrals so
          that they may be analytically evaluated, instead of the
          use of grid quadratures.  The grid free DFT computations
          in their present form have various numerical errors.

          DFTTYP = NONE     means ab initio computation (default)
                            exchange functionals:
                 = XALPHA   X-Alpha exchange (alpha=0.7)
                 = SLATER   Slater exchange (alpha=2/3)
                 = LOCAL    a synonym for SLATER
                 = LSDA     a synonym for SLATER
                 = BECKE    Becke's 1988 exchange
                 = DEPRISTO Depristo/Kress exchange
                 = CAMA     Handy et al's mods to Becke exchange
                 = HALF     50-50 mix of Becke and HF exchange
                            correlation functionals:
                 = VWN      Vosko/Wilke/Nusair correlation, formula 5
                 = PWLOC    Perdew/Wang local correlation
                 = LYP      Lee/Yang/Parr correlation
                            exchange/correlation functionals:
                 = BVWN     Becke exchange + VWN correlation
                 = BLYP     Becke exchange + LYP correlation
                 = BPWLOC   Becke exchange + Perdew/Wang correlation
                 = B3LYP    hybridized HF/Becke/LYP using VWN formula 5
                 = CAMB     CAMA exchange + Cambridge correlation
                 = XVWN     Xalpha exchange + VWN formula 5 correlation
                 = XPWLOC   Xalpha exchange + Perdew/Wang correlation
                 = SVWN     Slater exchange + VWN correlation
                 = SPWLOC   Slater exchange + PWLOC correlation
                 = WIGNER   Wigner exchange + correlation
                 = WS       Wigner scaled exchange + correlation
                 = WIGEXP   Wigner exponential exchange + correlation

          AUXFUN = AUX0  uses no auxiliary basis set for resolution
                         of the identity, limiting accuracy.
                 = AUX3  uses the 3rd generation of RI basis sets,
                         These are available for the elements H to
                         Ar, but have been carefully considered for
                         H-Ne only.  (DEFAULT)

          THREE  = a flag to use a resolution of the identity to
                   turn four center overlap integrals into three
                   center integrals.  This can be used only if
                   no auxiliary basis is employed. (default=.FALSE.)


            Do not use this input group without reading about the
           numerical limitations of the grid free code in REFS.DOC



$MP2 group (relevant to SCFTYP=RHF,UHF,ROHF if MPLEVL=2)

               Controls 2nd order Moller-Plesset perturbation runs,
          if requested by MPLEVL in $CONTRL.  See also the DIRSCF
          keyword in $SCF to select direct MP2.  MP2 is implemented
          for RHF, high spin ROHF, or UHF wavefunctions.  Analytic
          gradients and the first order correction to the wave-
          function (i.e. properties) are available only for RHF.
          The $MP2 group is usually not given.  See also $MCQDPT.

          NCORE = n  Omits the first n occupied orbitals from the
                     calculation.  The default for n is the number
                     of chemical core orbitals.

          MP2PRP=    a flag to turn on property computation for
                     RHF MP2 jobs with RUNTYP=ENERGY.  This is
                     appreciably more expensive than just evaluating
                     the 2nd order energy correction alone, so the
                     default is .FALSE.  Properties are always
                     computed during gradient runs, when they are
                     an almost free byproduct. (default=.FALSE.)

                     This parameter applies only to the serial MP2
                     program.  To see properties using the parallel
                     DDI code, use RUNTYP=GRADIENT.

          LMOMP2=    a flag to analyze the closed shell MP2 energy
                     in terms of localized orbitals.  Any type of
                     localized orbital may be used.  This option
                     is implemented only for RHF, and its selection
                     forces use of the METHOD=3 transformation.
                     The default is .FALSE.

          OSPT=      selects open shell spin-restricted perturbation.
                     This parameter applies only when SCFTYP=ROHF.
                     Please see the 'further information' section for
                     more information about this choice.
              = ZAPT picks Z-averaged perturbation theory. (default)
              = RMP  picks RMP (aka ROHF-MBPT) perturbation theory.

          CUTOFF=    transformed integral retention threshold, the
                     default is 1.0d-9.

          The last 3 input variables apply only UHF+MP2 or ROHF+MP2
          using OSPT=RMP, or to runs on one compute node only.

          NWORD =    controls memory usage.  The default uses all
                     available memory.  (default=0)

          METHOD= n  selects transformation method, 2 being the
                     segmented transformation, and 3 being a more
                     conventional two phase bin sort implementation.
                     3 requires more disk, but less memory.  The
                     default is to attempt method 2 first, and
                     method 3 second.

          AOINTS=    defines AO integral storage during conventional
                     integral transformations, during parallel runs.
                  DUP stores duplicated AO lists on each node, and
                     is the default for parallel computers with slow
                     interprocessor communication, e.g. ethernet.
                  DIST distributes the AO integral file across all
                     nodes, and is the default for parallel
                     computers with high speed communications.



$GUESS group (optional, relevant for all SCFTYP's)

              This group controls the selection of initial molecular

          GUESS = Selects type of initial orbital guess.
                = HUCKEL   Carry out an extended Huckel calculation
                           using a Huzinaga MINI basis set, and
                           project this onto the current basis.
                           This is implemented for atoms up to Rn,
                           and will work for any all electron or
                           ECP basis set.  (default for most runs)
                = HCORE    Diagonalize the one electron Hamiltonian
                           to obtain the initial guess orbitals.
                           This method is applicable to any basis
                           set, but does not work as well as the
                           HUCKEL guess.
                = MOREAD   Read in formatted vectors punched by an
                           earlier run.  This requires a $VEC group,
                           and you MUST pay attention to NORB below.
                = MOSAVED  (default for restarts)  The initial
                           orbitals are read from the DICTNRY file
                           of the earlier run.
                = SKIP     Bypass initial orbital selection.  The
                           initial orbitals and density matrix are
                           assumed to be in the DICTNRY file.  Mostly
                           used for RUNTYP=HESSIAN when the hessian
                           is being read in from the input.

              All GUESS types except 'SKIP' permit reordering of the
          orbitals, carry out an orthonormalization of the orbitals,
          and generate the correct initial density matrix, for RHF,
          UHF, ROHF, and GVB, but note that correct computation of
          the GVB density requires also CICOEF in $SCF.  The density
          matrix cannot be generated from the orbitals alone for MP2,
          CI, or MCSCF, so property evaluation for these should be
          RUNTYP=ENERGY rather than RUNTYP=PROP using GUESS=MOREAD.

          PRTMO = a flag to control printing of the initial guess.

          PUNMO = a flag to control punching of the initial guess.

          MIX    = rotate the alpha and beta HOMO and LUMO orbitals
                   so as to generate inequivalent alpha and beta
                   orbital spaces.  This pertains to UHF singlets
                   only.  This may require use of NOSYM=1 in $CONTRL
                   depending on your situation.  (default=.FALSE.)

          NORB   = The number of orbitals to be read in the $VEC
                   group.  This applies only to GUESS=MOREAD.

          For -RHF-, -UHF-, -ROHF-, and -GVB-, NORB defaults to the
          number of occupied orbitals.  NORB must be given for -CI-
          and -MCSCF-.  For -UHF-, if NORB is not given, only the
          occupied alpha and beta orbitals should be given, back to
          back.  Otherwise, both alpha and beta orbitals must
          consist of NORB vectors.
          NORB may be larger than the number of occupied MOs, if you
          wish to read in the virtual orbitals.  If NORB is less
          than the number of atomic orbitals, the remaining orbitals
          are generated as the orthogonal complement to those read.

          NORDER = Orbital reordering switch.
                 = 0  No reordering (default)
                 = 1  Reorder according to IORDER and JORDER.

          IORDER = Reordering instructions.
                   Input to this array gives the new molecular
                   orbital order.  For example, IORDER(3)=4,3 will
                   interchange orbitals 3 and 4, while leaving the
                   other MOs in the original order.  This parameter
                   applies to all orbitals (alpha and beta) except
                   for -UHF-, where it only affects the alpha MOs.
                   (default is IORDER(i)=i )

          JORDER = Reordering instructions.
                   Same as IORDER, but for the beta MOs of -UHF-.

         INSORB = the first INSORB orbitals specified in the $VEC
                  group will be inserted into the Huckel guess,
                  making the guess a hybrid of HUCKEL/MOREAD.  This
                  keyword is meaningful only when GUESS=HUCKEL, and
                  it is useful mainly for QM/MM runs where some
                  orbitals (buffer) are frozen and need to be
                  transferred to the initial guess vector set,
                  see $MOFRZ.  (default=0)

            * * * the next are 3 ways to clean up orbitals * * *

          PURIFY = flag to symmetrize starting orbitals.
                   This is the most soundly based of the possible
                   procedures.  (default=.FALSE.)

          TOLZ   = level below which MO coefficients will be set
                   to zero.  (default=1.0E-7)

          TOLE   = level at which MO coefficients will be equated.
                   This is a relative level, coefficients are set
                   equal if one agrees in magnitude to TOLE times
                   the other.  (default=5.0E-5)

          SYMDEN = project the totally symmetric from the density.
                   Maybe useful if the HUCKEL or HCORE give orbitals
                   of inexact symmetry.  Since the density matrix is
                   not idempotent, this can generate a non-variational
                   energy on the first iteration.  For the same
                   reason, this should never be used with orbitals
                   of MOREAD quality.  (default=.FALSE.)



$VEC group (optional, relevant for all SCFTYP's; required if GUESS=MOREAD)

                This group consists of formatted vectors, as written
          onto file PUNCH in a previous run.  It is considered good
          form to retain the titling comment cards punched before
          the $VEC card, as a reminder to yourself of the origin of
          the orbitals.

                For Morokuma decompositions, the names of this group
          are $VEC1, $VEC2, ... for each monomer, computed in the
          identical orientation as the supermolecule.  For transition
          moment or spin-orbit coupling runs, orbitals for states
          one and possibly two are $VEC1 and $VEC2.



$MOFRZ group (optional, relevant for RHF, ROHF, GVB)

              This group controls freezing the molecular orbitals
          of your choice during the SCF procedure.  If you choose
          this option, select DIIS in $SCF since SOSCF will not
          converge as well.  GUESS=MOREAD is required in $GUESS.

          FRZ   = flag which triggers MO freezing. (default=.FALSE.)

          IFRZ  = an array of MOs in the input $VEC set which are
                  to be frozen.  There is no default for this.



$STATPT group (optional, for RUNTYP=OPTIMIZE or SADPOINT)

              This group controls the search for stationary points.
          Note that NZVAR in $CONTRL determines if the geometry
          search is conducted in Cartesian or internal coordinates.

          METHOD = optimization algorithm selection.  Pick from

                   NR   Straight Newton-Raphson iterate. This will
                        attempt to locate the nearest stationary
                        point, which may be of any order. There
                        is no steplength control. RUNTYP can be
                        either OPTIMIZE or SADPOINT

                   RFO  Rational Function Optimization. This is
                        one of the augmented Hessian techniques
                        where the shift parameter(s) is(are) chosen
                        by a rational function approximation to
                        the PES. For SADPOINT searches it involves
                        two shift parameters. If the calculated
                        stepsize is larger than DXMAX the step is
                        simply scaled down to size.

                   QA   Quadratic Approximation. This is another
                        version of an augmented Hessian technique
                        where the shift parameter is chosen such
                        that the steplength is equal to DXMAX.
                        It is completely equivalent to the TRIM
                        method. (default)

                   SCHLEGEL The quasi-NR optimizer by Schlegel.

                   CONOPT, CONstrained OPTimization. An algorithm
                        which can be used for locating TSs.
                        The starting geometry MUST be a minimum!
                        The algorithm tries to push the geometry
                        uphill along a chosen Hessian mode (IFOLOW)
                        by a series of optimizations on hyperspheres
                        of increasingly larger radii.
                        Note that there currently are no restart
                        capabilitites for this method, not even

          OPTTOL = gradient convergence tolerance, in Hartree/Bohr.
                   Convergence of a geometry search requires the
                   largest component of the gradient to be less
                   than OPTTOL, and the root mean square gradient
                   less than 1/3 of OPTTOL.  (default=0.0001)

          NSTEP  = maximum number of steps to take.  Restart data
                   is punched if NSTEP is exceeded. (default=20)

                --- the next four control the step size ---

          DXMAX  = initial trust radius of the step, in Bohr.
                   For METHOD=RFO, QA, or SCHLEGEL, steps will
                   be scaled down to this value, if necessary.
                   (default=0.3 for OPTIMIZE and 0.2 for SADPOINT)
                   For METHOD=NR, DXMAX is inoperative.
                   For METHOD=CONOPT, DXMAX is the step along the
                   previous two points to increment the hypersphere
                   radius between constrained optimizations.

              the next three apply only to METHOD=RFO or QA:

          TRUPD  = a flag to allow the trust radius to change as
                   the geometry search proceeds.  (default=.TRUE.)

          TRMAX  = maximum permissible value of the trust radius.
                   (default=0.5 for OPTIMIZE and 0.3 for SADPOINT)

          TRMIN  = minimum permissible value of the trust radius.

               --- the next three control mode following ---

          IFOLOW = Mode selection switch, for RUNTYP=SADPOINT.
                   For METHOD=RFO or QA, the mode along which the
                   energy is maximized, other modes are minimized.
                   Usually refered to as "eigenvector following".
                   For METHOD=SCHLEGEL, the mode whose eigenvalue
                   is (or will be made) negative.  All other
                   curvatures will be made positive.
                   For METHOD=CONOPT, the mode along which the
                   geometry is initially perturbed from the minima.
                   (default is 1)
                   In Cartesian coordinates, this variable doesn't
                   count the six translation and rotation degrees.
                   Note that the "modes" aren't from mass-weighting.

          STPT   = flag to indicate whether the initial geometry
                   is considered a stationary point. If .true.
                   the initial geometry will be perturbed by
                   a step along the IFOLOW normal mode with
                   stepsize STSTEP. (default=.false.)
                   The positive direction is taken as the one where
                   the largest component of the Hessian mode is
                   positive. If there are more than one largest
                   component (symmetry), the first is taken as
                   Note that STPT=.TRUE. has little meaning with
                   HESS=GUESS as there will be many degenerate

          STSTEP = Stepsize for jumping off a stationary point.
                   Using values of 0.05 or more may work better.

          IFREEZ = array of coordinates to freeze.  These may be
                   internal or Cartesian coordinates.  For example,
                   IFREEZ(1)=1,3 freezes the two bond lengths in
                   the $ZMAT example, while optimizing the angle.
                   If NZVAR=0, so that this value applies to the
                   Cartesian coordinates instead, the input of
                   IFREEZ(1)=4,7 means to freeze the x coordinates
                   if the 2nd and 3rd atoms in the molecule.

                   See also IFZMAT and FVALUE in $ZMAT, and IFCART
                   below, as IFREEZ does not apply to DLC internals.

                   In a numerical Hessian run, IFREEZ specifies
                   Cartesian displacements to be skipped for a
                   Partial Hessian Analysis.  For more information:
                   J.D. Head, Int. J. Quantum Chem. 65, 827, 1997
                   H. Li, J.H. Jensen, manuscript in preparation.

          IFCART = array of Cartesian coordinates to freeze during
                   a geometry optimization using delocalized internal

           --- The next two control the hessian matrix quality ---

          HESS   = selects the initial hessian matrix.
                 = GUESS chooses a positive definite diagonal
                         hessian. (default for RUNTYP=OPTIMIZE)
                 = READ  causes the hessian to be read from a $HESS
                         group. (default for RUNTYP=SADPOINT)
                 = RDAB  reads only the ab initio part of the
                         hessian, and approximates the effective
                         fragment blocks.
                 = RDALL reads the full hessian, then converts
                         any fragment blocks to 6x6 T+R shape.
                         (this option is seldom used).
                 = CALC  causes the hessian to be computed, see
                         the $FORCE group.

          IHREP  = the number of steps before the hessian is
                   recomputed.  If given as 0, the hessian will
                   be computed only at the initial geometry if
                   you choose HESS=CALC, and never again.  If
                   nonzero, the hessian is recalculated every
                   IHREP steps, with the update formula used on
                   other steps.  (default=0)

          HSSEND = a flag to control automatic hessian evaluation
                   at the end of a successful geometry search.

             --- the next two control the amount of output ---
              Let 0 mean the initial geometry, L mean the last
              geometry, and all mean every geometry.
              Let INTR mean the internuclear distance matrix.
              Let HESS mean the approximation to the hessian.
              Note that a directly calculated hessian matrix
              will always be punched, NPUN refers only to the
              updated hessians used by the quasi-Newton step.

          NPRT   =  1  Print INTR at all, orbitals at all
                    0  Print INTR at all, orbitals at 0+L (default)
                   -1  Print INTR at all, orbitals never
                   -2  Print INTR at 0+L, orbitals never

          NPUN   =  3  Punch all orbitals and HESS at all
                    2  Punch all orbitals at all
                    1  same as 0, plus punch HESS at all
                    0  Punch all orbitals at 0+L, otherwise only
                       occupied orbitals (default)
                   -1  Punch occ orbitals at 0+L only
                   -2  Never punch orbitals

           ---- the following parameters are quite specialized ----

          PURIFY = a flag to help eliminate the rotational and
                   translational degrees of freedom from the
                   initial hessian (and possibly initial gradient).
                   This is much like the variable of the same name
                   in $FORCE, and will be relevant only if internal
                   coordinates are in use.  (default=.FALSE.)

          PROJCT = a flag to eliminate translation and rotational
                   degrees of freedom from Cartesian optimizations.
                   The default is .TRUE. since this normally will
                   reduce the number of steps, except that this
                   variable is set false when POSITION=FIXED is
                   used during EFP runs.

          ITBMAT = number of micro-iterations used to compute the
                   step in Cartesians which corresponds to the
                   desired step in internals.  The default is 5.

          UPHESS = SKIP     do not update Hessian (not recommended)
                   BFGS     default for OPTIMIZE using RFO or QA
                   POWELL   default for OPTIMIZE using NR or CONOPT
                   POWELL   default for SADPOINT
                   MSP      mixed Murtagh-Sargent/Powell update
                   SCHLEGEL only choice for METHOD=SCHLEGEL

          MOVIE  = a flag to create a series of structural data
                   which can be show as a movie by the MacIntosh
                   program Chem3D.  The data is written to the
                   file IRCDATA.  (default=.FALSE.)

           ---- NNEG, RMIN, RMAX, RLIM apply only to SCHLEGEL ----

          NNEG   = The number of negative eigenvalues the force
                   constant matrix should have. If necessary the
                   smallest eigenvalues will be reversed. The
                   default is 0 for RUNTYP=OPTIMIZE, and 1 for

          RMIN   = Minimum distance threshold. Points whose root
                   mean square distance from the current point is
                   less than RMIN are discarded. (default=0.0015)

          RMAX   = Maximum distance threshold. Points whose root
                   mean square distance from the current point is
                   greater than RMAX are discarded. (default=0.1)

          RLIM   = Linear dependence threshold. Vectors from the
                   current point to the previous points must not
                   be colinear.  (default=0.07)


$TRUDGE group (optional, required for RUNTYP=TRUDGE)

              This group defines the parameters for a non-gradient
          optimization of exponents or the geometry.  The TRUDGE
          package is a modified version of the same code from Michel
          Dupuis' HONDO 7.0 system, origially written by H.F.King.
          Presently the program allows for the optimization of 10

              Exponent optimization works only for uncontracted
          primitives, without enforcing any constraints.  Two
          non-symmetry equivalent H atoms would have their p
          function exponents optimized separately, and so would two
          symmetry equivalent atoms!  A clear case of GIGO.

              Geometry optimization works only in HINT internal
          coordinates (see $CONTRL and $DATA groups).  The total
          energy of all types of SCF wavefunctions can be optimized,
          although this would be extremely stupid as gradient
          methods are far more efficient.  The main utility is for
          open shell MP2 or CI geometry optimizations, which may
          not be done in any other way with GAMESS.  If your run
          requires NOSYM=1 in $CONTRL, you must be sure to use only
          C1 symmetry in the $DATA group.

          OPTMIZ = a flag to select optimization of either geometry
                   or exponents of primitive gaussian functions.
                 = BASIS    for basis set optimization.
                 = GEOMETRY for geometry optimization (default).
                   This means minima search only, there is no saddle
                   point capability.

          NPAR   = number of parameters to be optimized.

          IEX    = defines the parameters to be optimized.

                   If OPTMIZ=BASIS, IEX declares the serial number
              of the Gaussian primitives for which the exponents
              will be optimized.

                   If OPTMIZ=GEOMETRY, IEX define the pointers to
              the HINT internal coordinates which will be optimized.
              (Note that not all internal coordinates have to be
              optimized.) The pointers to the internal coordinates
              are defined as:  (the number of atom on the input
              list)*10 + (the number of internal coordinate for that
              atom).  For each atom, the HINT internal coordinates
              are numbered as 1, 2, and 3 for BOND, ALPHA, and BETA,

          P  =  Defines the initial values of the parameters to be
                optimized.  You can use this to reset values given
                in $DATA.  If omitted, the $DATA values are used.
                If given here, geometric data must be in Angstroms
                and degrees.

          A complete example is a TCSCF multireference 6-31G
          geometry optimization for methylene,

                   COORD=HINT $END
           $BASIS  GBASIS=N31 NGAUSS=6 $END
          Methylene TCSCF+CISD geometry optimization
          Cnv 2

          C    6.     LC  0.00  0.0  0.00  -  O  K
          H    1.    PCC  1.00  53.  0.00  +  O  K  I
           $SCF    NCO=3 NPAIR=1 $END
                   IEX(1)=21,22   P(1)=1.08 $END
                   NEXT=-1 $END

          using GVB-PP(1), or TCSCF orbitals in the CI.  The starting
          bond length is reset to 1.09, while the initial angle will
          be 106 (twice 53).  Result after 17 steps is R=1.1283056,
          half-angle=51.83377, with a CI energy of -38.9407538472

              Note that you may optimize the geometry for an excited
          CI state, just specify
                    $GUGDIA   NSTATE=5  $END
                    $GUGDM    IROOT=3   $END
          to find the equilibrium geometry of the third state (of
          five total states) of the symmetry implied by your $CIDRT.



$TRURST group (optional, relevant for RUNTYP=TRUDGE)

                This  group  specifies restart parameters for TRUDGE
          runs and accuracy thresholds.

          KSTART indicates the conjugate gradient direction in which
          the optimization will proceed. ( default = -1 )
               -1 .... indicates that this is a non-restart run.
                0 .... corresponds to a restart run.

          FNOISE accuracy of function values.
          Variation smaller than FNOISE are not considered to be
          significant (Def. 0.0005)

          TOLF accuracy required of the function (Def. 0.001)

          TOLR accuracy required of conjugate directions (Def. 0.05)

              For geometry optimization, the values which give
          better results (closer to the ones obtained with gradient
          methods) are:  TOLF=0.0001, TOLR=0.001, FNOISE=0.00001



$FORCE group (optional, relevant for RUNTYP=HESSIAN,OPTIMIZE,SADPOINT)

              This group controls the computation of the hessian
          matrix (the energy second derivative tensor, also known
          as the force constant matrix), and an optional harmonic
          vibrational analysis.  This can be a very time consuming
          calculation.  However, given the force constant matrix,
          the vibrational analysis for an isotopically substituted
          molecule is very cheap.  Related input is HESS= in
          $STATPT, and the $MASS, $HESS, $GRAD, $DIPDR, $VIB groups.

          METHOD = chooses the computational method.
                 = ANALYTIC is implemented only for SCFTYPs RHF,
                            ROHF, and GVB (when NPAIR is 0 or 1).
                            This is the default for these cases.
                 = NUMERIC  is the default for all other cases:
                            UHF, MCSCF, and all MP2 or CI runs.

          RDHESS = a flag to read the hessian from a $HESS group,
                   rather than computing it.  This variable pertains
                   only to RUNTYP=HESSIAN.  See also HESS= in the
                   $STATPT group.  (default is .FALSE.)

          PURIFY = controls cleanup
                   Given a $ZMAT, the hessian and dipole derivative
                   tensor can be "purified" by transforming from
                   Cartesians to internals and back to Cartesians.
                   This effectively zeros the frequencies of the
                   translation and rotation "modes", along with
                   their IR intensities.  The purified quantities
                   are punched out.  Purification does change the
                   Hessian slightly, frequencies at a stationary
                   point can change by a wave number or so.  The
                   change is bigger at non-stationary points.
                   (default=.FALSE. if $ZMAT is given)

          PRTIFC = prints the internal coordinate force constants.
                   You MUST have defined a $ZMAT group to use this.

            --- the next four apply only to METHOD=NUMERIC ----

          NVIB   =    Number of displacements in each Cartesian
                      direction for force field computation.
                 = 1  Move one VIBSIZ unit in each positive
                      Cartesian direction.  This requires 3N+1
                      evaluations of the wavefunction, energy, and
                      gradient, where N is the number of SYMMETRY
                      UNIQUE atoms given in $DATA.  (default)
                 = 2  Move one VIBSIZ unit in the positive direction
                      and one VIBSIZ unit in the negative direction.
                      This requires 6N+1 evaluations of the
                      wavefunction and gradient, and gives a small
                      improvement in accuracy.  In particular, the
                      frequencies will change from NVIB=1 results by
                      no more than 10-100 wavenumbers, and usually
                      much less.  However, the normal modes will be
                      more nearly symmetry adapted, and the residual
                      rotational and translational "frequencies"
                      will be much closer to zero.

          VIBSIZ =    Displacement size (in Bohrs). Default=0.01

                 Let 0 mean the Vib0 geometry, and
                 D mean all the displaced geometries

          NPRT   = 1  Print orbitals at 0 and D
                 = 0  Print orbitals at 0 only (default)

          NPUN   = 2  Punch all orbitals at 0 and D
                 = 1  Punch all orbitals at 0 and occupied orbs at D
                 = 0  Punch all orbitals at 0 only (default)

            ----- the rest control normal coordinate analysis ----

          VIBANL = flag to activate vibrational analysis.
                   (the default is .TRUE. for RUNTYP=HESSIAN, and
                   otherwise is .FALSE.)

          SCLFAC = scale factor for vibrational frequencies, used
                   in calculating the zero point vibrational energy.
                   Some workers correct for the usual overestimate
                   in SCF frequencies by a factor 0.89.  ZPE or other
                   methods might employ other factors, see A.P.Scott,
                   L.Radom  J.Phys.Chem.  100, 16502-16513 (1996).
                   The output always prints unscaled frequencies, so
                   this value is used only during the thermochemical
                   analysis.  (Default is 1.0)

          TEMP   = an array of up to ten temperatures at which the
                   thermochemistry should be printed out.  The
                   default is a single temperature, 298.15 K.  To
                   use absolute zero, input 0.001 degrees.

          FREQ   = an array of vibrational frequencies.  If the
                   frequencies are given here, the hessian matrix
                   is not computed or read.  You enter any imaginary
                   frequencies as negative numbers, omit the
                   zero frequencies corresponding to translation
                   and rotation, and enter all true vibrational
                   frequencies.  Thermodynamic properties will be
                   printed, nothing else is done by the run.

          PRTSCN = flag to print contribution of each vibrational
                   mode to the entropy.  (Default is .FALSE.)

          DECOMP = activates internal coordinate analysis.
                   Vibrational frequencies will be decomposed into
                   "intrinsic frequencies", by the method of
                   J.A.Boatz and M.S.Gordon, J.Phys.Chem., 93,
                   1819-1826(1989).  If set .TRUE., the $ZMAT group
                   may define more than 3N-6 (3N-5) coordinates.

          PROJCT = controls the projection of the hessian matrix.
                   The projection technique is described by
                   W.H.Miller, N.C.Handy, J.E.Adams in J. Chem.
                   Phys. 1980, 72, 99-112.  At stationary points,
                   the projection simply eliminates rotational and
                   translational contaminants.  At points with
                   non-zero gradients, the projection also ensures
                   that one of the vibrational modes will point
                   along the gradient, so that there are a total of
                   7 zero frequencies.  The other 3N-7 modes are
                   constrained to be orthogonal to the gradient.
                   Because the projection has such a large effect on
                   the hessian, the hessian punched is the one
                   BEFORE projection.  For the same reason, the
                   default is .FALSE. to skip the projection, which
                   is mainly of interest in dynamical calculations.


          There is a set of programs for the calculation of kinetic
          or equilibrium isotope effects from the group of Piotr
          Paneth at the University of Lodz.  This ISOEFF package will
          accept data computed by GAMESS, and can be downloaded at



$CPHF group (relevant for analytic RUNTYP=HESSIAN)

              This group controls the solution of the response
          equations, also known as coupled Hartree-Fock.

          POLAR = a flag to request computation of the static
                  polarizability, alpha.  Because this property
                  needs 3 additional response vectors, beyond those
                  needed for the hessian, the default is to skip the
                  property.  (default = .FALSE.)

          NWORD = controls memory usage for this step.  The default
                  uses all available memory.  (default=0)



$HESS group (relevant for RUNTYP=HESSIAN if RDHESS=.TRUE., relevant for RUNTYP=IRC if FREQ,CMODE not given & relevant for RUNTYP=OPTIMIZE,SADPOINT if HESS=READ)

              Formatted force constant matrix (FCM), i.e. hessian
          matrix.  This data is punched out by a RUNTYP=HESSIAN job,
          in the correct format for subsequent runs.  The first card
          in the group must be a title card.

              A $HESS group is always punched in Cartesians.  It
          will be transformed into internal coordinate space if a
          geometry search uses internals.  It will be mass weighted
          (according to $MASS) for IRC and frequency runs.

              The initial FCM is updated during the course of a
          geometry optimization or saddle point search, and will be
          punched if a run exhausts its time limit.  This allows
          restarts where the job leaves off.  You may want to read
          this FCM back into the program for your restart, or you
          may prefer to regenerate a new initial hessian.  In any
          case, this updated hessian is absolutely not suitable for
          frequency prediction!



$GRAD group (relevant for RUNTYP=OPTIMIZE or SADPOINT & relevant for RUNTYP=HESSIAN when RDHESS=.TRUE.)

              Formatted gradient vector at the $DATA geometry.  This
          data is read in the same format it was punched out.

              For RUNTYP=HESSIAN, this information is used to
          determine if you are at a stationary point, and possibly
          for projection.  If omitted, the program pretends the
          gradient is zero, and otherwise proceeds normally.

              For geometry searches, this information (if known) can
          be read into the program so that the first step can be
          taken instantly.



$DIPDR group (relevant for RUNTYP=HESSIAN if RDHESS=.T.)

          Formatted dipole derivative tensor, punched in a previous
          RUNTYP=HESSIAN job.  If this group is omitted, then a
          vibrational analysis will be unable to predict the IR
          intensities, but the run can otherwise proceed.




              Formatted card image -restart- data.  This data is
          read in the format it was punched by a previous HESSIAN
          job to the file IRCDATA.  Just add a " $END" card, and if
          the final gradient was punched as zero, delete the last
          set of data.  Normally, IREST in $CONTRL will NOT be used
          in conjunction with a HESSIAN restart.  The mere presence
          of this deck triggers the restart from cards.  This deck
          can also be used to turn a single point differencing run
          into double differencing, as well as recovering from time
          limits, or other bombouts.



$MASS group (relevant for RUNTYP=HESSIAN, IRC, or DRC)

              This group permits isotopic substitution during the
          computation of mass weighted Cartesian coordinates.  Of
          course, the masses affect the frequencies and normal modes
          of vibration.

          AMASS = An array giving the atomic masses, in amu. The
                  default is to use the mass of the most abundant
                  isotope.  Masses through element 104 are stored.

          example - $MASS AMASS(3)=2.0140 $END
          will make the third atom in the molecule a deuterium.



$IRC group (relevant for RUNTYP=IRC)

              This group governs the location of the intrinsic
          reaction coordinate, a steepest descent path in mass
          weighted corrdinates, that connects the saddle point to
          reactants and products.

          ----- there are five integration methods chosen by PACE.

          PACE = GS2    selects the Gonzalez-Schlegel second order
                        method.  This is the default method.
                        Related input is:

            GCUT   cutoff for the norm of the mass-weighted gradient
                   tangent (the default is chosen in the range from
                   0.00005 to 0.00020, depending on the value for
                   STRIDE chosen below.
            RCUT   cutoff for Cartesian RMS displacement vector.
                   (the default is chosen in the range 0.0005 to
                   0.0020 Bohr, depending on the value for STRIDE)
            ACUT   maximum angle from end points for linear
                   interpolation (default=5 degrees)
            MXOPT  maximum number of contrained optimization steps
                   for each IRC point (default=20)
            IHUPD  is the hessian update formula.  1 means Powell,
                   2 means BFGS (default=2)
            GA     is a gradient from the previous IRC point, and is
                   used when restarting.
            OPTTOL is a gradient cutoff used to determine if the IRC
                   is approaching a minimum.  It has the same meaning
                   as the variable in $STATPT.  (default=0.0001)

          PACE = LINEAR selects linear gradient following (Euler's
                        method).  Related input is:

            STABLZ switches on Ishida/Morokuma/Komornicki reaction
                   path stabilization.  The default is .TRUE.
            DELTA  initial step size along the unit bisector, if
                   STABLZ is on.  Default=0.025 Bohr.
            ELBOW  is the collinearity threshold above which the
                   stabilization is skipped.  If the mass weighted
                   gradients at QB and QC are almost collinear, the
                   reaction path is deemed to be curving very little,
                   and stabilization isn't needed.  The default is
                   175.0 degrees.  To always perform stabilization,
                   input 180.0.
            READQB,EB,GBNORM,GB are energy and gradient data
                   already known at the current IRC point.  If it
                   happens that a run with STABLZ on decides to skip
                   stabilization because of ELBOW, this data will be
                   punched to speed the restart.

          PACE = QUAD   selects quadratic gradient following.
                        Related input is:

            SAB    distance to previous point on the IRC.
            GA     gradient vector at that historical point.

          PACE = AMPC4  selects the fourth order Adams-Moulton
                        variable step predictor-corrector.
                        Related input is:

            GA0,GA1,GA2 which are gradients at previous points.

          PACE = RK4    selects the 4th order Runge-Kutta variable
                        step method.  There is no related input.

          ----- The next two are used by all PACE choices -----

          STRIDE = Determines how far apart points on the reaction
                   path will be.  STRIDE is used to calculate the
                   step taken, according to the PACE you choose.
                   The default is good for the GS2 method, which is
                   very robust.  Other methods should request much
                   smaller step sizes, such as 0.10 or even 0.05.
                   (default = 0.30 sqrt(amu)-Bohr)
          NPOINT = The number of IRC points to be located in this
                   run. The default is to find only the next point.
                   (default = 1)

          ----- The next two let you choose your output volume -----

              Let F mean the first IRC point found in this run,
              and L mean the final IRC point of this run.
              Let INTR mean the internuclear distance matrix.

          NPRT   =  1  Print INTR at all, orbitals at all IRC points
                    0  Print INTR at all, orbitals at F+L (default)
                   -1  Print INTR at all, orbitals never
                   -2  Print INTR at F+L, orbitals never

          NPUN   =  1  Punch all orbitals at all IRC points
                    0  Punch all orbitals at F+L, only occupied
                       orbitals at IRC points between (default)
                   -1  Punch all orbitals at F+L only
                   -2  Never punch orbitals

          ----- The next two tally the reaction path results.  The
                defaults are appropriate for starting from a saddle
                point, restart values are automatically punched out.

          NEXTPT = The number of the next point to be computed.
          STOTAL = Total distance along the reaction path to next
                   IRC point, in mass weighted Cartesian space.

          ----- The following controls jumping off the saddle point.
                If you give a $HESS group, FREQ and CMODE will be
                generated automatically.

          SADDLE = A logical variable telling if the coordinates
                   given in the $DATA deck are at a saddle point
                   (.TRUE.) or some other point lying on the IRC
                   (.FALSE.).  If SADDLE is true, either a $HESS
                   group or else FREQ and CMODE must be given.
                   (default = .FALSE.)  Related input is:

          TSENGY = A logical variable controlling whether the energy
                   and wavefunction are evaluated at the transition
                   state coordinates given in $DATA.  Since you
                   already know the energy from the transition state
                   search and force field runs, the default is .F.
          FORWRD = A logical variable controlling the direction to
                   proceed away from a saddle point. The forward
                   direction is defined as the direction in which
                   the largest magnitude component of the imaginary
                   normal mode is positive. (default =.TRUE.)
          EVIB   = Desired decrease in energy when following the
                   imaginary normal mode away from a saddle point.
                   (default=0.0005 Hartree)
          FREQ   = The magnitude of the imaginary frequency, given
                   in cm**-1.
          CMODE  = An array of the components of the normal mode
                   whose frequency is imaginary, in Cartesian
                   coordinates.  Be careful with the signs!

             You must give FREQ and CMODE if you don't give a $HESS
             group, when SADDLE=.TRUE.  The option of giving these
             two variables instead of a $HESS does not apply to the
             GS2 method, which must have a hessian input, even for
             restarts.  Note also that EVIB is ignored by GS2 runs.



$VSCF group (optional, relevant to RUNTYP=VSCF)

              This group governs the computation of frequencies
          including anharmonic effects.  Besides the values shown
          below, the input file must contain a $HESS group and
          perhaps a $DIPDR group, to start with previously obtained
          harmonic vibrational information.  Energies are sampled
          along the directions of harmonic normal modes, and along
          pairs of harmonic normal modes, after which vibrational
          nuclear wavefunctions are obtained at an SCF-like level,
          termed VSCF, using product nuclear wavefunctions.  An
          MP2-like correction to the vibrational energy, termed
          correlation corrected (cc-VSCF), is also obtained.  By
          default, the dipole is computed at every grid point to
          give improved IR intensity values.  See also the restart
          group $VIBSCF.

          NGRID  = number of grid points to be computed along each
                   harmonic normal mode, and if NCOUP=2, along each
                   pair of modes.  Reasonable values are 8 or 16,
                   with 16 considered significantly more accurate.

          NCOUP  = the order of mode couplings included.
                 = 1 computes 1-D grids along each harmonic mode
                 = 2 adds additionally, 2-D grids along each pair
                     of normal modes. (default)

             The total number of energy and dipole evaluations
             for NCOUP=2 is M*NGRID + M*(M-1)/2*NGRID**2, where
             M is the number of normal modes: M = 3N-6 or 3N-5.

          IEXC   = 1 obtain fundamental frequencies (default)
                 = 2 instead, obtain first overtones
                 = 3 instead, obtain second overtones

                     IEXC higher than 1 may be speedily obtained
                     using the next parameter to restart with a
                     completed $VIBSCF group.

          READV  = flag to indicate restart data $VIBSCF should be
                   read in to resume an interrupted calculation, or
                   to obtain overtones in follow-on runs.
                   (default is .FALSE.)

          The next two relate to simplified intensity computation.
          These simplifications are aimed at speeding up MP2 runs,
          if one cares not so much about intensities, and so would
          like to reduce CPU for computing dipoles.  It is pointless
          to select DMDR for SCF electronic structure, where the
          dipoles are easily obtainable.  DMDR must not be used if
          overtones are being computed.

          DMDR   = if true, indicates that the harmonic dipole
                   derivative tensor $DIPDR is input, rather than
                   computing the dipoles.  (default is .FALSE.)

          MPDIP  = for MP2 electronic structure, a value of .FALSE.
                   uses SCF level dipoles in order to save the time
                   needed to obtain the MP2 density at every grid
                   point.  It is more accurate to use the DMDR flag
                   instead of this option, if $DIPDR is available.
                   Obviously this variable is irrelevant for SCF
                   level electronic structure.  (default=.TRUE.)

          VCFCT  = scaling factor for pair-coupling potential.
                   Sometimes when pair-coupling potential values
                   are larger than the corresponding single mode
                   values, they must be scaled down.  (Default=1.0)

             G.M.Chaban, J.O.Jung, R.B.Gerber
             J.Chem.Phys. 111, 1823-1829(1999)



$VIBSCF group (optional, relevant to RUNTYP=VSCF)

          This is restart data, as written to file IRCDATA in a
          partially completed previous run.  Append a " $END" line,
          and select READV=.TRUE. to read the data.



$DRC group (relevant for RUNTYP=DRC)

              This group governs the dynamical reaction coordinate,
          a classical trajectory method based on quantum chemical
          potential energy surfaces.  In GAMESS these may be either
          ab initio or semi-empirical.  Because the vibrational
          period of a normal mode with frequency 500 wavenumbers is
          67 fs, a DRC needs to run for many steps in order to
          sample a representative portion of phase space.  Almost
          all DRCs break molecular symmetry, so build your molecule
          with C1 symmetry in $DATA, or specify NOSYM=1 in $CONTRL.
          Restart data can be found in the job's OUTPUT file, with
          important results summarized to the IRCDATA file.

          NSTEP  = The number of DRC points to be calculated, not
                   including the initial point.  (default = 1000)

          DELTAT = is the time step.  (default = 0.1 fs)

          TOTIME = total duration of the DRC computed in a previous
                   job, in fs.  The default is the correct value
                   when initiating a DRC.  (default=0.0 fs)

                                     * * *

                In general, a DRC can be initiated anywhere,
                so $DATA might contain coordinates of the
                equilibrium geometry, or a nearby transition
                state, or something else.  You must also
                supply an initial kinetic energy, and the
                direction of the initial velocity, for which
                there are a number of options:

	  EKIN   = The initial kinetic energy (default = 0.0 kcal/mol)
                   See also ENM, NVEL, and VIBLVL regarding alternate
                   ways to specify the initial value.

          VEL    = an array of velocity components, in Bohr/fs.
                   When NVEL is false, this is simply the direction
                   of the velocity vector.  Its magnitude will be
                   automatically adjusted to match the desired initial
                   kinetic energy, and it will be projected so that
                   the translation of the center of mass is removed.
                   Give in the order vx1, vy1, vz1, vx2, vy2, ...

	  NVEL   = a flag to compute the initial kinetic energy from
                   the input VEL using the sum of mass*VEL*VEL/2.
                   This flag is usually selected only for restarts.

                                     * * *

                   The next two allow the kinetic energy to be
                   partitioned over all normal modes.  The
                   coordinates in $DATA are likely to be from
                   a stationary point!  You must also supply a
                   $HESS group.

          VIBLVL = a flag to turn this option on (default=.FALSE.)

          VIBENG = an array of energies (in units of multiples of
                   the hv of each mode) to be imparted along each
                   normal mode.  The default is to assign the zero
                   point energy only, VIBENG(1)=0.5, 0.5, ..., 0.5.
                   If given as a negative number, the initial
                   direction of the velocity vector is along the
                   reverse direction of the mode.  "Reverse" means
                   the phase of the normal mode is chosen such that
                   the largest magnitude component is a negative
                   value.  An example might be VIBENG(4)=2.5 to add
                   two quanta to mode 4, along with zero point
                   energy in all modes.

                                     * * *

                   The next three pertain to initiating the DRC
                   along a single normal mode of vibration.  No
                   kinetic energy is assigned to the other modes.
                   You must also supply a $HESS group.

	  NNM    = The number of the normal mode to which the initial
		   kinetic energy is given.  The absolute value of NNM
                   must be in the range 1, 2, ..., 3N-6.  If NNM is a
                   positive/negative value, the initial velocity will
                   lie in the forward/reverse direction of the mode.
                   "Forward" means the largest component of the normal
                   mode is a positive value.  (default=0)

	  ENM    = the initial kinetic energy given to mode NNM,
                   in units of vibrational quanta hv, so the amount
                   depends on mode NNM's vibrational frequency, v.
                   If you prefer to impart an arbitrary initial
                   kinetic energy to mode NNM, specify EKIN instead.
                   (default = 0.0 quanta)

                                     * * *

          To summarize, there are five different ways to specify the
          DRC trajectory:

             1. VEL vector with NVEL=.TRUE.  This is difficult to
                specify at your initial point, and so this option
                is mainly used when restarting your trajectory.
                The restart information is always in this format.
             2. VEL vector and EKIN with NVEL=.FALSE.  This will
                give a desired amount of kinetic energy in the
                direction of the velocity vector.
             3. VIBLVL and VIBENG selected, to give initial kinetic
                energy to all of the normal modes.
             4. NNM and ENM to give quanta to a single normal mode.
             5. NNM and EKIN to give arbitrary kinetic energy to
                a single normal mode.

                                     * * *

                 The most common use of the next two is to
                 analyze a trajectory with respect to the
                 minimum energy geometry the trajectory is
                 traveling around.

          NMANAL = a flag to select mapping of the mass-weighted
                   Cartesian DRC coordinates and velocity (conjugate
                   momentum) in terms of normal modes step by step.
                   If you choose this option, you must supply both
                   C0 and a $HESS group from the stationary point.

          C0     = an array of the coordinates of the stationary
                   point (the coordinates in $DATA might well be
                   some other coordinates).  Give in the order
                                     * * *

               The next option applies to all input paths which
               read a hessian: NMANAL, NNM, or VIBLVL.  After
               the translations and rotations have been dropped,
               the normal modes are renumbered 1, 2, ..., 3N-6.

          HESSTS = a flag to say if the hessian corresponds to a
                   transition state or a minimum.  This parameter
                   controls deletion of the translation and rotation
                   degrees of freedom, i.e. the default is to drop
                   the first six "modes", while setting this flag
                   on drops modes 2 to 7 instead. (default=.FALSE.)

              The next variables can cause termination of a run, if
          molecular fragments get too far apart or close together.

          NFRGPR = Number of atom pairs whose distance will be
                   checked.  (default is 0)

          IFRGPR = Array of the atom pairs.  2 times NFRGPR values.

          FRGCUT = Array for a boundary distance (in Bohr) for atom
                   pairs to end DRC calculations.  The run will
                   stop if any distance exceeds the tolerance, or if
                   a value is given as a negative number, if the
                   distance becomes shorter than the absolute value.
                   In case the trajectory starts outside the bounds
                   specified, they do not apply until after the
                   trajectory reaches a point where the criteria
                   are satisfied, and then goes outside again.
                   Give NFRGPR values.

                                     * * *

              The final variables control the volume of output.
              Let F mean the first DRC point found in this run,
              and L mean the last DRC point of this run.

          NPRTSM = summarize the DRC results every NPRTSM steps,
                   to the file IRCDATA.  (default = 1)

          NPRT   =  1  Print orbitals at all DRC points
                    0  Print orbitals at F+L (default)
                   -1  Never print orbitals

          NPUN   =  2  Punch all orbitals at all DRC points
                    1  Punch all orbitals at F+L, and occupied
                       orbitals at DRC points between
                    0  Punch all orbitals at F+L only (default)
                   -1  Never punch orbitals



          J.J.P.Stewart, L.P.Davis, L.W.Burggraf,
              J.Comput.Chem. 8, 1117-1123 (1987)
          S.A.Maluendes, M.Dupuis,  J.Chem.Phys. 93, 5902-5911(1990)
          T.Taketsugu, M.S.Gordon,  J.Phys.Chem. 99, 8462-8471(1995)
          T.Taketsugu, M.S.Gordon,  J.Phys.Chem. 99, 14597-604(1995)
          T.Taketsugu, M.S.Gordon,  J.Chem.Phys. 103, 10042-9(1995)
          T.Taketsugu, M.S.Gordon,  J.Chem.Phys. 104, 2834-40(1996)
          M.S.Gordon, G.Chaban, T.Taketsugu
              J.Phys.Chem. 100, 11512-11525(1996)



$GLOBOP group (optional, relevant to RUNTYP=GLOBOP)

              This controls a search for the global minimum energy
          when effective fragment are being used.  It has options
          for a single temperature Monte Carlo search, or a multi-
          temperature simulated annealing.  Local minimization of
          some or all of the structures selected by the Monte Carlo
          is optional.  See REFS.DOC for an overview of this RUNTYP.

          TEMPI  =  initial temperature used in the simulation.
                    (default = 40000 K)

          TEMPF  =  final temperature. If TEMPF is not given and
                    NTEMPS is greater than 1, TEMPF will be
                    calculated based on a cooling factor of 0.95.

          NTEMPS =  number of temperatures used in the simulation.
                    If NTEMPS is not given but TEMPF is given,
                    NTEMP will be calculated based on a cooling
                    factor of 0.95. If neither NTEMP nor TEMPF is
                    given, the job defaults to a single temperature
                    Monte Carlo calculation.

          NFRMOV =  number of fragments to move on each step.

          NGEOPT =  number of geometries to be evaluated at each
                    temperature.  (default = 10)

          NTRAN  =  number of translational steps in each block.

          NROT   =  number of rotational steps in each block.

          NBLOCK =  the number of blocks of steps can be set directly
                    with this variable, instead of being calculated
                    from NGEOPT, NTRAN, and NROT.  If NBLOCK is input,
                    the number of geometries at each temperature will
                    be NBLOCK*(NTRAN+NROT).  Each block has NTRAN
                    translational steps followed by NROT rotational

          MCMIN  =  flag to enable geometry optimization to minimize
                    the energy is carried out every NSTMIN steps.

          NSTMIN =  After this number of geometry steps are taken, a
                    local (Newton-Raphson) optimization will be
                    carried out.  If this variable is set to 1, a
                    local minimization is carried out on every step,
                    reducing the MC space to the set of local minima.
                    Irrelevant if MCMIN is false.  (default=10)

          OPTN   =  if set to .TRUE., at the end of the run local
                    minimizations are carried out on the final
                    geometry and on the minimum-energy geometry.

          SCALE  =  an array of length two.  The first element is the
                    initial maximum step size for the translational
                    coordinates (Angstroms).  The second element is
                    the initial maximum stepsize for the rotational
                    coordinates (pi-radians). (defaults = 1,1)

          SMODIF =  scale factor for moving ab initio atoms in the
                    MC simulation.  If set to zero, the ab initio
                    atoms do not move. (default=0.1)

          ALPHA  =  controls the rate at which information from
                    successful steps is folded into the maximum step
                    sizes for each of the 6*(number of fragments)
                    coordinates.  ALPHA varies between 0 and 1.
                    ALPHA=0 means do not change the maximum step
                    sizes, and ALPHA=1 throws out the old step sizes
                    whenever there is a successful step and uses the
                    successful step sizes as the new maxima.  This
                    update scheme was used with the Parks method
                    where all fragments are moved on every step.  It
                    is normally not used with the Metropolis method.
                    (default = 0)

          DACRAT =  the desired acceptance ratio, the program tries
                    to achieve this by adjusting the maximum step
                    size.  (default = 0.5)

          UPDFAC =  the factor used to update the maximum step size
                    in the attempt to achive the desired acceptance
                    ratio (DACRAT).  If the acceptance ratio at the
                    previous temperature was below DACRAT, the step
                    size is decreased by multiplying it by UPDFAC.
                    If the acceptance ratio was above DACRAT, the
                    step size is increased by dividing it by DACRAT
                    It should be between 0 and 1. (default = 0.95)

          SEPTOL =  the separation tolerence between atoms in the ab
                    initio piece and atoms in the fragments, as well
                    as between atoms in different fragments.  If a
                    step moves atoms closer than this tolerence, the
                    step is rejected. (default = 1.5 Angstroms)

          XMIN, XMAX, YMIN, YMAX, ZMIN, ZMAX = mimimum and maximum
                    values for the Cartesian coordinates of the
                    fragment.  If the first point in a fragment steps
                    outside these boundaries, periodic boundary
                    conditions are used and the fragment re-enters on
                    the opposite side of the box.  The defaults of
                    -10 for minima and +10 for maxima should usually
                    be changed.

          BOLTWT =  method for calculating the Boltzmann factor,
                    which is used as the probability of accepting a
                    step that increases the energy.
                 =  STANDARD = use the standard Boltzmann factor,
                    exp(-delta(E)/kT)  (default)
                 =  AVESTEP = scale the temperature by the average
                    step size, as recommended in the Parks reference
                    when using values of ALPHA greater than 0.

          NPRT   =  controls the amount of output, with
                 = -2 reduces output below that of -1
                 = -1 reduces output further, needed for MCMIN=.true.
                 =  0 gives minimal output (default)
                 =  1 gives the normal GAMESS amount of output
                 =  2 gives maximum output
                    For large simulations, even IOUT=0 may produce
                    a log file too large to work with easily.

          RANDOM =  controls the choice of random number generator.
                 =  DEBUG uses a simple random number generator with
                    a constant seed. Since the same sequence of
                    random numbers is generated during each job, it
                    is useful for debugging.
                 =  RAND1 uses the simple random number generator
                    used in DEBUG, but with a variable seed.
                 =  RAND3 uses a more sophisticated random number
                    generator described in Numerical Recipes, with a
                    variable seed (default).

          IFXFRG =  array whose length is the number of fragments.
                    It allows one or more fragments to be fixed
                    during the simulation.
                 =0 allows the fragment to move during the run
                 =1 fixes the fragment
                    For example, IFXFRG(3)=1 would fix the third
                    fragment, the default is IFXFRG(1)=0,0,0,...,0

          MOVIE2 = a flag to create a series of structural data
                   which can be shown as a movie by the MacIntosh
                   program Chem3D.  The coordinates of each accepted
                   geometry are written.  The data is written to the
                   file IRCDATA.  (default=.FALSE.)



$GRADEX group (optional, for RUNTYP=GRADEXTR)

             This group controls the gradient extremal following
          algorithm.  The GEs leave stationary points parallel to
          each of the normal modes of the hessian.  Sometimes a GE
          leaving a minimum will find a transition state, and thus
          provides us with a way of finding that saddle point.  GEs
          have many unusual mathematical properties, and you should
          be aware that they normally differ a great deal from IRCs.

             The search will always be performed in cartesian
          coordinates, but internal coordinates along the way may
          be printed by the usual specification of NZVAR and $ZMAT.

          METHOD = algorithm selection.
                   SR   A predictor-corrector method due to Sun
                        and Ruedenberg (default).
                   JJH  A method due to Jorgensen, Jensen and

          NSTEP  = maximum number of predictor steps to take.

          DPRED  = the stepsize for the predictor step.
                   (default = 0.10)

          STPT   = a flag to indicate whether the initial geometry
                   is considered a stationary point. If .TRUE.,
                   the geometry will be perturbed by STSTEP along
                   the IFOLOW normal mode.
                   (default = .TRUE.)

          STSTEP = the stepsize for jumping away from a stationary
                   point. (default = 0.01)

          IFOLOW = Mode selection option.  (default is 1)
                   If STPT=.TRUE., the intial geometry will be
                   perturbed by STSTEP along the IFOLOW normal mode.
                   Note that IFOLOW can be positive or negative,
                   depending on the direction the normal mode
                   should be followed in. The positive direction
                   is defined as the one where the largest component
                   of the Hessian eigenvector is positive.

                   If STPT=.FALSE. the sign of IFOLOW determines
                   which direction the GE is followed in. A positive
                   value will follow the GE in the uphill direction.
                   The value of IFOLOW should be set to the Hessian
                   mode which is parallel to the gradient to avoid
                   miscellaneous warning messages.

          GOFRST = a flag to indicate whether the algorithm should
                   attempt to locate a stationary point.  If .TRUE.,
                   a straight NR search is performed once the NR
                   step length drops below SNRMAX.  10 NR step are
                   othen allowed, a value which cannot be changed.
                   (default = .TRUE.)

          SNRMAX = upper limit for switching to straight NR search
                   for stationary point location.
                   (default = 0.10 or DPRED, whichever is smallest)

          OPTTOL = gradient convergence tolerance, in Hartree/Bohr.
                   Used for optimizing to a stationary point.
                   Convergence of a geometry search requires the
                   rms gradient to be less than OPTTOL.

          HESS   = selection of the initial hessian matrix, if
                 = READ causes the hessian to be read from a $HESS
                 = CALC causes the hessian to be computed. (default)

          ---- parameters on this page apply only to METHOD=SR ----

          DELCOR = the corrector step should be smaller than this
                   value before the next predictor step is taken.
                   (default = 0.001)

          MYSTEP = maximum number of micro iteration allowed to
                   bring the corrector step length below DELCOR.

          SNUMH  = stepsize used in the numerical differentiation
                   of the Hessian to produce third derivatives.
                   (default = 0.0001)

          HSDFDB = flag to select determination of third derivatives.
                   At the current geometry we need the gradient, the
                   Hessian, and the partial third derivative matrix
                   in the gradient direction.

                   If .TRUE., the gradient is calculated at the
                   current geometry, and two Hessians are calculated
                   at SNUMH distance to each side in the gradient
                   direction.  The Hessian at the geometry is formed
                   as the average of the two displaced Hessians.

                   If .FALSE., both the gradient and Hessian are
                   calculated at the current geometry, and one
                   additional Hessian is calculated at SNUMH in the
                   gradient direction.

                   The default double-sided differentiation produces
                   a more accurate third derivative matrix, at the
                   cost of an additional wave function and gradient.
                   (default = .TRUE.)



$SURF group (relevant for RUNTYP=SURFACE)

              This group allows you to probe a potential energy
          surface along a small grid of points.  Note that there is
          no option to vary angles, only distances.  The scan can
          be made for any SCFTYP, or for the MP2 or CI surface.

          IVEC1  = an array of two atoms, defining a coordinate from
                   the first atom given, to the second.

          IGRP1  = an array specifying a group of atoms, which must
                   include the second atom given in IVEC1.  The
                   entire group will be translated (rigidly) along
                   the vector IVEC1, relative to the first atom
                   given in IVEC1.

          ORIG1  = starting value of the coordinate, which may be
                   positive or negative.  Zero corresponds to the
                   distance given in $DATA.

          DISP1  = step size for the coordinate.

          NDISP1 = number of steps to take for this coordinate.

                  There are no reasonable defaults for these
                  keywords, so you should input all of them.
                  ORIG1 and DISP1 should be given in Angstrom.

          IVEC2, IGRP2, ORIG2, DISP2, NDISP2 = have the identical
          meaning as their "1" counterparts, and permit you to make
          a two dimensional map along two displacement coordinates.
          If the "2" data are not input, the surface map proceeds in
          only one dimension.

          Note that properties are not computed at these points,
          other than the energy.



$LOCAL group (relevant for LOCAL=RUEDNBRG, BOYS, or POP)

              This group allows input of additional data to control
          the localization methods.  If no input is provided, the
          valence orbitals will be localized as much as possible,
          while still leaving the wavefunction invariant.

          PRTLOC = a flag to control supplemental printout.  The
                   extra output is the rotation matrix to the
                   localized orbitals, and, for the Boys method,
                   the orbital centroids, for the Ruedenberg
                   method, the coulomb and exchange matrices,
                   for the population method, atomic populations.

          MAXLOC = maximum number of localization cycles.  This
                   applies to BOYS or POP methods only.  If the
                   localization fails to converge, a different
                   order of 2x2 pairwise rotations will be tried.

          CVGLOC = convergence criterion.  The default provides
                   LMO coefficients accurate to 6 figures.

          SYMLOC = a flag to restrict localization so that
                   orbitals of different symmetry types are not
                   mixed.  This option is not supported in all
                   possible point groups.  The purpose of this
                   option is to give a better choice for the
                   starting orbitals for GVB-PP or MCSCF runs,
                   without destroying the orbital's symmetry.
                   This option is compatible with each of the
                   3 methods of selecting the orbitals to be
                   included.  (default=.FALSE.)

              These parameters select the orbitals which are
              to be included in the localization.  You may
              select from FCORE, NOUTA/NOUTB, or NINA/NINB,
              but may choose only one of these.

          FCORE  = flag to freeze all the chemical core orbitals
                   present.   All the valence orbitals will be
                   localized.  (default=.TRUE.)

                                 * * *

          NOUTA  = number of alpha orbitals to hold fixed in the
                   localization.  (default=0)

          MOOUTA = an array of NOUTA elements giving the numbers of
                   the orbitals to hold fixed.  For example, the
                   input NOUTA=2 MOOUTA(1)=8,13 will freeze only
                   orbitals 8 and 13.  You must enter all the
                   orbitals you want to freeze, including any cores.
                   This variable has nothing to do with cows.

          NOUTB =  number of beta orbitals to hold fixed in -UHF-
                   localizations.  (default=0)

          MOOUTB = same as MOOUTA, except that it applies to the
                   beta orbitals, in -UHF- wavefunctions only.

                                 * * *

          NINA   = number of alpha orbitals which are to be
                   included in the localization.  (default=0)

          MOINA  = an array of NINA elements giving the numbers of
                   the orbitals to be included in the localization.
                   Any orbitals not mentioned will be frozen.

          NINB   = number of -UHF- beta MOs in the localization.

          MOINB  = same as MOINA, except that it applies to the
                   beta orbitals, in -UHF- wavefunctions only.

          N.B.  Since Boys localization needs the dipole integrals,
                do not turn off dipole moment calculation in $ELMOM.

          ----- The following keywords are used for the localized
                charge distribution (LCD) energy decomposition.

          EDCOMP = flag to turn on LCD energy decomposition.
                   Note that this method is currently implemented
                   for SCFTYP=RHF and ROHF and LOCAL=RUEDNBRG only.
                   The SCF LCD forces all orbitals to be localized,
                   overriding input on the previous page.  See also
                   LMOMP2 in the $MP2 group.  (default = .FALSE.)

          MOIDON = flag to turn on LMO identification and subsequent
                   LMO reordering, and assign nuclear LCD automat-
                   ically.  (default = .FALSE.)

          DIPDCM = flag for LCD molecular dipole decomposition.
                   (default = .FALSE.)

          QADDCM = flag for LCD molecular quadrupole decomposition.
                   (default = .FALSE.)

          POLDCM = flag to turn on LCD polarizability decomposition.
                   This method is implemented for SCFTYP=RHF or ROHF
                   and LOCAL=BOYS or RUEDNBRG. (default=.FALSE.)

          POLNUM = flag to forces numerical rather than analytical
                   calculation of the polarizabilities.  This may be
                   useful in larger molecules.  The numerical
                   polarizabilities of bonds in or around aromatic
                   rings sometimes are unphysical. (default=.FALSE.)
                   See D.R.Garmer, W.J.Stevens
                       J.Phys.Chem. 93, 8263-8270 (1989).

          POLAPP = flag to force calculation of the polarizabilities
                   using a perturbation theory expression.  This may
                   be useful in larger molecules. (default=.FALSE.)
                   See R.M. Minikis, V. Kairys, J.H. Jensen
                       J.Phys.Chem.A 2000, submitted.

          POLANG = flag to choose units of localized polarizability
                   output. The default is Angstroms**3, while false
                   will give Bohr**3.  (default=.TRUE.)

          ZDO    = flag for LCD analysis of a composite wave function,
                   given in a $VEC group of a van der Waals complex,
                   within the zero differential overlap approximation.
                   The MOs are not orthonormalized and the inter-
                   molecular electron exchange energy is neglected.
                   In addition, the molecular overlap matrix is printed
                   out.  This is a very specialized option.
                   (default = .FALSE.)

          ----- The remaining keywords can be used to define the
                nuclear part of an LCD.  They are usually used to
                rectify mistakes in the automatic definition
                made when MOIDON=.TRUE.  The index defining the
                LMO number then refers to the reordered list of LMOs.

          NNUCMO = array giving the number of nuclei assigned to a
                   particular LMO.

          IJMO   = is an array of pairs of indices (I,J), giving
                   the row (nucleus I) and column (orbital J)
                   index of the entries in ZIJ and MOIJ.

          MOIJ   = arrays of integers K, assigning nucleus K as the
                   site of the Ith charge of LCD J.

          ZIJ    = array of floating point numbers assigning a
                   charge to the Ith charge of LCD J.

          IPROT  = array of integers K, defining nucleus K as a

          DEPRNT = a flag for additional decomposition printing,
                   such as pair contributions to various energy
                   terms, and centroids of the Ruedenberg orbitals.
                   (default = .FALSE.)



$TWOEI group (relevant for EDCOMP=.TRUE. in $LOCAL)

             Formatted transformed two-electron Coulomb and Exchange
          integrals as punched during a LOCAL=RUEDNBRG run.  If this
          group is present it will automaticall be read in during
          such a run and the two-electron integrals do not have to
          be re-transformed.  This group is especially useful for
          EDCOMP=.TRUE. runs when the localization has to be repeated
          for different definitions of nuclear LCDs.



$TRUNCN group (optional, relevant for RHF)

              This group controls the truncation of some of the
          localized orbitals to just the AOs on a subset of the
          atoms.  This option is particularly useful to generate
          localized orbitals to be frozen when the effective
          fragment potential is used to partition a system across a
          chemical bond.  In other words, this group prepares the
          frozen buffer zone orbitals.  This group should be used in
          conjunction with RUNTYP=ENERGY (or PROP if the orbitals
          are available) and either LOCAL=RUEDNBRG or BOYS, with
          MOIDON set in $LOCAL.

          DOPROJ = flag to activate MO projection/truncation, the
                   default is to skip this (default=.FALSE.)

          AUTOID = forces identification of MOs (analogous to MOIDON
                   in $LOCAL).  This keyword is provided in case the
                   localized orbitals are already present in $VEC,
                   in which case this is a faster RUNTYP=PROP with
                   LOCAL=NONE job.  Obviously, GUESS=MOREAD.

          PLAIN  = flag to control the MO tail truncation.  A value
                   of .FALSE. uses corresponding orbital projections,
                   H.F.King, R.E.Stanton, H.Kim, R.E.Wyatt, R.G.Parr
                   J. Chem. Phys. 47, 1936-1941(1967) and generates
                   orthogonal orbitals.  A value of .TRUE. just sets
                   the unwanted AOs to zero, so the resulting MOs
                   need to go through the automatic orthogonalization
                   step when MOREAD in the next job. (default=.FALSE.)

          IMOPR  = an array specifying which MOs to be truncated. In
                   most cases involving normal bonding, the options
                   MOIDON or AUTOID will correctly identify all
                   localized MOs belonging to the atoms in the zone
                   being truncated.  However, you can inspect the
                   output, and give a list of all MOs which you want
                   to be truncated in this array, in case you feel
                   the automatic assignment is incorrect.
                   Any orbital not in the truncation set, whether
                   this is chosen automatically or by IMOPR, is left
                   completely unaltered.

                                  - - -

          There are now two ways to specify what orbitals are to
          be truncated.  The most common usage is for preparation of
          a buffer zone for QM/MM computations, with an Effective
          Fragment Potential representing the non-quantum part of
          the system.  This input is NATAB, NATBF, ICAPFR, ICAPBF,
          in which case the $DATA input must be sorted into three
          zones.  The first group of atoms are meant to be treated
          in later runs by full quantum mechanics, the second
          group by frozen localized orbitals as a 'buffer', and the
          third group is to be substituted later by an effective
          fragment potential (multipoles, polarizabilities, ...).
          Note that in the DOPROJ=.TRUE. run, all atoms are still
          quantum atoms.

          NATAB  = number of atoms to be in the 'ab initio' zone.

          NATBF  = number of atoms to be in the 'buffer' zone.
                   The program can obtain the number of atoms in
                   the remaining zone by subtraction, so it need
                   not be input.

          In case the MOIDON or AUTOID options lead to confused
          assignments (unlikely in ordinary bonding situations
          around the buffer zone), there are two fine tuning values.

          ICAPFR = array indicating the identity of "capping atoms"
                   which are on the border between the ab initio and
                   buffer zones (in the ab initio zone).

          ICAPBK = array indicating the identity of "capping atoms"
                   which are on the border between the buffer and EFP
                   zones (in the effective fragment zone).

          See also IXCORL and IXLONE below.

                                  - - -

          In case truncation seems useful for some other purpose,
          you can specify the atoms in any order within the $DATA
          group, by the IZAT/ILAT approach.  You are supposed to
          give only one of these two lists, probably whichever is

          IZAT   = an array containing the atoms which are NOT in
                   the buffer zone.

          ILAT   = an array containing the atoms which are in
                   the buffer zone.

          The AO coefficients of the localized orbitals present in
          the buffer zone which lie on atoms outside the buffer will
          be truncated.

          See also IXCORL and IXLONE below.

                                  - - -

          The next two values let you remove additional orbitals
          within the buffer zone from the truncation process, if that
          is desirable.  These arrays can only include atoms that are
          already in the buffer zone, whether this was defined by
          NATBF, or IZAT/ILAT.  The default is to include all core
          and lone pair orbitals, not just bonding orbitals, as the
          buffer zone orbitals.

          IXCORL = an array of atoms whose core and lone pair
                   orbitals are to be considered as not belonging
                   to the buffer zone orbitals.

          IXLONE = an array of atoms for which only the lone pair
                   orbitals are to be considered as not belonging
                   to the buffer zone orbitals.

                                  - - -

          The final option controls output of the truncated orbitals
          to file PUNCH for use in later runs:

          NPUNOP =    punch out option for the truncated orbitals
                 = 1  the MOs are not reordered.
                 = 2  punch the truncated MOs as the first vectors
                      in the $VEC MO set, with untransformed vectors
                      following immediately after. (default)



$ELMOM group (not required)

          This group controls electrostatic moments calculation.

          IEMOM  = 0 - skip this property
                   1 - calculate monopole and dipole (default)
                   2 - also calculate quadrupole moments
                   3 - also calculate octupole moments

          WHERE  = COMASS   - center of mass (default)
                   NUCLEI   - at each nucleus
                   POINTS   - at points given in $POINTS.

          OUTPUT = PUNCH, PAPER, or BOTH (default)

          IEMINT = 0 - skip printing of integrals (default)
                   1 - print dipole integrals
                   2 - also print quadrupole integrals
                   3 - also print octupole integrals
                  -2 - print quadrupole integrals only
                  -3 - print octupole integrals only

              The quadrupole and octupole tensors on the printout
          are formed according to the definition of Buckingham.
          Caution: only the first nonvanishing term in the multi-
          ipole charge expansion is independent of the coordinate
          origin chosen, which is normally the center of mass.



$ELPOT group (not required)

          This group controls electrostatic potential calculation.

          IEPOT = 0 skip this property (default)
                  1 calculate electric potential

          WHERE  = COMASS   - center of mass
                   NUCLEI   - at each nucleus (default)
                   POINTS   - at points given in $POINTS
                   GRID     - at grid given in $GRID
                   PDC      - at points controlled by $PDC.

          OUTPUT = PUNCH, PAPER, or BOTH (default)

              This property is the electrostatic potential V(a) felt
          by a test positive charge, due to the molecular charge
          density.  A nucleus at the evaluation point is ignored.
          If this property is evaluated at the nuclei, it obeys the
               sum on nuclei(a)   Z(a)*V(a) = 2*V(nn) + V(ne).
          The electronic portion of this property is called the
          diamagnetic shielding.



$ELDENS group (not required)

          This group controls electron density calculation.

          IEDEN  = 0 skip this property (default)
                 = 1 compute the electron density.
                 = 2 compute the orbital values rather than the
                     density (TU Graz extension).

          MORB   = The molecular orbital whose electron density is
                   to be computed.  If zero, the total density is
                   computed.  (default=0)

          WHERE  = COMASS   - center of mass
                   NUCLEI   - at each nucleus (default)
                   POINTS   - at points given in $POINTS
                   GRID     - at grid given in $GRID

          OUTPUT = PUNCH, PAPER, or BOTH (default)

          IEDINT = 0 - skip printing of integrals (default)
                   1 - print the electron density integrals



$ELFLDG group (not required)

              This group controls electrostatic field and electric
          field gradient calculation.

          IEFLD  = 0 - skip this property (default)
                   1 - calculate field
                   2 - calculate field and gradient

          WHERE  = COMASS   - center of mass
                   NUCLEI   - at each nucleus (default)
                   POINTS   - at points given in $POINTS

          OUTPUT = PUNCH, PAPER, or BOTH (default)

          IEFINT = 0 - skip printing these integrals (default)
                   1 - print electric field integrals
                   2 - also print field gradient integrals
                  -2 - print field gradient integrals only

          The Hellman-Feynman force on a nucleus is the nuclear
          charge multiplied by the electric field at that nucleus.
          The electric field is the gradient of the electric
          potential, and the field gradient is the hessian of the
          electric potential.  The components of the electric field
          gradient tensor are formed in the conventional way, i.e.
          see D.Neumann and J.W.Moskowitz.



$POINTS group (not required)

              This group is used to input points at which properties
          will be computed.  This first card in the group must
          contain the string ANGS or BOHR, followed by an integer
          NPOINT, the number of points to be used.  The next NPOINT
          cards are read in free format, containing the X, Y, and Z
          coordinates of each desired point.



$GRID group (not required)

              This group is used to input a grid (plane through the
          molecule) on which properties will be calculated.

          ORIGIN(i) = coordinates of the lower left corner of
                      the plot.
          XVEC(i)   = coordinates of the lower right corner of
                      the plot.
          YVEC(i)   = coordinates of the upper left corner of
                      the plot.
          SIZE      = grid increment, default is 0.25.
          UNITS     = units of the above four values, it can be
                      either BOHR or ANGS (the default).

          Note that XVEC and YVEC are not necessarily parallel to
          the X and Y axes, rather they are the axes which you
          desire to see plotted by the MEPMAP contouring program.



$PDC group (relevant if WHERE=PDC in $ELPOT)

               This group determines the points at which to compute
          the electrostatic potential, for the purpose of fitting
          atomic charges to this potential.  Constraints on the fit
          which determines these "potential determined charges" can
          include the conservation of charge, the dipole, and the

          PTSEL  = determines the points to be used, choose from
	           GEODESIC to use a set of points on several fused
		         sphere van der Waals surfaces, with points
			 selected using an algorithm due to Mark
			 Spackman.  The results are similar to those
			 from the Kollman/Singh method, but are
			 less rotation dependent. (default)
	           CONNOLLY to use a set of points on several fused
		         sphere van der Waals surfaces, with points
			 selected using an algorithm due to Michael
			 Connolly.  This is identical to the method
			 used by Kollman & Singh (see below)
                   CHELPG to use a modified version of the CHELPG
                         algorithm, which produces a symmetric
                         grid of points for a symmetric molecule.

          CONSTR = NONE   - no fit is performed.  The potential at
                            the points is instead output according
                            to OUTPUT in $ELPOT.
                   CHARGE - the sum of fitted atomic charges is
                            constrained to reproduce the total
                            molecular charge. (default)
                   DIPOLE - fitted charges are constrained to
                            exactly reproduce the total charge
                            and dipole.
                   QUPOLE - fitted charges are constrained to
                            exactly reproduce the charge, dipole,
                            and quadrupole.

              Note: the number of constraints cannot exceed
              the number of parameters, which is the number
              of nuclei.  Planar molecules afford fewer
              constraint equations, namedly two dipole
              constraints and three quadrupole constraints,
              instead of three and five, repectively.

          * * * the next 5 pertain to PTSEL=GEODESIC or CONNOLLY * * *

          VDWSCL = scale factor for the first shell of VDW spheres.
                   The default of 1.4 seems to be an empirical best
                   value. Values for VDW radii for most elements up
                   to Z=36 are internally stored.

          VDWINC = increment for successive shells (default = 0.2).
		   The defaults for VDWSCL and VDWINC will result
		   in points chosen on layers at 1.4, 1.6, 1.8 etc
		   times the VDW radii of the atoms.

	  LAYER  = number of layers of points chosen on successive
		   fused sphere VDW surfaces (default = 4)

          NFREQ  = flag for particular geodesic tesselation of
		   points.  Only relevant if PTSEL=GEODESIC.
		   Options are:
                    (10*h + k)  for   {3,5+}h,k tesselations
                   -(10*h + k)  for   {5+,3}h,k tesselations
                   (of course both nh and nk must be less than 10,
		   so NFREQ must lie within the range -99 to 99)
   		   The default value is NFREQ=30 (=03)

	  PTDENS = density of points on the surface of each scaled
		   VDW sphere (in points per square au).  Only relevant
		   if PTSEL=CONNOLLY.  Default is 0.28 per au squared,
		   which corresponds to 1.0 per square Angstrom, the
		   default recommended by Kollman & Singh.

             * * * the next two pertain to PTSEL=CHELPG * * *

          RMAX   = maximum distance from any point to the closest
                   atom.  (default=3.0 Angstroms)

          DELR   = distance between points on the grid.
                   (default=0.8 Angstroms)

          MAXPDC = an estimate of the total number of points whose
                   electrostatic potential will be included in the
                   fit. (default=10000)

                                 * * *

          CENTER = an array of coordinates at which the moments were

          DPOLE  = the molecular dipole.

          QPOLE  = the molecular quadrupole.

          PDUNIT = units for the above values.  ANGS (default) will
                   mean that the coordinates are in Angstroms, the
                   dipole in Debye, and quadrupole in Buckinghams.
                   BOHR implies atomic units for all 3.

            Note: it is easier to compute the moments in the
            current run, by setting IEMOM to at least 2 in
            $ELMOM.  However, you could fit experimental data,
            for example, by reading it in here.


               There is no unique way to define fitted atomic
          charges.  Smaller numbers of points at which the electro-
          static potential is fit, changes in VDW radii, asymmetric
          point location, etc. all affect the results.  A useful
          bibliography is

          U.C.Singh, P.A.Kollman, J.Comput.Chem. 5, 129-145(1984)
          L.E.Chirlain, M.M.Francl, J.Comput.Chem. 8, 894-905(1987)
          R.J.Woods, M.Khalil, W.Pell, S.H.Moffatt, V.H.Smith,
             J.Comput.Chem. 11, 297-310(1990)
          C.M.Breneman, K.B.Wiberg, J.Comput.Chem. 11, 361-373(1990)
          K.M.Merz, J.Comput.Chem. 13, 749(1992)
	  M.A.Spackman, J.Comput.Chem. 17, 1-18(1996)



$MOLGRF group (relevant only if you have MOLGRAPH)

             This option provides an interface for viewing orbitals
          through a commercial package named MOLGRAPH, from Daikin
          Industries.  Note that this option uses three disk files
          which are not defined in the GAMESS execution scripts we
          provide, since we don't use MOLGRAPH ourselves.  You will
          need to define files 28, 29, 30, as generic names PRGRID,
          COGRID, MOGRID, of which the latter is passed to MOLGRAPH.

          GRID3D = a flag to generate 3D grid data.
                   (default is .false.).

          TOTAL  = a flag to generate a total density grid data.
                   "Total" means the sum of the orbital densities
                   given by NPLT array.  (default is .false.).

          MESH   = numbers of grids.  You can use different numbers
                   for three axes.  (default is MESH(1)=21,21,21).

          BOUND  = boundary coordinates of a 3D graphical cell.  The
                   default is that the cell is larger than the
                   molecular skeleton by 3 bohr in all directions.
                   E.g., BOUND(1)=xmin,xmax,ymin,ymax,zmin,zmax

          NPLOTS = number of orbitals to be used to generate 3D grid
                   data. (default is NPLOTS=1).

          NPLT   = orbital IDs.  The default is 1 orbital only, the
                   HOMO or SOMO.  If the LOCAL option is given in
                   $CONTRL, localized orbital IDs should be given.
                   For example, NPLT(1)=n1,n2,n3,...

          CHECK  = debug option, printing some of the grid data.

          If you are interested in graphics, look at the WWW page
          for information about other graphics packages with GAMESS.



$STONE group (optional)

              This group defines the expansion points for Stone's
          distributed multipole analysis (DMA) of the electrostatic

              The DMA takes the multipolar expansion of each overlap
          charge density defined by two gaussian primitives, and
          translates it from the center of charge of the overlap
          density to the nearest expansion point.  Some references
          for the method are

              Stone, Chem.Phys.Lett. 83, 233 (1981)
              Price and Stone, Chem.Phys.Lett. 98, 419 (1983)
              Buckingham and Fowler, J.Chem.Phys. 79, 6426 (1983)
              Stone and Alderton, Mol.Phys. 56, 1047 (1985)

              The existence of a $STONE group in the input is what
          triggers the analysis.  Enter as many lines as you wish,
          in any order, terminated by a $END record.


          ATOM i name, where

                ATOM     is a keyword indicating that a particular
                         atom is selected as an expansion center.
                i        is the number of the atom
                name     is an optional name for the atom. If not
                         entered the name will be set to the name
                         used in the $DATA input.


          ATOMS          is a keyword selecting all nuclei in the
                         molecule as expansion points.  No other
                         input on the line is necessary.


          BONDS          is a keyword selecting all bond midpoints
                         in the molecule as expansion points.  No
                         other input on the line is necessary.


          BOND i j name, where

                BOND     is a keyword indicating that a bond mid-
                         point is selected as an expansion center.
                i,j      are the indices of the atoms defining the
                         bond, corresponding to two atoms in $DATA.
                name     an optional name for the bond midpoint.
                         If omitted, it is set to 'BOND'.


          CMASS          is a keyword selecting the center of mass
                         as an expansion point.  No other input on
                         the line is necessary.


          POINT x y z name, where

                POINT    is a keyword indicating that an arbitrary
                         point is selected as an expansion point.
                x,y,z    are the coordinates of the point, in Bohr.
                name     is an optional name for the expansion
                         point.  If omitted, it is set to 'POINT'.


          While making the EFPs for QM/MM run, a single keyword
          QMMMBUF is necessary.  Adding additional keywords may lead
          to meaningless results.  The program will automatically
          select atoms and bond midpoints which are outside the
          buffer zone as the multipole expansion points.

          QMMMBUF  nmo, where

                QMMMBUF  is a keyword specifying the number of QM/MM
                         buffer molecular orbitals, which must be the
                         first NMO orbitals in the MO set.  These
                         orbitals must be frozen in the buffer zone,
                         so this is useful only if $MOFRZ is given.
                NMO      is the number of buffer MO-s
                         (if NMO is omitted, it will be set to the
                         number of frozen MOs in $MOFRZ)


          The second and third moments on the printout can be
          converted to Buckingham's tensors by formula 9 of
            A.D.Buckingham, Quart.Rev. 13, 183-214 (1959)
          These can in turn be converted to spherical tensors
          by the formulae in the appendix of
            S.L.Price, et al.  Mol.Phys. 52, 987-1001 (1984)



$RAMAN group (relevant for all SCFTYPs)

              This input controls the computation of Raman intensity
          by the numerical differentiation produre of Komornicki and
          others.  It is applicable to any wavefunction for which
          the analytic gradient is available, including some MP2 and
          CI cases.  The calculation involves the computation of 19
          nuclear gradients, one without applied electric fields,
          plus 18 no symmetry runs with electric fields applied in
          various directions.  The numerical second differencing
          produces intensity values with 2-3 digits of accuracy.

              This run must follow an earlier RUNTYP=HESSIAN job,
          and the $GRAD and $HESS groups from that first job must be
          given as input.  If the $DIPDR is computed analytically
          by this Hessian job, it too may be read in, if not, the
          numerical Raman job will evaluate $DIPDR.  Once the data
          from the 19 applied fields is available, the $ALPDR tensor
          is evaluated.  Then the nuclear derivatives of the dipole
          moment and alpha polarizability will be combined with the
          normal coordinate information to produce the IR and Raman
          intensity of each mode.

              To study isotopic substitution speedily, input the
          $GRAD, $HESS, $DIPDR, and $ALPDR groups along with the
          desired atomic masses in $MASS.

             The code does not permit any semi-empirical or solvation
          models to be used.

          EFIELD = applied electric field strenth.  The literature
                   suggests values in the range 0.001 to 0.005.
                   (default = 0.002 a.u.)



$ALPDR group (relevant for RUNTYP=RAMAN or HESSIAN)

          Formatted alpha derivative tensor, punched by a previous
          RUNTYP=RAMAN job.  If both $DIPDR and this group are found
          in the input file, the applied field computation will be
          skipped, to immediately evaluate IR and Raman intensities.  

          If this group is found during a Hessian job, the Raman
          intensities will be added to the output.  You might want
          to run as RUNTYP=HESSIAN instead of RUNTYP=RAMAN in order
          to have access to PROJCT or the other options available in
          the $FORCE group.



$MOROKM group (relevant for RUNTYP=MOROKUMA)

              This group controls how the supermolecule input in the
          $DATA group is divided into two or more monomers.  Both
          the supermolecule and its constituent monomers must be
          well described by RHF wavefunctions.

          MOROKM = a flag to request Morokuma-Kitaura decomposition.
                   (default is .TRUE.)

          RVS    = a flag to request "reduced variation space"
                   decomposition.  This differs from the Morokuma
                   option, and one or the other or both may be
                   requested in the same run.  (default is .FALSE.)

          BSSE   = a flag to request basis set superposition error
                   be computed.  You must ensure that CTPSPL is
                   selected.  This option applies only to MOROKM
                   decompositions, as a basis superposition error is
                   automatically generated by the RVS scheme.  This
                   is not the full Boys counterpoise correction, as
                   explained in the reference.  (default is .FALSE.)

                                     * * *

          IATM   = An array giving the number of atoms in each of
                   the monomer.  Up to ten monomers may be defined.
                   Your input in $DATA must have all the atoms in
                   the first monomer defined before the atoms in the
                   second monomer, before the third monomer...  The
                   number of atoms belonging to the final monomer
                   can be omitted.  There is no sensible default for
                   IATM, so don't omit it from your input.

          ICHM   = An array giving the charges of the each monomer.
                   The charge of the final monomer may be omitted,
                   as it is fixed by ICH in $CONTRL, which is the
                   total charge of the supermolecule.  The default
                   is neutral monomers, ICHM(1)=0,0,0,...

          EQUM   = a flag to indicate all monomers are equivalent
                   by symmetry (in addition to containing identical
                   atoms). If so, which is not often true, then only
                   the unique computations will be done.
                   (default is .FALSE.)

          CTPSPL = a flag to decompose the interaction energy into
                   charge transfer plus polarization terms.  This
                   is most appropriate for weakly interacting
                   monomers. (default is .TRUE.)

          CTPLX  = a flag to combine the CT and POL terms into a
                   single term.  If you select this, you might want
                   to turn CTPSPL off to avoid the extra work that
                   that decomposition entails, or you can analyze
                   both ways in the same run (default=.FALSE.)

          RDENG  = a flag to enable restarting, by reading the
                   lines containing "FINAL ENERGY" from a previous
                   run.  The $ENERGY group is single lines read
                   under format A16,F20.10 containing the E, and a
                   card $END to complete.  The 16 chars = anything.
                   (default is .FALSE.)


          The present implementation has some quirks:

          1. The initial guess of the monomer orbitals is not
             controlled by $GUESS.  The program first looks for a
             $VEC1, $VEC2, ... group for each monomer.  If they
             are found, they will be MOREAD.  If any of these are
             missing, the guess for that monomer will be constructed
             by HCORE.   Check your monomer energies carefully!  The
             initial guess orbitals for the supermolecule are formed
             by a block diagonal matrix of the monomer orbitals.
          2. The use of symmetry is turned off internally.
          3. There is no direct SCF option.  File ORDINT will be a
             full C1 list of integrals.  File AOINTS will contain
             whatever subset of these is needed for each particular
             decomposition step.  So extra disk space is needed
             compared to RUNTYP=ENERGY.
          4. This kind of run applies only to ab initio cases.
          5. This kind of run will work in parallel.
          6. Spherical harmonics may not be used.


          C.Coulson  in "Hydrogen Bonding", D.Hadzi, H.W.Thompson,
             Eds., Pergamon Press, NY, 1957, pp 339-360.
          C.Coulson  Research, 10, 149-159 (1957).
          K.Morokuma  J.Chem.Phys. 55, 1236-44 (1971).
          K.Kitaura, K.Morokuma  Int.J.Quantum Chem. 10, 325 (1976).
          K.Morokuma, K.Kitaura  in "Chemical Applications of
             Electrostatic Potentials", P.Politzer,D.G.Truhlar, Eds.
             Plenum Press, NY, 1981, pp 215-242.
          The method coded is the newer version described in the
          latter two papers.  Note that the CT term is computed
          separately for each monomer, as described in the words
          below equation 16 of the 1981 paper, not simultaneously.

          Reduced Variational Space:
          W.J.Stevens, W.H.Fink, Chem.Phys.Lett. 139, 15-22(1987).

          A comparison of the RVS and Morokuma decompositions can
          be found in the review article: "Wavefunctions and
          Chemical Bonding" M.S.Gordon, J.H.Jensen in "Encyclopedia
          of Computational Chemistry", volume 5, P.V.R.Schleyer,
          editor, John Wiley and Sons, Chichester, 1998.

          BSSE during Morokuma decomposition:
          R.Cammi, R.Bonaccorsi, J.Tomasi
          Theoret.Chim.Acta 68, 271-283(1985).

          The present implementation:
          "Energy decomposition analysis for many-body interactions,
           and application to water complexes"
          W.Chen, M.S.Gordon   J.Phys.Chem. 100, 14316-14328(1996)



$FFCALC group (relevant for RUNTYP=FFIELD)

              This group permits the study of the influence of an
          applied electric field on the wavefunction.  The most
          common finite field calculation applies a sequence of
          fields to extract the linear polarizability and first and
          second order hyperpolarizability.  The method is general,
          and so works for all ab initio wavefunctions in GAMESS.

          EFIELD      = applied electric field strength
                        (default=0.001 a.u.)

          IAXIS and JAXIS specify the orientation of the applied
                          field.  1,2,3 mean x,y,z respectively.
                          The default is IAXIS=3 and JAXIS=0.

            If IAXIS=i and JAXIS=0, the program computes alpha(ii),
            beta(iii), and gamma(iiii) from the energy changes, and
            a few more components from the dipole changes.  Five
            wavefunction evaluations are performed.

            If IAXIS=i and JAXIS=j, the program computes the cross
            terms beta(ijj), beta(iij), and gamma(iijj) from the
            energy changes, and a few more components from the
            dipole changes.  This requires nine evaluations of the

          AOFF        = a flag to permit evaluation of alpha(ij)
                        when the dipole moment is not available.
                        This is necessary only for MP2, and means
                        the off-axial calculation will do 13, not
                        9 energy evaluations.  Default=.FALSE.

          SYM         = a flag to specify when the fields to be
                        applied along the IAXIS and/or JAXIS (or
                        according to EONE below) do not break the
                        molecular symmetry.  Since most fields do
                        break symmetry, the default is .FALSE.

          ONEFLD      = a flag to specify a single applied field
                        calculation will be performed.  Only the
                        energy and dipole moment under this field
                        are computed.  If this option is selected,
                        only SYM and EONE input is heeded.  The
                        default is .FALSE.

          EONE        = an array of the three x,y,z components of
                        the single applied field.

          There are notes on RUNTYP=FFIELD on the next page.

              Finite field calculations require large basis sets,
          and extraordinary accuracy in the wavefunction.  To
          converge the SCF to many digits is sometimes problematic,
          but we suggest you use the input to increase integral
          accuracy and wavefunction convergence, for example

             $SCF    NCONV=10 FDIFF=.FALSE. $END

              In many cases, the applied fields will destroy the
          molecular symmetry.  This means the integrals are
          calculated once with point group symmetry to do the
          initial field free wavefunction evaluation, and then again
          with point group symmetry turned off.  If the fields
          applied do not destroy symmetry, you can avoid this second
          calculation of the integrals by SYM=.TRUE.  This option
          also permits use of symmetry during the applied field
          wavefunction evaluations.

              Examples of fields that do not break symmetry are a
          Z-axis field for an axial point group which is not
          centrosymmetric (i.e. C2v).  However, a second field in
          the X or Y direction does break the C2v symmetry.
          Application of a Z-axis field for benzene breaks D6h
          symmetry.  However, you could enter the group as C6v in
          $DATA while using D6h coordinates, and regain the prospect
          of using SYM=.TRUE.  If you wanted to go on to apply a
          second field for benzene in the X direction, you might
          want to enter Cs in $DATA, which will necessitate the
          input of two more carbon and hydrogen atom, but recovers
          use of SYM=.TRUE.

          Reference: H.A.Kurtz, J.J.P.Stewart, K.M.Dieter
                     J.Comput.Chem.  11, 82-87 (1990).

              For analytic computation of static and also frequency
          dependent NLO proerties, for closed shell cases, see the
          $TDHF group.



$TDHF group (relevant for SCFTYP=RHF if RUNTYP=TDHF)

              This group permits the analytic calculation of various
          static and/or frequency dependent polarizabilities, with
          an emphasis on important NLO properties such as second and
          third harmonic generation.  The method is programmed only
          for closed shell wavefunctions, at the semi-empirical or
          ab initio level.  Ab initio calculations may be direct SCF,
          or parallel, if desired.

              Because the Fock matrices computed during the time-
          dependent Hartree-Fock CPHF are not symmetric, you may not
          use symmetry.  You must enter NOSYM=1 in $CONTRL!

              For a more general numerical approach to the static
          properties, see $FFCALC.

          NFREQ  = Number of frequencies to be used. (default=1)

          FREQ   = An array of energy values in atomic units.  For
                   example: if NFREQ=3 then FREQ(1)=0.0,0.1,0.25.
                   By default, only the static polarizabilities are
                   computed.  (default is freq(1)=0.0)

              The conversion factor from Hartree to wave
              numbers is 219,474.6, and the wavelength is
              given (in nm) by 45.56/FREQ.

          MAXITA = Maximum number of iterations for an alpha
                   computation. (default=100)

          MAXITU = Maximum number of iterations in the second order
                   correction calculation.  This applies to iterative
                   beta values and all gammas. (default=100)

          ATOL   = Tolerance for convergence of first-order results.

          BTOL   = Tolerance for convergence of second-order results.

          RETDHF = a flag to choose starting points for iterative
                   calculations from best previous results.

          * * * the following NLO properties are available  * * *

          INIB   = 0 turns off all beta computation (default)
                 = 1 calculates only noniterative beta
                 = 2 calculate iterative and noniterative beta
                     The next flags allow further BETA tuning

          BSHG   = Calculate beta for second harmonic generation.

          BEOPE  = Calculate beta for electrooptic Pockels effect.

          BOR    = Calculate beta for optical rectification.

          INIG   = 0 turns off all gamma computation (default)
                 = 1 calculates only noniterative gamma
                 = 2 calculate iterative and noniterative gamma
                     The next flags allow further GAMMA tuning

          GTHG   = Calculate gamma for third harmonic generation.

          GEFISH = Calculate gamma for electric-field induced
                   second harmonic generation.

          GIDRI  = Calculate gamma for intensity dependent
                   refractive index.

          GOKE   = Calculate gamma for optical Kerr effect.

              These will be computed only if a nonzero energy is
          requested.  The default for each flag is .TRUE., and they
          may be turned off individually by setting some .FALSE.
          Note however that the program determines the best way to
          calculate them.  For example, if you wish to have the SHG
          results but no gamma results are needed, the SHG beta will
          be computed in a non-iterative way from alpha(w) and
          alpha(2w).  However if you request the computation of the
          THG gamma, the second order U(w,w) results are needed and
          an iterative SHG calculation will be performed whether
          you request it or not, as it is a required intermediate.

          S.P.Karna, M.Dupuis J.Comput.Chem.  12, 487-504 (1991).
          P.Korambath, H.A.Kurtz, in "Nonlinear Optical Materials",
          ACS Symposium Series 628, S.P.Karna and A.T.Yeates, Eds.
          pp 133-144, Washington DC, 1996.

          Review: D.P.Shelton, J.E.Rice, Chem.Rev. 94, 3-29(1994).



$EFRAG group (optional)

             This group gives the name and position of one or more
          effective fragment potentials.  It consists of a series of
          free format card images, which may not be combined onto a
          single line!  The position of a fragment is defined by
          giving any three points within the fragment, relative to
          the ab initio system defined in $DATA, since the effective
          fragments have a frozen internal geometry.  All other
          atoms within the fragment are defined by information in
          the $FRAGNAME group.


          -1-   a line containing one or more of these options:

               COORD   =CART     selects use of Cartesians coords
                                 to define the fragment position at
                                 line -3-.  (default)
                       =INT      selects use of Z-matrix internal
                                 coordinates at line -3-.
               POLMETHD=SCF      indicates the induced dipole for
                                 each fragment due to the ab initio
                                 electric field and other fragment
                                 fields is updated only once during
                                 each SCF iteration.
                       =FRGSCF   requests microiterations during
                                 each SCF iteration to make induced
                                 dipoles due to ab initio and other
                                 fragment fields self consistent
                                 amoung the fragments.  (default)
                                 Both methods converge to the same
                                 dipolar interaction.
               POSITION=OPTIMIZE Allows full optimization within the
                                 ab initio part, and optimization of
                                 the rotational and translational
                                 motions of each fragment. (default)
                       =FIXED    Allows full optimization of the
                                 ab initio system, but freezes the
                                 position of the fragments.  This
                                 makes sense only with two or more
                                 fragments, as what is frozen is the
                                 fragments' relative orientation.
                       =EFOPT    the same as OPTIMIZE, but if the
                                 fragment gradient is large, up to
                                 5 geometry steps in which only the
                                 fragments move may occur, before
                                 the geometry of the ab initio piece
                                 is relaxed.  This may save time by
                                 reusing the two electron integrals
                                 for the ab initio system.

               NBUFFMO = n       First n orbitals in the MO matrix
                                 are deemed to belong to the QM/MM
                                 buffer and will be excluded from
                                 the interaction with the EFP region.
                                 This makes sense only if these first
                                 MOs are frozen via the $MOFRZ group.

          Input a blank line if all the defaults are acceptable.

          -2-  FRAGNAME=XXX

          XXX is the name of the fragment whose coordinates are to
          be given next.  All other information defining the
          fragment is given in a supplemental $XXX group, which is
          referred to below as a $FRAGNAME group.

          A RHF/DZP EFP for water is internally stored in GAMESS.
          Choose FRAGNAME=H2OEF2 to look up this numerical data,
          and then skip the input of $H2OEF2 and $FRGRPL groups.

          -3-   NAME, X, Y, Z                           (COORD=CART)

          NAME     = the name of a fragment point.  The name used
                     here must match one of the points in $FRAGNAME.

          X, Y, Z  = Cartesian coordinates defining the position of
                     this fragment point RELATIVE TO THE COORDINATE
                     ORIGIN used in $DATA.  The choice of units is
                     controlled by UNITS in $CONTRL.

          I, DISTANCE, J, BEND, K, TORSION = the usual Z-matrix
                     connectivity internal coordinate definition.
                     The atoms I, J, K must be atoms in the ab
                     initio system from in $DATA, or fragment points
                     already defined in the current fragment or
                     previously defined fragments.

          Line -3- must be given a total of three times to define
          this fragment's position.

          Repeat lines -2- and -3- to enter as many fragments as you
          desire, and then end the group with a $END line.

          Note that it is quite typical to repeat the same fragment
          name at line -2-, to use the same fragment system at many
          different positions.



$FRAGNAME group (required for each FRAGNAME given in $EFRAG)

             This group gives all pertinent information for a given
          effective fragment potential (EFP).  This information
          falls into three categories:
               electrostatic (distributed multipoles, screening)
               distributed polarizabilities
               exchange repulsion
          It is input using several different subgroups, which
          should be given in the order shown below.  Each subgroup
          is specified by a particular name, and is terminated by
          the word STOP.  You may omit any of the subgroups to omit
          that term from the EFP.  All values are given in atomic

          To input monopoles,             follow input sequence -EM-
          To input dipoles,               follow input sequence -ED-
          To input quadrupoles,           follow input sequence -EQ-
          To input octupoles,             follow input sequence -EO-
          To input screening parameters,  follow input sequence -ES-
          To input polarizable points,    follow input sequence -P-
          To input repulsive points,      follow input sequence -R-


          -1-   a single descriptive title card

          -2-   COORDINATES

          COORDINATES signals the start of the subgroup containing
          the multipolar expansion terms (charges, dipoles, ...).
          Optionally, one can also give the coordinates of the
          polarizable points, or centers of exchange repulsion.

          -3-   NAME, X, Y, Z, WEIGHT, ZNUC

          NAME is a unique string identifying the point.
          X, Y, Z are the Cartesian coordinates of the point.
          WEIGHT and ZNUC are the atomic mass and nuclear charge,
          and are given only for the points which are nuclei.

          Typically the true nuclei will appear twice, once for
          defining the positive nuclear charge and its screening,
          and a second time for defining the electronic distributed

          Repeat line -3- for each expansion point, and terminate
          the list with a "STOP".

          -EM1-  MONOPOLES

          MONOPOLES signals the start of the subgroup containing
          the electronic and nuclear monopoles.

          -EM2-  NAME, CHARGE

          NAME must match one given in the COORDINATES subgroup.
          CHARGE = nuclear or electronic monopole at this point.

          Repeat -EM2- to define all desired charges.
          Terminate this subgroup with a "STOP".

          -ED1-  DIPOLES

          DIPOLES signals the start of the subgroup containing the
          dipolar part of the multipolar expansion.

          -ED2-  NAME, MUX, MUY, MUZ

          NAME must match one given in the COORDINATES subgroup.
          MUX, MUY, MUZ are the components of the electronic dipole.

          Repeat -ED2- to define all desired dipoles.
          Terminate this subgroup with a "STOP".

          -EQ1-  QUADRUPOLES

          QUADRUPOLES signals the start of the subgroup containing
          the quadrupolar part of the multipolar expansion.

          -EQ2-  NAME, XX, YY, ZZ, XY, XZ, YZ

          NAME must match one given in the COORDINATES subgroup.
          XX, YY, ZZ, XY, XZ, and YZ are the components of the
          electronic quadrupole moment.

          Repeat -EQ2- to define all desired quadrupoles.
          Terminate this subgroup with a "STOP".

          -EO1-  OCTUPOLES

          OCTUPOLES signals the start of the subgroup containing
          the octupolar part of the multipolar expansion.

          -EO2-  NAME, XXX, YYY, ZZZ, XXY, XXZ,
                       XYY, YYZ, XZZ, YZZ, XYZ

          NAME must match one given in the COORDINATES subgroup.
          XXX, ...  are the components of the electronic octupole.

          Repeat -EO2- to define all desired octupoles.
          Terminate this subgroup with a "STOP".

          -ES1-  SCREEN

          SCREEN signals the start of the subgroup containing the
          screening terms (A*exp[-B*r**2]) for the distributed
          multipoles, which account for charge penetration effects.

          -ES2-  NAME, A, B

          NAME must match one given in the COORDINATES subgroup.
          A, B are the parameters of the Gaussian screening term.

          Repeat -ES2- to define all desired screening points.
          Terminate this subgroup with a "STOP".


          POLARIZABLE POINTS signals the start of the subgroup
          containing the distributed polarizability tensors, and
          their coordinates.

          -P2-  NAME, X, Y, Z

          NAME gives a unique identifier to the location of this
          polarizability tensor.  It might match one of the points
          already defined in the COORDINATES subgroup, but often
          does not.  Typically the distributed polarizability
          tensors are located at the centroids of localized MOs.

          X, Y, Z are the coordinates of the polarizability point.
          They should be omitted if NAME did appear in COORDINATES.
          The units are controlled by UNITS= in $CONTRL.

          -P3-  XX, YY, ZZ, XY, XZ, YZ, YX, ZX, ZY

          XX, ... are components of the distributed polarizability,
          which is not a symmetric tensor.  XY means dMUx/dFy, where
          MUx is a dipole component, and Fy is a component of an
          applied field.

          Repeat -P2- and -P3- to define all desired polarizability
          tensors, and terminate this subgroup with a "STOP".


          REPULSIVE POTENTIAL signals the start of the subgroup
          containing the fitted exchange repulsion potential, for
          the interaction between the fragment and the ab initio
          part of the system.  This term also accounts for charge
          transfer effects.  The term has the form

                sum   C * exp[-D  * r**2]
                 i     i        i

          -R2-  NAME, X, Y, Z, N

          NAME may match one given in the COORDINATES subgroup,
          but need not.  If NAME does not match one of the
          known points, you must give its coordinates X, Y, and
          Z, otherwise omit these three values.  N is the total
          number of terms in the fitted repulsive potential.

          -R3-  C, D

          These two values define the i-th term in the repulsive
          potential.  Repeat line -R3- for all N terms.

          Repeat -R2- and -R3- to define all desired repulsive
          potentials,  and terminate this subgroup with a "STOP".


          The entire $FRAGNAME group is terminated by a " $END".



$FRGRPL group

          This group defines the inter-fragment repulsive potential,
          which consists primarily of exchange repulsions but also
          includes charge transfer.  Note that the functional form
          used for the fragment-fragment repulsion differs from
          that used for the ab initio-fragment repulsion, which is
          defined in the $FRAGNAME group.  The form of the potential
                sum   A * exp[-B * r]
                 i     i        i


          -1-  PAIR=FRAG1 FRAG2

          specifies which two fragment repulsions are being defined.
          $FRAGNAME input for the two names FRAG1 and FRAG2 must
          have been given.

          -2-  NAME1 NAME2 A B
               NAME1 NAME2 'EQ' NAME3 NAME4

          NAME1 must be one of the "NAME" points defined in the
          $FRAG1 group's COORDINATE section.  Similarly NAME2 must
          be a point from the $FRAG2 group.  In addition, NAME1 or
          NAME2 could be the keyword CENTER, indicating the center
          of mass of the fragment.

          A and B are the parameters of the fitted repulsive

          The second form of the input allows equal potential fits
          to be used.  The syntax implies that the potential between
          the points NAME1 and NAME2 should be taken the same as the
          potential previously given in this group for the pair of
          points NAME3 and NAME4.

          If there are NPT1 points in FRAG1, and NPT2 points in
          FRAG2, input line -2- should be repeated NPT1*NPT2 times.
          Terminate the pairs of potentials with a "STOP" card.
          Any pairs which you omit will be set to zero interaction.

          Typically the number of points on which fitted potentials
          might be taken to be all the nuclei in a fragment, plus
          the center of mass.

          Repeat lines -1- and -2- for all pairs of fragments, then
          terminate the group with a $END line.



$PCM group (optional)

             This group controls solvent effect computations using
          the Polarizable Continuum Method.  If this group is found
          in the input file, a PCM computation is performed.  The
          default calculation, chosen by selecting only the SOLVNT
          keyword, is to compute the electrostatic free energy.
          Appropriate numerical constants are provided for a wide
          range of solvents.  Additional keywords allow for more
          sophisticated computations, namely cavitation, repulsion,
          and dispersion free energies.  The methodology for these
          is general, but only numerical constants for water are
          provided.  There is additional information on PCM in the
          References chapter of this manual.

              PCM is programmed only for RHF and MCSCF wavefunctions.

              Geometry optimization with PCM will probably not be able
          to converge to the default OPTTOL in $STATPT due to some
          numerical inaccuracies.  You will likely need to raise this
          value if you attempt solvent optimizations, by a factor of
          two to five.

          --- the first set of parameters controls the computation:
                  IEF, ICOMP, ICAV, IDISP, IREP, IDP, and IFIELD.

          IEF      switch to choose boundary element method or the
                   integral equation formula as the PCM solver.
                 = 0 isotropic dielectrics using BEM PCM
                 = 1 anisotropic dielectrics using IEF PCM, see $IEFPCM
                 = 2 ionic solutions using IEF PCM, see $IEFPCM
                 = 3 isotropic dielectrics using IEF PCM (default)

             *** at the present time, there is a bug with IEF=1 or 2.

          ICOMP  = Compensation procedure for induced charges.
                   Gradient runs require ICOMP be 0 or 2 only.
                 = 0 No. (default)
                 = 1 Yes, each charge is corrected in proportion
                     to the area of the tessera to which it belongs.
                 = 2 Yes, using the same factor for all tesserae.
                 = 3 Yes, with explicit consideration of the
                     portion of solute electronic charge outside
                     the cavity, by the method of Mennucci and
                     Tomasi.  See the $NEWCAV group.

          The behaviour of PCM prior to Oct. 2000 can be recovered
          by selecting IEF=0 and ICOMP=2.  Options IEF=1 or 2 are
          incompatible with gradients and also must choose ICOMP=0.
          IEF=3 may not choose ICOMP=3, but if diffuse functions
          are in use, this default choice may benefit from ICOMP=2.
          The BEM method (IEF=0) should normally choose ICOMP=2.

          ICAV   = At the end of the run, calculate the cavitation
                   energy, by the method of Pierotti and Claverie:
                 = 0 skip the computation (default)
                 = 1 perform the computation.

            If ICAV=1, the following parameter is relevant:

          TABS   = the absolute temperature, in units K.

               There are two procedures for the calculation
               of the repulsion and dispersion free energy.
               IDISP is incompatible with IREP and IDP.

          IDISP  = Calculation of both dispersion and repulsion
                   free energy through the empirical method of
                   Floris and Tomasi.
                 = 0 skip the computation (default)
                 = 1 perform the computation.  See $DISREP group.

            The next two options add repulsive and dispersive terms
            to the solute hamiltonian, in an ab initio manner, by
            the method of Amovilli and Mennucci.

          IREP   = Calculation of repulsion free energy
                 = 0 skip the computation (default)
                 = 1 perform the computation.  See $NEWCAV group.

          IDP    = Calculation of dispersion free energy
                 = 0 skip the computation (default)
                 = 1 perform the computation.  See $DISBS group.

            If IDP=1, then three additional parameters must be
            defined.  The two solvent values correspond to water,
            and therefore these must be input for other solvents.

          WA     = solute average transition energy.  This is
                   computed from the orbital energies for RHF,
                   but must be input for MCSCF runs.
          WB     = ionization potential of solvent, in Hartrees.
          ETA2   = square of the zero frequency refractive index
                   of the solvent.  (default=1.75)

          IFIELD = At run end, calculate the electric potential
                   and electric field generated by the apparent
                   surface charges.
                 = 0 skip the computation (default)
                 = 1 on nuclei
                 = 2 on a planar grid

            If IFIELD=2, the following data must be input:

          AXYZ,BXYZ,CXYZ = each defines three components of the
                           vertices of the plane where the reaction
                           field is to be computed (in Angstroms)
                A ===> higher left corner of the grid
                B ===> lower left corner of the grid
                C ===> higher right corner of the grid
          NAB = vertical subdivision (A--B edge) of the grid
          NAC = horizontal subdivision (A--C edge) of the grid.

          IPRINT = 0 normal printing (default)
                 = 1 turns on debugging printout

          --- the next group of keywords defines the solvent

          SOLVNT = keyword naming the solvent of choice.  The eight
                   numerical constants defining the solvent are
                   internally stored for the following:
                       WATER (or H2O)
                       CH3OH                      C2H5OH
                       CLFORM (or CHCl3)          CTCL (or CCl4)
                       METHYCL (or CH2Cl2)        12DCLET (or C2H4Cl2)
                       BENZENE (or C6H6)          TOLUENE (or C6H5CH3)
                       CLBENZ (or C6H5Cl)         NITMET (or CH3NO2)
                       NEPTANE (or C7H16)         CYCHEX (or C6H12)
                       ANILINE (or C6H5NH2)       ACETONE (or CH3COCH3)
                       THF                        DMSO (or DMETSOX)
                   The default solvent name is
                   which indicates you will specify your solvent by
                   giving the following 8 numerical values:

          RSOLV  = the solvent radius, in units Angstrom
          EPS    = the dielectric constant
          EPSINF = the dielectric constant at infinite frequency.
                   This value must be given only for RUNTYP=TDHF,
                   if the external field frequency is in the optical
                   range and the solvent is polar; in this case the
                   solvent response is described by the electronic
                   part of its polarization.  Hence the value of the
                   dielectric constant to be used is that evaluated
                   at infinite frequency, not the static one (EPS).
                   For nonpolar solvents, the difference between
                   the two is almost negligible.
          TCE    = the thermal expansion coefficient, in units 1/K
          VMOL   = the molar volume, in units ml/mole
          STEN   = the surface tension, in units dyne/cm
          DSTEN  = the thermal coefficient of log(STEN)
          CMF    = the cavity microscopic coefficient

          Values for TCE, VMOL, STEN, DSTEN, CMF need to be given
          only for the case ICAV=1.  Input of any or all of these
          values will override the internally stored value.

          --- the next set of keywords defines the molecular cavity

          NESFP  = the number of initial spheres.
                   (default = number of atoms in solute molecule)

          ICENT  = option for definition of initial spheres.
                 = 0 centers spheres on each nucleus.  (default)
                 = 1 sphere centers XE, YE, ZE and radii RIN will be
                     specified explicitly in $PCMCAV.

             The cavity generation algorithm may use additional
             spheres to smooth out sharp grooves, etc.  The
             following parameters control how many extra spheres
             are generated:

          OMEGA and FRO = GEPOL parameters for the creation of the
                   `added spheres' defining the solvent accessible
                   surface. When an excessive number of spheres is
                   created, which may cause problems of convergence,
                   the value of OMEGA and/or FRO must be increased.
                   For example, OMEGA from 40 to 50 ... up to 90,
                                FRO from 0.2 ... up to 0.7.
                   (defaults are OMEGA=40.0, FRO=0.7)

          RET    = minimum radius (in A) of the added spheres.
                   Increasing RET decreases the number of added
                   spheres.  A value of 100.0 inhibits the addition
                   of spheres.  (default=0.2)



$PCMCAV group (optional)

             This group controls generation of the cavity holding
          the solute during Polarizable Continuum Method runs.
          The cavity is a union of spheres, according to ICENT and
          associated input values given in $PCM.  The data given
          here must be given in Angstrom units.

          XE,YE,ZE = arrays giving the coordinates of the spheres.
              if ICENT=0, the atomic positions will be used.
              if ICENT=1, you must supply NESFP values here.

          RIN = an array giving the sphere radii.
              if ICENT=0, the program will look up the internally
                          stored van der Waals radius for:  H,He,
                             B,C,N,O,F,Ne,   Na,Al,Si,P,S,Cl,Ar,
                             K,As,Se,Br,Kr,  Rb,Sb,Te,I,  Cs,Bi
                          Data for other elements is not tabulated.
              if ICENT=1, give NESFP values.

          ALPHA = an array of scaling factors, for the definition of
                  the solvent accessible surface.  If only the first
                  value is given, all radii are scaled by the same
                  factor.  (default is ALPHA(1)=1.2)

          Example: Suppose the 4th atom in your molecule is Fe, but
                   all other atoms have van der Waals radii.  You
                   decide a good guess for Fe is twice the covalent
                   radius:  $PCMCAV RIN(4)=2.33 $END

          The source for the van der Waals radii is "The Elements",
          2nd Ed., John Emsley, Clarendon Press, Oxford, 1991,
          except that for C,N,O, the U.Pisa's experience with the
          best radii for PCM treatment of singly bonded C,N,O atoms
          is used instead.  The radii for a few transition metals
          are given by A.Bondi, J.Phys.Chem. 68, 441-451(1968).



$NEWCAV group (optional)

             This group controls generation of the "escaped charge"
          cavity, used when ICOMP=3 or IREP=1 in $PCM.  This cavity
          is used only to calculate the fraction of the solute
          electronic charge escapes from the original cavity.

          IPTYPE = choice for tessalation of the cavity's spheres.
                 = 1 uses a tetrahedron
                 = 2 uses a pentakisdodecahedron (default)

          ITSNUM = m, the number of tessera to use on each sphere.
                 if IPTYPE=1, input m=30*(n**2), with n=1,2,3 or 4
                 if IPTYPE=2, input m=60*(n**2), with n=1,2,3 or 4
                 (default is 60)

             *** the next three parameters pertain to IREP=1 ***

          RHOW   = density, relative to liquid water (default = 1.0)

          PM     = molecular weight (default = 18.0)

          NEVAL  = number of valence electrons on solute (default=8)

          The defaults for RHOW, PM, and NEVAL correspond to water,
          and therefore must be correctly input for other solvents.



$IEFPCM group (optional)

              This group defines data for the integral equation
          formalism version of PCM solvation.  It includes special
          options for ionic or anisotropic solutions.

          The next two sets are relevant only for anisotropic
          solvents, namely IEF=1:

          EPS1, EPS2, EPS3 =
                  diagonal values of the dielectric permittivity
                  tensor with respect to the laboratory frame.
                  The default is EPS in $PCM

          EUPHI, EUTHE, EUPSI =
                  Eulerian angles which give the rotation of the
                  solvent orientation with respect to the lab frame.
                  The term lab frame means $DATA orientation.
                  The default for each is zero degrees.

          The next two are relevant to ionic solvents, namely IEF=2:

          EPSI = the ionic solutions's dielectric, the default is
                 EPS from $PCM.

          DISM = the ionic strength, in Molar units (mol/dm**3)
                 The default is 0.0



$DISBS group (optional)

             This group defines auxiliary basis functions used to
          evaluate the dispersion free energy by the method of
          Amovilli and Mennucci.  These functions are used only for
          the dispersion calculation, and thus have nothing to do
          with the normal basis given in $BASIS or $DATA.  If the
          input group is omitted, only the normal basis is used for
          the IDP=1 dispersion energy.

          NADD   = the number of added shells

          XYZE   = an array giving the x,y,z coordinates (in bohr)
                   of the center, and exponent of the added shell,
                   for each of the NADD shells.

          NKTYPE = an array giving the angular momenta of the shells

          An example placing 2s,2p,2d,1f on one particular atom,

           $DISBS  NADD=7 NKTYP(1)= 0 0 1 1 2 2 3
                   XYZE(1)=2.9281086   0.0  .0001726   0.2
                           2.9281086   0.0  .0001726   0.05
                           2.9281086   0.0  .0001726   0.2
                           2.9281086   0.0  .0001726   0.05
                           2.9281086   0.0  .0001726   0.75
                           2.9281086   0.0  .0001726   0.2
                           2.9281086   0.0  .0001726   0.2  $END



$DISREP group (optional)

             This group controls evaluation of the dispersion and
          repulsion energies by the empirical method of Floris and
          Tomasi.  The group must be given with IDISP=1 in $PCM.
          The two options are controlled by ICLAV and ILJ, only one
          of which should be selected.

          ICLAV = selects Claverie's disp-rep formalism.
                = 0 skip computation.
                = 1 Compute the solute-solvent disp-rep interaction
                    as a sum over atom-atom interactions through a
                    Buckingham-type formula (R^-6 for dispersion,
                    exp for repulsion).  (default)
                    Ref: Pertsin-Kitaigorodsky "The atom-atom
                         potential method", page 146.

          ILJ   = selects a Lennard-Jones formalism.
                = 0 skip computation. (default)
                = 1 solute atom's-solvent molecule interaction is
                    modeled by Lennard-Jones type potentials, R^-6
                    for dispersion, R^-12 for repulsion).

          ---- the following data must given for ICLAV=1:

          RHO   = solvent numeral density
          N     = number of atom types in the solvent molecule
          NT    = an array of the number of atoms of each type in a
                  solvent molecule
          RDIFF = distances between the first atoms of each type
                  and the cavity
          DK,RW = parameters of atom-atom interactions

          The defaults are chosen for water,

          ---- the following data must given for ILJ=1:

          RHO   = solvent numeral density
          EPSI  = an array of energy constants referred to each atom
                  of the solute molecule.
          SIGMA = an array of typical distances, relative to each
                  solute atom



$COSGMS group (optional)

              The presence of this group in the input turns on the
          use of the conductor-like screening model with molecular
          shaped cavity for RHF and closed shell MP2.  For RHF, the
          energy and gradient can be computed, while MP2 is limited
          to the energy only.

          EPSI   = the dielectric constant, 80 is often used for H2O
                   This parameter must be given.

          RSOLV  = the multiplicative factor for the van der Waals
                   radius used for cavity construction.

          NSPA   = the number of surface points on each atomic
                   sphere that form the cavity.  (default=92)

               Additional information on the COSMO model can be
                found in the References chapter of this manual.



$SCRF group (optional)

              The presence of this group in the input turns on the
          use of the Kirkwood-Onsager spherical cavity model for the
          study of solvent effects.  The method is implemented for
          RHF, UHF, ROHF, GVB and MCSCF wavefunctions and gradients,
          and so can be used with any RUNTYP involving the gradient.
          The method is not implemented for MP2, CI, any of the
          semiempirical models, or for analytic hessians.

          DIELEC = the dielectric constant, 80 is often used for H2O

          RADIUS = the spherical cavity radius, in Angstroms

          G      = the proportionality constant relating the solute
                   molecule's dipole to the strength of the reaction
                   field.  Since G can be calculated from DIELEC and
                   RADIUS, do not give G if they were given.

               Additional information on the SCRF model can be
               found in the References chapter of this manual.



$ECP group (required if ECP=READ in $CONTRL)

              This group lets you read in effective core potentials,
          for some or all of the atoms in the molecule.  You can
          use built in potentials for some of the atoms if you like.
          This is a free format (positional) input group.

          *** Give a card set -1-, -2-, and -3- for each atom ***

          -card 1-    PNAME, PTYPE, IZCORE, LMAX+1

          PNAME is a 8 character descriptive tag for this potential.
                If it is repeated for a subsequent atom, no other
                information need be given on this card, and cards
                -2- and -3- may also be skipped.  The information
                will be copied from the first atom by this PNAME.
                Do not use the option to repeat the previously read
                ECP for an atom with PTYPE=NONE, instead type "none".
          PTYPE = GEN    a general potential should be read.
                = SBKJC  look up the Stevens/Basch/Krauss/Jasien/
                         Cundari potential for this type of atom.
                = HW     look up the Hay/Wadt built in potential
                         for this type of atom.
                = NONE   treat all electrons on this atom.
          IZCORE is the number of core electrons to be removed.
          LMAX   is the maximum angular momentum occupied in the
                 core orbitals being removed (usually).  Give
                 IZCORE and LMAX only if PTYPE is GEN.

          *** For the first occurence of PNAME, if PTYPE is GEN, ***
          *** then give cards -2- and -3-.  Otherwise go to -1-. ***

          *** Card sets -2- and -3- are repeated LMAX+1 times    ***

              The potential U(LMAX+1) is given first,
              followed by U(L)-U(LMAX+1), for L=1,LMAX.

          -card 2-    NGPOT

          NGPOT is the number of Gaussians in this part of the
                local effective potential.

          -card 3-    CLP,NLP,ZLP   (repeat this card NGPOT times)

          CLP is the coefficient of this Gaussian in the potential.
          NLP is the power of r for this Gaussian.
          ZLP is the exponent of this Gaussian.

                                * * *

          By far the easiest way to use the SBKJC potential for all
          atoms in the formic acid molecule is to request ECP=SBKJC
          in $CONTRL.  But the next page shows two alternatives.

          The first way is to look up the program's internally
          stored SBKJC potentials one atom at a time:

          C-ECP SBKJC
          H-ECP NONE
          O-ECP SBKJC
          H-ECP NONE

          The second oxygen duplicates the first, no core electrons
          are removed for hydrogen.  The order of the atoms must
          follow that generated by $DATA.  Note PTYPE allows you to
          type in one or more atoms explicitly, while using built in
          data for some other atoms.

          The second example reads all SBKJC potentials explicitly:

          C-ECP GEN 2 1
          1      ----- CARBON U(P) -----
           -0.89371  1  8.56468
          2      ----- CARBON U(S)-U(P) -----
            1.92926  0  2.81497
           14.88199  2  8.11296
          H-ECP NONE
          O-ECP GEN 2 1
          1      ----- OXYGEN U(P) -----
           -0.92550  1 16.11718
          2      ----- OXYGEN U(S)-U(P) -----
            1.96069  0  5.05348
           29.13442  2 15.95333
          H-ECP NONE

          Again, the 2nd oxygen copies from the first.  It is handy
          to use the rest of card -2- as a descriptive comment.

          As a final example, for antimony we have LMAX+1=3 (there
          are core d's).  One must first enter U(f), followed by
          U(s)-U(f), U(p)-U(f), U(d)-U(f).

          At the present time, there are some numerical problems in
          the ECP code, which has been programed to use spdfg basis
          sets, and core potentials up to g.  If one is using a g
          basis function, or a g potential (bottom row elements
          beyond the lanthanide series), there are small errors in
          the ECP integrals.  A tight optimization (OPTTOL=1D-05)
          will usually result in the energy rising slightly during
          the last few geometry steps.  The error seems to be about
          0.000002 Hartree, so be cautious about using a g potential
          or basis function.  When both are used in the same run,
          the error in the energy is about 0.0002, which means the
          use of both is too inaccurate to be trusted.  The use of
          f functions or f potentials (or lower) seems to be free of
          any errors.



$RELWFN group (optional)

              This group is relevant if RELWFN in $CONTRL chose the
          NESC or RESC option for elimination of small components
          from relativistic wavefunctions, to produce a corrected
          single component wavefunction.  In case of RESC, only the
          one electron integral corrections are added, whereas for
          NESC, corrections to two electron integrals are accounted
          for by means of a relativistically averaged basis set.

              For NESC, you must provide three basis sets, for the
          large and small components and an averaged one, which are
          given in $DATAL, $DATAS, $DATA, respectively.  The only
          possible choice for these basis sets is due to Dyall, and
          these are available from
          Their names are similar to cc-pVnZ(pt/sf/lc), pt=point or
          fi=finite nucleus, sf for spin-free and the final field is
          lc=large component ($DATAL), sc=small component ($DATAS),
          and wf is a typo for Foldy-Wouthuysen 2e- basis ($DATA).
          In GAMESS you can only use point nucleus approximation.
          The need to input three basis sets means that you cannot
          use a $BASIS group, and you must use COORD=UNIQUE style
          input in the various $DATA's.  The three $DATA groups must
          contain identical information except for the primitive
          expansion coefficients, as the three basis sets must have
          the same exponents.  In case the option to treat only some
          atoms relativistically is chosen, all non-relativistic
          atoms must have identical basis input in all three groups.

              For RESC, ordinary basis sets are used.

              Analytic gradients are programmed for both RESC and
          NESC computations.  For NESC, the one electron part of
          the spin-orbit operator can be corrected, while for RESC,
          one can compute spin-orbit coupling with relativistic
          corrections to both one and two electron SOC integrals.

          NESOC  =   relativistic corrections to SOC integrals.
                     Choose only if RELWFN=RESC or NESC, and if
                     OPERAT=HSO1, HSO2P, or HSO2, for RUNTYP=TRANSITN
                 = 0 no corrections
                 = 1 one-electron spin-orbit integrals (NESC default)
                 = 2 one and two-electron integrals (RESC default)

          The remaining parameters pertain only to RELWFN=NESC:

          NRATOM the number of different elements to be treated
                 nonrelativistically.  For example, in Pb3O4, to
                 treat only lead relativistically, enter NRATOM=1.

          CHARGE array containing charges of atoms to be treated
                 nonrelativistically.  (e.g. CHARGE(1)=8.0, to drop
                 all oxygen atoms)



$EFIELD group (not required)

              This group permits the study of the influence of an
          external electric field on the molecule.  The method is
          general, and so works for all ab initio SCFTYPs.

          EVEC        = an array of the three x,y,z components of
                        the applied electric field.

          SYM         = a flag to specify when the field to be
                        applied breaks the molecular symmetry.
                        Since most fields break symmetry, the
                        default is .FALSE.


          Restrictions: analytic hessians are not available, but
          numerical hessians are.  Because an external field causes
          a molecule with a dipole to experience a torque, geometry
          optimizations must be done in Cartesian coordinates only.
          Internal coordinates eliminate the rotational degrees of
          freedom, which are no longer free.

          Notes: a hessian calculation will have two rotational
          modes with non-zero "frequency", caused by the torque.
          A gas phase molecule will rotate so that the dipole
          moment is anti-parallel to the applied field.  To carry
          out this rotation during geometry optimization will take
          many steps, and you can help save much time by inputting
          a field opposite the molecular dipole.  There is also
          a stationary point at higher energy with the dipole
          parallel to the field, which will have two imaginary
          frequencies in the hessian.  Careful, these will appear
          as the first two modes in a hessian run, but will not
          have the i for imaginary included on the printout since
          they are rotational modes.



$INTGRL group (optional)

              This group controls AO integral formats.  It should
          probably never be given, as the program always picks
          sensible values.

           SCHWRZ = a flag to activate use of the Schwarz inequality
                    to predetermine small integrals.  There is no
                    loss of accuracy when choosing this option, and
                    there are appreciable time savings for bigger
                    molecules.  Default=.TRUE. for over 5 atoms, or
                    for direct SCF, and is .FALSE. otherwise.

           NOPK   = 0 PK integral option on, which is permissible
                      for RHF, UHF, ROHF, GVB energy/gradient runs.
                  = 1 PK option off (default for all jobs).
                      Must be off for anything with a transformation.

           NORDER = 0 (default)
                  = 1 Sort integrals into canonical order.  There
                      is little point in selecting this option, as
                      no part of GAMESS requires ordered integrals.
                      See also NSQUAR.

           NINTMX =   Maximum no. of integrals in a record block.
                      (default=15000 for J or P file, =10000 for PK)

                The following parameters control the integral sort.
                (values given are defaults)

           NSQUAR = 0 Sorted integrals will be in triangular
                      canonical order (default)
                  = 1 instead sort to square canonical order.
           NDAR   = Number of direct access logical records to be
                    used for the integral sort (default=2000)
           LDAR   = Length of direct access records (site dependent)
           NBOXMX =  200   Maximum number of bins.
           NWORD  =    0   Memory to be used (default=all of it).
           NOMEM  =    0   If non-zero, force external sort.

                The following parameters control integral restarts
                (values given are defaults)
           IST=    1      JST=    1    KST=    1    LST=    1
           NREC=   1      INTLOC= 1



$TRANS group (optional for -CI- or -MCSCF-, relevant to analytic hessians & relevant to energy localization)

               This group controls the integral tranformation.  MP2
          integral transformations are controlled instead by the
          $MP2 input group.  There is little reason to give any but
          the first variable.

           DIRTRF = a flag to recompute AO integrals rather than
                    storing them on disk.  The default is .FALSE.
                    for MCSCF and CI runs.  If your job reads $SCF,
                    and you select DIRSCF=.TRUE. in that group, a
                    direct transformation will be done, no matter
                    how DIRTRF is set.

              Note that the transformation may do many passes over
              the AO integrals for large basis sets, and thus the
              direct recomputation of AO integrals can be very time

           MPTRAN = method to use for the integral transformation.
                    the default is try 0, then 1, then 2.
                    0 means use the incore method
                    1 means use the segmented method.  This is the
                      only method that works in parallel.
                    2 means use the alternate method, which uses
                      less memory than 2, but requires an extra
                      large disk file.

           NWORD  = Number of words of fast memory to allow.  Zero
                    uses all available memory. (default=0)

           CUTTRF = Threshold cutoff for keeping transformed two
                    electron integrals.  (default= 10**(-9))

           AOINTS = defines AO integral storage during conventional
                    integral transformations, during parallel runs.
                    DUP stores duplicated AO lists on each node, and
                    is the default for parallel computers with slow
                    interprocessor communication, e.g. ethernet.
                    DIST distributes the AO integral file across
                    all nodes, and it is the default for parallel
                    computers with high speed communications.



$CIINP group (optional, relevant for CITYP=GUGA or ALDET)

              This group is the control box for Graphical Unitary
          Group Approach (GUGA) CI calculations, or Ames Laboratory
          determinant (ALDET) full CI.  Each step which is executed
          potentially requires a further input group described later.

          NRNFG = An array of 10 switches controlling which steps of
                  a CI computation are performed.
                  1 means execute the module, 0 means don't.

            NRNFG(1) = Generate the configurations.  See either
                       $CIDRT or $CIDET input.  (default=1)
            NRNFG(2) = Transform the integrals. See $TRANS.
            NRNFG(3) = Sort integrals and calculate the Hamiltonian
                       matrix. See $CISORT and $GUGEM. (default=1)
                       This does not apply to ALDET.
            NRNFG(4) = Diagonalize the Hamiltonian matrix.
                       See $GUGDIA or $CIDET. (default=1)
            NRNFG(5) = Construct the one electron density matrix,
                       and generate NO's. See $GUGDM or $CIDET.
            NRNFG(6) = Construct the two electron density matrix.
                       See $GUGDM2 or $CIDET.
                       (default=0 normally, but 1 for CI gradients)
            NRNFG(7) = Construct the Lagrangian of the CI function.
                       Requires DM2 matrix exists.  See $LAGRAN.
                       (default=0 normally, but 1 for CI gradients)
                       This does not apply to ALDET.
            NRNFG(8-10) are not used.

          Users are not encouraged to change these values, as the
          defaults are quite reasonable ones.

          NPFLG = An array of 10 switches to produce debug printout.
                  There is a one to one correspondance to NRNFG, set
                  to 1 for output. (default = 0,0,0,0,0,0,0,0,0,0)
                  The most interesting is NPFLG(2)=1 to see the
                  transformed 1e- integrals, NPFLG(2)=2 adds the
                  very numerous transformed 2e- integrals to this.

          IREST = n    Restart the -CI- at stage NRNFG(n).


$DET group (required for SCFTYP=MCSCF if CISTEP=ALDET)

$CIDET group (required if CITYP=ALDET)

             This group describes the determinants to be used in a
          full MCSCF active space, or full CI wavefunction.

             Determinants contain several spin states, in contrast
          to configuration state functions.  The Sz quantum number
          of each determinant is the same, but the Hamiltonian
          eigenvectors will have various spins S=Sz, Sz+1, Sz+2, ...
          so NSTATE may need to account for state of other spin
          symmetry.  In Abelian groups, you can specify the exact
          spatial symmetry you desire.

             The first two determine the symmetry of the states:

          GROUP  = name of the point group.  The default is to copy
                   this from $DATA, if that group is Abelian (C2,
                   Ci, Cs, C2v, C2h, D2, or D2h).  If not, the
                   group is set to C1 (no symmetry used).

          ISTSYM = specifies the spatial symmetry of the state.
                   This table is exactly the same as in $DRT input.
                     ISTSYM= 1   2   3   4   5   6   7   8
                        C1   A
                        Ci   Ag  Au
                        Cs   A'  A''
                        C2   A   B
                        C2v  A1  A2  B1  B2
                        C2h  Ag  Bu  Bg  Au  <- differs from $MCQDPT!
                        D2   A   B1  B2  B3
                        D2h  Ag  B1g B2g B3g Au  B1u B2u B3u
                   Default is ISTSYM=1, the totally symmetric state.

             The next four define the filled and active orbital space.
          There is no default for NCORE, NACT, and NELS:

          NCORE  = total number of orbitals doubly occupied in all

          NACT   = total number of active orbitals.

          NELS   = total number of active electrons.

          SZ     = azimuthal spin quantum number for each of the
                   determinants, two times SZ is therefore the
                   number of excess alpha spins in each determinant.
                   The default is SZ=S, extracted from the MULT=2S+1
                   given in $CONTRL.

           * * * the following control the diagonalization * * *

          NSTATE = Number of CI states to be found, the default is
                   1.  The maximum number of states is 100.

          PRTTOL = Printout tolerance for CI coefficients, the
                   default is to print any larger than 0.05.

          ITERMX = Maximum number of Davidson iterations per root.
                   The default is 100.

          CVGTOL = Convergence criterion for Davidson eigenvector
                   routine.  This value is proportional to the
                   accuracy of the coeficients of the eigenvectors
                   found.  The energy accuracy is proportional to
                   its square.  The default is 1.0E-5.

          NHGSS  = dimension of the Hamiltonian submatrix which
                   is diagonalized to obtain the initial guess
                   eigenvectors.  The determinants forming the
                   submatrix are chosen on the basis of a low
                   diagonal energy, or if needed to complete a
                   spin eigenfunction.  The default is 300.

          NSTGSS = Number of eigenvectors from the initial guess
                   Hamiltonian to be included in the Davidson's
                   iterative scheme.  It is seldom necessary to
                   include extra states to obtain convergence to
                   the desired states.  The default is the value
                   for NSTATE.

          MXXPAN = Maximum number of expansion basis vectors in the
                   iterative subspace during the Davidson iterations
                   before the expansion basis is truncated.  The
                   default is the larger of 10 or 2*NSTGSS.  Larger
                   values might help convergence, do not decrease
                   this parameter below 2*NSTGSS.

          * * * the following control the 1st order density * * *
          These are ignored during MCSCF, but are used during a CI.

          IROOT  = the root whose density is saved on the disk file
                   for subsequent property analysis.  Only one root
                   can be saved, and the default value of 1 means
                   the ground state.  Be sure to set NFLGDM to form
                   the density of the state you are interested in!

          NFLGDM = Controls each state's density formation.
                   0 -> do not form density for this state.
                   1 -> form density and natural orbitals for this
                        state, print and punch occ.nums. and NOs.
                   2 -> same as 1, plus print density over MOs.
                   The default is NFLGDM(1)=1,0,0,...,0 meaning
                   only ground state NOs are generated.

              * * * the following control the state averaged
              * * * 1st and 2nd order density matrix computation
          Usually ignored by CI runs, these are relevant to MCSCF.

          PURES  = a flag controlling the spin purity of the state
                   avaraging.  If true, the WSTATE array pertains
                   to the lowest states of the same S value as is
                   given by the MULT keyword in $CONTRL.  In this
                   case the value of NSTATE will need to be bigger
                   than the total number of weights given by WSTATE
                   if there are other spin states present at low
                   energies.  If false, it is possible to state
                   average over more than one S value, which might
                   be of interest in spin-orbit coupling jobs.
                   The default is .TRUE.

          WSTATE = An array of up to 100 weights to be given to the
	   densities of each state in forming the average.
	   The default is to optimize a pure ground state,
	   A small amount of the ground state can help the
	   convergence of excited states greatly.
	   Gradient runs are possible only with pure states.
	   Be sure to set NSTATE above appropriately!



$DRT group (required for SCFTYP=MCSCF if CISTEP=GUGA)

$CIDRT group (required if CITYP=GUGA)

              This group describes the -MCSCF- or -CI- wavefunction.
          The distinct row table is the means by which the Graphical
          Unitary Group Approach (GUGA) names the configurations.

             The group is spelled DRT for MCSCF runs, and CIDRT for
          CI runs.  The main difference in these is NMCC vs. NFZC.

              There is no default for GROUP, and you must choose one
          of FORS, FOCI, SOCI, or IEXCIT.

          GROUP = the name of the point group to be used.  This is
                  usually the same as that in $DATA, except for
                  RUNTYP=HESSIAN, when it must be C1.  Choose from
                  the following: C1, C2, CI, CS, C2V, C2H, D2, D2H,
                  C4V, D4, D4H.  If your $DATA group is not listed,
                  choose only C1 here.

          FORS  = flag specifying the Full Optimized Reaction Space
                  set of configuration should be generated.  This
                  is usually set true for MCSCF runs, but if it is
                  not, see FORS in $MCSCF.  (Default=.FALSE.)

          FOCI  = flag specifying first order CI.  In addition to
                  the FORS configurations, all singly excited CSFs
                  from the FORS reference are included.

          SOCI  = flag specifying second order CI.  In addition to
                  the FORS configurations, all singly and doubly
                  excited configurations from the FORS reference
                  are included.  (Default=.FALSE.)

          IEXCIT= electron excitation level, for example 2 will
                  lead to a singles and doubles CI.  This variable
                  is computed by the program if FORS, FOCI, or
                  SOCI is chosen, otherwise it must be entered.

            * * the next variables define the single reference * *

              The single configuration reference is defined by
          filling in the orbitals by each type, in the order shown.
          The default for each type is 0.

                 Core orbitals, which are always doubly occupied:
          NMCC = number of MCSCF core MOs (in $DRT only).
          NFZC = number of CI frozen core MOs (in $CIDRT only).

                 Internal orbitals, which are partially occupied:
          NDOC = number of doubly occupied MOs in the reference.
          NAOS = number of alpha occupied MOs in the reference,
                 which are singlet coupled with a corresponding
                 number of NBOS orbitals.
          NBOS = number of beta spin singly occupied MOs.
          NALP = number of alpha spin singly occupied MOs in the
                 reference, which are coupled high spin.
          NVAL = number of empty MOs in the reference.

                 External orbitals, occupied only in FOCI or SOCI:
          NEXT = number of external MOs.  If given as -1, this will
                 be set to all remaining orbitals (apart from any
                 frozen virtual orbitals).
          NFZV = number of frozen virtual MOs, never occupied.

                  * * the next two help with state symmetry * *

          ISTSYM= irreducible representation for GUGA wavefunction.
                  This option overwrites whatever symmetry is implied
                  by NALP/NAOS/NBOS.  Default=0 means the symmetry of
                  the reference configuration is to be used.
                     ISTSYM= 1   2   3   4   5   6   7   8
                        C1   A
                        Ci   Ag  Au
                        Cs   A'  A''
                        C2   A   B
                        C2v  A1  A2  B1  B2
                        C2h  Ag  Bu  Bg  Au  <- differs from $MCQDPT!
                        D2   A   B1  B2  B3
                        D2h  Ag  B1g B2g B3g Au  B1u B2u B3u

          NOIRR= controls labelling of the CI state symmetries.
               = 1 no labelling (default)
               = 0 usual labelling.  This can be very time consuming
                   if the group is non-Abelian.
               =-1 fast labelling, in which all CSFs with small CI
                   coefficients are ignored.  This can produce weights
                   quite different from one, due to ignoring the small
                   coefficients, but overall seems to work OK.
                   Note that it is normal for the weights not to sum
                   to 1 even for NOIRR=0 because for simplicity the
                   weight determination is focused on the relative
                   weights rather than absolute.  However weight do
                   not sum to one only for row-mixed MOs.
               = -2,-3... fast labelling and sets SYMTOL=10**NOIRR
                   for runs other than TRANSITN.  All irreps with
                   weights greater than SYMTOL are considered.

                 * * * the final choices are seldom used * * *

          INTACT= flag to select the interacting space option.
                  The CI will include only those spin couplings
                  which have a nonvanishing matrix element with
                  the reference configuration.

          MXNINT = Buffer size for sorted integrals. (default=20000)

          MXNEME = Buffer size for energy matrix.  (default=10000)

          NPRT   = Configuration printout control switch.
                   This can consume a HUMUNGUS amount of paper!
                   0 = no print (default)
                   1 = print electron occupancies, one per line.
                   2 = print determinants in each CSF.
                   3 = print determinants in each CSF (for Ms=S-1).



$MCSCF group (optional for -MCSCF-)

              This group controls the MCSCF orbital optimization
          step.  The difference between the four convergence methods
          is outlined in Chapter Four of this manual, which you must
          carefully study before attempting MCSCF computations.

           --- the next choose the configuration basis ---

          CISTEP = ALDET chooses the Ames Lab. determinant CI, and
                         requires $DET input. (default)
                 = GUGA  chooses the graphical unitary group CSFs,
                         and requires $DRT input.  This is the
                         only value usable with the QUAD converger.

           --- the next four choose the orbital optimizer ---

          FOCAS  = a flag to select a method with a first order
                   convergence rate.  (default=.FALSE.)

          SOSCF  = a flag selecting an approximately second order
                   convergence method.  (default=.TRUE.)

          FULLNR = a flag selecting a second order method, with an
                   exact orbital hessian.  (default=.FALSE.)

          QUAD   = a flag to pick a fully quadratic (orbital and
                   CI coefficient) optimization method, which is
                   applicable to FORS or non-FORS wavefunctions.
                   QUAD may not be used with state-averaging.
                   (default = .FALSE.)

          Note that FOCAS must be used only with FORS=.TRUE. in $DRT.
          The other convergers are usable for either FORS or non-FORS
          wavefunctions, although convergence is always harder in the
          latter case, when FORS below must be set .FALSE.

             --- the next apply to all convergence methods ---

          FORS   = a flag to specify that the MCSCF function is of
                   the Full Optimized Reaction Space type, which is
                   sometimes known as CAS-SCF.  .TRUE. means omit
                   act-act rotations from the optimization.  Since
                   convergence is usually better for FULLNR with
                   these rotations included, the default is sensible
                   for the case FORS=.TRUE. in $DRT.  (default is
                   .TRUE. for FOCAS/SOSCF, .FALSE. for FULLNR/QUAD)

          ACURCY = the major convergence criterion, the maximum
                   permissible asymmetry in the Lagrangian matrix.

          ENGTOL = a secondary convergence criterion, the run is
                   considered converged when the energy change is
                   smaller than this value. (default=1.0E-10)

          MAXIT  = Maximum number of iterations (default=100 for
                   FOCAS, 60 for SOSCF, 30 for FULLNR or QUAD)

          MICIT  = Maximum number of microiterations within a
                   single MCSCF iteration. (default=5 for FOCAS
                   or SOSCF, or 1 for FULLNR or QUAD)

          NWORD  = The maximum memory to be used, the default is
                   to use all available memory.  (default=0)

          CANONC = a flag to cause formation of the closed shell
                   Fock operator, and generation of canonical core
                   orbitals.  This will order the MCC core by their
                   orbital energies.  (default=.TRUE.)

          EKT    = a flag to cause generation of extended Koopmans'
                   theorem orbitals and energies.  (Default=.FALSE.)
              For this option, see R.C.Morrison and G.Liu,
              J.Comput.Chem., 13, 1004-1010 (1992).  Note that
              the process generates non-orthogonal orbitals, as
              well as physically unrealistic energies for the
              weakly occupied MCSCF orbitals.  The method is
              meant to produce a good value for the first I.P.

          NPUNCH = MCSCF punch option (analogous to $SCF NPUNCH)
                   0  do not punch out the final orbitals
                   1  punch out the occupied orbitals
                   2  punch out occupied and virtual orbitals
                       The default is NPUNCH = 2.

          NPFLG  = an array of debug print control.  This is
                   analagous to the same variable in $CIINP.
                   Elements 1,2,3,4,6,8 make sense, the latter
                   controls debugging the orbital optimization.

            --- the next refers to SOSCF optimizations ---

          NOFO   = set to 1 to skip use of FOCAS for one iteration
                   during SOSCF.  This is a testing parameter, at
                   present NOFO defaults to 0 to do one FOCAS iter.

             --- the next three refer to FOCAS optimizations ---

          CASDII = threshold to start DIIS (default=0.05)

          CASHFT = level shift value (default=1.0)

          NRMCAS = renormalization flag, 1 means do Fock matrix
                   renormalization, 0 skips (default=1)

             --- the next applies to the QUAD method ---
              (note that FULLNR input is also relevant)

          QUDTHR = threshold on the orbital rotation parameter,
                   SQCDF, to switch from the initial FULLNR
                   iterations to the fully quadratic method.
                   (default = 0.05)

             --- all remaining input applies only to FULLNR ---

          DAMP   = damping factor, this is adjusted by the program
                   as necessary.  (default=0.0)

          METHOD = DM2 selects a density driven construction of the
                   Newton-Raphson matrices.  (default).
                 = TEI selects 2e- integral driven NR construction.
                   See the 'further information' section for more
                   details concerning these methods.  TEI is slow!

          LINSER = a flag to activate a method similar to direct
                   minimization of SCF.  The method is used if
                   the energy rises between iterations.  It may in
                   some circumstances increase the chance of
                   converging excited states.  (default=.FALSE.)

          FCORE  = a flag to freeze optimization of the MCC core
                   orbitals, which is useful in preparation for
                   RUNTYP=TRANSITN jobs.  Setting this flag will
                   automatically force CANONC false.  This option
                   is incompatible with gradients, so can only be
                   used with RUNTYP=ENERGY.  It is a good idea to
                   decrease TOLZ and TOLE in $GUESS by two orders
                   of magnitude to ensure the core orbitals are
                   unchanged during input.  (default=.FALSE.)

                    --- the next four are seldom used ---

          DROPC  = a flag to include MCC core orbitals during the
                   CI computation.  The default is to drop them
                   during the CI, instead forming Fock operators
                   which are used to build the correct terms in
                   the orbital hessian. (default = .TRUE.)

          NORB   = the number of orbitals to be included in the
                   optimization, the default is to optimize with
                   respect to the entire basis.  This option is
                   incompatible with gradients, so can only be used
                   with RUNTYP=ENERGY.  (default=number of AOs
                   given in $DATA).

          MOFRZ  = an array of orbitals to be frozen out of the
                   orbital optimization step (default=none frozen).

          NOROT  = an array of up to 250 pairs of orbital rotations
                   to be omitted from the NR optimization process.
                   The program automatically deletes all core-core
                   rotations, all act-act rotations if FORS=.T.,
                   and all core-act and core-virt rotations if
                   FCORE=.T.  Additional rotations are input as
                   I1,J1,I2,J2... to exclude rotations between
                   orbital I running from 1 to NORB, and J running
                   up to the smaller of I or NVAL in $TRANS.



$MCQDPT group (relevant to SCFTYP=MCSCF if MPLEVL=2)

               Controls 2nd order MCQDPT (multiconfiguration quasi-
          degenerate perturbation theory) runs, if requested by
          MPLEVL=2 in $CONTRL.  MCQDPT2 is implemented only for
          FORS (aka CASSCF) wavefunctions.   The MCQDPT method is a
          multistate, as well as multireference perturbation theory.
          The implementation is a separate program, interfaced to
          GAMESS, with its own procedures for determination of the
          canonical MOs, CSF generation, integral transformation,
          CI in the reference CAS, etc.  Therefore some of the input
          in this group repeats data given elsewhere, particularly
          the $DET/$DRT.  A more complete discussion may be found in
          the 'Further Information' chapter.  Analytic gradients are
          not available.

                 *** MCSCF reference wavefunction ***

          NEL    =   total number of electrons, including core.
                     (default from $DATA and ICHARG in $CONTRL)

          MULT   =   spin multiplicity (default from $CONTRL)

          NMOACT =   Number of orbitals in FORS active space
                     (default is the active space in $DET or $DRT)

          NMOFZC =   number of frozen core orbitals, NOT correlated
                     in the perturbation calculation.  (default is
                     number of chemical cores)

          NMODOC =   number of orbitals which are doubly occupied in
                     every MCSCF configuration, that is, not active
                     orbitals, which are to be included in the
                     perturbation calculation.  (The default is all
                     valence orbitals between the chemical core and
                     the active space)

          NMOFZV =   number of frozen virtuals, NOT occupied during
                     the perturbation calculation.  The default is
                     to use all virtuals in the MP2.  (default=0)

          NOSYM  = 0 use CSF symmetry (see the ISTSYM keyword).
                     off diagonal perturbations vanish if states are
                     of different symmetry, so the most efficient
                     computation is a separate run for every space
                     symmetry. (default)
                   1 turn off CSF state symmetry so that all states
                     are treated at once.  ISTSYM is ignored.

          ISTSYM =   the state symmetry of the target state(s).
                     This is given as an integer, note that only
                     Abelian groups in $DATA are supported:
                       ISTSYM= 1   2   3   4   5   6   7   8
                          C1   A
                          Ci   Ag  Au
                          Cs   A'  A''
                          C2   A   B
                          C2v  A1  A2  B1  B2
                          C2h  Ag  Bg  Au  Bu  <- differs $DRT/$DET!
                          D2   A   B1  B2  B3
                          D2h  Ag  B1g B2g B3g Au  B1u B2u B3u
                     (The default is 1, the totally symmetric state)

                 *** perturbation specification ***

          KSTATE=    state is used (1) or not (0) in the MCQDPT2.
                     Maximum of 20 elements, including zeros.
                     For example, if you want the perturbation
                     correction to the second and the fourth roots,

                 *** MO input and flow control ***

          INORB  = 0 optimize the MCSCF wavefunction in this run.
                 = 1 read the converged orbitals from a $VEC group,
                     and skip immediately to the MCQDPT computation.
                     A complete $VEC including virtuals must be given.

                 *** Canonical Fock orbitals ***

          IFORB  = 1 determine the canonical Fock orbitals
                 = 0 omit this step.  (default=1)

          WSTATE =   weight of each CAS-CI state in computing the
                     closed shell Fock matrix.  You must enter 0.0
                     whenever the same element in KSTATE is 0.
                     (default is WSTATE(1)=1.0,0.0,0.0,...)

                 *** Miscellaneous options ***

          REFWGT =   a flag to request decomposition of the second
                     order energy into internal, semi-internal, and
                     external contributions, and to obtain the weight
                     of the MCSCF reference in the 1st order wave
                     function.  This option significantly increases
                     the run time!  (default = .FALSE.)

          THRGEN =   threshold for one-, two-, and three-body
                     density matrix elements in the perturbation
                     calculation.  If you want to obtain energies,
                     for instance, to 6 figures after point, choose
                     THRGEN=1.0D-08 or 1.0D-09.  (default=1.0D-08)

          THRENE =   threshold for the energy convergence in the
                     Davidson's method CAS-CI.  (default=-1.0D+00)

          THRCON =   threshold for the vector convergence in the
                     Davidson's method CAS-CI.  (default=1.0D-06)

          LPOUT =    print option, 0 gives normal printout, while
                     <0 gives debug print (e.g. -1, -5, -10, -100)

          Finally, there are additional very specialized input
          options, described in the source code routine MQREAD:




$CISORT group (optional, relevant for -CI- and -MCSCF-)

               This group provides further control over the sorting
          of the transformed integrals.

          NDAR   = Number of direct access records.
                   (default = 2000)

          LDAR   = Length of direct access record (site dependent)

          NBOXMX = Maximum number of boxes in the sort.
                   (default = 200)

          NWORD  = Number of words of fast memory to use in this
                   step.  A value of 0 results in automatic use of
                   all available memory.  (default = 0)

          NOMEM  = 0 (set to one to force out of memory algorithm)



$GUGEM group (optional, relevant for -CI- or -MCSCF-)

              This group provides further control over the
          calculation of the energy (Hamiltonian) matrix.

          CUTOFF = Cutoff criterion for the energy matrix.

          NWORD  = not used.



$GUGDIA group (optional, relevant for -CI- or -MCSCF-)

               This group provides control over the Davidson method
          diagonalization step.

          NSTATE = Number of CI states to be found. (default=1)
                   You can solve for any number of states, but only
                   100 can be saved for subsequent sections, such
                   as state averaging.

          PRTTOL = Printout tolerance for CI coefficients
                   (default = 0.05)

          MXXPAN = Maximum no. of expansion basis vectors used
                   before the expansion basis is truncated.

          ITERMX = Maximum number of iterations (default=50)

          CVGTOL = Convergence criterion for Davidson eigenvector
                   routine.  This value is proportional to the
                   accuracy of the coeficients of the eigenvector(s)
                   found.  The energy accuracy is proportional to
                   its square.  (default = 1.0E-5)

          NWORD  = Number of words of fast memory to use in this
                   step.  A value of zero results in the use of all
                   available memory.  (default = 0)

          MAXHAM = specifies dimension of Hamiltonian to try to
                   store in memory.  The default is to use all
                   remaining memory to store this matrix in memory,
                   if it fits, to reduce disk I/O to a minimum.

          MAXDIA = maximum dimension of Hamiltonian to send to an
                   incore diagonalization.  If the number of CSFs
                   is bigger than MAXDIA, an iterative Davidson
                   procedure is invoked.  Default=100

          NIMPRV = Maximum no. of eigenvectors to be improved every
                   iteration. (default = nstate)

          NSELCT = Determines initial guess to eigenvectors.
                   = 0 ->  Unit vectors corresponding to the NSTATE
                           lowest diagonal elements and any diagonal
                           elements within SELTHR of them. (default)
                   < 0 ->  First abs(NSELCT) unit vectors.
                   > 0 ->  use NSELCT unit vectors corresponding to
                            the NSELCT lowest diagonal elements.

          SELTHR = Guess selection threshold when NSELCT=0.

          NEXTRA = Number of extra expansion basis vectors to be
                   included on the first iteration.  NEXTRA is
                   decremented by one each iteration.  This may be
                   useful in "capturing" vectors for higher states.

          KPRINT = Print flag bit vector used when
                   NPFLG(4)=1 in the $CIINP group       (default=8)
                   value  1 bit 0 print final eigenvalues
                   value  2 bit 1 print final tolerances
                   value  4 bit 2 print eigenvalues and tolerances
                                  at each truncation
                   value  8 bit 3 print eigenvalues every iteration
                   value 16 bit 4 print tolerances every iteration



$GUGDM group (optional, relevant for -CI-)

               This group provides further control over formation of
          the one electron density matrix.  See NSTATE in $GUGDIA.

          NFLGDM = Controls each state's density formation.
                   0 -> do not form density for this state.
                   1 -> form density and natural orbitals for this
                        state, print and punch occ.nums. and NOs.
                   2 -> same as 1, plus print density over MOs.
                   (default=1,99*0, means ground state NOs only)
                   Note that forming the 1-particle density for a
                   state is negligible against the diagonalization
                   time for that state.

          IROOT  = The -CI- root whose density matrix is saved on
                   the direct access dictionary file for later
                   computation of properties.  You may save only
                   one state's density for property evaluation.

          WSTATE = An array of up to 100 weights to be given to the
                   1 body density of each state in forming the DM1.
                   It is not physically reasonable to average over
                   any CI states that are not degenerate, but it
                   may be useful to use WSTATE to produce a totally
                   symmetric density when the states are degenerate.
                   The averaged density will be used for property
                   computations, as well as to generate natural
                   orbitals.  The default is to use NFLGDM/IROOT,
                   unless WSTATE information is given, in which case
                   NFLGDM/IROOT are ignored.

          IBLOCK = Density blocking switch. If nonzero, the off
                   diagonal block of the density below row IBLOCK
                   will be set to zero before the (approximate)
                   natural orbitals are found.  One use for this is
                   to keep the internal and external orbitals in a
                   FOCI or SOCI calculation from mixing, in which
                   case IBLOCK is the highest occupied internal
                   orbital.  (default=0)

          NWORD  = Number of words of memory to use.  Zero means use
                   all available memory (default=0).



$GUGDM2 group (optional, relevant for -CI- or -MCSCF-)

               This group provides control over formation of the
          2-particle density matrix.

          WSTATE = An array of up to 100 weights to be given to the
                   2 body density of each state in forming the DM2.
                   The default is to optimize a pure ground state.
                   A small amount of the ground state can help the
                   convergence of excited states greatly.
                   Gradient runs are possible only with pure states.

                   Be sure to set NSTATE in $GUGDIA appropriately!

          CUTOFF = Cutoff criterion for the 2nd-order density.
                   (default = 1.0E-9)

          NWORD  = Number of words of fast memory to use in sorting
                   the DM2.  The default uses all available memory.

          NOMEM  = 0 uses in memory sort, if possible.
                 = 1 forces out of memory sort.

          NDAR   = Number of direct access records. (default=4000)

          LDAR   = Length of direct access record (site dependent)

          NBOXMX = Maximum no. of boxes in the sort. (default=200)



$LAGRAN group (optional, relevant for -CI- gradient)

               This group provides further control over formation of
          the CI Lagrangian, a quantity which is necessary for the
          computation of CI gradients.

           NOMEM =   0 form in core, if possible
                 =   1 forces out of core formation
           NWORD =   0 (0=use all available memory)
           NDAR  = 4000
           LDAR  = Length of each direct access record
                   (default is NINTMX from $INTGRL)



$TRFDM2 group (optional, relevant for -CI- gradient)

               This group provides further control over the back
          transformation of the 2 body density to the AO basis.

           NOMEM =   0 transform and sort in core, if possible
                 =   1 transform in core, sort out of core, if poss.
                 =   2 transform out of core, sort out of core
           NWORD =   0 (0=use all available memory)
           CUTOFF= 1.0D-9, threshold for saving DM2 values
           NDAR  = 2000
           LDAR  = Length of each direct access record
                   (default is system dependent)
           NBOXMX= 200


          Usually neither of these two groups is given.  Since these
          groups are normally used only for CI gradient runs, we
          list here some of the restrictions on the CI gradients:

            a) SCFTYP=RHF, only
            b) no FZV orbitals in $CIDRT, all MOs must be used.
            c) the derivative integrals are computed in the 2nd
               derivative code, which is limited to spd basis sets.
            d) the code does not run in parallel.
            e) Use WSTATE in $GUGDM2 to specify the state whose
               gradient is to be found.  Use IROOT in $GUGDM to
               specify the state whose other properties will be
               found.  These must be the same state!
            f) excited states often have different symmetry than the
               ground state, so think about GROUP in $CIDRT.
            g) the gradient can probably be found for any CI for
               which you have sufficient disk to do the CI itself.
               Time is probably about 2/3 additional.



$TRANST group (relevant for RUNTYP=TRANSITN & only for CITYP=GUGA)

              This group controls the evaluation of the radiative
          transition moment, or spin orbit coupling (SOC).  Defaults
          assume that there is one common set of orbitals, all of
          which are occupied.  The program can use two separately
          optimized MO sets, provided certain conditions are met.
          If relativistic corrections were included in the underlying
          spin-free wavefunctions, it is possible to either include
          or neglect the corrections to SOC, see NESOC in $RELWFN.

          OPERAT selects the type of transition being computed.
                 = DM      calculates radiative transition moment
                           between states of same spin, using
                           the dipole moment operator. (default)
                 = HSO1    one-electron Spin-Orbit Coupling (SOC)
                 = HSO2P   partial two electron and full 1e- SOC,
                           namely core-active 2e- contributions are
                           computed, but active-active 2e- terms
                           are ignored.  This generally captures
                           >90% of the full HSO2 computation, but
                           with spin-orbit matrix element time
                           similar to the HSO1 calculation.
                 = HSO2    one and two-electron SOC, this is the
                           full Pauli-Breit operator.
                 = HSO2FF  one and two-electron SOC, the form factor
                           method gives the same result as HSO2, but
                           is more efficient in the case of small
                           active spaces, small numbers of CI states,
                           and large atomic basis sets.

              * * * next define orbitals and wavefunctions * * *

          NUMVEC = the number of different MO sets. This can be
                   either 1 or 2, but 2 can be chosen only for
                   FORS/CASSCF or FCI wavefunctions. (default=1)
                   Note that $GUESS is not read by this RUNTYP!
                   If you set NUMVEC=2 and you use symmetry in any
                   of the $DRTx groups, you may have to use ISTSYM
                   in the $DRT groups since the order of orbitals
                   from the corresponding orbital transformation
                   is unpredictable.

          NFZC   = When NUMVEC=2, this is the number of identical
                   core orbitals in the two vector sets, and must
                   be equal to NFZC in all $DRTx groups.  When
                   NUMVEC=1, it is the value of NFZC in the $DRTx.
                   The default is the number of AOs given in $DATA,
                   again this is not very reasonable.

          NOCC   = the number of occupied orbitals.  The default is
                   the number of AOs given in $DATA, which is not
                   usually correct.

          NUMCI  = the number of different CI calculations to do.
                   You may wish to define one $DRTx group for each
                   irreducible representation within each spin
                   multiplicity.  NUMCI may not exceed 64.  IVEX,
                   IROOTS, NSTATE, ENGYST below will all have NUMCI
                   values.  (default=1)

          IVEX   = array of indices of $VECx groups to be used for
                   each CI calculation.  The default for NUMVEC=2 is
                   IVEX(1)=1,2,1,1,1,1,1..., and of course for
                   NUMVEC=1, it is IVEX(1)=1,1,1,1,1...

          IROOTS = array containing the number of CI states for which
                   the transition moments are to be found.  The
                   default is 1 for each CI, which is probably not a
                   correct choice for OPERAT=DM runs, but is quite
                   reasonable for the HSO operators.  See also ETOL.

          ETOL   = energy tolerance for CI state elimination.  After
                   a CI computation finding up to NSTATE(i) CI roots
                   within each $DRTx, the number of states kept in
                   the computation is limited to the number given
                   by IROOTS, with a further limitation that the
                   state be within ETOL of the lowest energy state
                   found for any of the $DRTx.  The default is 100.0
                   Hartree, so that IROOTS is the only limitation.
                   This applies only to OPERAT=HSO1,2,2P.

          NSTATE = array of CI states to be found when diagonalising
                   the CI Hamiltonians.  Of these, the first IROOTS(i)
                   states will be used to find transition moments.
                   The default for NSTATE(i) is IROOTS(i).

          ENGYST = energy values to replace the CI spin-free energies.
                   A possible use for this is to use SOCI energies
                   on the diagonal of the Hamiltonian (obtained in
                   previous runs) but to use only FORS wavefunctions
                   to evaluate off diagonal HSO matrix elements.  The
                   CI runs are still conducted to obtain CI coefs,
                   needed to evaluate the off diagonal elements.
                   Enter MXRT*NUMCI values as a square array, by the
                   usual FORTRAN convention (that is, MXRT roots of
                   $DRT1, MXRT roots of $DRT2 etc), in hartrees, with
                   zeros added to fill each column to MXRT values.
                   MXRT is the maximum value in the IROOTS array.
                   (the default is the computed CI energies)

                  * * * some control tolerances * * *

          NOSYM= -1 forces use of symmetry-contaminated orbitals
                    symmetry analysis, otherwise the same as NOSYM=0
               =  0 fully use symmetry
               =  1 do not use point group symmetry, but still use
                    other symmetries (Hermiticity, spin).
               =  2 use no symmetry.   Also, include all CSFs for
                    HSO1, 2, 2P.
               =  3 force the code to assume the symmetry specified
                    in $DATA is the same as in all $DRT groups, but
                    is otherwise identical to NOSYM=-1.  This option
                    saves CPU time and money(memory).  Since the $DRT
                    works by mapping non-Abelian groups into their
                    highest Abelian subgroup, do not use NOSYM=3 for
                    non-Abelian groups.

          SYMTOL = relative error for the matrix elements.  This
                   parameter has a great impact upon CPU time, and
                   the default has been chosen to obtain nearly
                   full accuracy while still getting good speedups.

             * * * the next pertain only to spin-orbit runs * * *

          ZEFTYP          specifies effective nuclear charges to use
                 = TRUE   uses true nuclear charge of each atom,
                          except protons are removed if an ECP basis
                          is being used (default).
                 = 3-21G  selects values optimized for the 3-21G
                          basis, but these are probably appropriate
                          for any all electron basis set.  Rare gases,
                          transition metals, and Z>54 will use the
                          true nuclear charges.
                 = SBKJC  selects a set obtained for the SBKJC ECP
                          basis set, specifically.  It may not be
                          sensible to use this for other ECP sets.
                          Rare gases, lanthanides, and Z>86 will use
                          the true nuclear charges.

          ZEFF   = an array of effective nuclear charges, overriding
                   the charges chosen in ZEFTYP.

              Note that effective nuclear charges can be used for
              any HSO type OPERAT, but traditionally these are used
              mainly for HSO1 as an empirical correction to the
              omission of the 2e- term, or to compensate for missing
              core orbitals in ECP runs.

          JZ       controls the calculation of Jz eigenvalues
                 = 0 do not perform the calculation
                 = 1 do the calculation
                   By default, jz is set to 1 for molecules that are
                   recognised as linear (this includes atoms!).
                   The matrix of Jz=Lz+Sz operator is constructed
                   between spin-mixed states (eigenvalues of Hso).
                   Setting Jz to 1 can enforce otherwise avoided (by
                   symmetry) calculations of SOC matrix elements.
                   Interpretation of Jz eigenvalues for systems other
                   than linear molecules aligned along z-axis is left
                   at your discretion.  JZ applies only to HSO1,2,2P.

          TMOMNT = flag to control computation of the transition
                   dipole moment between spin-mixed wavefunctions
                   (that is, betweeen eigenvectors of the Pauli-Breit
                   Hamiltonian).  Applies only to HSO1,2,2P.
                   (default is .FALSE.)

          SKIPDM = flag to omit(.TRUE.) or include(.FALSE.) dipole
                   moment matrix elements during spin-orbit coupling.
                   Usually it takes almost no addition effort to
                   calculate  excluding some cases when the
                   calculation of forbidden by symmetry spin-orbit
                   coupling matrix elements  may have to be
                   performed since  and  are computed
                   simultaneously.  Applies only to HSO1,2,2P
                   (default is .TRUE.)

          IPRHSO = controls output style for matrix elements (HSO*)
                 =-1 do not output individual matrix elements
                 otherwise these are accumulative:
                 = 0 term-symbol like kind of labelling:
                     labels contain full symmetry info (default)
                 = 1 all states are numbered consequently within each
                     spin multiplicity (ye olde style)
                 = 2 output only nonzero (>=1e-4) matrix elements

          PRTPRM = flag to provide detailed information about the
                   composition of the spin-mixed states in terms of
                   adiabatic states. This flag also provides similar
                   information about Jz (if JZ set).
                   (default is .FALSE.)

               * * * expert mode HSO control options * * *

          MODPAR =    parallel options, which are independent bit
                      options, 0=off, 1=on.  Bit 1 refers only to
                      HSO2FF, bit 2 to HSO1,2,2P.  Enter a decimal
                      value 0, 1, 2, 3 meaning binary 00, 01, 10, 11.
           bit 1 = 0/1 (HSO2FF) uses static/dynamic load balancing in
                      parallel if available, otherwise use static
                      load balancing.  Dynamic algorithm is usually
                      faster but may utilize memory less efficiently,
                      and I/O can slow it down.  Also, dynamical
                      algorithm forces SAVDSK=.F. since its
                      unique distribution of FFs among nodes implies
                      no savings from precalculating form factors.
           bit 2 = 0/1 (HSO1,2,2P) duplicate/distribute SOC integrals
                      in parallel.  If set, 2e AO integrals and the
                      four-index transformation are divided over
                      nodes (distributed), and SOC MO integrals are
                      then summed over nodes.
           The default is 3, meaning both bits are set on (11)

          PHYSRC = flag to force the size of the physical record to
                   be equal to the size of the sorting buffers.
                   This option can have a dramatic effect on the
                   efficency.  Usually, setting PHYSRC=.t. is helpful
                   if the code complains that low memory enforces
                   SLOWFF=.TRUE., or you set it yourself. For large
                   active spaces and large memory (more precisely, if
                   reclen is larger than the physical record size)
                   PHYSRC=.TRUE. can slow the code down.  Setting
                   PHYSRC to .true. forces SLOWFF to be .false.
                   See MODPAR. (default .FALSE.) (only with HSO2FF)

          RECLEN = specifies the size of the record on file 40,
                   where form factors are stored. This parameter
                   significantly affects performance.
                   If not specified, RECLEN have to be guessed,
                   and the guess will usually be either an
                   overestimate or underestimate. If the former
                   you waste disk space, if the latter the program
                   aborts. Note that RECLEN will be different for
                   each pair of multiplicities and you must specify
                   the maximum for all pairs.  The meaning of this
                   number is how many non-zero form factors are
                   present given four MO indices.  You can decrease
                   RECLEN if you are getting a message "predicted
                   sorting buffer length is greater than needed..."
                   Default depends on active space. (only HSO2FF)

          SAVDSK = flag to repeat the form factor calculation twice.
                   This avoids wasting disk space as the actually
                   required record size is found during the 1st run.
                   (default=.FALSE.) (only with HSO2FF)

          SLOWFF = flag to choose a slower FF write-out method.
                   By default .FALSE., but this is turned on if:
                   1) not enough memory for the fast way is available
                   2) the maximum usable memory is available, as when
                      the buffer is as large as the maximum needed,
                      then the "slow FF" algorythm is faster.
                   Generally SLOWFF=.true. saves up to 50% or so of
                   disk space.  See PHYSRC.  (only with HSO2FF)

          ACTION          controls disk file DAFL30 reuse.
                 = NORMAL calculate the form factors in this run.
                 = SAVE   calculate, and store the form factors on
                          disk for future runs with the same active
                          space characteristics.
                 = READ   read the form factors from disk from an
                          earlier run which used SAVE.
                   (default=NORMAL) (only with HSO2FF)
                   Note that currently in order to use ACTION =
                   SAVE or READ you should specify MS= -1, 0, or 1

          MS     = refers to the difference between ket and bra's Ms
                   It pertains only to HSO2FF.
                  -1 do matrix elements for ms=-1 only
                   0 do matrix elements for ms=0 only
                   1 do matrix elements for ms=1 only
                  -2 do matrix elements for all ms (0, 1, and -1),
                     which is the default.
                  -3 calculates form factors so they can be saved.

          * * * the remaining parameters are not so important * * *

          PRTCMO = flag to control printout of the corresponding
                   orbitals.  (default is .FALSE.)

          HSOTOL = HSO matrix elements are considered zero if they
                   are smaller than HSOTOL.  This parameter is used
                   only for print-out and statistics.
                   (default=1.0E-1 cm-1)

          TOLZ   = MO coefficient zero tolerance (as for $GUESS).

          TOLE   = MO coefficient equating tolerance (as for
                   $GUESS).  (default=1.0E-5)


Informations required by Austrian law (Offenlegung gem. §25 MedienG): Dr. Michael Ramek, Graz