Документ взят из кэша поисковой машины. Адрес оригинального документа : http://erg.biophys.msu.ru/~comcon1/qmdata/gamess-us-input-2014dec.pdf
Дата изменения: Sun Apr 5 09:23:30 2015
Дата индексирования: Sat Apr 9 22:28:09 2016
Кодировка:
Input Description

2-1 (5 December 2014) ********************************* * * * Section 2 - Input Description * * * *********************************

This section of the manual describes the input to GAMESS. The section is written in a reference, rather than tutorial fashion. However, there are frequent reminders that more information can be found on a particular input group, or type of calculation, in the 'Further Information' section of this manual. Numerous complete input files are shown in the 'Input Examples' section. Note that this chapter of the manual can be searched online by means of the "gmshelp" command, if your computer runs Unix. A command such as gmshelp scf will display the $SCF input group. With no arguments, the gmshelp command will show you all of the input group names. Type "" to see the next screen, "b" to back up to the previous screen, and "q" to exit the pager. If gmshelp does not work, ask the person who installed GAMESS to fix the 'gmshelp' script, as it is extremely useful. The order of this section is chosen to approximate the order in which most people prepare their input ($CONTRL, $BASIS/$DATA, $GUESS, and so on). The next few pages contain a list of all possible input groups, grouped in this way. The PDF version of this file contains an index of all group names in alphabetical order.


Input Description name function module:routine -----------------------Molecule, basis set, wavefunction specification: $CONTRL chemical control data INPUTA:START $SYSTEM computer related options INPUTA:START $BASIS basis set INPUTB:BASISS $DATA molecule, geometry, basis set INPUTB:MOLE $ZMAT internal coordinates ZMATRX:ZMATIN $LIBE linear bend coordinates ZMATRX:LIBE $SCF HF-SCF wavefunction control SCFLIB:SCFIN $SCFMI SCF-MI input control data SCFMI :MIINP $DFT density functional theory DFT :DFTINP $TDDFT time-dependent DFT TDDFT :TDDINP $CIS singly excited CI CISGRD:CISINP $CISVEC vectors for CIS CISGRD:CISVRD $MP2 2nd order Moller-Plesset MP2 :MP2INP $RIMP2 resolution of the identity MP2 RIMP2 :RIDRVR $AUXBAS RI-MP2's basis set specification RIMP2 :RIDRVR $CCINP coupled cluster input CCSDT :CCINP $EOMINP equation of motion CC EOMCC :EOMINP $MOPAC semi-empirical specification MPCMOL:MOLDAT $GUESS initial orbital selection GUESS :GUESMO $VEC orbitals (formatted) GUESS :READMO $MOFRZ freezes MOs during SCF runs EFPCOV:MFRZIN $DFTB DFTB input DFTBLB:INPUT $DFTBSK Slater-Koster table input DFTBSK:SKTAB Note that MCSCF and CI input is listed below. Potential energy surface options: $STATPT $TRUDGE $TRURST $FORCE $CPHF $CPMCHF $MASS $HESS $GRAD $DIPDR $ALPDR $VIB $VIB2 $VSCF $VIBSCF $GAMMA $EQGEOM geometry search control nongradient optimization restart data for TRUDGE hessian, normal coordinates coupled-Hartree-Fock options coupled-MR-Hartree-Fock options isotope selection force constant matrix (formatted) gradient vector (formatted) dipole deriv. matrix (formatted) alpha polar. der. (formatted) HESSIAN restart data (formatted) num GRAD/HESS restart (formatted) vibrational anharmonicity VSCF restart data (formatted) 3rd nuclear derivatives equilibrium geometry data STATPT:SETSIG TRUDGE:TRUINP TRUDGE:TRUDGX HESS :HESSX CPHF :CPINP MCPCGX:MCPCGX VIBANL:RAMS HESS :FCMIN HESS :EGIN HESS :DDMIN RAMAN :ADMIN HESS :HSSNUM HESS :HSSFUL VSCF :VSCFIN VSCF :VGRID HESS :GAMMXX HESS :FFCARX

2-2 *


Input Description $HLOWT $GLOWT $IRC $DRC $MEX $CONICL $MD $RDF $GLOBOP $GLBFRG $GRADEX $SURF hessian data from equilibrium 3rd derivatives at equilibrium intrinsic reaction coordinate dynamic reaction path minimum energy crossing point conical intersection search molecular dynamics trajectory radial dist. functions for MD Monte Carlo global optimization Monte Carlo atom groups gradient extremal path potential surface scan HESS :FFCARX HESS :FFCARX RXNCRD:IRCX DRC :DRCDRV MEXING:MEXINP MDEFP :MDX MDEFP :RDFX GLOBOP:GLOPDR GLOBOP:GLOPDR GRADEX:GRXSET SURF :SRFINP

2-3

Interpretation, properties: $LOCAL $TRUNCN $ELMOM $ELPOT $ELDENS $ELFLDG $POINTS $GRID $PDC $MGC $RADIAL $MOLGRF $STONE $COMP $RAMAN $NMR $MOROKM $LMOEDA $QMEFP $FFCALC $TDHF $TDHFX localized molecular orbitals localized orbital truncations electrostatic moments electrostatic potential electron density electric field/gradient property calculation points property calculation mesh MEP fitting mesh mean gradient charges atomic orbital radial data orbital plots distributed multipole analysis thermochemical calculation Raman intensity NMR shielding tensors Morokuma energy decomposition LMO-based energy decomposition QM/EFP energy decomposition finite field polarizabilities time dependent HF of NLO props TDHF for NLO, Raman, hyperRaman LOCAL :LMOINP EFPCOV:TRNCIN PRPLIB:INPELM PRPLIB:INPELP PRPLIB:INPELD PRPLIB:INPELF PRPLIB:INPPGS PRPLIB:INPPGS PRPLIB:INPPDC PRPPOP:RADWFN PARLEY:PLTMEM PRPPOP:STNRD COMP :COMPX RAMAN :RAMANX NMR :NMRX MOROKM:MOROIN MOROKM:MMOEDIN EFINP :QMEFPAX FFIELD:FFLDX TDHF :TDHFX TDX:FINDTDHFX

Solvation models: $EFRAG use effective fragment potential $FRAGNAME specifically named fragment pot. $FRGRPL inter-fragment repulsion $EWALD Ewald sums for EFP electrostatics $MAKEFP generate effective fragment pot. $PRTEFP simplified EFP generation $DAMP EFP multipole screening fit $DAMPGS initial guess screening params $PCM polarizable continuum model EFINP :EFINP EFINP :RDSTFR EFINP :RDDFRL EWALD :EWALDX EFINP :EFPX EFINP :PREFIN CHGPEN:CGPINP CHGPEN:CGPINP PCM :PCMINP


Input Description $PCMGRD $PCMCAV $TESCAV $NEWCAV $IEFPCM $PCMITR $DISBS $DISREP $SVP $SVPIRF $COSGMS $SCRF PCM gradient control PCM cavity generation PCM cavity tesselation PCM escaped charge cavity PCM integral equation form. data PCM iterative IEF input PCM dispersion basis set PCM dispersion/repulsion Surface Volume Polarization model reaction field points (formatted) conductor-like screening model self consistent reaction field PCMCV2:PCMGIN PCM :MAKCAV PCMCV2:TESIN PCM :DISREP PCM :IEFDAT PCMIEF:ITIEFIN PCMDIS:ENLBS PCMVCH:MORETS SVPINP:SVPINP SVPINP:SVPIRF COSMO :COSMIN SCRF :ZRFINP

2-4

Integral, and integral modification options: $ECP $MCP $RELWFN $EFIELD $INTGRL $FMM $TRANS effective core potentials model core potentials scalar relativistic integrals external electric field 2e- integrals fast multipole method integral transformation ECPLIB:ECPPAR MCPINP:MMPRED INPUTB:RWFINP PRPLIB:INPEF INT2A :INTIN QMFM :QFMMIN TRANS :TRFIN

Fragment Molecular Orbital method: $FMO $FMOPRP $FMOXYZ $OPTFMO $FMOHYB $FMOBND $FMOENM $FMOEND $OPTRST $GDDI define FMO fragments FMO properties and convergers atomic coordinates for FMO input for special FMO optimizer localized MO for FMO boundaries FMO bond cleavage definition monomer energies for FMO restart dimer energies for FMO restart OPTFMO restart data group DDI definition FMOIO :FMOMIN FMOIO :FMOPIN FMOIO :FMOXYZ FMOGRD:OPTFMO FMOIO :FMOLMO FMOIO :FMOBON FMOIO :EMINOU FMOIO :EDIN FMOGRD:RSTOPT INPUTA:GDDINP

Polymer model: $ELG polymer elongation method ELGLIB:ELGINP

Divide and conquer model: $DANDC $DCCORR $SUBSCF $SUBCOR $MP2RES $CCRES DC SCF input DC correlation method input subsystem definition for SCF subsystem definition for MP2/CC restart data for DC-MP2 restart data for DC-CC DCLIB DCLIB DCLIB DCLIB DCMP2 DCCC :DCINP :DCCRIN :DFLCST :DFLCST :RDMPDC :RDCCDC


Input Description clusters in molecules $CIMINP $CIMATM $CIMFRG controls clusters in molecules fine tune calculation level fine tune atomic fragmentation CIMINF:CIMINP CIMINF:CIMINP CIMINF:CIMPRT

2-5

quantum mechanics/molecular mechanics model: $FFDATA $FFPDB QuanPol calculation for molecules QUANPO:QUANPOL QuanPol calculation for proteins QUANPO:QUANPOL

MCSCF and CI wavefunctions, and their properties: $CIINP $DET $CIDET $GEN $CIGEN $ORMAS $CEEIS $CEDATA $GCILST $GMCPT $PDET $ADDDET $REMDET $SODET $DRT $CIDRT $MCSCF $MRMP $DETPT $MCQDPT $EXCORR $CASCI $IVOORB $CISORT $GUGEM $GUGDIA $GUGDM $GUGDM2 $LAGRAN $TRFDM2 $DIABAT $TRANST control over CI calculation determinant full CI for MCSCF determinant full CI determinant general CI for MCSCF determinant general CI determinant multiple active space CI energy extrapolation restart data for CEEIS general MCSCF/CI determinant list general MCSCF/CI determinant list parent determinant list add determinants to reference remove determinants from ref. determinant second order CI GUGA distinct row table for MCSCF GUGA CI (CSF) distinct row table control over MCSCF calculation MRPT selection det. multireference pert. theory CSF multireference pert. theory interface to MPQC's R12 program IVO-CASCI input fine tuning of IVO-CASCI GUGA CI integral sorting GUGA CI Hamiltonian matrix GUGA CI diagonalization GUGA CI 1e- density matrix GUGA CI 2e- density matrix GUGA CI Lagrangian GUGA CI 2e- density backtransform diabatic states transition moments, spin-orbit GAMESS:WFNCI ALDECI:DETINP ALDECI:DETINP ALGNCI:GCIINP ALGNCI:GCIINP ORMAS :FCINPT CEEIS :CEEISIN CEEIS :RDCEEIS ALGNCI:GCIGEN GMCPT :OSRDDAT GMCPT :OSMKREF GMCPT :OSMKREF GMCPT :OSMKREF FSODCI:SOCINP GUGDRT:ORDORB GUGDRT:ORDORB MCSCF :MCSCF MP2 :MRMPIN DEMRPT:DMRINP MCQDPT:MQREAD EXCORR:GETEXC IVOCAS:IVODRV IVOCAS:ORBREAD GUGSRT:GUGSRT GUGEM :GUGAEM GUGDGA:GUGADG GUGDM :GUGADM GUGDM2:GUG2DM LAGRAN:CILGRN TRFDM2:TRF2DM DIAB:DIABINP TRNSTN:TRNSTX

* this column is more useful to programmers than to users.


Input Description

$CONTRL

2-6

==========================================================

$CONTRL

group

(note:

only one "oh"!)

This group specifies the type of wavefunction, the type of calculation, use of core potentials, spherical harmonics, coordinate choices, and similar fundamental job options. Because this is a very long input group, here is a short list of its most important keywords: SCFTYP, MPLEVL, CITYP, CCTYP, DFTTYP, TDDFT RUNTYP, ICHARG, MULT, RELWFN/PP, NZVAR, ISPHER SCFTYP = RHF = UHF = ROHF = GVB = MCSCF = NONE specifies the self-consistent field wavefunction. You may choose from Restricted Hartree Fock calculation (default) Unrestricted Hartree Fock calculation Restricted open shell Hartree-Fock. (high spin, see GVB for low spin) Generalized valence bond wavefunction, or low spin ROHF. (needs $SCF input) Multiconfigurational SCF wavefunction (this requires $DET or $DRT input) indicates a single point computation, rereading a converged SCF function. This option requires that you select CITYP=ALDET, ORMAS, FSOCI, GENCI, or GUGA, requesting only RUNTYP=ENERGY or TRANSITN, and using GUESS=MOREAD.

The treatment of electron correlation for the above SCF wavefunctions is controlled by the keywords DFTTYP, VBTYP, MPLEVL, CITYP, and CCTYP contained in this group. No more than one of these may be chosen in a single run (except as part of RUNTYP=SURFACE). Scalar relativistic effects may be incorporated using RELWFN for any of these wavefunction choices, correlated or not.


Input Description

$CONTRL

2-7

DFTTYP = NONE = XXXXXX

ab initio computation (default) perform density functional theory run, using the functional specified. Many choices for XXXXXX are listed in the $DFT and $TDDFT input groups. no excited states (default) generate time-dependent DFT excitation energies, using the DFTTYP= functional, for RHF or UHF references. Analytic nuclear gradients are available for RHF. See $TDDFT. spin-flip TD-DFT, for either UHF or ROHF references. Nuclear gradients and solvent effects are coded. See $TDDFT. (hyper)polarizability calculation, for RHF only. See $TDDFT. *****

TDDFT

= NONE = EXCITE

= SPNFLP = POL

VBTYP

= NONE = VB2000

no valence bond calculation (default) use the VB2000 program to generate VB wavefunctions, for SCFTYP=RHF or ROHF. Analytic nuclear gradients are not available. A $VB2000 input group is required. See ~/gamess/vb2000/DOC/readme.GAMESS for info about $VB2000, and see also http://www.scinetec.com/~vb *****

MPLEVL = =0 =2

chooses Moller-Plesset perturbation theory level, after the SCF. See $MP2, or $MRMP for MCSCF. skip the MP computation (default) perform second order energy correction.

MP2 (a.k.a. MBPT(2)) is implemented for RHF, UHF, ROHF, and MCSCF wavefunctions, but not GVB. Gradients are available for RHF, UHF, or ROHF based MP2, but for MCSCF, you must choose numerical derivatives to use any RUNTYP other than ENERGY, TRUDGE, SURFACE, or FFIELD.


Input Description

$CONTRL *****

2-8

CITYP

= = NONE = CIS

= SFCIS = ALDET = ORMAS = FSOCI = GENCI = GUGA

chooses CI computation after the SCF, for any SCFTYP except UHF. skips the CI. (default) single excitations from a SCFTYP=RHF reference, only. This is for excited states, with analytic nuclear gradients available. See the $CIS input group. spin-flip style CIS, see $CIS input. runs the Ames Laboratory determinant full CI package, requiring $CIDET. runs an Occupation Restricted Multiple Active Space determinant CI. The input is $CIDET and $ORMAS. runs a full second order CI using determinants, see $CIDET and $SODET. runs a determinant CI program that permits arbitrary specification of the determinants, requiring $CIGEN. runs the Unitary Group CI package, which requires $CIDRT input. Analytic gradients are available only for RHF, so for other SCFTYPs, you may choose only RUNTYP=ENERGY, TRUDGE, SURFACE, FFIELD, TRANSITN.

PMTD1

= For CITYP=ALDET or ORMAS, or for these two CI steps in MCSCF runs, for EFP solvent calculations, this flag enables use of "polarization method 1" for the effective fragments. See also FSTATE in $CIDET or $DET = .TRUE. The EFP dipoles will not be re-polarized to the CITYP wavefunction (default) = .FALSE. The EFP dipoles will be re-polarized to the CITYP wavefunction *****

CCTYP

chooses a Coupled-Cluster (CC calculation for the ground state and, optionally, Equation of Motion Coupled-Cluster (EOMCC) computation for excited states, both performed after the SCF (RHF or ROHF). See also $CCINP and $EOMINP.


Input Description

$CONTRL

2-9

Only CCSD and CCSD(T) for RHF can run in parallel. For ROHF, you may choose only CCSD and CR-CCL. = NONE = LCCD = CCD = CCSD = CCSD(T) = R-CC = CR-CC = CR-CCL skips CC computation (default). perform a coupled-cluster calculation using the linearized coupled-cluster method with double excitations. perform a CC calculation using the coupled-cluster method with doubles. perform a CC calculation with both single and double excitations. in addition to CCSD, the non-iterative triples corrections are computed, giving standard CCSD[T] and CCSD(T) energies. in addition to all CCSD(T) calculations, compute the renormalized R-CCSD[T] and R-CCSD(T) energies. in addition to all R-CC calculations, the completely renormalized CR-CCSD[T] and CR-CCSD(T) energies are computed. in addition to a CCSD ground state, the non-iterative triples energy correction defining the rigorously size extensive completely renormalized CR-CC(2,3), also called CR-CCSD(T)_L theory, is computed. Ground state only (zero NSTATE vector) CCTYP=CR-EOM type CR-EOMCCSD(T) energies and CCSD properties are also generated. For further information about accuracy, and A to D CR-CC(2,3) energy types, see REFS.DOC. in addition to all R-CC calculations, non-iterative triple and quadruple corrections are used, to give CCSD(TQ) and various R-CCSD(TQ) energies. in addition to all CR-CC and CCSD(TQ) calculations, the CR-CCSD(TQ) energies are obtained.

= CCSD(TQ)

= CR-CC(Q)

excited state options, note that EOM-CCSD is available for RHF or ROHF references, but triples corrections only for RHF cases. = EOM-CCSD in addition to a CCSD ground state, excited states are calculated using the equation of motion coupled-cluster


Input Description

$CONTRL

2-10

= CR-EOM

= CR-EOML

method with singles and doubles. in addition to the CCSD and EOM-CCSD, noniterative triples corrections to CCSD ground-state and EOM-CCSD excited-state energies are found, using completely renormalized CR-EOMCCSD(T) approaches. in addition to printing all results that CR-EOM obtains, this solves the lambda equations, and gives triples corrections analogous to ground state CR-CCL.

ionization processes, = IP-EOM2 ionized EOMCC with up to 2h1p excitations (i.e., IP-EOMCCSD) = IP-EOM3A ionized EOMCC with all 1h and 2h1p, and active-space 3h2p excitations (i.e., IP-EOMCCSDt) = EA-EOM2 electron-attached EOMCC with up to 2p1h excitations (i.e., EA-EOMCCSD) = EA-EOM3A electron-attached EOMCC with all 1p and 2p1h, and active-space 3p2h excitations (i.e., EA-EOMCCSDt). Labels "p" and "h" in the description of IP and EA EOMCC methods refer to particles (unoccupied correlated orbitals) and holes (occupied correlated orbitals). EA and IP runs produce both ground and excited states of systems obtained by attaching an electron to or removing an electron from the underlying CCSD reference ground state, using the EOMCC formalism. Thus, EA and IP runs read $CCINP as well as $EOMINP inputs. Any publication describing the results of CC calculations obtained using GAMESS should reference the appropriate papers, which are listed on the output of every run, and in chapter 4 of this manual. Analytic gradients are not available, so use CCTYP only for RUNTYP=ENERGY, TRUDGE, SURFACE, or maybe FFIELD, or request numerical derivatives. Generally speaking, the Renormalized energies are obtained at similar cost to the standard values, while Completely Renormalized energies cost twice the time. For usage tips and more information about resources on the various Coupled Cluster methods, see Section 4, 'Further Information'.


Input Description

$CONTRL

2-11

CIMTYP

chooses a Cluster-In-Molecule (CIM) calculation. = NONE skip CIM computation, i.e., perform a canonical calculation (default). = SECIM perform a single-environment CIM (SECIM) computation. = DECIM perform a dual-environment CIM (DECIM) computation. = GSECIM perform a generalized SECIM (GSECIM) computation. The $CIMFRG must be included as well. See also $CIMINP and, optionally, $CIMFRG and $CIMATM. If CIMTYP is given, SUBMTD in $CIMINP is required. Only RUNTYP=ENERGY and SCFTYP=RHF or ROHF work when CIMTYP is given. See SUBMTD in $CIMINP for more details. ***** RELWFN = NONE (default) See also the $RELWFN input group. = IOTC infinite-order two-component method of M. Barysz and A.J. Sadlej. This is the most sound scalar relativity choice, being at the infinite order instead of just third. = DK Douglas-Kroll transformation, available at the 1st, 2nd, or 3rd order. = RESC relativistic elimination of small component, the method of T. Nakajima and K. Hirao, available at 2nd order only. = NESC normalised elimination of small component, the method of K. Dyall, 2nd order only. ***** RUNTYP = ENERGY = GRADIENT = HESSIAN specifies the type of computation, for example at a single geometry point: Molecular energy. (default) Molecular energy plus gradient. Molecular energy plus gradient plus second derivatives, including harmonic harmonic vibrational analysis. See the $FORCE and $CPHF input groups. For FMO, use HESSFMO instead of HESSIAN. the same as HESSIAN, for FMO runs, supported only for RHF, R-DFT, UHF,

= FMOHESS


Input Description

$CONTRL

2-12

= GAMMA

U-DFT, and ROHF. Evaluate up to 3rd nuclear derivatives, by finite differencing of Hessians. See $GAMMA, and also NFFLVL in $CONTRL. multiple geometry options:

= OPTIMIZE = TRUDGE = SADPOINT = MEX = CONICAL = IRC = VSCF = DRC = MD = GLOBOP = OPTFMO = GRADEXTR = SURFACE

Optimize the molecular geometry using analytic energy gradients. See $STATPT. Non-gradient total energy minimization. See $TRUDGE and $TRURST. Locate saddle point (transition state). See $STATPT. Locate minimum energy crossing point on the intersection seam of two potential energy surfaces. See $MEX. Locate conical intersection point on the intersection seam of two potential energy surfaces. See $CONICL. Follow intrinsic reaction coordinate. See $IRC. anharmonic vibrational corrections. See $VSCF. Follow dynamic reaction coordinate. See $DRC. molecular dynamics trajectory, see $MD. Monte Carlo-type global optimization. See $GLOBOP. genuine FMO geometry optimization using nearly analytic gradient. See $OPTFMO. Trace gradient extremal. See $GRADEX. Scan linear cross sections of the potential energy surface. See $SURF. single geometry property options:

= COMP = G3MP2 = PROP

composite thermochemistry calculation, including G3MP2. See $COMP input. evaluate heat of formation using the G3(MP2,CCSD(T)) methodology. See test example exam43.inp for more information. Molecular properties will be calculated. Orbital localization can be requested as well. See $ELPOT, etc. Converged orbitals must be input in a


Input Description

$CONTRL

2-13

= RAMAN = NACME

= NMR = EDA = QMEFPEA = TRANSITN = FFIELD = TDHF = TDHFX = MAKEFP = FMO0

$VEC input, which suffice to reproduce the wavefunction only for simple SCF: RHF, UHF, ROHF, or DFT counterparts. GVB also works (CICOEF may be needed). All other calculations must instead use RUNTYP=ENERGY to regenerate the density matrix. computes Raman intensities, see $RAMAN. non-adiabatic coupling matrix element between two or more state averaged MCSCF wavefunctions. The calculation has no specific input group, but must use only SCFTYP=MCSCF with CISTEP=ALDET or ORMAS. NMR shielding tensors for closed shell molecules by the GIAO method. See $NMR. Perform energy decomposition analysis. Give one of $MOROKM or $LMOEDA inputs. QM/EFP solvent energy analysis, see $QMEFP. Compute radiative transition moment or spin-orbit coupling. See $TRANST. applies finite electric fields, most commonly to extract polarizabilities. See $FFCALC. analytic computation of time dependent polarizabilities. See $TDHF. extended TDHF package, including nuclear polarizability derivatives, and Raman and Hyper-Raman spectra. See $TDHFX. creates an effective fragment potential, for SCFTYP=RHF or ROHF only. See $MAKEFP, $DAMP, $DAMPGS, $STONE, ... performs the free state FMO calculation. See $FMO.

***************************** Note that RUNTYPs which require the nuclear gradient are GRADIENT, HESSIAN, OPTIMIZE, SADPOINT, GLOBOP, IRC, GRADEXTR, DRC, and RAMAN These are efficient with analytic gradients, which are available only for certain CI or MP2 calculations, but no CC calculations, as indicated above. See NUMGRD. ***************************** NUMGRD Flag to allow numerical differentiation


Input Description

$CONTRL

2-14

of the energy. Each gradient requires the energy be computed twice (forward and backward displacements) along each totally symmetric modes. It is thus recommended only for systems with just a few symmetry unique atoms in $DATA. The default is .FALSE. EXETYP = RUN = CHECK Actually do the run. (default) Wavefunction and energy will not be evaluated. This lets you speedily check input and memory requirements. See the overview section for details. Note that you must set PARALL=.TRUE. in $SYSTEM to test distributed memory allocations. Massive amounts of output are printed, useful only if you hate trees. Maximum output is generated by the routine named. Check the source for the routines this applies to. ******* ICHARG = MULT = =1 = 2,3,... Molecular charge. (default=0, neutral)

= DEBUG = routine

Multiplicity of the electronic state singlet (default) 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. * * * 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. = PRINAXIS Cartesian coordinates will be input, and transformed to principal axes. Please read the warning just below!!!


Input Description = ZMT = ZMTMPC = FRAGONLY

$CONTRL

2-15

GAUSSIAN style internals will be input. MOPAC style internals will be input. means no part of the system is treated by ab initio means, hence $DATA is not given. The system is defined by $EFRAG.

Note: the choices PRINAXIS, ZMT, ZMTMPC require input of all atoms in the molecule. They also orient the molecule, and then determine which atoms are unique. The reorientation is likely to change the order of the atoms from what you input. When the point group contains a 3fold 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: the choices PRINAXIS, ZMT, ZMTMPC require the use of a group named $BASIS 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 =0 =M Use Cartesian coordinates (default). If COORD=ZMT or ZMTMPC, and $ZMAT is not given: the internal coordinates will be those defining the molecule in $DATA. In this case, $DATA may not contain any dummy atoms. M is usually 3N-6, or 3N-5 for linear. 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

=M

NZVAR refers mainly to the coordinates used by OPTIMIZE or SADPOINT runs, but may also print the internal's


Input Description

$CONTRL

2-16

values for other run types. You can use internals to define the molecule, but Cartesians during optimizations! ******* Pseudopotentials may be of two types: ECP (effective core potentials) which generate nodeless valence orbitals, and MCP (model core potentials) producing valence orbitals with the correct radial nodal structure. At present, ECPs have analytic nuclear gradients and Hessians, while MCPs have analytic nuclear gradients. PP = = NONE = READ = SBKJC = HW = MCP pseudopotential selection. all electron calculation (default). read ECP potentials in the $ECP input. use Stevens, Basch, Krauss, Jasien, Cundari ECP potentials for all heavy atoms (Li-Rn are available). use Hay, Wadt ECP potentials for heavy atoms (Na-Xe are available). use Huzinaga's Model Core Potentials. The correct MCP potential will be chosen to match the requested MCP valence basis set (see $BASIS). ******* LOCAL controls orbital localization. NONE Skip localization (default). BOYS Do Foster-Boys-like localization. RUEDNBRG Do Edmiston-Ruedenberg localization. POP Do Pipek-Mezey population localization. SVD Do single value decomposition, to project the molecular orbitals onto atoms. This is available only for SCFTYP=RHF, ROHF, and MCSCF (full space or ORMAS). The ORIENT keyword in $LOCAL is pertinent. See the related $LOCAL input. Localization is not available for SCFTYP=GVB. DFTB only works with LOCAL=POP (and NONE). ******* ISPHER = = -1 Spherical Harmonics option Use Cartesian basis functions to construct = = = = = =


Input Description

$CONTRL symmetry-adapted linear combination (SALC) of basis functions. The SALC space is the linear variation space used. (default) 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. 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.

2-17

=0

= +1

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. MAXIT = Maximum number of SCF iteration cycles. This pertains only to RHF, UHF, ROHF, or GVB runs. See also MAXIT in $MCSCF. (default = 30) * * * interfaces to other programs * * * MOLPLT = flag that produces an input deck for a molecule drawing program distributed with GAMESS. (default is .FALSE.)


Input Description

$CONTRL

2-18

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 http://www.chemistry.mcmaster.ca/aimpac 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 OPTIMIZE, SAPDOINT, IRC... 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. *** NFFLVL used to determine energies and gradients away from equilibrium structures, at the coordinates given in $DATA. The method will use a Taylor expansion of the potential surface around the stationary point. See $EQGEOM, $HLOWT, $GLOWT. This may be used with RUNTYP=ENERGY or GRADIENT. = 2 uses only Hessian information, which gives a reasonable energy, but not such a good gradient. = 3 uses Hessian and 3rd nuclear derivatives in the Taylor expansion, producing more accurate values for the energy and for the gradient.


Input Description

$CONTRL

2-19

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

=0 =1

ETOLLZ

= threshold to label molecular orbitals by Lz values. Small matrices of the Lz operator are diagonalized for the sets of MOs whose orbital


Input Description

$CONTRL

2-20

energies are degenerate to within ETOLLZ. This option may be used in molecules with distorted linear symmetry for approximate labelling. Default: 1.0d-6 for linear, 0 (disable) if not. INTTYP selects the integral package(s) used, all of which produce equally accurate results. This is therefore used only for debugging purposes. = BEST use the fastest integral code available for any particular shell quartet (default): s,p,L or s,p,d,L rotated axis code first. ERIC s,p,d,f,g precursor transfer equation code second, up to 5 units total ang. mom. Rys quadrature for general s,p,d,f,g,L, or for uncontracted quartets. = ROTAXIS means don't use ERIC at all, e.g. rotated axis codes, or else Rys quadrature. = ERIC means don't use rotated axis codes, e.g. ERIC code, or else Rys quadrature. = RYSQUAD means use Rys quadrature for everything. GRDTYP = BEST use Schlegel routines for spL gradient blocks, and Rys quadrature for all other gradient integrals. (default) = RYSQUAD use Rys quadrature for all gradient integrals. This option is only slightly more accurate, but is rather slower. normalize the basis functions (default) no normalization input contraction coefficients refer to normalized Gaussian primitives. (default) the opposite. primitive cutoff factor (default=20) products of primitives whose exponential factor is less than 10**(-n) are skipped. integrals less than 10.0**(-n) are not saved on disk. (default = 9). Direct SCF will calculate to a cutoff 1.0d-10 or 5.0d-11 depending on FDIFF=.F. or .T. proceed as usual

NORMF NORMP

=0 =1 =0 =1

ITOL

= =n =n

ICUT

ISKPRP = 0


Input Description 1

$CONTRL

2-21

skip computation of some properties which are not well parallelised. This includes bond orders and virial theorem, and can help parallel scalability if many CPUs are used. Note that NPRINT=-5 disables most property computations as well, so ISKPRP=1 has no effect in that case. (default: 0) * * * restart options * * *

IREST

= = -1 = = = = = 0 1 2 3 4

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


Input Description

$SYSTEM

2-22

==========================================================

$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! MWORDS = the maximum replicated memory which your job can use, on every core. This is given in units of 1,000,000 words (as opposed to 1024*1024 words), where a word is defined as 64 bits. (default=1)

In case finer control over the replicated memory is needed, this value can be given in units of words, with the old 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 processor core for a run using p cores is therefore MEMDDI/p + MWORDS. The parallel runs that currently require MEMDDI are: SCFTYP=RHF MPLEVL=2 energy or gradient SCFTYP=UHF MPLEVL=2 energy or gradient SCFTYP=ROHF MPLEVL=2 OSPT=ZAPT energy or gradient SCFTYP=MCSCF MPLEVL=2 energy SCFTYP=MCSCF using the FULLNR or JACOBI convergers SCFTYP=MCSCF analytic hessian SCFTYP=any CITYP=ALDET, ORMAS, GUGA SCFTYP=any energy localization SCFTYP=RHF CCTYP=CCSD or CCSD(T) All other parallel runs should enter MEMDDI=0, for they use only replicated memory. Some serial runs execute the parallel code (on just 1 CPU), for there is only a parallel code. These serial runs must give MEMDDI as a result: SCFTYP=ROHF MPLEVL=2 OSPT=ZAPT gradient/property run SCFTYP=MCSCF analytic hessian Two kinds of runs (RI-MP2 and parallel CCSD(T)) use an additional type of memory, for which there is no input


Input Description

$SYSTEM

2-23

keyword. Please read EXETYP=CHECK output carefully to learn the total memory/node requirements for these two! TIMLIM = time limit, in minutes. Set to about 95 percent of the time limit given to the batch job (if you use a queueing system) so that GAMESS can stop itself gently. (default=525600.0 minutes) a flag to cause the program to execute the parallel algorithm, in cases where different serial and parallel codes exist, if you happen to be running on only one core. The default is .FALSE. if you are running on one core. The main purpose of this keyword is to allow you to do EXETYP=CHECK runs on only one core, when your intent is perform the actual calculation in parallel. PARALL is ignored for runs on more than one core, when of course parallel algorithms are executed. diagonalization control switch use a vectorized diagonalization routine if one is available on your machine, else use EVVRSP. (default) use EVVRSP diagonalization. This may be more accurate than KDIAG=0. use GIVEIS diagonalization (not as fast or reliable as EVVRSP) use JACOBI diagonalization (this is the slowest method)

PARALL =

KDIAG

= =0 =1 =2 =3

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.)

