*********************************
*
*
* 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. There are
also a number
of examples shown in the 'Input Examples'
section.
It is useful to note that this chapter of the manual
can be searched
online by means of the "gmshelp" command,
if your computer
is of the Unix type. A command such as
"gmshelp scf"
will display the $SCF input group. With
no arguments,
the gmshelp command will show you all input
group names.
Type "q" to exit the pager, and note that
some pagers
will let you back up by means of "b".
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). After that comes run
type related
input, then properties input, input for two
different solvation
models, integral related input, and
finally CI/MCSCF
input. The next page contains a list of
all possible
input groups, in the order in which they can
be found in
this section.
1
*
name
function
module:routine
----
--------
--------------
Molecule, basis, wavefunction specification:
$CONTRL
chemical control data
INPUTA:START
$SYSTEM
computer related control data INPUTA:START
$BASIS
basis set
INPUTB:BASISS
$DATA
molecule, basis set
INPUTB:MOLE
$ZMAT
coded z-matrix
ZMATRX:ZMATIN
$LIBE
linear bend data
ZMATRX:LIBE
$SCF
HF-SCF wavefunction control SCFLIB:SCFIN
$SCFMI
SCF-MI input control data
SCFMI :MIINP
$DFT
density functional input
DFT :DFTINP
$MP2
2nd order Moller-Plesset
MP2 :MP2INP
$GUESS
initial orbital selection
GUESS :GUESMO
$VEC
orbitals
(formatted) GUESS :READMO
$MOFRZ
freezes MOs during SCF runs EFPCOV:MFRZIN
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
$HESS
force constant matrix (formatted) HESS :FCMIN
$GRAD
gradient vector (formatted) HESS
:EGIN
$DIPDR
dipole deriv. matrix (formatted) HESS :DDMIN
$VIB
HESSIAN restart data (formatted) HESS :HSSNUM
$MASS
isotope selection
VIBANL:RAMS
$IRC
intrinsic reaction path
RXNCRD:IRCX
$VSCF
vibrational SCF and MP2
VSCF :VSCFIN
$VIBSCF
VSCF restart data (formatted) VSCF :VGRID
$DRC
dynamic reaction path
DRC :DRCDRV
$GLOBOP
Monte Carlo global fragment opt GLOBOP:GLOPDR
$GRADEX
gradient extremal path
GRADEX:GRXSET
$SURF
potential surface scan
SURF :SRFINP
continued on the next page...
1
*
name
function
module:routine
----
--------
--------------
Interpretation, properties:
$LOCAL
orbital localization control LOCAL :LMOINP
$TWOEI
J,K integrals (formatted)
LOCCD :TWEIIN
$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
$MOLGRF
orbital plots
PARLEY:PLTMEM
$STONE
distributed multipole analysis PRPPOP:STNRD
$RAMAN
Raman intensity
RAMAN :RAMANX
$ALPDR
alpha polar. der. (formatted) RAMAN :ADMIN
$MOROKM
Morokuma energy decomposition MOROKM:MOROIN
$FFCALC
finite field polarizabilities FFIELD:FFLDX
$TDHF
time dependent HF NLO properties TDHF :TDHFX
Solvation models:
$EFRAG
effective fragment potentials EFINP :EFINP
$FRAGNAME specific
named fragment pot. EFINP :RDSTFR
$FRGRPL
inter-fragment repulsion
EFINP :RDDFRL
$PCM
polarizable continuum model PCM
:PCMINP
$PCMCAV
PCM cavity generation
PCM :MAKCAV
$NEWCAV
PCM escaped charge cavity
PCM :DISREP
$IEFPCM
PCM integral equation form. data PCM :IEFDAT
$DISBS
PCM dispersion basis set
PCMDIS:ENLBS
$DISREP
PCM dispersion/repulsion
PCMVCH:MORETS
$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
$RELWFN
relativistic correction
INPUTB:RWFINP
$EFIELD
external electric field
PRPLIB:INPEF
$INTGRL
format for 2e- integrals
INT2A :INTIN
$FMM
fast multipole method
QMFM :QFMMIN
$TRANS
integral transformation
TRANS :TRFIN
continued on the next page...
1
*
name
function
module:routine
----
--------
--------------
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
$GCILST
general determinant list
ALGNCI:GCIGEN
$DRT
distinct row table for MCSCF GUGDRT:ORDORB
$CIDRT
distinct row table for CI
GUGDRT:ORDORB
$MCSCF
parameters for MCSCF
MCSCF :MCSCF
$MCQDPT
multireference pert. theory MCQDPT:MQREAD
$CISORT
integral sorting
GUGSRT:GUGSRT
$GUGEM
Hamiltonian matrix formation GUGEM :GUGAEM
$GUGDIA
Hamiltonian eigenvalues/vectors GUGDGA:GUGADG
$GUGDM
1e- density matrix
GUGDM :GUGADM
$GUGDM2
2e- density matrix
GUGDM2:GUG2DM
$LAGRAN
CI lagrangian matrix
LAGRAN:CILGRN
$TRFDM2
2e- density backtransformation TRFDM2:TRF2DM
$TRANST
transition moments, spin-orbit TRNSTN:TRNSTX
* this column is more useful to programmers than to users.
1
$CONTRL
==========================================================
$CONTRL group
(optional)
This is a free format group specifying global switches.
SCFTYP
together with MPLEVL or CITYP specifies
the wavefunction. You may choose from
= RHF Restricted Hartree Fock calculation
(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 OCBSE type ROHF. (needs $SCF input)
= MCSCF Multiconfigurational SCF wavefunction
(this requires $DET or $DRT input)
= NONE indicates a single point computation,
rereading a converged SCF function.
This option requires that you select
CITYP=GUGA, ALDET, or GENCI with
only RUNTYP=ENERGY or TRANSITN
and with GUESS=MOREAD.
MPLEVL =
chooses Moller-Plesset perturbation
theory level, after the SCF.
See $MP2 and $MCQDPT input groups.
= 0 skips the MP computation
(default)
= 2 performs a second order
energy
correction. MP2 is implemented only
for RHF, UHF, ROHF, and MCSCF wave
functions. Gradients are available
only for RHF, so for the others you
may pick from RUNTYP=ENERGY, TRUDGE,
SURFACE, or FFIELD only.
1
CITYP =
chooses CI computation after the SCF,
for any SCFTYP except UHF.
= NONE skips the CI. (default)
= GUGA runs the Unitary Group CI package,
which requires $CIDRT input.
Gradients are available only for RHF,
so for other SCFTYPs, you may choose
only RUNTYP=ENERGY, TRUDGE, SURFACE,
FFIELD, TRANSITN.
= ALDET runs the Ames Laboratory determinant
full CI package, requiring $CIDET
input. RUNTYP=ENERGY only.
= GENCI runs a determinant CI program that
permits arbitrary specification of
the determinants, requiring $CIGEN
input. RUNTYP=ENERGY only.
Obviously, at
most one of MPLEVL or CITYP may be chosen.
Likewise, you
should not mix either of these with DFT.
1
$CONTRL
RUNTYP
specifies the type of computation, for
example at a single geometry point:
= ENERGY Molecular energy. (default)
= GRADIENT Molecular energy plus gradient.
= HESSIAN Molecular energy plus gradient plus
second derivatives, including harmonic
harmonic vibrational analysis. See the
$FORCE and $CPHF input groups.
multiple geometry options:
= OPTIMIZE Optimize the molecular geometry using
analytic energy gradients. See $STATPT.
= TRUDGE Non-gradient total energy minimization.
See groups $TRUDGE and $TRURST.
= SADPOINT Locate saddle point (transition state).
See the $STATPT group.
= IRC Follow intrinsic reaction coordinate.
See the $IRC group.
= VSCF Compute anharmonic vibrational
corrections (see $VSCF)
= DRC Follow dynamic reaction coordinate.
See the $DRC group.
= GLOBOP global optimization of effective fragment
positions via Monte Carlo. See $GLOBOP.
= GRADEXTR Trace gradient extremal.
See the $GRADEX group.
= SURFACE Scan linear cross sections of the
potential energy surface. See $SURF.
single geometry property options:
= PROP Properties will be calculated.
A $DATA
deck and converged $VEC group should be
input. Optionally, orbital localization
can be done. See $ELPOT, etc.
= RAMAN computes Raman intensities, see $RAMAN.
= MOROKUMA Performs monomer energy decomposition.
See the $MOROKM group.
= TRANSITN Compute radiative transition moment or
spin-orbit coupling. See $TRANST group.
= FFIELD applies finite electric fields, most
commonly to extract polarizabilities.
See the $FFCALC group.
= TDHF analytic computation of time dependent
polarizabilities. See the $TDHF group.
= MAKEFP creates an effective fragment potential.
* * * * * * * * * * * * * * * * * * * * * * * * * *
Note that RUNTYPs involving the energy gradient,
namely GRADIENT, HESSIAN, OPTIMIZE, SADPOINT,
GLOBOP, IRC, GRADEXTR, and DRC, cannot be used for
any CI or MP2 computation, except when SCFTYP=RHF.
* * * * * * * * * * * * * * * * * * * * * * * * * *
1
$CONTRL
EXETYP = RUN
Actually do the run. (default)
= CHECK Wavefunction and energy will not be
evaluated. This lets you speedily
check input and memory requirements.
See the overview section for details.
= DEBUG Massive amounts of output are printed,
useful only if you hate trees.
= routine Maximum output is generated by the
routine named. Check the source for
the routines this applies to.
MAXIT =
Maximum number of SCF iteration cycles.
Pertains only to RHF, UHF, ROHF, or
GVB runs. See also MAXIT in $MCSCF.
(default = 30)
* * * * * * *
ICHARG = Molecular charge. (default=0, neutral)
MULT
= Multiplicity
of the electronic state
= 1 singlet (default)
= 2,3,... doublet, triplet, and so on.
ICHARG and MULT are used directly for RHF, UHF, ROHF.
For GVB, these are implicit in the $SCF input, while
for MCSCF or CI, these are implicit in $DRT/$CIDRT or
$DET/$CIDET input. You must still give them correctly.
* * * * * * *
ECP
= effective
core potential control.
= NONE all electron calculation (default).
= READ read the potentials in $ECP group.
= SBKJC use Stevens, Basch, Krauss, Jasien,
Cundari potentials for all heavy
atoms (Li-Rn are available).
= HW use Hay, Wadt potentials
for all the
heavy atoms (Na-Xe are available).
* * * * * * *
RELWFN = NONE
(default) See also $RELWFN input group.
= NESC normalised elimination of small component,
the method of K. Dyall
= RESC relativistic elimination of small component,
the method of T. Nakajima and K. Hirao.
1
$CONTRL
* * * the next three control molecular geometry * * *
COORD =
choice for molecular geometry in $DATA.
= UNIQUE only the symmetry unique atoms will be
given, in Cartesian coords (default).
= HINT only the symmetry unique atoms will
be
given, in Hilderbrandt style internals.
= CART Cartesian coordinates will be input.
Please read the warning just below!!!
= ZMT GAUSSIAN style internals will
be input.
= ZMTMPC MOPAC style internals will be input.
= FRAGONLY means no part of the system is treated
by ab initio means, hence $DATA is not
given. The system is specified by $EFRAG.
Note
that the CART, ZMT, ZMTMPC choices require input of
all atoms in the molecule. These three also orient the
molecule, and then determine which atoms are unique. The
reorientation is very likely to change the order of the
atoms from what you input. When the point group contains
a 3-fold or higher rotation axis, the degenerate moments
of inertia often cause problems choosing correct symmetry
unique axes, in which case you must use COORD=UNIQUE
rather than Z-matrices.
Warning:
The reorientation into principal axes is done
only for atomic coordinates, and is not applied to the
axis dependent data in the following groups: $VEC, $HESS,
$GRAD, $DIPDR, $VIB, nor Cartesian coords of effective
fragments in $EFRAG. COORD=UNIQUE avoids reorientation,
and thus is the safest way to read these.
Note
that the choices CART, ZMT, ZMTMPC require the use
of a $BASIS group to define the basis set. The first
two choices might or might not use $BASIS, as you wish.
UNITS =
distance units, any angles must be in degrees.
= ANGS Angstroms (default)
= BOHR Bohr atomic units
NZVAR =
0 Use Cartesian coordinates (default).
= M If COORD=ZMT or ZMTMPC and a $ZMAT is not given:
the internal coordinates will be those defining
the molecule in $DATA. In this case, $DATA must
not contain any dummy atoms. M is usually 3N-6,
or 3N-5 for linear.
= M For other COORD choices, or if $ZMAT is given:
the internal coordinates will be those defined
in $ZMAT. This allows more sophisticated
internal coordinate choices. M is ordinarily
3N-6 (3N-5), unless $ZMAT has linear bends.
NZVAR
refers mainly to the coordinates used by OPTIMIZE
or SADPOINT runs, but may also print the internal's
values for other run types. You can use internals to
define the molecule, but Cartesians during optimizations!
1
$CONTRL
LOCAL =
controls orbital localization.
= NONE Skip localization (default).
= BOYS Do Foster-Boys localization.
= RUEDNBRG Do Edmiston-Ruedenberg localization.
= POP Do Pipek-Mezey population localization.
See the $LOCAL group. Localization
does not work for SCFTYP=GVB or CITYP.
ISPHER =
Spherical Harmonics option
= -1 Use Cartesian basis functions to construct
symmetry-adapted linear combination (SALC)
of basis functions. The SALC space is the
linear variation space used. (default)
= 0 Use spherical harmonic functions to create
SALC functions, which are then expressed
in terms of Cartesian functions. The
contaminants are not dropped, hence this
option has EXACTLY the same variational
space as ISPHER=-1. The only benefit to
obtain from this is a population analysis
in terms of pure s,p,d,f,g functions.
= +1 Same as ISPHER=0, but the function space
is truncated to eliminate all contaminant
Cartesian functions [3S(D), 3P(F), 4S(G),
and 3D(G)] before constructing the SALC
functions. The computation corresponds
to the use of a spherical harmonic basis.
QMTTOL = linear
dependence threshhold
Any functions in the SALC variational space whose
eigenvalue of the overlap matrix is below this
tolerence is considered to be linearly dependent.
Such functions are dropped from the variational
space. What is dropped is not individual basis
functions, but rather some linear combination(s)
of the entire basis set that represent the linear
dependent part of the function space. The default
is a reasonable value for most purposes, 1.0E-6.
When many diffuse functions are used, it is common
to see the program drop some combinations. On
occasion, in multi-ring molecules, we have raised
QMTTOL to 3.0E-6 to obtain SCF convergence, at the
cost of some energy.
1
$CONTRL
* * * 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
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.
The NBO program of Frank Weinhold's group can be
attached to
GAMESS. The input to control the natural
bond order analysis
is read by the add in code, so is
not described
here. The NBO program is available by
anonymous FTP
to ftp.osc.edu, in the directory
pub/chemistry/software/SOURCES/FORTRAN/nbo
1
$CONTRL
* * * computation control switches * * *
For the most part, the default is the only sensible
value, and unless
you are sure of what you are doing,
these probably
should not be touched.
NPRINT =
Print/punch control flag
See also EXETYP for debug info.
(options -7 to 5 are primarily debug)
= -7 Extra printing from Boys
localization.
= -6 debug for geometry searches
= -5 minimal output
= -4 print 2e-contribution to
gradient.
= -3 print 1e-contribution to
gradient.
= -2 normal printing, no punch
file
= 1 extra printing for
basis,symmetry,ZMAT
= 2 extra printing for
MO guess routines
= 3 print out property
and 1e- integrals
= 4 print out 2e- integrals
= 5 print out SCF data
for each cycle.
(Fock and density matrices, current MOs
= 6 same as 7, but wider
132 columns output.
This option isn't perfect.
= 7 normal printing and
punching (default)
= 8 more printout than
7. The extra output
is (AO) Mulliken and overlap population
analysis, eigenvalues, Lagrangians, ...
= 9 everything in 8 plus
Lowdin population
analysis, final density matrix.
NOSYM =
0 the symmetry specified in $DATA is used
as much as possible in integrals, SCF,
gradients, etc. (this is the default)
= 1 the symmetry specified in the $DATA group
is used to build the molecule, then
symmetry is not used again. Some GVB
or MCSCF runs (those without a totally
symmetric charge density) require you
request no symmetry.
INTTYP = POPLE
use fast Pople-Hehre routines for sp integral
blocks, and HONDO Rys polynomial code for
all other integrals. (default)
= HONDO use HONDO/Rys integrals for all integrals.
This option produces slightly more accurate
integrals but is also slower.
When diffuse functions are used, the inaccuracy in
Pople/Hehre sp integrals shows up as inaccurate
LCAO coefficients in virtual orbitals. This means
the error in SCF (meaning RHF to MCSCF) energies is
expected to be about 5d-8 Hartree, but the error in
computations that OCCUPY the virtual orbitals may
be much larger. We have seen an energy error of
1d-4 in an MP2 energy when diffuse functions were
used. We recommend that all MP2 or CI jobs with
diffuse functions select INTTYP=HONDO.
1
NORMF =
0 normalize the basis functions (default)
= 1 no normalization
NORMP =
0 input contraction coefficients refer to
normalized Gaussian primitives. (default)
= 1 the opposite.
$CONTRL
ITOL
= primitive cutoff factor (default=20)
= n products of primitives whose exponential
factor is less than 10**(-n) are skipped.
ICUT
= n integrals less than 10.0**(-n) are not
saved on disk. (default = 9)
* * * restart options * * *
IREST =
restart control options
(for OPTIMIZE run restarts, see $STATPT)
Note that this option is unreliable!
= -1 reuse dictionary file from previous run,
useful with GEOM=DAF and/or GUESS=MOSAVED.
Otherwise, this option is the same as 0.
= 0 normal run (default)
= 1 2e restart (1-e integrals and MOs saved)
= 2 SCF restart (1-,2-e integrls and MOs saved)
= 3 1e gradient restart
= 4 2e gradient restart
GEOM
= select where to obtain molecular
geometry
= INPUT from $DATA input (default for IREST=0)
= DAF read from DICTNRY file (default otherwise)
As noted in the first chapter, binary file restart is
not a well tested
option!
==========================================================
1
$SYSTEM
==========================================================
$SYSTEM group (optional)
This group provides global control information for
your computer's
operation. This is system related input,
and will not
seem particularly chemical to you!
TIMLIM =
time limit, in minutes. Set to about 95 percent
of the time limit given to the batch job so that
GAMESS can stop itself gently. (default=600.0)
MWORDS =
the maximum replicated memory which your job can
use, on every node. This is given in units of
1,000,000 words (as opposed to 1024*1024 words),
where a word is always a 64 bit quantity. Most
systems allocate this memory at run time, but
some more primitive systems may have an upper
limit chosen at compile time. (default=1)
In case finer control over the memory is needed,
this value can be given in units of words by
using the keyword MEMORY instead of MWORDS.
MEMDDI =
the grand total memory needed for the distributed
data interface (DDI) storage, given in units of
1,000,000 words. See Chapter 5 of this manual for
an extended explanation of running with MEMDDI.
note: the memory
required on each node for a run using
p processors is therefore MWORDS + MEMDDI/p.
PARALL =
a flag to cause the distributed data parallel
MP2 program to execute the parallel algorithm
even if you are running on only one node.
The main purpose of this is to allow you to
do EXETYP=CHECK runs to learn what the correct
value of MEMDDI needs to be.
KDIAG =
diagonalization control switch
= 0 use a vectorized diagonalization routine
if one is available on your machine,
else use EVVRSP. (default)
= 1 use EVVRSP diagonalization. This may
be more accurate than KDIAG=0.
= 2 use GIVEIS diagonalization
(not as fast or reliable as EVVRSP)
= 3 use JACOBI diagonalization
(this is the slowest method)
1
COREFL =
a flag to indicate whether or not GAMESS
should produce a "core" file for debugging
when subroutine ABRT is called to kill
a job. This variable pertains only to
UNIX operating systems. (default=.FALSE.)
* * * the next three refer to parallel GAMESS * * *
The next three
apply only to parallel runs, and as they
are more or
less obsolete, their use is discourged.
BALTYP =
Parallel load balence scheme
LOOP turns off dynamic load balancing (DLB)
NXTVAL uses dynamic load balancing
(default = LOOP)
XDR
= a flag to indicate whether or not messages
should be converted into a generic format
known as external data representation.
If true, messages can exchange between
machines of different vendors, at the cost
of performing the data type conversions.
(default=.FALSE.) --inactive at present--
PTIME =
a logical flag to print extra timing info
during parallel runs. This is not currently
implemented.
==========================================================
1
$BASIS
==========================================================
$BASIS group (optional)
This group allows certain standard basis sets to be
easily given.
If this group is omitted, the basis set
must be given
instead in the $DATA group.
GBASIS =
Name of the Gaussian basis set.
= MINI - Huzinaga's 3 gaussian minimal basis set.
Available H-Rn.
= MIDI - Huzinaga's 21 split valence basis set.
Available H-Rn.
= STO - Pople's STO-NG minimal basis set.
Available H-Xe, for NGAUSS=2,3,4,5,6.
= N21 - Pople's N-21G split valence basis set.
Available H-Xe, for NGAUSS=3.
Available H-Ar, for NGAUSS=6.
= N31 - Pople's N-31G split valence basis set.
Available H-Ne,P-Cl for NGAUSS=4.
Available H-He,C-F for NGAUSS=5.
Available H-Ar, for NGAUSS=6.
For Ga-Kr, N31 selects the BC basis.
= N311 - Pople's "triple split" N-311G basis set.
Available H-Ne, for NGAUSS=6.
Selecting N311 implies MC for Na-Ar.
= DZV - "double zeta valence" basis set.
a synonym for DH for H,Li,Be-Ne,Al-Cl.
(14s,9p,3d)/[5s,3p,1d] for K-Ca.
(14s,11p,5d/[6s,4p,1d] for Ga-Kr.
= DH - Dunning/Hay "double zeta" basis set.
(3s)/[2s] for H.
(9s,4p)/[3s,2p] for Li.
(9s,5p)/[3s,2p] for Be-Ne.
(11s,7p)/[6s,4p] for Al-Cl.
= TZV - "triple zeta valence" basis set.
(5s)/[3s] for H.
(10s,3p)/[4s,3p] for Li.
(10s,6p)/[5s,3p] for Be-Ne.
a synonym for MC for Na-Ar.
(14s,9p)/[8s,4p] for K-Ca.
(14s,11p,6d)/[10s,8p,3d] for Sc-Zn.
= MC - McLean/Chandler "triple split" basis.
(12s,9p)/[6s,5p] for Na-Ar.
Selecting MC implies 6-311G for H-Ne.
additional values for GBASIS are on the next page.
1
* * * the next two are ECP bases only * * *
GBASIS = SBKJC-
Stevens/Basch/Krauss/Jasien/Cundari
valence basis set, for Li-Rn. This choice
implies an unscaled -31G basis for H-He.
= HW - Hay/Wadt valence basis.
This is a -21 split, available Na-Xe,
except for the transition metals.
This implies a 3-21G basis for H-Ne.
* * * semiempirical basis sets * * *
The elements for which these exist can be found
in the 'further information' section of this
manual. If you pick one of these, all other data
in this group is ignored. Semi-empirical runs
actually use valence-only STO bases, not GTOs.
GBASIS = MNDO - selects MNDO model hamiltonian
= AM1 - selects AM1 model hamiltonian
= PM3 - selects PM3 model hamiltonian
NGAUSS = the
number of Gaussians (N). This parameter
pertains only to GBASIS=STO, N21, N31, or N311.
NDFUNC = number
of heavy atom polarization functions to
be used. These are usually d functions, except
for MINI/MIDI. The term "heavy" means Na on up
when GBASIS=STO, HW, or N21, and from Li on up
otherwise. The value may not exceed 3. The
variable POLAR selects the actual exponents to
be used, see also SPLIT2 and SPLIT3. (default=0)
NFFUNC = number
of heavy atom f type polarization
functions to be used on Li-Cl. This may only
be input as 0 or 1. (default=0)
NPFUNC = number
of light atom, p type polarization
functions to be used on H-He. This may not
exceed 3, see also POLAR. (default=0)
DIFFSP = flag
to add diffuse sp (L) shell to heavy atoms.
Heavy means Li-F, Na-Cl, Ga-Br, In-I, Tl-At.
The default is .FALSE.
DIFFS =
flag to add diffuse s shell to hydrogens.
The default is .FALSE.
Warning: if you
use diffuse functions, please read QMTTOL
and INTTYP in
the $CONTRL group for numerical concerns.
1
$BASIS
POLAR =
exponent of polarization functions
= POPLE (default for all other cases)
= POPN311 (default for GBASIS=N311, MC)
= DUNNING (default for GBASIS=DH, DZV)
= HUZINAGA (default for GBASIS=MINI, MIDI)
= HONDO7 (default for GBASIS=TZV)
SPLIT2 = an array
of splitting factors used when NDFUNC
or NPFUNC is 2. Default=2.0,0.5
SPLIT3 = an array
of splitting factors used when NDFUNC
or NPFUNC is 3. Default=4.00,1.00,0.25
EXTFIL = a flag
to read basis sets from an external file,
defined by EXTBAS, instead of $DATA.
No external file is provided with GAMESS, instead
you would supply your own. The GBASIS keyword
must give an 8 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 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 GBASIS string allows
you to have several families of bases in the same
file, identified by different strings.
(default=.false.)
==========================================================
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 group. Thus,
if $BASIS does
not quite cover the basis you want, you can
obtain this
full $DATA group from EXETYP=CHECK, and then
change polarization
exponents, add Rydbergs, etc.
1
$DATA
==========================================================
$DATA group
(required)
$DATAS group
(if NESC chosen, gives small component basis)
$DATAL group
(if NESC chosen, gives large component basis)
This group describes the global molecular data such as
point group
symmetry, nuclear coordinates, and possibly
the basis set.
It consists of a series of free format
card images.
See $RELWFN for more 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", Chapter 3.
Armed with only the name of the group, GAMESS is able
to exploit the
molecular symmetry throughout almost all of
the program,
and thus save a great deal of computer time.
GAMESS does
not require that you know very much else about
group theory,
although a deeper knowledge (character
tables, irreducible
representations, term symbols, and so
on) is useful
when dealing with the more sophisticated
wavefunctions.
1
$DATA
Cards -3- and
-4- are quite complicated, and are rarely
given.
A *SINGLE* blank card may replace both cards -3-
and -4-, to
select the 'master frame', which is defined on
the next page.
If you choose to enter a blank card, skip
to the bottom
of the next page.
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.
----------------------------------------------------------
1
$DATA
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.
1
$DATA
----------------------------------------------------------
-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-.
1
$DATA
----------------------------------------------------------
-6U- GBASIS,
NGAUSS, (SCALF(i),i=1,4)
GBASIS has exactly
the same meaning as in $BASIS. You may
choose from
MINI, MIDI, STO, N21, N31, N311, DZV, DH, BC,
TZV, MC, SBKJC,
or HW. In addition, you may choose S, P,
D, F, G, or
L to enter an explicit basis set. Here, L
means both an
s and p shell with a shared exponent.
NGAUSS is the
number of Gaussians (N) in the Pople style
basis, or user
input general basis. It has meaning only
for GBASIS=STO,
N21, N31, or N311, and S,P,D,F,G, or L.
Up to four scale
factors may be entered. If omitted,
standard values
are used. They are not documented as
every GBASIS
treats these differently. Read the source
code if you
need to know more. They are seldom given.
----------------------------------------------------------
* * * If GBASIS
is not S,P,D,F,G, or L, either add more
shells by repeating card -6U-, or go on to -8U-.
* * * If GBASIS=S,P,D,F,G,
or L, enter NGAUSS cards -7U-.
----------------------------------------------------------
-7U- IG, ZETA,
C1, C2
IG = a counter, IG takes values 1, 2, ..., NGAUSS.
ZETA = Gaussian exponent of the IG'th primitive.
C1 = Contraction coefficient for S,P,D,F,G shells,
and for the s function of L shells.
C2 = Contraction coefficient for the p in L shells.
----------------------------------------------------------
* * * For more
shells on this atom, go back to card -6U-.
* * * If there
are no more shells, go on to card -8U-.
----------------------------------------------------------
-8U-
A blank card ends the basis set for this atom.
----------------------------------------------------------
Continue entering
atoms with -5U- through -8U- until all
are given, then
terminate the group with a " $END " card.
--- this is the end of card sequence U ---
1
$DATA
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 ---
1
$DATA
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!
1
$DATA
o
In -6G- to -8G-, symbolic strings may be given in
place of numeric values for BLENGTH, ALPHA, and BETA.
The same string may be repeated, which is handy in
enforcing symmetry. If the string is preceeded by a
minus sign, the numeric value which will be used is
the opposite, of course. Any mixture of numeric data
and symbols may be given. If any strings were given
in -6G- to -8G-, you must provide cards -9G- and
-10G-, otherwise you may terminate the group now with
a " $END " card.
----------------------------------------------------------
-9G- A blank line terminates the Z-matrix section.
----------------------------------------------------------
-10G- STRING VALUE
STRING is a symbolic
string used in the Z-matrix.
VALUE
is the numeric value to substitute for that string.
----------------------------------------------------------
Continue entering
-10G- until all STRINGs are defined.
Note that any
blank card encountered while reading -10G-
will be ignored.
GAMESS regards all STRINGs as variables
(constraints
are sometimes applied in $STATPT). It is not
necessary to
place constraints to preserve point group
symmetry, as
GAMESS will never lower the symmetry from
that given at
-2-. When you have given all STRINGs a
VALUE, terminate
the group with a " $END " card.
--- this is the end of card sequence G ---
* * * *
The documentation for sequence G above and sequence M
below presumes
you are reasonably familiar with the input
to GAUSSIAN
or MOPAC. It is probably too terse to be
understood very
well if you are unfamiliar with these. A
good tutorial
on both styles of Z-matrix input can be
found in Tim
Clark's book "A Handbook of Computational
Chemistry",
published by John Wiley & Sons, 1985.
Both Z-matrix input styles must generate a molecule
which possesses
the symmetry you requested at -2-. If
not, your job
will be terminated automatically.
1
$DATA
COORD=ZMTMPC input: (MOPAC style internals)
----------------------------------------------------------
-5M- ATOM
Only the name
of the first atom is required.
See -8M- for
a description of this information.
----------------------------------------------------------
-6M- ATOM BLENGTH
Only a name and
a bond distance is required for atom 2.
See -8M- for
a description of this information.
----------------------------------------------------------
-7M- ATOM BLENGTH j1 ALPHA j2
Only a bond distance
from atom 2, and an angle with repect
to atom 1 is
required for atom 3. If you prefer to hook
atom 3 to atom
1, you must give connectivity as in -8M-.
See -8M- for
a description of this information.
----------------------------------------------------------
-8M- ATOM BLENGTH j1 ALPHA j2 BETA j3 i1 i2 i3
ATOM, BLENGTH,
ALPHA, BETA, i1, i2 and i3 are as described
at -8G-.
However, BLENGTH, ALPHA, and BETA must be given
as numerical
values only. In addition, BETA is always a
dihedral angle.
i1, i2, i3 must be integers only.
The j1, j2 and
j3 integers, used in MOPAC to signal
optimization
of parameters, must be supplied but are
ignored here.
You may give them as 0, for example.
----------------------------------------------------------
Continue entering
atoms 3, 4, 5, ... with -8M- cards until
all are given,
and then terminate the group by giving a
" $END " card.
--- this is the end of card sequence M ---
==========================================================
This is the end of $DATA!
If you have any
doubt about what molecule and basis set
you are defining,
or what order the atoms will be
generated in,
simply execute an EXETYP=CHECK job to find
out!
1
$ZMAT
==========================================================
$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.
1
$ZMAT
**** the next define delocalized coordinates ****
DLC
is a flag to request delocalized coordinates.
(default is .FALSE.)
AUTO
is a flag to generate all redundant coordinates,
automatically. The DLC space will consist of all
non-redundant combinations of these which can be
found. The list of redundant coordinates will
consist of bonds, angles, and torsions only.
(default is .FALSE.)
NONVDW is an
array of atom pairs which are to be joined
by a bond, but might be skipped by the routine
that automatically includes all distances shorter
than the sum of van der Waals radii. Any angles
and torsions associated with the new bond(s) are
also automatically included.
The format for IXZMAT, IRZMAT, IFZMAT is that of IZMAT:
IXZMAT is an
extra array of simple internal coordinates
which you want to have added to the list generated
by AUTO. Unlike NONVDW, IXZMAT will add only the
coordinate(s) you specify.
IRZMAT is an
array of simple internal coordinates which
you would like to remove from the AUTO list of
redundant coordinates. It is sometimes necessary
to remove a torsion if other torsions around a bond
are being frozen, to obtain a nonsingular G matrix.
IFZMAT is an
array of simple internal coordinates which
you would like to freeze. See also FVALUE below.
Note that IFZMAT/FVALUE work only with DLC, see the
IFREEZ option in $STATPT to freeze coordinates if
you wish to freeze simple or natural coordinates.
FVALUE is an
array of values to which the internal
coordinates should be constrained. It is not
necessary to input $DATA such that the initial
values match these desired final values, but it is
helpful if the initial values are not too far away.
1
$ZMAT $LIBE
**** 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.
==========================================================
1
$SCF
==========================================================
$SCF group
relevant if SCFTYP = RHF, UHF, or ROHF,
required if SCFTYP = GVB)
This group of parameters provides additional control
over the RHF,
UHF, ROHF, or GVB SCF steps. It must be
used for GVB
open shell or perfect pairing wavefunctions.
DIRSCF = a flag
to activate a direct SCF calculation,
which is implemented for all the Hartree-Fock
type wavefunctions: RHF, ROHF, UHF, and GVB.
This keyword also selects direct MP2 computation.
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 sometimes oscillate at the end, rather than
converging. Turning this parameter off will
normally give convergence.
---- The next flags affect convergence rates.
EXTRAP = controls
Pople extrapolation of the Fock matrix.
DAMP
= controls Davidson damping of the Fock matrix.
SHIFT
= controls level shifting of the Fock matrix.
RSTRCT = controls
restriction of orbital interchanges.
DIIS
= controls Pulay's DIIS interpolation.
SOSCF
= controls second order SCF orbital optimization.
(default=.TRUE. for RHF, Abelian group ROHF, GVB)
(default=.FALSE. for UHF, non-Abelian group ROHF)
DEM
= controls direct energy minimization, which is
implemented only for RHF. (default=.FALSE.)
defaults for
EXTRAP DAMP SHIFT RSTRCT DIIS SOSCF
ab initio:
T F F
F T T/F
semiempirical:
T F F
F F F
The above parameters are implemented for all SCF
wavefunction
types, except that DIIS will work for GVB
only for those
cases with NPAIR=0 or NPAIR=1. If both
DIIS and SOSCF
are chosen, SOSCF is stronger than DIIS,
and so DIIS
will not be used.
Once either DIIS or SOSCF are initiated, any other
accelerator
in effect is put in abeyance.
1
$SCF
---- 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 or MP2 gradients use 1.0d-06.
SOGTOL = second
order gradient tolerance. SOSCF will be
initiated when the orbital gradient falls below
this threshold. (default=0.25 au)
ETHRSH = energy
error threshold for initiating DIIS. The
DIIS error is the largest element of e=FDS-SDF.
Increasing ETHRSH forces DIIS on sooner.
(default = 0.5 Hartree)
MAXDII = Maximum
size of the DIIS linear equations, so
that at most MAXDII-1 Fock matrices are used
in the interpolation. (default=10)
DEMCUT = Direct
energy minimization will not be done
once the density matrix change falls below
this threshold. (Default=0.5)
DMPCUT = Damping
factor lower bound cutoff. The damping
damping factor will not be allowed to drop
below this value. (default=0.0)
note: The damping factor need not be zero to achieve
valid convergence (see Hsu, Davidson, and
Pitzer, J.Chem.Phys., 65, 609 (1976), see
especially the section on convergence control),
but it should not be astronomical either.
* * * * * * * * * * * * * * * * * * * * *
For more info on the convergence methods,
see the 'Further Information' section.
* * * * * * * * * * * * * * * * * * * * *
1
----- miscellaneous options -----
UHFNOS = flag
controlling generation of the natural
orbitals of a UHF function. (default=.FALSE.)
MVOQ
= 0 Skip MVO generation (default)
= n Form modified virtual orbitals, using a cation
with n electrons removed. Implemented for
RHF, ROHF, and GVB. If necessary to reach a
closed shell cation, the program might remove
n+1 electrons. Typically, n will be about 6.
= -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.
NPUNCH = SCF
punch option
= 0 do not punch out the final orbitals
= 1 punch out the occupied orbitals
= 2 punch out occupied and virtual orbitals
The default is NPUNCH = 2.
----- options for virial scaling -----
VTSCAL =
A flag to request that the virial theorem be
satisfied. An analysis of the total energy
as an exact sum of orbital kinetic energies
is printed. The default is .FALSE.
This option is implemented for RHF, UHF, and ROHF,
for RUNTYP=ENERGY, OPTIMIZE, or SADPOINT. Related
input is as follows:
SCALF =
initial exponent scale factor when VTSCAL is
in use, useful when restarting. The default
is 1.0.
MAXVT =
maximum number of iterations (at a single
geometry) to satisfy the energy virial theorem.
The default is 20.
VTCONV =
convergence criterion for the VT, which is
satisfied when 2<T> + <V> + R x dE/dR is less
than VTCONV. The default is 1.0D-6 Hartree.
For more information
on this option, which is most
economically
employed during a geometry search, see
M.Lehd and F.Jensen,
J.Comput.Chem. 12, 1089-1096(1991).
1
$SCF
The next parameters define the GVB wavefunction. Note
that ALPHA and
BETA also have meaning for ROHF. See also
MULT in the
$CONTRL group. The GVB wavefunction assumes
orbitals are
in the order core, open, pairs.
NCO
= The number of closed shell orbitals. The
default almost certainly should be changed!
(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. For example, a two pair
case for water, say, might be
CICOEF(1)=0.95,-0.05,0.95,-0.05. If not
normalized, as in the default, they will be.
This parameter is useful in restarting a GVB
run, with the current CI coefficients.
(default = 0.90,-0.20,0.90,-0.20,...)
COUPLE =
A switch controlling the input of F, ALPHA,
and BETA. The default is to use internally
stored values for these variables. Note
ALPHA and BETA can be given for -ROHF-, as
well as -GVB-. (Default=.FALSE.)
F = An vector of fractional occupations.
ALPHA =
An array of A coupling coefficients given in
lower triangular order.
BETA
= An array of B coupling coefficients given in
lower triangular order.
Note: The default for F, ALPHA, and BETA depends on
the state chosen.
Defaults for the most commonly occuring
cases are internally
stored.
* * * * * * * * * * * * * * * * * * *
For more discussion of GVB/ROHF input
see the 'further information' section
* * * * * * * * * * * * * * * * * * *
==========================================================
1
$SCFMI
==========================================================
$SCFMI group
(optional, relevant if SCFTYP=RHF)
The SCF-MI method is a modification of the Roothaan
equations that
avoids basis set superposition error (BSSE)
in intermolecular
interaction calculations, by expanding
each monomer's
orbitals using only its own basis set.
Thus, the resulting
orbitals are not orthogonal. The
presence of
a $SCFMI group in the input triggers the use
of this option.
The implementation is limited to two monomers, treated
at the RHF level.
The energy, gradient, and therefore
numerical hessian
are available. The SCF step may be run
in direct SCF
mode. The first 4 parameters must be given.
All atoms of
monomer A must be given in $DATA before the
atoms of monomer
B.
NA
= number of doubly occupied MOs on fragment A.
NB
= number of doubly occupied MOs on fragment B.
MA
= number of basis functions on fragment A.
MB
= number of basis functions on fragment B.
ITER
= maximum number of SCF-MI cycles, overriding
the usual MAXIT value. (default is 50).
DTOL
= SCF-MI density convergence criteria.
(default is 1.0d-10)
ALPHA
= possible level shift parameter.
(default is 0.0, meaning shifting is not used)
IOPT
= prints additional debug information.
= 0 standard outout (default)
= 1 print for each SCF-MI cycle MOs, overlap
between the MOs, CPU times.
= 2 print some extra informations in secular
systems solution.
MSHIFT
= debugging option that permits to shift all
the memory pointer of the SCF-MI section
of code of the quantity MSHIFT (default is 0).
==========================================================
"Modification of Roothan Equations to Exclude BSSE
from Molecular Interaction Calculations"
E. Gianinetti, M. Raimondi, E. Tornaghi
Int. J. Quantum Chem. 60, 157 (1996)
A. Famulari, E. Gianinetti, M. Raimondi, and M. Sironi
Int. J. Quantum Chem. (1997), submitted.
1
$DFT
==========================================================
$DFT group (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 and other associated
inputs are documented
separately below. Every functional
that has the
same name in both lists is the 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.
Achieving a self-consistent field with DFT is somewhat
more difficult
than for normal HF, so DIIS is the default
converger.
Since the DFT iterations are also more time
consuming, the
use of GUESS=MOREAD may be very helpful.
Both DFT programs will run in parallel.
1
$DFT
DFTTYP = NONE means no DFT is performed (default)
METHOD = selects
grid based DFT or grid free DFT.
= GRID Grid based DFT (default)
= GRIDFREE Grid free DFT
----- options for METHOD=GRID -----
DFTTYP = specifies
exchange and correlation functionals.
= SLATER Slater exchange
= BECKE Becke 1988 exchange
= GILL Gill 1996 exchange
= PBE Perdew-Burke-Ernzerhof (PBE) exchange
Note that the PBE correlation functional
is not implemented.
= SVWN SLATER + Vosko-Wilk-Nusair correlation,
using their electron gas formula 5 (VWN5)
Also known as LDA/LSDA for RHF/UHF.
= SLYP SLATER + Lee-Yang-Parr (LYP) correlation
= SOP SLATER + One-parameter Progressive
corr.
= BVWN BECKE exchange + VWN5 correlation
= BLYP BECKE exchange + LYP correlation
= BOP BECKE exchange + OP correlation
= GVWN GILL exchange + VWN5 correlation
= GLYP GILL exchange + LYP correlation
= GOP GILL exchange + OP correlation
= PBEVWN PBE exchange + VWN5 correlation
= PBELYP PBE exchange + LYP correlation
= PBEOP PBE exchange + OP correlation
= HVWN Hartree-Fock exchange + VWN5 correlation
= HLYP Hartree-Fock exchange + LYP correlation
= HOP Hartree-Fock exchange + OP/B88 correlation
= BHHLYP HF and BECKE exchange + LYP correlation
= B3LYP this is a hybrid method combining five
functionals, namely Becke + Slater + HF
exchange and LYP + VWN5 correlation.
An extensive
bibliography for these functionals can be
found in the
'Further References' section of this manual.
1
$DFT
NRAD
= number of radial grids in Euler-Maclaurin
quadrature. (default=96)
NTHE
= number of angle theta grids in Gauss-Legendre
quadrature. (default=12)
NPHI
= number of angle phi grids in Gauss-Legendre
quadrature. NPHI should be double NTHE so that
points are spherically distributed. (default=24)
NRAD*NTHE*NPHI
grid points will be constructed around each
atom.
Time is linear in the number of grid points, so be
careful.
Energies can be compared only when the identical
grid density
has been used, analogous to needing to compare
with the identical
basis set expansions. A very accurate
"army grade"
grid capable of producing an integration error
less than a
microHartree/atom is NRAD=96 NTHE=36 NPHI=72.
NRAD0, NTHE0,
NPHI0 define a smaller grid used during the
SCF iterations
before some initial convergence is reached.
After that,
the full grid defined by NRAD, NTHE, NPHI will
be used.
This can save considerable CPU time in the early
SCF iterations.
SWITCH = when
the change in the density matrix between
iterations falls below this threshhold, switch
to use of the desired full grid (default=3.0E-4)
NRAD0
= same as NRAD, but defines initial (smaller) grid.
NTHE0 =
same as NTHE, but defines initial (smaller) grid.
NPHI0
= same as NPHI, but defines initial (smaller) grid.
Default values
for the initial grid depend upon NRAD, NTHE,
and NPHI.
For the default full grid settings, the initial
grid is NRAD0=24,
NTHE0=8, NPHI0=16, for other values the
formula is NRAD0
the larger of NRAD/4 or 24, for NTHE0 the
larger of NTHE/3
or 8, and for NPHI0 the larger of NPHI/3
or 16.
In case of slow convergence of the SCF or if using
the "army grade
grid", NRAD0=48 NTHE0=12 NPHI0=24 and
SWITCH=1.0E-4
may be better. Numerical hessian runs set
the coarse grid
to the same size as the full grid, by
default.
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
1
$DFT
----- 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 + VWN correlation
= BLYP Becke exchange + LYP correlation
= BPWLOC Becke exchange + Perdew/Wang correlation
= B3LYP hybridized HF/Becke/LYP using VWN formula 5
= CAMB CAMA exchange + Cambridge correlation
= XVWN Xalpha exchange + VWN formula 5 correlation
= XPWLOC Xalpha exchange + Perdew/Wang correlation
= SVWN Slater exchange + VWN correlation
= SPWLOC Slater exchange + PWLOC correlation
= WIGNER Wigner exchange + correlation
= WS Wigner scaled exchange + correlation
= WIGEXP Wigner exponential exchange + correlation
AUXFUN = AUX0
uses no auxiliary basis set for resolution
of the identity, limiting accuracy.
= AUX3 uses the 3rd generation of RI basis sets,
These are available for the elements H to
Ar, but have been carefully considered for
H-Ne only. (DEFAULT)
THREE =
a flag to use a resolution of the identity to
turn four center overlap integrals into three
center integrals. This can be used only if
no auxiliary basis is employed. (default=.FALSE.)
==========================================================
1
$MP2
==========================================================
$MP2 group (relevant to SCFTYP=RHF,UHF,ROHF if MPLEVL=2)
Controls 2nd order Moller-Plesset perturbation runs,
if requested
by MPLEVL in $CONTRL. See also the DIRSCF
keyword in $SCF
to select direct MP2. MP2 is implemented
for RHF, high
spin ROHF, or UHF wavefunctions. Analytic
gradients and
the first order correction to the wave-
function (i.e.
properties) are available only for RHF.
The $MP2 group
is usually not given. See also $MCQDPT.
NCORE = n
Omits the first n occupied orbitals from the
calculation. The default for n is the number
of chemical core orbitals.
MP2PRP=
a flag to turn on property computation for
RHF MP2 jobs with RUNTYP=ENERGY. This is
appreciably more expensive than just evaluating
the 2nd order energy correction alone, so the
default is .FALSE. Properties are always
computed during gradient runs, when they are
an almost free byproduct. (default=.FALSE.)
LMOMP2=
a flag to analyze the closed shell MP2 energy
in terms of localized orbitals. Any type of
localized orbital may be used. This option
is implemented only for RHF, and its selection
forces use of the METHOD=3 transformation.
The default is .FALSE.
OSPT=
selects open shell spin-restricted perturbation.
This parameter applies only when SCFTYP=ROHF.
Please see the 'further information' section for
more information about this choice.
= ZAPT picks Z-averaged perturbation theory. (default)
= RMP picks RMP (aka ROHF-MBPT) perturbation theory.
CUTOFF=
transformed integral retention threshold, the
default is 1.0d-9.
1
The last 3 input
variables apply to UHF+MP2 or ROHF+MP2
using OSPT=RMP,
or to any serial MP2 run.
NWORD =
controls memory usage. The default uses all
available memory. (default=0)
METHOD= n
selects transformation method, 2 being the
segmented transformation, and 3 being a more
conventional two phase bin sort implementation.
3 requires more disk, but less memory. The
default is to attempt method 2 first, and
method 3 second.
AOINTS=
defines AO integral storage during conventional
integral transformations, during parallel runs.
DUP stores duplicated AO lists on each node, and
is the default for parallel computers with slow
interprocessor communication, e.g. ethernet.
DIST distributes the AO integral file across all
nodes, and is the default for parallel
computers with high speed communications.
==========================================================
1
$GUESS
==========================================================
$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
ECP basis set. (default for most runs)
= HCORE Diagonalize the one electron Hamiltonian
to obtain the initial guess orbitals.
This method is applicable to any basis
set, but does not work as well as the
HUCKEL guess.
= MOREAD Read in formatted vectors punched by an
earlier run. This requires a $VEC group,
and you MUST pay attention to NORB below.
= RDMINI Read in a $VEC group from a converged
calculation that used GBASIS=MINI and no
polarization functions, and project these
orbitals onto the current basis. Do not
use this option if the current basis
involve ECP basis sets.
= MOSAVED (default for restarts) The initial
orbitals are read from the DICTNRY file
of the earlier run.
= SKIP Bypass initial orbital selection.
The
initial orbitals and density matrix are
assumed to be in the DICTNRY file. Mostly
used for RUNTYP=HESSIAN when the hessian
is being read in from the input.
All GUESS types except 'SKIP' permit reordering of the
orbitals, carry
out an orthonormalization of the orbitals,
and generate
the correct initial density matrix, for RHF,
UHF, ROHF, and
GVB, but note that correct computation of
the GVB density
requires also CICOEF in $SCF. The density
matrix cannot
be generated from the orbitals alone for MP2,
CI, or MCSCF,
so property evaluation for these should be
RUNTYP=ENERGY
rather than RUNTYP=PROP using GUESS=MOREAD.
PRTMO = a flag
to control printing of the initial guess.
(default=.FALSE.)
PUNMO = a flag
to control punching of the initial guess.
(default=.FALSE.)
1
$GUESS
MIX
= rotate the alpha and beta HOMO and LUMO orbitals
so as to generate inequivalent alpha and beta
orbital spaces. This pertains to UHF singlets
only. This may require use of NOSYM=1 in $CONTRL
depending on your situation. (default=.FALSE.)
NORB
= The number of orbitals to be read in the $VEC
group. This applies only to GUESS=MOREAD.
For -RHF-, -UHF-,
-ROHF-, and -GVB-, NORB defaults to the
number of occupied
orbitals. NORB must be given for -CI-
and -MCSCF-.
For -UHF-, if NORB is not given, only the
occupied alpha
and beta orbitals should be given, back to
back.
Otherwise, both alpha and beta orbitals must
consist of NORB
vectors.
NORB may be
larger than the number of occupied MOs, if you
wish to read
in the virtual orbitals. If NORB is less
than the number
of atomic orbitals, the remaining orbitals
are generated
as the orthogonal complement to those read.
NORDER = Orbital
reordering switch.
= 0 No reordering (default)
= 1 Reorder according to IORDER and JORDER.
IORDER = Reordering
instructions.
Input to this array gives the new molecular
orbital order. For example, IORDER(3)=4,3 will
interchange orbitals 3 and 4, while leaving the
other MOs in the original order. This parameter
applies to all orbitals (alpha and beta) except
for -UHF-, where it only affects the alpha MOs.
(default is IORDER(i)=i )
JORDER = Reordering
instructions.
Same as IORDER, but for the beta MOs of -UHF-.
INSORB = the first
INSORB orbitals specified in the $VEC
group will be inserted into the Huckel guess,
making the guess a hybrid of HUCKEL/MOREAD. This
keyword is meaningful only when GUESS=HUCKEL, and
it is useful mainly for QM/MM runs where some
orbitals (buffer) are frozen and need to be
transferred to the initial guess vector set,
see $MOFRZ. (default=0)
1
$GUESS
* * * 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 totally symmetric from the density.
Maybe useful if the HUCKEL or HCORE give orbitals
of inexact symmetry. Since the density matrix is
not idempotent, this can generate a non-variational
energy on the first iteration. For the same
reason, this should never be used with orbitals
of MOREAD quality. (default=.FALSE.)
==========================================================
1
$VEC
==========================================================
$VEC group
(optional, relevant for all SCFTYP's)
(required if GUESS=MOREAD)
This group consists of formatted vectors, as written
onto file PUNCH
in a previous run. It is considered good
form to retain
the titling comment cards punched before
the $VEC card,
as a reminder to yourself of the origin of
the orbitals.
For Morokuma decompositions, the names of this group
are $VEC1, $VEC2,
... for each monomer, computed in the
identical orientation
as the supermolecule. For transition
moment or spin-orbit
coupling runs, orbitals for states
one and possibly
two are $VEC1 and $VEC2.
==========================================================
$MOFRZ group (optional, relevant for RHF, ROHF, GVB)
This group controls freezing the molecular orbitals
of your choice
during the SCF procedure. If you choose
this option,
select DIIS in $SCF since SOSCF will not
converge as
well. GUESS=MOREAD is required in $GUESS.
FRZ = flag which triggers MO freezing. (default=.FALSE.)
IFRZ =
an array of MOs in the input $VEC set which are
to be frozen. There is no default for this.
==========================================================
1
$STATPT
==========================================================
$STATPT group (optional, for RUNTYP=OPTIMIZE or SADPOINT)
This group controls the search for stationary points.
Note that NZVAR
in $CONTRL determines if the geometry
search is conducted
in Cartesian or internal coordinates.
METHOD = optimization algorithm selection. Pick from
NR Straight Newton-Raphson iterate. This will
attempt to locate the nearest stationary
point, which may be of any order. There
is no steplength control. RUNTYP can be
either OPTIMIZE or SADPOINT
RFO Rational Function Optimization. This is
one of the augmented Hessian techniques
where the shift parameter(s) is(are) chosen
by a rational function approximation to
the PES. For SADPOINT searches it involves
two shift parameters. If the calculated
stepsize is larger than DXMAX the step is
simply scaled down to size.
QA Quadratic Approximation. This is another
version of an augmented Hessian technique
where the shift parameter is chosen such
that the steplength is equal to DXMAX.
It is completely equivalent to the TRIM
method. (default)
SCHLEGEL The quasi-NR optimizer by Schlegel.
CONOPT, CONstrained OPTimization. An algorithm
which can be used for locating TSs.
The starting geometry MUST be a minimum!
The algorithm tries to push the geometry
uphill along a chosen Hessian mode (IFOLOW)
by a series of optimizations on hyperspheres
of increasingly larger radii.
Note that there currently are no restart
capabilitites for this method, not even
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. (default=20)
1
$STATPT
--- 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 refered to as "eigenvector following".
For METHOD=SCHLEGEL, the mode whose eigenvalue
is (or will be made) negative. All other
curvatures will be made positive.
For METHOD=CONOPT, the mode along which the
geometry is initially perturbed from the minima.
(default is 1)
In Cartesian coordinates, this variable doesn't
count the six translation and rotation degrees.
Note that the "modes" aren't from mass-weighting.
STPT
= flag to indicate whether the initial geometry
is considered a stationary point. If .true.
the initial geometry will be perturbed by
a step along the IFOLOW normal mode with
stepsize STSTEP. (default=.false.)
The positive direction is taken as the one where
the largest component of the Hessian mode is
positive. If there are more than one largest
component (symmetry), the first is taken as
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)
1
$STATPT
IFREEZ = array
of coordinates to freeze. These may be
internal or Cartesian coordinates. For example,
IFREEZ(1)=1,3 freezes the two bond lengths in
the $ZMAT example, while optimizing the angle.
If NZVAR=0, so that this value applies to the
Cartesian coordinates instead, the input of
IFREEZ(1)=4,7 means to freeze the x coordinates
if the 2nd and 3rd atoms in the molecule.
See also IFZMAT and FVALUE in $ZMAT, and IFCART
below, as IFREEZ does not apply to DLC internals.
In a numerical Hessian run, IFREEZ specifies
Cartesian displacements to be skipped for a
Partial Hessian Analysis. For more information:
J.D. Head, Int. J. Quantum Chem. 65, 827, 1997
H. Li, J.H. Jensen, manuscript in preparation.
IFCART = array
of Cartesian coordinates to freeze during
a geometry optimization using delocalized internal
coordinates.
--- The next two control the hessian matrix quality ---
HESS
= selects the initial hessian matrix.
= GUESS chooses a positive definite diagonal
hessian. (default for RUNTYP=OPTIMIZE)
= READ causes the hessian to be read from a $HESS
group. (default for RUNTYP=SADPOINT)
= RDAB reads only the ab initio part of the
hessian, and approximates the effective
fragment blocks.
= RDALL reads the full hessian, then converts
any fragment blocks to 6x6 T+R shape.
(this option is seldom used).
= CALC causes the hessian to be computed, see
the $FORCE group.
IHREP =
the number of steps before the hessian is
recomputed. If given as 0, the hessian will
be computed only at the initial geometry if
you choose HESS=CALC, and never again. If
nonzero, the hessian is recalculated every
IHREP steps, with the update formula used on
other steps. (default=0)
HSSEND = a flag
to control automatic hessian evaluation
at the end of a successful geometry search.
(default=.FALSE.)
1
$STATPT
--- the next two control the amount of output ---
Let 0 mean the initial geometry, L mean the last
geometry, and all mean every geometry.
Let INTR mean the internuclear distance matrix.
Let HESS mean the approximation to the hessian.
Note that a directly calculated hessian matrix
will always be punched, NPUN refers only to the
updated hessians used by the quasi-Newton step.
NPRT
= 1 Print INTR at all, orbitals at all
0 Print INTR at all, orbitals at 0+L (default)
-1 Print INTR at all, orbitals never
-2 Print INTR at 0+L, orbitals never
NPUN
= 3 Punch all orbitals and HESS at all
2 Punch all orbitals at all
1 same as 0, plus punch HESS at all
0 Punch all orbitals at 0+L, otherwise only
occupied orbitals (default)
-1 Punch occ orbitals at 0+L only
-2 Never punch orbitals
---- the following parameters are quite specialized ----
PURIFY = a flag
to help eliminate the rotational and
translational degrees of freedom from the
initial hessian (and possibly initial gradient).
This is much like the variable of the same name
in $FORCE, and will be relevant only if internal
coordinates are in use. (default=.FALSE.)
PROJCT = a flag
to eliminate translation and rotational
degrees of freedom from Cartesian optimizations.
The default is .TRUE. since this normally will
reduce the number of steps, except that this
variable is set false when POSITION=FIXED is
used during EFP runs.
ITBMAT = number
of micro-iterations used to compute the
step in Cartesians which corresponds to the
desired step in internals. The default is 5.
UPHESS = SKIP
do not update Hessian (not recommended)
BFGS default for OPTIMIZE using RFO or QA
POWELL default for OPTIMIZE using NR or CONOPT
POWELL default for SADPOINT
MSP mixed Murtagh-Sargent/Powell update
SCHLEGEL only choice for METHOD=SCHLEGEL
MOVIE =
a flag to create a series of structural data
which can be show as a movie by the MacIntosh
program Chem3D. The data is written to the
file IRCDATA. (default=.FALSE.)
1
---- 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
* * * * * * * * * * * * * * * * * * * * *
1
$TRUDGE
==========================================================
$TRUDGE group (optional, required for RUNTYP=TRUDGE)
This group defines the parameters for a non-gradient
optimization
of exponents or the geometry. The TRUDGE
package is a
modified version of the same code from Michel
Dupuis' HONDO
7.0 system, origially written by H.F.King.
Presently the
program allows for the optimization of 10
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 groups). The total
energy of all
types of SCF wavefunctions can be optimized,
although this
would be extremely stupid as gradient
methods are
far more efficient. The main utility is for
open shell MP2
or CI geometry optimizations, which may
not be done
in any other way with GAMESS. If your run
requires NOSYM=1
in $CONTRL, you must be sure to use only
C1 symmetry
in the $DATA group.
OPTMIZ = a flag
to select optimization of either geometry
or exponents of primitive gaussian functions.
= BASIS for basis set optimization.
= GEOMETRY for geometry optimization (default).
This means minima search only, there is no saddle
point capability.
NPAR = number of parameters to be optimized.
IEX = defines the parameters to be optimized.
If OPTMIZ=BASIS, IEX declares the serial number
of the Gaussian primitives for which the exponents
will be optimized.
If OPTMIZ=GEOMETRY, IEX define the pointers to
the HINT internal coordinates which will be optimized.
(Note that not all internal coordinates have to be
optimized.) The pointers to the internal coordinates
are defined as: (the number of atom on the input
list)*10 + (the number of internal coordinate for that
atom). For each atom, the HINT internal coordinates
are numbered as 1, 2, and 3 for BOND, ALPHA, and BETA,
respectively.
1
$TRUDGE
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.
==========================================================
1
$TRURST
==========================================================
$TRURST group (optional, relevant for RUNTYP=TRUDGE)
This group specifies restart parameters for TRUDGE
runs and accuracy
thresholds.
KSTART indicates
the conjugate gradient direction in which
the optimization
will proceed. ( default = -1 )
-1 .... indicates that this is a non-restart run.
0 .... corresponds to a restart run.
FNOISE accuracy
of function values.
Variation smaller
than FNOISE are not considered to be
significant
(Def. 0.0005)
TOLF accuracy required of the function (Def. 0.001)
TOLR accuracy required of conjugate directions (Def. 0.05)
For geometry optimization, the values which give
better results
(closer to the ones obtained with gradient
methods) are:
TOLF=0.0001, TOLR=0.001, FNOISE=0.00001
==========================================================
1
$FORCE
==========================================================
$FORCE group
(optional, relevant for RUNTYP=HESSIAN,OPTIMIZE,SADPOINT)
This group controls the computation of the hessian
matrix (the
energy second derivative tensor, also known
as the force
constant matrix), and an optional harmonic
vibrational
analysis. This can be a very time consuming
calculation.
However, given the force constant matrix,
the vibrational
analysis for an isotopically substituted
molecule is
very cheap. Related input is HESS= in
$STATPT, and
the $MASS, $HESS, $GRAD, $DIPDR, $VIB groups.
METHOD = chooses
the computational method.
= ANALYTIC is implemented only for SCFTYPs RHF,
ROHF, and GVB (when NPAIR is 0 or 1).
This is the default for these cases.
= NUMERIC is the default for all other cases:
UHF or MCSCF, RESC or NESC relativistic
correction, and all MP2, CI, or DFT runs.
RDHESS = a flag
to read the hessian from a $HESS group,
rather than computing it. This variable pertains
only to RUNTYP=HESSIAN. See also HESS= in the
$STATPT group. (default is .FALSE.)
PURIFY = controls
cleanup
Given a $ZMAT, the hessian and dipole derivative
tensor can be "purified" by transforming from
Cartesians to internals and back to Cartesians.
This effectively zeros the frequencies of the
translation and rotation "modes", along with
their IR intensities. The purified quantities
are punched out. Purification does change the
Hessian slightly, frequencies at a stationary
point can change by a wave number or so. The
change is bigger at non-stationary points.
(default=.FALSE. if $ZMAT is given)
PRTIFC = prints
the internal coordinate force constants.
You MUST have defined a $ZMAT group to use this.
(Default=.FALSE.)
1
$FORCE
--- the next four apply only to METHOD=NUMERIC ----
NVIB
= Number of displacements in each Cartesian
direction for force field computation.
= 1 Move one VIBSIZ unit in each positive
Cartesian direction. This requires 3N+1
evaluations of the wavefunction, energy, and
gradient, where N is the number of SYMMETRY
UNIQUE atoms given in $DATA. (default)
= 2 Move one VIBSIZ unit in the positive direction
and one VIBSIZ unit in the negative direction.
This requires 6N+1 evaluations of the
wavefunction and gradient, and gives a small
improvement in accuracy. In particular, the
frequencies will change from NVIB=1 results by
no more than 10-100 wavenumbers, and usually
much less. However, the normal modes will be
more nearly symmetry adapted, and the residual
rotational and translational "frequencies"
will be much closer to zero.
VIBSIZ = Displacement size (in Bohrs). Default=0.01
Let 0 mean the Vib0 geometry, and
D mean all the displaced geometries
NPRT
= 1 Print orbitals at 0 and D
= 0 Print orbitals at 0 only (default)
NPUN
= 2 Punch all orbitals at 0 and D
= 1 Punch all orbitals at 0 and occupied orbs at D
= 0 Punch all orbitals at 0 only (default)
----- the rest control normal coordinate analysis ----
VIBANL = flag
to activate vibrational analysis.
(the default is .TRUE. for RUNTYP=HESSIAN, and
otherwise is .FALSE.)
SCLFAC = scale
factor for vibrational frequencies, used
in calculating the zero point vibrational energy.
Some workers correct for the usual overestimate
in SCF frequencies by a factor 0.89. ZPE or other
methods might employ other factors, see A.P.Scott,
L.Radom J.Phys.Chem. 100, 16502-16513 (1996).
The output always prints unscaled frequencies, so
this value is used only during the thermochemical
analysis. (Default is 1.0)
1
$FORCE
TEMP
= an array of up to ten temperatures at which the
thermochemistry should be printed out. The
default is a single temperature, 298.15 K. To
use absolute zero, input 0.001 degrees.
FREQ
= an array of vibrational frequencies. If the
frequencies are given here, the hessian matrix
is not computed or read. You enter any imaginary
frequencies as negative numbers, omit the
zero frequencies corresponding to translation
and rotation, and enter all true vibrational
frequencies. Thermodynamic properties will be
printed, nothing else is done by the run.
PRTSCN = flag
to print contribution of each vibrational
mode to the entropy. (Default is .FALSE.)
DECOMP = activates
internal coordinate analysis.
Vibrational frequencies will be decomposed into
"intrinsic frequencies", by the method of
J.A.Boatz and M.S.Gordon, J.Phys.Chem., 93,
1819-1826(1989). If set .TRUE., the $ZMAT group
may define more than 3N-6 (3N-5) coordinates.
(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 punched is the one
BEFORE projection. For the same reason, the
default is .FALSE. to skip the projection, which
is mainly of interest in dynamical calculations.
==========================================================
There is a set
of programs for the calculation of kinetic
or equilibrium
isotope effects from the group of Piotr
Paneth at the
University of Lodz. This ISOEFF package will
accept data
computed by GAMESS, and can be downloaded at
http://ck-sg.p.lodz.pl/isoeff/isoeff.html
1
$CPHF
==========================================================
$CPHF group (relevant for analytic RUNTYP=HESSIAN)
This group controls the solution of the response
equations, also
known as coupled Hartree-Fock.
POLAR = a flag
to request computation of the static
polarizability, alpha. Because this property
needs 3 additional response vectors, beyond those
needed for the hessian, the default is to skip the
property. (default = .FALSE.)
NWORD = controls
memory usage for this step. The default
uses all available memory. (default=0)
==========================================================
1
$HESS $GRAD $DIPDR
==========================================================
$HESS group (relevant
for RUNTYP=HESSIAN if RDHESS=.TRUE.)
(relevant for RUNTYP=IRC if FREQ,CMODE not given)
(relevant for RUNTYP=OPTIMIZE,SADPOINT if HESS=READ)
Formatted force constant matrix (FCM), i.e. hessian
matrix.
This data is punched out by a RUNTYP=HESSIAN job,
in the correct
format for subsequent runs. The first card
in the group
must be a title card.
A $HESS group is always punched in Cartesians. It
will be transformed
into internal coordinate space if a
geometry search
uses internals. It will be mass weighted
(according to
$MASS) for IRC and frequency runs.
The initial FCM is updated during the course of a
geometry optimization
or saddle point search, and will be
punched if a
run exhausts its time limit. This allows
restarts where
the job leaves off. You may want to read
this FCM back
into the program for your restart, or you
may prefer to
regenerate a new initial hessian. In any
case, this updated
hessian is absolutely not suitable for
frequency prediction!
==========================================================
$GRAD group
(relevant for RUNTYP=OPTIMIZE or SADPOINT)
(relevant for RUNTYP=HESSIAN when RDHESS=.TRUE.)
Formatted gradient vector at the $DATA geometry. This
data is read
in the same format it was punched out.
For RUNTYP=HESSIAN, this information is used to
determine if
you are at a stationary point, and possibly
for projection.
If omitted, the program pretends the
gradient is
zero, and otherwise proceeds normally.
For geometry searches, this information (if known) can
be read into
the program so that the first step can be
taken instantly.
==========================================================
$DIPDR group (relevant for RUNTYP=HESSIAN if RDHESS=.T.)
Formatted dipole
derivative tensor, punched in a previous
RUNTYP=HESSIAN
job. If this group is omitted, then a
vibrational
analysis will be unable to predict the IR
intensities,
but the run can otherwise proceed.
==========================================================
1
$VIB $MASS
==========================================================
$VIB group (relevant for RUNTYP=HESSIAN, METHOD=NUMERIC)
Formatted card image -restart- data. This data is
read in the
format it was punched by a previous HESSIAN
job to the file
IRCDATA. Just add a " $END" card, and if
the final gradient
was punched as zero, delete the last
set of data.
Normally, IREST in $CONTRL will NOT be used
in conjunction
with a HESSIAN restart. The mere presence
of this deck
triggers the restart from cards. This deck
can also be
used to turn a single point differencing run
into double
differencing, as well as recovering from time
limits, or other
bombouts.
==========================================================
$MASS group (relevant for RUNTYP=HESSIAN, IRC, or DRC)
This group permits isotopic substitution during the
computation
of mass weighted Cartesian coordinates. Of
course, the
masses affect the frequencies and normal modes
of vibration.
AMASS = An array
giving the atomic masses, in amu. The
default is to use the mass of the most abundant
isotope. Masses through element 104 are stored.
example - $MASS
AMASS(3)=2.0140 $END
will make the
third atom in the molecule a deuterium.
==========================================================
1
$IRC
==========================================================
$IRC group (relevant for RUNTYP=IRC)
This group governs the location of the intrinsic
reaction coordinate,
a steepest descent path in mass
weighted corrdinates,
that connects the saddle point to
reactants and
products.
----- there are five integration methods chosen by PACE.
PACE = GS2
selects the Gonzalez-Schlegel second order
method. This is the default method.
Related input is:
GCUT
cutoff for the norm of the mass-weighted gradient
tangent (the default is chosen in the range from
0.00005 to 0.00020, depending on the value for
STRIDE chosen below.
RCUT cutoff for Cartesian RMS displacement vector.
(the default is chosen in the range 0.0005 to
0.0020 Bohr, depending on the value for STRIDE)
ACUT maximum angle from end points for linear
interpolation (default=5 degrees)
MXOPT maximum number of contrained optimization steps
for each IRC point (default=20)
IHUPD is the hessian update formula. 1 means Powell,
2 means BFGS (default=2)
GA is a gradient from the previous IRC point, and
is
used when restarting.
OPTTOL is a gradient cutoff used to determine if the IRC
is approaching a minimum. It has the same meaning
as the variable in $STATPT. (default=0.0001)
PACE = LINEAR
selects linear gradient following (Euler's
method). Related input is:
STABLZ
switches on Ishida/Morokuma/Komornicki reaction
path stabilization. The default is .TRUE.
DELTA initial step size along the unit bisector, if
STABLZ is on. Default=0.025 Bohr.
ELBOW is the collinearity threshold above which the
stabilization is skipped. If the mass weighted
gradients at QB and QC are almost collinear, the
reaction path is deemed to be curving very little,
and stabilization isn't needed. The default is
175.0 degrees. To always perform stabilization,
input 180.0.
READQB,EB,GBNORM,GB are energy and gradient data
already known at the current IRC point. If it
happens that a run with STABLZ on decides to skip
stabilization because of ELBOW, this data will be
punched to speed the restart.
1
$IRC
PACE = QUAD
selects quadratic gradient following.
Related input is:
SAB
distance to previous point on the IRC.
GA gradient vector at that historical point.
PACE = AMPC4
selects the fourth order Adams-Moulton
variable step predictor-corrector.
Related input is:
GA0,GA1,GA2
which are gradients at previous points.
PACE = RK4
selects the 4th order Runge-Kutta variable
step method. There is no related input.
----- The next two are used by all PACE choices -----
STRIDE = Determines
how far apart points on the reaction
path will be. STRIDE is used to calculate the
step taken, according to the PACE you choose.
The default is good for the GS2 method, which is
very robust. Other methods should request much
smaller step sizes, such as 0.10 or even 0.05.
(default = 0.30 sqrt(amu)-Bohr)
NPOINT = The
number of IRC points to be located in this
run. The default is to find only the next point.
(default = 1)
----- The next two let you choose your output volume -----
Let F mean the first IRC point found in this run,
and L mean the final IRC point of this run.
Let INTR mean the internuclear distance matrix.
NPRT
= 1 Print INTR at all, orbitals at all IRC points
0 Print INTR at all, orbitals at F+L (default)
-1 Print INTR at all, orbitals never
-2 Print INTR at F+L, orbitals never
NPUN
= 1 Punch all orbitals at all IRC points
0 Punch all orbitals at F+L, only occupied
orbitals at IRC points between (default)
-1 Punch all orbitals at F+L only
-2 Never punch orbitals
1
$IRC
----- The next
two tally the reaction path results. The
defaults are appropriate for starting from a saddle
point, restart values are automatically punched out.
NEXTPT = The
number of the next point to be computed.
STOTAL = Total
distance along the reaction path to next
IRC point, in mass weighted Cartesian space.
----- The following
controls jumping off the saddle point.
If you give a $HESS group, FREQ and CMODE will be
generated automatically.
SADDLE = A logical
variable telling if the coordinates
given in the $DATA deck are at a saddle point
(.TRUE.) or some other point lying on the IRC
(.FALSE.). If SADDLE is true, either a $HESS
group or else FREQ and CMODE must be given.
(default = .FALSE.) Related input is:
TSENGY = A logical
variable controlling whether the energy
and wavefunction are evaluated at the transition
state coordinates given in $DATA. Since you
already know the energy from the transition state
search and force field runs, the default is .F.
FORWRD = A logical
variable controlling the direction to
proceed away from a saddle point. The forward
direction is defined as the direction in which
the largest magnitude component of the imaginary
normal mode is positive. (default =.TRUE.)
EVIB
= Desired decrease in energy when following the
imaginary normal mode away from a saddle point.
(default=0.0005 Hartree)
FREQ
= The magnitude of the imaginary frequency, given
in cm**-1.
CMODE
= An array of the components of the normal mode
whose frequency is imaginary, in Cartesian
coordinates. Be careful with the signs!
You must give FREQ and CMODE if you don't give a $HESS
group, when SADDLE=.TRUE. The option of giving these
two variables instead of a $HESS does not apply to the
GS2 method, which must have a hessian input, even for
restarts. Note also that EVIB is ignored by GS2 runs.
==========================================================
* * * * * * * * * * * * * * * * * *
For hints about IRC tracking, see
the 'further information' section.
* * * * * * * * * * * * * * * * * *
1
$VSCF
==========================================================
$VSCF group (optional, relevant to RUNTYP=VSCF)
This group governs the computation of frequencies
including anharmonic
effects. Besides the values shown
below, the input
file must contain a $HESS group and
perhaps a $DIPDR
group, to start with previously obtained
harmonic vibrational
information. Energies are sampled
along the directions
of harmonic normal modes, and along
pairs of harmonic
normal modes, after which vibrational
nuclear wavefunctions
are obtained at an SCF-like level,
termed VSCF,
using product nuclear wavefunctions. An
MP2-like correction
to the vibrational energy, termed
correlation
corrected (cc-VSCF), is also obtained. In
addition, degenerate
pertubation theory is performed,
based on a CI
reference that includes linear combinations
of degenerate
states. By default, the dipole is computed
at every grid
point to give improved IR intensity values.
See also the
restart group $VIBSCF.
NGRID =
number of grid points to be computed along each
harmonic normal mode, and if NCOUP=2, along each
pair of modes. Reasonable values are 8 or 16,
with 16 considered significantly more accurate.
(default=16)
NCOUP =
the order of mode couplings included.
= 1 computes 1-D grids along each harmonic mode
= 2 adds additionally, 2-D grids along each pair
of normal modes. (default)
The total number of energy and dipole evaluations
for NCOUP=2 is M*NGRID + M*(M-1)/2*NGRID**2, where
M is the number of normal modes: M = 3N-6 or 3N-5.
The next five
relate to the solver for the vibrational
states.
The default is a degenerate perturbation theory
treatment including
the ground and every singly excited
vibrational
level.
VDPT
= flag to use 2nd order degenerate perturbation
theory to find vibrational energys. Turning this
off causes only a CI singles and doubles treatment
to be made. (default=.true.)
ICASX
= vibrational excitation level to include in the
solver's basis. 1,2... mean first, second...
excitations will be included. The default, 1,
includes single quantum excited states only.
ICAS1, ICAS2
= starting and ending vibrations whose quanta
are included, according to ICASX. The default is
all modes, ICAS1=1 and ICAS2=3N-6.
SFACT
= a numerical cutoff for small contributions in
the solver. The default is 1d-5.
1
$VSCF
CASMIN = a flag,
largely redundant, that ensures default
settings for
The solver finds
the ground state (v=0) by default, but
will readily
find excited levels (such as all v=1) if
restarted.
Note that IEXC is one greater than the sum of
the vibrational
quantum numbers.
IEXC
= 1 obtain fundamental frequencies (default)
= 2 instead, obtain first overtones
= 3 instead, obtain second overtones
IEXC higher than 1 may be speedily obtained
using the next parameter to restart with a
completed $VIBSCF group.
READV =
flag to indicate restart data $VIBSCF should be
read in to resume an interrupted calculation, or
to obtain overtones in follow-on runs.
(default is .FALSE.)
The next two
relate to simplified intensity computation.
These simplifications
are aimed at speeding up MP2 runs,
if one cares
not so much about intensities, and so would
like to reduce
CPU for computing dipoles. It is pointless
to select DMDR
for SCF electronic structure, where the
dipoles are
easily obtainable. DMDR must not be used if
overtones are
being computed.
DMDR
= if true, indicates that the harmonic dipole
derivative tensor $DIPDR is input, rather than
computing the dipoles. (default is .FALSE.)
MPDIP =
for MP2 electronic structure, a value of .FALSE.
uses SCF level dipoles in order to save the time
needed to obtain the MP2 density at every grid
point. It is more accurate to use the DMDR flag
instead of this option, if $DIPDR is available.
Obviously this variable is irrelevant for SCF
level electronic structure. (default=.TRUE.)
VCFCT =
scaling factor for pair-coupling potential.
Sometimes when pair-coupling potential values
are larger than the corresponding single mode
values, they must be scaled down. (Default=1.0)
==========================================================
$VIBSCF group (optional, relevant to RUNTYP=VSCF)
This is restart
data, as written to file IRCDATA in a
partially completed
previous run. Append a " $END" line,
and select READV=.TRUE.
to read the data.
==========================================================
1
$DRC
==========================================================
$DRC group (relevant for RUNTYP=DRC)
This group governs the dynamical reaction coordinate,
a classical
trajectory method based on quantum chemical
potential energy
surfaces. In GAMESS these may be either
ab initio or
semi-empirical. Because the vibrational
period of a
normal mode with frequency 500 wavenumbers is
67 fs, a DRC
needs to run for many steps in order to
sample a representative
portion of phase space. Almost
all DRCs break
molecular symmetry, so build your molecule
with C1 symmetry
in $DATA, or specify NOSYM=1 in $CONTRL.
Restart data
can be found in the job's OUTPUT file, with
important results
summarized to the IRCDATA file.
NSTEP =
The number of DRC points to be calculated, not
including the initial point. (default = 1000)
DELTAT = is the time step. (default = 0.1 fs)
TOTIME = total
duration of the DRC computed in a previous
job, in fs. The default is the correct value
when initiating a DRC. (default=0.0 fs)
* * *
In general, a DRC can be initiated anywhere,
so $DATA might contain coordinates of the
equilibrium geometry, or a nearby transition
state, or something else. You must also
supply an initial kinetic energy, and the
direction of the initial velocity, for which
there are a number of options:
EKIN = The initial kinetic energy (default
= 0.0 kcal/mol)
See also ENM, NVEL, and VIBLVL regarding alternate
ways to specify the initial value.
VEL
= an array of velocity components, in Bohr/fs.
When NVEL is false, this is simply the direction
of the velocity vector. Its magnitude will be
automatically adjusted to match the desired initial
kinetic energy, and it will be projected so that
the translation of the center of mass is removed.
Give in the order vx1, vy1, vz1, vx2, vy2, ...
NVEL = a flag to compute the initial kinetic
energy from
the input VEL using the sum of mass*VEL*VEL/2.
This flag is usually selected only for restarts.
(default=.FALSE.)
1
$DRC
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 a
$HESS group, 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 a $HESS group at 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 component of the normal
mode is a positive value. (default=0)
ENM = the initial kinetic energy given
to mode NNM,
in units of vibrational quanta hv, so the amount
depends on mode NNM's vibrational frequency, v.
If you prefer to impart an arbitrary initial
kinetic energy to mode NNM, specify EKIN instead.
(default = 0.0 quanta)
1
$DRC
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 a $HESS2 group 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.
1
$DRC
The next variables can cause termination of a run, if
molecular fragments
get too far apart or close together.
NFRGPR = Number
of atom pairs whose distance will be
checked. (default is 0)
IFRGPR = Array of the atom pairs. 2 times NFRGPR values.
FRGCUT = Array
for a boundary distance (in Bohr) for atom
pairs to end DRC calculations. The run will
stop if any distance exceeds the tolerance, or if
a value is given as a negative number, if the
distance becomes shorter than the absolute value.
In case the trajectory starts outside the bounds
specified, they do not apply until after the
trajectory reaches a point where the criteria
are satisfied, and then goes outside again.
Give NFRGPR values.
* * *
The final variables control the volume of output.
Let F mean the first DRC point found in this run,
and L mean the last DRC point of this run.
NPRTSM = summarize
the DRC results every NPRTSM steps,
to the file IRCDATA. (default = 1)
NPRT
= 1 Print orbitals at all DRC points
0 Print orbitals at F+L (default)
-1 Never print orbitals
NPUN
= 2 Punch all orbitals at all DRC points
1 Punch all orbitals at F+L, and occupied
orbitals at DRC points between
0 Punch all orbitals at F+L only (default)
-1 Never punch orbitals
==========================================================
References: see REFS.DOC.
1
$GLOBOP
==========================================================
$GLOBOP group (optional, relevant to RUNTYP=GLOBOP)
This controls a search for the global minimum energy
when effective
fragment are being used. It has options
for a single
temperature Monte Carlo search, or a multi-
temperature
simulated annealing. Local minimization of
some or all
of the structures selected by the Monte Carlo
is optional.
See REFS.DOC for an overview of this RUNTYP.
TEMPI =
initial temperature used in the simulation.
(default = 40000 K)
TEMPF =
final temperature. If TEMPF is not given and
NTEMPS is greater than 1, TEMPF will be
calculated based on a cooling factor of 0.95.
NTEMPS =
number of temperatures used in the simulation.
If NTEMPS is not given but TEMPF is given,
NTEMP will be calculated based on a cooling
factor of 0.95. If neither NTEMP nor TEMPF is
given, the job defaults to a single temperature
Monte Carlo calculation.
NFRMOV =
number of fragments to move on each step.
(default=1)
NGEOPT =
number of geometries to be evaluated at each
temperature. (default = 10)
NTRAN =
number of translational steps in each block.
(default=1)
NROT
= number of rotational steps in each block.
(default=1)
NBLOCK =
the number of blocks of steps can be set directly
with this variable, instead of being calculated
from NGEOPT, NTRAN, and NROT. If NBLOCK is input,
the number of geometries at each temperature will
be NBLOCK*(NTRAN+NROT). Each block has NTRAN
translational steps followed by NROT rotational
steps.
MCMIN =
flag to enable geometry optimization to minimize
the energy is carried out every NSTMIN steps.
(default=.true.)
NSTMIN =
After this number of geometry steps are taken, a
local (Newton-Raphson) optimization will be
carried out. If this variable is set to 1, a
local minimization is carried out on every step,
reducing the MC space to the set of local minima.
Irrelevant if MCMIN is false. (default=10)
1
$GLOBOP
OPTN
= if set to .TRUE., at the end of the run local
minimizations are carried out on the final
geometry and on the minimum-energy geometry.
(default=.FALSE.)
SCALE =
an array of length two. The first element is the
initial maximum step size for the translational
coordinates (Angstroms). The second element is
the initial maximum stepsize for the rotational
coordinates (pi-radians). (defaults = 1,1)
SMODIF =
scale factor for moving ab initio atoms in the
MC simulation. If set to zero, the ab initio
atoms do not move. (default=0.1)
ALPHA =
controls the rate at which information from
successful steps is folded into the maximum step
sizes for each of the 6*(number of fragments)
coordinates. ALPHA varies between 0 and 1.
ALPHA=0 means do not change the maximum step
sizes, and ALPHA=1 throws out the old step sizes
whenever there is a successful step and uses the
successful step sizes as the new maxima. This
update scheme was used with the Parks method
where all fragments are moved on every step. It
is normally not used with the Metropolis method.
(default = 0)
DACRAT =
the desired acceptance ratio, the program tries
to achieve this by adjusting the maximum step
size. (default = 0.5)
UPDFAC =
the factor used to update the maximum step size
in the attempt to achive the desired acceptance
ratio (DACRAT). If the acceptance ratio at the
previous temperature was below DACRAT, the step
size is decreased by multiplying it by UPDFAC.
If the acceptance ratio was above DACRAT, the
step size is increased by dividing it by DACRAT
It should be between 0 and 1. (default = 0.95)
SEPTOL =
the separation tolerence between atoms in the ab
initio piece and atoms in the fragments, as well
as between atoms in different fragments. If a
step moves atoms closer than this tolerence, the
step is rejected. (default = 1.5 Angstroms)
XMIN, XMAX, YMIN,
YMAX, ZMIN, ZMAX = mimimum and maximum
values for the Cartesian coordinates of the
fragment. If the first point in a fragment steps
outside these boundaries, periodic boundary
conditions are used and the fragment re-enters on
the opposite side of the box. The defaults of
-10 for minima and +10 for maxima should usually
be changed.
1
$GLOBOP
BOLTWT =
method for calculating the Boltzmann factor,
which is used as the probability of accepting a
step that increases the energy.
= STANDARD = use the standard Boltzmann factor,
exp(-delta(E)/kT) (default)
= AVESTEP = scale the temperature by the average
step size, as recommended in the Parks reference
when using values of ALPHA greater than 0.
NPRT
= controls the amount of output, with
= -2 reduces output below that of -1
= -1 reduces output further, needed for MCMIN=.true.
= 0 gives minimal output (default)
= 1 gives the normal GAMESS amount of output
= 2 gives maximum output
For large simulations, even IOUT=0 may produce
a log file too large to work with easily.
RANDOM =
controls the choice of random number generator.
= DEBUG uses a simple random number generator with
a constant seed. Since the same sequence of
random numbers is generated during each job, it
is useful for debugging.
= RAND1 uses the simple random number generator
used in DEBUG, but with a variable seed.
= RAND3 uses a more sophisticated random number
generator described in Numerical Recipes, with a
variable seed (default).
IFXFRG =
array whose length is the number of fragments.
It allows one or more fragments to be fixed
during the simulation.
=0 allows the fragment to move during the run
=1 fixes the fragment
For example, IFXFRG(3)=1 would fix the third
fragment, the default is IFXFRG(1)=0,0,0,...,0
MOVIE2 = a flag
to create a series of structural data
which can be shown as a movie by the MacIntosh
program Chem3D. The coordinates of each accepted
geometry are written. The data is written to the
file IRCDATA. (default=.FALSE.)
==========================================================
1
$GRADEX
==========================================================
$GRADEX group (optional, for RUNTYP=GRADEXTR)
This group controls the gradient extremal following
algorithm.
The GEs leave stationary points parallel to
each of the
normal modes of the hessian. Sometimes a GE
leaving a minimum
will find a transition state, and thus
provides us
with a way of finding that saddle point. GEs
have many unusual
mathematical properties, and you should
be aware that
they normally differ a great deal from IRCs.
The search will always be performed in cartesian
coordinates,
but internal coordinates along the way may
be printed by
the usual specification of NZVAR and $ZMAT.
METHOD = algorithm
selection.
SR A predictor-corrector method due to Sun
and Ruedenberg (default).
JJH A method due to Jorgensen, Jensen and
Helgaker.
NSTEP =
maximum number of predictor steps to take.
(default=50)
DPRED =
the stepsize for the predictor step.
(default = 0.10)
STPT
= a flag to indicate whether the initial geometry
is considered a stationary point. If .TRUE.,
the geometry will be perturbed by STSTEP along
the IFOLOW normal mode.
(default = .TRUE.)
STSTEP = the
stepsize for jumping away from a stationary
point. (default = 0.01)
IFOLOW = Mode
selection option. (default is 1)
If STPT=.TRUE., the intial geometry will be
perturbed by STSTEP along the IFOLOW normal mode.
Note that IFOLOW can be positive or negative,
depending on the direction the normal mode
should be followed in. The positive direction
is defined as the one where the largest component
of the Hessian eigenvector is positive.
If STPT=.FALSE. the sign of IFOLOW determines
which direction the GE is followed in. A positive
value will follow the GE in the uphill direction.
The value of IFOLOW should be set to the Hessian
mode which is parallel to the gradient to avoid
miscellaneous warning messages.
1
$GRADEX
GOFRST = a flag
to indicate whether the algorithm should
attempt to locate a stationary point. If .TRUE.,
a straight NR search is performed once the NR
step length drops below SNRMAX. 10 NR step are
othen allowed, a value which cannot be changed.
(default = .TRUE.)
SNRMAX = upper
limit for switching to straight NR search
for stationary point location.
(default = 0.10 or DPRED, whichever is smallest)
OPTTOL = gradient
convergence tolerance, in Hartree/Bohr.
Used for optimizing to a stationary point.
Convergence of a geometry search requires the
rms gradient to be less than OPTTOL.
(default=0.0001)
HESS
= selection of the initial hessian matrix, if
STPT=.TRUE.
= READ causes the hessian to be read from a $HESS
group.
= CALC causes the hessian to be computed. (default)
1
$GRADEX
---- parameters on this page apply only to METHOD=SR ----
DELCOR = the
corrector step should be smaller than this
value before the next predictor step is taken.
(default = 0.001)
MYSTEP = maximum
number of micro iteration allowed to
bring the corrector step length below DELCOR.
(default=20)
SNUMH =
stepsize used in the numerical differentiation
of the Hessian to produce third derivatives.
(default = 0.0001)
HSDFDB = flag
to select determination of third derivatives.
At the current geometry we need the gradient, the
Hessian, and the partial third derivative matrix
in the gradient direction.
If .TRUE., the gradient is calculated at the
current geometry, and two Hessians are calculated
at SNUMH distance to each side in the gradient
direction. The Hessian at the geometry is formed
as the average of the two displaced Hessians.
If .FALSE., both the gradient and Hessian are
calculated at the current geometry, and one
additional Hessian is calculated at SNUMH in the
gradient direction.
The default double-sided differentiation produces
a more accurate third derivative matrix, at the
cost of an additional wave function and gradient.
(default = .TRUE.)
==========================================================
* * * * * * * * * * * * * * * * * * *
See the 'further information' section
for some help with GRADEXTR runs.
* * * * * * * * * * * * * * * * * * *
1
$SURF
==========================================================
$SURF group (relevant for RUNTYP=SURFACE)
This group allows you to probe a potential energy
surface along
a small grid of points. Note that there is
no option to
vary angles, only distances. The scan can
be made for
any SCFTYP, or for the MP2 or CI surface. You
may specify
two rather different calculations to be done
at each point
on the grid, through the RUNTYPn, SCFTYPn,
and CITYPn keywords.
* * * below, 1 and 2 refer to different calculations * * *
RUNTP1,RUNTYP2
= some RUNTYP supported in $CONTRL
First RUNTYP=RUNTP1 and then RUNTYP=RUNTP2 will be
performed, for each point on the grid. The second
run is omitted if RUNTP2 is set to NONE.
default: RUNTP1=ENERGY RUNTP2=NONE
SCFTP1,SCFTYP2
= some SCFTYP supported in $CONTRL
default: SCFTYP in $CONTRL
CITYP1,CITYP2
= some CITYP supported in $CONTRL
default: CITYP in $CONTRL
MPLEV1,MPLEV2
= some MPLEVL supported in $CONTRL
default: MPLEVL in $CONTRL
You may need
to help by giving values in $CONTRL that will
permit the program
to estimate what is coming in the values
here.
For example, if you want to request hessians here,
it may be good
to give RUNTYP=HESSIAN in $CONTRL so that
in its earliest
stages of a job, the program can initialize
for 2nd derivatives.
There is less checking here than on
$CONTRL input,
so don't request something impossible such
as CI and MP2
simultaneously, or analytic hessians for MP2,
or other things
that are impossible.
* * * below, 1 and 2 refer to different coordinates * * *
IVEC1 =
an array of two atoms, defining a coordinate from
the first atom given, to the second.
IGRP1 =
an array specifying a group of atoms, which must
include the second atom given in IVEC1. The
entire group will be translated (rigidly) along
the vector IVEC1, relative to the first atom
given in IVEC1.
ORIG1 =
starting value of the coordinate, which may be
positive or negative. Zero corresponds to the
distance given in $DATA.
1
$SURF
DISP1 =
step size for the coordinate. If DISP1 is set
to zero, then the keyword GRID1 is read.
NDISP1 = number of steps to take for this coordinate.
GRID1 =
an array of grid points at which to compute the
energy. This option is an alternative to the
ORIG1, DISP1 input which produces an equidistant
grid. To use GRID1, one has to set DISP1=0.0.
The number of grid points is given in NDISP1, and
is limite to at most 100 grid points. The input
of GRID1(1)=ORIG1,ORIG1+DISP1,ORIG1+DISP1*2,...
would reproduce an equidistant grid given by ORIG1
and DISP1.
ORIG1, DISP1, and GRID1 should be given in Angstrom.
There are no reasonable defaults for these keywords.
IVEC2, IGRP2,
ORIG2, DISP2, NDISP2, GRID2 have the same
meaning as their
"1" counterparts, and permit you to make
a two dimensional
map along two displacement coordinates.
If the "2" data
are not input, the surface map proceeds in
only one dimension.
==========================================================
1
$LOCAL
==========================================================
$LOCAL group (relevant for LOCAL=RUEDNBRG, BOYS, or POP)
This group allows input of additional data to control
the localization
methods. If no input is provided, the
valence orbitals
will be localized as much as possible,
while still
leaving the wavefunction invariant. There are
many specialized
options for Localized Charge Distribution
analysis, and
for EFP generation.
N.B. Since
Boys localization needs the dipole integrals,
do not turn off dipole moment calculation in $ELMOM.
MAXLOC = maximum
number of localization cycles. This
applies to BOYS or POP methods only. If the
localization fails to converge, a different
order of 2x2 pairwise rotations will be tried.
(default=250)
CVGLOC = convergence
criterion. The default provides
LMO coefficients accurate to 6 figures.
(default=1.0E-6)
SYMLOC = a flag
to restrict localization so that orbitals
of different symmetry types are not mixed. This
option is not supported in all possible point
groups. The purpose of this option is to give a
better choice for the starting orbitals for GVB-PP
or MCSCF runs, without destroying the orbital's
symmetry. This option is compatible with each of
the 3 methods of selecting the orbitals to be
included. (default=.FALSE.)
ORIENT = a flag
to request orientation of the localized
orbitals for bond-order analysis. After the
localization, the orbitals on each atom are
rotated only among themselves, in order to direct
the orbitals towards neighboring atom's orbitals,
to which they are bonded. The density matrix,
or bond-order matrix, of these Oriented LMOs is
readily interpreted as atomic populations and
bond orders. This option can be used only for
SCFTYP=MCSCF and LOCAL=RUEDENBRG.
(default=.FALSE.)
PRTLOC = a flag
to control supplemental printout. The
extra output is the rotation matrix to the
localized orbitals, and, for the Boys method,
the orbital centroids, for the Ruedenberg
method, the coulomb and exchange matrices,
for the population method, atomic populations.
(default=.FALSE.)
1
$LOCAL
----- The
following keywords select the orbitals which
are to be included in the localization. You may
select from FCORE, NOUTA/NOUTB, or NINA/NINB,
but may choose only one of these three groups.
FCORE =
flag to freeze all the chemical core orbitals
present. All the valence orbitals will be
localized. You must explicitly turn this
option off to choose one of the other two
orbital selection options. (default=.TRUE.)
* * *
NOUTA =
number of alpha orbitals to hold fixed in the
localization. (default=0)
MOOUTA = an array
of NOUTA elements giving the numbers of
the orbitals to hold fixed. For example, the
input NOUTA=2 MOOUTA(1)=8,13 will freeze only
orbitals 8 and 13. You must enter all the
orbitals you want to freeze, including any cores.
This variable has nothing to do with cows.
NOUTB =
number of beta orbitals to hold fixed in -UHF-
localizations. (default=0)
MOOUTB = same
as MOOUTA, except that it applies to the
beta orbitals, in -UHF- wavefunctions only.
* * *
NINA
= number of alpha orbitals which are to be
included in the localization. (default=0)
MOINA =
an array of NINA elements giving the numbers of
the orbitals to be included in the localization.
Any orbitals not mentioned will be frozen.
NINB
= number of -UHF- beta MOs in the localization.
(default=0)
MOINB =
same as MOINA, except that it applies to the
beta orbitals, in -UHF- wavefunctions only.
----- The following
keywords are used for the localized
charge distribution (LCD) energy decomposition.
EDCOMP = flag
to turn on LCD energy decomposition.
Note that this method is currently implemented
for SCFTYP=RHF and ROHF and LOCAL=RUEDNBRG only.
The SCF LCD forces all orbitals to be localized,
overriding input on the previous page. See also
LMOMP2 in the $MP2 group. (default = .FALSE.)
1
$LOCAL
MOIDON = flag
to turn on LMO identification and subsequent
LMO reordering, and assign nuclear LCD automat-
ically. (default = .FALSE.)
DIPDCM = flag
for LCD molecular dipole decomposition.
(default = .FALSE.)
QADDCM = flag
for LCD molecular quadrupole decomposition.
(default = .FALSE.)
POLDCM = flag
to turn on LCD polarizability decomposition.
This method is implemented for SCFTYP=RHF or ROHF
and LOCAL=BOYS or RUEDNBRG. (default=.FALSE.)
POLNUM = flag
to forces numerical rather than analytical
calculation of the polarizabilities. This may be
useful in larger molecules. The numerical
polarizabilities of bonds in or around aromatic
rings sometimes are unphysical. (default=.FALSE.)
See D.R.Garmer, W.J.Stevens
J.Phys.Chem. 93, 8263-8270 (1989).
POLAPP = flag
to force calculation of the polarizabilities
using a perturbation theory expression. This may
be useful in larger molecules. (default=.FALSE.)
See R.M. Minikis, V. Kairys, J.H. Jensen
J.Phys.Chem.A 2000, submitted.
POLANG = flag
to choose units of localized polarizability
output. The default is Angstroms**3, while false
will give Bohr**3. (default=.TRUE.)
ZDO
= flag for LCD analysis of a composite wave function,
given in a $VEC group of a van der Waals complex,
within the zero differential overlap approximation.
The MOs are not orthonormalized and the inter-
molecular electron exchange energy is neglected.
In addition, the molecular overlap matrix is printed
out. This is a very specialized option.
(default = .FALSE.)
----- The following
keywords can be used to define the
nuclear part of an LCD. They are usually used to
rectify mistakes in the automatic definition
made when MOIDON=.TRUE. The index defining the
LMO number then refers to the reordered list of LMOs.
NMOIJ =
array giving the number of nuclei assigned to a
particular LMO.
IJMO
= is an array of pairs of indices (I,J), giving
the row (nucleus I) and column (orbital J)
index of the entries in ZIJ and MOIJ.
1
$LOCAL
MOIJ
= arrays of integers K, assigning nucleus K as the
site of the Ith charge of LCD J.
ZIJ
= array of floating point numbers assigning a
charge to the Ith charge of LCD J.
IPROT =
array of integers K, defining nucleus K as a
proton.
DEPRNT = a flag
for additional decomposition printing,
such as pair contributions to various energy
terms, and centroids of the Ruedenberg orbitals.
(default = .FALSE.)
----- The following
keywords are used to build large EFPs
from several RUNTYP=MAKEFP runs on smaller molecular
fragments, by excluding common regions of overlap.
For example, an EFP for n-octanol can be build from
two MAKEFP runs, on n-pentane and n-pentanol,
CH3CH2CH2CH2-CH2CH2CH2CH2OH
CH3CH2CH2CH2[-CH3]
[CH3]-CH2CH2CH2CH2OH
by excluding operlapping regions shown in brackets
from the two EFPs. See J.Phys.Chem.A 105, 3829-3837,
(2001) for more information.
NOPATM = array
of atoms that define an area to be excluded
from a DMA ($STONE) during a RUNTYP=MAKEFP run.
All atomic centers specified, and the midpoints
of any bonds to them, are excluded as expansion
points. The density due to all LMOs primarily
centered on these atoms are excluded from the DMA
(see also KMIDPT). Furthermore, polarizability
tensors for these LMOs are excluded.
KPOINT = array
of "boundary atoms", those atoms that are
covalently bonded to the atoms given in NOATM.
KMIDPT = flag
to indicate whether the density due to bond
LMOs (and associated expansion points) between
the NOPATM atoms and the KPOINT atoms are to be
included in the DMA. (default = .TRUE.)
==========================================================
* * * * * * * * * * * * * * * * * *
For hints about localizations, and
the LCD energy decomposition, see
the 'further information' section.
* * * * * * * * * * * * * * * * * *
1
$TWOEI
==========================================================
$TWOEI group (relevant for EDCOMP=.TRUE. in $LOCAL)
Formatted transformed two-electron Coulomb and Exchange
integrals as
punched during a LOCAL=RUEDNBRG run. If this
group is present
it will automaticall be read in during
such a run and
the two-electron integrals do not have to
be re-transformed.
This group is especially useful for
EDCOMP=.TRUE.
runs when the localization has to be repeated
for different
definitions of nuclear LCDs.
1
$TRUNCN
==========================================================
$TRUNCN group (optional, relevant for RHF)
This group controls the truncation of some of the
localized orbitals
to just the AOs on a subset of the
atoms.
This option is particularly useful to generate
localized orbitals
to be frozen when the effective
fragment potential
is used to partition a system across a
chemical bond.
In other words, this group prepares the
frozen buffer
zone orbitals. This group should be used in
conjunction
with RUNTYP=ENERGY (or PROP if the orbitals
are available)
and either LOCAL=RUEDNBRG or BOYS, with
MOIDON set in
$LOCAL.
DOPROJ = flag
to activate MO projection/truncation, the
default is to skip this (default=.FALSE.)
AUTOID = forces
identification of MOs (analogous to MOIDON
in $LOCAL). This keyword is provided in case the
localized orbitals are already present in $VEC,
in which case this is a faster RUNTYP=PROP with
LOCAL=NONE job. Obviously, GUESS=MOREAD.
(default=.FALSE.)
PLAIN =
flag to control the MO tail truncation. A value
of .FALSE. uses corresponding orbital projections,
H.F.King, R.E.Stanton, H.Kim, R.E.Wyatt, R.G.Parr
J. Chem. Phys. 47, 1936-1941(1967) and generates
orthogonal orbitals. A value of .TRUE. just sets
the unwanted AOs to zero, so the resulting MOs
need to go through the automatic orthogonalization
step when MOREAD in the next job. (default=.FALSE.)
IMOPR =
an array specifying which MOs to be truncated. In
most cases involving normal bonding, the options
MOIDON or AUTOID will correctly identify all
localized MOs belonging to the atoms in the zone
being truncated. However, you can inspect the
output, and give a list of all MOs which you want
to be truncated in this array, in case you feel
the automatic assignment is incorrect.
Any orbital not in the truncation set, whether
this is chosen automatically or by IMOPR, is left
completely unaltered.
1
- - -
There are now
two ways to specify what orbitals are to
be truncated.
The most common usage is for preparation of
a buffer zone
for QM/MM computations, with an Effective
Fragment Potential
representing the non-quantum part of
the system.
This input is NATAB, NATBF, ICAPFR, ICAPBF,
in which case
the $DATA input must be sorted into three
zones.
The first group of atoms are meant to be treated
in later runs
by full quantum mechanics, the second
group by frozen
localized orbitals as a 'buffer', and the
third group
is to be substituted later by an effective
fragment potential
(multipoles, polarizabilities, ...).
Note that in
the DOPROJ=.TRUE. run, all atoms are still
quantum atoms.
NATAB = number of atoms to be in the 'ab initio' zone.
NATBF =
number of atoms to be in the 'buffer' zone.
The program can obtain the number of atoms in
the remaining zone by subtraction, so it need
not be input.
In case the MOIDON
or AUTOID options lead to confused
assignments
(unlikely in ordinary bonding situations
around the buffer
zone), there are two fine tuning values.
ICAPFR = array
indicating the identity of "capping atoms"
which are on the border between the ab initio and
buffer zones (in the ab initio zone).
ICAPBK = array
indicating the identity of "capping atoms"
which are on the border between the buffer and EFP
zones (in the effective fragment zone).
See also IXCORL and IXLONE below.
- - -
In case truncation
seems useful for some other purpose,
you can specify
the atoms in any order within the $DATA
group, by the
IZAT/ILAT approach. You are supposed to
give only one
of these two lists, probably whichever is
shorter:
IZAT
= an array containing the atoms which are NOT in
the buffer zone.
ILAT
= an array containing the atoms which are in
the buffer zone.
The AO coefficients
of the localized orbitals present in
the buffer zone
which lie on atoms outside the buffer will
be truncated.
See also IXCORL and IXLONE below.
1
- - -
The next two
values let you remove additional orbitals
within the buffer
zone from the truncation process, if that
is desirable.
These arrays can only include atoms that are
already in the
buffer zone, whether this was defined by
NATBF, or IZAT/ILAT.
The default is to include all core
and lone pair
orbitals, not just bonding orbitals, as the
buffer zone
orbitals.
IXCORL = an array
of atoms whose core and lone pair
orbitals are to be considered as not belonging
to the buffer zone orbitals.
IXLONE = an array
of atoms for which only the lone pair
orbitals are to be considered as not belonging
to the buffer zone orbitals.
- - -
The final option
controls output of the truncated orbitals
to file PUNCH
for use in later runs:
NPUNOP =
punch out option for the truncated orbitals
= 1 the MOs are not reordered.
= 2 punch the truncated MOs as the first vectors
in the $VEC MO set, with untransformed vectors
following immediately after. (default)
==========================================================
1
$ELMOM $ELPOT
==========================================================
$ELMOM group
(not required)
This group controls
electrostatic moments calculation.
The symmetry
properties of multipoles are discussed in
A.Gelessus, W.Thiel, W.Weber
J.Chem.Ed. 72, 505-508(1995)
IEMOM =
0 - skip this property
1 - calculate monopole and dipole (default)
2 - also calculate quadrupole moments
3 - also calculate octupole moments
WHERE =
COMASS - center of mass (default)
NUCLEI - at each nucleus
POINTS - at points given in $POINTS.
OUTPUT = PUNCH, PAPER, or BOTH (default)
IEMINT = 0 -
skip printing of integrals (default)
1 - print dipole integrals
2 - also print quadrupole integrals
3 - also print octupole integrals
-2 - print quadrupole integrals only
-3 - print octupole integrals only
The quadrupole and octupole tensors on the printout
are formed according
to the definition of Buckingham.
Caution: only
the first nonvanishing term in the multi-
ipole charge
expansion is independent of the coordinate
origin chosen,
which is normally the center of mass.
==========================================================
$ELPOT group (not required)
This group controls electrostatic potential calculation.
IEPOT = 0 skip
this property (default)
1 calculate electric potential
WHERE =
COMASS - center of mass
NUCLEI - at each nucleus (default)
POINTS - at points given in $POINTS
GRID - at grid given in $GRID
PDC - at points controlled by $PDC.
OUTPUT = PUNCH, PAPER, or BOTH (default)
This property is the electrostatic potential V(a) felt
by a test positive
charge, due to the molecular charge
density.
A nucleus at the evaluation point is ignored.
If this property
is evaluated at the nuclei, it obeys the
equation
sum on nuclei(a) Z(a)*V(a) = 2*V(nn) + V(ne).
The electronic
portion of this property is called the
diamagnetic
shielding.
==========================================================
1
$ELDENS $ELFLDG
==========================================================
$ELDENS group (not required)
This group controls electron density calculation.
IEDEN =
0 skip this property (default)
= 1 compute the electron density.
MORB
= The molecular orbital whose electron density is
to be computed. If zero, the total density is
computed. (default=0)
WHERE =
COMASS - center of mass
NUCLEI - at each nucleus (default)
POINTS - at points given in $POINTS
GRID - at grid given in $GRID
OUTPUT = PUNCH, PAPER, or BOTH (default)
IEDINT = 0 -
skip printing of integrals (default)
1 - print the electron density integrals
==========================================================
$ELFLDG group (not required)
This group controls electrostatic field and electric
field gradient
calculation.
IEFLD =
0 - skip this property (default)
1 - calculate field
2 - calculate field and gradient
WHERE =
COMASS - center of mass
NUCLEI - at each nucleus (default)
POINTS - at points given in $POINTS
OUTPUT = PUNCH, PAPER, or BOTH (default)
IEFINT = 0 -
skip printing these integrals (default)
1 - print electric field integrals
2 - also print field gradient integrals
-2 - print field gradient integrals only
The Hellman-Feynman
force on a nucleus is the nuclear
charge multiplied
by the electric field at that nucleus.
The electric
field is the gradient of the electric
potential, and
the field gradient is the hessian of the
electric potential.
The components of the electric field
gradient tensor
are formed in the conventional way, i.e.
see D.Neumann
and J.W.Moskowitz.
==========================================================
1
$POINTS $GRID
==========================================================
$POINTS group (not required)
This group is used to input points at which properties
will be computed.
This first card in the group must
contain the
string ANGS or BOHR, followed by an integer
NPOINT, the
number of points to be used. The next NPOINT
cards are read
in free format, containing the X, Y, and Z
coordinates
of each desired point.
==========================================================
$GRID group (not required)
This group is used to input a grid (plane through the
molecule) on
which properties will be calculated.
ORIGIN(i) = coordinates
of the lower left corner of
the plot.
XVEC(i)
= coordinates of the lower right corner of
the plot.
YVEC(i)
= coordinates of the upper left corner of
the plot.
SIZE
= grid increment, default is 0.25.
UNITS
= units of the above four values, it can be
either BOHR or ANGS (the default).
Note that XVEC
and YVEC are not necessarily parallel to
the X and Y
axes, rather they are the axes which you
desire to see
plotted by the MEPMAP contouring program.
==========================================================
* * * * * * * * * * * * * * * * * * * *
For conversion factors, and references
see the 'further information' section.
* * * * * * * * * * * * * * * * * * * *
1
$PDC
==========================================================
$PDC group (relevant if WHERE=PDC in $ELPOT)
This group determines the points at which to compute
the electrostatic
potential, for the purpose of fitting
atomic charges
to this potential. Constraints on the fit
which determines
these "potential determined charges" can
include the
conservation of charge, the dipole, and the
quadrupole.
PTSEL =
determines the points to be used, choose from
GEODESIC to use a set of points on several fused
sphere
van der Waals surfaces, with points
selected using an algorithm due to Mark
Spackman. The results are similar to those
from the Kollman/Singh method, but are
less rotation dependent. (default)
CONNOLLY to use a set of points on several fused
sphere
van der Waals surfaces, with points
selected using an algorithm due to Michael
Connolly. This is identical to the method
used by Kollman & Singh (see below)
CHELPG to use a modified version of the CHELPG
algorithm, which produces a symmetric
grid of points for a symmetric molecule.
CONSTR = NONE
- no fit is performed. The potential at
the points is instead output according
to OUTPUT in $ELPOT.
CHARGE - the sum of fitted atomic charges is
constrained to reproduce the total
molecular charge. (default)
DIPOLE - fitted charges are constrained to
exactly reproduce the total charge
and dipole.
QUPOLE - fitted charges are constrained to
exactly reproduce the charge, dipole,
and quadrupole.
Note: the number of constraints cannot exceed
the number of parameters, which is the number
of nuclei. Planar molecules afford fewer
constraint equations, namedly two dipole
constraints and three quadrupole constraints,
instead of three and five, repectively.
1
* * * the next 5 pertain to PTSEL=GEODESIC or CONNOLLY * * *
VDWSCL = scale
factor for the first shell of VDW spheres.
The default of 1.4 seems to be an empirical best
value. Values for VDW radii for most elements up
to Z=36 are internally stored.
VDWINC = increment
for successive shells (default = 0.2).
The defaults for VDWSCL and VDWINC will result
in points chosen on layers at 1.4, 1.6, 1.8
etc
times the VDW radii of the atoms.
LAYER = number of layers of points chosen on successive
fused sphere VDW surfaces (default = 4)
NFREQ =
flag for particular geodesic tesselation of
points. Only relevant if PTSEL=GEODESIC.
Options are:
(10*h + k) for {3,5+}h,k tesselations
-(10*h + k) for {5+,3}h,k tesselations
(of course both nh and nk must be less than 10,
so NFREQ must lie within the range -99 to
99)
The default value is NFREQ=30
(=03)
PTDENS = density of points on the surface of each scaled
VDW sphere (in points per square au).
Only relevant
if PTSEL=CONNOLLY. Default is 0.28 per
au squared,
which corresponds to 1.0 per square Angstrom,
the
default recommended by Kollman & Singh.
* * * the next two pertain to PTSEL=CHELPG * * *
RMAX
= maximum distance from any point to the closest
atom. (default=3.0 Angstroms)
DELR
= distance between points on the grid.
(default=0.8 Angstroms)
1
MAXPDC = an estimate
of the total number of points whose
electrostatic potential will be included in the
fit. (default=10000)
* * *
CENTER = an array
of coordinates at which the moments were
computed.
DPOLE = the molecular dipole.
QPOLE = the molecular quadrupole.
PDUNIT = units
for the above values. ANGS (default) will
mean that the coordinates are in Angstroms, the
dipole in Debye, and quadrupole in Buckinghams.
BOHR implies atomic units for all 3.
Note:
it is easier to compute the moments in the
current run, by setting IEMOM to at least 2 in
$ELMOM. However, you could fit experimental data,
for example, by reading it in here.
==========================================================
There is no unique way to define fitted atomic
charges.
Smaller numbers of points at which the electro-
static potential
is fit, changes in VDW radii, asymmetric
point location,
etc. all affect the results. A useful
bibliography
is
U.C.Singh, P.A.Kollman,
J.Comput.Chem. 5, 129-145(1984)
L.E.Chirlain,
M.M.Francl, J.Comput.Chem. 8, 894-905(1987)
R.J.Woods, M.Khalil,
W.Pell, S.H.Moffatt, V.H.Smith,
J.Comput.Chem. 11, 297-310(1990)
C.M.Breneman,
K.B.Wiberg, J.Comput.Chem. 11, 361-373(1990)
K.M.Merz, J.Comput.Chem.
13, 749(1992)
M.A.Spackman, J.Comput.Chem. 17, 1-18(1996)
1
$MOLGRF
==========================================================
$MOLGRF group (relevant only if you have MOLGRAPH)
This option provides an interface for viewing orbitals
through a commercial
package named MOLGRAPH, from Daikin
Industries.
Note that this option uses three disk files
which are not
defined in the GAMESS execution scripts we
provide, since
we don't use MOLGRAPH ourselves. You will
need to define
files 28, 29, 30, as generic names PRGRID,
COGRID, MOGRID,
of which the latter is passed to MOLGRAPH.
GRID3D = a flag
to generate 3D grid data.
(default is .false.).
TOTAL =
a flag to generate a total density grid data.
"Total" means the sum of the orbital densities
given by NPLT array. (default is .false.).
MESH
= numbers of grids. You can use different numbers
for three axes. (default is MESH(1)=21,21,21).
BOUND =
boundary coordinates of a 3D graphical cell. The
default is that the cell is larger than the
molecular skeleton by 3 bohr in all directions.
E.g., BOUND(1)=xmin,xmax,ymin,ymax,zmin,zmax
NPLOTS = number
of orbitals to be used to generate 3D grid
data. (default is NPLOTS=1).
NPLT
= orbital IDs. The default is 1 orbital only, the
HOMO or SOMO. If the LOCAL option is given in
$CONTRL, localized orbital IDs should be given.
For example, NPLT(1)=n1,n2,n3,...
CHECK =
debug option, printing some of the grid data.
If you are interested
in graphics, look at the WWW page
for information
about other graphics packages with GAMESS.
==========================================================
1
$STONE
==========================================================
$STONE group (optional)
This group defines the expansion points for Stone's
distributed
multipole analysis (DMA) of the electrostatic
potential.
The DMA takes the multipolar expansion of each overlap
charge density
defined by two gaussian primitives, and
translates it
from the center of charge of the overlap
density to the
nearest expansion point. Some references
for the method
are
Stone, Chem.Phys.Lett. 83, 233 (1981)
Price and Stone, Chem.Phys.Lett. 98, 419 (1983)
Buckingham and Fowler, J.Chem.Phys. 79, 6426 (1983)
Stone and Alderton, Mol.Phys. 56, 1047 (1985)
The existence of a $STONE group in the input is what
triggers the
analysis. Enter as many lines as you wish,
in any order,
terminated by a $END record.
----------------------------------------------------------
ATOM i name, where
ATOM is a keyword indicating that a particular
atom is selected as an expansion center.
i is the number of the atom
name is an optional name for the atom. If not
entered the name will be set to the name
used in the $DATA input.
----------------------------------------------------------
ATOMS
is a keyword selecting all nuclei in the
molecule as expansion points. No other
input on the line is necessary.
----------------------------------------------------------
BONDS
is a keyword selecting all bond midpoints
in the molecule as expansion points. No
other input on the line is necessary.
----------------------------------------------------------
1
$STONE
----------------------------------------------------------
BOND i j name, where
BOND is a keyword indicating that a bond mid-
point is selected as an expansion center.
i,j are the indices of the atoms defining
the
bond, corresponding to two atoms in $DATA.
name an optional name for the bond midpoint.
If omitted, it is set to 'BOND'.
----------------------------------------------------------
CMASS
is a keyword selecting the center of mass
as an expansion point. No other input on
the line is necessary.
----------------------------------------------------------
POINT x y z name, where
POINT is a keyword indicating that an arbitrary
point is selected as an expansion point.
x,y,z are the coordinates of the point, in Bohr.
name is an optional name for the expansion
point. If omitted, it is set to 'POINT'.
----------------------------------------------------------
While making
the EFPs for QM/MM run, a single keyword
QMMMBUF is necessary.
Adding additional keywords may lead
to meaningless
results. The program will automatically
select atoms
and bond midpoints which are outside the
buffer zone
as the multipole expansion points.
QMMMBUF nmo, where
QMMMBUF is a keyword specifying the number of QM/MM
buffer molecular orbitals, which must be the
first NMO orbitals in the MO set. These
orbitals must be frozen in the buffer zone,
so this is useful only if $MOFRZ is given.
NMO is the number of buffer MO-s
(if NMO is omitted, it will be set to the
number of frozen MOs in $MOFRZ)
==========================================================
The second and
third moments on the printout can be
converted to
Buckingham's tensors by formula 9 of
A.D.Buckingham, Quart.Rev. 13, 183-214 (1959)
These can in
turn be converted to spherical tensors
by the formulae
in the appendix of
S.L.Price, et al. Mol.Phys. 52, 987-1001 (1984)
1
$RAMAN
==========================================================
$RAMAN group (relevant for all SCFTYPs)
This input controls the computation of Raman intensity
by the numerical
differentiation produre of Komornicki and
others.
It is applicable to any wavefunction for which
the analytic
gradient is available, including some MP2 and
CI cases.
The calculation involves the computation of 19
nuclear gradients,
one without applied electric fields,
plus 18 no symmetry
runs with electric fields applied in
various directions.
The numerical second differencing
produces intensity
values with 2-3 digits of accuracy.
This run must follow an earlier RUNTYP=HESSIAN job,
and the $GRAD
and $HESS groups from that first job must be
given as input.
If the $DIPDR is computed analytically
by this Hessian
job, it too may be read in, if not, the
numerical Raman
job will evaluate $DIPDR. Once the data
from the 19
applied fields is available, the $ALPDR tensor
is evaluated.
Then the nuclear derivatives of the dipole
moment and alpha
polarizability will be combined with the
normal coordinate
information to produce the IR and Raman
intensity of
each mode.
To study isotopic substitution speedily, input the
$GRAD, $HESS,
$DIPDR, and $ALPDR groups along with the
desired atomic
masses in $MASS.
The code does not permit any semi-empirical or solvation
models to be
used.
EFIELD = applied
electric field strenth. The literature
suggests values in the range 0.001 to 0.005.
(default = 0.002 a.u.)
==========================================================
$ALPDR group (relevant for RUNTYP=RAMAN or HESSIAN)
Formatted alpha
derivative tensor, punched by a previous
RUNTYP=RAMAN
job. If both $DIPDR and this group are found
in the input
file, the applied field computation will be
skipped, to
immediately evaluate IR and Raman intensities.
If this group
is found during a Hessian job, the Raman
intensities
will be added to the output. You might want
to run as RUNTYP=HESSIAN
instead of RUNTYP=RAMAN in order
to have access
to PROJCT or the other options available in
the $FORCE group.
==========================================================
1
$MOROKM
==========================================================
$MOROKM group (relevant for RUNTYP=MOROKUMA)
This group controls how the supermolecule input in the
$DATA group
is divided into two or more monomers. Both
the supermolecule
and its constituent monomers must be
well described
by RHF wavefunctions.
MOROKM = a flag
to request Morokuma-Kitaura decomposition.
(default is .TRUE.)
RVS
= a flag to request "reduced variation space"
decomposition. This differs from the Morokuma
option, and one or the other or both may be
requested in the same run. (default is .FALSE.)
BSSE
= a flag to request basis set superposition error
be computed. You must ensure that CTPSPL is
selected. This option applies only to MOROKM
decompositions, as a basis superposition error is
automatically generated by the RVS scheme. This
is not the full Boys counterpoise correction, as
explained in the reference. (default is .FALSE.)
* * *
IATM
= An array giving the number of atoms in each of
the monomer. Up to ten monomers may be defined.
Your input in $DATA must have all the atoms in
the first monomer defined before the atoms in the
second monomer, before the third monomer... The
number of atoms belonging to the final monomer
can be omitted. There is no sensible default for
IATM, so don't omit it from your input.
ICHM
= An array giving the charges of the each monomer.
The charge of the final monomer may be omitted,
as it is fixed by ICH in $CONTRL, which is the
total charge of the supermolecule. The default
is neutral monomers, ICHM(1)=0,0,0,...
EQUM
= a flag to indicate all monomers are equivalent
by symmetry (in addition to containing identical
atoms). If so, which is not often true, then only
the unique computations will be done.
(default is .FALSE.)
CTPSPL = a flag
to decompose the interaction energy into
charge transfer plus polarization terms. This
is most appropriate for weakly interacting
monomers. (default is .TRUE.)
1
CTPLX =
a flag to combine the CT and POL terms into a
single term. If you select this, you might want
to turn CTPSPL off to avoid the extra work that
that decomposition entails, or you can analyze
both ways in the same run (default=.FALSE.)
RDENG =
a flag to enable restarting, by reading the
lines containing "FINAL ENERGY" from a previous
run. The $ENERGY group is single lines read
under format A16,F20.10 containing the E, and a
card $END to complete. The 16 chars = anything.
(default is .FALSE.)
==========================================================
The present implementation has some quirks:
1. The initial
guess of the monomer orbitals is not
controlled by $GUESS. The program first looks for a
$VEC1, $VEC2, ... group for each monomer. If they
are found, they will be MOREAD. If any of these are
missing, the guess for that monomer will be constructed
by HCORE. Check your monomer energies carefully! The
initial guess orbitals for the supermolecule are formed
by a block diagonal matrix of the monomer orbitals.
2. The use of
symmetry is turned off internally.
3. There is
no direct SCF option. File ORDINT will be a
full C1 list of integrals. File AOINTS will contain
whatever subset of these is needed for each particular
decomposition step. So extra disk space is needed
compared to RUNTYP=ENERGY.
4. This kind
of run applies only to ab initio cases.
5. This kind
of run will work in parallel.
6. Spherical
harmonics may not be used.
References:
C.Coulson
in "Hydrogen Bonding", D.Hadzi, H.W.Thompson,
Eds., Pergamon Press, NY, 1957, pp 339-360.
C.Coulson
Research, 10, 149-159 (1957).
K.Morokuma
J.Chem.Phys. 55, 1236-44 (1971).
K.Kitaura, K.Morokuma
Int.J.Quantum Chem. 10, 325 (1976).
K.Morokuma,
K.Kitaura in "Chemical Applications of
Electrostatic Potentials", P.Politzer,D.G.Truhlar, Eds.
Plenum Press, NY, 1981, pp 215-242.
The method coded
is the newer version described in the
latter two papers.
Note that the CT term is computed
separately for
each monomer, as described in the words
below equation
16 of the 1981 paper, not simultaneously.
Reduced Variational
Space:
W.J.Stevens,
W.H.Fink, Chem.Phys.Lett. 139, 15-22(1987).
1
A comparison
of the RVS and Morokuma decompositions can
be found in
the review article: "Wavefunctions and
Chemical Bonding"
M.S.Gordon, J.H.Jensen in "Encyclopedia
of Computational
Chemistry", volume 5, P.V.R.Schleyer,
editor, John
Wiley and Sons, Chichester, 1998.
BSSE during Morokuma
decomposition:
R.Cammi, R.Bonaccorsi,
J.Tomasi
Theoret.Chim.Acta
68, 271-283(1985).
The present implementation:
"Energy decomposition
analysis for many-body interactions,
and application
to water complexes"
W.Chen, M.S.Gordon
J.Phys.Chem. 100, 14316-14328(1996)
1
$FFCALC
==========================================================
$FFCALC group (relevant for RUNTYP=FFIELD)
This group permits the study of the influence of an
applied electric
field on the wavefunction. The most
common finite
field calculation applies a sequence of
fields to extract
the linear polarizability and first and
second order
hyperpolarizability. The method is general,
and so works
for all ab initio wavefunctions in GAMESS.
EFIELD
= applied electric field strength
(default=0.001 a.u.)
IAXIS and JAXIS
specify the orientation of the applied
field. 1,2,3 mean x,y,z respectively.
The default is IAXIS=3 and JAXIS=0.
If
IAXIS=i and JAXIS=0, the program computes alpha(ii),
beta(iii), and gamma(iiii) from the energy changes, and
a few more components from the dipole changes. Five
wavefunction evaluations are performed.
If
IAXIS=i and JAXIS=j, the program computes the cross
terms beta(ijj), beta(iij), and gamma(iijj) from the
energy changes, and a few more components from the
dipole changes. This requires nine evaluations of the
wavefunction.
AOFF
= a flag to permit evaluation of alpha(ij)
when the dipole moment is not available.
This is necessary only for MP2, and means
the off-axial calculation will do 13, not
9 energy evaluations. Default=.FALSE.
SYM
= a flag to specify when the fields to be
applied along the IAXIS and/or JAXIS (or
according to EONE below) do not break the
molecular symmetry. Since most fields do
break symmetry, the default is .FALSE.
ONEFLD
= a flag to specify a single applied field
calculation will be performed. Only the
energy and dipole moment under this field
are computed. If this option is selected,
only SYM and EONE input is heeded. The
default is .FALSE.
EONE
= an array of the three x,y,z components of
the single applied field.
There are notes on RUNTYP=FFIELD on the next page.
1
$FFCALC
Finite field calculations require large basis sets,
and extraordinary
accuracy in the wavefunction. To
converge the
SCF to many digits is sometimes problematic,
but we suggest
you use the input to increase integral
accuracy and
wavefunction convergence, for example
$CONTRL ICUT=20 ITOL=30 INTTYP=HONDO $END
$SCF CONV=1.0d-10 FDIFF=.FALSE. $END
In many cases, the applied fields will destroy the
molecular symmetry.
This means the integrals are
calculated once
with point group symmetry to do the
initial field
free wavefunction evaluation, and then again
with point group
symmetry turned off. If the fields
applied do not
destroy symmetry, you can avoid this second
calculation
of the integrals by SYM=.TRUE. This option
also permits
use of symmetry during the applied field
wavefunction
evaluations.
Examples of fields that do not break symmetry are a
Z-axis field
for an axial point group which is not
centrosymmetric
(i.e. C2v). However, a second field in
the X or Y direction
does break the C2v symmetry.
Application
of a Z-axis field for benzene breaks D6h
symmetry.
However, you could enter the group as C6v in
$DATA while
using D6h coordinates, and regain the prospect
of using SYM=.TRUE.
If you wanted to go on to apply a
second field
for benzene in the X direction, you might
want to enter
Cs in $DATA, which will necessitate the
input of two
more carbon and hydrogen atom, but recovers
use of SYM=.TRUE.
Reference: H.A.Kurtz,
J.J.P.Stewart, K.M.Dieter
J.Comput.Chem. 11, 82-87 (1990).
For analytic computation of static and also frequency
dependent NLO
proerties, for closed shell cases, see the
$TDHF group.
==========================================================
1
$TDHF
==========================================================
$TDHF group (relevant for SCFTYP=RHF if RUNTYP=TDHF)
This group permits the analytic calculation of various
static and/or
frequency dependent polarizabilities, with
an emphasis
on important NLO properties such as second and
third harmonic
generation. The method is programmed only
for closed shell
wavefunctions, at the semi-empirical or
ab initio level.
Ab initio calculations may be direct SCF,
or parallel,
if desired.
Because the Fock matrices computed during the time-
dependent Hartree-Fock
CPHF are not symmetric, you may not
use symmetry.
You must enter NOSYM=1 in $CONTRL!
For a more general numerical approach to the static
properties,
see $FFCALC.
NFREQ = Number of frequencies to be used. (default=1)
FREQ
= An array of energy values in atomic units. For
example: if NFREQ=3 then FREQ(1)=0.0,0.1,0.25.
By default, only the static polarizabilities are
computed. (default is freq(1)=0.0)
The conversion factor from Hartree to wave
numbers is 219,474.6, and the wavelength is
given (in nm) by 45.56/FREQ.
MAXITA = Maximum
number of iterations for an alpha
computation. (default=100)
MAXITU = Maximum
number of iterations in the second order
correction calculation. This applies to iterative
beta values and all gammas. (default=100)
ATOL
= Tolerance for convergence of first-order results.
(default=1.0d-05)
BTOL
= Tolerance for convergence of second-order results.
(default=1.0d-05)
RETDHF = a flag
to choose starting points for iterative
calculations from best previous results.
(default=.true.)
1
* * * the following NLO properties are available * * *
INIB
= 0 turns off all beta computation (default)
= 1 calculates only noniterative beta
= 2 calculate iterative and noniterative beta
The next flags allow further BETA tuning
BSHG = Calculate beta for second harmonic generation.
BEOPE = Calculate beta for electrooptic Pockels effect.
BOR = Calculate beta for optical rectification.
INIG
= 0 turns off all gamma computation (default)
= 1 calculates only noniterative gamma
= 2 calculate iterative and noniterative gamma
The next flags allow further GAMMA tuning
GTHG = Calculate gamma for third harmonic generation.
GEFISH = Calculate
gamma for electric-field induced
second harmonic generation.
GIDRI =
Calculate gamma for intensity dependent
refractive index.
GOKE = Calculate gamma for optical Kerr effect.
These will be computed only if a nonzero energy is
requested.
The default for each flag is .TRUE., and they
may be turned
off individually by setting some .FALSE.
Note however
that the program determines the best way to
calculate them.
For example, if you wish to have the SHG
results but
no gamma results are needed, the SHG beta will
be computed
in a non-iterative way from alpha(w) and
alpha(2w).
However if you request the computation of the
THG gamma, the
second order U(w,w) results are needed and
an iterative
SHG calculation will be performed whether
you request
it or not, as it is a required intermediate.
Reference:
S.P.Karna, M.Dupuis
J.Comput.Chem. 12, 487-504 (1991).
P.Korambath,
H.A.Kurtz, in "Nonlinear Optical Materials",
ACS Symposium
Series 628, S.P.Karna and A.T.Yeates, Eds.
pp 133-144,
Washington DC, 1996.
Review: D.P.Shelton, J.E.Rice, Chem.Rev. 94, 3-29(1994).
==========================================================
1
$EFRAG
==========================================================
$EFRAG group (optional)
This group gives the name and position of one or more
effective fragment
potentials. It consists of a series of
free format
card images, which may not be combined onto a
single line!
The position of a fragment is defined by
giving any three
points within the fragment, relative to
the ab initio
system defined in $DATA, since the effective
fragments have
a frozen internal geometry. All other
atoms within
the fragment are defined by information in
the $FRAGNAME
group.
----------------------------------------------------------
-1- a line containing one or more of these options:
COORD =CART selects use of Cartesians
coords
to define the fragment position at
line -3-. (default)
=INT selects use of Z-matrix internal
coordinates at line -3-.
POLMETHD=SCF indicates the induced dipole
for
each fragment due to the ab initio
electric field and other fragment
fields is updated only once during
each SCF iteration.
=FRGSCF requests microiterations during
each SCF iteration to make induced
dipoles due to ab initio and other
fragment fields self consistent
amoung the fragments. (default)
Both methods converge to the same
dipolar interaction.
POSITION=OPTIMIZE Allows full optimization within the
ab initio part, and optimization of
the rotational and translational
motions of each fragment. (default)
=FIXED Allows full optimization of the
ab initio system, but freezes the
position of the fragments. This
makes sense only with two or more
fragments, as what is frozen is the
fragments' relative orientation.
=EFOPT the same as OPTIMIZE, but if the
fragment gradient is large, up to
5 geometry steps in which only the
fragments move may occur, before
the geometry of the ab initio piece
is relaxed. This may save time by
reusing the two electron integrals
for the ab initio system.
1
$EFRAG
NBUFFMO = n First n orbitals in the
MO matrix
are deemed to belong to the QM/MM
buffer and will be excluded from
the interaction with the EFP region.
This makes sense only if these first
MOs are frozen via the $MOFRZ group.
Input a blank
line if all the defaults are acceptable.
----------------------------------------------------------
-2- FRAGNAME=XXX
XXX is the name
of the fragment whose coordinates are to
be given next.
All other information defining the
fragment is
given in a supplemental $XXX group, which is
referred to
below as a $FRAGNAME group.
A RHF/DZP EFP
for water is internally stored in GAMESS.
Choose FRAGNAME=H2OEF2
to look up this numerical data,
and then skip
the input of $H2OEF2 and $FRGRPL groups.
----------------------------------------------------------
-3-
NAME, X, Y, Z
(COORD=CART)
NAME, I, DISTANCE, J, BEND, K, TORSION (COORD=INT)
NAME
= the name of a fragment point. The name used
here must match one of the points in $FRAGNAME.
X, Y, Z
= Cartesian coordinates defining the position of
this fragment point RELATIVE TO THE COORDINATE
ORIGIN used in $DATA. The choice of units is
controlled by UNITS in $CONTRL.
I, DISTANCE,
J, BEND, K, TORSION = the usual Z-matrix
connectivity internal coordinate definition.
The atoms I, J, K must be atoms in the ab
initio system from in $DATA, or fragment points
already defined in the current fragment or
previously defined fragments.
Line -3- must
be given a total of three times to define
this fragment's
position.
----------------------------------------------------------
Repeat lines
-2- and -3- to enter as many fragments as you
desire, and
then end the group with a $END line.
Note that it
is quite typical to repeat the same fragment
name at line
-2-, to use the same fragment system at many
different positions.
==========================================================
* * * * * * * * * * * * * * * * * * * * *
For tips on effective fragment potentials
see the 'further information' section
* * * * * * * * * * * * * * * * * * * * *
1
$FRAGNAME
==========================================================
(required for each FRAGNAME given in $EFRAG)
$FRAGNAME group
This group gives all pertinent information for a given
effective fragment
potential (EFP). This information
falls into three
categories:
electrostatic (distributed multipoles, screening)
distributed polarizabilities
exchange repulsion
It is input
using several different subgroups, which
should be given
in the order shown below. Each subgroup
is specified
by a particular name, and is terminated by
the word STOP.
You may omit any of the subgroups to omit
that term from
the EFP. All values are given in atomic
units.
To input monopoles,
follow input sequence -EM-
To input dipoles,
follow input sequence -ED-
To input quadrupoles,
follow input sequence -EQ-
To input octupoles,
follow input sequence -EO-
To input screening
parameters, follow input sequence -ES-
To input polarizable
points, follow input sequence -P-
To input repulsive
points, follow input sequence -R-
----------------------------------------------------------
-1-
a single descriptive title card
----------------------------------------------------------
-2- COORDINATES
COORDINATES signals
the start of the subgroup containing
the multipolar
expansion terms (charges, dipoles, ...).
Optionally,
one can also give the coordinates of the
polarizable
points, or centers of exchange repulsion.
-3- NAME, X, Y, Z, WEIGHT, ZNUC
NAME is a unique
string identifying the point.
X, Y, Z are
the Cartesian coordinates of the point.
WEIGHT and ZNUC
are the atomic mass and nuclear charge,
and are given
only for the points which are nuclei.
Typically the
true nuclei will appear twice, once for
defining the
positive nuclear charge and its screening,
and a second
time for defining the electronic distributed
multipoles.
Repeat line -3-
for each expansion point, and terminate
the list with
a "STOP".
----------------------------------------------------------
1
$FRAGNAME
-EM1-
MONOPOLES
MONOPOLES signals
the start of the subgroup containing
the electronic
and nuclear monopoles.
-EM2- NAME, CHARGE1, CHARGE2
NAME must match
one given in the COORDINATES subgroup.
CHARGE1 = electronic
monopole at this point.
CHARGE2 = nuclear
monopole at this point. Omit or enter
zero if this is a bond midpoint or some other
expansion point that is not a nucleus.
Repeat -EM2-
to define all desired charges.
Terminate this
subgroup with a "STOP".
----------------------------------------------------------
-ED1-
DIPOLES
DIPOLES signals
the start of the subgroup containing the
dipolar part
of the multipolar expansion.
-ED2- NAME, MUX, MUY, MUZ
NAME must match
one given in the COORDINATES subgroup.
MUX, MUY, MUZ
are the components of the electronic dipole.
Repeat -ED2-
to define all desired dipoles.
Terminate this
subgroup with a "STOP".
----------------------------------------------------------
-EQ1-
QUADRUPOLES
QUADRUPOLES signals
the start of the subgroup containing
the quadrupolar
part of the multipolar expansion.
-EQ2- NAME, XX, YY, ZZ, XY, XZ, YZ
NAME must match
one given in the COORDINATES subgroup.
XX, YY, ZZ,
XY, XZ, and YZ are the components of the
electronic quadrupole
moment.
Repeat -EQ2-
to define all desired quadrupoles.
Terminate this
subgroup with a "STOP".
----------------------------------------------------------
-EO1-
OCTUPOLES
OCTUPOLES signals
the start of the subgroup containing
the octupolar
part of the multipolar expansion.
-EO2- NAME,
XXX, YYY, ZZZ, XXY, XXZ,
XYY, YYZ, XZZ, YZZ, XYZ
NAME must match
one given in the COORDINATES subgroup.
XXX, ...
are the components of the electronic octupole.
Repeat -EO2-
to define all desired octupoles.
Terminate this
subgroup with a "STOP".
----------------------------------------------------------
1
$FRAGNAME
-ES1- SCREEN
SCREEN signals
the start of the subgroup containing the
screening terms
(A*exp[-B*r**2]) for the distributed
multipoles,
which account for charge penetration effects.
-ES2- NAME, A, B
NAME must match
one given in the COORDINATES subgroup.
A, B are the
parameters of the Gaussian screening term.
Repeat -ES2-
to define all desired screening points.
Terminate this
subgroup with a "STOP".
----------------------------------------------------------
-P1- POLARIZABLE POINTS
POLARIZABLE POINTS
signals the start of the subgroup
containing the
distributed polarizability tensors, and
their coordinates.
-P2- NAME, X, Y, Z
NAME gives a
unique identifier to the location of this
polarizability
tensor. It might match one of the points
already defined
in the COORDINATES subgroup, but often
does not.
Typically the distributed polarizability
tensors are
located at the centroids of localized MOs.
X, Y, Z are the
coordinates of the polarizability point.
They should
be omitted if NAME did appear in COORDINATES.
The units are
controlled by UNITS= in $CONTRL.
-P3- XX, YY, ZZ, XY, XZ, YZ, YX, ZX, ZY
XX, ... are components
of the distributed polarizability,
which is not
a symmetric tensor. XY means dMUx/dFy, where
MUx is a dipole
component, and Fy is a component of an
applied field.
Repeat -P2- and
-P3- to define all desired polarizability
tensors, and
terminate this subgroup with a "STOP".
----------------------------------------------------------
1
$FRAGNAME
-R1- REPULSIVE
POTENTIAL
REPULSIVE POTENTIAL
signals the start of the subgroup
containing the
fitted exchange repulsion potential, for
the interaction
between the fragment and the ab initio
part of the
system. This term also accounts for charge
transfer effects.
The term has the form
N
sum C * exp[-D * r**2]
i i i
-R2- NAME, X, Y, Z, N
NAME may match
one given in the COORDINATES subgroup,
but need not.
If NAME does not match one of the
known points,
you must give its coordinates X, Y, and
Z, otherwise
omit these three values. N is the total
number of terms
in the fitted repulsive potential.
-R3- C, D
These two values
define the i-th term in the repulsive
potential.
Repeat line -R3- for all N terms.
Repeat -R2- and
-R3- to define all desired repulsive
potentials,
and terminate this subgroup with a "STOP".
==========================================================
The entire $FRAGNAME group is terminated by a " $END".
1
$FRGRPL
==========================================================
$FRGRPL group
This group defines
the inter-fragment repulsive potential,
which consists
primarily of exchange repulsions but also
includes charge
transfer. Note that the functional form
used for the
fragment-fragment repulsion differs from
that used for
the ab initio-fragment repulsion, which is
defined in the
$FRAGNAME group. The form of the potential
is
N
sum A * exp[-B * r]
i i i
----------------------------------------------------------
-1- PAIR=FRAG1 FRAG2
specifies which
two fragment repulsions are being defined.
$FRAGNAME input
for the two names FRAG1 and FRAG2 must
have been given.
----------------------------------------------------------
-2- NAME1
NAME2 A B
*or*
NAME1 NAME2 'EQ' NAME3 NAME4
NAME1 must be
one of the "NAME" points defined in the
$FRAG1 group's
REPULSION POTENTIAL section. Similarly
NAME2 must be
a point from the $FRAG2 group. In addition,
NAME1 or NAME2
could be the keyword CENTER, indicating the
center of mass
of the fragment.
A and B are the
parameters of the fitted repulsive
potential.
The second form
of the input allows equal potential fits
to be used.
The syntax implies that the potential between
the points NAME1
and NAME2 should be taken the same as the
potential previously
given in this group for the pair of
points NAME3
and NAME4.
If there are
NPT1 points in FRAG1, and NPT2 points in
FRAG2, input
line -2- should be repeated NPT1*NPT2 times.
Terminate the
pairs of potentials with a "STOP" card.
Any pairs which
you omit will be set to zero interaction.
Typically the
number of points on which fitted potentials
might be taken
to be all the nuclei in a fragment, plus
the center of
mass.
----------------------------------------------------------
Repeat lines
-1- and -2- for all pairs of fragments, then
terminate the
group with a $END line.
==========================================================
1
$PCM
==========================================================
$PCM group
(optional)
This group controls solvent effect computations using
the Polarizable
Continuum Method. If this group is found
in the input
file, a PCM computation is performed. The
default calculation,
chosen by selecting only the SOLVNT
keyword, is
to compute the electrostatic free energy.
Appropriate
numerical constants are provided for a wide
range of solvents.
Additional keywords allow for more
sophisticated
computations, namely cavitation, repulsion,
and dispersion
free energies. The methodology for these
is general,
but only numerical constants for water are
provided.
There is additional information on PCM in the
References chapter
of this manual.
PCM is programmed only for RHF and MCSCF wavefunctions.
Geometry optimization with PCM will probably not be able
to converge
to the default OPTTOL in $STATPT due to some
numerical inaccuracies.
You will likely need to raise this
value if you
attempt solvent optimizations, by a factor of
two to five.
--- the first
set of parameters controls the computation:
IEF, ICOMP, ICAV, IDISP, IREP, IDP, and IFIELD.
IEF
switch to choose boundary element method or the
integral equation formula as the PCM solver.
= 0 isotropic dielectrics using BEM PCM
= 1 anisotropic dielectrics using IEF PCM, see $IEFPCM
= 2 ionic solutions using IEF PCM, see $IEFPCM
= 3 isotropic dielectrics using IEF PCM (default)
*** at the present time, there is a bug with IEF=1 or 2.
ICOMP =
Compensation procedure for induced charges.
Gradient runs require ICOMP be 0 or 2 only.
= 0 No. (default)
= 1 Yes, each charge is corrected in proportion
to the area of the tessera to which it belongs.
= 2 Yes, using the same factor for all tesserae.
= 3 Yes, with explicit consideration of the
portion of solute electronic charge outside
the cavity, by the method of Mennucci and
Tomasi. See the $NEWCAV group.
The behaviour
of PCM prior to Oct. 2000 can be recovered
by selecting
IEF=0 and ICOMP=2. Options IEF=1 or 2 are
incompatible
with gradients and also must choose ICOMP=0.
IEF=3 may not
choose ICOMP=3, but if diffuse functions
are in use,
this default choice may benefit from ICOMP=2.
The BEM method
(IEF=0) should normally choose ICOMP=2.
1
$PCM
ICAV
= At the end of the run, calculate the cavitation
energy, by the method of Pierotti and Claverie:
= 0 skip the computation (default)
= 1 perform the computation.
If ICAV=1, the following parameter is relevant:
TABS
= the absolute temperature, in units K.
(default=298.0)
There are two procedures for the calculation
of the repulsion and dispersion free energy.
IDISP is incompatible with IREP and IDP.
IDISP =
Calculation of both dispersion and repulsion
free energy through the empirical method of
Floris and Tomasi.
= 0 skip the computation (default)
= 1 perform the computation. See $DISREP group.
The
next two options add repulsive and dispersive terms
to the solute hamiltonian, in an ab initio manner, by
the method of Amovilli and Mennucci.
IREP
= Calculation of repulsion free energy
= 0 skip the computation (default)
= 1 perform the computation. See $NEWCAV group.
IDP
= Calculation of dispersion free energy
= 0 skip the computation (default)
= 1 perform the computation. See $DISBS group.
If
IDP=1, then three additional parameters must be
defined. The two solvent values correspond to water,
and therefore these must be input for other solvents.
WA
= solute average transition energy. This is
computed from the orbital energies for RHF,
but must be input for MCSCF runs.
(default=1.10)
WB
= ionization potential of solvent, in Hartrees.
(default=0.451)
ETA2
= square of the zero frequency refractive index
of the solvent. (default=1.75)
IFIELD = At run
end, calculate the electric potential
and electric field generated by the apparent
surface charges.
= 0 skip the computation (default)
= 1 on nuclei
= 2 on a planar grid
1
$PCM
If IFIELD=2, the following data must be input:
AXYZ,BXYZ,CXYZ
= each defines three components of the
vertices of the plane where the reaction
field is to be computed (in Angstroms)
A ===> higher left corner of the grid
B ===> lower left corner of the grid
C ===> higher right corner of the grid
NAB = vertical
subdivision (A--B edge) of the grid
NAC = horizontal
subdivision (A--C edge) of the grid.
IPRINT = 0 normal
printing (default)
= 1 turns on debugging printout
--- the next group of keywords defines the solvent
SOLVNT = keyword
naming the solvent of choice. The eight
numerical constants defining the solvent are
internally stored for the following:
WATER (or H2O)
CH3OH
C2H5OH
CLFORM (or CHCl3)
CTCL (or CCl4)
METHYCL (or CH2Cl2) 12DCLET (or
C2H4Cl2)
BENZENE (or C6H6)
TOLUENE (or C6H5CH3)
CLBENZ (or C6H5Cl) NITMET
(or CH3NO2)
NEPTANE (or C7H16) CYCHEX
(or C6H12)
ANILINE (or C6H5NH2) ACETONE (or CH3COCH3)
THF
DMSO (or DMETSOX)
The default solvent name is
INPUT
which indicates you will specify your solvent by
giving the following 8 numerical values:
RSOLV =
the solvent radius, in units Angstrom
EPS
= the dielectric constant
EPSINF = the
dielectric constant at infinite frequency.
This value must be given only for RUNTYP=TDHF,
if the external field frequency is in the optical
range and the solvent is polar; in this case the
solvent response is described by the electronic
part of its polarization. Hence the value of the
dielectric constant to be used is that evaluated
at infinite frequency, not the static one (EPS).
For nonpolar solvents, the difference between
the two is almost negligible.
TCE
= the thermal expansion coefficient, in units 1/K
VMOL
= the molar volume, in units ml/mol
STEN
= the surface tension, in units dyne/cm
DSTEN
= the thermal coefficient of log(STEN)
CMF
= the cavity microscopic coefficient
Values for TCE,
VMOL, STEN, DSTEN, CMF need to be given
only for the
case ICAV=1. Input of any or all of these
values will
override the internally stored value.
1
$PCM $PCMCAV
--- the next set of keywords defines the molecular cavity
NESFP =
the number of initial spheres.
(default = number of atoms in solute molecule)
ICENT =
option for definition of initial spheres.
= 0 centers spheres on each nucleus. (default)
= 1 sphere centers XE, YE, ZE and radii RIN will be
specified explicitly in $PCMCAV.
The cavity generation algorithm may use additional
spheres to smooth out sharp grooves, etc. The
following parameters control how many extra spheres
are generated:
OMEGA and FRO
= GEPOL parameters for the creation of the
`added spheres' defining the solvent accessible
surface. When an excessive number of spheres is
created, which may cause problems of convergence,
the value of OMEGA and/or FRO must be increased.
For example, OMEGA from 40 to 50 ... up to 90,
FRO from 0.2 ... up to 0.7.
(defaults are OMEGA=40.0, FRO=0.7)
RET
= minimum radius (in A) of the added spheres.
Increasing RET decreases the number of added
spheres. A value of 100.0 inhibits the addition
of spheres. (default=0.2)
==========================================================
$PCMCAV group (optional)
This group controls generation of the cavity holding
the solute during
Polarizable Continuum Method runs.
The cavity is
a union of spheres, according to ICENT and
associated input
values given in $PCM. The data given
here must be
given in Angstrom units.
XE,YE,ZE = arrays
giving the coordinates of the spheres.
if ICENT=0, the atomic positions will be used.
if ICENT=1, you must supply NESFP values here.
RIN = an array
giving the sphere radii.
if ICENT=0, the program will look up the internally
stored van der Waals radius for: H,He,
B,C,N,O,F,Ne, Na,Al,Si,P,S,Cl,Ar,
K,As,Se,Br,Kr, Rb,Sb,Te,I, Cs,Bi
Data for other elements is not tabulated.
if ICENT=1, give NESFP values.
1
$PCMCAV $NEWCAV
ALPHA = an array
of scaling factors, for the definition of
the solvent accessible surface. If only the first
value is given, all radii are scaled by the same
factor. (default is ALPHA(1)=1.2)
Example: Suppose
the 4th atom in your molecule is Fe, but
all other atoms have van der Waals radii. You
decide a good guess for Fe is twice the covalent
radius: $PCMCAV RIN(4)=2.33 $END
The source for
the van der Waals radii is "The Elements",
2nd Ed., John
Emsley, Clarendon Press, Oxford, 1991,
except that
for C,N,O, the U.Pisa's experience with the
best radii for
PCM treatment of singly bonded C,N,O atoms
is used instead.
The radii for a few transition metals
are given by
A.Bondi, J.Phys.Chem. 68, 441-451(1968).
==========================================================
$NEWCAV group (optional)
This group controls generation of the "escaped charge"
cavity, used
when ICOMP=3 or IREP=1 in $PCM. This cavity
is used only
to calculate the fraction of the solute
electronic charge
escapes from the original cavity.
IPTYPE = choice
for tessalation of the cavity's spheres.
= 1 uses a tetrahedron
= 2 uses a pentakisdodecahedron (default)
ITSNUM = m, the
number of tessera to use on each sphere.
if IPTYPE=1, input m=30*(n**2), with n=1,2,3 or 4
if IPTYPE=2, input m=60*(n**2), with n=1,2,3 or 4
(default is 60)
*** the next three parameters pertain to IREP=1 ***
RHOW = density, relative to liquid water (default = 1.0)
PM = molecular weight (default = 18.0)
NEVAL = number of valence electrons on solute (default=8)
The defaults
for RHOW, PM, and NEVAL correspond to water,
and therefore
must be correctly input for other solvents.
==========================================================
1
==========================================================
$IEFPCM group
(optional)
This group defines data for the integral equation
formalism version
of PCM solvation. It includes special
options for
ionic or anisotropic solutions.
The next two
sets are relevant only for anisotropic
solvents, namely
IEF=1:
EPS1, EPS2, EPS3
=
diagonal values of the dielectric permittivity
tensor with respect to the laboratory frame.
The default is EPS in $PCM
EUPHI, EUTHE,
EUPSI =
Eulerian angles which give the rotation of the
solvent orientation with respect to the lab frame.
The term lab frame means $DATA orientation.
The default for each is zero degrees.
The next two are relevant to ionic solvents, namely IEF=2:
EPSI = the ionic
solutions's dielectric, the default is
EPS from $PCM.
DISM = the ionic
strength, in Molar units (mol/dm**3)
The default is 0.0
==========================================================
1
$DISBS
==========================================================
$DISBS group
(optional)
This group defines auxiliary basis functions used to
evaluate the
dispersion free energy by the method of
Amovilli and
Mennucci. These functions are used only for
the dispersion
calculation, and thus have nothing to do
with the normal
basis given in $BASIS or $DATA. If the
input group
is omitted, only the normal basis is used for
the IDP=1 dispersion
energy.
NADD = the number of added shells
XYZE
= an array giving the x,y,z coordinates (in bohr)
of the center, and exponent of the added shell,
for each of the NADD shells.
NKTYPE = an array giving the angular momenta of the shells
An example placing 2s,2p,2d,1f on one particular atom,
$DISBS
NADD=7 NKTYP(1)= 0 0 1 1 2 2 3
XYZE(1)=2.9281086 0.0 .0001726 0.2
2.9281086 0.0 .0001726 0.05
2.9281086 0.0 .0001726 0.2
2.9281086 0.0 .0001726 0.05
2.9281086 0.0 .0001726 0.75
2.9281086 0.0 .0001726 0.2
2.9281086 0.0 .0001726 0.2 $END
==========================================================
1
$DISREP
==========================================================
$DISREP group
(optional)
This group controls evaluation of the dispersion and
repulsion energies
by the empirical method of Floris and
Tomasi.
The group must be given with IDISP=1 in $PCM.
The two options
are controlled by ICLAV and ILJ, only one
of which should
be selected.
ICLAV = selects
Claverie's disp-rep formalism.
= 0 skip computation.
= 1 Compute the solute-solvent disp-rep interaction
as a sum over atom-atom interactions through a
Buckingham-type formula (R^-6 for dispersion,
exp for repulsion). (default)
Ref: Pertsin-Kitaigorodsky "The atom-atom
potential method", page 146.
ILJ
= selects a Lennard-Jones formalism.
= 0 skip computation. (default)
= 1 solute atom's-solvent molecule interaction is
modeled by Lennard-Jones type potentials, R^-6
for dispersion, R^-12 for repulsion).
---- the following data must given for ICLAV=1:
RHO
= solvent numeral density
N
= number of atom types in the solvent molecule
NT
= an array of the number of atoms of each type in a
solvent molecule
RDIFF = distances
between the first atoms of each type
and the cavity
DK,RW = parameters
of atom-atom interactions
The defaults
are chosen for water,
RHO=3.348D-02
N=2
NT(1)=2,1
RDIFF(1)=1.20,1.50
DKT(1)=1.0,1.36
RWT(1)=1.2,1.5
---- the following data must given for ILJ=1:
RHO
= solvent numeral density
EPSI =
an array of energy constants referred to each atom
of the solute molecule.
SIGMA = an array
of typical distances, relative to each
solute atom
==========================================================
1
$COSGMS $SCRF
==========================================================
$COSGMS group (optional)
The presence of this group in the input turns on the
use of the conductor-like
screening model with molecular
shaped cavity
for RHF and closed shell MP2. For RHF, the
energy and gradient
can be computed, while MP2 is limited
to the energy
only.
EPSI
= the dielectric constant, 80 is often used for H2O
This parameter must be given.
RSOLV =
the multiplicative factor for the van der Waals
radius used for cavity construction.
(default=1.2)
NSPA
= the number of surface points on each atomic
sphere that form the cavity. (default=92)
Additional information on the COSMO model can be
found in the References chapter of this manual.
==========================================================
$SCRF group (optional)
The presence of this group in the input turns on the
use of the Kirkwood-Onsager
spherical cavity model for the
study of solvent
effects. The method is implemented for
RHF, UHF, ROHF,
GVB and MCSCF wavefunctions and gradients,
and so can be
used with any RUNTYP involving the gradient.
The method is
not implemented for MP2, CI, any of the
semiempirical
models, or for analytic hessians.
DIELEC = the dielectric constant, 80 is often used for H2O
RADIUS = the spherical cavity radius, in Angstroms
G
= the proportionality constant relating the solute
molecule's dipole to the strength of the reaction
field. Since G can be calculated from DIELEC and
RADIUS, do not give G if they were given.
Additional information on the SCRF model can be
found in the References chapter of this manual.
==========================================================
1
$ECP
==========================================================
$ECP group
(required if ECP=READ in $CONTRL)
This group lets you read in effective core potentials,
for some or
all of the atoms in the molecule. You can
use built in
potentials for some of the atoms if you like.
This is a free
format (positional) input group.
*** Give a card set -1-, -2-, and -3- for each atom ***
-card 1- PNAME, PTYPE, IZCORE, LMAX+1
PNAME is a 8
character descriptive tag for this potential.
If it is repeated for a subsequent atom, no other
information need be given on this card, and cards
-2- and -3- may also be skipped. The information
will be copied from the first atom by this PNAME.
Do not use the option to repeat the previously read
ECP for an atom with PTYPE=NONE, instead type "none".
PTYPE = GEN
a general potential should be read.
= SBKJC look up the Stevens/Basch/Krauss/Jasien/
Cundari potential for this type of atom.
= HW look up the Hay/Wadt built in potential
for this type of atom.
= NONE treat all electrons on this atom.
IZCORE is the
number of core electrons to be removed.
LMAX
is the maximum angular momentum occupied in the
core orbitals being removed (usually). Give
IZCORE and LMAX only if PTYPE is GEN.
*** For the first
occurence of PNAME, if PTYPE is GEN, ***
*** then give
cards -2- and -3-. Otherwise go to -1-. ***
*** Card sets -2- and -3- are repeated LMAX+1 times ***
The potential U(LMAX+1) is given first,
followed by U(L)-U(LMAX+1), for L=1,LMAX.
-card 2- NGPOT
NGPOT is the
number of Gaussians in this part of the
local effective potential.
-card 3- CLP,NLP,ZLP (repeat this card NGPOT times)
CLP is the coefficient
of this Gaussian in the potential.
NLP is the power
of r for this Gaussian.
ZLP is the exponent
of this Gaussian.
* * *
By far the easiest
way to use the SBKJC potential for all
atoms in the
formic acid molecule is to request ECP=SBKJC
in $CONTRL.
But the next page shows two alternatives.
1
$ECP
The first way
is to look up the program's internally
stored SBKJC
potentials one atom at a time:
$ECP
C-ECP SBKJC
H-ECP NONE
O-ECP SBKJC
O-ECP
H-ECP NONE
$END
The second oxygen
duplicates the first, no core electrons
are removed
for hydrogen. The order of the atoms must
follow that
generated by $DATA. Note PTYPE allows you to
type in one
or more atoms explicitly, while using built in
data for some
other atoms.
The second example reads all SBKJC potentials explicitly:
$ECP
C-ECP GEN 2
1
1
----- CARBON U(P) -----
-0.89371
1 8.56468
2
----- CARBON U(S)-U(P) -----
1.92926 0 2.81497
14.88199
2 8.11296
H-ECP NONE
O-ECP GEN 2
1
1
----- OXYGEN U(P) -----
-0.92550
1 16.11718
2
----- OXYGEN U(S)-U(P) -----
1.96069 0 5.05348
29.13442
2 15.95333
O-ECP
H-ECP NONE
$END
Again, the 2nd
oxygen copies from the first. It is handy
to use the rest
of card -2- as a descriptive comment.
As a final example,
for antimony we have LMAX+1=3 (there
are core d's).
One must first enter U(f), followed by
U(s)-U(f), U(p)-U(f),
U(d)-U(f).
Caution:
-------
At the present
time, there are some numerical problems in
the ECP code,
which has been programed to use spdfg basis
sets, and core
potentials up to g. If one is using a g
basis function,
or a g potential (bottom row elements
beyond the lanthanide
series), there are small errors in
the ECP integrals.
A tight optimization (OPTTOL=1D-05)
will usually
result in the energy rising slightly during
the last few
geometry steps. The error seems to be about
0.000002 Hartree,
so be cautious about using a g potential
or basis function.
When both are used in the same run,
the error in
the energy is about 0.0002, which means the
use of both
is too inaccurate to be trusted. The use of
f functions
or f potentials (or lower) seems to be free of
any errors.
==========================================================
1
$RELWFN
==========================================================
$RELWFN group (optional)
This group is relevant if RELWFN in $CONTRL chose the
NESC or RESC
option for elimination of small components
from relativistic
wavefunctions, to produce a corrected
single component
wavefunction. In case of RESC, only the
one electron
integral corrections are added, whereas for
NESC, corrections
to two electron integrals are accounted
for by means
of a relativistically averaged basis set.
Analytic gradients are programmed for both RESC and
NESC computations.
For NESC, the one electron part of
the spin-orbit
operator can be corrected, while for RESC,
one can compute
spin-orbit coupling with relativistic
corrections
to both one and two electron SOC integrals,
unless internal
uncontraction is requested. In this case
only one electron
SOC integrals are modified. It should
be noted that
internally uncontracted basis functions with
very large exponents
have large SOC integrals, thus the
average asymmetry
due to RESC appears to be larger (before
contraction).
For NESC, you must provide three basis sets, for the
large and small
components and an averaged one, which are
given in $DATAL,
$DATAS, $DATA, respectively. The only
possible choice
for these basis sets is due to Dyall, and
these are available
from
http://www.emsl.pnl.gov:2080/forms/basisform.html
Their names
are similar to cc-pVnZ(pt/sf/lc), pt=point or
fi=finite nucleus,
sf for spin-free and the final field is
lc=large component
($DATAL), sc=small component ($DATAS),
and wf is a
typo for Foldy-Wouthuysen 2e- basis ($DATA).
In GAMESS you
can only use point nucleus approximation.
The need to
input three basis sets means that you cannot
use a $BASIS
group, and you must use COORD=UNIQUE style
input in the
various $DATA's. The three $DATA groups must
contain identical
information except for the primitive
expansion coefficients,
as the three basis sets must have
the same exponents.
In case the option to treat only some
atoms relativistically
is chosen, all non-relativistic
atoms must have
identical basis input in all three groups.
For RESC, ordinary basis sets are used. This however
is a misleading
statement, for while any basis set will
run, accurate
answers may be hard to obtain without the
use of basis
sets contracted using the RESC approximation.
Experience is
showing that large uncontracted basis sets
using non-relativistic
exponents are probably OK, but that
standard contractions
of these in NR calculations can lead
to spurious
results. Unfortunately, contractions using
the RESC approximation
are not yet available for ordinary
use.
1
$RELWFN
OPRESC gives
additive (bitwise) options, which pertain to
the RESC method:
= 0 original RESC implementation, reproduces the
results prior to June 2001. The accuracy of
the RI may be inadequate. (default)
= 1 to obtain more accurate integrals, use the
Gaussian primitives rather than the contracted
basis set to define the resolution of the
identity (RI), used to simplify the integrals in
order to evaluate them in closed form. This
internally uncontracted basis set can be large,
but produces considerably increased accuracy in
the integrals (see also NRATOM/CHARGE).
= 2 HONDO's implementation of the RI for RESC is
mimicked, namely that for ISPHER=+1 the space
used for the RI will not have spherical
contaminants (similarly to MO space).
No gradients for HONDO style are available.
= 4 split L-shells into s and p when generating the
internally uncontracted basis set. This is
necessary if you are using s or p primitives
with the same exponents declared as some L
shell. In such a case, the L shell must be
entered before the s or p. 4 requires 1.
These options are additive, for example OPRESC=5
is needed to select 1 as well as 4.
NESOC =
relativistic corrections to SOC integrals.
Choose only if RELWFN=RESC or NESC, and if
OPERAT=HSO1, HSO2P, or HSO2, for RUNTYP=TRANSITN
= 0 no corrections
= 1 one-electron spin-orbit integrals (NESC default)
= 2 one and two-electron integrals (RESC default)
For RESC and OPRESC=1, NESOC=2 is not implemented,
use NESOC=1 as the closest available possibility.
NRATOM the number
of different elements to be treated
nonrelativistically. For example, in Pb3O4, to
treat only lead relativistically, enter NRATOM=1.
For NESC, this parameter affects the choice of the
basis sets, you should use identical large, small,
and averaged basis set for such atoms.
For RESC, this parameter means that OPRESC=1 will
not cause uncontracting primitives for such atoms.
(default=0)
CHARGE array
containing charges of atoms to be treated
nonrelativistically. (e.g. CHARGE(1)=8.0, to drop
all oxygen atoms)
1
$RELWFN
* * * the next parameters are used only with RELWFN=RESC:
QMTTOL same as
in $CONTRL, but used for the preparation of
the RI space for RESC. (default: from $CONTRL).
RESCTO tolerance
for equating nearly degenerate eigenvalues
of the kinetic energy and overlaps, which is used
for evaluating RESC gradient. Values that are too
large (>1e-6) can cause numerical errors in the
gradient, approximately on the same order as RESCTO.
Too small values can add very large values to the
gradient due to division by numbers that are zero
within machine precision that are not avoided with
this tolerance filter. The recommended values for
OPRESC=1 are 1e-6 for gold to 1e-7 for silver. For
OPRESC=0, 1d-8 or smaller can be used.
==========================================================
1
$EFIELD
==========================================================
$EFIELD group (not required)
This group permits the study of the influence of an
external electric
field on the molecule. The method is
general, and
so works for all ab initio SCFTYPs.
EVEC
= an array of the three x,y,z components of
the applied electric field.
SYM
= a flag to specify when the field to be
applied breaks the molecular symmetry.
Since most fields break symmetry, the
default is .FALSE.
==========================================================
Restrictions:
analytic hessians are not available, but
numerical hessians
are. Because an external field causes
a molecule with
a dipole to experience a torque, geometry
optimizations
must be done in Cartesian coordinates only.
Internal coordinates
eliminate the rotational degrees of
freedom, which
are no longer free.
Notes: a hessian
calculation will have two rotational
modes with non-zero
"frequency", caused by the torque.
A gas phase
molecule will rotate so that the dipole
moment is anti-parallel
to the applied field. To carry
out this rotation
during geometry optimization will take
many steps,
and you can help save much time by inputting
a field opposite
the molecular dipole. There is also
a stationary
point at higher energy with the dipole
parallel to
the field, which will have two imaginary
frequencies
in the hessian. Careful, these will appear
as the first
two modes in a hessian run, but will not
have the i for
imaginary included on the printout since
they are rotational
modes.
For an application,
see
H.Kono, S.Koseki, M.Shiota, Y.Fujimura
J.Phys.Chem.A 105, 5627-5636(2001)
1
$INTGRL
==========================================================
$INTGRL group (optional)
This group controls AO integral formats. It should
probably never
be given, as the program always picks
sensible values.
SCHWRZ
= a flag to activate use of the Schwarz inequality
to predetermine small integrals. There is no
loss of accuracy when choosing this option, and
there are appreciable time savings for bigger
molecules. Default=.TRUE. for over 5 atoms, or
for direct SCF, and is .FALSE. otherwise.
QFMM
= a flag to use the quantum fast multipole method
for linear scaling Fock matrix builds. This is
available for RHF, UHF, and ROHF wavefunctions,
but should not be used with any correlation
treatment. You must select DIRSCF=.TRUE. in
$SCF if you use this option.
The Optimal Parameter FMM code will run at a
comparable speed to a ordinary run doing all
integrals for molecules about 15 Angstroms in
size, and should run faster for 20 Angtroms or
more. See also the $FMM group. (default=.FALSE.)
NINTMX
= Maximum no. of integrals in a record block.
(default=15000 for J or P file, =10000 for PK)
Various antiquated parameters follow:
NOPK
= 0 PK integral option on, which is permissible
for RHF, UHF, ROHF, GVB energy/gradient runs.
= 1 PK option off (default for all jobs).
Must be off for anything with a transformation.
NORDER
= 0 (default)
= 1 Sort integrals into canonical order. There
is little point in selecting this option, as
no part of GAMESS requires ordered integrals.
See also NSQUAR.
The following parameters control the integral sort.
NSQUAR
= 0 Sorted integrals will be in triangular
canonical order (default)
= 1 instead sort to square canonical order.
NDAR
= Number of direct access logical records to be
used for the integral sort (default=2000)
LDAR
= Length of direct access records (site dependent)
NBOXMX
= 200 Maximum number of bins.
NWORD
= 0 Memory to be used (default=all of it).
NOMEM
= 0 If non-zero, force external sort.
The following parameters control integral restarts.
IST=
1 JST= 1
KST= 1 LST= 1
NREC=
1 INTLOC= 1 (values given are
defaults)
==========================================================
1
$FMM
==========================================================
$FMM group (relevant if QFMM selected in $INTGRL)
This group controls the quantum fast multipole method
evaluation of
Fock matrices. The defaults are reasonable,
so there is
little need to give this input.
ITGERR = Target
error in final energy, to 10**-(ITGERR)
Hartree. The accuracy is usually better than
the setting of ITGERR, in fact QFMM runs should
suffer no loss of accuracy or be more accurate
than a conventional integral run (default=7).
QOPS
= a flag to use the Quantum Optimum Parameter
Searching technique, which finds an optimum FMM
parameter set. (Default=.TRUE.)
If QOPS=.FALSE.,
the ITGERR value is not used. In this
case the user
should specify the following parameters:
NP
= the highest multipole order for FMM (Default=15).
NS
= the highest subdivision level (Default=2).
IWS
= the minimum well-separateness (Default=2).
IDPGD
= point charge approximation error (10**(-IDPGD))
of the Gaussian products (Default=9).
IEPS
= very fast multipole method (vFMM) error,
(10**(-IEPS)) (Default=9)
==========================================================
1
$TRANS
==========================================================
$TRANS group
(optional for -CI- or -MCSCF-)
(relevant to analytic hessians)
(relevant to energy localization)
This group controls the integral tranformation. MP2
integral transformations
are controlled instead by the
$MP2 input group.
There is little reason to give any but
the first variable.
DIRTRF
= a flag to recompute AO integrals rather than
storing them on disk. The default is .FALSE.
for MCSCF and CI runs. If your job reads $SCF,
and you select DIRSCF=.TRUE. in that group, a
direct transformation will be done, no matter
how DIRTRF is set.
Note that the transformation may do many passes over
the AO integrals for large basis sets, and thus the
direct recomputation of AO integrals can be very time
consuming.
MPTRAN
= method to use for the integral transformation.
the default is try 0, then 1, then 2.
0 means use the incore method
1 means use the segmented method. This is the
only method that works in parallel.
2 means use the alternate method, which uses
less memory than 2, but requires an extra
large disk file.
NWORD
= Number of words of fast memory to allow. Zero
uses all available memory. (default=0)
CUTTRF
= Threshold cutoff for keeping transformed two
electron integrals. (default= 10**(-9))
AOINTS
= defines AO integral storage during conventional
integral transformations, during parallel runs.
DUP stores duplicated AO lists on each node, and
is the default for parallel computers with slow
interprocessor communication, e.g. ethernet.
DIST distributes the AO integral file across
all nodes, and it is the default for parallel
computers with high speed communications.
==========================================================
1
The remaining groups apply only to -CI- and -MCSCF- runs.
* * * * * * * * * * * * * * * * * * * *
For hints on how to do -CI- and -MCSCF-
see the 'further information' section
* * * * * * * * * * * * * * * * * * * *
$CIINP
==========================================================
$CIINP group (optional, relevant for CITYP=GUGA or ALDET)
This group is the control box for Graphical Unitary
Group Approach
(GUGA) CI calculations, or Ames Laboratory
determinant
(ALDET) full CI. Each step which is executed
potentially
requires a further input group described later.
NRNFG = An array
of 10 switches controlling which steps of
a CI computation are performed.
1 means execute the module, 0 means don't.
NRNFG(1)
= Generate the configurations. See either
$CIDRT or $CIDET input. (default=1)
NRNFG(2) = Transform the integrals. See $TRANS.
(default=1)
NRNFG(3) = Sort integrals and calculate the Hamiltonian
matrix. See $CISORT and $GUGEM. (default=1)
This does not apply to ALDET.
NRNFG(4) = Diagonalize the Hamiltonian matrix.
See $GUGDIA or $CIDET. (default=1)
NRNFG(5) = Construct the one electron density matrix,
and generate NO's. See $GUGDM or $CIDET.
(default=1)
NRNFG(6) = Construct the two electron density matrix.
See $GUGDM2 or $CIDET.
(default=0 normally, but 1 for CI gradients)
NRNFG(7) = Construct the Lagrangian of the CI function.
Requires DM2 matrix exists. See $LAGRAN.
(default=0 normally, but 1 for CI gradients)
This does not apply to ALDET.
NRNFG(8-10) are not used.
Users are not
encouraged to change these values, as the
defaults are
quite reasonable ones.
NPFLG = An array
of 10 switches to produce debug printout.
There is a one to one correspondance to NRNFG, set
to 1 for output. (default = 0,0,0,0,0,0,0,0,0,0)
The most interesting is NPFLG(2)=1 to see the
transformed 1e- integrals, NPFLG(2)=2 adds the
very numerous transformed 2e- integrals to this.
IREST = n
Restart the -CI- at stage NRNFG(n).
==========================================================
1
$DET/$CIDET/$GEN/$CIGEN
==========================================================
$DET group
(required for SCFTYP=MCSCF if CISTEP=ALDET)
$GEN group
(required for SCFTYP=MCSCF if CISTEP=GENCI)
$CIDET group
(required if CITYP=ALDET)
$CIGEN group
(required if CITYP=GENCI)
This group describes the determinants to be used in a
MCSCF or CI
wavefunction. For $DET or $CIDET, a full list
of determinants
is generated, i.e. the MCSCF is of FORS
(also known
as CAS) type, or this is to be a full CI.
For $GEN and
$CIGEN you must input an arbitrary set of
determinants,
according to the keyword GLIST.
Determinants contain several spin states, in contrast
to configuration
state functions. The Sz quantum number
of each determinant
is the same, but the Hamiltonian
eigenvectors
will have various spins S=Sz, Sz+1, Sz+2, ...
so NSTATE may
need to account for states of other spin
symmetry.
In Abelian groups, you can specify the exact
spatial symmetry
you desire.
GLIST =
general determinant list option
The keyword GLIST must not be given in a $DET or
$CIDET input group! These both generate full
determinant lists, automatically.
= INPUT means an input $GCILST group will be read.
= EXTRNL means the list will be read from a disk
file GCILIST generated in an earlier run.
= SACAS requests generation of sevaral CAS spaces
of different space symmetries, specified by
the input IRREPS. This option is intended
for state averaged calculations for cases
of high symmetry, where degenerate irreps
of the true group may fall into different
irreps of the Abelian subgroup used.
The next three determine the symmetry of the states:
GROUP =
name of the point group. The default is to copy
this from $DATA, if that group is Abelian (C2,
Ci, Cs, C2v, C2h, D2, or D2h). If not, the
group is set to C1 (no symmetry used).
ISTSYM = specifies
the spatial symmetry of the state.
This table is exactly the same as in $DRT input.
ISTSYM= 1 2 3 4 5
6 7 8
C1 A
Ci Ag Au
Cs A' A''
C2 A B
C2v A1 A2 B1 B2
C2h Ag Bu Bg Au
D2 A B1 B2 B3
D2h Ag B1g B2g B3g Au B1u B2u B3u
Default is ISTSYM=1, the totally symmetric state.
1
$DET/$CIDET
IRREPS = specifies
the symmetries of the GLIST=SACAS space
determinant list. This variable should always be
an array, as a single symmetry is more quickly
obtained by the regular full CI code. The values
given have the same meaning as the ISTSYM table.
The next four define the filled and active orbital space.
There is no default for NCORE, NACT, and NELS:
NCORE =
total number of orbitals doubly occupied in all
determinants.
NACT = total number of active orbitals.
NELS = total number of active electrons.
SZ
= azimuthal spin quantum number for each of the
determinants, two times SZ is therefore the
number of excess alpha spins in each determinant.
The default is SZ=S, extracted from the MULT=2S+1
given in $CONTRL.
* * * the following control the diagonalization * * *
NSTATE = Number
of CI states to be found, the default is
1. The maximum number of states is 100.
PRTTOL = Printout
tolerance for CI coefficients, the
default is to print any larger than 0.05.
ITERMX = Maximum
number of Davidson iterations per root.
The default is 100. A CI calculation will fail
if convergence is not obtained before reaching
the limit. MCSCF computations will not bomb
if the iteration limit is reached, instead the
last CI vector is used to proceed into the next
orbital update. In cases with very large active
spaces, it may be faster to input ITERMX=2 or 3
to allow the program to avoid fully converging
the CI eigenvalue problem during the early MCSCF
iterations. For small active spaces, it is
best to allow the CI step to be fully converged
on every iteration.
CVGTOL = Convergence
criterion for Davidson eigenvector
routine. This value is proportional to the
accuracy of the coeficients of the eigenvectors
found. The energy accuracy is proportional to
its square. The default is 1.0E-5.
NHGSS =
dimension of the Hamiltonian submatrix which
is diagonalized to obtain the initial guess
eigenvectors. The determinants forming the
submatrix are chosen on the basis of a low
diagonal energy, or if needed to complete a
spin eigenfunction. The default is 300.
1
$DET/$CIDET
NSTGSS = Number
of eigenvectors from the initial guess
Hamiltonian to be included in the Davidson's
iterative scheme. It is seldom necessary to
include extra states to obtain convergence to
the desired states. The default equals NSTATE.
MXXPAN = Maximum
number of expansion basis vectors in the
iterative subspace during the Davidson iterations
before the expansion basis is truncated. The
default is the larger of 10 or 2*NSTGSS. Larger
values might help convergence, do not decrease
this parameter below 2*NSTGSS.
* * * the following
control the 1st order density * * *
These are ignored
during MCSCF, but are used during a CI.
IROOT =
the root whose density is saved on the disk file
for subsequent property analysis. Only one root
can be saved, and the default value of 1 means
the ground state. Be sure to set NFLGDM to form
the density of the state you are interested in!
NFLGDM = Controls
each state's density formation.
0 -> do not form density for this state.
1 -> form density and natural orbitals for this
state, print and punch occ.nums. and NOs.
2 -> same as 1, plus print density over MOs.
The default is NFLGDM(1)=1,0,0,...,0 meaning
only ground state NOs are generated.
* * * the following control the state averaged
* * * 1st and 2nd order density matrix computation
Usually ignored
by CI runs, these are relevant to MCSCF.
PURES =
a flag controlling the spin purity of the state
avaraging. If true, the WSTATE array pertains
to the lowest states of the same S value as is
given by the MULT keyword in $CONTRL. In this
case the value of NSTATE will need to be bigger
than the total number of weights given by WSTATE
if there are other spin states present at low
energies. If false, it is possible to state
average over more than one S value, which might
be of interest in spin-orbit coupling jobs.
The default is .TRUE.
WSTATE = An array
of up to 100 weights to be given to the
densities of each state in forming the average.
The default is to optimize a pure ground state,
WSTATE(1)=1.0,0.0,...,0.0
A small amount of the ground state can help the
convergence of excited states greatly.
Gradient runs are possible only with pure states.
Be sure to set NSTATE above appropriately!
==========================================================
1
$GCILST
==========================================================
$GCILST group
(required for SCFTYP=MCSCF if CISTEP=GENCI)
(required if CITYP=GENCI)
This group defines space products to be used in the
general CI calculation,
or in a MCSCF wavefunction. The
input is free
format.
Line 1: NSPACE ISYM
The first line
gives the total number of space products to
be entered in
the second lines. The option ISYM can be
omitted, or
given as 0, in which case the program will
verify that
all space products typed in the second lines
indeed have
the spatial symmetry defined by ISTSYM in the
$GEN or $CIGEN
input groups. If ISYM is 1, the user is
indicating that
more than one space symmetry is known to
be in the list,
that this is intentional, and the program
should proceed
with the calculation. This might be of use
in state averaging
two representations in a group that has
more than two
total representations, and therefore faster
than turning
symmetry off completely by GROUP=C1. ISYM=2
has the same
meaning but turns on additional printing.
Line 2 is repeated
NSPACE times. Each line 2 contains NACT
integers, which
must be 0, 1, or 2, and therefore tells the
occupation of
each of the active orbitals in each space
product.
An example input is:
$GEN
GLIST=INPUT NELS=6 NACT=4 SZ=0.0 $END
$GCILST
5
2 2 2 0
2 1 2 1
2 0 2 2
2 2 0 2
0 2 2 2
$END
which generates
6 Ms=0 determinants, much less than the 16
determinants
in a C1 symmetry full list for 6 e- in 4 MOs.
The second space
product above generates two determinants.
All space products
with singly occupied orbitals are used
to form all
possible determinants, to ensure that the final
states are eigenfunctions
of the S**2 operator (meaning
they will be
pure spin states).
Note that there
is no way at present to generate lists
such as singles
and doubles from a single reference.
Convergence of
MCSCF calculations will depend on how well
chosen your
determinant list is, and may very well require
the use of the
JACOBI or FULLNR convergers.
==========================================================
1
$DRT/$CIDRT
==========================================================
$DRT group
(required for SCFTYP=MCSCF if CISTEP=GUGA)
$CIDRT group
(required if CITYP=GUGA)
This group describes the -MCSCF- or -CI- wavefunction.
The distinct
row table is the means by which the Graphical
Unitary Group
Approach (GUGA) names the configurations.
The group is spelled DRT for MCSCF runs, and CIDRT for
CI runs.
The main difference in these is NMCC vs. NFZC.
There is no default for GROUP, and you must choose one
of FORS, FOCI,
SOCI, or IEXCIT.
GROUP = the name
of the point group to be used. This is
usually the same as that in $DATA, except for
RUNTYP=HESSIAN, when it must be C1. Choose from
the following: C1, C2, CI, CS, C2V, C2H, D2, D2H,
C4V, D4, D4H. If your $DATA group is not listed,
choose only C1 here.
FORS =
flag specifying the Full Optimized Reaction Space
set of configuration should be generated. This
is usually set true for MCSCF runs, but if it is
not, see FORS in $MCSCF. (Default=.FALSE.)
FOCI =
flag specifying first order CI. In addition to
the FORS configurations, all singly excited CSFs
from the FORS reference are included.
Default=.FALSE.
SOCI =
flag specifying second order CI. In addition to
the FORS configurations, all singly and doubly
excited configurations from the FORS reference
are included. (Default=.FALSE.)
IEXCIT= electron
excitation level, for example 2 will
lead to a singles and doubles CI. This variable
is computed by the program if FORS, FOCI, or
SOCI is chosen, otherwise it must be entered.
1
$DRT/$CIDRT
* * the next variables define the single reference * *
The single configuration reference is defined by
filling in the
orbitals by each type, in the order shown.
The default
for each type is 0.
Core orbitals, which are always doubly occupied:
NMCC = number
of MCSCF core MOs (in $DRT only).
NFZC = number
of CI frozen core MOs (in $CIDRT only).
Internal orbitals, which are partially occupied:
NDOC = number
of doubly occupied MOs in the reference.
NAOS = number
of alpha occupied MOs in the reference,
which are singlet coupled with a corresponding
number of NBOS orbitals.
NBOS = number
of beta spin singly occupied MOs.
NALP = number
of alpha spin singly occupied MOs in the
reference, which are coupled high spin.
NVAL = number
of empty MOs in the reference.
External orbitals, occupied only in FOCI or SOCI:
NEXT = number
of external MOs. If given as -1, this will
be set to all remaining orbitals (apart from any
frozen virtual orbitals).
NFZV = number
of frozen virtual MOs, never occupied.
* * the next two help with state symmetry * *
ISTSYM= irreducible
representation for GUGA wavefunction.
This option overwrites whatever symmetry is implied
by NALP/NAOS/NBOS. Default=0 means ISTSYM will be
inferred from the symmetry of the reference, namely
from the symmetry of NALP/NAOS/NBOS orbitals.
ISTSYM= 1 2 3 4 5
6 7 8
C1 A
Ci Ag Au
Cs A' A''
C2 A B
C2v A1 A2 B1 B2
C2h Ag Bu Bg Au
D2 A B1 B2 B3
D2h Ag B1g B2g B3g Au B1u B2u B3u
It is no doubt easier to just select the desired
ISTSYM directly. Its computation from the singly
occupied orbitals is kept merely to preserve old
input files.
1
NOIRR= controls
labelling of the CI state symmetries.
= 1 no labelling (default)
= 0 usual labelling. This can be very time consuming
if the group is non-Abelian.
=-1 fast labelling, in which all CSFs with small CI
coefficients are ignored. This can produce weights
quite different from one, due to ignoring the small
coefficients, but overall seems to work OK.
Note that it is normal for the weights not to sum
to 1 even for NOIRR=0 because for simplicity the
weight determination is focused on the relative
weights rather than absolute. However weight do
not sum to one only for row-mixed MOs.
= -2,-3... fast labelling and sets SYMTOL=10**NOIRR
for runs other than TRANSITN. All irreps with
weights greater than SYMTOL are considered.
* * * the final choices are seldom used * * *
INTACT= flag
to select the interacting space option.
The CI will include only those spin couplings
which have a nonvanishing matrix element with
the reference configuration.
MXNINT = Buffer size for sorted integrals. (default=20000)
MXNEME = Buffer size for energy matrix. (default=10000)
NPRT
= Configuration printout control switch.
This can consume a HUMUNGUS amount of paper!
0 = no print (default)
1 = print electron occupancies, one per line.
2 = print determinants in each CSF.
3 = print determinants in each CSF (for Ms=S-1).
==========================================================
1
$MCSCF
==========================================================
$MCSCF group (optional for -MCSCF-)
This group controls the MCSCF orbital optimization
step.
The difference between the five convergence methods
is outlined
in Chapter Four of this manual, which you must
carefully study
before attempting MCSCF computations.
--- the next chooses the configuration basis ---
CISTEP = ALDET
chooses the Ames Lab. determinant full CI,
and requires $DET input. (default)
= GUGA chooses the graphical unitary group CSFs,
and requires $DRT input. This is the
only value usable with the QUAD converger.
= GENCI chooses the Ames Lab. general CI, and
requires $GEN input.
--- the next five choose the orbital optimizer ---
FOCAS =
a flag to select a method with a first order
convergence rate. (default=.FALSE.)
SOSCF =
a flag selecting an approximately second order
convergence method, using an approximate orbital
hessian. (default=.TRUE.)
FULLNR = a flag
selecting a second order method, with an
exact orbital hessian. (default=.FALSE.)
QUAD
= a flag to pick a fully quadratic (orbital and
CI coefficient) optimization method, which is
applicable to FORS or non-FORS wavefunctions.
QUAD may not be used with state-averaging.
(default = .FALSE.)
JACOBI = a flag
to pick a program that minimizes the
MCSCF energy by a sequence of 2x2 Jacobi
orbital rotations. This is very systematic in
forcing convergence, although the number of
iterations may be high and the time longer
than the other procedures. This option does
not compute the orbital Lagrangian, hence at
present nuclear gradients may not be computed.
(default = .FALSE.)
Note that FOCAS
must be used only with FORS=.TRUE. in $DRT.
The other convergers
are usable for either FORS or non-FORS
wavefunctions,
although convergence is always harder in the
latter case,
when FORS below must be set .FALSE.
1
$MCSCF
--- the next apply to all convergence methods ---
FORS
= a flag to specify that the MCSCF function is of
the Full Optimized Reaction Space type, which is
sometimes known as CAS-SCF. .TRUE. means omit
act-act rotations from the optimization. Since
convergence is usually better for FULLNR with
these rotations included, the default is sensible
for the case FORS=.TRUE. in $DRT. (default is
.TRUE. for FOCAS/SOSCF, .FALSE. for FULLNR/QUAD)
ACURCY = the
major convergence criterion, the maximum
permissible asymmetry in the Lagrangian matrix.
(default=1.0E-05)
ENGTOL = a secondary
convergence criterion, the run is
considered converged when the energy change is
smaller than this value. (default=1.0E-10)
MAXIT =
Maximum number of iterations (default=100 for
FOCAS, 60 for SOSCF, 30 for FULLNR or QUAD)
MICIT =
Maximum number of microiterations within a
single MCSCF iteration. (default=5 for FOCAS
or SOSCF, or 1 for FULLNR or QUAD)
NWORD =
The maximum memory to be used, the default is
to use all available memory. (default=0)
CANONC = a flag
to cause formation of the closed shell
Fock operator, and generation of canonical core
orbitals. This will order the MCC core by their
orbital energies. (default=.TRUE.)
EKT
= a flag to cause generation of extended Koopmans'
theorem orbitals and energies. (Default=.FALSE.)
For this option, see R.C.Morrison and G.Liu,
J.Comput.Chem., 13, 1004-1010 (1992). Note that
the process generates non-orthogonal orbitals, as
well as physically unrealistic energies for the
weakly occupied MCSCF orbitals. The method is
meant to produce a good value for the first I.P.
NPUNCH = MCSCF
punch option (analogous to $SCF NPUNCH)
0 do not punch out the final orbitals
1 punch out the occupied orbitals
2 punch out occupied and virtual orbitals
The default is NPUNCH = 2.
NPFLG =
an array of debug print control. This is
analagous to the same variable in $CIINP.
Elements 1,2,3,4,6,8 make sense, the latter
controls debugging the orbital optimization.
1
$MCSCF
--- the next refers to SOSCF optimizations ---
NOFO
= set to 1 to skip use of FOCAS for one iteration
during SOSCF. This is a testing parameter, at
present NOFO defaults to 0 to do one FOCAS iter.
--- the next three refer to FOCAS optimizations ---
CASDII = threshold to start DIIS (default=0.05)
CASHFT = level shift value (default=1.0)
NRMCAS = renormalization
flag, 1 means do Fock matrix
renormalization, 0 skips (default=1)
--- the next applies to the QUAD method ---
(note
that all FULLNR input is also relevant to QUAD)
QUDTHR = threshold
on the orbital rotation parameter,
SQCDF, to switch from the initial FULLNR
iterations to the fully quadratic method.
(default = 0.05)
--- The JACOBI converger accepts FULLNR options ---
--- NORB, NOROT, MOFRZ, and FCORE as input ---
--- all remaining input applies only to FULLNR ---
DAMP
= damping factor, this is adjusted by the program
as necessary. (default=0.0)
METHOD = DM2
selects a density driven construction of the
Newton-Raphson matrices. (default).
= TEI selects 2e- integral driven NR construction.
See the 'further information' section for more
details concerning these methods. TEI is slow!
LINSER = a flag
to activate a method similar to direct
minimization of SCF. The method is used if
the energy rises between iterations. It may in
some circumstances increase the chance of
converging excited states. (default=.FALSE.)
FCORE =
a flag to freeze optimization of the MCC core
orbitals, which is useful in preparation for
RUNTYP=TRANSITN jobs. Setting this flag will
automatically force CANONC false. This option
is incompatible with gradients, so can only be
used with RUNTYP=ENERGY. It is a good idea to
decrease TOLZ and TOLE in $GUESS by two orders
of magnitude to ensure the core orbitals are
unchanged during input. (default=.FALSE.)
1
$MCSCF
--- the last four FULLNR options are seldom used ---
DROPC =
a flag to include MCC core orbitals during the
CI computation. The default is to drop them
during the CI, instead forming Fock operators
which are used to build the correct terms in
the orbital hessian. (default = .TRUE.)
NORB
= the number of orbitals to be included in the
optimization, the default is to optimize with
respect to the entire basis. This option is
incompatible with gradients, so can only be used
with RUNTYP=ENERGY. (default=number of AOs
given in $DATA).
MOFRZ =
an array of orbitals to be frozen out of the
orbital optimization step (default=none frozen).
NOROT =
an array of up to 250 pairs of orbital rotations
to be omitted from the NR optimization process.
The program automatically deletes all core-core
rotations, all act-act rotations if FORS=.T.,
and all core-act and core-virt rotations if
FCORE=.T. Additional rotations are input as
I1,J1,I2,J2... to exclude rotations between
orbital I running from 1 to NORB, and J running
up to the smaller of I or NVAL in $TRANS.
==========================================================
1
$MCQDPT
==========================================================
$MCQDPT group (relevant to SCFTYP=MCSCF if MPLEVL=2)
Controls 2nd order MCQDPT (multiconfiguration quasi-
degenerate perturbation
theory) runs, if requested by
MPLEVL=2 in
$CONTRL. MCQDPT2 is implemented only for
FORS (aka CASSCF)
wavefunctions. The MCQDPT method is a
multistate,
as well as multireference perturbation theory.
The implementation
is a separate program, interfaced to
GAMESS, with
its own procedures for determination of the
canonical MOs,
CSF generation, integral transformation,
CI in the reference
CAS, etc. Therefore some of the input
in this group
repeats data given elsewhere, particularly
for $DET/$DRT.
Analytic gradients are not available. Spin-orbit
coupling may
be treated as a perturbation, included at
the same time
as the energy perturbation. If spin-
orbit calculations
are performed, the input groups for
each multiplicity
are named $MCQD1, $MCQD2, ... rather
than $MCQDPT.
Parallel calculation is implemented.
When applied to only one state, the theory is known as
multi-reference
Moller-Plesset (MRMP), so that the term
MCQDPT is more
appropriate when this program is used in its
multistate form.
Please note that this perturbation theory
is not exactly
the same as CASPT2, and should -NEVER- be
called that.
A more complete discussion may be found in
the 'Further
Information' chapter.
*** MCSCF reference wavefunction ***
NEL
= total number of electrons, including core.
(default from $DATA and ICHARG in $CONTRL)
MULT = spin multiplicity (default from $CONTRL)
NMOACT =
Number of orbitals in FORS active space
(default is the active space in $DET or $DRT)
NMOFZC =
number of frozen core orbitals, NOT correlated
in the perturbation calculation. (default is
number of chemical cores)
NMODOC =
number of orbitals which are doubly occupied in
every MCSCF configuration, that is, not active
orbitals, which are to be included in the
perturbation calculation. (The default is all
valence orbitals between the chemical core and
the active space)
NMOFZV =
number of frozen virtuals, NOT occupied during
the perturbation calculation. The default is
to use all virtuals in the MP2. (default=0)
1
$MCQDPT
ISTSYM =
the state symmetry of the target state(s).
This is given as an integer, note that only
Abelian groups in $DATA are supported:
ISTSYM= 1 2 3 4 5
6 7 8
C1 A
Ci Ag Au
Cs A' A''
C2 A B
C2v A1 A2 B1 B2
C2h Ag Bu Bg Au
D2 A B1 B2 B3
D2h Ag B1g B2g B3g Au B1u B2u B3u
(The default is inherited from $DET or $DRT)
NOSYM =
0 use CSF symmetry (see the ISTSYM keyword).
off diagonal perturbations vanish if states are
of different symmetry, so the most efficient
computation is a separate run for every space
symmetry. (default)
1 turn off CSF state symmetry so that all states
are treated at once. ISTSYM is ignored.
Presently this option does not seem to work!!
-1 Symmetry purify the orbitals. Since $GUESS is
not read by MCQDPT runs, this option can be used
as a substitute for its PURIFY. After cleaning
the orbitals, they are reorthogonalised within
each irrep and within each group (core, double,
active, virtual) separately. Since this occurs
after MCSCF optimization (see INORB), it is
*your* responsibility to verify that the changes
made to the orbitals are small enough that the
CAS energies for the original CASSCF and the
CAS-CI performed during MCQDPT give the same
energies!
*** perturbation specification ***
KSTATE=
state is used (1) or not (0) in the MCQDPT2.
Maximum of 20 elements, including zeros.
For example, if you want the perturbation
correction to the second and the fourth roots,
KSTATE(1)=0,1,0,1
See also WSTATE. (default=1,0,0,0,0,0,0,...)
*** MO input and flow control ***
INORB =
0 optimize the MCSCF wavefunction in this run.
= 1 read the converged orbitals from a $VEC group,
and skip immediately to the MCQDPT computation.
A complete $VEC including virtuals must be given.
(default=0)
1
$MCQDPT
*** Intruder State Removal ***
EDSHFT =
energy denominator shifts. (default=0.0,0.0)
Intruder State Free (ISF, aka ISA) calculations
can be made by changing the energy denominators
around poles (where the denominator is zero).
Each denominator x is replaced by x + EDSHFT/x,
so that far from the poles (when x is large) the
effect of such change is small. EDSHFT is an
array of two values, the first is used in spin-
free MCQDPT, and the second is for spin-orbit
MCQDPT. Both values are used if RUNTYP=TRNSTN,
only the first is used otherwise. A suggested
pair of values is 0.02,0.1, but experimentation
with your system is recommended. Setting these
values to zero is ordinary MCQDPT, and infinite
collapses to the MCSCF reference. Note that the
energy denominators (which are ket-dependent in
MCQDPT) are changed in a different way for each
ket-vector, that is, for each row in MCQDPT
Hamiltonian matrix. In other words, the zeroth
order energies are not "universal", but state
specific. This is strictly speaking some weak
inconsistency in defining zeroth order energies
that are usually chosen "universally".
In order to maintain continuity when studying
a PES, one usually uses the same EDSHFT values
for all points on PES. In order to study the
potential surface for any extended range of
geometries, it is recommended to use ISF, as it
is quite likely that one or more regions of the
PES will be unphysical due to intruder states.
For an example of how intruder states can appear
at some points on the PES, see Figures 1,2,7 of
K.R.Glaesemann, M.S.Gordon, H.Nakano
Phys.Chem.Chem.Phys. 1, 967-975(1999)
For a discussion of intruder state removal from
MCQDPT, see
Y.-K.Choe, H.A.Witek, J.P.Finley, K.Hirao
J.Chem.Phys. 114, 3913-3918(2001)
Note however that the formula given above for
the use of EDSHFT does not match the one used
in section V of this paper!
See also REFWGT.
1
$MCQDPT
REFWGT =
a flag to request decomposition of the second
order energy into internal, semi-internal, and
external contributions, and to obtain the weight
of the MCSCF reference in the 1st order wave
function. This option significantly increases
the run time! When you run in parallel, only
the transformation steps will speed up, as the
PT part of the reference weight calculation has
not been adapted for speedups (default = .FALSE.)
The EDSHFT option does not apply if REFWGT is
used. One purpose of using REFWGT is to try to
understand the nature of the intruder states.
*** Canonical Fock orbitals ***
IFORB =
0 omit this step.
= 1 determine the canonical Fock orbitals. (default)
= 3 canonicalise the Fock orbitals averaged over
all $MCQDx input groups. This option pertains
only to RUNTYP=TRANSITN. It is primarily meant
to include spin-orbit coupling perturbation into
the energy perturbation, but could also be used
in conjunction with OPERAT=DM to calculate only
the second order energy perturbation. IFORB=3
means that WSTATE is used as follows: In each
$MCQDx group, the WSTATE weights are divided by
the total number of states (sum(i) IROOTS(i)),
so the sum over all WSTATE values in all $MCQDx
groups is normalized to sum to 1. Thus there is
no normalisation to 1 within each $MCQDx group.
This option might be used to speed up an atomic
MCQDPT, e.g. if computing the 3-P ground state
of carbon, one would want to average over all
three spatial components of the P term, to be
sure of spatial degeneracy, but then run the
perturbation using symmetry, separately on the
B1g+B2g+B3g subspecies (within D2h) of a P term.
It is very important to give weights, appropriate
for the symmetry, the input requires care.
WSTATE =
weight of each CAS-CI state in computing the
closed shell Fock matrix. You must enter 0.0
whenever the same element in KSTATE is 0.
In most cases setting all WSTATE for states
to be included in the MCQDPT to an equal value
is the best.
(default is WSTATE(1)=1.0,0.0,0.0,...)
1
$MCQDPT
*** Miscellaneous options ***
ISELCT
is an option to select only the important CSFs
for inclusion into the CAS-CI reference states.
Set to 1 to select, or to 0 to avoid selection of
CSFs (default = 0)
All CSFs in a preliminary complete active space
CI whose CI coefficients exceed the square root
of THRWGT are kept in a smaller CI to determine
the zero-th order states. Note that the CSFs
with smaller coefficients, while excluded from
the reference states, are still used during the
perturbation calculation, so most of their energy
contribution is still retained. This can save
appreciable computer time in cases with large
active spaces.
THRWGT =
weight threshold for retaining CSFs in selected
configuration runs. In quantum mechanics, the
weight of a CSF is the square of its CI
coefficient. (default=1d-6)
THRGEN =
threshold for one-, two-, and three-body
density matrix elements in the perturbation
calculation. If you want to obtain energies,
for instance, to 6 figures after point, choose
THRGEN=1.0D-08 or 1.0D-09. (default=1.0D-08)
THRENE =
threshold for the energy convergence in the
Davidson's method CAS-CI. (default=-1.0D+00)
THRCON =
threshold for the vector convergence in the
Davidson's method CAS-CI. (default=1.0D-06)
tip: using THRCON=1E-8
THRGEN=1E-10 can help improve the
number of digits of equal energy in degenerate roots.
MDI
= dimension of small Hamiltonian diagonalized
to prepare the initial guess CI states.
MXBASE =
maximum number of expansion vectors in the
Davidson diagonalization (e.g. MXXPAN).
NSOLUT =
number of states to be solved for in the
Davidson's method, this might need to exceed
the number of states in the perturbation
treatment in order to "capture" the correct
roots.
NSTOP =
maximum number of iterations to permit in
the Davidson's diagonalization.
1
$MCQDPT
LPOUT =
print option, 0 gives normal printout, while
<0 gives debug print (e.g. -1, -5, -10, -100)
In particular, LPOUT=-1 gives more detailed
timing information. (default=0)
The next three parameters refer to parallel execution:
DOORD0 =
a flag to select reordering of AO integrals
which speeds the integral transformations.
This reduces disk writes, but increases disk
reads, so you can try turning it off if your
machine has slow writes. (default=.TRUE.)
PARAIO =
access 2e- integral file on every node, at
the same time. This affects only runs with
DOORD0 true, and it may be useful to turn
this off in the case of SMP nodes sharing
a common disk drive. (default=.TRUE.)
DELSCR =
a flag to delete file 56 containing half-
transformed integrals after it has been
used. This reduces total disk requirements
if this file is big. (default=.FALSE.)
Note that parallel
execution will be more effective
if you use distributed
memory, MEMDDI in $SYSTEM. Use
of AOINTS=DIST
in $TRANS is likely to be helpful in
situations with
relatively poor I/O rates compared to
communication,
e.g. SMP enclosures forced to share a
single scratch
disk system. See PROG.DOC for more
information
on parallel execution.
Finally, there
are additional very specialized input
options, described
in the source code routine MQREAD:
IROT, LENGTH,
MAXCSF, MAXERI, MAXROW, MXTRFR, THRERI,
MAINCS, NSTATE
==========================================================
1
$CISORT $GUGEM $GUGDIA
==========================================================
$CISORT group (optional, relevant for -CI- and -MCSCF-)
This group provides further control over the sorting
of the transformed
integrals.
NDAR
= Number of direct access records.
(default = 2000)
LDAR = Length of direct access record (site dependent)
NBOXMX = Maximum
number of boxes in the sort.
(default = 200)
NWORD =
Number of words of fast memory to use in this
step. A value of 0 results in automatic use of
all available memory. (default = 0)
NOMEM = 0 (set to one to force out of memory algorithm)
==========================================================
$GUGEM group (optional, relevant for -CI- or -MCSCF-)
This group provides further control over the
calculation
of the energy (Hamiltonian) matrix.
CUTOFF = Cutoff
criterion for the energy matrix.
(default=1.0E-8)
NWORD = not used.
==========================================================
$GUGDIA group (optional, relevant for -CI- or -MCSCF-)
This group provides control over the Davidson method
diagonalization
step.
NSTATE = Number
of CI states to be found. (default=1)
You can solve for any number of states, but only
100 can be saved for subsequent sections, such
as state averaging.
PRTTOL = Printout
tolerance for CI coefficients
(default = 0.05)
MXXPAN = Maximum
no. of expansion basis vectors used
before the expansion basis is truncated.
(default=30)
1
$GUGDIA
ITERMX = Maximum number of iterations (default=50)
CVGTOL = Convergence
criterion for Davidson eigenvector
routine. This value is proportional to the
accuracy of the coeficients of the eigenvector(s)
found. The energy accuracy is proportional to
its square. (default = 1.0E-5)
NWORD =
Number of words of fast memory to use in this
step. A value of zero results in the use of all
available memory. (default = 0)
MAXHAM = specifies
dimension of Hamiltonian to try to
store in memory. The default is to use all
remaining memory to store this matrix in memory,
if it fits, to reduce disk I/O to a minimum.
MAXDIA = maximum
dimension of Hamiltonian to send to an
incore diagonalization. If the number of CSFs
is bigger than MAXDIA, an iterative Davidson
procedure is invoked. Default=100
NIMPRV = Maximum
no. of eigenvectors to be improved every
iteration. (default = nstate)
NSELCT = Determines
initial guess to eigenvectors.
= 0 -> Unit vectors corresponding to the NSTATE
lowest diagonal elements and any diagonal
elements within SELTHR of them. (default)
< 0 -> First abs(NSELCT) unit vectors.
> 0 -> use NSELCT unit vectors corresponding to
the NSELCT lowest diagonal elements.
SELTHR = Guess
selection threshold when NSELCT=0.
(default=0.01)
NEXTRA = Number
of extra expansion basis vectors to be
included on the first iteration. NEXTRA is
decremented by one each iteration. This may be
useful in "capturing" vectors for higher states.
(default=5)
KPRINT = Print
flag bit vector used when
NPFLG(4)=1 in the $CIINP group (default=8)
value 1 bit 0 print final eigenvalues
value 2 bit 1 print final tolerances
value 4 bit 2 print eigenvalues and tolerances
at each truncation
value 8 bit 3 print eigenvalues every iteration
value 16 bit 4 print tolerances every iteration
==========================================================
1
$GUGDM
==========================================================
$GUGDM group (optional, relevant for -CI-)
This group provides further control over formation of
the one electron
density matrix. See NSTATE in $GUGDIA.
NFLGDM = Controls
each state's density formation.
0 -> do not form density for this state.
1 -> form density and natural orbitals for this
state, print and punch occ.nums. and NOs.
2 -> same as 1, plus print density over MOs.
(default=1,99*0, means ground state NOs only)
Note that forming the 1-particle density for a
state is negligible against the diagonalization
time for that state.
IROOT =
The -CI- root whose density matrix is saved on
the direct access dictionary file for later
computation of properties. You may save only
one state's density for property evaluation.
(default=1)
WSTATE = An array
of up to 100 weights to be given to the
1 body density of each state in forming the DM1.
It is not physically reasonable to average over
any CI states that are not degenerate, but it
may be useful to use WSTATE to produce a totally
symmetric density when the states are degenerate.
The averaged density will be used for property
computations, as well as to generate natural
orbitals. The default is to use NFLGDM/IROOT,
unless WSTATE information is given, in which case
NFLGDM/IROOT are ignored.
IBLOCK = Density
blocking switch. If nonzero, the off
diagonal block of the density below row IBLOCK
will be set to zero before the (approximate)
natural orbitals are found. One use for this is
to keep the internal and external orbitals in a
FOCI or SOCI calculation from mixing, in which
case IBLOCK is the highest occupied internal
orbital. (default=0)
NWORD =
Number of words of memory to use. Zero means use
all available memory (default=0).
==========================================================
1
$GUGDM2
==========================================================
$GUGDM2 group (optional, relevant for -CI- or -MCSCF-)
This group provides control over formation of the
2-particle density
matrix.
WSTATE = An array
of up to 100 weights to be given to the
2 body density of each state in forming the DM2.
The default is to optimize a pure ground state.
(Default=1.0,99*0.0)
A small amount of the ground state can help the
convergence of excited states greatly.
Gradient runs are possible only with pure states.
Be sure to set NSTATE in $GUGDIA appropriately!
CUTOFF = Cutoff
criterion for the 2nd-order density.
(default = 1.0E-9)
NWORD =
Number of words of fast memory to use in sorting
the DM2. The default uses all available memory.
(default=0).
NOMEM =
0 uses in memory sort, if possible.
= 1 forces out of memory sort.
NDAR = Number of direct access records. (default=4000)
LDAR = Length of direct access record (site dependent)
NBOXMX = Maximum no. of boxes in the sort. (default=200)
==========================================================
1
$LAGRAN $TRFDM2
==========================================================
$LAGRAN group (optional, relevant for -CI- gradient)
This group provides further control over formation of
the CI Lagrangian,
a quantity which is necessary for the
computation
of CI gradients.
NOMEM =
0 form in core, if possible
= 1 forces out of core formation
NWORD
= 0 (0=use all available memory)
NDAR
= 4000
LDAR
= Length of each direct access record
(default is NINTMX from $INTGRL)
==========================================================
$TRFDM2 group (optional, relevant for -CI- gradient)
This group provides further control over the back
transformation
of the 2 body density to the AO basis.
NOMEM =
0 transform and sort in core, if possible
= 1 transform in core, sort out of core, if poss.
= 2 transform out of core, sort out of core
NWORD
= 0 (0=use all available memory)
CUTOFF=
1.0D-9, threshold for saving DM2 values
NDAR
= 2000
LDAR
= Length of each direct access record
(default is system dependent)
NBOXMX=
200
==========================================================
Usually neither
of these two groups is given. Since these
groups are normally
used only for CI gradient runs, we
list here some
of the restrictions on the CI gradients:
a)
SCFTYP=RHF, only
b) no FZV orbitals in $CIDRT, all MOs must be used.
c) the derivative integrals are computed in the 2nd
derivative code, which is limited to spd basis sets.
d) the code does not run in parallel.
e) Use WSTATE in $GUGDM2 to specify the state whose
gradient is to be found. Use IROOT in $GUGDM to
specify the state whose other properties will be
found. These must be the same state!
f) excited states often have different symmetry than the
ground state, so think about GROUP in $CIDRT.
g) the gradient can probably be found for any CI for
which you have sufficient disk to do the CI itself.
Time is probably about 2/3 additional.
1
$TRANST
==========================================================
$TRANST group
(relevant for RUNTYP=TRANSITN)
(only for CITYP=GUGA or MPLEVL=2)
This group controls the evaluation of the radiative
transition moment,
or spin orbit coupling (SOC). An SOC
calculation
can be based on variational CI wavefunctions,
using GUGA CSFs,
or based on 2nd order perturbation theory
using the MCQDPT
multireference perturbation theory.
These are termed
SO-CI and SO-MCQDPT below. The orbitals
are typically
obtained by MCSCF computations, and since
the CI or MCQDPT
wavefunctions are based on those MCSCF
states, the
zero-th order states are referred to below as
the CAS-CI states.
SOC jobs prepare a model Hamiltonian
in the CAS-CI
basis, and diagonalize it to produce spin-
mixed states,
which are linear combinations of the CAS-CI
states.
If scalar relativistic corrections were included
in the underlying
spin-free wavefunctions, it is possible
either to include
or to neglect similar corrections to the
spin-orbit integrals,
see keyword NESOC in $RELWFN.
An input file to perform SO-CI will contain
SCFTYP=NONE CITYP=GUGA MPLEVL=0 RUNTYP=TRANSITN
while a SO-MCQDPT
calculation will have
SCFTYP=NONE CITYP=NONE MPLEVL=2 RUNTYP=TRANSITN
The SOC job
will compute a Hamiltonian matrix as the sum
of spin-free
terms and spin-orbit terms, H = H-sf + H-so.
For SO-CI, the
matrix H-sf is diagonal in the CAS-CI state
basis, with
the LS-coupled CAS-CI energies as the diagonal
elements, and
H-so contains only off-diagonal couplings
between these
LS states,
H-sf = CAS-CI spin-free E
H-so = CAS SOC Hamiltonian (e.g. HSO1, HSO2P, HSO2)
For SO-MCQDPT,
the additional input PARMP defines these
matrices differently.
For PARMP=0, the spin-free term
has diagonal
and off-diagonal MCQDPT perturbations:
H-sf - CAS-CI spin-free E + 2nd order spin-free MCQDPT
H-so - CAS SOC Hamiltonian
For PARMP not
equal to 0, the spin orbit operator is also
included into
the perturbing Hamiltonian of the MCQDPT:
H-sf - CAS-CI spin-free E + 2nd order spin-free MCQDPT
H-so - CAS SOC Hamiltonian + 2nd order SO-MCQDPT
Pure transition moment calculations (OPERAT=DM) are
presently limited
to CI wavefunctions, so please use only
CITYP=GUGA MPLEVL=0.
The transition moments computed by
SO-MCQDPT runs
(see TMOMNT flag) will form the transition
density for
the CAS-CI zeroth order states rather than the
1st order perturbed
wavefunctions.
Please see REFS.DOC for additional information on what
is actually
a fairly complex input file to prepare.
1
$TRANST
OPERAT selects
the type of transition being computed.
= DM calculates radiative transition moment
between states of same spin, using
the dipole moment operator. (default)
= HSO1 one-electron Spin-Orbit Coupling (SOC)
= HSO2P partial two electron and full 1e- SOC,
namely core-active 2e- contributions are
computed, but active-active 2e- terms
are ignored. This generally captures
>90% of the full HSO2 computation, but
with spin-orbit matrix element time
similar to the HSO1 calculation.
= HSO2 one and two-electron SOC, this is the
full Pauli-Breit operator.
= HSO2FF one and two-electron SOC, the form factor
method gives the same result as HSO2, but
is more efficient in the case of small
active spaces, small numbers of CAS-CI
states, and large atomic basis sets.
This final option applies only to SO-CI.
PARMP =
controls inclusion of the SOC terms in SO-MCQDPT,
for OPERAT=HSO1 (default=1) or for HSO2P/HSO2
(default=3) only.
0 - no SOC terms should be included in the MCQDPT
corrections at 2nd order, but they will be
included in the CAS states on which the MCQDPT
(i.e. up to 1st order)
1 - include the 1e- SOC perturbation in MCQDPT
-1 - defined under "3", read on...
3 - full 1-electron and partial 2-electron in the
form of the mean field perturbation (this is
very similar to HSO2P, but in the MCQDPT2
perturbation). Only doubly occupied orbitals
(NMODOC) are used for the core 2e contribution.
if the option is set to -1, then all core
orbitals (NMOFZC+NMODOC) are used. Neither
calculation includes extra diagrams including
filled orbitals, so both are "partial".
PARMP=3 (or
-1) has almost no extra cost compared to
PARMP=1, but
can only be used with OPERAT=HSO2 or HSO2P.
The options
-1 and 3 are not rigorously justified, contrary
to HOS2P for
a SO-CI, as 2e integrals with 2 core indices
appear in the
second order in two ways. There is a mean-
field addition
to 1e integrals, which is included when you
choose PARMP=3
or -1. But, there are separate terms from
additional diagrams
that are not implemented, so that there
is some imbalance
in including the partial 2e correction.
Nevetheless,
it may be better to include such "partial"
partial 2e contributions
than not to. Note that at first
order in the
energy (the CAS-CI states) the N-electron
terms are treated
exactly as specified by OPERAT.
1
$TRANST
It is advisable
to tighten up the convergence criteria in
the $MCQDx groups
since SOC is a fairly small effect, and
the spin-free
energies should be accurately computed, for
example THRCON=1e-8
THRGEN=1e-10.
PARMP has a rather
different meaning for OPERAT=HSO2FF:
It refers to
the difference between ket and bra's Ms,
-1 do matrix elements for ms=-1 only
0 do matrix elements for ms=0 only
1 do matrix elements for ms=1 only
-2 do matrix elements for all ms (0, 1, and -1),
which is the default.
-3 calculates form factors so they can be saved
* * * next defines the orbitals and wavefunctions * * *
NUMCI =
For SO-CI, this parameter tells how many CI
calculations to do, and therefore defines how
many $DRTx groups will be read in.
For SO-MCQDPT, this parameter tells how many
MCQDPT calculations to do, and therefore defines
how many $MCQDx groups will be read in.
(default=1)
IROOTS, IVEX, NSTATE, and ENGYST below will all
have NUMCI values. NUMCI may not exceed 64.
You may wish
to define one $DRTx or $MCQDx group for each
spatial symmetry
representation occuring within each spin
multiplicity,
as the use of symmetry during these separate
calculations
may make the entire job run much faster.
NUMVEC = the
meaning is different depending on the run:
a) spin-orbit CI (SO-CI),
Gives the number of different MO sets. This can
be either 1 or 2, but 2 can be chosen only for
FORS/CASSCF or FCI wavefunctions. (default=1)
If you set NUMVEC=2 and you use symmetry in any
of the $DRTx groups, you may have to use ISTSYM
in the $DRT groups since the order of orbitals
from the corresponding orbital transformation
is unpredictable.
b) spin-orbit perturbation (SO-MCQDPT),
The option to have different MOs for different
states is not implemented, so your job will have
only one $VEC1 group, and IVEX will not normally
be input. The absolute value of NUMVEC should be
be equal to the value of NUMCI above. If NUMVEC
positive, the orbitals in the $VEC1 will be used
exactly as given, whereas if NUMVEC is a negative
number, the orbitals will be canonicalized
according to IFORB in $MCQDx. Using NUMVEC=-NUMCI
and IFORB=3 in all $MCQDx to canonicalize over all
states is recommended.
Note that $GUESS
is not read by this RUNTYP! Orbitals must
be in $VEC1
and possibly $VEC2 input groups.
1
$TRANST
NFZC
= For SO-CI, this is equal to NFZC in each $DRTx
group. When NUMVEC=2, this is also the number of
identical core orbitals in the two vector sets.
For SO-MCQDPT, this should be NMOFZC+NMODOC given
in each of the $MCQDx groups.
The default is the number of AOs given in $DATA,
this is not very reasonable.
NOCC
= the number of occupied orbitals. For SO-CI this
should be NFZC+NDOC+NALP+NAOS+NBOS+NVAL, but
add the external orbitals if the CAS-CI states
are CI-SD or FOCI or SOCI type instead of CAS.
For SO-MCQDPT enter NUMFZC+NUMDOC+NUMACT.
The default is the number of AOs given in $DATA,
which is not usually correct.
Note: IROOTS, NSTATE, ENGYST, IVEX contain NUMCI values.
IROOTS = array
containing the number of CAS-CI states to
be used from each CI or MCQDPT calculation.
The default is 1 for every calculation, which is
probably not a correct choice for OPERAT=DM runs,
but is quite reasonable for the HSO operators.
The total number of states included in the SOC
Hamiltonian is the summation of the NUMCI values
of IROOTS times the multiplicity of each CI or
MCQDPT. See also ETOL.
NSTATE = array
containing the number of CAS-CI states to be
found by diagonalising the spin-free Hamiltonians.
Of these, the first IROOTS(i) states will be used
to find transition moments or SOC. Obviously,
enter NSTATE(i) >= IROOTS(i).
The default for NSTATE(i) is IROOTS(i), but might
be bigger if you are curious about the additional
energies, or to help the Davidson diagonalizer.
NSTATE is ignored by SO-MCQDPT runs, and you must
ensure that your IROOTS input corresponds to the
KSTATE option in $MCQDx.
ETOL
= energy tolerance for CI state elimination.
This applies only to SO-CI and OPERAT=HSO1,2,2P.
After each CI finds NSTATE(i) CI roots for each
$DRTx, the number of states kept in the run is
normally IROOTS(i), but ETOL applies the further
constraint that the states kept be within ETOL of
the lowest energy found for any of the $DRTx.
The default is 100.0 Hartree, so that IROOTS is
the only limitation.
IVEX
= Array of indices of $VECx groups to be used for
each CI calculation. The default for NUMVEC=2 is
IVEX(1)=1,2,1,1,1,1,1..., and of course for
NUMVEC=1, it is IVEX(1)=1,1,1,1,1...
This applies only to CITYP=GUGA jobs.
1
$TRANST
ENGYST = energy
values to replace the CI spin-free energies.
This parameter applies to SO-CI only.
A possible use for this is to use first or second
order CI energies (FOCI or SOCI in $DRT) on the
diagonal of the Hamiltonian (obtained in some
earlier runs) but to use only CAS wavefunctions
to evaluate off diagonal HSO matrix elements. The
CAS-CI runs are still conducted to obtain CI coefs,
needed to evaluate the off diagonal elements.
Enter MXRT*NUMCI values as a square array, by the
usual FORTRAN convention (that is, MXRT roots of
$DRT1, MXRT roots of $DRT2 etc), in hartrees, with
zeros added to fill each column to MXRT values.
MXRT is the maximum value in the IROOTS array.
(the default is the computed CAS-CI energies)
* * * the next pertain only to spin-orbit runs * * *
RSTATE = sets
the zero energy level
format: ndrt*1000+iroot for adiabatic state (CI root)
0000 sets zero energy to the lowest diabatic root
default: 1001 (1st root in $DRT1 or $MCQD1)
ZEFTYP specifies
effective nuclear charges to use.
= TRUE uses true nuclear charge of each atom,
except protons are removed if an ECP basis
is being used (default).
= 3-21G selects values optimized for the 3-21G
basis, but these are probably appropriate
for any all electron basis set. Rare gases,
transition metals, and Z>54 will use the
true nuclear charges.
= SBKJC selects a set obtained for the SBKJC ECP
basis set, specifically. It may not be
sensible to use these for other ECP sets.
Rare gases, lanthanides, and Z>86 will use
the true nuclear charges.
ZEFF
= an array of effective nuclear charges, overriding
the charges chosen in ZEFTYP.
Note that effective nuclear charges can be used for
any HSO type OPERAT, but traditionally these are used
mainly for HSO1 as an empirical correction to the
omission of the 2e- term, or to compensate for missing
core orbitals in ECP runs.
1
$TRANST
JZ
controls the calculation of Jz eigenvalues
= 0 do not perform the calculation
= 1 do the calculation
By default, Jz is set to 1 for molecules that are
recognised as linear (this includes atoms!).
Jz cannot be computed for nonlinear molecules.
The matrix of Jz=Lz+Sz operator is constructed
between spin-mixed states (eigenvalues of Hso).
Setting Jz to 1 can enforce otherwise avoided (by
symmetry) calculations of SOC matrix elements.
JZ applies only to HSO1,2,2P.
TMOMNT = flag
to control computation of the transition
dipole moment between spin-mixed wavefunctions
(that is, betweeen eigenvectors of the Pauli-Breit
Hamiltonian). Applies only to HSO1,2,2P.
(default is .FALSE.)
SKIPDM = flag
to omit(.TRUE.) or include(.FALSE.) dipole
moment matrix elements during spin-orbit coupling.
Usually it takes almost no addition effort to
calculate <R> excluding some cases when the
calculation of forbidden by symmetry spin-orbit
coupling matrix elements <Hso> may have to be
performed since <R> and <Hso> are computed
simultaneously. Applies only to HSO1,2,2P.
Since the lack of a MCQDPT density matrix means
there are no MCQDPT dipole moments at present,
SO-MCQDPT jobs will compute the dipole matrix
elements for the CAS-CI states only. However,
the dipole moments in the spin-mixed states will
be computed with the MCQDPT mixing coefficients.
(default is .TRUE.)
IPRHSO = controls
output style for matrix elements (HSO*)
=-1 do not output individual matrix elements
otherwise these are accumulative:
= 0 term-symbol like kind of labelling:
labels contain full symmetry info (default)
= 1 all states are numbered consequently within each
spin multiplicity (ye olde style)
= 2 output only nonzero (>=1e-4) matrix elements
PRTPRM = flag
to provide detailed information about the
composition of the spin-mixed states in terms of
adiabatic states. This flag also provides similar
information about Jz (if JZ set).
(default is .FALSE.)
1
$TRANST
* * * expert mode HSO control options * * *
MODPAR =
parallel options, which are independent bit
options, 0=off, 1=on. Bit 1 refers only to
HSO2FF, bit 2 to HSO1,2,2P. Enter a decimal
value 0, 1, 2, 3 meaning binary 00, 01, 10, 11.
bit 1
= 0/1 (HSO2FF) uses static/dynamic load balancing in
parallel if available, otherwise use static
load balancing. Dynamic algorithm is usually
faster but may utilize memory less efficiently,
and I/O can slow it down. Also, dynamical
algorithm forces SAVDSK=.F. since its
unique distribution of FFs among nodes implies
no savings from precalculating form factors.
bit 2
= 0/1 (HSO1,2,2P) duplicate/distribute SOC integrals
in parallel. If set, 2e AO integrals and the
four-index transformation are divided over
nodes (distributed), and SOC MO integrals are
then summed over nodes.
The default
is 3, meaning both bits are set on (11)
PHYSRC = flag
to force the size of the physical record to
be equal to the size of the sorting buffers.
This option can have a dramatic effect on the
efficency. Usually, setting PHYSRC=.t. is helpful
if the code complains that low memory enforces
SLOWFF=.TRUE., or you set it yourself. For large
active spaces and large memory (more precisely, if
reclen is larger than the physical record size)
PHYSRC=.TRUE. can slow the code down. Setting
PHYSRC to .true. forces SLOWFF to be .false.
See MODPAR. (default .FALSE.) (only with HSO2FF)
RECLEN = specifies
the size of the record on file 40,
where form factors are stored. This parameter
significantly affects performance.
If not specified, RECLEN have to be guessed,
and the guess will usually be either an
overestimate or underestimate. If the former
you waste disk space, if the latter the program
aborts. Note that RECLEN will be different for
each pair of multiplicities and you must specify
the maximum for all pairs. The meaning of this
number is how many non-zero form factors are
present given four MO indices. You can decrease
RECLEN if you are getting a message "predicted
sorting buffer length is greater than needed..."
Default depends on active space. (only HSO2FF)
SAVDSK = flag
to repeat the form factor calculation twice.
This avoids wasting disk space as the actually
required record size is found during the 1st run.
(default=.FALSE.) (only with HSO2FF)
1
$TRANST
SLOWFF = flag
to choose a slower FF write-out method.
By default .FALSE., but this is turned on if:
1) not enough memory for the fast way is available
2) the maximum usable memory is available, as when
the buffer is as large as the maximum needed,
then the "slow FF" algorythm is faster.
Generally SLOWFF=.true. saves up to 50% or so of
disk space. See PHYSRC. (only with HSO2FF)
ACTION
controls disk file DAFL30 reuse.
= NORMAL calculate the form factors in this run.
= SAVE calculate, and store the form factors on
disk for future runs with the same active
space characteristics.
= READ read the form factors from disk from an
earlier run which used SAVE.
(default=NORMAL) (only with HSO2FF)
Note that currently in order to use ACTION =
SAVE or READ you should specify MS= -1, 0, or 1
* * * some control tolerances * * *
NOSYM= -1 forces
use of symmetry-contaminated orbitals
symmetry analysis, otherwise the same as NOSYM=0
= 0 fully use symmetry
= 1 do not use point group symmetry, but still use
other symmetries (Hermiticity, spin).
= 2 use no symmetry. Also, include all CSFs for
HSO1, 2, 2P.
= 3 force the code to assume the symmetry specified
in $DATA is the same as in all $DRT groups, but
is otherwise identical to NOSYM=-1. This option
saves CPU time and money(memory). Since the $DRT
works by mapping non-Abelian groups into their
highest Abelian subgroup, do not use NOSYM=3 for
non-Abelian groups.
SYMTOL = relative
error for the matrix elements. This
parameter has a great impact upon CPU time, and
the default has been chosen to obtain nearly
full accuracy while still getting good speedups.
(default=1.0E-4)
1
$TRANST
* * * the remaining parameters are not so important * * *
PRTCMO = flag
to control printout of the corresponding
orbitals. (default is .FALSE.)
HSOTOL = HSO
matrix elements are considered zero if they
are smaller than HSOTOL. This parameter is used
only for print-out and statistics.
(default=1.0E-1 cm-1)
TOLZ
= MO coefficient zero tolerance (as for $GUESS).
(default=1.0E-8)
TOLE
= MO coefficient equating tolerance (as for
$GUESS). (default=1.0E-5)
==========================================================