(15 Jul 2024)
*********************************
* *
* 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.
*
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
$RICC resolution-of-the-identity CC RICCRHF:DRRICC_RHF
$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
$REKS REKS input REKS :REXINPUT
Note that MCSCF and CI input is listed below.
Potential energy surface options:
$STATPT geometry search control STATPT:SETSIG
$TRUDGE nongradient optimization TRUDGE:TRUINP
$TRURST restart data for TRUDGE TRUDGE:TRUDGX
$FORCE hessian, normal coordinates HESS :HESSX
$CPHF coupled-Hartree-Fock options CPHF :CPINP
$CPMCHF coupled-MR-Hartree-Fock options MCPCGX:MCPCGX
$MASS isotope selection VIBANL:RAMS
$HESS force constant matrix (formatted) HESS :FCMIN
$GRAD gradient vector (formatted) HESS :EGIN
$DIPDR dipole deriv. matrix (formatted) HESS :DDMIN
$ALPDR alpha polar. der. (formatted) RAMAN :ADMIN
$VIB HESSIAN restart data (formatted) HESS :HSSNUM
$VIB2 num GRAD/HESS restart (formatted) HESS :HSSFUL
$VSCF vibrational anharmonicity VSCF :VSCFIN
$VIBSCF VSCF restart data (formatted) VSCF :VGRID
$GAMMA 3rd nuclear derivatives HESS :GAMMXX
$EQGEOM equilibrium geometry data HESS :FFCARX
$HLOWT hessian data from equilibrium HESS :FFCARX
$GLOWT 3rd derivatives at equilibrium HESS :FFCARX
$IRC intrinsic reaction coordinate RXNCRD:IRCX
$DRC dynamic reaction path DRC :DRCDRV
$MEX minimum energy crossing point MEXING:MEXINP
$CONICL conical intersection search
$MD molecular dynamics trajectory MDEFP :MDX
$RDF radial dist. functions for MD MDEFP :RDFX
$GLOBOP Monte Carlo global optimization GLOBOP:GLOPDR
$GLBFRG Monte Carlo atom groups GLOBOP:GLOPDR
$GRADEX gradient extremal path GRADEX:GRXSET
$SURF potential surface scan SURF :SRFINP
$NEB nudged elastic band NEBPATH:NEBPX
Interpretation, properties:
$LOCAL localized molecular orbitals LOCAL :LMOINP
$TRUNCN localized orbital truncations EFPCOV:TRNCIN
$ELMOM electrostatic moments PRPLIB:INPELM
$ELPOT electrostatic potential PRPLIB:INPELP
$ELDENS electron density PRPLIB:INPELD
$ELFLDG electric field/gradient PRPLIB:INPELF
$POINTS property calculation points PRPLIB:INPPGS
$GRID property calculation mesh PRPLIB:INPPGS
$PDC MEP fitting mesh PRPLIB:INPPDC
$MGC mean gradient charges
$RADIAL atomic orbital radial data PRPPOP:RADWFN
$MOLGRF orbital plots PARLEY:PLTMEM
$STONE distributed multipole analysis PRPPOP:STNRD
$COMP thermochemical calculation COMP :COMPX
$RAMAN Raman intensity RAMAN :RAMANX
$NMR NMR shielding tensors NMR :NMRX
$MOROKM Morokuma energy decomposition MOROKM:MOROIN
$LMOEDA LMO-based energy decomposition MOROKM:MMOEDIN
$QMEFP QM/EFP energy decomposition EFINP :QMEFPAX
$FFCALC finite field polarizabilities FFIELD:FFLDX
$TDHF time dependent HF of NLO props TDHF :TDHFX
$TDHFX TDHF for NLO, Raman, hyperRaman TDX:FINDTDHFX
Solvation models:
$EFRAG use effective fragment potential EFINP :EFINP
$EFPARM additional parameters for EFP EFINP :EFPARMINP
$FRAGNAME specifically named fragment pot. EFINP :RDSTFR
$FRGRPL inter-fragment repulsion EFINP :RDDFRL
$EWALD Ewald sums for EFP electrostatics EWALD :EWALDX
$MAKEFP generate effective fragment pot. EFINP :EFPX
$PRTEFP simplified EFP generation EFINP :PREFIN
$DAMP EFP multipole screening fit CHGPEN:CGPINP
$DAMPGS initial guess screening params CHGPEN:CGPINP
$PCM polarizable continuum model PCM :PCMINP
$PCMGRD PCM gradient control PCMCV2:PCMGIN
$PCMCAV PCM cavity generation PCM :MAKCAV
$TESCAV PCM cavity tesselation PCMCV2:TESIN
$REORG solvent reorganization in IEF-PCM REORG :RORGIN
$NEWCAV PCM escaped charge cavity PCM :DISREP
$IEFPCM PCM integral equation form. data PCM :IEFDAT
$PCMITR PCM iterative IEF input PCMIEF:ITIEFIN
$DISBS PCM dispersion basis set PCMDIS:ENLBS
$DISREP PCM dispersion/repulsion PCMVCH:MORETS
$SVP Surface Volume Polarization model SVPINP:SVPINP
$SVPIRF reaction field points (formatted) SVPINP:SVPIRF
$COSGMS conductor-like screening model COSMO :COSMIN
$SCRF self consistent reaction field SCRF :ZRFINP
Integral, and integral modification options:
$ECP effective core potentials ECPLIB:ECPPAR
$MCP model core potentials MCPINP:MMPRED
$RELWFN scalar relativistic integrals INPUTB:RWFINP
$EFIELD external electric field PRPLIB:INPEF
$INTGRL 2e- integrals INT2A :INTIN
$FMM fast multipole method QMFM :QFMMIN
$TRANS integral transformation TRANS :TRFIN
Fragment Molecular Orbital method:
$FMO define FMO fragments FMOIO :FMOMIN
$FMOPRP FMO properties and convergers FMOIO :FMOPIN
$FMOXYZ atomic coordinates for FMO FMOIO :FMOXYZ
$AFOMOD capping atom input for FMO FMOLIB:MAKELMO
$OPTFMO input for special FMO optimizer FMOGRD:OPTFMO
$FMOHYB localized MO for FMO boundaries FMOIO :FMOLMO
$FMOBND FMO bond cleavage definition FMOIO :FMOBON
$FMOENM monomer energies for FMO restart FMOIO :EMINOU
$FMOEND dimer energies for FMO restart FMOIO :EDIN
$OPTRST OPTFMO restart data FMOGRD:RSTOPT
$GDDI group DDI definition INPUTA:GDDINP
Polymer model:
$ELG polymer elongation method ELGLIB:ELGINP
Divide and conquer model:
$DANDC DC SCF input DCLIB :DCINP
$DCCORR DC correlation method input DCLIB :DCCRIN
$SUBSCF subsystem definition for SCF DCLIB :DFLCST
$SUBCOR subsystem definition for MP2/CC DCLIB :DFLCST
$MP2RES restart data for DC-MP2 DCMP2 :RDMPDC
$CCRES restart data for DC-CC DCCC :RDCCDC
clusters in molecules
$CIMINP controls clusters in molecules CIMINF:CIMINP
$CIMATM fine tune calculation level CIMINF:CIMINP
$CIMFRG fine tune atomic fragmentation CIMINF:CIMPRT
Quantum mechanics/molecular mechanics model:
$QUANPO QuanPol calculation QUANPO:QUANPOL
$FFDATA QuanPol coordinates for molecules QUANPO:QUANPOL
$FFDATB QuanPol coordinates for molecules QUANPO:QUANPOL
$FFPDB QuanPol coordinates for proteins QUANPO:QUANPOL
MCSCF and CI wavefunctions, and their properties:
$CIINP control over CI calculation GAMESS:WFNCI
$DET determinant full CI for MCSCF ALDECI:DETINP
$CIDET determinant full CI ALDECI:DETINP
$GEN determinant general CI for MCSCF ALGNCI:GCIINP
$CIGEN determinant general CI ALGNCI:GCIINP
$ORMAS determinant multiple active space ORMAS :FCINPT
$CEEIS CI energy extrapolation CEEIS :CEEISIN
$CEDATA restart data for CEEIS CEEIS :RDCEEIS
$GCILST general MCSCF/CI determinant list ALGNCI:GCIGEN
$GMCPT general MCSCF/CI determinant list GMCPT :OSRDDAT
$PDET parent determinant list GMCPT :OSMKREF
$ADDDET add determinants to reference GMCPT :OSMKREF
$REMDET remove determinants from ref. GMCPT :OSMKREF
$SODET determinant second order CI FSODCI:SOCINP
$DRT GUGA distinct row table for MCSCF GUGDRT:ORDORB
$CIDRT GUGA CI (CSF) distinct row table GUGDRT:ORDORB
$MCSCF control over MCSCF calculation MCSCF :MCSCF
$MRMP MRPT selection MP2 :MRMPIN
$DETPT det. multireference pert. theory DEMRPT:DMRINP
$MCQDPT CSF multireference pert. theory MCQDPT:MQREAD
$EXCORR interface to MPQC's R12 program EXCORR:GETEXC
$CASCI IVO-CASCI input IVOCAS:IVODRV
$IVOORB fine tuning of IVO-CASCI IVOCAS:ORBREAD
$CISORT GUGA CI integral sorting GUGSRT:GUGSRT
$GUGEM GUGA CI Hamiltonian matrix GUGEM :GUGAEM
$GUGDIA GUGA CI diagonalization GUGDGA:GUGADG
$GUGDM GUGA CI 1e- density matrix GUGDM :GUGADM
$GUGDM2 GUGA CI 2e- density matrix GUGDM2:GUG2DM
$LAGRAN GUGA CI Lagrangian LAGRAN:CILGRN
$TRFDM2 GUGA CI 2e- density backtransform TRFDM2:TRF2DM
$DIABAT diabatic states DIAB:DIABINP
$TRANST transition moments, spin-orbit TRNSTN:TRNSTX
LibXC and related sections:
$LIBXC control LibXC variables LIBXC :LIBXC_INPUT
$LDA_K LDA kinetic functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$LDA_X LDA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$LDA_C LDA correlation functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$LDA_XC LDA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$HLDA_XC hybrid LDA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$GGA_K control GGA kinetic functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$GGA_X GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$GGA_C GGA correlation functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$GGA_XC GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$HGGA_X hybrid GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$HGGA_XC hybrid GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$MGGA_K m-GGA kinetic functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$MGGA_X m-GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$MGGA_C m-GGA correlation functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$MGGA_XC m-GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$HMGGA_XC hybrid m-GGA XC functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
$HMGGA_X hybrid m-GGA exchange functionals LIBXC :LIBXC_CHOOSE_FUNCTIONAL
RISM and related sections:
$RISMUN solvent effect using RISM RSM_UNINP :RSM_UNCN
$FITINP fitting of cSED RSM_SDDUTL:RSM_FITINP
$ABSINP specifies the weight in the fitting RSM_SDDUTL:RSM_ABSIN
OpenMP offload and accelerator sections:
$DEVICE Accelerator device settings RIMP2GPU:DEVICEINP
$OFFLOAD OpenMP target offload controls MOD_OFFLOAD:MOD_OFFLOAD
* this column is more useful to programmers than to users.
==========================================================
$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 specifies the self-consistent field
wavefunction. You may choose from
= RHF Restricted Hartree Fock calculation
(default)
= UHF Unrestricted Hartree Fock calculation
= ROHF Restricted open shell Hartree-Fock.
(high spin, see GVB for low spin)
= GVB Generalized valence bond wavefunction,
or low spin ROHF. (needs $SCF input)
= MCSCF Multiconfigurational SCF wavefunction
(this requires $DET or $DRT input)
= REKS Spin-restricted ensemble referenced
Kohn-Sham method.
(this requires $REKS input)
= NONE 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, or if CITYP=SFDET). Scalar relativistic
effects may be incorporated using RELWFN for any of these
wavefunction choices, correlated or not.
DFTTYP = NONE ab initio computation (default)
= XXXXXX perform density functional theory run,
using the functional specified. Many
choices for XXXXXX are listed in the
$DFT and $TDDFT input groups.
TDDFT = NONE no excited states (default)
= EXCITE generate time-dependent DFT excitation
energies, using the DFTTYP= functional,
for RHF or UHF references. Analytic
nuclear gradients are available for RHF.
See $TDDFT.
= SPNFLP spin-flip TD-DFT, for either UHF or ROHF
references. Nuclear gradients and
solvent effects are programmed. See $TDDFT.
= MRSF Mixed-Reference Spin-Flip (MRSF) TD-DFT
for ROHF reference. Nuclear gradients and
nonadiabatic MD are programmed for ROHF.
See $TDDFT. CUHF option is also utilized
to improve the convergence of MRSF.
= POL (hyper)polarizability calculation, for
RHF only. See $TDDFT.
* * * * *
PDFTYP = NONE ab initio computation (default)
= XXXXXX perform multiconfiguration pair-density
functional theory run, using the
functional specified. Many
choices for XXXXXX are TPBE, FTPBE,
TBLYP,FTBLYP,TREVPBE and FTREVPBE.
SCFTYP=MCSCF is required.
* * * * *
VBTYP = NONE no valence bond calculation (default)
= VB2000 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 = chooses Moller-Plesset perturbation
theory level, after the SCF. See $MP2,
or $MRMP for MCSCF.
= 0 skip the MP computation (default)
= 2 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.
MP2 can also be used on top of CITYP=SFDET to get the
SF-MPRMP2 energy correction.
* * * * *
CITYP = chooses CI computation after the SCF,
for any SCFTYP except UHF.
= NONE skips the CI. (default)
= CIS single excitations from a SCFTYP=RHF
reference, only. This is for excited
states, with analytic nuclear gradients
available. See the $CIS input group.
= SFCIS spin-flip style CIS, see $CIS input.
= ALDET runs the Ames Laboratory determinant
full CI package, requiring $CIDET.
= ORMAS runs an Occupation Restricted Multiple
Active Space determinant CI. The input
is $CIDET and $ORMAS.
= FSOCI runs a full second order CI using
determinants, see $CIDET and $SODET.
= GENCI runs a determinant CI program that
permits arbitrary specification of
the determinants, requiring $CIGEN.
= GUGA 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.
= SFORM performs SF-ORMAS calculation. SCFTYP
must be ROHF, and the MULSF must be
specified in $CIDET. The active space
may be defined as usual in $ORMAS.
Analytic gradients are available.
= SFDET performs SF-CI calculation using the
ALDET code. SCFTYP must be ROHF and
MULSF must be specified. This is mostly
useful when one wants to perform
an SF-MRMP2 calculation.
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 ground-state coupled-cluster (CC) or,
optionally, equation-of-motion coupled-cluster
(EOMCC) calculation for excited, electron-attached,
and ionized electronic states using either the RHF
or ROHF reference. See also $CCINP and $EOMINP.
Only CCSD and CCSD(T) for RHF can run in parallel
using DDI. All CC and EOMCC methods may run in
parallel on a single SMP node if GAMESS is compiled
with a parallel math library, such as MKL.
Only CCSD and CCSD(T) calculations can run using the
RI approximation for two-electron repulsion integrals
and a hybrid MPI/OpenMP parallelization model.
The ground-state CC options:
= NONE skips CC computation (default).
= LCCD perform a CC calculation
using the linearized coupled-cluster
method with double excitations.
= CCD perform a CC calculation using the
coupled-cluster method with doubles.
= CCSD perform a CC calculation with
single and double excitations.
= CCSD(T) in addition to CCSD, the non-iterative
triples corrections are computed, giving
standard CCSD[T] and CCSD(T) energies.
= R-CC in addition to all CCSD(T) calculations,
the renormalized R-CCSD[T] and
R-CCSD(T) energies are computed.
= CR-CC in addition to all R-CC calculations,
the completely renormalized CR-CCSD[T]
and CR-CCSD(T) energies are computed.
= CR-CCL in addition to a CCSD ground state and
the associated left CCSD calculations, the
non-iterative triples energy correction
defining the rigorously size-extensive
completely renormalized CR-CC(2,3) approach
(sometimes also called CR-CCSD(T)_L)
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 Section 4 of this manual.
Generally, CR-CC(2,3),D
provides the most accurate energetics.
= CCSD(TQ) in addition to all R-CC calculations,
non-iterative triples and quadruples
corrections are determined, to give CCSD(TQ)
and various R-CCSD(TQ) energies.
= CR-CC(Q) in addition to all CR-CC and CCSD(TQ)
calculations, the CR-CCSD(TQ) energies
are obtained.
= CCSD3A perform CC calculations with singles,
doubles, and active-space triples,
designated as CCSDt. When all orbitals
in the MO basis are set to be active,
CCSDt becomes full CCSDT.
= CCT3 in addition to CCSDt, compute the non-iterative
CC(P;Q)-based CC(t;3) energy correction due to
the missing triples not captured by CCSDt.
As for CR-CCL, variant D of CC(t;3) is
generally most accurate and closest to CCSDT.
= ACCSD perform ACCSD calculations, i.e., CCSD-type
calculations that select or scale T2^2 diagrams
in the equations projected on doubly excited
determinants in the spirit of approximate
coupled-pair (ACP) approaches.
= ACC23 in addition to ACCSD calculations, compute the
CR-CC(2,3)-type non-iterative triples correction
to the ACCSD energy.
= ACCSD3A perform ACCSDt calculations, i.e., CCSDt-type
calculations that select or scale T2^2 diagrams
in the equations projected on doubly excited
determinants in the spirit of ACP approaches.
= ACCT3 in addition to ACCSDt calculations, compute the
CC(t;3)-type non-iterative energy correction
due to the missing triples not captured by
ACCSDt.
ACCSD, ACC23, ACCSD3A, and ACCT3 options require the
ACP array with T2^2 diagram scaling factors other
than 1.0 to be provided in $CCINP.
CCSD, CR-CCL, CCSD3A, CCT3, ACCSD, ACC23, ACCSD3A,
and ACCT3 options work for both RHF and ROHF references;
the remaining ground-state CC calculations can only be
performed using SCFTYP=RHF.
Excited-state EOMCC options:
= EOM-CCSD in addition to a CCSD ground state,
excited states are calculated using the
equation-of-motion coupled-cluster
method with singles and doubles (EOMCCSD).
= CR-EOM in addition to CCSD and EOMCCSD,
non-iterative triples corrections to CCSD
ground-state and EOMCCSD excited-state
energies are determined using completely
renormalized CR-EOMCCSD(T) approaches.
In typical applications, delta-CR-EOMCCSD(T)
triples corrections to EOMCCSD excitation
energies are most useful.
= CR-EOML in addition to printing all results that
CR-EOM obtains and solving the left
eigenstate EOMCCSD problem used for
properties other than energy, the CR-EOMCC(2,3)
triples corrections analogous to those of the
ground-state CR-CCL approach are calculated.
In typical applications, delta-CR-EOMCC(2,3)
triples corrections to EOMCCSD excitation
energies are most useful.
CCTYP=EOM-CCSD is available for RHF or ROHF references;
the CR-EOM and CR-EOML triples corrections can only
be calculated using SCFTYP=RHF.
Single and double electron-attachment and ionization
potential EOMCC options for open-shell systems with
one or two electrons or holes outside closed-shell
cores (available for RHF or ROHF references,
although the EA- and IP-EOMCC calculations using
ROHF reference are discouraged):
= EA-EOM2 electron-attachment EOMCC calculations
with up to 2p1h excitations treated
fully (i.e., EA-EOMCCSD).
= EA-EOM3A electron-attachment EOMCC calculations
with all 1p and 2p1h and active-space
3p2h excitations [i.e., EA-EOMCCSDt or
EA-EOMCC(3p2h){Nu}].
= IP-EOM2 ionization potential EOMCC calculations
with up to 2h1p excitations treated fully
(i.e., IP-EOMCCSD).
= IP-EOM3A ionization potential EOMCC calculations
with all 1h and 2h1p and active-space
3h2p excitations [i.e., IP-EOMCCSDt or
IP-EOMCC(3h2p){No}].
= DEAEOM2 double electron-attachment EOMCC with up to
3p1h excitations treated fully [i.e.,
DEA-EOMCC(3p1h) or DEA-EOMCCSD].
= DEAEOM2A double electron-attachment EOMCC with all
2p and active-space 3p1h excitations
[i.e., DEA-EOMCC(3p1h){Nu}].
= DEAEOM3 double electron-attachment EOMCC with up to
4p2h excitations treated fully [i.e.,
DEA-EOMCC(4p2h); most accurate, but
very expensive!].
= DEAEOM3A double electron-attachment EOMCC with all 2p
and active-space 3p1h and 4p2h excitations
[i.e., DEA-EOMCC(3p1h,4p2h){Nu}].
= DEAEOM3B double electron-attachment EOMCC with all 2p
and 3p1h and active-space 4p2h excitations
[i.e., DEA-EOMCC(4p2h){Nu}].
= DIPEOM2 double ionization potential EOMCC with up to
3h1p excitations treated fully [i.e.,
DIP-EOMCC(3h1p) or DIP-EOMCCSD].
= DIPEOM3 double ionization potential EOMCC with up to
4h2p excitations treated fully [i.e.,
DIP-EOMCC(4h2p)].
= DIPEOM3A double ionization potential EOMCC with all 2h
and 3h1p and active-space 4h2p excitations
[i.e., DIP-EOMCC(4h2p){No}].
Labels "p" and "h" in the description of (D)EA- and (D)IP-EOMCC
methods refer to particles (correlated orbitals unoccupied
in the underlying reference) and holes (correlated orbitals
occupied in the underlying reference), respectively. Nu and No in
the description of (D)EA- and (D)IP-EOMCC methods indicate the numbers
of active unoccupied and active occupied orbitals used to select
the dominant 3p2h, 3h2p, 3p1h, 4p2h, and 4h2p contributions.
EA/IP-EOMCC runs produce ground and excited states of systems
obtained by attaching one electron to (EA) or removing one electron
from (IP) the underlying closed-shell cores described by CCSD.
DEA/DIP-EOMCC runs produce ground and excited states of systems
obtained by attaching two electrons to (DEA) or removing two electrons
from (DIP) the underlying closed-shell cores described by CCSD.
As all other EOMCC calculations, EA/IP- and DEA/DIP-EOMCC
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 in 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 CC energies are obtained
at similar cost to the standard values, while completely
renormalized ones cost about twice the time. For usage tips
and more information about resources on the various coupled
cluster methods, see Section 4, 'Further Information'.
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 Selects all-electron scalar relativity treatment.
See the $RELWFN input group for more information,
including nuclear derivative availability.
= NONE use the basic Schrodinger equation (default)
= LUT-IOTC local unitary transformation modification
of IOTC, due to H.Nakai, J.Seino, Y.Nakajima.
This is the fastest and most numerically
reliable scalar relativity method, so it is
preferred over RESC, DK, or IOTC.
= LUTIOTC2 local unitary transformation modification
of IOTC2E for one- and two-electron terms,
due to H. Nakai and J. Seino. This option
is not available to range separated functionals
in DFT.
= IOTC infinite-order two-component method of
M. Barysz and A.J. Sadlej.
= IOTC2E infinite-order two-component method including
two-electron term of J. Seino and M. Hada,
in addition to IOTC. This option is not available
to range separated functionals in DFT.
= 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 specifies the type of computation, for
example at a single geometry point:
= ENERGY Molecular energy. (default)
= GRADIENT Molecular energy plus gradient.
= HESSIAN Molecular energy plus gradient plus
second derivatives, including harmonic
harmonic vibrational analysis.
Analytic FMO Hessians and PAVE should use
RUNTYP=FMOHESS. Semi-numerical FMO
Hessians except for PAVE should use
RUNTYP=HESSIAN.
See the $FORCE and $CPHF input groups.
= FMOHESS used for FMO Hessians: analytic for DFTB,
RHF, DFT, UHF, and ROHF;
numerical runs also possible when using PAVE.
FMOHESS reads $FORCE as HESSIAN.
= GAMMA Evaluate up to 3rd nuclear derivatives,
by finite differencing of Hessians.
See $GAMMA, and also NFFLVL in $CONTRL.
multiple geometry options:
= OPTIMIZE Optimize the molecular geometry using
analytic energy gradients. See $STATPT.
= TRUDGE Non-gradient total energy minimization.
See $TRUDGE and $TRURST.
= SADPOINT Locate saddle point (transition state).
See $STATPT.
= MEX Locate minimum energy crossing point on
the intersection seam of two potential
energy surfaces. See $MEX.
= CONICAL Locate conical intersection point on
the intersection seam of two potential
energy surfaces. See $CONICL.
= IRC Follow intrinsic reaction coordinate.
See $IRC.
= VSCF anharmonic vibrational corrections.
See $VSCF.
= DRC Follow dynamic reaction coordinate.
See $DRC.
= MD molecular dynamics trajectory, see $MD.
= GLOBOP Monte Carlo-type global optimization.
See $GLOBOP.
= OPTFMO genuine FMO geometry optimization using
nearly analytic gradient. See $OPTFMO.
= GRADEXTR Trace gradient extremal. See $GRADEX.
= SURFACE Scan linear cross sections of the
potential energy surface. See $SURF.
single geometry property options:
= COMP composite thermochemistry calculation,
including G3MP2. See $COMP input.
= G3MP2 evaluate heat of formation using the
G3(MP2,CCSD(T)) methodology. See test
example exam43.inp for more information.
= PROP Molecular properties will be calculated.
Orbital localization can be requested as
well. See $ELPOT, etc.
Converged orbitals must be input in a
$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.
= RAMAN computes Raman intensities, see $RAMAN.
= NACME 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 NMR shielding tensors for closed shell
molecules by the GIAO method. See $NMR.
= EDA Perform energy decomposition analysis.
Give one of $MOROKM or $LMOEDA inputs.
= QMEFPEA QM/EFP solvent energy analysis,
see $QMEFP.
= TRANSITN Compute radiative transition moment or
spin-orbit coupling. See $TRANST.
= FFIELD applies finite electric fields, most
commonly to extract polarizabilities.
See $FFCALC.
= TDHF analytic computation of time dependent
polarizabilities. See $TDHF.
= TDHFX extended TDHF package, including nuclear
polarizability derivatives, and Raman
and Hyper-Raman spectra. See $TDHFX.
= MAKEFP creates an effective fragment potential,
for SCFTYP=RHF or ROHF only.
See $MAKEFP, $DAMP, $DAMPGS, $STONE, ...
= FMO0 performs the free state FMO calculation.
See $FMO.
= FCIDUMP Performs an ROHF calculation and generates
a FCIDUMP file.
= GRADNAC $GRADIENT + NonAdibatic Coupling Term.
NACT is designed in cases of using external
MD drivers and requires coordinates of
current and previous steps.
* * * * * * * * * * * * * * * * * * * * * * * * * * * * *
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
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.
FMO uses GDDI in the fragment-way, other
runs use GDDI for offloading for offsets.
The default is .FALSE.
EXETYP = RUN Actually do the run. (default)
= CHECK Wavefunction and energy will not be
evaluated. This lets you speedily
check input and memory requirements.
See the overview section for details.
Note that you must set PARALL=.TRUE.
in $SYSTEM to test distributed memory
allocations.
= EXPERT Declare your expert status to enjoy
special priveleges: some runs will be
allowed to proceed that normally stop
with the view of protecting novice users.
In particular, unnormalized basis sets
may be allowed.
Also some warnings may be suppressed.
= FMOHMO Normalize HMOs in $FMOHYB and quit.
= DEBUG Massive amounts of output are printed,
useful only if you hate trees.
= routine Maximum output is generated by the
routine named. Check the source for
the routines this applies to.
* * * * * * *
ICHARG = Molecular charge. (default=0, neutral)
MULT = Multiplicity of the electronic state
= 1 singlet (default)
= 2,3,... doublet, triplet, and so on.
ICHARG and MULT are used directly for RHF, UHF, ROHF.
For GVB, these are implicit in the $SCF input, while
for MCSCF or CI, these are implicit in $DRT/$CIDRT or
$DET/$CIDET input. You must still give them correctly.
* * * 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!!!
= ZMT GAUSSIAN style internals will be input.
= ZMTMPC MOPAC style internals will be input.
= FRAGONLY means no part of the system is treated
by ab initio means, hence $DATA is not
given. The system is 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 3-
fold or higher rotation axis, the degenerate moments of
inertia often cause problems choosing correct symmetry
unique axes, in which case you must use COORD=UNIQUE rather
than Z-matrices.
Warning: The reorientation into principal axes is done
only for atomic coordinates, and is not applied to the axis
dependent data in the following groups: $VEC, $HESS, $GRAD,
$DIPDR, $VIB, nor Cartesian coords of effective fragments
in $EFRAG. COORD=UNIQUE avoids reorientation, and thus is
the safest way to read these.
Note: 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
NZVAR = 0 Use Cartesian coordinates (default).
= M 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.
= M For other COORD choices, or if $ZMAT is given:
the internal coordinates will be those defined
in $ZMAT. This allows more sophisticated
internal coordinate choices. M is ordinarily
3N-6 (3N-5), unless $ZMAT has linear bends.
NZVAR refers mainly to the coordinates used by OPTIMIZE
or SADPOINT runs, but may also print the internal's
values for other run types. You can use internals to
define the molecule, but Cartesians during optimizations!
* * * * * * *
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 = pseudopotential selection.
= NONE all electron calculation (default).
= READ read ECP potentials in the $ECP input.
= SBKJC use Stevens, Basch, Krauss, Jasien,
Cundari ECP potentials for all heavy
atoms (Li-Rn are available).
= HW use Hay, Wadt ECP potentials for heavy
atoms (Na-Xe are available).
= MCP 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 = Spherical Harmonics option
= -1 Use Cartesian basis functions to construct
symmetry-adapted linear combination (SALC)
of basis functions. The SALC space is the
linear variation space used. (default)
= 0 Use spherical harmonic functions to create
SALC functions, which are then expressed
in terms of Cartesian functions. The
contaminants are not dropped, hence this
option has EXACTLY the same variational
space as ISPHER=-1. The only benefit to
obtain from this is a population analysis
in terms of pure s,p,d,f,g functions.
= +1 Same as ISPHER=0, but the function space
is truncated to eliminate all contaminant
Cartesian functions [3S(D), 3P(F), 4S(G),
and 3D(G)] before constructing the SALC
functions. The computation corresponds
to the use of a spherical harmonic basis.
QMTTOL = linear dependence threshhold
Any functions in the SALC variational space whose
eigenvalue of the overlap matrix is below this
tolerence is considered to be linearly dependent.
Such functions are dropped from the variational
space. What is dropped is not individual basis
functions, but rather some linear combination(s)
of the entire basis set that represent the linear
dependent part of the function space. The default
is a reasonable value for most purposes, 1.0E-6.
When many diffuse functions are used, it is common
to see the program drop some combinations. On
occasion, in multi-ring molecules, we have raised
QMTTOL to 3.0E-6 to obtain SCF convergence, at the
cost of some energy.
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.)
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
DGRID = flag to add extra digits in molecular orbitals to
the log file for use by Kohout's DGrid program:
http://www2.cpfs.mpg.de/~kohout/dgrid.html
This is one of the modern alternatives to the old
AIMPAC codes, in the QTAIM/ELF arena.
(default .FALSE.)
NUMCOR = an array setting up the number of core orbitals
for each element (up to Z=103). NUMCOR is only
used in MP2 currently.
Default: -1 for all atoms, which means use the
default values hardwired in the MP2 code.
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.
* * * computation control switches * * *
For the most part, the default is the only sensible
value, and unless you are sure of what you are doing,
these probably should not be touched.
NPRINT = Print/punch control flag
See also EXETYP for debug info.
(options -7 to 5 are primarily debug)
= -7 Extra printing from Boys localization.
= -6 debug for geometry searches
= -5 minimal output
= -4 print 2e-contribution to gradient.
= -3 print 1e-contribution to gradient.
= -2 normal printing, no punch file
= 1 extra printing for basis,symmetry,ZMAT
= 2 extra printing for MO guess routines
= 3 print out property and 1e- integrals
= 4 print out 2e- integrals
= 5 print out SCF data for each cycle.
(Fock and density matrices, current MOs
= 6 same as 7, but wider 132 columns output.
This option isn't perfect.
= 7 normal printing and punching (default)
= 8 more printout than 7. The extra output
is (AO) Mulliken and overlap population
analysis, eigenvalues, Lagrangians, ...
= 9 everything in 8 plus Lowdin population
analysis, final density matrix.
NOSYM = 0 the symmetry specified in $DATA is used
as much as possible in integrals, SCF,
gradients, etc. (this is the default)
= 1 the symmetry specified in the $DATA 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.
ETOLLZ = threshold to label molecular orbitals by Lz
values. Small matrices of the Lz operator are
diagonalized for the sets of MOs whose orbital
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.
NORMF = 0 normalize the basis functions (default)
= 1 no normalization
NORMP = 0 input contraction coefficients refer to
normalized Gaussian primitives. (default)
= 1 the opposite.
ITOL = primitive cutoff factor (default=20)
= n products of primitives whose exponential
factor is less than 10**(-n) are skipped.
ICUT = n integrals less than 10.0**(-n) are not
saved on disk. (default = 9). Direct
SCF will calculate to a cutoff 1.0d-10
or 5.0d-11 depending on FDIFF=.F. or .T.
ISKPRP = 0 proceed as usual
1 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)
MODDOS = a bit-additive option to calculate
density of states (DOS)
+ 1: calculate total DOS
+ 2: automatically set emin and emax to the whole
MO range (see DOSPAR(5))
+ 4: set the Fermi level to 0
+ 8: omit lines with no density
Default: 0
DOSPAR = an array of 5 parameters for DOS calculations
(1) minimum energy (eV)
(2) maximum energy (eV)
(3) energy grid step (eV)
(4) broadening Gaussian sigma (eV)
(5) (1) and (2) are set automatically using
MO energies padded with this value (MODDOS=2).
Default: 0,0,0,0,0
* * * restart options * * *
IREST = restart control options
(for OPTIMIZE run restarts, see $STATPT)
Note that this option is unreliable!
= -1 reuse dictionary file from previous run,
useful with GEOM=DAF and/or GUESS=MOSAVED.
Otherwise, this option is the same as 0.
= 0 normal run (default)
= 1 2e restart (1-e integrals and MOs saved)
= 2 SCF restart (1-,2-e integrals and MOs saved)
= 3 1e gradient restart
= 4 2e gradient restart
GEOM = select where to obtain molecular geometry
= INPUT from $DATA input (default for IREST=0)
= DAF read from DICTNRY file (default otherwise)
As noted in the first chapter, binary file restart is
not a well tested option!
==========================================================
==========================================================
$SYSTEM group (optional)
This group provides global control information for
your computer's operation. This is system related input,
and will not seem particularly chemical to you!
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 Section 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
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)
PARALL = 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.
KDIAG = diagonalization control switch
= 0 use a vectorized diagonalization routine
if one is available on your machine,
else use EVVRSP. (default)
= 1 use EVVRSP diagonalization. This may
be more accurate than KDIAG=0.
= 2 use GIVEIS diagonalization
(not as fast or reliable as EVVRSP)
= 3 use JACOBI diagonalization
(this is the slowest method)
= 4 use DSYEV diagonalization from LAPACK
= 5 use DSYEVD for DFTB Fock matrices and
DSYEV elsewhere (both from LAPACK).
To use LAPACK, you should link it.
You may be able to edit comp and lked to
take advantage of your LAPACK library.
Some math libraries such as MKL, OpenBLAS,
and Fujitsu SSL come with LAPACK.
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.
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 option should not be used on multicore
machines, and is really outdated now.
(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)
MODIO = 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;
however, in CPKS F22 is always opened).
16 - reduce I/O (skip various printing)
32 - do not open file F15 (in SOSCF) on workers/agents
and do not close on all.
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.
256 - reduce timing output.
512 - use XYZ file to store coordinates.
1024 - enforce minimal output
2048 - a single flag to turn on various accceleration
options
4096 - redirect OUTPUT and PUNCH output to /dev/null
on all parallel ranks except the global master.
Use only if OS permits it (a UNIX).
Mainly intended for GDDI runs to prevent masters
from writing anything. It is assumed that the
global master writes all important data.
16384 - force workers/agents to write error messages when aborting
a job due to an error; mostly useful in GDDI.
MEM10 = words used to store dictionary file F10 in memory.
Selecting this option will replace any I/O for F10
by memory access.
Default: 0 (disk-based F10)
MEM22 = words used to store grid file F22 for DFT in memory.
Selecting this option will skip any I/O for F22.
Mainly useful for DFT Hessian and analytic FMO-DFT
gradient, both using CPKS. MEM22 is not used for
F22 in DFT energy/gradient (use MODIO=8 for removing
I/O to F22).
Default: 0 (disk-based F22)
NALIGN = align dynamic memory blocks so that relative addresses
are multiples of NALIGN (absolute addresses are
given by OS, usually a multiple of 64 bits).
Default: 1 (effectively, no allignment)
==========================================================
==========================================================
$BASIS group (optional)
This group allows certain standard basis sets to be
easily requested. Basis sets are specified by:
a) GBASIS plus optional supplementations such as NDFUNC,
b) BASNAM to read custom basis sets from your input,
c) EXTFIL to read custom bases from an external file,
d) 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.
= 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.
= MINIX The basis set for HF-3C.
GBASIS=MINIX will not add three energy corrections
(to add them, use GBASIS=HF-3C).
Supported elements are H-Kr, and MINIX means
H-Ne, B-Ne: MINIS (a minimum basis set)
Li-Be, Na-Mg, Al-Ar: MINIS+p (plus polarization)
K-Zn: SV (a minimum basis set)
Ga-Kr: SVP (SV + polarization)
= HF-3C The minimal basis set MINIX and three
add-on energy corrections (3c). These
three corrections are D3(BJ), GCP and SRB.
See $DFT's dispersion corrections.
* * 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.
= CCnC - As CCn, but augmented with tight functions
for recovering core and core-valence
correlation, e.g. cc-pCVnZ.
Available H-Ar for n=D,T,Q, also n=5 for H-Ne.
= 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!
Sapporo valence basis sets:
= SPK-nZP - Sapporo family of non-relativistic
bases, n=D,T,Q, available H-Xe
= SPK-AnZP - diffuse augmentation of the above.
= SPKrnZP - 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=LUT-IOTC.
= SPKrAnZP - diffuse augmentation of the above.
See extended notes below!
Sapporo core/valence basis sets:
= SPK-nZC - 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-Rn.
To be used only with a relativistic
transformation of the integrals, such
as RELWFN=LUT-IOTC.
= 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 i-
functions. 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
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=RESC, DK,
IOTC, IOTC2E, LUT-IOTC, or LUTIOTC2,
while using NR sets for elements lighter
than Sc. DK versions of ACCD or ACCT are available for Sc-
Zn (but not for the rest of the row, 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 core/valence sets are available for all
atoms including the 6th row of the periodic table (H-Rn).
4. It is extremely illogical to use any of the all-electron
relativistic bases without turning on scalar relativity!
So, choose RELWFN=LUT-IOTC (or IOTC, IOTC2E, LUTIOTC2,
DK, RESC) in $CONTRL.
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 4s-
4f,5s-5d,6s for f block (lanthanides). This suggests you
should change the default number of core orbitals, such as
NACORE in $MP2 or NCORE in $CCINP, to correlate the
indicated semi-core orbitals (this is not automatic).
6. The relativistic sets ("r") are identical to the non-
relativistic choices ("-") for atoms H-Ar, where scalar
relativity has almost no effect on orbital shapes.
7. The relativistic bases were optimized at the 3rd order
of the Douglas-Kroll transformation, with a Gaussian nuclei
model. It should be fine to use them with any RESC, DK,
IOTC, IOTC2E, LUT-IOTC, or LUTIOTC2 calculation.
8. 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.
9. The SPK basis sets were extracted from the data base of
Segmented Gaussian Basis Sets, maintained by Takeshi Noro,
University of Hokkaido, Sapporo, Japan:
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
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-DZP
MCP/NOSeC-V-TZP MCP-TZP
MCP/NOSeC-V-QZP MCP-QZP
MCP/NOSeC-V-TZP+diffuse MCP-ATZP
MCP/NOSeC-V-QZP+diffuse MCP-AQZP
MCP/NOSeC-CV-DZP MCPCDZP
MCP/NOSeC-CV-TZP MCPCTZP
MCP/NOSeC-CV-QZP MCPCQZP
MCP/NOSeC-CV-DZP+diffuse MCPACDZP
MCP/NOSeC-CV-TZP+diffuse MCPACTZP
MCP/NOSeC-CV-QZP+diffuse MCPACQZP
GBASIS = IMCP-SR1 and IMCP-SR2 -
valence basis sets to be used with the improved
MCPs with scalar relativistic effects.
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 aug-
cc-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 all-
electron 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 - selects MNDO model Hamiltonian
= AM1 - selects AM1 model Hamiltonian
= PM3 - selects PM3 model Hamiltonian
= RM1 - selects RM1 model Hamiltonian
= DFTB - selects 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
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)
= DUNNING (default for GBASIS=DH, DZV)
= HUZINAGA (default for GBASIS=MINI, MIDI)
= HONDO7 (default for GBASIS=TZV)
SPLIT2 = an array of splitting factors used when NDFUNC
or NPFUNC is 2. Default=2.0,0.5
SPLIT3 = an array of splitting factors used when NDFUNC
or NPFUNC is 3. Default=4.00,1.00,0.25
The splitting factors are from the Pople school, and are
probably too far apart. See for example the Binning and
Curtiss paper. For example, the SPLIT2 value will usually
cause an INCREASE over the 1d energy at the HF level for
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:
$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
! unpolarized basis for hydrogens
$ligH
n31 6
$end
$ecp
Cr-ecp SBKJC
O-ecp none
...snipped... there must be 6 O's given here
O-ecp none
H-ecp none
...snipped... there must be 12 H's given here
H-ecp none
$end
* * *
EXTFIL = a flag to read basis sets from an external file,
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.
===========================================================
==========================================================
$DATA group (required)
$DATAS group (if NESC chosen, for small component basis)
$DATAL group (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.
----------------------------------------------------------
-1- TITLE a single descriptive title card.
----------------------------------------------------------
-2- GROUP, NAXIS
GROUP is the Schoenflies symbol of the symmetry group,
you may choose from
C1, Cs, Ci, Cn, S2n, Cnh, Cnv, Dn, Dnh, Dnd,
T, Th, Td, O, Oh.
NAXIS is the order of the highest rotation axis, and
must be given when the name of the group contains an N.
For example, "Cnv 2" is C2v. "S2n 3" means S6. 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", Section 3.
Armed with only the name of the group, GAMESS is able
to exploit the molecular symmetry throughout almost all of
the program, and thus save a great deal of computer time.
GAMESS does not require that you know very much else about
group theory, although a deeper knowledge (character
tables, irreducible representations, term symbols, and so
on) is useful when dealing with the more sophisticated
wavefunctions.
Cards -3- and -4- are quite complicated, and are rarely
given. A *SINGLE* blank card may replace both cards -3-
and -4-, to select the 'master frame', which is defined on
the next page. If you choose to enter a blank 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).
----------------------------------------------------------
-3- X1, Y1, Z1, X2, Y2, Z2
For C1 group, there is no card -3- or -4-.
For CI group, give one point, the center of inversion.
For CS group, any two points in the symmetry plane.
For axial groups, any two points on the principal axis.
For tetrahedral groups, any two points on a two-fold axis.
For octahedral groups, any two points on a four-fold axis.
----------------------------------------------------------
-4- X3, Y3, Z3, DIRECT
third point, and a directional parameter.
For CS group, one point of the symmetry plane,
noncollinear with points 1 and 2.
For CI group, there is no card -4-.
For other groups, a generator sigma-v plane (if any) is
the (x,z) plane of the local frame (CNV point groups).
A generator sigma-h plane (if any) is the (x,y) plane of
the local frame (CNH and dihedral groups).
A generator C2 axis (if any) is the x-axis of the local
frame (dihedral groups).
The perpendicular to the principal axis passing through
the third point defines a direction called D1. If
DIRECT='PARALLEL', the x-axis of the local frame coincides
with the direction D1. If DIRECT='NORMAL', the x-axis of
the local frame is the common perpendicular to D1 and the
principal axis, passing through the intersection point of
these two lines. Thus D1 coincides in this case with the
negative y axis.
----------------------------------------------------------
The 'master frame' is just a standard orientation for
the molecule. By default, the 'master frame' assumes that
1. z is the principal rotation axis (if any),
2. x is a perpendicular two-fold axis (if any),
3. xz is the sigma-v plane (if any), and
4. xy is the sigma-h plane (if any).
Use the lowest number rule that applies to your molecule.
Some examples of these rules:
Ammonia (C3v): the unique H lies in the XZ plane (R1,R3).
Ethane (D3d): the unique H lies in the YZ plane (R1,R2).
Methane (Td): the H lies in the XYZ direction (R2). Since
there is more than one 3-fold, R1 does not apply.
HP=O (Cs): the mirror plane is the XY plane (R4).
In general, it is a poor idea to try to reorient the
molecule. Certain sections of the program, such as the
orbital symmetry assignment, do not know how to deal with
cases where the 'master frame' has been changed.
Linear molecules (C4v or D4h) must lie along the z axis,
so do not try to reorient linear molecules.
You can use EXETYP=CHECK to quickly find what atoms are
generated, and in what order. This is typically necessary
in order to use the general $ZMAT coordinates.
* * * *
Depending on your choice for COORD in $CONTROL,
if COORD=UNIQUE, follow card sequence U
if COORD=HINT, follow card sequence U
if COORD=CART, follow card sequence C
if COORD=ZMT, follow card sequence G
if COORD=ZMTMPC, follow card sequence M
Card sequence U is the only one which allows you to define
a completely general basis here in $DATA.
Recall that UNIT in $CONTRL determines the distance units.
----------------------------------------------------------
-5U- Atom input. Only the symmetry unique atoms are
input, GAMESS will generate the symmetry equivalent atoms
according to the point group selected above.
if COORD=UNIQUE NAME, ZNUC, X, Y, Z
***************
NAME = 10 character atomic name, used only for printout.
Thus you can enter H or Hydrogen, or whatever.
ZNUC = nuclear charge. It is the nuclear charge which
actually defines the atom's identity.
X,Y,Z = Cartesian coordinates.
if COORD=HINT
*************
NAME,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
(origin and unit points on axes of master frame).
defaults: POINT1='O', POINT2='I', POINT3='J'
ref- R.L. Hilderbrandt, J.Chem.Phys. 51, 1654 (1969).
You cannot understand HINT input without reading this.
Note that if ZNUC is negative, the internally stored
basis for ABS(ZNUC) is placed on this center, but the
calculation uses ZNUC=0 after this. This is useful
for basis set superposition error (BSSE) calculations.
----------------------------------------------------------
* * * If you gave $BASIS, continue entering cards -5U-
until all the unique atoms have been specified.
When you are done, enter a " $END " card.
* * * If you did not, enter cards -6U-, -7U-, -8U-.
----------------------------------------------------------
-6U- GBASIS, NGAUSS, (SCALF(i),i=1,4)
GBASIS 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".
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-.
----------------------------------------------------------
-8U- A blank card ends the basis set for this atom.
----------------------------------------------------------
Continue entering atoms with -5U- through -8U- until all
are given, then terminate the group with a " $END " card.
--- this is the end of card sequence U ---
COORD=CART input:
----------------------------------------------------------
-5C- Atom input.
Cartesian coordinates for all atoms must be entered. They
may be arbitrarily rotated or translated, but must possess
the actual point group symmetry. GAMESS will reorient the
molecule into the 'master frame', and determine which
atoms are the unique ones. Thus, the final order of the
atoms may be different from what you enter here.
NAME, ZNUC, X, Y, Z
NAME = 10 character atomic name, used only for printout.
Thus you can enter H or Hydrogen, or whatever.
ZNUC = nuclear charge. It is the nuclear charge which
actually defines the atom's identity.
X,Y,Z = Cartesian coordinates.
----------------------------------------------------------
Continue entering atoms with card -5C- until all are
given, and then terminate the group with a " $END " card.
--- this is the end of card sequence C ---
COORD=ZMT input: (GAUSSIAN style internals)
----------------------------------------------------------
-5G- ATOM
Only the name of the first atom is required.
See -8G- for a description of this information.
----------------------------------------------------------
-6G- ATOM i1 BLENGTH
Only a name and a bond distance is required for atom 2.
See -8G- for a description of this information.
----------------------------------------------------------
-7G- ATOM i1 BLENGTH i2 ALPHA
Only a name, distance, and angle are required for atom 3.
See -8G- for a description of this information.
----------------------------------------------------------
-8G- ATOM i1 BLENGTH i2 ALPHA i3 BETA i4
ATOM is the chemical symbol of this atom. It can be
followed by numbers, if desired, for example Si3.
The chemical symbol implies the nuclear charge.
i1 defines the connectivity of the following bond.
BLENGTH is the bond length "this atom-atom i1".
i2 defines the connectivity of the following angle.
ALPHA is the angle "this atom-atom i1-atom i2".
i3 defines the connectivity of the following angle.
BETA is either the dihedral angle "this atom-atom i1-
atom i2-atom i3", or perhaps a second bond
angle "this atom-atom i1-atom i3".
i4 defines the nature of BETA,
If BETA is a dihedral angle, i4=0 (default).
If BETA is a second bond angle, i4=+/-1.
(sign specifies one of two possible directions).
----------------------------------------------------------
o Repeat -8G- for atoms 4, 5, ...
o The use of ghost atoms is possible, by using X or BQ
for the chemical symbol. Ghost atoms preclude the
option of an automatic generation of $ZMAT.
o The connectivity i1, i2, i3 may be given as integers,
1, 2, 3, 4, 5,... or as strings which match one of
the ATOMs. In this case, numbers must be added to the
ATOM strings to ensure uniqueness!
o In -6G- to -8G-, symbolic strings may be given in
place of numeric values for BLENGTH, ALPHA, and BETA.
The same string may be repeated, which is handy in
enforcing symmetry. If the string is 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.
----------------------------------------------------------
-9G- A blank line terminates the Z-matrix section.
----------------------------------------------------------
-10G- STRING VALUE
STRING is a symbolic string used in the Z-matrix.
VALUE is the numeric value to substitute for that string.
----------------------------------------------------------
Continue entering -10G- until all STRINGs are defined.
Note that any blank card encountered while reading -10G-
will be ignored. GAMESS regards all STRINGs as variables
(constraints are sometimes applied in $STATPT). It is not
necessary to place constraints to preserve point group
symmetry, as GAMESS will never lower the symmetry from
that given at -2-. When you have given all STRINGs a
VALUE, terminate the group with a " $END " card.
--- this is the end of card sequence G ---
* * * *
The documentation for sequence G above and sequence M
below presumes you are reasonably familiar with the input
to GAUSSIAN or MOPAC. It is probably too terse to be
understood very well if you are unfamiliar with these. A
good tutorial on both styles of Z-matrix input can be
found in Tim Clark's book "A Handbook of Computational
Chemistry", published by John Wiley & Sons, 1985.
Both Z-matrix input styles must generate a molecule
which possesses the symmetry you requested at -2-. If
not, your job will be terminated automatically.
COORD=ZMTMPC input: (MOPAC style internals)
----------------------------------------------------------
-5M- ATOM
Only the name of the first atom is required.
See -8M- for a description of this information.
----------------------------------------------------------
-6M- ATOM BLENGTH
Only a name and a bond distance is required for atom 2.
See -8M- for a description of this information.
----------------------------------------------------------
-7M- ATOM BLENGTH j1 ALPHA j2
Only a bond distance from atom 2, and an angle with 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.
----------------------------------------------------------
-8M- ATOM BLENGTH j1 ALPHA j2 BETA j3 i1 i2 i3
ATOM, BLENGTH, ALPHA, BETA, i1, i2 and i3 are as described
at -8G-. However, BLENGTH, ALPHA, and BETA must be given
as numerical values only. In addition, BETA is always a
dihedral angle. i1, i2, i3 must be integers only.
The j1, j2 and j3 integers, used in MOPAC to signal
optimization of parameters, must be supplied but are
ignored here. You may give them as 0, for example.
----------------------------------------------------------
Continue entering atoms 3, 4, 5, ... with -8M- cards until
all are given, and then terminate the group by giving a
" $END " card.
--- this is the end of card sequence M ---
==========================================================
This is the end of $DATA!
If you have any doubt about what molecule and basis set
you are defining, or what order the atoms will be
generated in, simply execute an EXETYP=CHECK job to find
out!
==========================================================
$ZMAT group (required if NZVAR is nonzero in $CONTRL)
This group lets you define the internal coordinates in
which the gradient geometry search is carried out. These
need not be the same as the internal coordinates used in
$DATA. The coordinates may be simple Z-matrix types,
delocalized coordinates, or natural internal coordinates.
You must input a total of M=3N-6 internal coordinates
(M=3N-5 for linear molecules). NZVAR in $CONTRL can be
less than M IF AND ONLY IF you are using linear bends. It
is also possible to input more than M coordinates if they
are used to form exactly M linear combinations for new
internals. These may be symmetry coordinates or natural
internal coordinates. If NZVAR > M, you must input IJS and
SIJ below to form M new coordinates. See DECOMP in $FORCE
for the only circumstance in which you may enter a larger
NZVAR without giving SIJ and IJS.
**** IZMAT defines simple internal coordinates ****
IZMAT is an array of integers defining each coordinate.
The general form for each internal coordinate is
code number,I,J,K,L,M,N
IZMAT =1 followed by two atom numbers. (I-J bond length)
=2 followed by three numbers. (I-J-K bond angle)
=3 followed by four numbers. (dihedral angle)
Torsion angle between planes I-J-K and J-K-L.
=4 followed by four atom numbers. (atom-plane)
Out-of-plane angle from bond I-J to plane J-K-L.
=5 followed by three numbers. (I-J-K linear bend)
Counts as 2 coordinates for the degenerate bend,
normally J is the center atom. See $LIBE.
=6 followed by five atom numbers. (dihedral angle)
Dihedral angle between planes I-J-K and K-L-M.
=7 followed by six atom numbers. (ghost torsion)
Let A be the midpoint between atoms I and J, and
B be the midpoint between atoms M and N. This
coordinate is the dihedral angle A-K-L-B. The
atoms I,J and/or M,N may be the same atom number.
(If I=J AND M=N, this is a conventional torsion).
Examples: N2H4, or, with one common pair, H2POH.
Example - a nonlinear triatomic, atom 2 in the middle:
$ZMAT IZMAT(1)=1,1,2, 2,1,2,3, 1,2,3 $END
This sets up two bonds and the angle between them.
The blanks between each coordinate definition are
not necessary, but improve readability mightily.
**** the next define delocalized coordinates ****
DLC is a flag to request delocalized coordinates.
(default is .FALSE.)
AUTO is a flag to generate all redundant coordinates,
automatically. The DLC space will consist of all
non-redundant combinations of these which can be
found. The list of redundant coordinates will
consist of bonds, angles, and torsions only.
(default is .FALSE.)
NONVDW is an array of atom pairs which are to be joined
by a bond, but might be skipped by the routine
that automatically includes all distances shorter
than the sum of van der Waals radii. Any angles
and torsions associated with the new bond(s) are
also automatically included.
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:
IXZMAT is an extra array of simple internal coordinates
which you want to have added to the list generated
by AUTO. Unlike NONVDW, IXZMAT will add only the
coordinate(s) you specify.
IRZMAT is an array of simple internal coordinates which
you would like to remove from the AUTO list of
redundant coordinates. It is sometimes necessary
to remove a torsion if other torsions around a bond
are being frozen, to obtain a nonsingular G matrix.
IFZMAT is an array of simple internal coordinates which
you would like to freeze. See also FVALUE below,
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
gives the matrix S= 1.0 1.0 0.0
0.0 0.0 1.0
1.0 -1.0 0.0
which defines the symmetric stretch, asymmetric stretch,
and bend of water.
references for natural internal coordinates:
P.Pulay, G.Fogarasi, F.Pang, J.E.Boggs
J.Am.Chem.Soc. 101, 2550-2560(1979)
G.Fogarasi, X.Zhou, P.W.Taylor, P.Pulay
J.Am.Chem.Soc. 114, 8191-8201(1992)
reference for delocalized coordinates:
J.Baker, A. Kessi, B.Delley
J.Chem.Phys. 105, 192-212(1996)
==========================================================
==========================================================
$LIBE group (required if linear bends are used in $ZMAT)
A degenerate linear bend occurs in two orthogonal planes,
which are specified with the help of a point A. The first
bend occurs in a plane containing the atoms I,J,K and the
user input point A. The second bend is in the plane
perpendicular to this, and containing I,J,K. One such
point must be given for each pair of bends used.
APTS(1)= x1,y1,z1,x2,y2,z2,... for linear bends 1,2,...
Note that each linear bend serves as two coordinates, so
that if you enter 2 linear bends (HCCH, for example), the
correct value of NZVAR is M-2, where M=3N-6 or 3N-5, as
appropriate.
==========================================================
==========================================================
$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.
DIIS = selects Pulay's DIIS interpolation.
SOSCF = selects second order SCF orbital optimization.
Only one of DIIS or SOSCF may be .TRUE. in any run.
Which is chosen by default depends on the run:
for RHF, GVB, UHF, or ROHF (if Abelian): SOSCF.TRUE.
for any DFT, or for non-Abelian groups: DIIS=.TRUE.
NOCONV = .TRUE. means neither SOSCF nor DIIS will be used,
specify NOCONV=.TRUE. DIIS=.FALSE. SOSCF=.FALSE.
in the rare case you don't wish to use either one.
NOCONV's default is .FALSE., meaning the program
will obey your choice of DIIS/SOSCF, or else pick
its default for them.
NOSCAL = a flag to avoid scaling rotations in SOSCF.
Once either DIIS or SOSCF are initiated, the following
less important accelerators are placed in abeyance:
EXTRAP = selects Pople extrapolation of the Fock matrix.
DAMP = selects Davidson damping of the Fock matrix.
SHIFT = selects level shifting of the Fock matrix.
RSTRCT = selects restriction of orbital interchanges.
DEM = selects direct energy minimization, which is
implemented only for RHF. (default=.FALSE.)
defaults for EXTRAP DAMP SHIFT RSTRCT DIIS SOSCF
ab initio: T F F F F/T T/F
semiempirical: T F F F F F
The above parameters are implemented for all SCF
wavefunction types, except that DIIS will work for GVB only
for those cases with NPAIR=0 or NPAIR=1.
FSHIFT = a level static shift parameter of the Fock matrix
in Hartree (default = 0.0).
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, UDFT or MRSF-TDDFT
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)
DIRTHR = energy threshold to force a recomputation of the
Fock matrix if energy does not decrease with DIRSCF.
(deafult = 0 hartree, disabling the option)
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 more info on the convergence methods,
see the 'Further Information' section.
* * * * * * * * * * * * * * * * * * * * *
---- 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 (non-relativistic) and H-Rn
(relativistic), for RHF, ROHF, and GVB wavefunctions, as
well as DFT runs (see also VVOS in $MCSCF). IVVTYP in
$LOCAL allows for the explicit selection non-relativistic
or relativistic orbitals. VVOs are restricted to all electron
basis sets and model core potentials; ECPs may not be used.
VVOS should be much better MCSCF starting orbitals than
either MVOQ or ACAVO type virtuals.
MVOQ = 0 Skip MVO generation (default)
= n Form modified virtual orbitals, using a cation
with n electrons removed. Implemented for
RHF, ROHF, and GVB. If necessary to reach a
closed shell cation, the program might remove
n+1 electrons. Typically, n will be about 6.
= -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.
ACAVO = Flag to request Approximate Correlation-Adapted
Virtual Orbitals. Implemented for RHF, ROHF,
and GVB (w/o direct SCF). Default is .FALSE.
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.
Of course, canonical virtuals are PACAVO(1)=1,1,2,2,-1,-1.
SFOVVO = Flag to request the substitution of virtual
canonical orbitals for VVOs in energy calculations
for SF-ORMAS and SF-ORMAS-PDFT.
* * *
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.
NHOMO = 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.
NLUMO = 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 and ROHF.
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.
note for MRSF-TDDFT: When there are nearby triplet states
the SCF process can choose a wrong triplet state. This can
cause a discontinuity in the potential energy surface.
MOM=.t. of ROHF SCF tries to keep the same reference ROHF
wavefunction of the initial MOREAD set to avoid the
complication.
However, one can also use this keyword more aggresively
for choosing a particular MO set to be optimized, which
would produce new respose states that could not have been
obtained without the option. This can be especially useful
when there are near-degenerate triplet states.
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
----- GVB wavefunction input -----
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).
NSETO = The number of sets of open shells in the
function. Maximum of 10. (default=0)
NO = An array giving the degeneracy of each open
shell set. Give NSETO values.
(default=0,0,0,...).
NPAIR = The number of geminal pairs in the -GVB-
function. Maximum of 12. The default
corresponds to open shell SCF (default=0).
CICOEF = An array of ordered pairs of CI coefficients
for the -GVB- pairs.
(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 = An vector of fractional shell occupations.
ALPHA = An array of A coupling coefficients, given in
lower triangular order.
BETA = An array of B coupling coefficients, given in
lower triangular order.
Note: The default for F, ALPHA, and BETA depends on
the state chosen. Defaults for the most commonly occurring
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 = option for output to the PUNCH file
= 0 do not punch out the final orbitals
= 1 punch out the occupied orbitals
= 2 punch out occupied and virtual orbitals
The default is NPUNCH = 2.
NPREO = energy and orbital printing options, applied
after other output options, for example NPRINT=-5
for no orbital output overrules NPREO. NPREO
affects only printing, see NPUNCH just above.
Default: 1,9999,2,1 (meaning print all orbitals,
but no separate list of orbital energies).
Orbitals from NPREO(1) to NPREO(2) and orbital
energies from NPREO(3) to NPREO(4) are printed.
Positive values are explicit orbitals, while
negative numbers are relative to the HOMO. Here,
HOMO means NE/2, ie RHF-like counting, no matter
what the SCFTYP actually is:
NPREO(1)=25,35 prints orbitals 25 to 35
NPREO(1)=-3,-4 prints 8 orbitals: from 3
below the HOMO, the HOMO,
to 4 above the HOMO.
NPREO(3) to NPREO(4) define separate print-out of
the orbital energies, by default this is skipped,
since the starting value is higher than the end.
To print only HOMO and "LUMO" LCAO coefficients,
and all orbital energies, enter:
NPREO(1)=0,-1,1,9999
----- 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:
SCALF = initial exponent scale factor when VTSCAL is
in use, useful when restarting. The default
is 1.0.
MAXVT = maximum number of iterations (at a single
geometry) to satisfy the energy virial theorem.
The default is 20.
VTCONV = convergence criterion for the VT, which is
satisfied when 2 + + R x dE/dR is less
than VTCONV. The default is 1.0D-6 Hartree.
For more information on this option, which is most useful
during a geometry search, see M.Lehd and F.Jensen,
J.Comput.Chem. 12, 1089-1096(1991).
* * * * * * * * * * * * * * * * * * *
For more discussion of GVB/ROHF input
see the 'further information' section
* * * * * * * * * * * * * * * * * * *
==========================================================
==========================================================
$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.
NF = an array containing the number of doubly
occupied
MOs for each fragment.
MF = an array containing the number of atomic basis
functions located on each fragment.
ITER = maximum number of SCF-MI cycles, overriding
the usual MAXIT value. (default is 50).
DTOL = SCF-MI density convergence criteria.
(default is 1.0d-10)
ALPHA = possible level shift parameter.
(default is 0.0, meaning shifting is not used)
DIISON = a flag to active the DIIS convergence.
(default is .TRUE.)
MXDIIS = the maximum number of previous effective Fock
and
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 = maximum number of iteration in the localization
step (Default is 50)
CNVLOC = convergence parameter for the localization.
(default is .01).
IOPT = prints additional debug information.
= 0 standard outout (default)
= 1 print for each SCF-MI cycle MOs, overlap
between the MOs, CPU times.
= 2 print some extra informations in secular
systems solution.
==========================================================
"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)
==========================================================
$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
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):
= SVWN SLATER exchange + VWN5 correlation
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
= HCTH Hamprecht/Cohen/Tozer/Handy's 1998 mod
to B97, omitting HF exchange (m=4)
= 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)
= HCTH407P also known as HCTH/407+, designed by
Boese/Chandra/Martin/Marx for
describing hydrogen bonds
= HCTH14 also known as HCTH p=1/4, made by
Menconi/Wilson/Tozer
= HCTH76 also known as HCTH p=7/6, made by
Menconi/Wilson/Tozer
= 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
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.
= APF mixed functional based on PBE0 and B3PW91
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, LC-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
= CAMQTP00 Verma/Bartlett's reparametrization of
CAMB3LYP for vertical ionization
potential
= CAMQTP01 Jin/Bartlett's reparametrization of
CAMB3LYP for vertical ionization
potential
M11 is also range-separated, but is listed below with the
other meta-GGAs.
"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
correct second order gradient expansion.
= M11 M11 range-separated hybrid
= M11-L M11 local (0% HF exchange) with
dual-range exchange
= MN12-L uses a nonseparable functional form aiming to provide
balanced peformance for both chemistry and solid-state
physics applications
= MN12-SX screened-exchanged (SX) hybrid functional with 20% HF
exchange for short-range and 0% HF exchange for
long-range
= MN15 uses a nonseparable functional form with 44% HF
exchange. This fuctional is a global hybride (no
range-separation).
= MN15-L a local functional with 0% HF exchange.
= REVM06 revised hybrid M06
= REVM06-L revised local M06-L
= REVM11 revised M11 with 40% HF exchange
More than one hundred functionals are available by:
= USELIBXC See more details in $LIBXC
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.
https://en.wikipedia.org/wiki/Minnesota_functionals
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, DC and LRDFLG. DC uses empirically chosen C6 and C8
coefficients, while LRDFLG 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.
Four 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 if DC is chosen and IDCVER isn't given.
= 4 means modified 3rd implementation DFT-D3(BJ).
(-4 is used for DFT-D3(BJ) for HF-3c).
Setting IDCVER will force DC=.TRUE.
GCP = a flag for the geometric counterpoise scheme
correction in HF-3c.
SRB = a flag for short-range basis set incompleteness
(SRB) correction in HF-3c.
DCCHG = a flag to use Chai-Head-Gordon damping function
instead of Grimme's 2006 function. Pertinent only
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) non-
additive 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)
and for DFT-D3(BJ):
S. Grimme, S. Ehrlich and L. Goerigk,
J.Comput.Chem. 32, 1456-1465 (2011)
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.
DCSR = 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.
DCS6 = 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.
DCS8 = 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.
DCA1 = a1 parameter appearing in the -D3(BJ) dispersion
model. Optimized values are automatically
selected for a set of known functionals,
otherwise the default is 0.50.
DCA2 = a2 parameter appearing in the -D3(BJ) dispersion
model. Optimized values are automatically
selected for a set of known functionals,
otherwise the default is 4.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.
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).
MU = 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)
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:
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 = amount of HF exchange. (default=0.53)
CMP2 = 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.
RADTYP = type of radial quadrature
= MHL : Murray, Handy and Laming radial grid (default)
= MK3 : Mura and Knowles Log3 radial grid
= TA : Treutler and Ahlrichs radial grid
Lebedev grid:
NRAD = number of radial points in the Euler-MacLaurin
quadrature. (default=96)
NLEB = 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...
The same switch can be used to select spherical grids
developed by A.S.Popov:
246, 264, 342, 432 : rotaional octahedral symmetry grids
350 and 398 : full octahedral symmetry grids (like Lebedev grids)
212 : rotational icosahedral grid
rotational icosahedral grids with 192 and 242 points
are currently disabled in the source code
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.
JANS = 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.
polar coordinate grid:
NRAD = number of radial points in the Euler-MacLaurin
quadrature. (96 is reasonable)
NTHE = number of angle theta grids in Gauss-Legendre
quadrature (polar coordinates). (12 is reasonable)
NPHI = number of angle phi grids in Gauss-Legendre
quadrature. NPHI should be double NTHE so points
are spherically distributed. (24 is reasonable)
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.
Setting to zero disables DFT grid switching.
NRAD0 = same as NRAD, but defines initial coarse grid.
default = smaller of 24 and NRAD/4
NLEB0 = same as NLEB, but defines initial coarse grid.
default = 110
NTHE0 = same as NTHE, but defines initial coarse grid.
default = smaller of 8, NTHE/3
NPHI0 = same as NPHI, but defines initial coarse grid.
default = smaller of 16, NPHI/3
molecular grid construction parameters:
BFCTYP = specify algorithm of molecular grid construction
from atomic grids (Becke's fuzzy cell algorithm)
= BECKE : use Becke's fuzzy cell method from the original
paper.
= SSF : use algorithm by Stratmann, Scuseria and Frisch.
It usually screens out more points than Becke's
algorithm
= LEGACY : use legacy DFT code (default)
BECKE and SSF switches allow multithreaded DFT run, but
does not support DC, LRD and analytical Hessian calculation
types, as well as polar spherical grid with non-C1 symmetry.
PARTFN = select grid partitioning function
relevant only if BFCTYP=SSF or BECKE
= BECKE : use polynomial proposed in Becke's original work.
Default when BFCTYP=BECKE
= SSF : use polynomial proposed by Stratmann, Scuseria
and Frisch. Default when BFCTYP=SSF
= ERF : use erf-based partitioning function based on
unpublished ERF1 algorithm in NWChem
= SMSTEP2 : use smoothstep-2 polynomial (very effective,
but aggressive partitioning)
= SMSTEP3 : use smoothstep-3 polynomial
= SMSTEP4 : use smoothstep-4 polynomial
= SMSTEP5 : use smoothstep-5 polynomial
technical parameters:
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.
WTDER = switch on/off grid weigths derivative contribution to
the nuclear gradient.
= .TRUE. : compute grid weights derivatives
= .FALSE. : do not compute grid weights derivatives
Default is .TRUE. if BFCTYP=LEGACY, and .FALSE.
otherwise.
Note: BFCTYP=BECKE and BFCTYP=SSF currently do not
support grid weight derivatives.
Grid weights are only important for low-quality grids.
The default grid settings are good enough to disable
grid weight derivatives calculation. This switch can
be turned off to improve convergence of geometry
optimization, if the optimized structure is not the
one with lowest energy.
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
AUXFUN = AUX0 uses no auxiliary basis set for resolution
of the identity, limiting accuracy.
= AUX3 uses the 3rd generation of RI basis sets,
These are available for the elements H to
Ar, but have been carefully considered for
H-Ne only. (DEFAULT)
THREE = a flag to use a resolution of the identity to
turn four center overlap integrals into three
center integrals. This can be used only if
no auxiliary basis is employed. (default=.FALSE.)
==========================================================
==========================================================
$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.
Like spin-flip TD-DFT, the orbitals for a reference
state of Mixed-Reference Spin-Flip (MRSF) TDDFT (TDDFT=MRSF
in $CONTRL) are obtained by either UHF or ROHF type, with
MULT in $CONTRL determining its Ms quantum number. Only
triplet reference (MULT=3 in $CONTRL) state is supported
now. While spin-flip TDDFT only utilizes Ms=+1 of triplet
state in the generation of response states, MRSF-TDDFT
takes both Ms=+1 and -1, which eliminates a bulk of the
notorious spin-contamination of SF-TDDFT. Accordingly,
MRSF-TDDFT utilizes both alpha to beta as well as beta to
alpha spin-filp excitations in the generation of response
determinants with Ms=0, 2. Using these, the spin state
specific response target states with singlet, triplet as
well as quintet spin states can be generated indivisually.
Thus, target states have S=0, 1 with Ms=0 and S=2 with
Ms=2 (MULT=1, 3, 5 in $TDDFT). Analytic nuclear energy
gradients, (RUNTYP=GRADIENT in $CONTRL), 2-state and
3-state conical-intersection search (see $CONICL),
ionization potential and electron affinity based on extended
Koopmans' theorem, nonadiabatic MD simulation (see $MD), and
some properties are available with MRSF-TDDFT. However,
gradient is only available for ROHF.
The respective papers are
Energy and Gradient of MRSF-TDDFT:
Lee S, Filatov M, Lee S, Choi CH.
J.Chem.Phys.2018 Sep 14;149(10):104101.
Lee S, Kim EE, Nakata H, Lee S, Choi CH.
J.Chem.Phys.2019 May 14;150(18):184111.
2- and 3 States Conical Intersection:
Baek YS, Lee S, Filatov M, Choi CH.
J.Phys.Chem.A2021 Mar 2;125(9):1994-2006.
Ionization Potential:
Pomogaev V, Lee S, Shaik S, Filatov M, Choi CH.
J.Phys.Chem.Lett.2021 Oct 7;12:9963-72.
Nonadiabatic MD:
Park W, Lee S, Huix-Rotllant M, Filatov M, Choi CH.
J.Phys.Chem.Lett.2021 Apr 30;12(18):4339-46.
===========================================================
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
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, B3LYPV1R, PBE, PBE0
CAMQTP00, CAMQTP01
For evaluation of just the excitation energies, you may use
many more functionals, notably including the metaGGAs in
the last three lines:
NONE
SVWN, SVWN1RPA, SPZ81, SP86, SOP, SLYP,
BVWN, BVWN1RPA, BPZ81, BP86, BOP, BLYP, OLYP,
B3LYP, CAMB3LYP, B3LYPV1R, B3PW91, X3LYP,
PW91, PBE, PBE0,
APF,
VS98, PKZB,
M05, M05-2X, M06, M06-HF, M06-L, M06-2X, M08-HX, M08-SO,
MN12-L, MN12-SX, MN15, MN15-L,
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. DFTB
(possibly combined with PCM) may be chosen as well, and
analytic gradients for singlet and triplet are available.
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 and MRSF:
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. Note that
spin-flip TD-DFT in the Tamm/Dancoff approximation using
DFTTYP=NONE is equivalent to spin-flip CIS.
In the case of SFNFLP, MULT below is ignored, since the
equation of motion of SFNFLP cannot restrict a particular
spin state. Therefore, the spin state of resulted response
states should be determined by inspection.
On the other hand, MRSF utilizes it to specify target spin
state, since MRSF formulates independent equation of motion
for each response spin states. Both spin-flip and MRSF codes
operate 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.
It is important to know that the first response
state corresponds to the ground state in the case
of TDDFT=SPNFLP or MRSF.
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!
MULT = Multiplicity of the excited states specifying spin
state of response states, i.e.,
1 (Singlet) or 3 (Triplet) in TDDFT=EXCITE and
1 (Singlet), 3 (Triplet) or 5 (Quintet) in
TDDFT=MRSF (default = 1). It should be emphasized
that the singlet response states by MRSF include
the ground singlet state.
However, this parameter is ignored when
TDDFT=SPNFLP. Therefore, the spin state of
resulted response states should be determined by
inspection in the case of TDDFT=SPNFLP.
TDPRP = 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. It includes
natural orbitals for TDDFT=MRSF.
UNRLX = a flag that requests the use of the unrelaxed
density in property computations for the state IROOT.
UNRLX=.TRUE. also requires TDPRP=.TRUE.
ALTRN = a flag that requests property computations on unrelaxed
densities of all NSTATE excited states. ALTRN=.TRUE. also
requires TDPRP=.TRUE. and UNRLX=.TRUE.
TRNSD = a flag that requests the use of the transition density
in property computations. The transition density is
taken for the transition between the ground state
and the state IROOT. All nuclear contributions to
properties are set to zero when using the transition
density. By default this flag is set to .FALSE.
TRNSD=.TRUE. also requires TDPRP=.TRUE. If this flag
is combined with ALTRN, properties are calculated
for transition densities of the transitions between
the ground state and all NSTATE excited states as well
as for all state-to-state transitions between the NSTATE
excited states.
IFEDAT = an array of four integers that define atom ranges of
donor/acceptor for FED. The format is: start of donor,
end of donor, start of acceptor, end of acceptor.
For the donor/acceptor definition to be valid,
the ranges must be valid and non overlapping. By
default IFEDAT(1)=0,0,0,0 is set. This means that
no donor/acceptor definition is given.
FED = a flag that requests the calculation of
fragment excitation difference (FED), fragment charge
difference (FCD), and Multi-FED-FCD couplings. All
three calculations will be done by setting
FED=.TRUE. . FED=.TRUE also requires UNRLX=.TRUE.
TPA = a flag requesting two-photon absorption cross-
sections. These are computed for each of the
NSTATE excited states, after first evaluating
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 and MRSF,
which is only coded in the Tamm/Dancoff
approximation.
NONEQ 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.
MODTD = bit-additive options
1 use fast sort of virt - occ MO energies
Default: 0
MREKT = a flag calling Dyson Orbitals (DOs) of Extended
Koopmans' Theorem (EKT) using MRSF-DFT Excited State
Response for Ionization Energies. (default, .FALSE.)
Note, that NPRINT=8 in $CONTRL prints out a final
MRSF Density Matrix of a selected electronic state
as well as overlaping between MRSF-DOs and initial
canonical MOs if there is a $VEC block.
MRDEA = a flag calling DOs of EKT using MRSF-DFT Excited State
Response for Electron Affinities. (default, .FALSE.)
It requares DIRSCF=.FALSE. in $SCF group.
* * * 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 = requests the polarizability. Default=.FALSE.
If BETA is not chosen, give just one PFREQ.
BETA = requests the hyperpolarizability. Default=.FALSE.
Two values are required for PFREQ.
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).
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)
NLEB = number of angular points in the Lebedev grid.
(default=110)
NTHE = number of theta grid points if a polar coordinate
grid is used.
NPHI = number of phi grid points if a polar coordinate
grid is used. NPHI should be twice NTHE.
SG1 = flag selecting "standard grid one".
(default=.FALSE.)
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:
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).
ZVTOL = the convergence criteria of Z-Vector routine.
(default is 1.0d-10. But it appears to be too
tight. So 1.0d-8 is recommended.)
IXCORE = an array of occupied orbitals that are included
in the excitations during MRSF-TDDFT calculation.
This is an option of REW (Restricted Energy
Window) technique, which can limit excitations of
response theory. MRSF-TDDFT utilizes a ROHF
reference, which is composed of doubly occupied
(C), two singly occupied (O) and virtual orbitals
(V). Only the elements of IXCORE array are allowed
in the excitations from C to O and V.
For example, IXCORE(1)=1,10 only allows
the transitions from 1 and 10th doubly occupied
orbitals to O and V. All the transitions that
starts from the other orbitals of C are excluded
in the response space.
==========================================================
==========================================================
$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 and FMO-DFTB. Excited state calculations for
full DFTB may be performed through the standard (linear-
response) time-dependent formalism (only closed shell). PCM
can be used for both ground and excited state calculations,
and energy and gradient can be evaluated.
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.
NDFTB 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.
LCDFTB = a flag to turn on long range corrected DFTB, as
derived from LC-DFT. At present, only DFTB2 can be
combined with this option. See EMU below.
Note that LC-DFTB requires parameters optimized
for it.
Default: .FALSE.
EMU = a long-range separation parameter for LC-DFTB (like
MU in $DFT). To use LC-DFTB properly, some non-zero
value (for example, 0.3) should be assigned to EMU.
The value of EMU can be found in the readme file
for the parameters used.
Default: 0.
DAMPXH = 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.
DAMPEX = an exponent used in the damping function for X-H
atomic pairs. The default value is 4.0 (taken
from the 3OB parameter set).
SRSCC = 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 Slater-
Koster tables.
Using .TRUE. enables the use of proper Hubbard
values for p and d orbitals, implemented only
for DFTB1 and DFTB2.
Default: .FALSE.
ITYPMX Convergence method of SCC calculations.
= -1 Use standard GAMESS convergence methods.
SOSCF and DIIS are supported, but DEM is not.
= 0 Broyden's method.
Interpolation is applied for atomic
(or shell-resolved when SRSCC=.TRUE.)
charges, but not Hamiltonian matrix.
= 1 (reserved)
= 2 DIIS for charges.
Default: 0.
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.
OCCTHR = threshold for identifying doubly occupied, partially
occupied or unoccupied molecular orbitals with ETEMP.
If an occupation number n is OCCTHR1 (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.
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.
The OMPRIMP2 program is an OpenMP threaded version of the
RIMP2 code. TO use the OMPRIMP2 program you must build
GAMESS with GMS_OPENMP set to true in the following files:
- $GMS_DIR/install.info
- $GMS_DIR/Makefile
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)
CUTOFF = transformed integral retention threshold, the
default is 1.0d-9 (1.0d-12 in FMO runs).
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 SCS-
MP2 potential surface. For CODE=IMS, the nuclear gradient
can be evaluated analytically. See NUMGRD in $CONTRL if
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.
CPHFBS = 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.
NWORD = controls memory usage. The default uses all
available memory. Applies to CODE=SERIAL.
(default=0)
METHOD= n selects transformation method, 2 being the
segmented transformation, and 3 being a more
conventional two phase bin sort implementation.
3 requires more disk, but less memory. The
default is to attempt method 2 first, and
method 3 second. Applies only to CODE=SERIAL.
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.
==========================================================
===========================================================
$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.
USEDM = 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.
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.
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)
IVMTD = selects the procedure for removing redundancies
when inverting the two-center, two-e- matrix.
= 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.
===========================================================
$HPCCHEM group (optional, relevent for libcchem v2.0)
This group controls the behavior of libcchem v2.0
NMODE = 0 do not use libcchem v2.0 (default)
= 1 enable libcchem v2.0
ETHRSH = energy SCF convergence threshold
defult value = 1.0d-06
DTHRSH = starting value for dynamic integrals screening
threshold. As calculation proceeds through
SCF iterations the value decreses
to the minimum value (10e-10).
NPRINT = print mode. Default = 1
higher numbers print less detail.
See LIBCCHEM 2.0 documentation for details
of print modes.
SCF = use libcchem's SCF routines (default false)
FOCK = use libcchem's fock build routines (default false)
===========================================================
$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!
In FMO/AP, CABNAM(2) is used to set the auxiliary basis set
for the larger basis in AP.
Note IAUXBF in $RIMP2 for selecting spherical harmonics
versus Cartesian Gaussians.
In the RI-G3(MP2,CCSD(T) (also called the PDG method), rimp2
energy and gradient are computed. The gradient level is rimp2/6-31g(d)/aux_grd,
the energy level is rimp2/g3l/aux_eng. Here, g3l is the large G3 AO
basis set, aux_grd and aux_eng are auxiliary basis sets customized
for gradient and energy computations, respectively. For good performance
and accuracy aux_grd is usually smaller than aux_eng. To input aux_grd
and aux_eng set:
cabnam(1)=aux_grd,aux_eng
One recommendation in our publication
[Phys. Chem. A 2021, 125, 42, 9421-9429] is to use:
cabnam(1)=ccd,cct
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.
===========================================================
==========================================================
$CCINP group (optional, relevant for any CCTYP)
This group controls the 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, only valence electrons are correlated.
This group is also necessary to execute ACP calculations.
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, CR-EOM, or
CR-EOML, and the EA, IP, DEA, and DIP-EOMCC computations
read $CCINP to define orbital ranges for the underlying
ground-state CCSD steps preceding EOMCC diagonalizations
under $EOMINP's control.
A number of CCTYP choices have been superseded by more
advanced or improved approaches. For example, R-CC and
CR-CC were developed prior to their CR-CCL replacement,
while CR-EOML supersedes CR-EOM. CR-CCL provides an
accurate approximation to the fully iterative CCSDT method,
superior to the widely used CCSD(T), especially in bond
breaking situations. If it is important to account for the
coupling of singly, doubly, and triply excited clusters,
which non-iterative corrections to CCSD neglect, or if a
highly accurate representation of full CCSDT energetics is
desired, it is best to use the CCT3 approach. If one is
interested in methods with up to connected triple
excitations, the most reasonable choices are:
RHF ROHF (high spin)
--- ----------------
ground states [properties]:
CCSD [CCPRP] CCSD [n/a]
CCSD(T) n/a
CR-CCL CR-CCL
CCSD3A CCSD3A
CCT3 CCT3
excited states [properties]:
EOM-CCSD [CCPRPE] EOM-CCSD
CR-EOML n/a
If one is interested in strongly correlated systems (e.g.,
models of metal-insulator transitions), it is worth considering
the ACP options, namely,
ACCSD
ACC23
ACCSD3A
ACCT3
which work for both RHF and ROHF references (see the ACP
keyword below for further information).
If one is interested in a spin-adapted and accurate description
of ground and excited states of an open-shell system obtained by
attaching one electron to or removing one electron from the parent
closed-shell core, such as a radical, the recommended choices,
which should be used with SCFTYP=RHF, are:
EA-EOM3A
IP-EOM3A
If one is interested in a spin-adapted and accurate description
of ground and excited states of an open-shell system obtained by
attaching two electrons to or removing two electrons from the
parent closed-shell core, such as a biradical, the recommended
choices, which can work with RHF or ROHF orbitals of the target
system or RHF orbitals of the underlying closed-shell core, are:
DEAEOM3A
DEAEOM3B
DIPEOM3A
although the less expensive DEAEOM2, DEAEOM2A, and DIPEOM2 options
may be useful too, for example, as part of composive schemes aimed
at larger systems that mix different DEA- and DIP-EOMCC levels.
For a brief explanation of the above acronyms, see CCTYP. See
Section 4 of this manual and papers cited in the outputs for
more information.
When reading CR-EOM and CR-EOML outputs, it is recommended
to focus on the DELTA-CR-EOMCCSD(T) excitation energies, their
rigorously size-intensive delta-CR-EOMCC(2,3) counterparts, and
the delta-CR-EOMCC(2,3) total energies.
Parallel computation using DDI 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.
See the "Further Information" section of this manual
for more details about coupled-cluster runs.
**** The first six keywords 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.
NFZV = the number of frozen virtual orbitals to be
omitted from the calculation. (default is 0)
MAXCC = defines the maximum number of CCSD (or LCCD, CCD,
CCSDt, ACCSD, ACCSDt) iterations. This parameter
also applies to ROHF's left CCSD vector solver,
but not RHF's left CCSD vector. See MAXCCL for RHF.
(default=30 for LCCD, CCD, and CCSD, and 50 for
CCSD3A, ACCSD, and ACCSD3A)
ICONV = 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.)
NACTO = number of active occupied orbitals to be used
in CCSDt (CCTYP=CCSD3A), CC(t;3) (CCTYP=CCT3),
and selected ACP (CCTYP=ACCSD3A, ACCT3)
calculations. This number corresponds to the beta
active occupied orbitals, which is the same as
the number of active orbitals that are doubly
occupied in the RHF or ROHF reference. The alpha
active occupied orbitals are automatically
calculated using the formula (MULT - 1) + NACTO,
which means that the singly occupied alpha orbitals
in open-shell systems are automatically treated as
active, even if NACTO=0 [CCSDt and CC(t;3) become
CCSD and CR-CC(2,3), respectively, by setting
NACTO=-1].
NACTU = number of active unoccupied orbitals to be used
in CCSDt (CCTYP=CCSD3A), CC(t;3) (CCTYP=CCT3),
and selected ACP (CCTYP=ACCSD3A, ACCT3)
calculations. This number corresponds to the alpha
active unnocupied orbitals, which is the same as
the number of active orbitals which are empty in
the RHF or ROHF reference. The beta active
unoccupied orbitals are automatically calculated
according to the formula (MULT - 1) + NACTU.
SHIFT = a small positive value (usually <1.0 Hartree) added
to the denominators entering the iterative amplitude
equations of the CCSDt (CCTYP=CCSD3A), CC(t;3)
(CCTYP=CCT3), and ACP (CCTYP=ACCSD, ACC23, ACCSD3A,
and ACCT3) systems to help convergence in cases of
more severe quasi-degeneracies. (default is 0)
ACP = an array of five real numbers specifying the scaling
factors multiplying each of the five T2^2 diagrams
projected on doubles [see, e.g., Figure 1 of
I. Magoulas, J. Shen, and P. Piecuch, Mol. Phys.
120, e2057365 (2022) for the numbering convention
of these T2^2 diagrams]. There is no default for
the ACP array. When one of the ACP options
(CCTYP=ACCSD, ACC23, ACCSD3A, and ACCT3) is selected,
an five-number array, excluding (1.0,1.0,1.0,1.0,1.0),
must be provided. For the study of strongly
correlated systems, the recommended options include
ACP(1)=0.0,0.0,0.0,1.0,1.0
ACP(1)=1.0,0.0,0.0,1.0,0.0
and
ACP(1)=1.0,0.0,no/(no+nu),nu/(no+nu),0.0
where no and nu are the numbers of correlated occupied
and unoccupied orbitals, respectively. This keyword
is ignored unless one of the ACP options is selected.
**** the next group of keywords pertains to SCFTYP=RHF only ****
CCPRP = a flag to select computation of the CCSD-level
ground-state density matrix (see, also, CCPRPE
in $EOMINP for EOMCCSD-level density matrices
involving excited states). The computation of
the required left eigenstates of CCSD takes
extra time, 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
EOMCCSD 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. Use of CCTYP=CR-EOM
with CCPRP=.TRUE. or CCTYP=CR-EOML will calculate triples
corrections and CCSD-level properties.
There is usually little reason to select any of these:
MAXCCL = maximum number of iterations for the left CCSD
calculations needed to determine CCSD properties
or CR-CCL energies in SCFTYP=RHF runs.
This is just a synonym for MAXEOM in $EOMINP.
If you want to alter the left CCSD state's convergence
threshold for RHF, 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.
IREST = defines the restart option for CC runs using SCFTYP=RHF.
If the value of IREST is greater or equal 3, program
will restart from the earlier CC SCFTYP=RHF run.
This option requires saving the disk file CCREST
from the previous CC run. Values of IREST between
0 and 3 should not be used. The value of IREST also
sets the iteration counter in the restarted run.
The default is 0, meaning no restart is attempted.
IREST does not apply to CCSD3A, CCT3, ACCSD, ACC23,
ACCSD3A, and ACCT3 runs using SCFTYP=RHF, which are
automatically redirected to SCFTYP=ROHF and MULT=1.
The CCSD3A, CCT3, ACCSD, ACC23, ACCSD3A, and ACCT3
runs can be restarted using KREST and, in the case of
CCT3, ACC23, and ACCT3, KREST and LREST described
below, independent of SCFTYP.
MXDIIS = defines the number of cluster amplitude vectors
from previous iterations to be included in the
DIIS extrapolation during the CCSD (or LCCD, CCD)
SCFTYP=RHF run. 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. MXDIIS does not
apply to CCSD3A, ACCSD, and ACCSD3A iterations.
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. AMPTSH does not apply to CCSD3A,
ACCSD, and ACCSD3A iterations.
CCERI = defines whether or not a tensor decomposition
technique should be applied to the two-electron
repulsion integrals (2ERIs) for CC calculations.
= STANDARD, no tensor decomposition, the standard
4-center 2ERIs are used (default).
= RI, the resolution-of-the-identity approximation
is employed and the 4-center 2ERIs are assembled
as products of 3-center RI integrals on the fly.
The RI approximation largely eliminates memory
bottlenecks of CC calculations and offers a better
parallel scaling and CPU usage.
Don't forget to see also the $RICC input group for
additional parameters.
**** the next group of keywords pertains to SCFTYP=ROHF only ****
There is little reason to select any of these.
MULT = spin multiplicity to use in the reference
determinant during the CC computation. 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 = number of correlated occupied spatial
orbitals and NU = number of correlated unoccupied
spatial orbitals)
= 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 or CCSDt iterations
(default)
= 1 restart from AMPROCC ($JOB.F70) file of a
previous CCSD run, or from the $JOB.MOE file of a
previous ACCSD run when CCTYP=ACCSD or ACC23, a
previous CCSDt run when CCTYP=CCSD3A or CCT3, or a
previous ACCSDt run when CCTYP=ACCSD3A or ACCT3.
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 iterations
(CCTYP=CCSD with CCPRP=.TRUE. or CR-CCL) or
the left CCSD-like iterations using singly and
doubly excited cluster amplitudes obtained in
CCSDt (CCTYP=CCT3) or ACP (CCTYP=ACC23, ACCT3)
calculations. (default)
= 1 restart from AMPROCC ($JOB.F70) file of a
previous left CCSD (CCTYP=CCSD with CCPRP=.TRUE.
or CR-CCL) or left CCSD-like (CCTYP=CCT3,
ACC23, and ACCT3) run.
LREST=1 must be combined with KREST=1.
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.
**** The IOSIZE input parameter applies to all CC ****
calculations.
IOSIZE = n enables "chunking" of direct access I/O during
CC computations. The direct access (DA) record
size of several scratch files generated during CC
computations can become quite large, but many
compilers will not permit DA record lengths in excess
of 2GB. If the record size L of a DA file exceeds n,
then each record will be written to/read from disk
in multiple chunks, each of which is at most n bytes
in length. For example, if the logical record length
L = 2,400,000,000 and n = 1,000,000,000, then each
logical record will be processed in 3 chunks, with
respective record lengths of 10^9, 10^9, and 4^8 bytes,
thereby avoiding the 2GB record length limit.
= 0 disables I/O chunking (default.)
An exetyp=check job with iosize=0 will indicate if
chunking needs to be enabled and will also list
recommended values of n.
==========================================================
==========================================================
$RICC group
(optional, for SCFTYP=RHF and CCTYP=CCSD/CCSD(T) in $CONTRL
and essential if CCERI=RI in $CCINP)
This group controls the calculation of ground state
energies for closed-shell molecules (SCFTYP=RHF) using the
coupled-cluster method within the singles and doubles (CCSD)
truncation scheme or CCSD augmented with an approximate
triples correction: CCSD(T). These calculations run in parellel
using a hybrid MPI/OpenMP parallelization model and employ
the resolution-of-the-identity (RI) approximation for two-
electron repulsion integrals (2ERIs).
The following three paremeters are related to the definition
of the auxiliary basis set to be used for the RI approximation
to 2ERIs.
IAUXBF = 0 uses Cartesian Gaussians
= 1 uses spherical harmonics
for the auxiliary basis set. The default is
inherited from ISPHER.
CABNAM specifies built-in auxiliary basis sets:
= 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!
In FMO/AP, CABNAM(2) is used to set the auxiliary basis set
for the larger basis in AP.
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.
*** Please note that the current RI integral evaluation
module in GAMESS can only deal with up to G shells.
Therefore, remember to discard basis functions with
higher angular momenta if they are present in your
external auxiliary basis.
*** The auxiliary basis sets optimized for RIMP2 calculations
(extension -RI) introduce smaller errors to the RI-CCSD
or RI-CCSD(T) energies compared to the RI-JK basis sets.
*** The common practice is to use the same class of auxiliary
basis set as the AO basis set, e.g., cc-pVDZ (AO) and
cc-pVDZ-RI (i.e., GBASIS=CCD in $BASIS and CABNAM=CCD in
$RICC). Please note that the chosen auxiliary basis set
controls the accuracy of the RI approximation. If a higher
accuracy is desired, we recommend using a larger auxiliary
basis set, e.g., if GBASIS=CCD in $BASIS, then CABNAM in
$RICC may be set to CCT or ACCT. Importantly, increasing
the size of the auxiliary basis set only marginally increases
the computational cost of CC calculations, as long as the
AO basis remains the same.
The following five paremeters are related to the evaluation and
storage of the three-center RI integrals.
OTHAUX = flag to orthogonalize the auxiliary basis set by
diagonalization of the overlap matrix. If there
is reason to suspect linear dependence may exist
in the auxiliary basis, set OTHAUX to .TRUE. to have
a more numerically stable result. Larger auxiliary
basis sets such as CCT and ACCT, in particular, may
benefit from selecting this. Since the CCn/ACCn
auxiliary basis sets are very commonly used for
coupled-cluster calculations, and they often entail
linear dependence, the default is set to .TRUE.
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)
IVMTD = selects the procedure for removing redundancies
when inverting the two-center, two-e- matrix.
= 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)
USEDM3 = flag to control the storage of the three-center RI
integrals in distributed memory or in the shared
memory of each MPI/DDI process.
= .FALSE. if stored in the shared memory of each MPI/
DDI process (default).
= .TRUE. if stored in distributed memory. For this
choice, the MEMDDI parameter in $SYSTEM is a
mandatory input.
The three-center RI integrals do not consume a lot of memory
and do not present the true memory bottleneck for RI-CC calculations.
This makes their storage in the shared memory of each MPI/DDI
process feasible. This is also a good choice, as this reduces
communication between the MPI/DDI processes. Hence, the default
value of USEDM3 is set to .FALSE.
For large molecules, however, choosing USEDM3=.TRUE. might be a
good option for making calculations feasible. For example,
if the shared memory requirement per MPI/DDI processor is nearly
equal to the available memory on a single node of your computer
cluster, choosing USEDM3=.TRUE. would reduce the shared memory
need per MPI/DDI process to some extent. Note that this helps
only a little and not a lot, since the real memory bottleneck
of RI-CC calculations do not come from the three-center RI
integrals.
The following paremeter is related to the choice of an algorithm
for RI-CCSD iterations.
USEDM4 = flag to control whether or not to temporarily store
blocks of preassembled four-center 2ERIs in the shared
memory of the MPI/DDI proccesses during CCSD iterations.
= .TRUE., needs larger memory per MPI/DDI process but the
algorithm is faster (default).
= .FALSE. a memory-economic algorithm but takes a longer
wall time/CCSD iteration. This algorithm completely avoids
the storage of four-center 2ERIs and assembles them from
the 3-center RI integrals on the fly.
Choosing USEDM4=.FALSE. is useful if the computational resources
are limited, e.g., for computer clusters with a small number
of compute nodes. In this case, selecting USEDM4=.FALSE. would
allow you to perform your desired calculation at the expense of
some additional wall time.
*** We strongly recommend running an EXETYP=CHECK job always
prior to submitting a real RI-CC calculation. This provides
a detailed information on the memory requirements for the
various options described above and a guideline to set
the parameters MWORDS and MEMDDI in the $SYSTEM input group.
The RI-CC program is OpenMP threaded. Hence, to use this program
you must build GAMESS with GMS_OPENMP set to true in the following
files:
- $GMS_DIR/install.info
- $GMS_DIR/Makefile
==========================================================
==========================================================
$EOMINP group
(optional, for CCTYP=EOM-CCSD, CR-EOM, or CR-EOML)
(optional, for CCTYP=EA-EOMx or IP-EOMx)
(optional, for CCTYP=DEAEOMx or DIPEOMx)
(optional for CCSD properties, or CCTYP=CR-CCL)
This group controls the calculation of excited electronic
states by the equation-of-motion (EOM) coupled-cluster (CC) with
single and double excitations, with optional triples corrections.
It also pertains to single and double electron attachment and
detachment processes, which may result in the system being
left in an excited state. EOM-CCSD, EA-EOMx, IP-EOMx, DEAEOMx,
and DIPEOMx can be selected for RHF and ROHF reference states,
however the use of EA-EOMx and IP-EOMx for open-shell references
is strongly discouraged. The remaining CCTYP options listed above,
namely CR-EOM and CR-EOML, can only be used with SCFTYP=RHF.
These EOM-type runs consist of an SCF calculation to obtain
the appropriate reference state entering a subsequent ground-state
CCSD computation (see the $CCINP group to control the ground-state
calculation and the correlated orbital range) followed by EOMCCSD,
IP-EOMCC, EA-EOMCC, DIP-EOMCC, or DEA-EOMCC calculations to
determine the relevant excited, attached, or ionized electronic
states of interest (see NSTATE below). In the case of EOMCCSD,
triples corrections based on the method of moments approach may
follow.
This input group permits the user to select how many states
are computed, however keep in mind that machine time is linear
in the number of requested states. Because the default state
selection is rather simplistic (only one state in the totally
symmetric representation), it is usually necessary to provide
inputs for this group, which includes selecting the NSTATE and
IROOT keywords sensibly.
Because this input group is used for several CCTYP
calculations, not all keywords are used in every case and
certain keywords may have slightly different meanings. In
particular,
a) if the reference type is RHF and CCTYP is EOM-CCSD, CR-
EOM, 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, CVGCI, MAXCI, and
MICCI are ignored. MULT is sometimes ignored.
c) if CCTYP is EA-EOMx or IP-EOMx (x=2,3A), keywords CCPRPE,
MTRIP, MEOM, MCI, MACT, CVGCI, MAXCI, MICCI are ignored.
d) if CCTYP is DEAEOMx (x=2,2A,3,3A,3B) or DIPEOMx (x=2,3,3A),
keywords GROUP, IROOT, CCPRPE, MTRIP, MEOM, MCI, MOACT,
MICEOM, CVGCI, NAXCI, MICCI are ignored.
Additional information on CC and EOMCC methods can be
found in the "Further Information" section of this manual.
--- spin and spatial symmetry and state selection
(relevant to CCTYP=EOM-CCSD, CR-EOM, CR-EOML, EA-EOMx,
and IP-EOMx)
MULT = Spin multiplicities for the target states.
The meaning depends on the particular calculation.
The default value for 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 routine will include these
states.
If the calculation is of the IP-EOM or EA-EOM type, the run
will pass through an open-shell code even when the reference
in $CONTRL is given as RHF. The resulting IP or EA states will
be spin-adapted. The IP-EOM and EA-EOM methods can also be used
with ROHF references and MULT=1 in $CONTRL and will similarly
result in spin-adapted states. Although IP-EOM and EA-EOM
calculations will run using genuine open-shell ROHF references
(MULT>1 in $CONTRL), the resulting ionized or attached states
will not be spin-adapted, so we do not recommend to do this.
Therefore, assuming that $CONTRL specifies a closed-shell
reference (accomplished by using either SCFTYP=RHF or ROHF
with MULT=1),
MULT = -1 means determine both doublet and quartet states
= 2 means determine only doublet states
= 4 means determine only quartet states.
In case the run is EOM-CCSD with SCFTYP=ROHF and $CONTRL
selects a closed-shell reference using MULT=1, the meaning
of MULT in $EOMINP is
MULT = -1 determine singlet, triplet, and quintet
excited states,
= 1 determine only singlet excited states,
= 3 consider only triplet excited states,
= 5 consider only quintet excited states.
For ROHF-based EOMCCSD calculations using an open-shell
reference (selected with MULT>1 in $CONTRL), 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. The MULT input will be
be ignored.
In case the run is RHF-based EOM-CCSD (possibly followed by
the triples corrections obtained using CR-EOM or CR-EOML),
the target states are singlets only. In this case, the MULT
keyword in $EOMINP will be ignored.
GROUP the name of the Abelian group to be used, which
may only be 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
state selection, for controlling the EOMCC
calculations, and for labeling the calculated
states in the output. Spatial symmetry is
not used not to speed up the calculations.
NSTATE an array of up to 8 integers that indicates
how many 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
EOMCCSD runs. For EA or IP runs, the NSTATE
input MUST include the ground state of the target
species and may also include excited states.
See also ISELCT below.
There is no particular reason the first excited state (or
ground state in EA/IP runs) should be totally symmetric, so
most calculations require the user to provide a sensible
NSTATE input. Up to 10 states can be found in any irrep.
Keep in mind that 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, using NSTATE(1)=0,0,0,0,0,0,0,0 for RHF
references will calculate the ground state only. This
can be used in CR-EOM runs to generate the type I, II,
or possibly III CR-CCSD(T) energies, provided that
MTRIP (described below) is given appropriately, which
are not otherwise available in standard ground-state
calculations.
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 defines the irrep number, obtained
from the same table as NSTATE, and the second integer
specifies the state number. For example, if
GROUP=C2V, then the input IROOT(1)=3,2 refers to
the second B1 state.
(default is IROOT(1)=1,0)
IROOT's default is moderately sensible for RHF-based
excited-state runs, corresponding to the ground state
(labeled as state 0), which must belong to the totally
symmetric representation. For ROHF-based excitations or EA/IP
runs, this is not generally the case and one should select
IROOT appropriately!
If degenerate EOMCCSD states are detected, only one such
state will be considered for subsequent triples corrections,
if they are requested using CCTYP=CR-EOM or CR-EOML. The
state chosen for possible triples corrections will be the
state with 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 and ISELECT
selects those that will carry into the calculation.
For example, NSTATE(1)=2,2,2,2 with ISELCT(1)=1,3,5,7
prepares two guesses in each irrep, but only iterates the
EOMCCSD equations for the lowest state in each
irrep (the guesses are counted serially).
The next keyword addresses triples corrections to
EOMCCSD. Note that non-iterative triples corrections
to EOMCCSD are only available for SCFTYP=RHF.
MTRIP selects the type of non-iterative triples
corrections to RHF-based EOMCCSD energies
(relevant 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,
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.
Note: The solver used in CISD calculations needed by
MTRIP=2,3,4 is selected by the keyword MCI described
later on.
IEOMTP = Type of the EOM-CC energy:
0 default choice (check the output)
1 CR-CC(2,3),D
2 delta-CR-CC(2,3),D
The energy chosen by IEOMTP will be used as the
final EOM-CC energy, affecting numerical gradients,
FMO runs etc. Note that IEOMTP= 1 or 2 is meaningful
only when the specified values are available (as
determined by CCTYP).
Default: 0
NACT pertains only to EA-EOM3A or IP-EOM3A runs
(for NACT pertaining to DEAEOMx and DIPEOMx, see the
subsection discussing DEA/DIP-EOMCC calculations below):
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)-electron system, NACT refers to the NACT
lowest unoccupied orbitals of the N-electron
reference system.
For CCTYP=IP-EOM3A, used to describe the
(N-1)-electron system, NACT refers to the NACT
highest occupied orbitals of the N-electron
reference system.
The default for NACT is 0, which does not allow any
3p2h or 3h2p operators, thus yielding only EA-EOMCC2
IP-EOMCC2 results.
In other words, you should input a nonzero value
for NACT!
CCPRPE = a flag to select computation of the EOMCCSD-level
excited-state density matrices (see, also, CCPRP in
$CCINP for the CCSD-level ground-state density matrix).
The computation of the required left eigenstates of
EOMCCSD takes extra time, so the default is .FALSE.,
except for CCTYP=CR-EOML, where the work
required for properties must be done anyway.
Notes: CCPRPE can only be used if SCFTYP=RHF and CCTYP=EOM-CCSD,
CR-EOM, or CR-EOML. The property printout includes ground-
and excited-state dipole moments, transition moments and
oscillator strengths involving pairs of calculated states,
occupation numbers for natural orbitals, and selected additional
properties for the state indicated by IROOT (cf. $ELMOM). The
CVGEOM convergence criterion for the left eigenstates of EOMCCSD,
which is set for excitation energies, may need tightening
if the calculated properties require higher precision.
The CC/EOMCC density matrices are square, but not symmetric,
which means that CCSD and EOMCCSD natural orbitals come in
left/right pairs. To reduce the amount of output, only the
right natural orbitals for the state indicated by IROOT
will be found in the output file. The complete set of the
calculated CCSD and EOMCCSD density matrices are printed
into the PUNCH file for possible use elsewhere. Use of
CCTYP=CR-EOM with CCPRPE=.TRUE. or CCTYP=CR-EOML will
calculate triples corrections to CCSD and EOMCCSD as well as
CCSD/EOMCCSD-level properties.
--- iterative solver selection
(relevant to CCTYP=EOM-CCSD, CR-EOM, CR-EOML, EA-EOMx,
and IP-EOMx)
MEOM selects the solver for the EOMCCSD calculations:
0 = one root at a time using a united iterative
space for all calculated roots (default)
1 = one root at a time using a separate iterative
space for each calculated root
2 = the Hirao-Nakatsuji multi-root solver
3 = one root at a time using a separate iterative space
for each computed right/left root. (compare to 1)
4 = one root at a time using a united iterative space
for all right/left roots (compare to 0).
For ROHF-based EOMCCSD or IP/EA runs, there is only one
iterative solver, so MEOM is ignored in these cases.
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 greater than 1.
1 = one root at a time using a separate iterative
space for each calculated root (default)
2 = the Hirao-Nakatsuji multi-root solver (slower).
--- initial guess for EOMCCSD, EA/IP, and possible CISD solvers:
The keywords MINIT and MACT below help select the type of
initial guess used to begin EOMCCSD and EA/IP-EOMCC iterations.
In the context of EOMCCSD runs, "S" and "D" stand for using all
singles and doubles, while "s" and "d" mean restricting those
excitations to and from a smaller number of active orbitals.
Meanwhile, in the context of IP-EOMCCSD or IP-EOMCCSDt and
EA-EOMCCSD or EA-EOMCCSDt calculations, "S" and "D" refer to
using all 1h and 2h1p excitations in the case of IP-EOMCC and
all 1p and 2p1h excitations in the case of EA-EOMCC, while
"s" and "d" denote restricting these kinds of excitations
to and from a smaller set of active 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 because low-lying states with non-negligible double
excitation or significant multi-configurational character
are missing in simplistic CIS, 1h, or 1p guesses and, thus,
may not appear in the final converged EOMCCSD or EA/IP-EOMCC
calculations.
MINIT selects the initial guess procedure for EOMCCSD,
EA-EOMCC, IP-EOMCC, and possibly CISD iterations,
in case MTRIP>1.
1 = Use EOMCCSd to start the EOMCCSD iterations,
use IP-EOMCCSd to start the IP-EOMCCSD or
IP-EOMCCSDt iterations, use EA-EOMCCSd to
start the EA-EOMCCSD or EA-EOMCCSDt iterations,
and CISd to start possible CISD iterations.
2 = Use CIS wave functions to start EOMCCSD as well
as any possible CISD calculations, use 1h wave
functions to start IP-EOMCC iterations, and use
1p wave functions to start EA-EOMCC iterations.
MINIT's default is 2, but MINIT=1 is HIGHLY RECOMMENDED!
MACT = fine tuning of MINIT's EOMCCSD initial guess.
For MINIT=1 MACT=0, use EOMCCSd guess.
For MINIT=1 MACT=1, use EOMCCsd guess.
For MINIT=2 MACT=0, use CIS guess.
For MINIT=2 MACT=1, use CIs guess.
The default for MACT is 0. MACT is ignored for EA/IP-EOMCC
runs.
The next three keywords define active orbitals
used in the initial guesses. 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 carefully to avoid splitting
any degenerate (e.g., pi) orbital shells.
NOACT the number of occupied MOs in the active space
for little "s" or little "d" initial guesses.
NUACT the number of unoccupied MOs in the active space
for little "s" or little "d" initial guesses.
MOACT an array allowing for explicit selection of the
active orbitals used to define the 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,
where the user has searched the virtual space
looking for valence-like orbitals, might be
MINIT=1 NOACT=3 NUACT=5
MOACT(1)=13,14,15, 19,20,24,25,30
---- restarts of EOMCCSD, EA-EOMCC, and IP-EOMCC runs:
JREST = 0 means no restart is used (default).
= 1 means that restart data is read from
AMPROCC file.
One use for JREST=1 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.
--- iteration control:
CVGEOM convergence criterion for the EOMCCSD, EA/IP-EOMCC,
and DEA/DIP-EOMCC excitation amplitudes, which also
applies to the RHF-based left EOMCCSD iterations.
(default=1.0d-4).
CVGEOM also serves as a convergence criterion for
left CCSD calculations using RHF, in which case
the default is 10**(-ICONV).
MAXEOM maximum number of iterations in the EOMCCSD,
EA/IP-EOMCC, and DEA/DIP-EOMCC calculations
(default=50). In the case of RHF-basd EOMCCSD
runs, this keyword takes on a different meaning
depending on the choice of solver
selected using MEOM. 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.
MICEOM maximum number of microiterations in the
EOMCCSD calculations (default=80). Rarely used.
Again, for RHF-based EOMCCSD runs, the meaning
of this keyword changes based on the value of
MEOM. 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., remain within a single iteration cycle).
If for some reason the EOMCCSD convergence is
very slow and the iterative space becomes
too 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.
The next three keywords only apply to RHF-based
CR-EOMCCSD(T) triples corrections (CCTYP=CR-EOM) if MTRIP
is greater than 1:
CVGCI convergence criterion for the CISD expansion
coefficients (default=1.0d-4).
MAXCI 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.
MICCI maximum number of microiterations in the
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).
---- DEA/DIP-EOMCC calculations:
DEA-EOMCC and DIP-EOMCC runs use the CVGEOM and MAXEOM
keywords described above, in addition to the following
keywords described below.
The keyword MULT specifies the spin multiplicities for the
targeted doubly electron attached or doubly ionized states.
If quintet states are sought in the DEA-EOMCC calculations,
be sure to set MINIT=1 and MACT=0 or 1.
MULT = -1 determine singlet, triplet, and quintet states (default)
= 1 determine only singlet states
= 3 determine only triplet states
= 5 determine only quintet states.
NSTATE = the number of states to be computed. The default is
NSTATE=1.
ISELCT = an array to specify the states to be determined. The
default is ISELCT(1)=1,...,NSTATE. For example,
MULT=1 NSTATE=2 ISELCT(1)=1,3 will use the first
and third roots of the singlet symmetry obtained in
the initial guess calculation.
NACT = the number of active MOs defining the active-space
DEA/DIP-EOMCC methods. In the case of DEA-EOMCC
calculations, NACT refers to the number of active
unoccupied MOs used to select the active 3p1h
excitations in DEA-EOMCC(3p1h){Nu} (CCTYP=DEAEOM2A),
the active 3p1h and 4p2h excitations in
DEA-EOMCC(3p1h,4p2h){Nu} (CCTYP=DEAEOM3A), or
active 4p2h excitations in DEA-EOMCC(4p2h){Nu}
(CCTYP=DEAEOM3B). For DIP-EOMCC methods, NACT
refers to the number of active occupied MOs used
to select the active 4h2p excitations in the
DIP-EOMCC(4h2p){No} (CCTYP=DIPEOM3A) calculations.
The default value for NACT is 0, which does not allow any 3p1h, 4p2h or
4h2p operators, and thus yields only DEA-EOMCC(2p) (when CCTYP=DEAEOM2A
or DEAEOM3A), DEA-EOMCC(3p1h) (when CCTYP=DEAEOM3B), or DIP-EOMCC(3h1p)
(when CCTYP=DIPEOM3A) results.
MINIT and MACT select the initial guess procedure used to start the
DEA-EOMCC iterations. These keywords are not applicable to DIP-EOMCC runs,
which are initiated using full 2h and, if NOACT and NUACT are given,
active 3h1p guesses.
The allowed choices of MINIT and MACT in DEA-EOMCC runs are as follows:
For MINIT=1 MACT=0, use full 2p and active 3p1h guess
For MINIT=1 MACT=1, use active 2p and active 3p1h guess
For MINIT=2 MACT=0, use full 2p guess (default)
For MINIT=2 MACT=1, use active 2p guess.
MINIT=2 MACT=0 is a default, but MINIT=1 MACT=0 with the value of NOACT and NUACT
specified is HIGHLY RECOMMENDED! In order for the DEA-EOMCC calculations
to find quintet states using MULT=-1 or 5, MINIT must be set at 1.
NOACT = the number of active occupied MOs defining active
3h1p (CCTYP=DIPEOMx) or active 3p1h (CCTYP=DEAEOMx) initial guesses.
NUACT = the number of active unoccupied MOs defining
active 3h1p (CCTYP=DIPEOMx), active 2p or active 3p1h (CCTYP=DEAEOMx)
initial guesses.
The default values of NOACT and NUACT are 0, however, values for NOACT and NUACT
must be provided unless MINIT=2 and MACT=0.
The keyword JREST below controls restarts of DEA-EOMCC and DIP-EOMCC runs.
The allowed values of JREST are as follows:
JREST = 0 no restart, initial guess procedure will be performed (default)
=-1 stop after completing the initial guess procedure
= 1 restart data is read from $SCR/$JOB/$JOB.R1 file
= 2 initial guess vectors are read from $SCR/$JOB/$JOB.R0 file.
Initial guess vectors are written into $SCR/$JOB/$JOB.R0 file when
JREST=0 or -1.
The reference determinant for the DEA-EOMCC and DIP-EOMCC calculations
defined in $CONTRL can be either of the RHF or the ROHF singlet or
triplet type. The keyword NREF below determines whether the DEA/DIP-EOMCC
calculation is carried out using the MOs of the underlying closed-shell
core or the MOs of the doubly electron attached or doubly ionized target
species.
NREF = 0 singlet RHF MOs (SCFTYP=RHF or SCFTYP=ROHF, MULT=1
in $CONTRL) of the (N-2)-electron (DEA-EOMCC) or
(N+2)-electron (DIP-EOMCC) closed-shell reference system
will be used. Set ICHARGE for the closed-shell
reference system in $CONTRL accordingly. (default)
= 1 singlet RHF (SCFTYP=RHF or SCFTYP=ROHF, MULT=1 in $CONTRL)
or triplet ROHF (SCFTYP=ROHF, MULT=3 in $CONTRL) MOs
of the N-electron target system will be used. In this case,
the input in $CONTRL should describe the N-electron target
system. This choice of NREF requires DIRSCF=.FALSE. in $SCF
group.
NSAVE = 0 temporary binary $SCR/$JOB/$JOB.F* files will be deleted
after all calculations are completed (default)
= 1 temporary binary files are kept; this requires modifying
the rungms file to keep all $SCR/$JOB/$JOB.F* files
= 2 CCSD amplitudes and H-BAR will be read from existing
binary files generated by previous DEA/DIP-EOMCC
calculations (with NSAVE=1). All procedures before
DEA/DIP-EOMCC, including SCF, CCSD, and the
generation of H-BAR, will be skipped.
==========================================================
==========================================================
$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.
==========================================================
==========================================================
$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:
= FMO Read orbitals from the DICTNRY file, from
previous FMO run with MODPRP=1.
= HUCSUB Perform a Huckel guess in each subsystem
of a Divide and Conquer run
= DMREAD 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.)
NORB = The number of orbitals to be read in the $VEC
group. This applies only to GUESS=MOREAD.
For -RHF-, -UHF-, -ROHF-, and -GVB-, NORB defaults to the
number of occupied orbitals. NORB must be given for -CI-
and -MCSCF-. For -UHF-, if NORB is not given, only the
occupied alpha and beta orbitals should be given, back to
back. Otherwise, both alpha and beta orbitals must
consist of NORB vectors.
NORB may be larger than the number of occupied MOs, if you
wish to read in the virtual orbitals. If NORB is less
than the number of atomic orbitals, the remaining orbitals
are generated as the orthogonal complement to those read.
NORDER = Orbital reordering switch.
= 0 No reordering (default)
= 1 Reorder according to IORDER and JORDER.
IORDER = Reordering instructions, 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 = level below which MO coefficients will be set
to zero. (default=1.0E-7)
TOLE = level at which MO coefficients will be equated.
This is a relative level, coefficients are set
equal if one agrees in magnitude to TOLE times
the other. (default=5.0E-5)
SYMDEN = project the initial density in order to generate
symmetric orbitals. This may be useful if the
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.)
==========================================================
==========================================================
$VEC group (optional, relevant for all SCFTYP's)
(required if GUESS=MOREAD)
This group consists of formatted vectors, as written
onto file PUNCH in a previous run. It is considered good
form to retain the titling comment cards punched before
the $VEC card, as a reminder to yourself of the origin of
the orbitals.
For Morokuma decompositions, the names of this group
are $VEC1, $VEC2, ... for each monomer, computed in the
identical orientation as the supermolecule. For transition
moment or spin-orbit coupling runs, orbitals for states
one and possibly two are $VEC1 and $VEC2.
==========================================================
$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 = flag which triggers MO freezing. (default=.FALSE.)
IFRZ = an array of MOs in the input $VEC set which are
to be frozen. There is no default for this.
==========================================================
==========================================================
$STATPT group (for RUNTYP=OPTIMIZE or SADPOINT)
This group controls the search for stationary points.
Note that NZVAR in $CONTRL determines if the geometry
search is conducted in Cartesian or internal coordinates.
METHOD = optimization algorithm selection. Pick from
NR Straight Newton-Raphson iterate. This will
attempt to locate the nearest stationary
point, which may be of any order. There
is no steplength control. RUNTYP can be
either OPTIMIZE or SADPOINT
RFO Rational Function Optimization. This is
one of the augmented Hessian techniques
where the shift parameter(s) is(are) chosen
by a rational function approximation to
the PES. For SADPOINT searches it involves
two shift parameters. If the calculated
stepsize is larger than DXMAX the step is
simply scaled down to size.
QA Quadratic Approximation. This is another
version of an augmented Hessian technique
where the shift parameter is chosen such
that the steplength is equal to DXMAX.
It is completely equivalent to the TRIM
method. (default)
SCHLEGEL The quasi-NR optimizer by Schlegel.
CONOPT, CONstrained OPTimization. An algorithm
which can be used for locating TSs.
The starting geometry MUST be a minimum!
The algorithm tries to push the geometry
uphill along a chosen Hessian mode (IFOLOW)
by a series of optimizations on hyperspheres
of increasingly larger radii.
Note that there currently are no restart
capabilitites for this method, not even
manually.
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.
KDIAGH = eigenvalue solver for the Hessian matrix.
Possible choices are listed in $SYSTEM.
The default is 2 for FMO or DC, and -1 otherwise
(-1 means use KDIAG from $SYSTEM).
NPRICO = Print coordinates every NPRICO optimization steps
to reduce the output file.
(default: 1)
--- 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 = a flag to allow the trust radius to change as
the geometry search proceeds. (default=.TRUE.)
TRMAX = maximum permissible value of the trust radius.
(default=0.5 for OPTIMIZE and 0.3 for SADPOINT)
TRMIN = minimum permissible value of the trust radius.
(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.
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,
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.
IHREP = the number of steps before the hessian is
recomputed. If given as 0, the hessian will
be computed only at the initial geometry if
you choose HESS=CALC, and never again. If
nonzero, the hessian is recalculated every
IHREP steps, with the update formula used on
other steps. (default=0)
HSSEND = a flag to control automatic hessian evaluation
at the end of a successful geometry search.
(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 Print INTR at all, orbitals at all
0 Print INTR at all, orbitals at 0+L (default)
-1 Print INTR at all, orbitals never
-2 Print INTR at 0+L, orbitals never
NPUN = 3 Punch all orbitals and HESS at all
2 Punch all orbitals at all
1 same as 0, plus punch HESS at all
0 Punch all orbitals at 0+L, otherwise only
occupied orbitals (default)
-1 Punch occ orbitals at 0+L only
-2 Never punch orbitals
---- the 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 allowed codes are: 1 (stretches), 2 (angles),
and 3 (dihedrals)
SHMCON = equilibrium constraint values for the distances
specified by IHMCON, given in Angstrom (bond
stretches) or degrees (for angles and dihedrals)
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 do not update Hessian (not recommended)
BFGS default for OPTIMIZE using RFO or QA
POWELL default for OPTIMIZE using NR or CONOPT
POWELL default for SADPOINT
MSP mixed Murtagh-Sargent/Powell update
SCHLEGEL only choice for METHOD=SCHLEGEL
---- 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.
RMIN = Minimum distance threshold. Points whose root
mean square distance from the current point is
less than RMIN are discarded. (default=0.0015)
RMAX = Maximum distance threshold. Points whose root
mean square distance from the current point is
greater than RMAX are discarded. (default=0.1)
RLIM = Linear dependence threshold. Vectors from the
current point to the previous points must not
be colinear. (default=0.07)
==========================================================
* * * * * * * * * * * * * * * * * * * * *
See the 'further information' section for
some help with OPTIMIZE and SADPOINT runs
* * * * * * * * * * * * * * * * * * * * *
==========================================================
$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 = number of parameters to be optimized.
IEX = defines the parameters to be optimized.
If OPTMIZ=BASIS, IEX declares the serial number
of the Gaussian primitives for which the exponents
will be optimized.
If OPTMIZ=GEOMETRY, IEX define the pointers to
the HINT internal coordinates which will be optimized.
(Note that not all internal coordinates have to be
optimized.) The pointers to the internal coordinates
are defined as: (the number of atom on the input
list)*10 + (the number of internal coordinate for that
atom). For each atom, the HINT internal coordinates
are numbered as 1, 2, and 3 for BOND, ALPHA, and BETA,
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.
==========================================================
==========================================================
$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
==========================================================
==========================================================
$DERIV group (optional for numerical derivatives)
This group sets options for numerical gradients. Do
not confuse with $FORCE ($FORCE is for second derivatives).
Note that numerical gradients for FMO cannot use GDDI
(the only option is to use DDI) whereas non-FMO runs can
use GDDI).
VIBSIZ = Displacement size (in bohr). If negative, do
not transform coordinates projecting out
translations+rotations; use Cartesian
coordiantes and no symmetry for offsets.
Otherwise, if positive, project out T+R and
use symmetry. For PBC, use a negative VIBSIZ.
Default=0.01
NVIB = 1 or 2, the number of offsets per coordinate,
similar to NVIB in $FORCE.
Default=2
STRESS = a logical flag; if set, then the numerical
gradient is computed with respect to cell
parameters (defining a stress tensor) in PBC.
Currently, only DFTB/PBC uses STRESS.
Default=false
MODGRN = Options for numerical gradient
0 Use the default general driver.
1 Use built-in FMO numerical gradient driver
(STRESS cannot be used).
Note that the use of GDDI differs in these two
methods; FMO gradients can benefit from the full
support of GDDI/2 in MODGRN=1, on par with any other
FMO run (no, GDDI/3 is not available).
In contrast, the general driver uses GDDI/2
to distribute offsets to groups.
Default: 0
==========================================================
==========================================================
$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, 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
$STATPT input group. (default is .FALSE.)
PURIFY = controls cleanup
Given a $ZMAT, the hessian and dipole derivative
tensor can be "purified" by transforming from
Cartesians to internals and back to Cartesians.
This effectively zeros the frequencies of the
translation and rotation "modes", along with
their IR intensities. The purified quantities
are punched out. Purification does change the
Hessian slightly, frequencies at a stationary
point can change by a wave number or so. The
change is bigger at non-stationary points.
(default=.FALSE. if $ZMAT is given)
PRTIFC = prints the internal coordinate force constants.
You MUST have provided $ZMAT input to use this.
(Default=.FALSE.)
HSSTYP selects the integral package used for integral
derivatives in analytic hessians. This is
used only for debugging purposes.
= BEST use experimental routines, not fully
tested.
= GENERAL use slower safe code (default).
--- 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.
= 1 Move one VIBSIZ unit in each positive
Cartesian direction. This requires 3N+1
evaluations of the wavefunction, energy, and
gradient, where N is the number of SYMMETRY
UNIQUE atoms given in $DATA.
= 2 Move one VIBSIZ unit in the positive direction
and one VIBSIZ unit in the negative direction.
This requires 6N+1 evaluations of the
wavefunction and gradient, and gives a small
improvement in accuracy. In particular, the
frequencies will change from NVIB=1 results by
no more than 10-100 wavenumbers, and usually
much less. However, the normal modes will be
more nearly symmetry adapted, and the residual
rotational and translational "frequencies"
will be much closer to zero. (default)
VIBSIZ = Displacement size (in Bohrs). This pertains to
both SEMINUM and FULLNUM. The second value of
VIBSIZ may be specified for Raman with PAVE,
with the meaning of EFIELD in $RAMAN.
Default=0.01
Let 0 mean the Vib0 geometry, and
D mean all the displaced geometries
NPRT = 1 Print orbitals at 0 and D
= 0 Print orbitals at 0 only (default)
NPUN = 2 Punch all orbitals at 0 and D
= 1 Punch all orbitals at 0 and occupied orbs at D
= 0 Punch all orbitals at 0 only (default)
NPRHSS = an option to reduce I/O for RUNTYP=FMOHESS
1 - stop printing the Hessian and dipole derivatives
2 - stop punching the Hessian and dipole derivatives
4 - stop printing/punching normal modes
16 - in PAVE, punch full dimension Hessian
default: 0
NZPHA = partition analysis for RUNTYP=FMOHESS
-2 do not do it
-1 define the number of normals modes to skip
automatically, i.e., ignore all complex-valued
frequencies.
>=0 the number of low normal modes to skip
(complex frequencies are always skept. For example,
if there are 4 complex values and NZPHA=3 then 4 are
skept; if NZPHA=6, then 6 are skept, 4 complex and 2
real).
Common values are -2, -1, 0, 3 or 6.
(3 is chosen to skip translational modes (usually lowest),
6 is for translational + rotational).
Default: -2
JACTAT = atoms in the partition analysis for RUNTYP=FMOHESS
The format of JACTAT is the same as IACTAT in STATPT.
VDOSPA = partial vibrational density for RUNTYP=FMOHESS.
If set to a non-zero value, this density is
printed discretized by the VDOSPAR value.
----- 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.
FREQ = an array of vibrational frequencies. If the
frequencies are given here, the hessian matrix
is not computed or read. You enter any imaginary
frequencies as negative numbers, omit the
zero frequencies corresponding to translation
and rotation, and enter all true vibrational
frequencies. Thermodynamic properties will be
printed, nothing else is done by the run.
PRTSCN = flag to print contribution of each vibrational
mode to the entropy. (Default is .FALSE.)
DECOMP = activates internal coordinate analysis.
Vibrational frequencies will be decomposed into
"intrinsic frequencies", by the method of
J.A.Boatz and M.S.Gordon, J.Phys.Chem., 93,
1819-1826(1989). If set .TRUE., the $ZMAT 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
==========================================================
$CPHF group (relevant for analytic RUNTYP=HESSIAN)
This group controls the solution of the response
equations, also known as coupled perturbed 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)
NCPHF = maximum number of iteractions in CPHF (currently,
applies to restricted methods only).
Note that DFTB uses MXCPIT in $DFTB instead of NCPHF.
(default=50)
TOLCP = tolerance for convergence in CPHF (currently,
applies to restricted methods only).
Note that DFTB uses CPCONV in $DFTB instead of TOLCP.
(default=5.0D-05)
NSCZV = maximum numbers of iteractions for SCZV in FMO:
NSCZV(1) refers to global iteractions over fragments
NSCZV(2) refers to local iteractions in a fragment
(default=50,50)
TOSCZV = tolerances for convergence for SCZV in FMO
TOSCZV(1) refers to global iteractions over fragments
TOSCZV(2) refers to local iteractions in a fragment
(default=1.0D-06,1.0D-08)
==========================================================
==========================================================
$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 = maximum iterations. (default=300)
CPTOL = accuracy tolerance for cpmchf equations Ax=b.
(compared to r/||b||, within orbital, CI, and
state averaged components) (default=1.0D-07)
PRCCHG = a flag to adjust the linear equation's
preconditioning. (Default=.FALSE.)
For ORMAS runs in particular, the standard
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)
NACST = 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)
--- 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.)
GCRO = a flag to select Generalized Conjugate Residual
with inner Orthogonalization. (default=.FALSE.)
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)
KICIT = 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, 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)
IPDIR = 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)
ITERB = the number of GCRO iterations to apply
before alternation to the (linear) CG converger.
(default=5)
==========================================================
==========================================================
$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.
==========================================================
==========================================================
$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.)
Formatted gradient vector at the $DATA geometry. This
data is read in the same format it was punched out.
For RUNTYP=HESSIAN, this information is used to
determine if you are at a stationary point, and possibly
for projection. If omitted, the program pretends the
gradient is zero, and otherwise proceeds normally.
For geometry searches, this information (if known) can
be read into the program so that the first step can be
taken instantly.
==========================================================
==========================================================
$DIPDR group (relevant for RUNTYP=HESSIAN if RDHESS=.T.)
Formatted dipole derivative tensor, punched in a
previous RUNTYP=HESSIAN job. If this group is omitted,
then a vibrational analysis will be unable to predict the
IR intensities, but the run can otherwise proceed.
==========================================================
$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.
==========================================================
==========================================================
$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.
===========================================================
==========================================================
$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.
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 Section 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.
= 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)
AMP = 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)
STPSZ = 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)
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.
* * * * *
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.
Selection of VDPT turns VCI on, as well.
(default=.FALSE.)
VCI = 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.)
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
IEXC2 = 0 skip combination bands (default)
= 1 add one additional quanta in other modes
= 2 add two other quanta in one mode at a time.
IEXC IEXC2 for H2O, which has only three modes:
0 0 only 000 ground state, no transitions
1 0 000, and 100, 010, 001 (fundamentals)
2 0 000, and 200, 020, 002 (1st overtones)
3 0 000, and 300, 030, 003 (2nd overtones)
1 1 000, and 100, 010, 001, 110, 101, 110
(1st overtones and combinations)
1 2 000, and 100, 010, 001, 210, 201, 021
2 1 000, and 200, 020, 002, 120, 102, 012
between them, 1st and 2nd overtones,
and all 2-1-0 combinations.
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.
VCFCT = scaling factor for pair-coupling potential.
Sometimes when pair-coupling potential values
are larger than the corresponding single mode
values, they must be scaled down. 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 MP2-
level 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.)
MPDIP = If .TRUE. the run will compute MP2 level dipoles
for the IR intensity evaluation.
Entering .FALSE. uses SCF level dipoles instead.
Default=.TRUE. for MP2 runs, except when using the
RI-MP2 program, which cannot compute MP2 dipoles,
and so chooses .FALSE. here.
It is more accurate to use the DMDR flag instead
instead of turning off MPDIP, if an MP2 level
$DIPDR is available from the MP2 hessian run.
* * * *
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.
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.
==========================================================
==========================================================
$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 group required if NFFLVL=2 or 3 in $CONTRL
$GLOWT group 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.
==========================================================
==========================================================
$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 cutoff for the norm of the mass-weighted gradient
tangent (the default is chosen in the range from
0.00005 to 0.00020, depending on the value for
STRIDE chosen below.
RCUT cutoff for Cartesian RMS displacement vector.
(the default is chosen in the range 0.0005 to
0.0020 Bohr, depending on the value for STRIDE)
ACUT maximum angle from end points for linear
interpolation (default=5 degrees)
MXOPT maximum number of constrained optimization steps
for each IRC point (default=20)
IHUPD is the hessian update formula. 1 means Powell,
2 means BFGS (default=2)
GA is a gradient from the previous IRC point, and is
used when restarting.
OPTTOL is a gradient cutoff used to determine if the IRC
is approaching a minimum. It has the same meaning
as the variable in $STATPT. (default=0.0001)
PACE = LINEAR selects linear gradient following (Euler's
method). Related input is:
STABLZ switches on Ishida/Morokuma/Komornicki reaction
path stabilization. The default is .TRUE.
DELTA initial step size along the unit bisector, if
STABLZ is on. Default=0.025 Bohr.
ELBOW is the collinearity threshold above which the
stabilization is skipped. If the mass weighted
gradients at QB and QC are almost collinear, the
reaction path is deemed to be curving very little,
and stabilization isn't needed. The default is
175.0 degrees. To always perform stabilization,
input 180.0.
READQB,EB,GBNORM,GB are energy and gradient data
already known at the current IRC point. If it
happens that a run with STABLZ on decides to skip
stabilization because of ELBOW, this data will be
punched to speed the restart.
PACE = QUAD selects quadratic gradient following.
Related input is:
SAB distance to previous point on the IRC.
GA gradient vector at that historical point.
PACE = AMPC4 selects the fourth order Adams-Moulton
variable step predictor-corrector.
Related input is:
GA0,GA1,GA2 which are gradients at previous points.
PACE = RK4 selects the 4th order Runge-Kutta variable
step method. There is no related input.
----- The next two are used by all PACE choices -----
STRIDE = Determines how far apart points on the reaction
path will be. STRIDE is used to calculate the
step taken, according to the PACE you choose.
The default is good for the GS2 method, which is
very robust. Other methods should request much
smaller step sizes, such as 0.10 or even 0.05.
(default = 0.30 sqrt(amu)-Bohr)
NPOINT = The number of IRC points to be located in this
run. The default is to find only the next point.
(default = 1)
----- 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 Print INTR at all, orbitals at all IRC points
0 Print INTR at all, orbitals at F+L (default)
-1 Print INTR at all, orbitals never
-2 Print INTR at F+L, orbitals never
NPUN = 1 Punch all orbitals at all IRC points
0 Punch all orbitals at F+L, only occupied
orbitals at IRC points between (default)
-1 Punch all orbitals at F+L only
-2 Never punch orbitals
----- The next two tally the reaction path results. The
defaults are appropriate for starting from a saddle
point, restart values are automatically punched out.
NEXTPT = The number of the next point to be computed.
STOTAL = Total distance along the reaction path to next
IRC point, in mass weighted Cartesian space.
----- The following controls jumping off the saddle point.
If you give $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
GS2 method, which must have a hessian input, even for
restarts. Note also that EVIB is ignored by GS2 runs.
* * * * * * * * * * * * * * * * * *
For hints about IRC tracking, see
the 'further information' section.
* * * * * * * * * * * * * * * * * *
==========================================================
==========================================================
$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)
DELTAT = is the time step. (default = 0.1 fs)
TOTIME = total duration of the DRC computed in a previous
job, in fs. The default is the correct value
when initiating a DRC. (default=0.0 fs)
* * *
In general, a DRC can be initiated anywhere,
so $DATA might contain coordinates of the
equilibrium geometry, or a nearby transition
state, or something else. You must also
supply an initial kinetic energy, and the
direction of the initial velocity, for which
there are a number of options:
EKIN = The initial kinetic energy (default = 0.0
kcal/mol)
See also ENM, NVEL, and VIBLVL regarding alternate
ways to specify the initial value.
VEL = an array of velocity components, in Bohr/fs.
When NVEL is false, this is simply the direction
of the velocity vector. Its magnitude will be
automatically adjusted to match the desired
initial
kinetic energy, and it will be projected so that
the translation of the center of mass is removed.
Give in the order vx1, vy1, vz1, vx2, vy2, ...
NVEL = a flag to compute the initial kinetic energy from
the input VEL using the sum of mass*VEL*VEL/2.
This flag is usually selected only for restarts.
(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.
* * *
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)
ENM = the initial kinetic energy given to mode NNM,
in units of vibrational quanta hv, so the amount
depends on mode NNM's vibrational frequency, v.
If you prefer to impart an arbitrary initial
kinetic energy to mode NNM, specify EKIN instead.
(default = 0.0 quanta)
To summarize, there are 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
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 = MIN indicates the hessian supplied for the initial
geometry corresponds to a minimum (default).
= TS indicates the hessian is for a saddle point.
HESS2 = 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
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 Print orbitals at all DRC points
0 Print orbitals at F+L (default)
-1 Never print orbitals
NPUN = 2 Punch all orbitals at all DRC points
1 Punch all orbitals at F+L, and occupied
orbitals at DRC points between
0 Punch all orbitals at F+L only (default)
-1 Never punch orbitals
==========================================================
===========================================================
$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 spin-
orbit 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 = maximum number of search steps (default=50)
STPSZ = 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)
NMOS1 = Number of orbitals for first state's $VEC1.
NMOS2 = Number of orbitals for second state's $VEC2.
NPRT = 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 = maximum of effective gradient between the two
states (default = 5.0D-04)
TGRMS = 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).
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)
===========================================================
$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 (three)
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
= PENALTT penalty-constrained optimization method for
three state conical intersection search
See J. Phys.Chem.A 125.9 (2021)
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.
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!
ITROOT = array of three states whose CI point is sought.
For example, this might be ITROOT(1)=2,3,4 The
roots are counted exactly the same as IROOT in
the $DET or $TDDFT input groups. There is no
default for ITROOT. Only OPTTYP=PENALTY is
supported.
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 = parameter ensuring a singularity free penalty,
default=0.02 Hartree
SIGMA = 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
T.W.Keal, A.Koslowski, W.Thiel
Theoret.Chem.Acc. 118, 837(2007)
===========================================================
===========================================================
$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.
DT = MD time step size, in seconds, default=1.0d-15,
which is a femtosecond.
NVTNH selects a integrator step appropriate to the
desired ensemble. This is only implemented for
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. In REMD, BATHT is an array with the
size given by the number of replicas.
* * *
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. default=.FALSE.
DTEMP = 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.
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.
DELR = spacing for radial bins in RDF calculations,
default=0.02 Angstroms.
NPROP = 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.
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).
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
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.