BALTYP = Parallel load balance scheme: = SLB uses static load balancing. = DLB uses dynamic load balancing (default). Dynamic load balancing attempts to spread out possibly unequal work assignments based on the rate at which different nodes complete tasks. For historical reasons, it is permissible to spell SLB as LOOP, and DLB as NXTVAL.


Input Description

$SYSTEM

2-24

MXSEQ2 = 300 (default) MXSEQ3 = 150 (default) Matrix/vector problem size in loops requiring either O(N**2) or O(N**3) work, respectively. Problems below these sizes are run purely serial, to avoid poor communication/computation ratios. NODEXT = array specifying node extensions in GDDI for each file. Non-zero values force no extension. E.g., NODEXT(40)=1 forces file 40 (file numbers are unit numbers used in GAMESS, see "rungms" or PROG.DOC) to have the name of $JOB.F40 on all nodes, rather than $JOB.F40, $JOB.F40.001, $JOB.F40.002 etc. This is convenient for FMO restart jobs, so that the file name need not be changed for each node, when copying the restart file. Note that on machines when several CPUs use the same directory (e.g., SMP) NODEXT should be zero. (default: all zeros) IOSMP = Parallelise I/O on SMP machines with multiple hard disks. Two parameters are specified, whose meaning should be clear from the example. iosmp(1)=2,6 2 refers to the number of HDDs per SMP box. 6 is the location of the character in the file names that switches HDDs, i.e. if HDDs are mounted as /work1 and /work2, then 6 refers to the position of the number 1 in /work1. The file system should permit disks attached with directory names differing by one symbol. (default: 0,0, disable the feature) = Global I/O options (bitwise additive) (default: 0) 1 - forbid flushing files 2 - do not close dictionary file in GDDI (only record indices are reset) 4 - do not print timings on each rank at run end. 8 - forbid grid data saving in DFT (prevent F22 from being opened) 16 - reduce I/O in MD. 32 - do not open file F15 (in SOSCF) on slaves and do not close on all.

MODIO


Input Description

$SYSTEM

2-25

64 - use in-memory F15 (in SOSCF). This also parallelizes one more step in SOSCF. 128 - always run EVVRSP sequentially. This is useful on mixed CPU type clusters. MODEIO=15 is recommended on computers w/o local disks. MEM10 = words used to store dictionary file F10 in memory. Selecting this option will skip any I/O for F10. Default: 0 (disk-based F10)

==========================================================


Input Description

$BASIS

2-26

==========================================================

$BASIS

group

(optional)

This group allows certain standard basis sets to be easily requested. Basis sets are specified by: a) b) c) d) GBASIS plus optional supplementations such as NDFUNC, BASNAM to read custom basis sets from your input, EXTFIL to read custom bases from an external file, or omit this group entirely, and give the basis set in the $DATA input, which is completely general.

GBASIS requests various Gaussian basis sets. These include options for effective core and model core potentials. Rather oddly, GBASIS also can select semi-empirical models, and in that case requests the Slater-type orbitals for the MOPAC-type calculation. Note: The first two groups of GBASIS keywords below (except G3L and G3LX) define only the basic functions, without any polarization functions and/or diffuse functions. For example, main group elements have the basic functions for their s,p valence orbitals. Polarization and/or diffuse supplements are added separately to these GBASIS values, with keywords NPFUNC, NDFUNC, NFFUNC, DIFFS, DIFFSP, POLAR, SPLIT2, and SPLIT3, which are defined at the end of this input group. GBASIS = 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-Kr, for NGAUSS=6, note that the bases for K,Ca,Ga-Kr were changed 9/2006. = N311 - Pople's "triple split" N-311G basis set. Available H-Ne, for NGAUSS=6. Selecting N311 implies MC for Na-Ar. = G3L - Pople's G3MP2Large basis set, for H-Kr.


Input Description

$BASIS

2-27

= G3LX - Pople's G3MP2LargeXP basis set, for H-Kr. NGAUSS = the number of Gaussians (N). This parameter pertains to GBASIS=STO, N21, N31, or N311. GBASIS = MINI - Huzinaga's 3 gaussian minimal basis set. Available H-Rn. = MIDI - Huzinaga's 21 split valence basis set. Available H-Rn. = 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. * * systematic basis set families * * * These four families provide a hierachy of basis sets approaching the complete basis set limits. These families include relevant polarization and diffuse augmentations, as indicated in their names. GBASIS = CCn - Dunning-type Correlation Consistent basis sets, officially called cc-pVnZ. Use n = D,T,Q,5,6 to indicate the level of polarization. Available for H-He, Li-Ne, Na-Ar, Ca, Ga-Kr and for Sc-Zn for n=T,Q. = ACCn - As CCn, but augmented with a set of diffuse functions, e.g. aug-cc-pVnZ. Availability is the same as CCn.


Input Description

$BASIS

2-28

= CCnC - As CCn, but augmented with tight functions for recovering core and core-valence correlation, e.g. cc-pCVnZ. = ACCnC- As CCn, augmented with diffuse as well as CCnC's tight functions, e.g. aug-cc-pCVnZ. Availability is the same as CCnC. = CCnWC the omega form of CCnC, e.g. cc-pwCVnZ, for H-Ar, for n=T only. CCnWC's tight functions are considered superior to CCnC's for recovery of core/valence correlation. = ACCnWC augmented form of CCnWC: aug-cc-pwCVnZ. See extended notes below! = PCseg-n Polarization Consistent basis sets. n = 0,1,2,3,4 indicates the level of polarization. (n=0 is unpolarized, n=1 is ~DZP, n=2 is ~TZ2P, etc.). These provide a hierarchy of basis sets suitable for DFT and HF calculations. Available for H-Kr. = APCseg-n - These are the PCseg-n bases, with diffuse augmentation. See extended notes below! = SPK-nZP Sapporo valence basis sets: Sapporo family of non-relativistic bases, n=D,T,Q, available H-Xe diffuse augmentation of the above. Sapporo family of relativistic bases n=D,T,Q, available H-Xe. These should be used only with a relativistic transformation of the integrals, such as RELWFN=IOTC or RELWFN=DK. These sets are identical to SPK-nZP up to Ar. diffuse augmentation of the above. See extended notes below!

Available H-Ar for n=D,T,Q, also n=5 for H-Ne.

= SPK-AnZP = SPKrnZP -

= SPKrAnZP -

Sapporo core/valence basis sets: - Sapporo family of non-relativistic bases, n=D,T,Q, available H-Xe = SPK-nZCD - diffuse augmentation of the above. = SPKrnZC - Sapporo family of relativistic bases n=D,T,Q, available H-Xe, and La-Lu. To be used only with a relativistic = SPK-nZC


Input Description

$BASIS

2-29

transformation of the integrals, such as RELWFN=IOTC or RELWFN=DK. These sets are identical to SPK-nZC up to Ar. = SPKrnZCD - diffuse augmentation of the above. See extended notes below! = KTZV - Karlsruhe valence triple zeta basis, as developed by Prof.Ahlrichs, see REFS.DOC. = KTZVP- Karlsruhe valence triple zeta basis with a set of single polarization (P). = KTZVPP-Karlsruhe valence triple zeta basis with a set of double polarization (PP). The Karlsruhe sets are provided for H-Ar. Normally these families are used as spherical harmonics, see ISPHER=1 in $CONTRL. Failure to set ISPHER=1 will result in discrepancies in energy values compared to the literature or other programs, difficulties in converging SCF/DFT, CC, CI, and/or response equation iterations, and longer run times due to retention of unimportant MOs. The calculations will refuse to run without ISPHER being set. Important note about the PCseg basis set family: 1. These should be used only in spherical harmonic form. 2. The PCn basis sets included in GAMESS versions prior to March 2014 were generally contracted, but were replaced by computationally more efficient segmented contractions, and renamed to PCseg-n. The segemented contractions have the same or slightly better accuracy (especially for n=0) as the original PCn bases, which are no longer available. Important notes about the CC basis set family: 1. These should be used only in spherical harmonic form. 2. The CC5 and CC6 basis sets (and corresponding augmented versions) contain h-functions, and CC6 also contains ifunctions. As of January 2013, GAMESS integral code can correctly use h & i functions, so these three call up the true basis sets. Prior to January 2013, GAMESS' integral codes was restricted to g-functions, so these three truncated away any h & i functions, to spdfg subsets, and therefore were not the true basis sets. 3. Note that the CC basis sets are generally contracted, which GAMESS can only handle by replicating the primitive


Input Description

$BASIS

2-30

basis functions, leading to a less than optimum performance in AO integral evaluation. 4. The implementation of the cc-pVnZ and cc-pCVnZ basis sets for Na-Ar include one additional tight d-function, producing the so-called cc-pV(n+d)Z and cc-pV(n+d)Z sets, which are known to improve results (see J.Chem.Phys. 114, 9244(2001) and Theoret.Chem.Acc. 120, 119(2008)). These tight d versions are invoked by GBASIS=CCn or CCnC (and also their augmented counterparts ACCn or ACC). This means the old (and less accurate) basis sets without the tight d's are not available for Na-Ar. 5. Alkali and alkali earth basis sets (Li,Be,Na,Mg) were changed April 2013 so that regular, diffuse, tight d (for Na/Mg), core/valence, and weighted core/valence sets agree with their official publication: Theor. Chem. Acc. 128, 69(2011). 6. In case you are interested in scalar relativistic effects, the CCT-DK and CCQ-DK sets optimized for use with Douglas/Kroll are available for Sc-Kr. These will be used if you type GBASIS=CCT or CCQ along with RELWFN=DK or IOTC, using NR sets for elements lighter than Sc. DK versions of ACCD or ACCT are available for Sc-Zn (but not Ga-Kr). Notes about the Sapporo basis set family: 1. SPK is the international airport city code for Sapporo. 2. These should be used only in spherical harmonic form. 3. The relativistic bases were optimized at the 3rd order of the Douglas-Kroll transformation, with a Gaussian nuclei model, but can also be used with the infinite order twocomponent scheme IOTC (see RELWFN in $CONTRL). It would be very illogical to use these all-electron relativistic bases without also turning on scalar relativity! 4. Because they are stored in an external file supplied with GAMESS, these can only be accessed via GBASIS in this group, not by using them in-line in $DATA. 5. The core/valence basis sets treat (n-1)s,(n-1)p,ns for s-block elements; (n-1)s,(n-1)p,ns,np for p block elements; (n-1)s,(n-1)p,(n-1)d,ns for d block elements, and the 4s4f,5s-5d,6s for f block (lanthanides). Usage of these suggests you should change the default number of core orbitals, such as NACORE in $MP2 or NCORE in $CCINP. 6. The basis sets were extracted from the data base of Segmented Gaussian Basis Sets, maintained by Takeshi Noro, Quantum Chemistry Group, Sapporo, Japan:


Input Description

$BASIS

2-31

http://setani.sci.hokudai.ac.jp/sapporo/Welcome.do The mapping between the data base names and the keywords used in GAMESS is (for n=D,T,Q): data base name keyword Sapporo-nZP SPK-nZP Sapporo-nZP+diffuse SPK-AnZP Sapporo-DK-nZP SPKrnZP Sapporo-DK-nZP+diffuse SPKrAnZP Sapporo-nZP-2012 SPK-nZC Sapporo-nZP-2012+diffuse SPK-nZCD Sapporo-DK-nZP-2012 SPKrnZC Sapporo-DK-nZP-2012+diffuse SPKrnZCD

* * * Effective Core Potential (ECP) bases * * * 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. * * * Model Core Potential (MCP) bases * * * Notes: Select PP=MCP in $CONTRL to automatically use the model core potential matching your basis choice below. References for these bases, and other information about MCPs can be found in the REFS.DOC chapter. Another family covering almost all elements is available in $DATA only. GBASIS = MCP-DZP, MCP-TZP, MCP-QZP a family of double, triple, and quadruple zeta quality valence basis sets, which are akin to the correlation consistent sets, in that these include increasing levels of polarization (and so do not require "supplements" like NDFUNC or DIFFSP) and must be used as spherical harmonics (see ISPHER). Availability: MCP-DZP: 56 elements Z=3-88, except V-Zn, Y-Cd, La, Hf-Hg MCP-TZP, MCP-QZP: 85 elements Z=3-88, except La


Input Description

$BASIS

2-32

The basis sets for hydrogen atoms will be the corresponding Dunning's cc-pVNZ (N=D,T,Q). = MCP-ATZP, MCP-AQZP MCP-TZP and MCP-QZP core potentials whose basis sets were augmented with diffuse functions Availability: same as for MCP-TZP, MCP-QZP = MCPCDZP, MCPCTZP, MCPCQZP based on MCP-DZP, MCP-TZP, MCP-QZP, with core-valence functions provided for the alkali and alkaline earth atoms Na through Ra. = MCPACDZP, MCPACTZP, MCPACQZP based on MCPCDZP, MCPCTZP, MCPCQZP, with core-valence functions provided for the alkali and alkaline earth atoms Na through Ra, and augmented with diffuse functions. The basis sets were extracted from the data base Segmented Gaussian Basis Sets, maintained by Takeshi Noro, Quantum Chemistry Group, Sapporo, Japan: http://setani.sci.hokudai.ac.jp/sapporo/Welcome.do The mapping between the data base names and the names used in GAMESS is data base name GAMESS keyword MCP/NOSeC-V-DZP MCP/NOSeC-V-TZP MCP/NOSeC-V-QZP MCP/NOSeC-V-TZP+diffuse MCP/NOSeC-V-QZP+diffuse MCP/NOSeC-CV-DZP MCP/NOSeC-CV-TZP MCP/NOSeC-CV-QZP MCP/NOSeC-CV-DZP+diffuse MCP/NOSeC-CV-TZP+diffuse MCP/NOSeC-CV-QZP+diffuse MCP-DZP MCP-TZP MCP-QZP MCP-ATZP MCP-AQZP MCPCDZP MCPCTZP MCPCQZP MCPACDZP MCPACTZP MCPACQZP

GBASIS = IMCP-SR1 and IMCP-SR2 valence basis sets to be used with the improved MCPs with scalar relativistic effects.


Input Description

$BASIS

2-33

These are available for transition metals except La, and the main group elements B-Ne, P-Ar, Ge, Kr, Sb, Xe, Rn. The 1 and 2 refer to addition of first and second polarization shells, so again don't use any of the "supplements" and do use spherical harmonics. = IMCP-NR1 and IMCP-NR2 closely related valence basis sets, but with nonrelativistic model core potentials. GBASIS = ZFK3-DK3, ZFK4-DK3, ZFK5-DK3, or ZFK3LDK3, ZFK4LDK3, ZFK5LDK3 These are a family of model core potential basis sets developed by Zeng/Fedorov/Klobukowski, for the p-block elements from 2p to 6p. The potentials were paramaterized taking into account both DK3 scalar relativistic and DK-SOC effects. The fundamental basis functions are from the Well-Tempered Basis Sets. The number after ZFK indicates the augmentation levels, e.g. ZFK3 means the diffuse functions from aug-cc-pVTZ are added, ZFK4 means from augcc-pVQZ, etc. The difference between ZFKn-DK3 and ZFKnLDK3 is that the common s and p exponents have been contracted as a single L-shell for the outermost s and p valence shells to save time in the "L" case. The s-block elements from 1s to 4s have also been put in the library. For H/He, all-electron aug-cc-pVnZ basis sets are used. For Li/Be, the relativistically contracted atomic natural orbital allelectron basis sets (ANO-RCC) are used. For Na/Mg, and K/Ca, unpublished MCP and basis sets based on ANO-RCC are available, although the potentials have not been extensively tested yet. No d-block elements can be used. * * * semiempirical basis sets * * * GBASIS = = = = = MNDO AM1 PM3 RM1 DFTB selects selects selects selects selects MNDO model Hamiltonian AM1 model Hamiltonian PM3 model Hamiltonian RM1 model Hamiltonian tight binding Hamiltonian

Note: 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 Slater type


Input Description

$BASIS

2-34

orbitals (STOs), not Gaussian GTOs, but the keyword remains GBASIS. Except for NGAUSS, all other keywords such as NDFUNC, etc. will be ignored for these. If you add NGAUSS, STO-NG expansions of the valence STO functions in terms of Gaussians will be added to the log file. Plotting programs such as MacMolPlt can pick up this approximation to the STOs used up from the ouput, in order to draw the orbitals. The default NGAUSS=0 suppresses this output, but values up to 6 may be given to control the accuracy of the STO-NG printing. --- supplementary functions --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 in the $CONTRL input group for numerical concerns. POLAR = = = = exponent of polarization functions COMMON (default for GBASIS=STO,N21,HW,SBKJC) POPN31 (default for GBASIS=N31) POPN311 (default for GBASIS=N311, MC)


Input Description = DUNNING = HUZINAGA = HONDO7

$BASIS (default for GBASIS=DH, DZV) (default for GBASIS=MINI, MIDI) (default for GBASIS=TZV)

2-35

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 hydrocarbons. 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 input group. Thus, if $BASIS does not quite cover the basis you want, you can obtain this full $DATA from EXETYP=CHECK, and then change polarization exponents, add Rydbergs, etc. *** This may only be used with COORD=UNIQUE or HINT! BASNAM = an array of names of customized basis set input groups. BASNAM should obey the rule of no more than six characters starting with a letter names, and must avoid using any GBASIS string. However, the individual basis inputs can use any of the GBASIS sets by its standard name. Basis supplementations such as DIFFS or NDFUNC may only be given by explicit numerical values. This is best explained by an example where a core potential and valence-only basis set is used on a transition metal, but not its ligands:


Input Description

$BASIS

2-36

$contrl scftyp=rohf icharg=+3 mult=4 runtyp=gradient pp=read ispher=1 $end $system mwords=1 $end $guess guess=huckel $end $basis basnam(1)=metal, ligO,ligO,ligO,ligO,ligO,ligO, ligH,ligH,ligH,ligH,ligH,ligH, ligH,ligH,ligH,ligH,ligH,ligH $end $data Cr+3(H2O)6 complex...SBKJC & 6-31G(d) geometry Th CHROMIUM 24.0 .0000000000 .0 .0000000000 OXYGEN 8.0 .0000000000 .0 2.0398916104 HYDROGEN 1.0 .7757887450 .0 2.6122732372 $end ! core potential basis for Chromium $metal sbkjc $end ! normal 6-31G(d) for oxygen ligands $ligO n31 6 d 1 ; 1 0.8 1.0 $end ! $ligH n31 6 unpolarized basis for hydrogens

$end $ecp Cr-ecp SBKJC O-ecp none ...snipped... O-ecp none H-ecp none ...snipped... H-ecp none $end

there must be 6 O's given here there must be 12 H's given here

*** EXTFIL = a flag to read basis sets from an external file,


Input Description

$BASIS

2-37

defined by EXTBAS, rather than from $DATA. (default=.false.) It may be easier to use BASNAM to create custom basis sets! BASNAM has the bonus that your input file contains all information about the calculation, explicitly. Except for MCP basis sets, no external file is provided with GAMESS, thus you must create your own. The GBASIS keyword must give an 8 or less character string, obviously not using any internally stored names. Every atom must be defined in the external file by a line giving the chemical symbol, and this chosen string. Following this header line, give the basis in free format $DATA style, containing only S, P, D, F, G, and L shells, and terminating each atom by the usual blank line. The external file may have several families of bases in the same file, identified by different GBASIS strings. ===========================================================


Input Description

$DATA

2-38

==========================================================

$DATA group $DATAS group $DATAL group

(required) (if NESC chosen, for small component basis) (if NESC chosen, for 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 information on large and small component basis sets. The input structure of $DATAS and $DATAL is identical to the COORD=UNIQUE $DATA input. ----------------------------------------------------------1TITLE a single descriptive title card.

----------------------------------------------------------2GROUP, 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. Use of NAXIS up to 8 is supported in each axial groups. 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.


Input Description

$DATA

2-39

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 wavefunctions. Cards -3- and -4- are quite complicated, and are rarely given. A *SINGLE* blank card may replace both cards -3and -4-, to select the 'master frame', which is defined on the next page. If you choose to enter a blank line, skip to one of the -5- input sequences. Note! If the point group is C1 (no symmetry), skip over cards -3- and -4- (which means no blank card). ----------------------------------------------------------3For For For For For For X1, Y1, Z1, X2, Y2, Z2 C1 group, there is no card -3- or -4-. CI group, give one point, the center of inversion. CS group, any two points in the symmetry plane. axial groups, any two points on the principal axis. tetrahedral groups, any two points on a two-fold axis. octahedral groups, any two points on a four-fold axis.

----------------------------------------------------------4X3, 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).


Input Description

$DATA

2-40

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,


Input Description if if if if if COORD=UNIQUE, COORD=HINT, COORD=CART, COORD=ZMT, COORD=ZMTMPC, follow follow follow follow follow

$DATA card card card card card sequence sequence sequence sequence sequence U U C G M

2-41

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. ----------------------------------------------------------5UAtom 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 NAME, ZNUC, X, Y, Z

= 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,ZNUC,CONX,R,ALPHA,BETA,SIGN,POINT1,POINT2,POINT3 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


Input Description

$DATA

2-42

(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 -5Uuntil 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 can have exactly the same meaning as the keyword in $BASIS. You may choose from STO, N21, N31, N311, ACCT, PC4, ... A few of these require NGAUSS below. In addition, GBASIS can be S, P (or L), D, F, G, H, or I to enter an explicit basis set. L means both an S and P shell with the same exponent. See NGAUSS below. In addition, GBASIS may be defined as MCP, to indicate that the current atom is represented by a model core potential, and valence basis set. An internally stored basis and potential will be applied (see REFS.DOC for the details). The MCP basis supplies only the occupied atomic orbitals, e.g. sp for a main group element, so please supplement with any desired polarization. In case the keyword MCP is followed by the keyword READ, everything will be taken from the input file, namely the basis functions are read using the sequence -6U-, -7U-, and -8U-, from lines following the "MCP READ" line. In addition, "MCP READ" implies that the parameters of the model core potentials, together with core basis functions are in the input stream, in a $MCP input group. Other MCP bases are available in the $BASIS input, but note that to locate the MCP, the atom name must be a chemical symbol, that is "P" instead of "Phosphorus".


Input Description

$DATA

2-43

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, or explicit GTO types such as S,P,D,F... Up to 4 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,... either add more shells by repeating card -6U-, or go on to -8U-. * * * If GBASIS=S,P,D,F,... 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-. ----------------------------------------------------------8UA 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


Input Description

$DATA

2-44

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)

----------------------------------------------------------5GATOM

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

Only a name and a bond distance is required for atom 2. See -8G- for a description of this information. ----------------------------------------------------------7GATOM i1 BLENGTH i2 ALPHA

Only a name, distance, and angle are required for atom 3. See -8G- for a description of this information. ----------------------------------------------------------8GATOM ATOM i1 BLENGTH i2 ALPHA i3 BETA i4 It can be

is the chemical symbol of this atom.


Input Description

$DATA

2-45

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 i1atom 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 o o Repeat -8G- for atoms 4, 5, ... 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. 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! 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 preceded 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.

o

----------------------------------------------------------9GA blank line terminates the Z-matrix section.

----------------------------------------------------------10GSTRING VALUE

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


Input Description

$DATA

2-46

---------------------------------------------------------Continue entering -10G- until all STRINGs are defined. Note that any blank card encountered while reading -10Gwill 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)

----------------------------------------------------------5MATOM

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

Only a name and a bond distance is required for atom 2. See -8M- for a description of this information. ----------------------------------------------------------7MATOM BLENGTH j1 ALPHA j2


Input Description

$DATA

2-47

Only a bond distance from atom 2, and an angle with respect 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. ----------------------------------------------------------8MATOM 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 out!


Input Description

$ZMAT

2-48

==========================================================

$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.


Input Description

$ZMAT

2-49

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 AUTO is a flag to request delocalized coordinates. (default is .FALSE.) 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. Cases where the AUTO generation of DLC coordinates fails to find the full set of 3N-6 coordinates typically fall 6 short of 3N-6. These cases are invariably due to the system being divided into pieces too far apart to have bonds detected (for example, system A might be H-bonded to system B, finding 3N-12 coordinates only). Adding NONVDW input for that H-bond will tie A and B together, and result in a correct AUTO generation of all 3N-6 coordinates. Falling short by an integer multiple of 6 indicates more than two pieces, requiring several NONVDW pairs. Falling short by 3 coordinates indicates one of the separate systems A or B is likely a single atom, which has no rotational degrees of freedom, again it should be attached by NONVDW. DLC coordinate generation can be fine tuned by IXZMAT, IRZMAT, IFZMAT whose format is the same as IZMAT:


Input Description IXZMAT is an extra array which you want to by AUTO. Unlike coordinate(s) you

$ZMAT

2-50

of simple internal coordinates have added to the list generated NONVDW, IXZMAT will add only the 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, which is --required-- input when IFZMAT is given. IFZMAT/FVALUE work with ordinary coordinate input using IZMAT, as well as with DLC, but in the former case be careful that IFZMAT specifies coordinates that were already given in IZMAT. In addition, IFZMAT works only for IZMAT=1,2,3 type coordinates. See IFREEZ in $STATPT you wish to freeze regular or natural internal 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


Input Description

$ZMAT

2-51

gives the matrix S=

1.0 0.0 1.0

1.0 0.0 -1.0

0.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) ==========================================================


Input Description

$LIBE

2-52

==========================================================

$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 appropriate. ==========================================================


Input Description

$SCF

2-53

==========================================================

$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 given to define GVB open shell or perfect pairing wavefunctions. See $MCSCF for multireference inputs. 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. The default of .FALSE. stores integrals on disk storage for a conventional SCF calculation. 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, or large molecules, may sometimes be "mushy" at the end, rather than converging. Increasing ICUT in $CONTRL by one may help this, or consider turning this parameter off. ---- The next flags affect convergence rates. NOCONV = .TRUE. means neither SOSCF nor DIIS will be used. The default is .FALSE., making the choice of the primary converger as follows: for RHF, GVB, UHF, or ROHF (if Abelian): SOSCF for any DFT, or for non-Abelian groups: DIIS. DIIS = selects Pulay's DIIS interpolation. SOSCF = selects second order SCF orbital optimization. Once either DIIS or SOSCF are initiated, the following less important accelerators are placed in abeyance: EXTRAP = selects Pople extrapolation of the Fock matrix.


Input Description DAMP SHIFT RSTRCT DEM = = = =

$SCF

2-54

selects Davidson damping of the Fock matrix. selects level shifting of the Fock matrix. selects restriction of orbital interchanges. selects direct energy minimization, which is implemented only for RHF. (default=.FALSE.) EXTRAP T T DAMP F F SHIFT RSTRCT F F F F DIIS F/T F SOSCF T/F F

defaults for ab initio: semiempirical:

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. CUHF = flag requesting Constrained UHF, which causes the occupied beta orbitals of a UHF run to lie entirely within the occupied alpha orbital space. This produces results identical to high spin ROHF! Obviously, this keyword pertains only when using SCFTYP=UHF. The default is .FALSE., meaning a spin-contaminated ordinary UHF solution is sought. Applicable to UHF or UDFT energy and gradients, or to UMP2 energy calculations.

---- These parameters fine tune the various convergers. CONV = SCF density convergence criteria. Convergence is reached when the density change between two consecutive SCF cycles is less than this in absolute value. One more cycle will be executed after reaching convergence. Less accuracy in CONV gives questionable gradients. The default is 1.0d-05, except runs involving CI, MP2, CC, or TDDFT use 1.0d-06 to obtain more crisply converged virtual orbitals.

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)


Input Description

$SCF

2-55

MAXDII = Maximum size of the DIIS linear equations, so that at most MAXDII-1 Fock matrices are used in the interpolation. (default=10) SWDIIS = density matrix convergence at which to switch from DIIS to SOSCF. A value of zero means to keep using DIIS at all geometries, which is the default. However, it may be useful to have DIIS work only at the first geometry, in the initial iterations, for example transition metal ECP runs which has a less good Huckel guess, and then use SOSCF for the final SCF iterations at the first geometry, and ever afterwards. A suggested usage might be DIIS=.TRUE. ETHRSH=2.0 SWDIIS=0.005. This option is not programmed for GVB. LOCOPT = When set to .TRUE., SCF options are locked and do not change during the following: SOSCF and DIIS switch (SWDIIS), DFT grid switch, RHF -> DFT switch (SWOFF). If .FALSE., any of these switches resets some SCF options, such as SHIFT or DAMP. (Default: .FALSE.) RESET = In UHF, reset SOSCF or DIIS if energy rises (Default: .TRUE. for UHF, .FALSE. otherwise).

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 the section on convergence control), but it should not be astronomical either. ** For see ** ******************* more info on the convergence methods, the 'Further Information' section. *******************


Input Description

$SCF

2-56

---- orbital modification options ---The four options UHFNOS, VVOS, MVOQ, and ACAVO are mutually exclusive. The latter 3 require RUNTYP=ENERGY, and should not be used with any correlation treatment. UHFNOS = flag controlling generation of the natural orbitals of a UHF function. (default=.FALSE.) VVOS = flag controlling generation of Valence Virtual Orbitals. See J.Chem.Phys. 120, 2629-2637(2004). The default is .FALSE.

VVOs are a quantitative realization of the concept of the "lowest unoccupied molecular orbital". The implementation allows any elements H-Xe, for RHF, ROHF, and GVB wavefunctions, as well as DFT runs (see also VVOS in $MCSCF). Core potentials may not be used. VVOS should be much better MCSCF starting orbitals than either MVOQ or ACAVO type virtuals. MVOQ =0 =n Skip MVO generation (default) 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. = -1 The cation used will have each valence orbital half filled, to produce MVOs with valence-like character in all regions of the molecule. Implemented for RHF and ROHF only. Flag to request Approximate Correlation-Adapted Virtual Orbitals. Implemented for RHF, ROHF, and GVB (w/o direct SCF). Default is .FALSE.

ACAVO

=

PACAVO =

Parameters used to define the ACAVO generating operator, which is defined as a*T + b*Vne + c*Jcore + d*Jval + e*Kcore + f*Kval The default, PACAVO(1)=0,0,0,0,0,-1, maximizes the exchange interaction with valence MOs (see for example J.L.Whitten, J.Chem.Phys. 56, 5458-5466(1972). The K-orbitals of D.Feller, E.R.Davidson J.Chem.Phys. 74, 3977-3979(1981) are 0.06*F-K(valence), which is PACAVO(1)= 0.06,0.06,0.12,0.12,-0.06,-1.06.


Input Description

$SCF

2-57

Of course, canonical virtuals are PACAVO(1)=1,1,2,2,-1,-1. *** UHFCHK = a flag to perform a RHF --> UHF wavefunction stability test (relaxation of RHF orbitals to unequal alpha and beta orbitals). This option is implemented only for RHF, and involves testing only pairs of orbitals at a single time: this is not foolproof! default is .FALSE. = if UHFCHK is selected, the number of orbitals at The top of the occupied orbital space to be checked for instabilities. Default= -1, meaning check HOMO and HOMO-1. = if UHFCHK is selected, the number of orbitals at the bottom of the virtual space checked. Basis sets with diffuse functions should check many more orbitals, to get past the diffuse MOs. Default=2, meaning check LUMO, LUMO+1, LUMO+2. *** MOM = flag enabling the Molecular Overlap Method (MOM). Currently works only with SCFTYP=UHF. MOM computes excitations by selecting orbitals at each iteration to resemble earlier iterations. See J.Phys.Chem. A 112, 13164-13171(2008). (default=.FALSE.) note: MOM typically requires an initial MOREAD set of ground-state MO's, reordered by NORDER=1 in conjunction with the IORDER/JORDER arrays. The MOM flag is an algorithm very similar to that enabled by the RSTRCT flag. Note: RSTRCT works with all SCF types. KPROJ determines the flavor of the MOM projections p_j obtained from the overlap matrix O(i,j). KPROJ pertains only if the MOM flag is enabled. = 0 p_j = |sum_i O(i,j)| as given in the original MOM article = 1 p_j = sum_i O(i,j) (default) = 2 p_j = sum_i O(i,j)*O(i,j), as implemented in Q-Chem

NHOMO

NLUMO


Input Description

$SCF ----- GVB wavefunction input -----

2-58

The next parameters define the GVB wavefunction. See also MULT in the $CONTRL input. 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! (default=0). The number of sets of open shells in the function. Maximum of 10. (default=0) An array giving the degeneracy of each open shell set. Give NSETO values. (default=0,0,0,...). The number of geminal pairs in the -GVBfunction. Maximum of 12. The default corresponds to open shell SCF (default=0).

NSETO NO

= =

NPAIR

=

CICOEF =

An array of ordered pairs of CI coefficients for the -GVB- pairs. (default = 0.90,-0.20,0.90,-0.20,...) 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, CICOEF will be. This parameter is useful in restarting a GVB run, with the current CI coefficients. COUPLE = A switch controlling the input of F, ALPHA, and BETA. (Default=.FALSE.) Input for F, ALPHA, BETA will be ignored unless you select this variable as .TRUE. F ALPHA BETA = = = An vector of fractional shell occupations. An array of A coupling coefficients, given in lower triangular order. 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 occurring


Input Description

$SCF

2-59

cases are internally stored. See "Further Information" for other cases, including degenerate open shells. Note: ALPHA and BETA can be given for -ROHF- orbital canonicalization control, see "Further Information". ----- miscellaneous options ----NPUNCH = = = = NPREO 0 1 2 option for output to the PUNCH file do not punch out the final orbitals punch out the occupied orbitals punch out occupied and virtual orbitals The default is NPUNCH = 2.

= energy and orbital printing options, applying after other output options, for example NPRINT=-5 for no orbital output overrules this keyword. Orbitals from NPREO(1) to NPREO(2) and orbital energies from NPREO(3) to NPREO(4) are printed. Positive values indicate plain ordinal numbers. Non-positive values are relative to HOMO. For NPREO(1) and (3), 0 is HOMO, -1 is HOMO+1 etc. For NPREO(2) and (4), 0 is HOMO, -1 is HOMO+1 etc. Numbers exceeding the total orbital count are automatically adjusted to the maximum value. Orbitals printed by NPREO(1) and NPREO(2) will always have the orbital energy labels attached, NPREO(3) to NPREO(4) define separate print-out of the orbital energies. HOMO here means the highest occupied orbital, assuming a singlet RHF orbital occupation, that is to say NE/2, no matter what SCFTYP is. To print only the HOMO and LUMO LCAO coefficients. and all orbital energies, enter: NPREO(1)=0,-1,1,9999 Default: 1,9999,2,1 (meaning print all orbitals, but no separate list of orbital energies). ----- 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:


Input Description

$SCF

2-60

SCALF

=

initial exponent scale factor when VTSCAL is in use, useful when restarting. The default is 1.0. maximum number of iterations (at a single geometry) to satisfy the energy virial theorem. The default is 20. 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.

MAXVT

=

VTCONV =

For more information on this option, which is most useful during a geometry search, see M.Lehd and F.Jensen, J.Comput.Chem. 12, 1089-1096(1991). ** For see ** ***************** more discussion of GVB/ROHF input the 'further information' section *****************

==========================================================


Input Description

$SCFMI

2-61

==========================================================

$SCFMI

group

(optional, relevant if SCFTYP=RHF)

The Self Consistent Field for Molecular Interactions (SCF-MI) method is a modification of the usual 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 ten monomers, treated at the RHF level. The energy, gradient, and therefore semi-numerical hessian are available. The SCF step may be run in direct SCF mode, and parallel calculation is also enabled. The calculation must use Cartesian Gaussian AOs only, not spherical harmonics. The SCF-MI driver differs from normal RHF calculations, so not all converger methods are available. Finally, this option is not compatible with electron correlation treatments (DFT, MP2, CI, or CC). The first 3 parameters must be given. All atoms of a fragment must appear consecutively in $DATA. NFRAGS = number of distinct fragments present. Both the supermolecule and its constituent monomers must be well described as closed shells by RHF wavefunctions. = an array containing the number of doubly MOs for each fragment. = an array containing the number of atomic basis functions located on each fragment. = maximum number of SCF-MI cycles, overriding the usual MAXIT value. (default is 50). = SCF-MI density convergence criteria. (default is 1.0d-10)

NF occupied MF ITER DTOL


Input Description ALPHA DIISON MXDIIS and

$SCFMI

2-62

= possible level shift parameter. (default is 0.0, meaning shifting is not used) = a flag to active the DIIS convergence. (default is .TRUE.) = the maximum number of previous effective Fock

overlap matrices to be used in DIIS (default=10) DIISTL = the density change value at which DIIS starts. (default=0.01)

A Huckel guess is localized by the Boys procedure onto each fragment to provide starting orbitals for each: ITLOC CNVLOC IOPT = maximum number of iteration in the localization step (Default is 50) = convergence parameter for the localization. (default is .01). = 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.

========================================================== "Modification of Roothan Equations to exclude BSSE from Molecular Interaction calculations" E. Gianinetti, M. Raimondi, E. Tornaghi Int. J. Quantum Chem. 60, 157-166 (1996) "Implementation of Gradient optimization algorithms and Force Constant computations in BSSE-free direct and conventional SCF approaches" A. Famulari, E. Gianinetti, M. Raimondi, M. Sironi Int. J. Quantum Chem. 69, 151-158 (1997)


Input Description

$DFT

2-63

==========================================================

$DFT

group

(relevant if DFTTYP is chosen) (relevant if SCFTYP=RHF,UHF,ROHF)

Note that if DFTTYP=NONE, an ab initio calculation will be performed, rather than density functional theory. This group permits the use of various one electron (usually empirical) operators instead of the true many electron Hamiltonian. Two programs are provided, METHOD= GRID or GRIDFREE. The programs have different functionals available, and so the keyword DFTTYP (which is entered in $CONTRL) and other associated inputs are documented separately below. Every functional that has the same name in both lists is an identical functional, but each METHOD has a few functionals that are missing in the other. The grid free implementation is based on the use of the resolution of the identity to simplify integrals so that they may be analytically evaluated, without using grid quadratures. The grid free DFT computations in their present form have various numerical errors, primarily in the gradient vectors. Please do not use the grid-free DFT program without reading the discussion in the 'Further References' section regarding the gradient accuracy. The grid based DFT uses a typical grid quadrature to compute integrals over the rather complicated functionals, using two possible angular grid types. Achieving a self-consistent field with DFT is rather more difficult than for normal HF, so DIIS is the default converger. Both DFT programs will run in parallel. See the two lists below for possible functionals in the two programs. See also the $TDDFT input group for excited states. METHOD = selects grid based DFT or grid free DFT. = GRID Grid based DFT (default) = GRIDFREE Grid free DFT


Input Description

$DFT

2-64

DFTTYP is given in $CONTRL, not here in $DFT! Possible values for the grid-based program are listed first, ----- options for METHOD=GRID ----DFTTYP = NONE means ab initio computation (default)

Many choices are given below, perhaps the most sensible are local DFT: SVWN pure DFT GGA: BLYP, PW91, B97-D, PBE/PBEsol hybrid DFT GGA: B3LYP, X3LYP, PBE0 pure DFT meta-GGA: revTPSS hybrid DFT meta-GGA: TPSSh, M06 but of course, everyone has their own favorite! pure exchange functionals: = SLATER Slater exchange = BECKE Becke 1988 exchange = GILL Gill 1996 exchange = OPTX Handy-Cohen exchange = PW91X Perdew-Wang 1991 exchange = PBEX Perdew-Burke-Ernzerhof exchange These will be used with no correlation functional at all. pure correlation functionals: = VWN Vosko-Wilk-Nusair correlation, using their electron gas formula 5 (aka VWN5) = VWN3 Vosko-Wilk-Nusair correlation, using their electron gas formula 3 = VWN1RPA Vosko-Wilke-Nusair correlation, using their e- gas formula 1, with RPA params. = PZ81 Perdew-Zener 1981 correlation = P86 Perdew 1986 correlation = LYP Lee-Yang-Parr correlation = PW91C Perdew-Wang 1991 correlation = PBEC Perdew-Burke-Ernzerhof correlation = OP One-parameter Progressive correlation These will be used with 100% HF exchange, if chosen. combinations (partial list): SLATER exchange + VWN5 correlation

= SVWN


Input Description

$DFT

2-65

Called LDA/LSDA in physics for RHF/UHF. = SVWN1RPA Slater exchange + VWN1RPA correlation = BLYP BECKE exchange + LYP correlation = BOP BECKE exchange + OP correlation = BP86 BECKE exchange + P86 correlation = GVWN GILL exchange + VWN5 correlation = GPW91 GILL exchange + PW91 correlation = PBEVWN PBE exchange + VWN5 correlation = PBEOP PBE exchange + OP correlation = OLYP OPTX exchange + LYP correlation = PW91 means PW91 exchange + PW91 correlation = PBE means PBE exchange + PBE correlation There's a nearly infinite set of pairings (well, 6*9), so we show only enough to give you the idea. In other words, pairs are formed by abbreviating the exchange functionals SLATER=S, BECKE=B, GILL=G, OPTX=O, PW91X=PW91, PBEX=PBE and matching them with any correlation functional, of which only two are abbreviated when used in combinations, PW91C==>PW91, PBEC==>PBE The pairings shown above only scratch the surface, but clearly, many possibilities, such as PW91PBE, are nonsense! pure DFT GGA functionals: = EDF1 empirical density functional #1, which is a modified BLYP from Adamson/Gill/Pople. = PW91 Perdew/Wang 1991 = PBE Perdew/Burke/Ernzerhof 1996 = revPBE PBE as revised by Zhang/Yang = RPBE PBE as revised by Hammer/Hansen/Norskov = PBEsol PBE as revised by Perdew et al for solids = HCTH93 Hamprecht/Cohen/Tozer/Handy's 1998 mod to B97, omitting HF exchange, fitting to 93 atoms and molecules = HCTH120 later fit to 120 systems = HCTH147 later fit to 147 systems = HCTH407 later fit to 407 systems (best) = SOGGA PBE revised by Zhao/Truhlar for solids = MOHLYP metal optimized OPTX, half LYP = B97-D Grimme's modified B97, with dispersion correction (this forces DC=.TRUE.) = SOGGA11 optimized with broad applicability for chemistry, by Peverati/Zhao/Truhlar


Input Description

$DFT

2-66

hybrid GGA functionals: = BHHLYP HF and BECKE exchange + LYP correlation = B3PW91 Becke's 3 parameter exchange hybrid, with PW91 correlation functional = B3LYP this is a hybrid method combining five functionals: Becke + Slater + HF exchange (B3), with LYP + VWN5 correlation. B3LYPV5 is a synonym for B3LYP. = B3LYPV1R use VWN1RPA in place of VWN5, matches the e- gas formula chosen by some programs. = B3LYPV3 use VWN3 in place of B3LYP's VWN5 = B3P86 B3-type exchange, P86 correlation, using VWN3 as the LDA part of the correlation. B3P86V3 is a synonym for B3P86. = B3P86V1R use VWN1RPA in place of VWN3 = B3P86V5 use VWN5 in place of VWN3 = B97 Becke's 1997 hybrid functional = B97-1 Hamprecht/Cohen/Tozer/Handy's 1998 reparameterization of B97 = B97-2 Wilson/Bradley/Tozer's 2001 mod to B97 = B97-3 Keal/Tozer's 2005 mod to B97 = B97-K Boese/Martin's 2004 mod for kinetics = B98 Schmider/Becke's 1998 mode to B97, using their best "2c" parameters. = PBE0 a hybrid made from PBE = X3LYP HF+Slater+Becke88+PW91 exchange, and LYP+VWN1RPA correlation. = SOGGA11X a hybrid based on SOGGA11, with 40.15% HF exchange. Each includes some Hartree-Fock exchange, and also may use a linear combination of many DFT parts. range separated functionals: These are also known as "long-range corrected functionals". LC-BVWN, LB-BOP, LC-BLYP, or LC-BPBE are available by selecting BVWN, BOP, BLYP, or BPBE and also setting the flag LC=.TRUE. (see LC and also MU below). Others are selected by their specific name, without using LC: = CAMB3LYP coulomb attenuated B3LYP = wB97 omega separated form of B97 = wB97X wB97 with short-range HF exchange = wB97X-D dispersion corrected wB97X M11 is also range-separated, but is listed below with the other meta-GGAs.


Input Description

$DFT

2-67

"double hybrid" GGA: = B2PLYP mixes BLYP, HF exchange, and MP2! See related inputs CHF and CMP2 below. "double hybrid" and "range separated": = wB97X-2 intended for use with GBASIS=CCT,CCQ,CC5 = wB97X-2L intended for use with GBASIS=N311 NGAUSS=6 NDFUNC=3 NFFUNC=1 NPFUNC=3 DIFFSP=.T. DIFFS=.T. Note: there are no analytic gradients for "double hybrids". Note: the B2PLYP family uses the conventional MP2 energy and may be used for closed shell or spin-unrestricted open shell cases. The wB97X-2 family uses the SCS-MP2 energy, and thus is limited to closed shell cases at present. meta-GGA functionals: These are not hybridized with HF exchange, unless that is explicitly stated below. = VS98 Voorhis/Scuseria, 1998 = PKZB Perdew/Kurth/Zupan/Blaha, 1999 = tHCTH Boese/Handy's 2002 metaGGA akin to HCTH = tHCTHhyb tHCTH's hybrid with 15% HF exchange = BMK Boese/Martin's 2004 parameterization of tHCTHhyb for kinetics = TPSS Tao/Perdew/Staroverov/Scuseria, 2003 = TPSSh TPSS hybrid with 10% HF exchange = TPSSm TPSS with modified parameter, 2007 = revTPSS revised TPSS, 2009 = dlDF a reparameterized M05-2X, reproducing interaction energies which have had all dispersion removed. This MUST be used with a special -D correction to recover dispersion. See 'Further References'. = M05 Minnesota exchange-correlation, 2005 a hybrid with 28% HF exchange. = M05-2X M05, with doubled HF exchange, to 56% = M06 Minnesota exchange-correlation, 2006 a hybrid with 27% HF exchange. = M06-L M06, with 0% HF exchange (L=local) = M06-2X M06, with doubled HF exchange, to 54% = M06-HF M06 correlation, using 100% HF exchange = M08-HX M08 with 'high HF exchange' = M08-SO M08 with parameters that enforce the


Input Description

$DFT

2-68

correct second order gradient expansion. = M11 M11 range-separated hybrid = M11-L M11 local (0% HF exchange) with dual-range exchange When the M06 family was created, Truhlar recommended M06 for the general situation, but see his "concluding remarks" in the M06 reference about which functional is best for what kind of test data set. The most recent M11 family is probably a better choice, and two functionals fit all the needs of the older M05/M06/M08 families. An extensive bibliography for all functionals can be found in the 'Further References' section of this manual. Note that only a subset of these functionals can be used for TD-DFT energy or gradients. These subsets are listed in the $TDDFT input group. * * * dispersion corrections * * * Many exchange-correlation functionals fail to compute intra- and inter-molecular dispersion interactions accurately. Two possible correction schemes are provided below. The first uses empirically chosen C6 and C8 coefficients, while the latter obtains these from the molecular DFT densities. At most, only one of the LRDFLG or DC options below may be chosen. DC = a flag to turn on Grimme's empirical dispersion correction, involving scaled R**(-6) terms. N.B. This empiricism may also be added to plain Hartree-Fock, by choosing DFTTYP=NONE with DC=.T. Three different versions exist, see IDCVER. (default=.FALSE., except if DFTTYP=B97-D, wB97X-D)

IDCVER = 1 means 1st 2004 implementation. = 2 means 2nd 2006 implementation DFT-D2, default for B97-D, wB97X-D. = 3 means 3rd 2010 implementation DFT-D3, default for all others. Setting IDCVER will force DC=.TRUE. DCCHG = a flag to use Chai-Head-Gordon damping function instead of Grimme's 2006 function. Pertinent only


Input Description

$DFT

2-69

for the DFT-D2 method. Forces DC=.TRUE. (default=.FALSE. except for wB97X-D) DCABC = a flag to turn on the computation of the E(3) nonadditive energy term. Pertinent only for DFT-D3, it forces DC=.TRUE. (default=.FALSE.)

The following parameters govern Grimme's semiempirical dispersion term. They are basis set and functional dependent, so they exist for only a few DFTTYP. Default values are automatically selected and printed out in the output file for many common density functionals. The following keywords are for entering non-standard values. For DFT-D2 values, see also: R.Peverati and K.K.Baldridge J.Chem.Theory Comput. 4, 2030-2048 (2008). For DFT-D3 values, and a detailed explanation of each parameter, see: S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J.Chem.Phys. 132, 154104/1-19(2010) DCALP = alpha parameter in the DFT-D damping function (same as alpha6 in Grimme's DFT-D3 notation). Note also that alpha8 and alpha10 in DFT-D3 have constrained values of: alpha8 = alpha6 + 2, alpha10 = alpha8 + 2. Default=14.0 for DFT-D3 =20.0 for DFT-D2 =23.0 for DFT-D1 =6.00 for DCCHG=.TRUE. = sR exponential parameter to scale the van der Waals radii (same as sR,6 in Grimme's DFT-D3 notation). Note also that sR,8 in DFT-D3 have a fixed value of 1.0. Optimized values are automatically selected for some of the more common functionals, otherwise, the default is 1.00 for DFT-D3, 1.10 for DFT-D2, and 1.22 for DFT-D1. = s6 linear parameter for scaling the C6 term. Optimized values are automatically selected for some of the more common functionals, otherwise, the default is 1.00.

DCSR

DCS6


Input Description DCS8

$DFT

2-70

= s8 linear parameter for scaling the C8 term of DFT-D3. Pertinent only for DFT-D3. Optimized values are automatically selected for some of the more common functionals, otherwise, the default is 1.00.

The old keywords DCPAR and DCEXP were replaced by DCS6 and DCSR in 2010. Similarly, DCOLD has morphed into IDCVER. --The Local Response Dispersion (LRD) correction includes atomic pair-wise -C6/R**6, -C8/R**8, and -C10/R**10 terms, whose coefficients are computed from the molecular system's electron density and its nuclear gradient. The nuclear gradient assumes the dispersion coefficients do not vary with geometry, which causes only a very small error in the gradient. Optionally, 3 and 4 center terms may be added, at the 1/R**6 level; in this case, nuclear gradients may not be computed at all. Since the three numerical parameters are presently known only for the long-range exchange corrected BOP functional, calculations may specify simply DFTTYP=LCBOPLRD. The "LCBOPLRD" functional will automatically select the following: DFTTYP=BOP LC=.TRUE. MU=0.47 LRDFLG=.TRUE. LAMBDA=0.232 KAPPA=0.600 RZERO=3.22 leaving only the choice for MLTINT up to you. References for LRD are T.Sato, H.Nakai J.Chem.Phys. 131, 224104/1-12(2009) T.Sato, H.Nakai J.Chem.Phys. 133, 194101/1-9(2010) LRDFLG = flag choosing the Local Response Dispersion (LRD) C6, C8, and C10 corrections. Default=.FALSE. MLTINT = flag to add the 3 and 4 center 6th order terms, the default=.FALSE. Note that nuclear gradients are not available if these multi-center terms are requested. Three numerical parameters may be input. The defaults shown are optimized for the BOP functional with the LC correction for long-range exchange.


Input Description

$DFT

2-71

LAMBDA = parameter adjusting the density gradient correction for the atomic and atomic pair polarizabilities. (default=0.232) KAPPA = parameter in the damping function (default=0.600) RZERO = parameter in the damping function (default=3.22) It may be interesting to see a breakdown of the total dispersion correction, using these keywords: PRPOL = print out atomic effective polarizabilities (default=.FALSE.) PRCOEF = N (default N=0) print out dispersion coefficient to N-th order. PRPAIR = print out atomic pair dispersion energies (default=.FALSE.) * * * range separation * * * LC = flag to turn on the long range correction (LC), which smoothly replaces the DFT exchange by the HF exchange at long inter-electron distances. (default=.FALSE.) This option can be used only with the Becke exchange functional (Becke) and a few correlation functionals: DFTTYP=BVWN, BOP, BLYP, BPBE only. For example, B3LYP has a fixed admixture of HF exchange, so it cannot work with the LC option. See H.Iikura, T.Tsuneda, T.Yanai, and K.Hirao, J.Chem.Phys. 115, 3540 (2001). = A parameter for the long range correction scheme. Increasing MU increases the HF exchange used, very small MU produces the DFT limit. (default=0.33)

MU

Other range-separated options exist, invoked by naming the functional, such as DFTTYP=CAMB3LYP (see the DFTTYP keyword for a full list). * * * B2x-PLYP double hybrid functionals * * * B2xPLYP Double Hybrid functionals have the general formula:


Input Description

$DFT

2-72

Exc = (1-cHF) * ExGGA + cHF * ExHF + (1-cMP2) * EcGGA + cMP2 * E(2) The next keywords allow the choice of cHF and cMP2. Both values must be between 0 and 1 (0-100%). CHF CMP2 = amount of HF exchange. (default=0.53) = amount of MP2. (default=0.27)

Some other common double hybrid functionals are available simply by choosing DFTTYP=B2PLYP, and changing the CHF and CMP2 parameters. Popular parametrizations are: CHF CMP2 -----------------------------------------B2-PLYP (default) | 0.53 | 0.27 | -----------------------------------------B2K-PLYP | 0.72 | 0.42 | -----------------------------------------B2T-PLYP | 0.60 | 0.31 | -----------------------------------------B2GP-PLYP | 0.65 | 0.36 | -----------------------------------------* * * Grid Input * * * Only one of the three grid types may be chosen for the run. The default (if no selection is made) is the Lebedev grid. In order to duplicate results obtained prior to April 2008, select the polar coordinate grid NRAD=96 NTHE=12 NPHI=24. Energies can be compared if and only if the identical grid type and density is used, analogous to needing to compare with the identical basis set expansions. See REFS.DOC for more information on grids. See similar inputs in $TDDFT. Lebedev grid: NRAD NLEB = number of radial points in the Euler-MacLaurin quadrature. (default=96) = number of angular points in the Lebedev grids. (default=302). Possible values are 86, 110, 146, 170, 194, 302, 350, 434, 590, 770, 974, 1202, 1454, 1730, 2030...


Input Description

$DFT

2-73

Meta-GGA functionals require a tighter grid to achieve the same accuracy. For this reason a tighter default grid of NRAD=99 and NLEB=590 is chosen by default with all meta-GGA functionals. The default for NLEB means that nuclear gradients will be accurate to about the default OPTTOL=0.00010 (see $STATPT), 590 approaches OPTTOL=0.00001, and 1202 is "army grade". The next two specify radial/angular in a single keyword: SG1 = a flag to select the "standard grid 1", which has 24 radial points, and various pruned Lebedev grids, from 194 down to 6. (default=.FALSE. This grid is very fast, but produces gradients whose accuracy reaches only OPTTOL=0.00050. This grid should be VERY USEFUL for the early steps of a geometry optimization. = two unpublished grids due to Curtis Janssen, implemented here differently than in MPQC: = 1 uses 95 radial points for all atoms, and prunes from a Lebedev grid whose largest size is 434, thus using about 15,000 grid points/atom. = 2 uses 155 radial points for all atoms, and prunes from a Lebedev grid whose largest size is 974, thus using about 71,000 grid points/atom. This is a very accurate grid, e.g. "army grade". The information for pruning exists only for H-Ar, so heavier elements will use the large radial/ Lebedev grid without any pruning.

JANS

polar coordinate grid: NRAD NTHE NPHI = number of radial points in the Euler-MacLaurin quadrature. (96 is reasonable) = number of angle theta grids in Gauss-Legendre quadrature (polar coordinates). (12 is reasonable) = number of angle phi grids in Gauss-Legendre quadrature. NPHI should be double NTHE so points are spherically distributed. (24 is reasonable)


Input Description

$DFT

2-74

The number of angular points will be NTHE*NPHI. The values shown give a gradient accuracy near the default OPTTOL of 0.00010, while NTHE=24 NPHI=48 approaches OPTTOL=0.00001, and "army grade" is NTHE=36 NPHI=72. * * * Grid Switching * * * At the first geometry of the run, pure HF iterations will be performed, since convergence of DFT is greatly improved by starting with the HF density matrix. After DFT engages, most runs (at all geometries, except for PCM or numerical Hessians) will use a coarser grid during the early DFT iterations, before reaching some initial convergence. After that, the full grid will be used. Together, these switchings can save considerable CPU time. SWOFF = turn off DFT, to perform pure SCF iterations, until the density matrix convergence falls below this threshold. This option is independent of SWITCH and can be used with or without it. It is reasonable to pick SWOFF > SWITCH > CONV in $SCF. SWOFF pertains only to the first geometry that the run computes, and is automatically disabled if you choose GUESS=MOREAD to provide initial orbitals. The default is 5.0E-3.

SWITCH = when the change in the density matrix between iterations falls below this threshhold, switch to the desired full grid (default=3.0E-4) This keyword is ignored if the SG1 grid is used. NRAD0 NLEB0 NTHE0 NPHI0 = same as NRAD, but defines initial coarse grid. default = smaller of 24 and NRAD/4 = same as NLEB, but defines initial coarse grid. default = 110 = same as NTHE, but defines initial coarse grid. default = smaller of 8, NTHE/3 = same as NPHI, but defines initial coarse grid. default = smaller of 16, NPHI/3


Input Description technical parameters:

$DFT

2-75

THRESH = threshold for ignoring small contributions to the Fock matrix. The default is designed to produce no significant energy loss, even when the grid is as good as "army grade". If for some reason you want to turn all threshhold tests off, of course requiring more CPU, enter 1.0e-15. default: 1.0e-4/Natoms/NRAD/NTHE/NPHI GTHRE = threshold applied to gradients, similar to THRESH. < 1 assign this value to all thresholds = 1 use the default thresholds (default). > 1 divide default thresholds by this value. If you wish to increase accuracy, set GTHRE=10. The default introduces an error of roughly 1e-7 (a.u./bohr) in the gradient.


Input Description

$DFT

2-76

The keyword $DFTTYP is given in $CONTRL, and may have these values if the grid-free program is chosen: ----- options for METHOD=GRIDFREE ----DFTTYP = NONE means ab initio computation (default) exchange functionals: = XALPHA X-Alpha exchange (alpha=0.7) = SLATER Slater exchange (alpha=2/3) = 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 + VWN5 correlation = BLYP Becke exchange + LYP correlation = BPWLOC Becke exchange + Perdew/Wang correlation = B3LYP hybrid HF/Becke/LYP using VWN formula 5 = CAMB CAMA exchange + Cambridge correlation = XVWN Xalpha exchange + VWN5 correlation = XPWLOC Xalpha exchange + Perdew/Wang correlation = SVWN Slater exchange + VWN5 correlation = SPWLOC Slater exchange + PWLOC correlation = WIGNER Wigner exchange + correlation = WS Wigner scaled exchange + correlation = WIGEXP Wigner exponential exchange + correlation uses no auxiliary basis set for resolution of the identity, limiting accuracy. 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)

AUXFUN = AUX0 = AUX3

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.) ==========================================================


Input Description

$TDDFT

2-77

==========================================================

$TDDFT

group (relevant if TDDFT chosen in $CONTRL)

This group generates molecular excitation energies by time-dependent density functional theory computations (or time-dependent Hartree-Fock, also known as the Random Phase Approximation). The functional used for the excited states is necessarily the same one that is used for the reference state, specified by DFTTYP in $CONTRL. For conventional TD-DFT (TDDFT=EXCITE in $CONTRL), the orbitals are optimized for RHF or UHF type reference states. Analytic nuclear gradients are available for singlet excited states, while the energy of excited states of other multiplicities can be computed. Two-photon absorption cross-sections may be predicted for singlet excited states. Ground state hyperpolarizabilities may be computed with the TDDFT module. For spin-flip TD-DFT (TDDFT=SPNFLP in $CONTRL), the calculation obtains orbitals for a reference state of either UHF or ROHF type, with MULT in $CONTRL determining the Ms quantum number of the reference. The reference state's Ms is set equal to the S value implied by $CONTRL's MULT=2S+1. The SF-TD-DFT then uses only determinants with Ms=S-1, due to the flip of one alpha spin into a beta spin. This means that target states (which are spin contaminated) will have multiplicities around the range S-1 to S, only. It is quite possible for some of the target states to have a lower energy than the reference!!! Nuclear gradients and properties are available. See just below for "limitations" below regarding the two different TD-DFT types. TD-DFT is a single excitation theory. All of the caveats listed in the $CIS input group about states with double excitation character, need for Rydberg basis sets, greatly different topology of excited state surfaces, and so on apply here as well. Please read the introduction to the $CIS input group! If you use very large or very small Gaussian exponents, you may need to increase the number of


Input Description

$TDDFT

2-78

radial grid points (the program prints advice in such cases). TDHF, TDDFT, and CIS are related in the following way: -- Tamm/Dancoff approximation --> | TDHF CIS DFT | V TDDFT TDDFT/TDA Here TDHF means absorption of photons, to produce excited states (TDHF is called RPA in the physics community). This meaning of TDHF should not be confused with the photon scattering processes computed by RUNTYP=TDHF or TDHFX, which generate polarizabilities. Note, in particular, that CITYP=CIS is equivalent to using TDDFT=EXCITE DFTTYP=NONE TAMMD=.TRUE., provided the former is run with no frozen cores. Solvent effects for CIS calculations are therefore available via the TDDFT codes. Excited state properties are calculated using the TDDFT excited state electronic density only during gradient runs, or by setting TDPRP below. The TD-DFT codes excite all electrons, that is, there is no frozen core concept. Please see the 4th chapter of this manual for more information on both types of TD-DFT. "limitations" for TDDFT=EXCITE: Permissible values for DFTTYP are shown below. These include "NONE" which uses TDHF (i.e. the Random Phase Approximation), noting that extra states may need to be solved for in order to be sure of getting the first few states correctly. If nuclear gradients are needed, you may choose any of the following functionals: NONE SVWN, SOP, SLYP, OLYP, BVWN, BOP, BLYP (and their LC=.TRUE. versions) B3LYP, CAMB3LYP, B3LYP1, PBE, PBE0 For evaluation of just the excitation energies, you may use many more functionals, notably including the metaGGAs in the last three lines: NONE SVWN, SVWN1, SPZ81, SP86, SOP, SLYP, BVWN, BVWN1, BPZ81, BP86, BOP, BLYP, OLYP,


Input Description

$TDDFT

2-79

B3LYP, CAMB3LYP, B3LYP1, B3PW91, X3LYP, PW91, PBE, PBE0 VS98, PKZB, M05, M05-2X, M06, M06-HF, M06-L, M06-2X, M08-HX, M08-SO TPSS, TPSSm, TPSSh, and revTPSS The LC flag in $DFT automatically carries over to TDDFT runs. The LC option may be used with the "B" functionals, and (like the similar range-separated CAMB3LYP) is useful in obtaining better descriptions for charge-transfer excitations or Rydberg excitation energies than are the conventional exchange correlation functionals (whether pure or hybrid). The LC flag is also available for excited state gradient computation. Limits specific to the references for TDDFT=EXCITE are: For SCFTYP=RHF, excitation energies can be found for singlet or triplet coupled excited states. For singlet excited states only, analytic gradients and properties can be found, for either full TD-DFT or in the Tamm/Dancoff approximation. For RHF references, solvent effects can be included by EFP1 or PCM (or both together), for both TD-DFT excitation energies and their nuclear gradients. For SCFTYP=UHF, excited states with the same spin projection as the ground state are found. MULT in $CONTRL governs the number of alpha and beta electrons, hence Ms=(MULT-1)/2 is the only good quantum number for either the ground or excited states. Since U-TDDFT is a single excitation theory, excited states with values near Ms and near Ms+1 will appear in the calculation. There are no properties other than the excitation energy, nor gradients, nor solvent effects, at present. "limitations" for TDDFT=SPNFLP: Spin-flip TDDFT is programmed in the "collinear approximation" which means only the HF exchange term carries a large impact on the excitation energies. Pure DFT functionals may be used, but normally hybrids with large HF exchange fractions are used. The LC option for range-separation hybrids cannot be used, which also removes CAMB3LYP. Finally, no meta-GGA may be used. Note that spin-flip TD-DFT in the Tamm/Dancoff approximation using DFTTYP=NONE is equivalent to spin-flip CIS.


Input Description

$TDDFT

2-80

MULT below is ignored, as the Ms of target states is fixed solely by MULT in $CONTRL. The spin-flip code operates only in the Tamm/Dancoff approximation, so TAMMD below is automatically .TRUE. Nuclear gradients and/or excited state properties are available only in the gas phase. Solvation effects are available for both energy and gradient calculations, for EFP1, C-PCM, or both. --------NSTATE = Number of states to be found (excluding the reference state). The default is 1 more state. IROOT = State used for geometry optimization and property evaluation. (default=1) TDDFT=EXCITE counts the reference as 0, and this should be the lowest state. Hence IROOT=1 means the 1st excited state, just as you might guess. TDDFT=SPNFLP labels the reference state as 0, but this might not be the lowest state overall. The meaning of IROOT=1 is the lowest state, omitting the reference state from consideration. Hence IROOT=1 might specify the ground state! = Multiplicity (1 or 3) of the singly excited states. This keyword applies only when the reference is a closed shell. (default is 1) This parameter is ignored when TDDFT=SPNFLP. = a flag to request property computation for the state IROOT. Properties can only be obtained when the nuclear gradient is computable, see gradient restrictions noted in the introduction to this group. Properties require significant extra computer time, compared to the excitation energy alone, so the default is .FALSE. Properties are always evaluated during nuclear gradient runs, when they are a free by-product. = a flag requesting two-photon absorption crosssections. These are computed for each of the NSTATE excited states, after first evaluating

MULT

TDPRP

TPA


Input Description

$TDDFT

2-81

their excitation energies. The TPA calculation is only available for closed shell references, only for singlet excited states (MULT=1), and may not be used with the Tamm/Dancoff approximation. Solvent effects may be treated by EFP. TAMMD is a flag selecting the Tamm/Dancoff approximation be used. This may be used with closed shell excitation energies or gradients, or open shell excitation energies. Default = .FALSE. This parameter is ignored by TDDFT=SPNFLP, which is only coded in the Tamm/Dancoff approximation. is a flag controlling PCM's solvent behavior: .TRUE. splits the dielectric constant into a bulk value (EPS in $PCM) and a fast component (EPSINF), see Cossi and Barone, 2001. The idea is that NONEQ=.t. is appropriate for vertical excitations, and .f. for adiabatic. (the default is .TRUE.) This keyword is ignored by TDDFT=SPNFLP.

NONEQ

* * * ground state polarizability calculation * * * (use TDDFT=HPOL option in $CONTRL) These two frequency dependent polarizability calculations may be requested in the same run (more efficient). These properties are available only for closed shell references, require the default MULT=1 value in this input group, and may not be used with the Tamm/Dancoff approximation. Solvent effects may be treated by EFP. ALPHA BETA PFREQ = requests the polarizability. Default=.FALSE. If BETA is not chosen, give just one PFREQ. = requests the hyperpolarizability. Default=.FALSE. Two values are required for PFREQ. = an array of one or two input frequencies, omega1 and omega2, to yield the polarizability alpha(omega1;omega1) [if BETA=.F.] alpha(omega2;omega2) [if BETA=.T.] alpha(omega3;omega3) [if BETA=.T.] and/or to yield the hyperpolarizability beta(omega3;omega1,omega2).


Input Description

$TDDFT

2-82

The output photon frequency is determined from omega3=-(omega1+omega2). Useful special cases second harmonic generation beta(-2W;W,W), electro-optic Pockels effect beta(-W;W,0), and optical rectification beta(0;W,-W) are among the possibilities. The default is the static polarizability and/or static hyperpolarizability: PFREQ(1)=0.0,0.0 PFREQ is given in atomic units: PFREQ=45.56/lamda, for wavelength lambda in nm. * * * Grid Selection * * * The grid type and point density used in $TDDFT may be chosen independently of the values in $DFT. Excitation energies accurate to 0.01 eV may be obtained with grids that are much sparser than those needed for the ground state, and this is reflected in the defaults. Prior to April 2008, the default grid was NRAD=24 NTHE=8 NPHI=16. NRAD = number of radial grid points in Euler-Maclaurin quadrature, used in calculations of the second or third derivatives of density functionals. (default=48) = number of angular points in the Lebedev grid. (default=110) = number of theta grid points if a polar coordinate grid is used. = number of phi grid points if a polar coordinate grid is used. NPHI should be twice NTHE. = flag selecting "standard grid one". (default=.FALSE.)

NLEB NTHE NPHI SG1

See both $DFT and REFS.DOC for more information on grids. The "army grade" standard for $TDDFT is NRAD=96 combined with either NLEB=302 or NTHE=12/NPHI=24. the remaining parameters are technical in nature:


Input Description

$TDDFT

2-83

CNVTOL = convergence tolerance in the iterative TD-DFT step. (default=1.0E-7) MAXVEC = the maximum number of expansion vectors used by the solver's iterations, per state (default=50). The total size of the expansion space will be NSTATE*MAXVEC. NTRIAL = the number of initial expansion vectors used. (default is the larger of 5 and NSTATE). ==========================================================


Input Description

$DFTB

2-84

==========================================================

$DFTB

group

(relevant for GBASIS=DFTB)

Density-functional tight-binding (DFTB) is turned on by selecting GBASIS=DFTB in $BASIS. $DFTB controls optional parameters for a DFTB calculation. DFTB is formulated in a two-center approximation utilizing implicitly a minimal pseudoatomic orbital basis set with corresponding, pretabulated one- and two-center integrals. Because of this, many properties (for instances, multipoles higher than dipoles) and many options are ignored or not available in the current implementations of DFTB. DFTB also uses an independent SCF driver (SCF in DFTB is also called SCC, see below), so most SCF options are not available for DFTB. Only SCFTYP=RHF and UHF are implemented. SCFTYP=ROHF is available, only when all SPNCST values are zero. DFTB does not explicitly use symmetry (C1 throughout) since integrals are never computed during the calculations. Slater-Koster tables are only defined for spherical functions (5d) so DFTB sets ISPHER=1. Most $GUESS options do not work for DFTB (DFTB does not use initial orbitals in the usual sense). Other than the default (METHOD=HUCKEL, which is ignored), only METHOD=MOREAD works (note that SCC-DFTB can use initial charges on atoms, derived from the orbitals). RUNTYP=OPTIMIZE, HESSIAN and RAMAN are available for full (non-FMO) DFTB, whereas RUNTYP=OPTIMIZE and OPTFMO are available for FMO-DFTB. At present, DFTB may not be combined with any solvation model, and excited states (i.e., from time-dependent DFTB), are also not available. In DFTB calculation, the atom type is determined by its name, not its nuclear charge as elsewhere in GAMESS. The nuclear charge (the second column in $DATA) is used only in population analysis, but not in SCF. DFTB uses a notion of "species", which means an atomic type. The species are numbered according to the order in which atoms appear in $DATA. For instances, in water there are two species, O and H. An atomic type of each species needs MAXANG, which for most but not all atoms is set automatically.


Input Description NDFTB

$DFTB

2-85

order of the Taylor expansion of the total energy around a reference density in the DFTB model. = 1 NCC-DFTB, also called DFTB1. NCC stands for non-charge-consistent, i.e., no explicity charge-charge interaction term is included in the energy calculation. = 2 SCC-DFTB, also called DFTB2. SCC means a self-charge-consistent approach, and SCC implies that SCF iterations are carried out that converge monopolar charges towards self-consistency. = 3 DFTB3, including 3rd order correction using Hubbard derivatives (HUBDER). In order to reproduce the published DFTB3 approach, it is necessary to also specify DAMPXH=.TRUE. to add other terms. Gaus, M. et al. J. Chem. Theory Comput. 2011, 7, 931-948 is referred to as Gaus2011 below. Default: 2. a flag to include the damping function for X-H atomic pair in DFTB3. See also DAMPEX, and eq 21 in Gaus2011. The damping function is used when at least one atom in a pair is "H". "HYDROGEN" and any other name will turn off the damping. Default: .false. an exponent used in the damping function for X-H atomic pairs. The default value is 4.2 (see Table 2 in Gaus2011 for more details). a flag to perform shell-resolved SCC calculation. If set to .false., the code uses the Hubbard value for an s orbital for p and d orbitals, ignoring their Hubbard values defined in SlaterKoster tables. Using .true. enables the use of proper Hubbard values for p and d orbitals, implemented only for DFTB1 and DFTB2. Default: .false.

DAMPXH =

DAMPEX =

SRSCC

=

ITYPMX

Convergence method of SCC calculations. = -1 Use standard GAMESS convergence methods. SOSCF and DIIS are supported, but DEM is not.


Input Description =0

$DFTB Broyden's method. Interpolation is applied for atomic (or shell-resolved when SRSCC=.TRUE.) charges, but not Hamiltonian matrix. (reserved) DIIS for charges. Default: 0.

2-86

=1 =2 ETEMP

= electronic temperature in Kelvin. Non-zero values of ETEMP help SCF convergence of nearly-degenerate systems by smearing occupation numbers around the Fermi-level. Only the Fermi-Dirac distribution function is available as a smearing function. The default value is 0 Kelvin, meaning the smearing function is not used. ETEMP is implemented only for SCFTYP=RHF and when FMO is not used. dispersion model for DFTB. = NONE no Dispersion correction. = UFF UFF-type dispersion correction. Parameters for atomic numbers up to 36 are available internally or can be supplied in DISPPR for any atom. Built-in parameters are taken from Rappe et al. J. Am. Chem. Soc. 1992, 114, 10024. = SK The Slater-Kirkwood type dispersion correction omitting the change polarizability depending on the number of bonds. No default values of DISPPR are available. Some are listed in the manual of the DFTB+ program. an array of parameters used for dispersion correction, listed in sets for each species. For DISP=UFF, DISPPR(1) and DISPPR(2) define the non-bonded distance (Angs.) and energy (kcal/mol) for the first species, respectively, and so on. For DISP=SK, a set for a species has 3 parameters, the polarizability (Angstrom^3), cutoff length (Angstrom), and atomic charge. Default: see DISP. an array of Hubbard derivatives for each species (1 per species) used only for DFTB3 calculations.

DISP

DISPPR

HUBDER


Input Description

$DFTB

2-87

Default values are set only for C, H, N, O, and P using the final row of Table 2 in Gaus2011 (see the paper for other choices). MAXANG array of maximum angular momentum of each species, which determines the number of basis functions. DFTB uses only valence orbitals and electrons! Most elements have proper default values, but for some atomic types (i.e., species) you need to manually define the values. an array of spin constants used in unrestricted (UHF) DFTB calculation. Provide 6 spin constants, W_{ss}, W_{sp}, W_{pp}, W_{sd}, W_{pd}, & W_{dd}, for each species in a continuous array. Constants for some elements can be found in the manual of the DFTB+ program. *** The following options are FMO-DFTB specific (Nishimoto, Y. et al. J. Chem. Theory Comput. 10, 4801-4812, 2014). FMO-DFTB has many limitations and some FMO options are not supported (for instance, AFO, multilayer FMO etc). Only single layer, restricted (RHF) FMO2-DFTB1 and FMO2-DFTB2 are implemented at present. DFTB3, SRSCC, ETEMP etc are not available. FMO-DFTB is available for the energy and approximate gradient (no SCZV terms). MODESD = controls the behavior of ES-DIM (electrostatic dimer) approximation (bit additive). 1 Calculate interfragment repulsive energy for ES dimers (almost never used). 2 Add up all ES-DIM energies. This means that individual ES dimer energies are not calculated, but only their total lump sum, computed with the dynamic load balancing. 4 Lump ES-DIM routine with static load balancing. The bits of 2 or 4 are mutually exclusive. Default: 0 (i.e., individual ES dimer energies). MODGAM = controls the calculation of gamma values (interatomic 1/R-like function) in FMO-DFTB2.

SPNCST


Input Description

$DFTB

2-88

0 Calculate gamma values on the fly. (default) 1 Calculate once and prestore gamma values. This option is faster but takes more memory. ==========================================================


Input Description

$DFTSBK

2-89

==========================================================

$DFTBSK

group

(required if GBASIS=DFTB)

This group is required for all DFTB calculations. It defines the file names for Slater-Koster parameters describing pairwise interactions for each pair of species (see $DFTB for the introduction). The group is free-formatted, and consists of N*N lines, where N is the number of species. Each line has the format: species-i species-j "i-j parameter file name" Note that (species-i species-j) is different from (speciesj species-i), and both must be given! For example, C O "/home/user/DFTB/mio-0-1/C-O.skf" which defines carbon-oxygen parameters to be in that file. Naturally, this file must be readable on the computer where calculations are actually run. The names of species should be the same as in $DATA (for FMO, in $FMOXYZ). In FMO-DFTB calculations, these names must be atomic symbols (H, He, ...), otherwise names may be used (hydrogen, helium, ...) truncated to 8 symbols. The names may include numbers, as in (H.1). Slater-Koster parameters are available from http://www.dftb.org/ (for free) and from some other sources. There are many different parameter sets, for instance, the above example uses the "mio-0-1" set. Follow the requirements on proper citations when using these parameters. ==========================================================


Input Description

$CIS

2-90

==========================================================

$CIS

group

required when CITYP=CIS required when CITYP=SFCIS

The CIS method (singly excited CI) is the simplest way to treat excited states. By Brillouin's Theorem, a single determinant reference such as RHF will have zero matrix elements with singly substituted determinants. The ground state reference therefore has no mixing with the excited states treated with singles only. Reading the references given in Section 4 of this manual will show the CIS method can be thought of as a non-correlated method, rigorously so for the ground state, and effectively so for the various excited states. Some issues making CIS rather less than a black box method are: a) any states characterized by important doubles are simply missing from the calculation. b) excited states commonly possess Rydberg (diffuse) character, so the AO basis used must allow this. c) excited states often have different point group symmetry than the ground state, so the starting geometries for these states must reflect their actual symmetry. d) excited state surfaces frequently cross, and thus root flipping may very well occur. The normal CIS implementation allows the use of only RHF references, but can pick up both singlet and triplet excited states. Nuclear gradients are available, as are properties. The CIS run automatically includes computation of the dipole moments of all states, and all pairwise transition dipoles and oscillator strengths. The spin-flip type of CIS is very similar to spin-flip TDDFT (the $TDDFT input contains more information about how spin-flip runs select the target state's Ms by $CONTRL's MULT value). The reference state must be UHF or ROHF, with MULT in $CONTRL at least 3. The target states of the CIS have one lower Ms, after one alpha spin in the reference is flipped to beta. Nuclear gradients are possible. Solvent effects are not available for either CIS or SFCIS.


Input Description

$CIS

2-91

It is worthwhile to look at the $TDDFT input, which is a very similar calculation. The TD-DFT program offers the possibility of recovering some of the correlation energy, permits some solvent models, and can be used for MEX/CONICL surface intersection searches. The first six keywords are chemically important, while the remainder are mostly technical. NACORE = n Omits the first n occupied orbitals from the calculation (frozen core approximation). For CITYP=CIS, the default for n is the number of chemical core orbitals. For CITYP=SFCIS, the default, which is also the only possibility, is 0. NSTATE = IROOT = Number of states to be found (excluding the reference state). No default is provided. State for which properties and/or gradient will be calculated. Only one state can be chosen. The reference state is referred to as 0, and in the case of CITYP=SFCIS, might have a higher energy than some of the NSTATE target states. Flag to request the determination of CIS level properties, using the relaxed density. Relevant to RUNTYP=ENERGY jobs, although the default is .FALSE. because additional CPHF calculation will be required. Properties are an automatic byproduct of runs involving the CIS or SFCIS nuclear gradient. Type of CI Hamiltonian to use, if CITYP=CIS. SAPS spin-adapted antisymmetrized product of the desired MULT will be used (default) DETS determinant based, so both singlets and triplets will be obtained. Multiplicity (1 or 3) of the singly excited SAPS (the reference can only be singlet RHF). Only relevant for SAPS-based CITYP=CIS run, as SFCIS controls the Ms for target states by the value of MULT in $CONTRL.

CISPRP =

HAMTYP = = = MULT =


Input Description

$CIS ------------

2-92

DIAGZN = = =

Hamiltonian diagonalization method. DAVID use Davidson diagonalization. (default) FULL construct the full matrix in memory and diagonalize, thus determining all states (not recommended except for small cases). Flag to control whether approximate diagonal elements of the CIS Hamiltonian (based only on the orbital energies) are used in the Davidson algorithm. Note, this only affects the rate of convergence, not the resulting final energies. If set .FALSE., the exact diagonal elements are determined and used. Default=.TRUE. Dimension of the Hamiltonian submatrix that is diagonalized to form the initial CI vectors. The default is the greater of NSTATE*2 and 10. Maximum number of expansion basis vectors in the iterative subspace during Davidson iterations, before the expansion basis is truncated. The default is the larger of 8*NSTATE and NGSVEC. Maximum number of Davidson iterations. Default=50. Convergence criterion for Davidson eigenvectors. Eigenvector accuracy is proportional to DAVCVG, while the energy accuracy is proportional to its square. The default is 1.0E-05. Chooses type of CPHF solver to use. CONJG selects an ordinary preconditioned conjugate gradient solver. (default) DIIS selects a diis-like iterative solver. Flag to read CIS vectors from a $CISVEC input group in the input file. Default is .FALSE. Flag to force the use of the minimal amount of memory in construction of the CIS Hamiltonian diagonal elements. This is only relevant when DGAPRX=.FALSE., and is meant for debug purposes.

DGAPRX =

NGSVEC =

MXVEC

=

NDAVIT = DAVCVG =

CHFSLV = = = RDCISV = MNMEDG =


Input Description

$CIS

2-93

The default is .FALSE. MNMEOP = Flag to force the use of the minimal amount of memory during the Davidson iterations. This is for debug purposes. The default is .FALSE.

==========================================================

$CISVEC

group

required if RDCISV in $CIS is chosen

This is formatted data generated by a previous CIS run, to be read back in as starting vectors. Sometimes molecular orbital phase changes make these CI vectors problematic. ==========================================================


Input Description

$MP2

2-94

==========================================================

$MP2

group

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

Controls 2nd order Moller-Plesset perturbation runs, if requested by MPLEVL in $CONTRL. MP2 is implemented for RHF, high spin ROHF, or UHF wavefunctions, but see also $MRMP for MCSCF. Analytic gradients and the first order correction to the wavefunction (i.e. properties) are available for RHF, ROHF (if OSPT=ZAPT), and UHF. The $MP2 input group is not usually given. See also the DIRSCF keyword in $SCF to select AO integral direct MP2. The spin-component-scaled MP2 (SCS-MP2) energy of Grimme is printed for SCFTYP=RHF references during energy runs. See also the keyword SCSPT below. Only the CODE=IMS program is able to do analytic gradients for SCS-MP2. Special serial codes exist for RHF or UHF MP2 energy or gradient, or the ROHF MP2 energy. Parallel codes using distributed memory are available for RHF, ROHF, or UHF MP2 gradients. In fact, the only way that ROHF MP2 gradients can be computed on one node is with the parallel code, using MEMDDI! MP2 energy values using solution models are computed by using the solvated SCF orbitals in the perturbation step. All of the MP2 nuclear gradient programs contain additional terms required for EFP, PCM, EFP plus PCM, or COSMO solvation models. NACORE = n Omits the first n occupied orbitals from the calculation. The default for n is the number of chemical core orbitals. NBCORE = MP2PRP= Same as NACORE, for the beta orbitals of UHF. It is almost always the same value as NACORE. a flag to turn on property computation for jobs jobs with RUNTYP=ENERGY. This is appreciably more expensive than just evaluating the second order energy correction alone, so the default is to skip properties. Properties are always computed during gradient runs, when they are an almost free byproduct. (default=.FALSE.)


Input Description

$MP2

2-95

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. = the program implementation to use, choose from SERIAL, DDI, or IMS according to the following chart, depending on SCFTYP and whether the run involves gradients, UHF UHF ROHF ROHF ROHF energy gradient energy gradient energy OSPT=ZAPT ZAPT RMP SERIAL SERIAL SERIAL SERIAL DDI DDI DDI DDI RIMP2 -

CODE

RHF RHF energy gradient SERIAL DDI IMS RIMP2 SERIAL DDI IMS -

The default for serial runs (p=1) is CODE=IMS for RHF, and CODE=SERIAL for UHF or ROHF (provided PARALL is .FALSE. in $SYSTEM). When p>1 (or PARALL=.TRUE.), the default becomes CODE=DDI. However, if FMO is in use, the default for closed shell parallel runs is CODE=IMS. The "SERIAL" code for OSPT=RMP will run with modest scalability when p>1. The many different MP2 programs are written for different hardware situations. Here N is the number of atomic basis functions, and O is the number of correlated orbitals in the run: The original SERIAL programs use N**3 memory, and have larger disk files and generally takes longer than CODE=IMS. The IMS program uses N*O**2 memory, and places most of its data on local disks (so you must have good disk access), and will run in parallel...ideal for small clusters. Using this program on a node where the disks are of poor quality (SATA-type) and with many cores accessing that single disk may be very I/O bound. Adding more memory can make this program run more efficiently. Network traffic is modest when running in parallel.


Input Description

$MP2

2-96

The DDI program uses N**4 memory, but this is distributed across all nodes, and there is essentially no I/O...ideal for large parallel machines where the manufacturer has forgotten to include disk drives. MEMDDI must be given in $SYSTEM for these codes, so large problems may require many nodes to aggregate enough MEMDDI. The network traffic is high, so an Infiniband quality network or better preferred. Scalability is very good, for example, this program has been used up to 4,000 cores on Altix/ICE equipment. All of the programs just mentioned should generate the same numerical results, so select which one best matches your hardware. The RIMP2 program is an approximation to the true MP2 energy, using the "resolution of the identity" to reduce the amount of data stored (in memory and/or on disk), and also the total amount of computation. See the paper on this program for its reduced CPU and memory requirements. Network traffic is modest. The code has options within the $RIMP2 input to govern the use of replicated memory versus shared memory, as well as the use of disk storage versus distributed memory, so you can tune this to your hardware. References for the various programs are given in REFS.DOC. NOSYM = disables the orbital symmetry test completely. This is not recommended, as loss of orbital symmetry is likely to mean a bad calculation. It has the same meaning as the keyword in $CONTRL, but just for the MP2 step. (Default=0) transformed integral retention threshold, the default is 1.0d-9 (1.0d-12 in FMO runs).

CUTOFF =

The following keyword applies only to RHF references: SCSPT = spin component scaled MP2 energy selection. = NONE - the energy will be the normal MP2 value. This is the default. = SCS - the energy used for the potential surface will be the SCS energy value. Use of SCSPT=SCS causes gradients to be those for the SCSMP2 potential surface. For CODE=IMS, the nuclear gradient can be evaluated analytically. See NUMGRD in $CONTRL if


Input Description

$MP2

2-97

for some reason you wish to use the other two closed shell codes for SCS-MP2 gradients. The following keywords apply to any CODE=SERIAL MP2 run, or to parallel ROHF+MP2 runs using OSPT=RMP: 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, in serial runs only. The default is .FALSE. BASISMO solves the response equations during gradient computations in the MO basis. This is programmed only for RHF references without frozen core orbitals, when it is the default. BASISAO solves the response equations using AO integrals, for frozen core MP2 with a RHF reference, or for ROHF or UHF based MP2. controls memory usage. The default uses all available memory. Applies to CODE=SERIAL. (default=0) 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. Applies only to CODE=SERIAL.

CPHFBS =

=

NWORD =

METHOD= n

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. Applies only to parallel OSPT=RMP runs.

==========================================================


Input Description

$RIMP2

2-98

===========================================================

$RIMP2

group

(optional, relevant if CODE=RIMP2 in $MP2)

This group controls the resolution of the identity MP2 program, which approximately evaluates the MP2 energy. The RI approximation greatly reduces the computer resources required, while suffering only a small error in the energies. Thus, very large atomic basis sets may be used. The input below controls both utilization of the computer resources, and the accuracy of the calculation. See also $AUXBAS, regarding the auxiliary basis set, whose choice also affects the accuracy of the calculation. The program is enabled for parallel calculation, and is tuned to today's SMP nodes. It is limited to energy calculations only, without any solvent effects, for RHF or UHF references. IAUXBF = 0 uses Cartesian Gaussians = 1 uses spherical harmonics for the auxiliary basis set used to expand the MP2 energy expression into products of 3-index matrices. The default is inherited from ISPHER. The next two control computer resources, trading memory for disk storage. GOSMP = flag requesting shared memory use. The default is .TRUE. in multi-core nodes, but .FALSE. in a uniprocessor. This option means only one copy of certain large matrices is stored per node. = a flag to store two and three center repulsion integrals in distributed memory (.TRUE.), or in disk files (.FALSE., which is the default). Selection of this flag requires MEMDDI in $SYSTEM. The default is .TRUE.

USEDM

The RI approximation reduces CPU time, memory requirements, and total disk storage requirements compared to exact calculation. Experimentation with these two keywords will let you tune the program to your hardware situation. For example, choosing GOSMP=.TRUE. and USEDM=.TRUE. will run without any extra disk files, while setting GOSMP=.TRUE.


Input Description

$RIMP2

2-99

and USEDM .FALSE. will minimize memory usage (and network usage) at the expense of doing disk I/O. Total memory usage per node can be obtained by running EXETYP=CHECK. Note the largest replicated memory printed during the RIMP2's output, dividing by 1000000 to get the correct input for MWORDS (round up a bit). Note the largest shared memory requirement printed, also dividing by 100000, and rounding up a bit. Note the distributed memory requirement, which is already in megawords, and is the correct input for MEMDDI. Then, assuming you use p total compute process on multiple n-way nodes, the memory per node is GBytes/node= 8(n*MWORDS + shared + n*MEMDDI/p)/1024 Turning off GOSMP reduces the shared memory to 0 but increases MWORDS, which is multiplied by the number of cores per node! Turning off USEDM leads to MEMDDI=0 by using disk storage instead. If additional memory is available, increasing MWORDS can lead to a reduction in the level of the occupied orbital batch, or "LV". Larger MWORDS permits a smaller LV, which will in turn reduce the required computational time, and the required network traffic or disk I/O. The value of LV used is the last line appearing after "CHECKING SIZE OF OCCUPIED ORBITAL BATCH". The next four control numerical accuracy, but see $AUXBAS which is even more influential in regards the accuracy! OTHAUX = flag to orthogonalize the RI basis set by diagonalization of the overlap matrix. If there is reason to suspect linear dependence may exist in the RI basis, select this option to have a more numerically stable result. Larger RI basis sets such as CCT and ACCT, in particular, may benefit from selecting this. (default=.FALSE.) STOL = threshold at which to remove small overlap matrix eigenvectors, ignored if OTHAUX=.FALSE. This keyword is analogous to QMTTOL in $CONTRL for the true AO basis. (default= 1.0d-6) = selects the procedure for removing redundancies when inverting the two-center, two-e- matrix.

IVMTD


Input Description

$RIMP2

2-100

= 0 use Cholesky decomposition (default) = 2 use diagonalization VTOL = threshold at which to remove redundancies. This is ignored unless IVMTD=2 (default= 1.0d-6)

Don't forget to see also the $AUXBAS input group! An example of this program follows. The molecule is taxol, with 1032 AOs and MOs in the 6-31G(d) basis, correlating 164 valence orbitals. The RI basis set used is SVP, which matches the true basis set in quality. There are 4175 AOs in the RI basis. The job was run on a single 8-way node (n=8, p=1,2,4,8), using MWORDS=50 (leading to LV=6), MEMDDI=580, and the largest shared memory needed is 95 million words. The total node memory is thus (8 bytes/word)*(8*50 + 95 + 8*580/ 8)/1024 = 8.4 GBytes easily fitting into a modern 16 GByte node. It reduces to (8 bytes/word)*(8*50 + 95 + 8*580/16)/1024 = 6.1 GB/node if two 8-way nodes are used. Scaling is p SCF RI-MP2 job total 1 7391 7919 15366 2 3718 4131 7860 4 1857 2290 4174 8 952 1488 2479 16 486 758 1276 using two 8-way nodes. numerical results are E(RI-MP2)= -2920.607512 versus the exact E(MP2)= -2920.606231 The 0.0013 error should be measured against the total 2nd order correlation energy, which is -8.7855, while noting the time for the 2nd order E is similar to the SCF time. ===========================================================


Input Description

$AUXBAS

2-101

===========================================================

$AUXBAS

group

(required if CODE=RIMP2 in $MP2)

This group specifies the auxiliary basis set used to define the resolution of the identity in the RI-MP2 method. The RI methods are formally exact if the RI basis set is complete, so selecting larger bases improves the results. However, this also increases the computational cost of the run! It is reasonable to use smaller RI basis sets when the AO basis is modest, and increase the RI basis when you use very large AO bases. CABNAM specifies built-in basis sets for the RI: = SVP Ahlrich's SVP basis, available H-Kr = TZVP Ahlrich's TZVP basis, available H-Ar = TZVPP Ahlrich's TZVPP basis, available H-Ar = CCD cc-pVDZ basis, available H-Ar = ACCD aug-cc-pVDZ basis, available H-Ar = CCT cc-pVTZ basis, available H-Ar = ACCT aug-cc-pVTZ basis, available H-Ar = XXXXX externally defined: see EXTCAB. CABNAM has no default, this is a required input! Note IAUXBF in $RIMP2 for selecting spherical harmonics versus Cartesian Gaussians. EXTCAB = flag to read the basis from an external file. (default is .FALSE.) This is analogous to EXTBAS in $BASIS: no external files are provided with GAMESS. The value for CABNAM=XXXX must be 8 or fewer letters: avoid the name of any built in auxiliary basis. Your XXXX bases will be read from a file defined by environment variable EXTCAB, in the execution script. Every atom present in your molecule must be defined in the external file by a line giving its chemical symbol and then XXXX. Following this header, give the basis in free format $DATA style, containing only S, P, D, F, and G shells, terminating each atom by the usual blank line. ===========================================================


Input Description

$CCINP

2-102

==========================================================

$CCINP

group

(optional, relevant for any CCTYP)

This group controls a coupled-cluster calculation specified by CCTYP in $CONTRL. The reference orbitals may be RHF or high spin ROHF. If this group is not given, as is often the case, all valence electrons are correlated. Several ground state CCTYP choices obey at least a few of the keywords from $EOMINP, so please see that group too. Excited state runs such as CCTYP=EOM-CCSD or CR-EOML read $CCINP to define orbital ranges for the ground state CCSD prior to generating excitations under $EOMINP's control. A number of CCTYP choices have been superceded by more advanced equations. For example, R-CC and CR-CC were developed prior to their CR-CCL replacement, while CR-EOML supercedes CR-EOM. CR-CCL provides a good approximation to the fully iterated CCSDT method, and so is superior to the familiar CCSD(T). A reasonable menu is: RHF ROHF (high spin) -----------------ground states [properties]: CCSD [CCPRP] CCSD [n/a] CCSD(T) n/a CR-CCL CR-CCL excited states [properties]: EOM-CCSD [CCPRPE] EOM-CCSD MULT=1 MULT=1,3,5 or spin-contaminated CR-EOML n/a ionization processes: EA-EOM3a n/a IP-EOM3a n/a MULT=2,4 CR-CCL =left-CCSD + CR-CC(2,3) perturbative triples. CR-EOML=left-EOM-CCSD + CR-CC(2,3) perturbative triples. Parallel computation is possible for RHF references only, and only for CCTYP=CCSD or CCSD(T). Memory use in parallel runs is exotic: use EXETYP=CHECK with PARALL in $SYSTEM set on prints the per node memory requirements.


Input Description

$CCINP

2-103

See the "Further Information" section of this manual for more details about coupled-cluster runs. **** The first four pertain to both RHF and ROHF **** NCORE = gives the number of frozen core orbitals to be omitted from the CC calculation. The default is the number of chemical core orbitals. = the number of frozen virtual orbitals to be omitted from the calculation. (default is 0) = defines the maximum number of CCSD (or LCCD, CCD) iterations. This parameter also applies to ROHF's left CC vector solver, but not RHF's left vector. See MAXCCL for RHF. (default=30) = defines the convergence criterion for the cluster amplitudes, as 10**(-ICONV). The ROHF reference also uses this for its left eigenstate solver, but see CVGEOM in $EOMINP for RHF references. (default is 7, but it tightens to 8 for FMO-CC.)

NFZV MAXCC

ICONV

**** the next group pertains to RHF reference only **** CCPRP = a flag to select computation of the CCSD level ground state density matrix (see also CCPRPE in $EOMINP for EOM-CCSD level excited states). The computation takes significant extra time, to obtain left eigenstates, so the default is .FALSE. except for CCTYP=CR-CCL or CR-EOML, where the work required for properties must be done anyway. This keyword is only available in serial runs.

Notes: CCSD is the only level at which properties can be obtained. Therefore this option can only be chosen for CCTYP=CCSD, CR-CCL, EOM-CCSD, CR-EOM, or CR-EOML. A CCSD run requesting CCPRP=.TRUE. will internally change itself to EOM-CCSD to run the left CCSD, but since NSTATE of $EOMINP will still be zero, this remains a ground state calculation. Note that the convergence criterion for left eigenstates is CVGEOM in $EOMINP, which is set to obtain excitation energies, and may need tightening.


Input Description

$CCINP

2-104

There is little reason to select any of these: MAXCCL = iteration limit on the left eigenstate needed by CCSD properties, or CR-CCL energies. This is just a synonym for MAXEOM in $EOMINP. If you want to alter the left state's convergence tolerance, use CVGEOM in $EOMINP. The right CCSD state's convergence is set by MAXCC and ICONV. NWORD = a limit on memory to be used in the CC steps. The default is 0, meaning all memory available will be used. = defines the restart option. If the value of IREST is greater or equal 3, program will restart from the earlier CC run. This requires saving the disk file CCREST from the previous CC run. Values of IREST between 0 and 3 should not be used. In general, the value of IREST is used by the program to set the iteration counter in the restarted run. The default is 0, meaning no restart is attempted.

IREST

MXDIIS = defines the number of cluster amplitude vectors from previous iterations to be included in the DIIS extrapolation during the CCSD (or LCCD, CCD) iterative process. The default value of MXDIIS is 5 for all but small problems. The DIIS solver can be disengaged by entering MXDIIS = 0. It is not necessary to change the default value of MXDIIS, unless the CC equations do not converge in spite of increasing the value of MAXCC. AMPTSH = defines a threshold for eliminating small cluster amplitudes from the CC calculations. Amplitudes with absolute values smaller than AMPTSH are set to zero. The default is to retain all small amplitudes, meaning fully accurate CC iterations. Default = 0.0. **** the next group pertains to ROHF reference only **** There is little reason to select any of these. MULT = spin multiplicity to use in the reference determinant during the CC computation.


Input Description

$CCINP

2-105

The value of MULT given in the $CONTRL input determines the spin state for the ROHF orbital optimization, and is the default for the CC. It would be quite unusual to use a different spin in the SCF than in the CC. The MULT keyword in $EOMINP is of greater physical interest. IOPMET = method for the CR-CC(2,3) triples correction. = 0 means try 1 and then try 2 (default) = 1, the high memory option This option uses the most memory, but the least disk storage and the least CPU time. = 2, the high disk option This option uses least memory, by storing a large disk file. Time is slightly more than IOPMET=1, but the disk file is (NO**3 * NU**3)/6 words, where NO = correlated orbitals, and NU= virtuals. = 3, the high I/O option This option requires slightly more memory than 2, and slightly more disk than 1, but does much I/O. It is also the slowest of the three choices. Check runs will print memory needed by all three options. KREST = 0 fresh start of the CCSD equations (default) = 1 restart from AMPROCC file of a previous run

KMICRO = n performs DIIS extrapolation of the open shell CCSD, every n iterations (default is 6) Enter 0 to avoid using the DIIS converger. LREST = 0 fresh start of the left CCSD equations (default) = 1 restart from AMPROCC file of a previous run

LMICRO = n performs DIIS extrapolation of the open shell left equations, every n iterations (default is 5) Enter 0 to avoid using the DIIS converger. KMICRO and LMICRO are ignored for trivial problem sizes. ==========================================================


Input Description

$EOMINP

2-106

==========================================================

$EOMINP

group (optional, for CCTYP=EOM-CCSD, CR-EOM, or CR-EOML) (optional, for CCTYP=EA-EOMx or IP-EOMx) (optional for CCSD properties, or CCTYP=CR-CCL)

This group controls the calculation of excited states by the equation of motion coupled cluster with single and double excitations, with optional triples corrections. It also pertains to electron attachment and detachment processes, which may result in the system being left in an excited state. EOM-CCSD can be selected for RHF or ROHF reference states, while all other CCTYP listed above can be used only with SCFTYP=RHF. These EOM-type runs consist of an SCF calculation on the reference state, followed by a ground state CCSD (see the $CCINP input to control the ground state calculation, and the orbital range correlated), followed by an EOM-CCSD, IP-EOMCC, or EA-EOMCC calculations on the target states (see NSTATE below). In some cases, triples corrections based on the method of moments approach may follow. The input group permits selection of how many states are computed (machine time is linear in the number of states). Since the target state default is simplistic (only one excited state in the totally symmetric representation), it is usually necessary to give this group, to select NSTATE and IROOT sensibly. Because this input group is used for several CCTYP calculations, not all keywords are used in every case, or the keywords may have slightly different meanings: a) if the reference type is RHF, and CCTYP is EOM-CCSD, CREOM, or CR-EOML, keywords MULT and NACT are ignored. b) if the reference type is ROHF, and CCTYP is EOM-CCSD, keywords CCPRPE, NACT, MTRIP, MEOM, MCI, MINIT, CVGCI, MAXCI, MICCI are ignored. MULT is sometimes ignored. c) if the reference type if RHF, and CCTYP is an ionization process (EA/IP), keywords CCPRPE, MTRIP, MEOM, MCI, MINIT, CVGCI, MAXCI, MICCI are ignored. Additional information on CC and EOM methods can be found in the "Further Information" section of this manual.


Input Description

$EOMINP

2-107

--- spin and space symmetry, and state selection: MULT = Spin multiplicities for the target states. The meaning depends on the particular calculation. The default for cases using $EOMINP's MULT is -1. If any doubly excited states or EA/IP quartets are sought, be sure to select MINIT=1 so that the initial guess states include these.

In case the run is IP-EOM or EA-EOM type, the run will pass through open shell code, although the reference in $CONTRL must be given as RHF and MULT=1. The IP or EA states will be spin adapted. MULT = -1 means target both doublet and quartet states = 2 means consider only doublet states = 4 means consider only quartet states, which can be produced at the EOM-CCSD level by a double that unpairs two electrons, and attaches (or detaches) a third electron. In case the run is EOM-CCSD, SCFTYP=ROHF, but $CONTRL selects a closed shell reference by MULT=1: MULT = -1 means target singlet, triplet, and pentuplets. = 1 consider only singlet excited states = 3 consider only triplet excited states = 5 consider only pentuplet excited states. In case the run is EOM-CCSD, SCFTYP=ROHF, and $CONTRL selects a genuinely open shell reference, the EOM states will not be perfectly spin adapted. The spin projection of Ms from $CONTRL's (MULT-1)/2 is the only good quantum number. The excited states will have S values near to Ms, Ms+1, or Ms+2 since single and double excitations are treated. Note that target states with spins LOWER than Ms are not generated, even if they exist in nature. The output will not print approximate S or values, but will label spatial symmetry. MULT = input will be ignored In case the run is an RHF reference EOM-CCSD (or triples correction), the target states are singlets, only. MULT = input will be ignored


Input Description GROUP

$EOMINP

2-108

the name of the Abelian group to be used, which may be only one of the groups shown in the table below. The default is taken from $DATA, and is reset to C1 if the group is non-Abelian. The purpose is to let the Abelian symmetry be turned off by setting GROUP=C1, if desired. Symmetry is used to help with the initial excited state selection, for controlling the EOMCC calculations, and for labeling the calculated states in the output (not to speed up the calculations). an array of up to 8 integers telling how many excited states of each symmetry type should be computed. The default is NSTATE(1)=1,0,0,0,0,0,0,0 meaning 1 totally symmetric state is to be found. The ground state is always computed, and MUST NOT be included in NSTATE's input, for excited state runs. For EA or IP runs, the NSTATE input MUST include the target ion's ground state, and may include excited states of the ion. See also ISELCT below.

NSTATE

There is no particular reason the first excited state (or ionic ground state) should be totally symmetric, so most runs should give a sensible NSTATE input. Up to 10 states can be found in any irrep. Machine time is linear in the number of states to be found, so be realistic! NSTATE uses this order for irreducible representations: irrep 1 2 3 4 5 6 7 8 C1 A C2 A B Cs A' A'' Ci Ag Au C2v A1 A2 B1 B2 C2h Ag Au Bg Bu D2 A B1 B2 B3 D2h Ag Au B1g B1u B2g B2u B3g B3u As an aside, NSTATE(1)=0,0,0,0,0,0,0,0 for RHF references will calculate the ground state only, generating the type I, II, or possibly III CR-CCSD(T) energies, which aren't otherwise available in a direct ground state calculation.


Input Description

$EOMINP

2-109

IROOT

selects the state whose energy is to be saved for further calculations, such as numerical gradients, or whose properties are evaluated, see CCPRPE below. The first integer lists the irrep number, from the same table as NSTATE, and the second lists the number of the state. Thus, IROOT(1)=3,2 means the second B1 state, if GROUP=C2V. (default IROOT(1)=1,0)

IROOT's default is moderately sensible for a RHF-based excited state run, corresponding to the ground state (labeled as state 0), as this state must lie in the totally symmetric representation. ROHF-based excitations, or EA/IP runs should select something appropriate! If degenerate EOM-CCSD states are detected, only one such state will be triples-corrected. The state chosen for possible triples will be the lower irrep number, so make sure IROOT matches this. ISELCT = an array allowing experts to reduce the number of states that are actually solved for. When given, NSTATE determines the number of states generated by the initial guess procedures, with ISELECT selecting those which carry into the calculation. NSTATE(1)=2,2,2,2 with ISELCT(1)=1,3,5,7 prepares two guesses in each irrep, but only iterates the EOM-CCSD equations for the lowest state in each irrep (the guesses are counted serially).

The next two keywords address triples corrections. Note that non-iterative triples corrections are not presently available for any SCFTYP=ROHF reference. MTRIP selects the type of noniterative triples corrections to SCFTYP=RHF EOM-CCSD energies. MINIT applies only to CCTYP=CR-EOM or CR-EOML: 1 = compute the CR-EOMCCSD(T) triples corrections termed type I and II in the output. This is the default, which skips the iterative CISD calculations needed to construct the CR-EOMCCSD(T) triples corrections of type III. 2 = after performing an additional CISD calculation,


Input Description

$EOMINP

2-110

evaluate all types of the CR-EOMCCSD(T) triples corrections, including types I, II, and III. This choice of MTRIP uses approximately 50 % more memory, but less CPU time than MTRIP=4. 3 = evaluate the CR-EOMCCSD(T) corrections of type III only. As with MTRIP=2, this calculation includes the iterative CISD calculation, which is needed to construct the type III triples corrections, in addition to the EOMCCSD and CR-EOMCCSD(T) calculations. 4 = carry out MTRIP=1 calculations, followed by MTRIP=3 calculations, thus evaluating all types of the CR-EOMCCSD(T) corrections (types I, II, and III in the output). As with MTRIP=2, this calculation includes the CISD iterations, which are needed to construct the type III triples corrections, in addition to the EOMCCSD and CR-EOMCCSD(T) calculations. NACT pertains only to EA-EOM3A or IP-EOM3A runs: NACT = the number of active MOs used to select the 3p2h or 3h2p excitations in EA-EOMCCSDt (EA-EOM3A) or IP-EOMCCSDt (IP-EOM3A) calculations. For CCTYP=EA-EOM3A, used to describe the (N+1) esystem, NACT refers to the NACT lowest unoccupied orbitals of the N e- reference system. For CCTYP=IP-EOM3A, used to describe the (N-1) esystem, NACT refers to the NACT highest occupied orbitals of the N e- reference system. The default for NACT is 0, which allows no three particle or three hole operators, and thus yields only EA-EOMCC2 or IP-EOMCC2 results. In other words, you should input a value for NACT!

CCPRPE = a flag to select computation of the EOM-CCSD level excited state density matrices (see also CCPRP in $CCINP for ground states). The computation takes extra time, to obtain left eigenstates, so the default is .FALSE. CCPRPE can be used only if SCFTYP=RHF and CCTYP=EOM-CCSD, CR-EOM, or CR-EOML. The property printout includes transition moments and oscillator strengths between all pairs of states, as well as the full range of Gaussian


Input Description

$EOMINP

2-111

properties (see $ELMOM, etc), for state IROOT only. CC density matrices are square, not symmetric, which means that CCSD natural orbitals come in left/right pairs. To minimize the amount of output, only left natural orbitals for excited state IROOT will be found in the log file. --- iterative solver selection: MEOM selects the solver for the EOMCCSD calculations: 0 = one EOMCCSD root at a time, united iterative space for all calculated roots (default) 1 = one root at a time, separate iterative space for each calculated root 2 = the Hirao-Nakatsuji multi-root solver 3 = one root at a time, separate iterative space for all computed right/left roots. (compare to 1) 4 = one root at a time, united iterative spaces for each right/left root (compare to 0). For open shell references, or IP/EA runs, there is only one EOM-CCSD solver, so MEOM is ignored. MEOM=0,1,2 obtain all the right eigenvectors first, and then if properties are being computed, proceed to compute the left eigenvectors. MEOM=3,4 obtain right and left eigenvectors simultaneously, and therefore should only be chosen if you are computing properties (see CCPRP/CCPRPE). MCI selects the solver for the CISD step, which is irrelevant unless MTRIP is bigger than 1. 1 = one root at a time, separate iterative space for each calculated root (default) 2 = the Hirao-Nakatsuji multi-root solver (slower)

--- initial guess for EOM-CCSD (and possible CISD) solvers: For both MINIT and MACT below, S and D stand for using all singles or doubles, while s and d mean restricting those excitations, both from and to a smaller number of orbitals. Of course, to define the range of orbitals "active" in the initial guess, inputs NOACT and NUACT (and perhaps MOACT) below must be given. The reason that MINIT=1 is preferred is that low-lying states with non-negligible double excitation character, or significant multi-configurational


Input Description

$EOMINP

2-112

character are missing in a simple CIS guess, and thus may not appear in the final converged calculations. MINIT selects the initial guess procedure for EOM-CCSD, and possibly CISD iterations, in case MTRIP>1. MINIT applies to all runs reading $EOMINP. 1 = Use EOMCCSd to start the EOMCCSD iterations, and CISd to start possible CISD iterations. 2 = Use CIS wave functions to start both EOMCCSD, or any possible CISD calculations. MINIT's default is 2, but MINIT=1 is HIGHLY RECOMMENDED! = fine tuning For MINIT=1 For MINIT=1 For MINIT=2 For MINIT=2 The default of MINIT's EOM-CCSD initial guess MACT=0, use EOMCCSd guess MACT=1, use EOMCCsd guess MACT=0, use CIS guess MACT=1, use CIs guess for MACT is 0.

MACT

MINIT applies to all calculations reading $EOMINP, while MACT applies only if SCFTYP=ROHF. the next three define the initial guess active space: There are no default values of NOACT and NUACT, so the user MUST provide NOACT and NUACT values if they are needed. NOACT and NUACT are usually small (5 or so), but should be chosen to avoid splitting any degenerate orbital shells. NOACT NUACT MOACT the number of occupied MOs in the active space for little s or little d initial guesses. the number of unoccupied MOs in the active space for little s or little d initial guesses. array allows explicit selection of the active orbitals used to define the EOMCCSd and CISd initial guesses. If not provided, the MOACT array is filled such that the NOACT highest occupied and NUACT lowest unoccupied orbitals are selected. If MOACT is given, the number of values provided must be NOACT+NUACT. MOACT is most useful in the virtual space, where the lowest orbitals might be diffuse in nature. An example with 15 occupied orbitals, and where the user has searched the virtual space looking for valence-like orbitals, might be


Input Description

$EOMINP MINIT=1 NOACT=3 NUACT=5 MOACT(1)=13,14,15, 19,20,24,25,30

2-113

--- iteration control: CVGEOM MAXEOM convergence criterion on the EOMCCSD excitation amplitudes R1 and R2 (default=1.0d-4). maximum number of iterations in the EOMCCSD calculations (default=50). For MEOM=0 or 1, this is the maximum number of iterations per each calculated state. For MEOM=2, this is the maximum number of iterations for all states of the EOMCCSD multi-root procedure. maximum number of microiterations in the EOMCCSD calculations (default=80). Rarely used. For MEOM=1 (separate iterative space for each root), this is the maximum number of microiterations for each calculated state. For MEOM=0 or 2 (united iterative space for all calculated roots), this is the maximum number of microiterations for all calculated states. It is much better to perform calculations with MICEOM > MAXEOM (i.e., in a single iteration cycle). If for some reason the EOMCCSD convergence is very slow and the iterative space becomes very large, it may be worth changing the default MICEOM value to MICEOM < MAXEOM to reduce the disk usage. This is not going to happen too often and normally there is no need to change the default MICEOM value.

MICEOM

The next three apply only to closed shell reference triples, if the triples method MTRIP is greater than 1: CVGCI MAXCI convergence criterion for the CISD expansion coefficients (default=1.0d-4). maximum number of iterations in the CISD calculation (default=50). For MCI=1, this is the maximum number of iterations per each calculated CISD state. For MCI=2, this is the maximum number of iterations for all states of the CISD multi-root procedure. maximum number of microiterations in the

MICCI


Input Description

$EOMINP

2-114

CISD calculation (default=80). Rarely used. For MCI=1 (separate iterative space for each root), this is the maximum number of microiterations for each calculated state. For MCI=2 (united iterative space for all calculated roots), this is the maximum number of microiterations for all calculated states. In analogy to MICEOM, it is much better to perform the CISD calculations with MICCI > MAXCI (i.e., in a single iteration cycle). ---- restarts: JREST = 0 this is not a restart = 1 restart data is read from AMPROCC file One use for this is to request additional states, with the restart taking any converged roots from disk, and doing an initial guess for additional states. You must not change MULT when restarting.

==========================================================


Input Description

$MOPAC

2-115

==========================================================

$MOPAC

group

(relevant if GBASIS=PM3, AM1, or MNDO)

This group affects only semi-empirical jobs, which are selected in $BASIS by keyword GBASIS. PEPTID = flag for peptide bond correction. By default a molecular mechanics-style torsion potential term is added for every peptide bond linkage found. The intent is to correct these torsions to be closer to planar than they would otherwise be in the semi-empirical model. Here, the peptide bond means any O H \\ / C----N / \ X One such torsion is added for O-C-N-H and one for O-C-N-X. This term is parameterized as in MOPAC6. Default=.TRUE. ==========================================================


Input Description

$GUESS

2-116

==========================================================

$GUESS

group

(optional, relevant for all SCFTYP's)

This group controls the selection of initial molecular orbitals. 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 core potential 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 deck, and you MUST pay attention to NORB below. = RDMINI Read in a $VEC deck from a converged SCF calculation using GBASIS=MINI, to project the MINI orbitals onto the current basis. The option improves upon the Huckel guess because it involves SCF orbitals, which are typically easily obtained in the small MINI basis. This option doesn't work if the current basis uses core potentials. potentials. The $VEC from the MINI run must contain all virtual orbitals. = 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. The next options are less general, being for Fragment Molecular Orbital runs, or Divide and Conquer runs:


Input Description = FMO = HUCSUB = DMREAD

$GUESS

2-117

Read orbitals from the DICTNRY file, from previous FMO run with MODPRP=1. Perform a Huckel guess in each subsystem of a Divide and Conquer run Read a density matrix from a formatted $DM group, produced by a previous Divide and Conquer run, see NDCPRT in $DANDC.

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. (default=.FALSE.) PUNMO = a flag to control punching of the initial guess. (default=.FALSE.) 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.) = The number of orbitals to be read in the $VEC group. This applies only to GUESS=MOREAD.

NORB

For -RHF-, -UHF-, -ROHF-, and -GVB-, NORB defaults to the number of occupied orbitals. NORB must be given for -CIand -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.


Input Description

$GUESS

2-118

IORDER = Reordering instructions, giving the new molecular orbital order. This parameter applies to the common orbitals (both alpha and beta) except for UHF, where IORDER only affects the alpha MOs. Examples (let there be 10 occupied orbitals): transposition of HOMO and LUMO: IORDER(10)=11,10 a different transposition: IORDER(10)=15 IORDER(15)=10 a more general permutation: IORDER(8)=11,8,9,10 so the new orbital 10 is the original 9th. The 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. However it may fail in complicated groups when the orbitals are very unsymmetric. (default=.FALSE.) TOLZ TOLE = level below which MO coefficients will be set to zero. (default=1.0E-7) = 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 initial density in order to generate symmetric orbitals. This may be useful if the


Input Description

$GUESS

2-119

HUCKEL or HCORE guess types give orbitals of impure symmetry (?'s present). The procedure will generate a fairly high starting energy, and thus its use may not be a good idea for orbitals of the quality of MOREAD. (default=.FALSE.) ==========================================================


Input Description

$VEC $DM $MOFRZ

2-120

==========================================================

$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 are $VEC1, $VEC2, ... for each monomer, computed identical orientation as the supermolecule. For moment or spin-orbit coupling runs, orbitals for one and possibly two are $VEC1 and $VEC2. this group in the transition states

==========================================================

$DM

group

(relevant in Divide and Conquer runs)

This group consists of a formatted density matrix, read in exactly the format it was written. See GUESS=DM, and NDCPR in $DANDC. ==========================================================

$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 IFRZ = flag which triggers MO freezing. (default=.FALSE.) = an array of MOs in the input $VEC set which are to be frozen. There is no default for this.

==========================================================


Input Description

$STATPT

2-121

==========================================================

$STATPT

group

(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. NR Pick from

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 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. 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)

RFO

QA

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 manually.


Input Description

$STATPT

2-122

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. The default is 50 steps for a minimum search, but only 20 for a transition state search, which benefit from relatively frequent Hessian re-evaluations. --- 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. (default=0.1)

the next three apply only to METHOD=RFO or QA: TRUPD TRMAX TRMIN = a flag to allow the trust radius to change as the geometry search proceeds. (default=.TRUE.) = maximum permissible value of the trust radius. (default=0.5 for OPTIMIZE and 0.3 for SADPOINT) = minimum permissible value of the trust radius. (default=0.05)

--- 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 referred to as "eigenvector following". For METHOD=SCHLEGEL, the mode whose eigenvalue is (or will be made) negative. All other curvatures will be made positive.


Input Description

$STATPT

2-123

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 positive. Note that STPT=.TRUE. has little meaning with HESS=GUESS as there will be many degenerate eigenvalues.

STSTEP = Stepsize for jumping off a stationary point. Using values of 0.05 or more may work better. (default=0.01) 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, which was for a triatomic $CONTRL NZVAR=3 $END $ZMAT IZMAT(1)=1,1,2, 2,1,2,3, 1,2,3 $END while optimizing the angle. If NZVAR=0, so that this value applies to the Cartesian coordinates instead, the input of IFREEZ(1)=4,8 means to freeze the x coordinate of the 2nd and y coordinate of the 3rd atom. 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. IFREEZ can pertain to EFP particles, but only during RUNTYP=HESSIAN,


Input Description

$STATPT

2-124

where the 6 translational and rotational degrees of freedom of each EFP come AFTER the QM atom coordinates. For more information: J.D.Head, Int.J.Quantum Chem. 65, 827, 1997 H.Li, J.H.Jensen Theoret. Chem. Acc. 107, 211-219(2002) IFCART = array of Cartesian coordinates to freeze during a geometry optimization using delocalized internal coordinates. This probably works less well than IFREEZ when it freezes Cartesians. Only one of IFREEZ or IFCART may be chosen in a single run. IACTAT = array of "active atoms", which is a complimentary input to IFREEZ. Any atom *not* included in the list has its Cartesian coordinates frozen. Thus IACTAT(1)=3,-5,107,144,202,-211 allows 15 atoms, namely 3-5, 107, 144, and 202-211 to be optimized, while all other atoms are frozen. NZVAR in $CONTRL must be 0 when this option is chosen. IFREEZ and IACTAT are mutually exclusive. The latter acts by generating a IFREEZ for all atom coordinates not defined as "active", so users can input whichever list is shorter. --- The next two control the hessian matrix quality --HESS = selects the initial hessian matrix. = GUESS chooses an initial guess for the 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 compute the hessian, see $FORCE input. = 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


Input Description

$STATPT

2-125

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. (default=.FALSE.) --- 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 0 -1 -2 3 2 1 0 -1 -2 Print Print Print Print INTR INTR INTR INTR at at at at all, all, all, 0+L, orbitals orbitals orbitals orbitals at all at 0+L (default) never never

NPUN

=

Punch all orbitals and HESS at all Punch all orbitals at all same as 0, plus punch HESS at all Punch all orbitals at 0+L, otherwise only occupied orbitals (default) Punch occ orbitals at 0+L only Never punch orbitals

---- the next parameters control harmonic constraints --Harmonic constraints can be added to the current geometry by setting ALL the keywords below. For instance, to harmonically constrain the distance between atom 3 and 12 to a distance of 2.0 Angstrom and a force constant of 500 kcal/mol, the following example can be used: IHMCON(1)=1,3,12 SHMCON(1)=2.0 FHMCON(1)=500.0 The default is all zeros which means do not do this. IHMCON = array of coordinates to constrain. The input is similar to IZMAT in $ZMAT, a code integer, and the atoms involved in the coordinate. The code integer may only be 1, for stretches. SHMCON = equilibrium constraint values for the distances


Input Description

$STATPT

2-126

specified by IHMCON, given in Angstrom. FHMCON = array of force constants for the distances specified by IHMCON, given in kcal/mol. ---- 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 BFGS POWELL POWELL MSP SCHLEGEL do not update Hessian (not recommended) default for OPTIMIZE using RFO or QA default for OPTIMIZE using NR or CONOPT default for SADPOINT mixed Murtagh-Sargent/Powell update only choice for METHOD=SCHLEGEL

---- 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 RUNTYP=SADPOINT. = Minimum distance threshold. Points whose root mean square distance from the current point is less than RMIN are discarded. (default=0.0015)

RMIN


Input Description

$STATPT

2-127

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) ========================================================== ********************* See the 'further information' section for some help with OPTIMIZE and SADPOINT runs *********************


Input Description

$TRUDGE

2-128

==========================================================

$TRUDGE

group

(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 parameters. 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 inputs). 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 input. 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 IEX = number of parameters to be optimized. = 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.


Input Description

$TRUDGE

2-129

(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, respectively. 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, $CONTRL SCFTYP=GVB CITYP=GUGA RUNTYP=TRUDGE COORD=HINT $END $BASIS GBASIS=N31 NGAUSS=6 $END $DATA 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 $END $SCF NCO=3 NPAIR=1 $END $TRUDGE OPTMIZ=GEOMETRY NPAR=2 IEX(1)=21,22 P(1)=1.08 $END $CIDRT GROUP=C2V SOCI=.TRUE. NFZC=1 NDOC=3 NVAL=1 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. ==========================================================


Input Description

$TRURST

2-130

==========================================================

$TRURST

group

(optional, relevant if 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 ==========================================================


Input Description

$FORCE

2-131

==========================================================

$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 inputs. Calculation of the hessian automatically yields the dipole derivative tensor, giving IR frequencies. Raman intensities are obtained by following with RUNTYP=RAMAN. METHOD = chooses the computational method: = ANALYTIC is a fully analytic calculation. This is implemented for SCFTYP=RHF, UHF, ROHF, GVB (for NPAIR=0 or 1, only), and MCSCF (for CISTEP=ALDET or ORMAS, only). R-DFT and U-DFT are also analytic. This is the default for these cases. = SEMINUM does numerical differentiation of analytically computed first derivatives. This is the default for UHF, MCSCF using other CISTEPs, DFT, all solvent, models, relativistic corrections, and most MP2 or CI runs. = FULLNUM numerically differentiates the energy twice, which can be used by all other cases. It requires many energies (a check run will tell how many) and so it is mainly useful for systems with only very few symmetry unique atoms. The default for METHOD is to pick ANALYTIC over SEMINUM if that is programmed, and SEMINUM otherwise. FULLNUM will never be chosen unless you specifically request it. RDHESS = a flag to read the hessian from a $HESS input, rather than computing it. This variable pertains only to RUNTYP=HESSIAN. See also HESS= in the


Input Description

$FORCE (default is .FALSE.)

2-132

$STATPT input group.

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 provided $ZMAT input to use this. (Default=.FALSE.) --- the next four apply to numeric differentiation ---NVIB = The number of displacements in each Cartesian direction for force field computation. This pertains only to METHOD=SEMINUM, as FULLNUM always uses double difference formulae. 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. 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. (default) Displacement size (in Bohrs). This pertains to Both SEMINUM and FULLNUM. Default=0.01

=1

=2

VIBSIZ =


Input Description

$FORCE

2-133

Let 0 mean the Vib0 geometry, and D mean all the displaced geometries NPRT NPUN =1 =0 =2 =1 =0 Print orbitals at 0 and D Print orbitals at 0 only (default) Punch all orbitals at 0 and D Punch all orbitals at 0 and occupied orbs at D 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 J.P.Merrick, D.Moran, L.Radom J.Phys.Chem.A 111, 11683-11700 (2007). 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. = 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.

FREQ

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


Input Description

$FORCE

2-134

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 input may define more than 3N-6 (3N-5) coordinates. (default=.FALSE.) 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 is punched both before and after projection. For the same reason, the default is .FALSE. to skip the projection, which is mainly of interest in dynamical calculations. ========================================================== There is a program ISOEFF for the calculation of kinetic and equilibrium isotope effects from the group of Piotr Paneth at the Technical University of Lodz. This program will accepts data computed by GAMESS (and other programs), and can be requested from paneth@p.lodz.pl


Input Description

$CPHF

2-135

==========================================================

$CPHF

group

(relevant for analytic RUNTYP=HESSIAN)

This group controls the solution of the response equations, also known as coupled Hartree-Fock (CPHF). 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.) CPHF = nature of integrals driving the formation of the response equation. The defaults are intelligent, so this is meant mostly for experts/debugging. = AO forms response equations from AO integrals, which usually takes less memory, and is more parallel. This is the default for RHF, UHF, ROHF, R-DFT, or U-DFT. AO-driven mode is not available for any other cases. = MO forms response equations from transformed MO integrals. This is the default for GVB or MCSCF. This is available forRHF or ROHF. = AODDI forms response equations from AO integrals, using distributed memory (see MEMDDI). This does AO integrals about 2x more than AO, but spreads the CPHF memory requirement out across multiple nodes. Coded only for RHF.

SOLVER = linear equation solver choice. This is primarily a debugging option. For RHF analytic Hessians, choose from CONJG (default), DIIS, ONDISK, not all of which will work for all CPHF= choices. For imaginary frequency dependent polarizability responses (MAKEFP jobs), choose GMRES (default), biconjugate gradient stabilized BCGST, DODIIS, or an explicit solver GAUSS. Most response equations have only one solver programmed, and thus ignore this keyword. NWORD = controls memory usage for this step. The default uses all available memory. (default=0) ==========================================================


Input Description

$CPMCHF

2-136

==========================================================

$CPMCHF

group (relevant for analytic RUNTYP=HESSIAN,NACME,CONICAL)

This group controls the solution of the response equations, also called coupled perturbed multiconfiguration Hartree-Fock, for MCSCF wavefunctions. These are needed for analytic hessians (CISTEP=ALDET or ORMAS), or for state-averaged gradients or non-adiabatic coupling matrix element calculations. The full response equations are solved for hessians, while Z-vector equations are used for NACME and SA-gradients. The default converger is a (linear) conjugate gradient (CG) method, but three others may be chosen. Difficult cases might work upwards from the default CG method by: $cpmchf $end $cpmchf ipdir=50 $end $cpmchf gcro=.t. micit=5 kicit=10 $end $cpmchf gcrodr=.t. micit=10 kicit=5 $end $cpmchf gcrodr=.t. micit=30 kicit=10 reclin=.false. $end $cpmchf gcrodr=.t. micit=20 kicit=10 reclin=.false. prcchg=.true. prctol=1.0 $end $cpmchf gcro=.t. micit=50 kicit=100 prcchg=.true. prctol=1.0 $end where the last one is "the sledgehammer". The options shown in the next to last case very often work, and will be considerably faster than the last set, which should always work. GCRO will typically take many fewer iterations than CG, a measure of its robustness, but will use more machine time due to its microiterations. --- the next set apply to any CP-MCHF converger --MAXIT CPTOL = maximum iterations. (default=300)

= accuracy tolerance for cpmchf equations Ax=b. (compared to r/||b||, within orbital, CI, and state averaged components) (default=1.0D-07) = a flag to adjust the linear equation's preconditioning. (Default=.FALSE.) For ORMAS runs in particular, the standard

PRCCHG


Input Description

$CPMCHF

2-137

preconditioner might lead to ill-conditioning with very small elements. If selected, any preconditioner elements below PRCTOL will be reset to a value of 1.0 PRCTOL = a tolerance to keep preconditioner elements, if PRCCHG is chosen. Default=1.0D-6 --- for RUNTYP=NACME --NAPICK = a flag to select running through z-vector setup and choose which linear equations to solve. .TRUE. leads to z-vector linear equations. .FALSE. leads to non-z-vector linear equations. The z-vector equations are advantageous when the degrees of freedom exceeds the no. of electronic states involved in the state-averaging. (The defaults are set to enforce this option.) The z-vector equations are also advantageous when only a few NA couplings out of the total total possible couplings are of interest. If CISTEP=ORMAS, NAPICK=.TRUE. is forced. Selecting NAPICK=.TRUE. requires the choice of NA couplings in the NACST array (see below). (default=varies... see ROUTINE NACMEX for notes) = an array that indicates which NA couplings to Calculate, if NAPICK is chosen. For example, if WSTATE in $DET contains at least four states, the NACME can be limited to state pairs 1<->2 and 3<->4 by NACST(1)=3,4, 1,2. Note that you should pick increasing order within any pair of states! The program always generates the state-specific gradient of every state with a non-zero WSTATE. (default=none)

NACST

--- the next three choose the other CPMCHF convergers --If none is selected, conjugate gradient (CG) is used. GCRODR = a flag to select Generalized Conjugate Residual with inner Orthogonalization with Deflated Restarting. (default=.FALSE.) = a flag to select Generalized Conjugate Residual

GCRO


Input Description

$CPMCHF

2-138 (default=.FALSE.)

with inner Orthogonalization. CGGCRO

= a flag to alternate between GCRO and CG solvers. (default=.FALSE.)

--- next apply only to GCRODR: RECLIN = a flag to select recycling of Krylov subspaces from the first linear equation. (recycle linear). The first CP-MCHF linear equation is solved and a recyclable subspace is generated. Then, after a projection of approximate solution across the subspace from the first system, the rest of the linear equations are solved. The recycled subspace may or may not give rapid convergence with fewer iterations. See MICIT and KICIT. (default=.TRUE.) --- next apply only to GCRODR,GCRO,CGGCRO: MICIT = total size of the Krylov expansion space, namely the number of micro-iterations within an overall iteration. While the MICIT variable has no limit, fifty or more micro-iterations start to become computationally unmanagable for larger systems. The default often must be increased for large systems or for geometries far from equilibrium. In addition, the GCRODR converger has a slightly modified scheme for the micro-iterations. In the first iteration, MICIT micro-iterations are performed in a GMRES(MICIT) iteration. However, for subsequent iterations, (MICIT-KICIT) micro-iterations are performed. (default=5, often increased to 20 or 30) = the size of the recyclable Krylov basis saved and modified from iteration to iteration, created from eigenvectors with small eigenvalues. If this space is too small, the run will experience ill-conditioning, but if too large, the search space includes ineffective parts. In case PRCCHG is chosen, make sure the number of vectors reset from small values does not exceed

KICIT


Input Description

$CPMCHF

2-139

KICIT, if it does, increase KICIT. (default=5, usually small, e.g. 5-10 for GCRODR) (MICIT> KICIT for GCRODR) (KICIT>=MICIT for GCRO and CGGCRO) --- next apply only to (linear) CG: NACFAC = number of iterations before softening the current convergence tolerance by 1/3, for CG only. Use NACFAC=0 (or MAXIT) to prevent raising the from the initial CPTOL. (default=50) number of iterations before resetting the residual from a pseudo-residual to the true residual. This reset also resets the search directions, since keeping the old 'ill-rounded' directions is not very beneficial. If a run almost convergences but struggles in later iterations, IPDIR=50 is recommended. This option appears more useful for NACME rather than Hessian runs. (default=MAXIT) --- next apply only to CGGCRO: ITERA = the number of (linear) CG iterations to apply before alternation to the GCRO converger. (default=5)

IPDIR

=

ITERB

= the number of GCRO iterations to apply before alternation to the (linear) CG converger. (default=5) ==========================================================


Input Description

$MASS

2-140

==========================================================

$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. ==========================================================


Input Description

$HESS $GRAD

2-141

==========================================================

$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. $HESS information 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.) This

Formatted gradient vector at the $DATA geometry. 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. ==========================================================


Input Description

$DIPDR $ALPDR

2-142

==========================================================

$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. ==========================================================

$ALPDR

group

(relevant for RUNTYP=RAMAN or HESSIAN)

Formatted alpha polarizability derivative tensor, punched by a previous RUNTYP=RAMAN job. If both $DIPDR and $ALPDR group are found in the input file, the applied electric field computation will be skipped, to immediately evaluate IR and Raman intensities. This restart may be most relevant for isotopic substitution. If this group is found during RUNTYP=HESSIAN, the Raman intensities will be added to the output. You might want to restart as RUNTYP=HESSIAN instead of RUNTYP=RAMAN in order to have access to PROJCT or the other options available in the $FORCE input group. ==========================================================


Input Description

$VIB $VIB2

2-143

==========================================================

$VIB

group

(relevant for RUNTYP=HESSIAN, METHOD=SEMINUM)

Formatted restart data, consisting of energies, gradients, and dipole moments. This data is read in the same format by which is was written to the RESTART file. Just add a " $END" card, and place this group into the input file to effect a restart. If the final displacement's gradient was written as zero, delete the entire last data set (energy, gradient, and dipole). In case the numerical hessian was run in groups (see $GDDI), the $VIB entries will be out of order. The unsorted data can be read back in, if and only if the new run is also using processor groups. Note that assembling a complete $VIB could be sent to one core in one group, but use NGROUP=1 to make it look like a "group run". This group can be used to turn a less accurate single differencing run into a more accurate double differencing run (NVIB in $HESS). The mere presence of this group triggers the restart. ==========================================================

$VIB2

group

(relevant for hessians, METHOD=FULLNUM) (relevant for gradients, with NUMGRD=.TRUE.)

Formatted restart information, consisting of energy values, as written to the RESTART file. Just add a " $END" line at the bottom, and place this group into the input file to effect a restart. This group has the same name ($VIB2), but different contents, depending on whether you are restarting a numerical gradient or a fully numerical hessian job. The mere presence of this group triggers the restart. ===========================================================


Input Description

$VSCF

2-144

==========================================================

$VSCF

group

(optional, relevant to RUNTYP=VSCF)

This group governs the computation of vibrational frequencies including anharmonic effects. Besides the keywords shown below, the input file must contain a $HESS input (and perhaps a $DIPDR input), to start with previously obtained harmonic vibrational information. The VSCF method requires only energies, so any energy type in GAMESS may be used, perhaps with fully numerical harmonic vibrational information. Energies are sampled along the directions of the harmonic normal modes, and usually along pairs of harmonic normal modes, after which the nuclear vibrational wavefunctions are obtained. The dipole on the grid points may be used to give improved IR intensities. The most accurate calculation computes the potential surface directly, on all grid points, but this involves many energy evaluations. An attractive alternative is the Quartic Force Field approximation of Yagi et al., which computes a fit to the derivatives up to fourth order by computing a specialized set of points, after which this fit is used to generate the full grid of points for the solver. Since there are a great many independent energy evaluations, no matter which type of surface is computed, the VSCF method allows for computations in subgroups (much like the FMO method). Thus any $GDDI input group will be read and acted upon, if found. Vibrational wavefunctions are obtained at an SCF-like level, termed VSCF, using product nuclear wavefunctions, along with an MP2-like correction to the vibrational energy, which is termed correlation corrected (cc-VSCF). In addition, vibrational energy levels based on second order degenerate pertubation theory (see VDPT) or a CI analog (see VCI) may be obtained. Most VSCF applications have been carried out with an electronic structure level of MP2 with triple zeta basis sets. This is thought to give accuracy to 50 wavenumbers for the larger fundamentals. Use of internal coordinates is known to give improved accuracy for lower frequencies, particularly in weakly bound clusters.


Input Description

$VSCF

2-145

Restarts involve the $VIBSCF input (which has different formats for each PETYP), and the READV keyword. Restarts are safest on the same machine, where normal mode phases are reproducible. References for the VSCF method, the QFF approximation, and the solvers are given in Chapter 4 of this manual, along with a number of sample applications. ***** The first input variables control the generation of the potential surface on which the nuclear vibrations occur: PETYP = DIRECT computes the full potential energy surface, according to NCOUP/NGRID. The total number of energy/dipole calculations for NCOUP=2 will be M*NGRID + (M*(M-1)/2)*NGRID*NGRID, where M is the number of normal modes. This is the default. = QFF the Quartic Force Field approximation to the potential surface is obtained. This is usually only slightly less accurate, but has a greatly reduced computational burden, namely 6*M + 12*M*(M-1)/2 energy/dipoles.

INTCRD = flag setting the coordinate system used for the grids. Any internal coordinates to be used must be defined in $ZMAT, using 3N-6 simple, DLC, or natural internal coordinates. Of course, you must enter NZVAR in $CONTRL as well. The default is to use Cartesians (default .FALSE.) INTTYP = 0 default if INTCRD=.FALSE. (ignore this keyword) = 1 implies that the $ZMAT contains only stretches, bends, and torsions. It also selects an approximate transformation between Cartesian and internal coords. = 2 the other $ZMAT coordinates may be used, and the coordinate transformation will be iterated to convergence. (default if INTCRD=.TRUE.) NCOUP = the order of mode couplings included.


Input Description

$VSCF

2-146

= 1 computes 1-D grids along each harmonic mode = 2 adds additionally, 2-D grids along each pair of normal modes. (default=2) = 3 adds additionally, 3-D grids for mode triples, for PETYP=DIRECT only. NGRID = number of grid points to be used in solving for the anharmonic vibrational levels. In the case of PETYP=DIRECT, each of these grid points must be explicitly computed. For PETYP=QFF these grid points are obtained from a fitted quartic force field. Reasonable values are 8 or 16 for DIRECT, with 16 considered significantly more accurate. For PETYP=QFF, the generation of the solver grid is very fast, so use 16 always. (default=16) = step size for PETYP=DIRECT displacements. The maximum distance along each mode is a function of its frequency, amplitude(i)=sqrt(2*(AMP+1/2)/freq(i)) so that AMP resembles a vibrational quantum number. The default goes far enough past the classical turning points of the fundamentals to capture the relevant part of the surface. (default = 7.0) = step size for PETYP=QFF displacements. The step along each mode depends on the harmonic frequency, as well as this parameter, whose default is usually satisfactory (default=0.5)

AMP

STPSZ

In case the user wants to control each normal mode with a separate parameter, arrays of values may be given, using the keywords AMPX(1)=xx,yy,... or STPSZX(1)=xx,yy,zz... IMODE = array of modes for which anharmonic effects will be computed. IMODE(1)=10,19 computes anharmonic energies and wavefunctions for modes 10 and 19, only. In the current implementation, pairs of modes cannot be coupled, so NCOUP is forced to 1 if this option is specified. This approximation is intended for larger molecules, where the whole VSCF calculation is prohibitive. *****


Input Description

$VSCF

2-147

The next set of keywords relates to the solver step which finds the vibrational states. The results always include VSCF and cc-VSCF (SCF and non-degenerate MP2-like solutions). Use of the restart option makes comparing the solvers very fast, compared to the time to generate the electronic potential energy surface's points. VDPT = option to use 2nd order degenerate perturbation theory, based on the ground and singly excited vibrational levels. Results for virtual CI within the same singly excited space will also be given. (default=.TRUE.) = option to use the virtual CI solver within a space of the ground and both singly and doubly excited vibrational levels. Selection of VCI turns VDPT off. (default=.FALSE.)

VCI

The solver always finds the ground vibrational state (v=0) by default, and defaults to finding the fundamentals (v=1 in every mode). It can rapidly find excited levels (such as all v=2) if restarted (see READV) from $VIBSCF, using the following to control the excitation levels: IEXC = 1 obtain fundamental frequencies (default) = 2 instead, obtain first overtones = 3 instead, obtain second overtones = 0 skip combination bands (default) = 1 add one additional quanta in other modes = 2 add two other quanta in one mode at a time. IEXC2 0 0 0 0 1 2 1 for H2O, which has only three modes: only 000 ground state, no transitions 000, and 100, 010, 001 (fundamentals) 000, and 200, 020, 002 (1st overtones) 000, and 300, 030, 003 (2nd overtones) 000, and 100, 010, 001, 110, 101, 110 (1st overtones and combinations) 000, and 100, 010, 001, 210, 201, 021 000, and 200, 020, 002, 120, 102, 012 between them, 1st and 2nd overtones, and all 2-1-0 combinations.

IEXC2

IEXC 0 1 2 3 1 1 2


Input Description

$VSCF

2-148

ICAS1, ICAS2 = starting and ending vibrations whose quanta are included. The default is all modes, ICAS1=1 and ICAS2=3N-6 (or 3N-5). SFACT = a numerical cutoff for small contributions in the solver. The default is 1d-4: 5d-3 or 1d-3 may affect accuracy of results, 1d-4 is safer, and 1d-5 might not converge. = 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. It is seldom necessary to select a scaling other than unity. (Default=1.0) ***** The next two relate to simplified intensity computation. These simplifications are aimed at speeding up MP2 runs, if one does not care so much about intensities, and would like to eliminate the considerable extra time to compute MP2level dipoles. DMDR must not be used if overtones are being computed. DMDR = if true, indicates that the harmonic dipole derivative tensor $DIPDR will be read and used, rather than computing dipoles. (default=.FALSE.) = 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 an MP2 level $DIPDR is available. (default=.TRUE.) **** These relate to the initial harmonic mode generation. Normally, a $HESS is provided, from which harmonic modes are obtained. It is possible to give the harmonic data explicitly with the first two: RDFRQ = array of harmonic frequencies, starting from the smallest.

VCFCT

MPDIP


Input Description

$VSCF

2-149

CMODE

= array of normal mode displacements given in the same order as the frequencies read in RDFRQ. The data should be the x,y,z displacement of the first atom of the first mode, then x,y,z for the second atom, then going on to give each additional mode.

PROJCT = controls the projection of the hessian matrix (same meaning as in $FORCE). Default is .TRUE. which removes small mixings between rotations or translations and the harmonic modes. **** 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.)

GEONLY = option to generate all points on the potential energy surface needed by the VSCF routine, without energy evaluations. The purpose of this is to prepare a set of geometries at which the energy is needed. A possible use for this is to obtain energies from a different program package, which might have an energy unavailable in GAMESS, but which lacks its own VSCF program. (default=.false.) ==========================================================

$VIBSCF

group

(optional, relevant to RUNTYP=VSCF)

This is restart data, as written to the disk file RESTART in a complete or partially completed previous run. Append a " $END", and also select READV=.TRUE. to read the data. $VIBSCF's contents are different for PETYP=DIRECT or QFF. The format of this group changed in December 2006, so that old groups can no longer be used. ==========================================================


Input Description

$GAMMA $EQGEOM $HLOWT $GLOWT

2-150

==========================================================

$GAMMA

group

required if RUNTYP=GAMMA

This group governs evaluation of the 3rd derivative of the energy with respect to nuclear coordinates, by finite differentiation of Hessians (see $FORCE options). NFCM = n describes the amount of restart data provided. The default is n=-1, to evaluate everything. A value of n means that n+1 $FCM groups are to be read from the file (hessian #0 means the equilibrium geometry). Restart data is read from a .gamma file, created by an earlier run.

DELTA = step size, default=0.01 Bohr PRTALL = flag to print full Hessian and Gamma matrix, the default is .FALSE. PRTSYM = flag to print unsymmetrical Gamma elements, the default is .FALSE. PRTBIG = flag to print large Gamma elements, default = .F. ==========================================================

$EQGEOM

group

required if NFFLVL=2 or 3 in $CONTRL

The coordinates of the stationary point, where the hessian and possibly 3rd derivative information was evaluated, in exactly the format it was printed by RUNTYP=GAMMA. ==========================================================

$HLOWT $GLOWT

group group

required if NFFLVL=2 or 3 in $CONTRL required if NFFLVL=3 in $CONTRL

These are the lower triangular parts of the hessian and 3rd derivative matrices, read in the same format as printed by an earlier RUNTYP=GAMMA. ==========================================================


Input Description

$IRC

2-151

==========================================================

$IRC

group

(relevant for RUNTYP=IRC)

This group governs the location of the intrinsic reaction coordinate (also called the minimum energy path, MEP), a steepest descent path in mass weighted coordinates, that connects the saddle point to reactants and products. The IRC serves a proof of the mechanism for a reaction, and is a starting point for reaction path dynamics. The IRC may be found for systems with QM atoms, EFP particles, or the combinations of QM and EFP particles, or QM plus the optional SIMOMM plug-in MM atoms. Restart data for RUNTYP=IRC is written into the PUNCH file. Information summarizing the reaction path is written to the TRAJECT file, which should be saved, appending these as various restarts are done. The graphics program MacMolPlt can display a movie of the entire mechanism, if you join the entire forward and entire backwards trajectory files, while changing the path distance parameter in the reverse part to a negative value. ----- 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

RCUT ACUT MXOPT IHUPD

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. 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) maximum angle from end points for linear interpolation (default=5 degrees) maximum number of constrained optimization steps for each IRC point (default=20) is the hessian update formula. 1 means Powell, 2 means BFGS (default=2)


Input Description GA

$IRC

2-152

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 SAB GA selects quadratic gradient following. Related input is:

distance to previous point on the IRC. gradient vector at that historical point. selects the fourth order Adams-Moulton variable step predictor-corrector. Related input is:

PACE = AMPC4

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.


Input Description

$IRC

2-153

----- 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) ----- constraint ----Of course, applying a constraint to the saddle point search and the reaction path means that you are not locating the true saddle, nor following the true reaction path. IFREEZ = array of Cartesian coordinates to freeze. The IRC stepper works in mass-weighted Cartesian space, making it impossible to freeze internal coordinates. An input of IFREEZ(1)=4,8 means to freeze the x coordinate of the 2nd atom and the y coordinate of the 3rd atom, that is, we count coordinates x1,y1,z1,x2,y2,z2,x3,y3,z3,... ----- 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 0 -1 -2 1 0 -1 -2 Print Print Print Print INTR INTR INTR INTR at at at at all, all, all, F+L, orbitals orbitals orbitals orbitals at all IRC points at F+L (default) never never

NPUN

=

Punch all orbitals at all IRC points Punch all orbitals at F+L, only occupied orbitals at IRC points between (default) Punch all orbitals at F+L only Never punch orbitals


Input Description

$IRC

2-154

----- 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 $HESS input, 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


Input Description

$IRC

2-155

GS2 method, which must have a hessian input, even for restarts. Note also that EVIB is ignored by GS2 runs. ** For the ** **************** hints about IRC tracking, see 'further information' section. ****************

==========================================================


Input Description

$DRC

2-156

==========================================================

$DRC

group

(relevant for RUNTYP=DRC)

This group governs "direct dynamics", following the dynamical reaction coordinate, which is a classical trajectory based on quantum chemistry potential energy surfaces. These may be either ab initio or semi-empirical, and are computed "on the fly" as the trajectory proceeds. 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. Restart data can be found in the job's OUTPUT file, with important results summarized to the TRAJECT file. Almost all DRCs break molecular symmetry, so build your molecule with C1 symmetry in $DATA, or specify NOSYM=1 in $CONTRL. RUNTYP=DRC may not be used with EFP particles. NSTEP = The number of DRC points to be calculated, not including the initial point. (default = 1000) (default = 0.1 fs)

DELTAT = is the time step.

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


Input Description

$DRC

2-157

initial

of the velocity vector. Its magnitude will be automatically adjusted to match the desired 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. (default=.FALSE.) The next three 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 $HESS input, which is the nuclear force constant matrix at the starting geometry.

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 when HESS=MIN, and 0.0, 0.5, ..., 0.5 if HESS=TS. 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. RCENG = reaction coordinate energy, in kcal/mol. This is the initial kinetic energy given to the imaginary frequency normal mode when HESS=TS. If this is given as a negative value, the direction of the velocity vector will be the "reverse direction", meaning the phase of the normal mode will be chosen so its largest component is negative. ***


Input Description

$DRC

2-158

The next two 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 $HESS input for the initial geometry. 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 normal mode component is a positive value. (default=0) = 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)

ENM

To summarize, there are 5 ways to initiate a 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 and possibly RCENG, to give some initial kinetic energy to all 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 normal modes of a minimum energy geometry it travels around. NMANAL = a flag to select mapping of the mass-weighted Cartesian DRC coordinates and velocity (conjugate momentum) in terms of normal modes at a nearby


Input Description

$DRC

2-159

reference stationary point (which can be either a minimum or transition state). This reference geometry could in fact be the same as the initial point of the DRC, but does not need to be. If you choose this option, you must supply C0, HESS2, and $HESS2 input corresponding to the reference stationary point. (default=.FALSE.) C0 = an array of the coordinates of the stationary reference point (the coordinates in $DATA might well be some other coordinates). Give in the order x1,y1,z1,x2,y2,... in Angstroms. *** The next options apply to input choices which may read a $HESS at the initial DRC point, namely NNM or VIBLVL, or to those that read a $HESS2 at some reference geometry (NMANAL). HESS HESS2 = MIN indicates the hessian supplied for the initial geometry corresponds to a minimum (default). = TS indicates the hessian is for a saddle point. = MIN (default) or TS, the same meaning, for the reference geometry. These are used to decide if modes 1-6 (minimum) or modes 2-7 (TS) are to be excluded from the hessian as the translational and rotational contaminants. If the initial and reference geometries are the same, these two hessians will be duplicates of each other. 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


Input Description

$DRC

2-160

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 TRAJECT file. (default = 1) NPRT = 1 0 -1 2 1 0 -1 Print orbitals at all DRC points Print orbitals at F+L (default) Never print orbitals Punch all orbitals at all DRC points Punch all orbitals at F+L, and occupied orbitals at DRC points between Punch all orbitals at F+L only (default) Never punch orbitals

NPUN

=

==========================================================


Input Description

$MEX

2-161

===========================================================

$MEX

group

(relevant if RUNTYP=MEX)

This group governs a search for the lowest energy on the 3N-7 dimensional "seam" of intersection of two different electronic potential energy surfaces. Such Minimum Energy Crossing Points are important for processes such as spinorbit coupling that involve transfer from one surface to another, and thus are analogous to transition states on a single surface. The present program requires that the two surfaces differ in spin quantum number, or space symmetry, or both. Analytic gradients are used in the search. In case the two potential surfaces have identical spin and space symmetry, this kind of intersection point is referred to as a Conical Intersection. See $CONICL using RUNTYP=CONICAL instead. SCF1, SCF2 = define the molecular wavefunction types, possibly in conjunction with the usual MPLEVL and DFTTYP keywords.

MULT1, MULT2 = give the spin multiplicity of the states. Permissible combinations of wavefunctions are RHF with ROHF/UHF ROHF with ROHF UHF with UHF as well as their MP2 and DFT counterparts, and GVB with ROHF/UHF MCSCF with MCSCF (CISTEP=ALDET or GUGA only) NSTEP STPSZ = maximum number of search steps (default=50) = Step size during the search (default = 0.1D+00)

NRDMOS = = = = =

Initial orbitals can be read in 0 No initial orbitals (default) 1 Read in orbitals for first state (in $VEC1) 2 Read in orbitals for second state (in $VEC2) 3 Read in orbitals for both ($VEC1 and $VEC2)


Input Description NMOS1 NMOS2 NPRT

$MEX

2-162

= Number of orbitals for first state's $VEC1. = Number of orbitals for second state's $VEC2. = Printing orbitals = 0 No orbital printed out except at the first geometry (default) = 1 Orbitals are printed each geometry. If MCSCF is used, CI expansions are also printed.

Finer control of the convergence criterion: TDE = energy difference between two states (default = 1.0D-05)

TDXMAX = maximum displacement of coordinates (default = 2.0D-03) TDXRMS = root mean square displacement (default = 1.5D-03) TGMAX TGRMS = maximum of effective gradient between the two states (default = 5.0D-04) = root mean square effective gradient tolerance (default = 3.0D-04)

=========================================================== Usage notes: 1. Normally $CONTRL will not give SCFTYP or MULT keywords. SCF1 and SCF2 can be given in any order. The combinations permitted ensure roughly equal sophistication in the treatment of electron correlation. 2. After reading $MEX, SCFTYP and MULT will be set to the more complex of the two choices, which is considered to be RHF < ROHF < UHF < GVB < MCSCF. This permits the $SCF input defining a GVB wavefunction to be read and tested for correctness, in a GVB+ROHF run. Since only one SCFTYP is stored while reading the input, you might need to provide some keywords that are normally set by default for the other (such as ensuring DIIS is selected in $SCF if either of the states is UHF).


Input Description

$MEX

2-163

3. It is safest by far to prepare and read $VEC1 and $VEC2 groups so that you know what electronic states you start with. It is a good idea to regenerate both states at the end of the MEX search, to be sure that they remain as you began. 4. It is your responsibility to make sure that the states have a different space symmetry, or a different spin symmetry (or both). That is why note 3 is so important. 5. $GRAD1 and/or $GRAD2 groups containing gradients may be given to speed up the first geometry of the MEX search. 6. The search is even trickier than a saddle point search, for it involves the peaks and valleys of BOTH surfaces being generated. Starting geometries may be guessed as lying between the minima of the two surfaces, but the lowest energy on the crossing seam may turn out to be somewhere else. Be prepared to restart! 7. The procedure is a Newton-Raphson search, conducted in Cartesian coordinates, with a Lagrange multiplier imposing the constraint of equal energy upon the two states. The hessian matrices in the search are guessed at, and subjected to BFGS updates. Internal coordinates will be printed (for monitoring purposes) if you define $ZMAT, but the stepper operates in Cartesian coordinates only. No geometry constraints can be applied, apart from the point group in $DATA. A good paper to read about this kind of search is A.Farazdel, M.Dupuis J.Comput.Chem. 12, 276-282(1991)


Input Description

$CONICL

2-164

===========================================================

$CONICL

group

(relevant if RUNTYP=CONICAL)

This group governs a search for the lowest energy on the 3N-7 dimensional "seam" of intersection of two electronic potential energy surfaces of the same spin and space symmetry. Such Conical Intersections (CI) are important in photochemistry, where they serve as "funnels" for the transfer from an excited state to a lower state. See RUNTYP=MEX and the $MEX input for the simpler case where the two surfaces differ by either space or spin symmetry. Three search procedures are given, one of which requires the non-adiabatic coupling matrix element (NACME), and two others which do not require NACME information. The conical intersection search is available only for MCSCF (for which NACME are available) or for TD-DFT potential surfaces (where NACME are not available). The TD-DFT must be used in the Tamm/Dancoff approximation (see TAMMD in $TDDFT), but can be either conventional or spin-flip. The search utilizes some of the options of $STATPT, but note that the Schlegel stepper and HESS=CALC are not permitted. It may be reasonable to try the RFO stepper sometimes. The search can only be run in Cartesian coordinates. Restarts are possible only by updating the coordinates in $DATA. At present, the only solvation model that is supported is conventional TD-DFT with EFP1. OPTTYP = search procedure choice, see references below! = GPWNAC Gradient Projection with NACME, so this is only available for MCSCF. = BPUPD branching plane updating method (default) = PENALTY penalty-constrained optimization method Note that for MCSCF surfaces, if state-averaging is used, the program executes the code needed to produce NACME vectors, to producing the state-averaged gradients. There is essentially no extra time required to produce also the NACME, hence the GPWNAC stepper might as well be used.


Input Description

$CONICL

2-165

IXROOT = array of two states whose CI point is sought. For example, this might be IXROOT(1)=2,3 The roots are counted exactly the same as IROOT in the $DET or $TDDFT input groups. For the latter case, set IXROOT to 0 if you want the ground state to be one of the two surfaces searched on. There is no default for IXROOT! SYMOFF = flag to switch off point group symmetry, the default is .TRUE. DEBUG = flag to print debugging info, default is .FALSE.

The following are meaningful only for OPTTYP=PENALTY: TOLSTP = energy difference tolerance default=1d-6 Hartree TOLGRD = gradient convergence tolerance default=5d-3 Hartree/Bohr ALPHA SIGMA = parameter ensuring a singularity free penalty, default=0.02 Hartree = Lagrange multiplier for the penalty term. In case the energy gap between the states is not acceptable at the CI point, increase the value. default = 3.5 (unitless)

An understanding of the search procedures can be gained by reading the following papers: Gradient Projection with NACME: M.J.Bearpark, M.A.Robb, H.B.Schlegel Chem.Phys.Lett. 223, 269(1994) Branching Plane Updating method: S.Maeda, K.Ohno, K.Morokuma J.Chem.Theor Comput. 6, 1538(2010) Penalty constrained update method: B.G.Levine, C.Ko, J.Quenneville, T.J.Martinez Mol.Phys. 104, 1039(2006) B.G.Levine, J.D.Coe, T.J.Martinez J.Phys.Chem.B 112, 405(2008) A comparative study of the first two procedures is


Input Description

$CONICL

2-166

T.W.Keal, A.Koslowski, W.Thiel Theoret.Chem.Acc. 118, 837(2007) ===========================================================


Input Description

$MD

2-167

===========================================================

$MD

group

(relevant if RUNTYP=MD)

This group controls the molecular dynamics trajectory for a collection of quantum mechanical atoms and/or Effective Fragment Potential particles. A typical MD simulation starts with an equilibration phase, running long enough to produce a randomized structure and velocity distribution. Typically equilibration is done with an NVT ensemble, allowing the system to equilibrate to a desired temperature. A production run restarts with the positions and the velocity and quaternion data from the equilibration run, might use either a NVE or NVT ensemble, and collects radial distribution functions and other properties. Only a few properties are computed from the MD trajectory, apart from correct radial distribution functions. In particular, the pressures, diffusion constants, and heats of vaporization that appear on the printout (presently only for pure EFP runs) are from a preliminary code, which has not yet been verified. If the system contains only EFP particles, it may be placed in a periodic box, according to the minimum image convention. The optional periodic boundary conditions, along with cut-offs, are given in the $EFRAG input. See also the $EWALD input group for long-range electrostatic treatment if PBC is used. The first keywords relate to the steps: MDINT = MD integrator selection. = FROG (leapfrog). This is less accurate, and lacks the special ensemble stepper option NVTNH. = VVERLET (velocity Verlet) - default. = MD time step size, in seconds, default=1.0d-15, which is a femtosecond. selects a integrator step appropriate to the desired ensemble. This is only implemented for

DT NVTNH


Input Description

$MD

2-168

velocity Verlet. = 0 means use NVE Verlet stepping = 1 means use NVT Verlet stepping = 2 means use Nose/Hoover chain NVT Verlet stepping The default is 2 if either NVT option RSTEMP or RSRAND is chosen, but is 0 otherwise. NSTEPS = number of MD time steps to be found in this run, default=10000. TTOTAL = total time elapsed in the previous part of a MD trajectory which is being restarted (READ=.TRUE.). The default means this trajectory is a new one, or perhaps the start of a production phase of the MD. (default=0.0 seconds) *** BATHT = bath temperature, in Kelvin (default=300.0) This value is used during NVT runs, or if the MD is initialized to a Maxwell-Boltzmann velocity distribution. *** Two options exist to create NVT runs, to bring the system to a desired bath temperature. If neither is chosen, the ensemble is NVE: RSTEMP = flag to rescale the temperature. DTEMP default=.FALSE.

= temperature range for the RSTEMP option. The velocities are rescaled to the bath temperature if T < (BATHT-DTEMP) or T > (BATHT+DTEMP). The default is DTEMP=100.0 degrees.

RSRAND = flag to reset to Maxwell-Boltzmann distribution, using random numbers (same algorithm as MBT and MBR) to choose individual velocity magnitudes and directions. default=.FALSE. NRAND = number of steps for the RSRAND option. Reassign velocities (translational and rotational) every NRAND time steps. Default=1000.


Input Description

$MD

2-169

NVTOFF = step number at which to turn off either NVT thermostat, and switch to NVE. At this point, the NVTNH parameter will be reset to 0, and the PROD flag will be turned on, so that the production run will start (gathering and printing the RDF information to .log file). This keyword is also useful in NVE runs to postpone the accumulation of production information. The default means no switch to NVE (default=0). JEVERY = report simulation quantities (write info such as energies, temps, etc. to .log file) and collect RDF info each JEVERY time step. Default=10 KEVERY = write coordinates (to log and TRAJECT files), velocity/quaternion restart info (to the TRAJECT file and RDFs (to log file) at each KEVERY step. default=100 PROD = production run, at present this means only that information for radial distribution functions is collected, and printed. default=.FALSE. = spacing for radial bins in RDF calculations, default=0.02 Angstroms. = step number at which to begin collecting data for the other properties, such as pressure and diffusion constants. This should be a value between 1 and NSTEPS, as it counts off the current run's steps. Default=0.

DELR NPROP

PBCOUT = print PBC coordinates in the end of simulation (i.e. all molecules will be contained in one box) Default=.FALSE. *** SSBP = flag to add spherical solvent boundary condition using the harmonic restraint potential(V) V=0.5*SFORCE*(R-DROFF)**2. (Default=.FALSE)

SFORCE = the force constant for SSBP in kcal/mol-A**2 (default: 0.0).


Input Description

$MD

2-170

CCMS

= flag to add a harmonic potential to constrain the center of mass of the QM subsystem. This will keep the QM subsystem in QM/MM-MD with spherical boundary conditions near the center of sphere. (Default=.FALSE)

CFORCE = the force constant for CCMS in kcal/mol-A**2 (default: 0.0). DROFF = is an array of distances such that V=0 if R
USAMP

UFORCE = the force constant for USAMP in kcal/mol-A**2 (default: 0.0). RZERO = array of reference values for USAMP. Distances and angles are in Angstrom and Degree, respectively. (default: all 0.0). specifies =0 =1 =2 =3 =4 the type of USAMP bond constraint (default value) asymmetric bond constraint angle constraint torsion constraint asymmetric and bond constraint. This is for 2 dimensional umbrella sampling of an atom shift reaction (atom 1 is moving from atom 2 to 3), and a bond distance between atom 4 and 5. So it requires 5 elements in IPAIR! = 5 asymmetric coordinate + Restraint IPAIR contains three atoms, where atom 1 is moving from atom 2 to 3 by asymmetric coordinate. This restraint forces atom 1 to stay between atom 2 and 3. = 6 Normal constraining for up to 2 bonds based on flexible constraints. IPAIR must contain 4 values for the two bonds.

IUSTYP


Input Description

$MD

2-171

= 7 Asymmtric coordinate based on flexible constraints. IPAIR should contains all the corresponding atoms. = 8 The center of mass coordinates of two QM parts. IPAIR should include only the first part of QM atoms. = 9 The center of mass coordinates with respect to the center of sphere (0,0,0). IPAIR should not be used, since GAMESS will automatically generate the center of mass of your QM part. IPAIR is an array of atom numbers to which the umbrella sampling potentials are applied. For example, a bond constraint for atom 1 and 2 is IPAIR(1)=1,2. An asymmetric bond constraint requires three atoms, in which the atom 1 migrates from atom 2 to 3. IUSTYP=4 requires 5 atoms: the first three atoms are for asymmetric bond constraint. Default: 0. *** MREMD specifies the type of REMD (Replica Exchange MD). NGROUPS in $GDDI must be greater than 1 in order for this option to be effective. Note that currently REMD works only for NVT. = 0 no REMD = 1 one dimension REMD Temperature exchange at LEVERY steps during REMD runs. (default: 1000) If JEVERY.gt.LEVERY, JEVERY is reset to LEVERY. *** The following keywords control starting MD conditions. Normally an MD trajectory is initiated with both MBT and MBR chosen, while restarts would select only READ. The restart data is written to the TRAJECT file. To restart requires merging particle coordinates into $DATA and/or $EFRAG, and placing the new $MD input below your existing $MD input, thus keeping your choices for the variables above (both $MD input groups will be read).

LEVERY


Input Description

$MD

2-172

MBT MBR QRAND

= get translational velocities from a random Maxwell-Boltzmann ensemble. Default=.FALSE. = get rotational velocities from a random MaxwellBoltzmann ensemble. Default=.FALSE. = if .TRUE., generate random quaternions, an option that is not normally chosen. if .FALSE., use EFP particle coordinates and the initial MBT/MBR assigned velocities to set correct quaternion data (default is .FALSE.) = read velocities (translational and rotational) and quaternions and their first and second derivatives from input file. Default is .FALSE. Set the other three values MBT/MBR/QRAND off if you choose restarting with READ.

READ

For READ=.TRUE., the following restart data is required. This data may be copied from the TRAJECT file, in exactly the format it was written out. The required data depends on your choice for the integrator, see MDINT above. In addition, you will need to update the particle coordinates in $DATA and/or $EFRAG, using data from the TRAJECT file. TVELQM(1)= quantum atom's translational velocities (both). TVEL(1)= array of EFP translational velocities (both). RVEL(1)= array of EFP rotational velocities (VVERLET). RMOM(1)= array of EFP rotational momenta (FROG). QUAT(1)= array of EFP quaternions (both). QUAT1D(1)= EFP quaternion first derivatives (VVERLET). QUAT2D(1)= EFP quaternion second derivatives (VVERLET). extra reading: "Computer Simulation of Liquids" M.P.Allen, D.J.Tildesley Oxford Science, 1987 "Understanding Molecular Simulation" D.Frenkel, B.Smit Academic Press, 2002 ===========================================================


Input Description

$RDF

2-173

==========================================================

$RDF

group

(relevant for RUNTYP=MD)

This group defines the pairs of atoms for which the radial distribution functions are to be computed, at the end of a molecular dynamics trajectory. The input is similar in style to $EFRAG, consisting of separate lines, with the word STOP ending each particular pair. Line 1. NRDF= gives the number of RDFs which should be computed. Line 2. gives a string for the printout (a good choice involves both atoms, such as ClCl), the name of the $FRAGNAME containing the first atom of the pair, the name of the $FRAGNAME input group with the second atom of the pair, and how many such pairs exist. Line 3.