Gamma(qp,KM), where

are matrix elements of the one-body property
operator W in a basis set of molecular spin-orbitals used
in the calculations. The calculation of reduced density
matrices provides the most convenient way of calculating CC
and EOMCC properties of ground and excited states. In
addition, by having reduced density matrices, one can
calculate CC and EOMCC electron densities,
rhoK(x) = Sum_pq Gamma(qp,K) (phi_q(x))* phi_p(x),
where phi_p(x) and phi_q(x) are molecular spin-orbitals and
x represents the electronic (spatial and spin) coordinates.
By diagonalizing Gamma(qp,K), one can determine the natural
occupation numbers and natural orbitals for the CC or EOMCC
state |PsiK>.
The above strategy of handling molecular properties
analytically by determining one-body reduced density
matrices was implemented in the CC/EOMCC programs
incorporated in GAMESS. At this time, the calculations of
reduced density matrices and selected properties are
possible at the CCSD (ground states) and EOMCCSD (ground
and excited states) levels of theory (T=T1+T2, R=R1+R2,
L=L0+L1+L2). Currently, in the main output the program
prints the CCSD and EOMCCSD electric multipole (dipole,
quadrupole, etc.) moments and several other one-electron
properties that one can extract from the CCSD/EOMCCSD
density matrices, the EOMCCSD transition dipole moments and
the corresponding dipole and oscillator strengths, and the
natural occupation numbers characterizing the CCSD/EOMCCSD
wave functions. In addition, the complete CCSD/EOMCCSD one-
body reduced density matrices and transition density
matrices in the RHF molecular orbital basis and the CCSD
and EOMCCSD natural orbital occupation numbers are printed
in the PUNCH output file. The eigenvalues of the density
matrix (natural occupation numbers) are ordered such that
the corresponding eigenvectors (CCSD or EOMCCSD natural
orbitals) have the largest overlaps with the consecutive
ground-state RHF MOs. Thus, the first eigenvalue of the
density matrix corresponds to the CCSD or EOMCCSD natural
orbital that has the largest overlap with the RHF MO 1, the
second with RHF MO 2, etc. This ordering is particularly
useful for analyzing excited states, since in this way one
can easily recognize orbital excitations that define a
given excited state.
One has to keep in mind that the reduced density
matrices originating from CC and EOMCC calculations are not
symmetric. Thus, if we, for example, want to calculate the
dipole strength between states K and M for the x component
of the dipole mu_x, |~~ values for
unrestricted DFT are smaller than for unrestricted HF.
2. GAMESS computes the ~~~~ quantity in an approximate
way, namely it pretend that the Kohn-Shan orbitals can be
used to form a determinant (WRONG, WRONG, WRONG, there is
no wavefunction in DFT!!!) and then uses the same formula
that UHF uses to evaluate that determinant's spin
expectation value. See
G.J.Laming, N.C.Handy, R.D.Amos
Mol.Phys. 80, 1121-1134(1993)
J.Baker, A.Scheiner, J.Andzelm
Chem.Phys.Lett. 216, 380-388(1993)
C.Adamo, V.Barone, A.Fortunelli
J.Chem.Phys. 98, 8648-8652(1994)
J.A.Pople, P.M.W.Gill, N.C.Handy
Int.J.Quantum Chem. 56, 303-305(1995)
J.Wang, A.D.Becke, V.H.Smith
J.Chem.Phys. 102, 3477-3480(1995)
J.M.Wittbrodt, H.B.Schlegel
J.Chem.Phys. 105, 6574-6577(1996)
J.Grafenstein, D.Cremer
Mol.Phys. 99, 981-989(2001)
and commentary in Koch & Holthausen, pp 52-54.
Orbital energies:
The discussion on page 49-50 of Koch and Holthausen shows
that although the highest occupied orbital's eigenvalue
should be the ionization potential for exact Kohn-Sham
calculations, the functionals we actually have greatly
underestimate IP values. The 5th reference below shows how
inclusion of HF exchange helps this, and provides a linear
correction formula for IPs. The first two papers below
connect the HOMO eigenvalue to the IP, and the third shows
that while the band gap is underestimated by existing
functionals, the gap's center is correctly predicted.
However, the 5th paper shows that DFT is actually pretty
hopeless at predicting these gaps. The 4th paper uses SCF
densities to generate exchange-correlation potentials that
actually give fairly good IP values:
J.F.Janak Phys.Rev.B 18, 7165-7168(1978)
M.Levy, J.P.Perdew, V.Sahni
Phys.Rev.A 30, 2745-2748(1984)
J.P.Perdew, M.Levy Phys.Rev.Lett. 51, 1884-1887(1983)
A.Nagy, M.Levy Chem.Phys.Lett. 296, 313-315(1998)
G.Zhang, C.B.Musgrave J.Phys.Chem.A 111, 1554-1561(2007)
Summary of excited state methods
This is not a "how to" section, as the actual calculations
will be carried out by means described elsewhere in this
chapter. Instead, a summary of methods that can treat
excited states is given.
The simplest possibility is SCFTYP. For example, a closed
shell molecule's first triplet state can always be treated
by SCFTYP=ROHF MULT=3. Assuming there is some symmetry
present, the GVB program may be able to do excited singlets
variationally, provided they are of a different space
symmetry than the ground state. The MCSCF program gives a
general entree into excited states, since upper roots of a
Hamiltonian are always variational: see for example NSTATE
and WSTATE and IROOT in $DET. Of course, 2nd order
perturbation theory can include correlation energy into
these SCF level calculations. Note in particular the
usefulness of quasi-degenerate multireference perturbation
theory when electronic states have similar energies.
CI calculations also give a simple entree into excitated
states. There are a variety of programs, selected by CITYP
in $CONTRL. Note in particular CITYP=CIS, programmed for
closed shell ground states, with gradient capability for
singlet excited states, and for the calculation of triplet
state energies. The other CI programs can generate very
flexible wavefunctions for the evaluation of the excitation
energy, and property values. Note that the GUGA program
will do nuclear gradients provided the reference is RHF.
The TD-DFT method treats singly excited states, including
correlation effects, and is a popular alternative to CIS.
The program allows for excitation energies from a UHF
reference, but is much more powerful for RHF references:
nuclear gradients and/or properties may be computed. Use
of a "long range corrected" or "range separated" functional
(the two terms are synonymous) is often thought to be
important when treating charge transfer or Rydberg states:
see the LC=.TRUE. flag or CAMB3LYP. Spin-flip TDDFT allows
the users to select as the reference state something more
appropriate to the orbital optimization stage. See $TDDFT
for details.
Equation of Motion (EOM) coupled cluster can give accurate
estimates of excitation energies. There are no gradients,
and properties exist only for the EOM-CCSD level, but
triples corrections to the energy are available. See
$EOMINP for more details.
Most of the runs will predict oscillator strengths, or
Einstein coefficients, or similar data regarding the
electronic transition moments. Full prediction of UV-vis
spectra is not possible without Franck-Condon information.
Excited states frequently come close together, and
crossings between them are of great interest.
See RUNTYP=TRANSITION for spin-orbit coupling, which is
responsible for InterSystem Crossing (ISC) between states
of different spin multiplicity. See RUNTYP=NACME for the
computation of the non-adiabatic coupling matrix elements
that cause Internal Conversion (IC) between states of the
same spin multiplicity. Alternatively, diabatic potential
surfaces may be generated at the MCSCF or MCQDPT levels:
see DIABAT in the $MCSCF group.
It is possible to search for the lowest energy on the
crossing seam between two surfaces. In case those surfaces
have different spins, or different space symmetries (or
both), see RUNTYP=MEX. When the surfaces have the same
symmetry, see RUNTYP=CONINT for location of conical
intersections.
Solvent effects (EFP and/or PCM) can easily be incorporated
when using SCFTYP to generate the states, and nuclear
gradients are available. It is now possible to assess
solvent effects on TD-DFT excitation energies from closed
shell references, using either EFP or PCM.
Excited states often possess Rydberg character, so diffuse
functions in the basis set are likely to be important.
Geometry Searches and Internal Coordinates
Stationary points are places on the potential energy
surface with a zero gradient vector (first derivative of
the energy with respect to nuclear coordinates). These
include minima (whether relative or global), better known
to chemists as reactants, products, and intermediates; as
well as transition states (also known as saddle points).
The two types of stationary points have a precise
mathematical definition, depending on the curvature of the
potential energy surface at these points. If all of the
eigenvalues of the hessian matrix (second derivative
of the energy with respect to nuclear coordinates) are
positive, the stationary point is a minimum. If there is
one, and only one, negative curvature, the stationary
point is a transition state. Points with more than one
negative curvature do exist, but are not important in
chemistry. Because vibrational frequencies are basically
the square roots of the curvatures, a minimum has all
real frequencies, and a saddle point has one imaginary
vibrational "frequency".
GAMESS locates minima by geometry optimization, as
RUNTYP=OPTIMIZE, and transition states by saddle point
searches, as RUNTYP=SADPOINT. In many ways these are
similar, and in fact nearly identical FORTRAN code is used
for both. The term "geometry search" is used here to
describe features which are common to both procedures.
The input to control both RUNTYPs is found in the $STATPT
group.
As will be noted in the symmetry section below, an
OPTIMIZE run does not always find a minimum, and a
SADPOINT run may not find a transtion state, even though
the gradient is brought to zero. You can prove you have
located a minimum or saddle point only by examining the
local curvatures of the potential energy surface. This
can be done by following the geometry search with a
RUNTYP=HESSIAN job, which should be a matter of routine.
quasi-Newton Searches
Geometry searches are most effectively done by what is
called a quasi-Newton-Raphson procedure. These methods
assume a quadratic potential surface, and require the
exact gradient vector and an approximation to the hessian.
It is the approximate nature of the hessian that makes the
method "quasi". The rate of convergence of the geometry
search depends on how quadratic the real surface is, and
the quality of the hessian. The latter is something you
have control over, and is discussed in the next section.
GAMESS contains different implementations of quasi-
Newton procedures for finding stationary points, namely
METHOD=NR, RFO, QA, and the seldom used SCHLEGEL. They
differ primarily in how the step size and direction are
controlled, and how the Hessian is updated. The CONOPT
method is a way of forcing a geometry away from a minimum
towards a TS. It is not a quasi-Newton method, and is
described at the very end of this section.
The NR method employs a straight Newton-Raphson step.
There is no step size control, the algorithm will simply
try to locate the nearest stationary point, which may be a
minimum, a TS, or any higher order saddle point. NR is
not intended for general use, but is used by GAMESS in
connection with some of the other methods after they have
homed in on a stationary point, and by Gradient Extremal
runs where location of higher order saddle points is
common. NR requires a very good estimate of the geometry
in order to converge on the desired stationary point.
The RFO and QA methods are two different versions of
the so-called augmented Hessian techniques. They both
employ Hessian shift parameter(s) in order to control the
step length and direction.
In the RFO method, the shift parameter is determined by
approximating the PES with a Rational Function, instead of
a second order Taylor expansion. For a RUNTYP=SADPOINT,
the TS direction is treated separately, giving two shift
parameters. This is known as a Partitioned Rational
Function Optimization (P-RFO). The shift parameter(s)
ensure that the augmented Hessian has the correct eigen-
value structure, all positive for a minimum search, and
one negative eigenvalue for a TS search. The (P)-RFO step
can have any length, but if it exceeds DXMAX, the step is
simply scaled down.
In the QA (Quadratic Approximation) method, the shift
parameter is determined by the requirement that the step
size should equal DXMAX. There is only one shift
parameter for both minima and TS searches. Again the
augmented Hessian will have the correct structure. There
is another way of describing the same algorithm, namely as
a minimization on the "image" potential. The latter is
known as TRIM (Trust Radius Image Minimization). The
working equation is identical in these two methods.
When the RFO steplength is close to DXMAX, there is
little difference between the RFO and QA methods. However,
the RFO step may in some cases exceed DXMAX significantly,
and a simple scaling of the step will usually not produce
the best direction. The QA step is the best step on the
hypersphere with radius DXMAX. For this reason QA is the
default algorithm.
Near a stationary point the straight NR algorithm is
the most efficient. The RFO and QA may be viewed as
methods for guiding the search in the "correct" direction
when starting far from the stationary point. Once the
stationary point is approached, the RFO and QA methods
switch to NR, automatically, when the NR steplength drops
below 0.10 or DXMAX, whichever is the smallest.
The QA method works so well that we use it exclusively,
and so the SCHLEGEL method will probably be omitted from
some future version of GAMESS.
You should read the papers mentioned below in order to
understand how these methods are designed to work. The
first 3 papers describe the RFO and TRIM/QA algorithms. A
good but slightly dated summary of search procedures is
given by Bell and Crighton, and see also the review by
Schlegel. Most of the FORTRAN code for geometry searches,
and some of the discussion in this section was written by
Frank Jensen of the University of Aarhus, whose paper
compares many of the algorithms implemented in GAMESS:
1. J.Baker J.Comput.Chem. 7, 385-395(1986)
2. T.Helgaker Chem.Phys.Lett. 182, 503-510(1991)
3. P.Culot, G.Dive, V.H.Nguyen, J.M.Ghuysen
Theoret.Chim.Acta 82, 189-205(1992)
4. H.B.Schlegel J.Comput.Chem. 3, 214-218(1982)
5. S.Bell, J.S.Crighton
J.Chem.Phys. 80, 2464-2475(1984).
6. H.B.Schlegel Advances in Chemical Physics (Ab Initio
Methods in Quantum Chemistry, Part I), volume 67,
K.P.Lawley, Ed. Wiley, New York, 1987, pp 249-286.
7. F.Jensen J.Chem.Phys. 102, 6706-6718(1995).
the nuclear Hessian
Although quasi-Newton methods require only an
approximation to the true hessian, the quality of this
matrix has a great affect on convergence of the geometry
search.
There is a procedure contained within GAMESS for
guessing a positive definite hessian matrix, HESS=GUESS.
If you are using Cartesian coordinates, the guess hessian
is based on pairwise atom stretches. The guess is more
sophisticated when internal coordinates are defined, as
empirical rules will be used to estimate stretching and
bending force constants. Other angular force constants are
set to 1/4. The guess often works well for minima, but
cannot possibly find transition states (because it is
positive definite). Therefore, GUESS may not be selected
for SADPOINT runs.
Two options for providing a more accurate hessian are
HESS=READ and CALC. For the latter, the true hessian is
obtained by direct calculation at the initial geometry, and
then the geometry search begins, all in one run. The READ
option allows you to feed in the hessian in a $HESS group,
as obtained by a RUNTYP=HESSIAN job. The second procedure
is actually preferable, as you get a chance to see the
frequencies. Then, if the local curvatures look good, you
can commit to the geometry search. Be sure to include a
$GRAD group (if the exact gradient is available) in the
HESS=READ job so that GAMESS can take its first step
immediately.
Note also that you can compute the hessian at a lower
basis set and/or wavefunction level, and read it into a
higher level geometry search. In fact, the $HESS group
could be obtained at the semiempirical level. This trick
works because the hessian is 3Nx3N for N atoms, no matter
what atomic basis is used. The gradient from the lower
level is of course worthless, as the geometry search must
work with the exact gradient of the wavefunction and basis
set in current use. Discard the $GRAD group from the lower
level calculation!
You often get what you pay for. HESS=GUESS is free, but
may lead to significantly more steps in the geometry
search. The other two options are more expensive at the
beginning, but may pay back by rapid convergence to the
stationary point.
The hessian update frequently improves the hessian for a
few steps (especially for HESS=GUESS), but then breaks
down. The symptoms are a nice lowering of the energy or
the RMS gradient for maybe 10 steps, followed by crazy
steps. You can help by putting the best coordinates into
$DATA, and resubmitting, to make a fresh determination of
the hessian.
The default hessian update for OPTIMIZE runs is BFGS,
which is likely to remain positive definite. The POWELL
update is the default for SADPOINT runs, since the hessian
can develop a negative curvature as the search progresses.
The POWELL update is also used by the METHOD=NR and CONOPT
since the Hessian may have any number of negative
eigenvalues in these cases. The MSP update is a mixture of
Murtagh-Sargent and Powell, suggested by Josep Bofill,
(J.Comput.Chem., 15, 1-11, 1994). It sometimes works
slightly better than Powell, so you may want to try it.
coordinate choices
Optimization in cartesian coordinates has a reputation
of converging slowly. This is largely due to the fact
that translations and rotations are usually left in the
problem. Numerical problems caused by the small eigen-
values associated with these degrees of freedom are the
source of this poor convergence. The methods in GAMESS
project the hessian matrix to eliminate these degrees of
freedom, which should not cause a problem. Nonetheless,
Cartesian coordinates are in general the most slowly
convergent coordinate system.
The use of internal coordinates (see NZVAR in $CONTRL
as well as $ZMAT) also eliminates the six rotational and
translational degrees of freedom. Also, when internal
coordinates are used, the GUESS hessian is able to use
empirical information about bond stretches and bends.
On the other hand, there are many possible choices for the
internal coordinates, and some of these may lead to much
poorer convergence of the geometry search than others.
Particularly poorly chosen coordinates may not even
correspond to a quadratic surface, thereby ending all hope
that a quasi-Newton method will converge.
Internal coordinates are frequently strongly coupled.
Because of this, Jerry Boatz has called them "infernal
coordinates"! A very common example to illustrate this
might be a bond length in a ring, and the angle on the
opposite side of the ring. Clearly, changing one changes
the other simultaneously. A more mathematical definition
of "coupled" is to say that there is a large off-diagonal
element in the hessian. In this case convergence may be
unsatisfactory, especially with a diagonal GUESS hessian,
where a "good" set of internals is one with a diagonally
dominant hessian. Of course, if you provide an accurately
computed hessian, it will have large off-diagonal values
where those are truly present. Even so, convergence may
be poor if the coordinates are coupled through large 3rd
or higher derivatives. The best coordinates are therefore
those which are the most "quadratic".
One very popular set of internal coordinates is the
usual "model builder" Z-matrix input, where for N atoms,
one uses N-1 bond lengths, N-2 bond angles, and N-3 bond
torsions. The popularity of this choice is based on its
ease of use in specifying the initial molecular geometry.
Typically, however, it is the worst possible choice of
internal coordinates, and in the case of rings, is not
even as good as Cartesian coordinates.
However, GAMESS does not require this particular mix
of the common types. GAMESS' only requirement is that you
use a total of 3N-6 coordinates, chosen from these 3 basic
types, or several more exotic possibilities. (Of course,
we mean 3N-5 throughout for linear molecules). These
additional types of internal coordinates include linear
bends for 3 collinear atoms, out of plane bends, and so on.
There is no reason at all why you should place yourself in
a straightjacket of N-1 bonds, N-2 angles, and N-3
torsions.
If the molecule has symmetry, be sure to use internals
which are symmetrically related.
For example, the most effective choice of coordinates
for the atoms in a four membered ring is to define all
four sides, any one of the internal angles, and a dihedral
defining the ring pucker. For a six membered ring, the
best coordinates seem to be 6 sides, 3 angles, and 3
torsions. The angles should be every other internal
angle, so that the molecule can "breathe" freely. The
torsions should be arranged so that the central bond of
each is placed on alternating bonds of the ring, as if
they were pi bonds in Kekule benzene. For a five membered
ring, we suggest all 5 sides, 2 internal angles, again
alternating every other one, and 2 dihedrals to fill in.
The internal angles of necessity skip two atoms where the
ring closes. Larger rings should generalize on the idea
of using all sides but only alternating angles. If there
are fused rings, start with angles on the fused bond, and
alternate angles as you go around from this position.
Rings and more especially fused rings can be tricky.
For these systems, especially, we suggest the Cadillac of
internal coordinates, the "natural internal coordinates"
of Peter Pulay. For a description of these, see
P.Pulay, G.Fogarosi, 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).
These are linear combinations of local coordinates, except
in the case of rings. The examples given in these two
papers are very thorough.
An illustration of these types of coordinates is given
in the example job EXAM25.INP, distributed with GAMESS.
This is a nonsense molecule, designed to show many kinds
of functional groups. It is defined using standard bond
distances with a classical Z-matrix input, and the angles
in the ring are adjusted so that the starting value of
the unclosed OO bond is also a standard value.
Using Cartesian coordinates is easiest, but takes a very
large number of steps to converge. This however, is better
than using the classical Z-matrix internals given in $DATA,
which is accomplished by setting NZVAR to the correct 3N-6
value. The geometry search changes the OO bond length to
a very short value on the 1st step, and the SCF fails to
converge. (Note that if you have used dummy atoms in the
$DATA input, you cannot simply enter NZVAR to optimize in
internal coordinates, instead you must give a $ZMAT which
involves only real atoms).
The third choice of internal coordinates is the best set
which can be made from the simple coordinates. It follows
the advice given above for five membered rings, and because
it includes the OO bond, has no trouble with crashing this
bond. It takes 20 steps to converge, so the trouble of
generating this $ZMAT is certainly worth it compared to the
use of Cartesians.
Natural internal coordinates are defined in the final
group of input. The coordinates are set up first for the
ring, including two linear combinations of all angles and
all torsions withing the ring. After this the methyl is
hooked to the ring as if it were a NH group, using the
usual terminal methyl hydrogen definitions. The H is
hooked to this same ring carbon as if it were a methine.
The NH and the CH2 within the ring follow Pulay's rules
exactly. The amount of input is much greater than a normal
Z-matrix. For example, 46 internal coordinates are given,
which are then placed in 3N-6=33 linear combinations. Note
that natural internals tend to be rich in bends, and short
on torsions.
The energy results for the three coordinate systems
which converge are as follows:
NSERCH Cart good Z-mat nat. int.
0 -48.6594935049 -48.6594935049 -48.6594935049
1 -48.6800538676 -48.6806631261 -48.6838361406
2 -48.6822702585 -48.6510215698 -48.6874045449
3 -48.6858299354 -48.6882945647 -48.6932811528
4 -48.6881499412 -48.6849667085 -48.6946836332
5 -48.6890226067 -48.6911899936 -48.6959800274
6 -48.6898261650 -48.6878047907 -48.6973821465
7 -48.6901936624 -48.6930608185 -48.6987652146
8 -48.6905304889 -48.6940607117 -48.6996366016
9 -48.6908626791 -48.6949137185 -48.7006656309
10 -48.6914279465 -48.6963767038 -48.7017273728
11 -48.6921521142 -48.6986608776 -48.7021504975
12 -48.6931136707 -48.7007305310 -48.7022405019
13 -48.6940437619 -48.7016095285 -48.7022548935
14 -48.6949546487 -48.7021531692 -48.7022569328
15 -48.6961698826 -48.7022080183 -48.7022570260
16 -48.6973813002 -48.7022454522 -48.7022570662
17 -48.6984850655 -48.7022492840
18 -48.6991553826 -48.7022503853
19 -48.6996239136 -48.7022507037
20 -48.7002269303 -48.7022508393
21 -48.7005379631
22 -48.7008387759
...
50 -48.7022519950
from which you can see that the natural internals are
actually the best set. The $ZMAT exhibits upward burps
in the energy at step 2, 4, and 6, so that for the
same number of steps, these coordinates are always at a
higher energy than the natural internals.
The initial hessian generated for these three columns
contains 0, 33, and 46 force constants. This assists
the natural internals, but is not the major reason for
its superior performance. The computed hessian at the
final geometry of this molecule, when transformed into the
natural internal coordinates is almost diagonal. This
almost complete uncoupling of coordinates is what makes
the natural internals perform so well. The conclusion
is of course that not all coordinate systems are equal,
and natural internals are the best. As another example,
we have run the ATCHCP molecule, which is a popular
geometry optimization test, due to its two fused rings:
H.B.Schlegel, Int.J.Quantum Chem., Symp. 26, 253-264(1992)
T.H.Fischer and J.Almlof, J.Phys.Chem. 96, 9768-9774(1992)
J.Baker, J.Comput.Chem. 14, 1085-1100(1993)
Here we have compared the same coordinate types, using a
guess hessian, or a computed hessian. The latter set of
runs is a test of the coordinates only, as the initial
hessian information is identical. The results show clearly
the superiority of the natural internals, which like the
previous example, give an energy decrease on every step:
HESS=GUESS HESS=READ
Cartesians 65 41 steps
good Z-matrix 32 23
natural internals 24 13
A final example is phosphinoazasilatrane, with three rings
fused on a common SiN bond, in which 112 steps in Cartesian
space became 32 steps in natural internals. The moral is:
"A little brain time can save a lot of CPU time."
In late 1998, a new kind of internal coordinate method
was included into GAMESS. This is the delocalized
internal
coordinate (DLC) of
J.Baker, A. Kessi, B.Delley
J.Chem.Phys. 105, 192-212(1996)
although as is the usual case, the implementation is not
exactly the same. Bonds are kept as independent
coordinates,
while angles are placed in linear combination by the DLC
process. There are some interesting options for applying
constraints, and other options to assist the automatic DLC
generation code by either adding or deleting coordinates.
It is simple to use DLCs in their most basic form:
$contrl nzvar=xx $end
$zmat dlc=.true. auto=.true. $end
Our initial experience is that the quality of DLCs is
not as good as explicitly constructed natural internals,
which benefit from human chemical knowledge, but are almost
always better than carefully crafted $ZMATs using only the
primitive internal coordinates (although we have seen a few
exceptions). Once we have more numerical experience with
the use of DLC's, we will come back and revise the above
discussion of coordinate choices. In the meantime, they
are quite simple to choose, so give them a go.
the role of symmetry
At the end of a succesful geometry search, you will
have a set of coordinates where the gradient of the energy
is zero. However your newly discovered stationary point
is not necessarily a minimum or saddle point!
This apparent mystery is due to the fact that the
gradient vector transforms under the totally symmetric
representation of the molecular point group. As a direct
consequence, a geometry search is point group conserving.
(For a proof of these statements, see J.W.McIver and
A.Komornicki, Chem.Phys.Lett., 10,303-306(1971)). In
simpler terms, the molecule will remain in whatever point
group you select in $DATA, even if the true minimum is in
some lower point group. Since a geometry search only
explores totally symmetric degrees of freedom, the only
way to learn about the curvatures for all degrees of
freedom is RUNTYP=HESSIAN.
As an example, consider disilene, the silicon analog
of ethene. It is natural to assume that this molecule is
planar like ethene, and an OPTIMIZE run in D2h symmetry
will readily locate a stationary point. However, as a
calculation of the hessian will readily show, this
structure is a transition state (one imaginary frequency),
and the molecule is really trans-bent (C2h). A careful
worker will always characterize a stationary point as
either a minimum, a transition state, or some higher order
stationary point (which is not of great interest!) by
performing a RUNTYP=HESSIAN.
The point group conserving properties of a geometry
search can be annoying, as in the preceeding example, or
advantageous. For example, assume you wish to locate the
transition state for rotation about the double bond in
ethene. A little thought will soon reveal that ethene is
D2h, the 90 degrees twisted structure is D2d, and
structures in between are D2. Since the saddle point is
actually higher symmetry than the rest of the rotational
surface, you can locate it by RUNTYP=OPTIMIZE within D2d
symmetry. You can readily find this stationary point with
the diagonal guess hessian! In fact, if you attempt to do
a RUNTYP=SADPOINT within D2d symmetry, there will be no
totally symmetric modes with negative curvatures, and it
is unlikely that the geometry search will be very well
behaved.
Although a geometry search cannot lower the symmetry,
the gain of symmetry is quite possible. For example, if
you initiate a water molecule optimization with a trial
structure which has unequal bond lengths, the geometry
search will come to a structure that is indeed C2v (to
within OPTTOL, anyway). However, GAMESS leaves it up to
you to realize that a gain of symmetry has occurred.
In general, Mother Nature usually chooses more
symmetrical structures over less symmetrical structures.
Therefore you are probably better served to assume the
higher symmetry, perform the geometry search, and then
check the stationary point's curvatures. The alternative
is to start with artificially lower symmetry and see if
your system regains higher symmetry. The problem with
this approach is that you don't necessarily know which
subgroup is appropriate, and you lose the great speedups
GAMESS can obtain from proper use of symmetry. It is good
to note here that "lower symmetry" does not mean simply
changing the name of the point group and entering more
atoms in $DATA, instead the nuclear coordinates themselves
must actually be of lower symmetry.
practical matters
Geometry searches do not bring the gradient exactly to
zero. Instead they stop when the largest component of the
gradient is below the value of OPTTOL, which defaults to
a reasonable 0.0001. Analytic hessians usually have
residual frequencies below 10 cm**-1 with this degree of
optimization. The sloppiest value you probably ever want
to try is 0.0005.
If a geometry search runs out of time, or exceeds
NSTEP, it can be restarted. For RUNTYP=OPTIMIZE, restart
with the coordinates having the lowest total energy
(do a string search on "FINAL"). For RUNTYP=SADPOINT,
restart with the coordinates having the smallest gradient
(do a string search on "RMS", which means root mean
square).
These are not necessarily at the last geometry!
The "restart" should actually be a normal run, that is
you should not try to use the restart options in $CONTRL
(which may not work anyway). A geometry search can be
restarted by extracting the desired coordinates for $DATA
from the printout, and by extracting the corresponding
$GRAD group from the PUNCH file. If the $GRAD group is
supplied, the program is able to save the time it would
ordinarily take to compute the wavefunction and gradient
at the initial point, which can be a substantial savings.
There is no input to trigger reading of a $GRAD group: if
found, it is read and used. Be careful that your $GRAD
group actually corresponds to the coordinates in $DATA, as
GAMESS has no check for this.
Sometimes when you are fairly close to the minimum, an
OPTIMIZE run will take a first step which raises the
energy, with subsequent steps improving the energy and
perhaps finding the minimum. The erratic first step is
caused by the GUESS hessian. It may help to limit the size
of this wrong first step, by reducing its radius, DXMAX.
Conversely, if you are far from the minimum, sometimes you
can decrease the number of steps by increasing DXMAX.
When using internals, the program uses an iterative
process to convert the internal coordinate change into
Cartesian space. In some cases, a small change in the
internals will produce a large change in Cartesians, and
thus produce a warning message on the output. If these
warnings appear only in the beginning, there is probably
no problem, but if they persist you can probably devise
a better set of coordinates. You may in fact have one of
the two problems described in the next paragraph. In
some cases (hopefully very few) the iterations to find
the Cartesian displacement may not converge, producing a
second kind of warning message. The fix for this may
very well be a new set of internal coordinates as well,
or adjustment of ITBMAT in $STATPT.
There are two examples of poorly behaved internal
coordinates which can give serious problems. The first
of these is three angles around a central atom, when
this atom becomes planar (sum of the angles nears 360).
The other is a dihedral where three of the atoms are
nearly linear, causing the dihedral to flip between 0 and
180. Avoid these two situations if you want your geometry
search to be convergent.
Sometimes it is handy to constrain the geometry search
by freezing one or more coordinates, via the IFREEZ array.
For example, constrained optimizations may be useful while
trying to determine what area of a potential energy surface
contains a saddle point. If you try to freeze coordinates
with an automatically generated $ZMAT, you need to know
that the order of the coordinates defined in $DATA is
y
y x r1
y x r2 x a3
y x r4 x a5 x w6
y x r7 x a8 x w9
and so on, where y and x are whatever atoms and molecular
connectivity you happen to be using.
saddle points
Finding minima is relatively easy. There are large
tables of bond lengths and angles, so guessing starting
geometries is pretty straightforward. Very nasty cases
may require computation of an exact hessian, but the
location of most minima is straightforward.
In contrast, finding saddle points is a black art.
The diagonal guess hessian will never work, so you must
provide a computed one. The hessian should be computed at
your best guess as to what the transition state (T.S.)
should be. It is safer to do this in two steps as outlined
above, rather than HESS=CALC. This lets you verify you
have guessed a structure with one and only one negative
curvature. Guessing a good trial structure is the hardest
part of a RUNTYP=SADPOINT!
This point is worth iterating. Even with sophisticated
step size control such as is offered by the QA/TRIM or RFO
methods, it is in general very difficult to move correctly
from a region with incorrect curvatures towards a saddle
point. Even procedures such as CONOPT or RUNTYP=GRADEXTR
will not replace your own chemical intuition about where
saddle points may be located.
The RUNTYP=HESSIAN's normal coordinate analysis is
rigorously valid only at stationary points on the surface.
This means the frequencies from the hessian at your trial
geometry are untrustworthy, in particular the six "zero"
frequencies corresponding to translational and rotational
(T&R) degrees of freedom will usually be 300-500 cm**-1,
and possibly imaginary. The Sayvetz conditions on the
printout will help you distinguish the T&R "contaminants"
from the real vibrational modes. If you have defined a
$ZMAT, the PURIFY option within $STATPT will help zap out
these T&R contaminants).
If the hessian at your assumed geometry does not have
one and only one imaginary frequency (taking into account
that the "zero" frequencies can sometimes be 300i!), then
it will probably be difficult to find the saddle point.
Instead you need to compute a hessian at a better guess
for the initial geometry, or read about mode following
below.
If you need to restart your run, do so with the
coordinates which have the smallest RMS gradient. Note
that the energy does not necessarily have to decrease in a
SADPOINT run, in contrast to an OPTIMIZE run. It is often
necessary to do several restarts, involving recomputation
of the hessian, before actually locating the saddle point.
Assuming you do find the T.S., it is always a good
idea to recompute the hessian at this structure. As
described in the discussion of symmetry, only totally
symmetric vibrational modes are probed in a geometry
search. Thus it is fairly common to find that at your
"T.S." there is a second imaginary frequency, which
corresponds to a non-totally symmetric vibration. This
means you haven't found the correct T.S., and are back to
the drawing board. The proper procedure is to lower the
point group symmetry by distorting along the symmetry
breaking "extra" imaginary mode, by a reasonable amount.
Don't be overly timid in the amount of distortion, or the
next run will come back to the invalid structure.
The real trick here is to find a good guess for the
transition structure. The closer you are, the better. It
is often difficult to guess these structures. One way
around this is to compute a linear least motion (LLM)
path. This connects the reactant structure to the product
structure by linearly varying each coordinate. If you
generate about ten structures intermediate to reactants
and products, and compute the energy at each point, you
will in general find that the energy first goes up, and
then down. The maximum energy structure is a "good" guess
for the true T.S. structure. Actually, the success of
this method depends on how curved the reaction path is.
A particularly good paper on the symmetry which a
saddle point (and reaction path) can possess is by
P.Pechukas, J.Chem.Phys. 64, 1516-1521(1976)
mode following
In certain circumstances, METHOD=RFO and QA can walk
from a region of all positive curvatures (i.e. near a
minimum) to a transition state. The criteria for whether
this will work is that the mode being followed should be
only weakly coupled to other close-lying Hessian modes.
Especially, the coupling to lower modes should be almost
zero. In practise this means that the mode being followed
should be the lowest of a given symmetry, or spatially far
away from lower modes (for example, rotation of methyl
groups at different ends of the molecule). It is certainly
possible to follow also modes which do not obey these
criteria, but the resulting walk (and possibly TS location)
will be extremely sensitive to small details such as the
stepsize.
This sensitivity also explain why TS searches often
fail, even when starting in a region where the Hessian has
the required one negative eigenvalue. If the TS mode is
strongly coupled to other modes, the direction of the mode
is incorrect, and the maximization of the energy along
that direction is not really what you want (but what you
get).
Mode following is really not a substitute for the
ability to intuit regions of the PES with a single local
negative curvature. When you start near a minimum, it
matters a great deal which side of the minima you start
from, as the direction of the search depends on the sign
of the gradient. We strongly urge that you read before
trying to use IFOLOW, namely the papers by Frank Jensen
and Jon Baker mentioned above, and see also Figure 3 of
C.J.Tsai, K.D.Jordan, J.Phys.Chem. 97, 11227-11237 (1993)
which is quite illuminating on the sensitivity of mode
following to the initial geometry point.
Note that GAMESS retains all degrees of freedom in its
hessian, and thus there is no reason to suppose the lowest
mode is totally symmetric. Remember to lower the symmetry
in the input deck if you want to follow non-symmetric
modes. You can get a printout of the modes in internal
coordinate space by a EXETYP=CHECK run, which will help
you decide on the value of IFOLOW.
* * *
CONOPT is a different sort of saddle point search
procedure. Here a certain "CONstrained OPTimization" may
be considered as another mode following method. The idea
is to start from a minimum, and then perform a series of
optimizations on hyperspheres of increasingly larger
radii. The initial step is taken along one of the Hessian
modes, chosen by IFOLOW, and the geometry is optimized
subject to the constraint that the distance to the minimum
is constant. The convergence criteria for the gradient
norm perpendicular to the constraint is taken as 10*OPTTOL,
and the corresponding steplength as 100*OPTTOL.
After such a hypersphere optimization has converged, a
step is taken along the line connecting the two previous
optimized points to get an estimate of the next hyper-
sphere geometry. The stepsize is DXMAX, and the radius of
hyperspheres is thus increased by an amount close (but not
equal) to DXMAX. Once the pure NR step size falls below
DXMAX/2 or 0.10 (whichever is the largest) the algorithm
switches to a straight NR iterate to (hopefully) converge
on the stationary point.
The current implementation always conducts the search
in cartesian coordinates, but internal coordinates may be
printed by the usual specification of NZVAR and ZMAT. At
present there is no restart option programmed.
CONOPT is based on the following papers, but the actual
implementation is the modified equations presented in
Frank Jensen's paper mentioned above.
Y. Abashkin, N. Russo,
J.Chem.Phys. 100, 4477-4483(1994).
Y. Abashkin, N. Russo, M. Toscano,
Int.J.Quant.Chem. 52, 695-704(1994).
There is little experience on how this method works in
practice, experiment with it at your own risk!
Intrinsic Reaction Coordinate Methods
The Intrinsic Reaction Coordinate (IRC) is defined as
the minimum energy path connecting the reactants to
products via the transition state. In practice, the IRC is
found by first locating the transition state for the
reaction. The IRC is then found in halves, going forward
and backwards from the saddle point, down the steepest
descent path in mass weighted Cartesian coordinates. This
is accomplished by numerical integration of the IRC
equations, by a variety of methods to be described below.
The IRC is becoming an important part of polyatomic
dynamics research, as it is hoped that only knowledge of
the PES in the vicinity of the IRC is needed for prediction
of reaction rates, at least at threshhold energies. The
IRC has a number of uses for electronic structure purposes
as well. These include the proof that a certain transition
structure does indeed connect a particular set of reactants
and products, as the structure and imaginary frequency
normal mode at the saddle point do not always unambiguously
identify the reactants and products. The study of the
electronic and geometric structure along the IRC is also of
interest. For example, one can obtain localized orbitals
along the path to determine when bonds break or form.
The accuracy to which the IRC is determined is dictated
by the use one intends for it. Dynamical calculations
require a very accurate determination of the path, as
derivative information (second derivatives of the PES at
various IRC points, and path curvature) is required later.
Thus, a sophisticated integration method (such as GS2), and
small step sizes (STRIDE=0.05, 0.01, or even smaller) may
be needed. In addition to this, care should be taken to
locate the transition state carefully (perhaps decreasing
OPTTOL by a factor of 10, to OPTTOL=1D-5), and in the
initiation of the IRC run. The latter might require a
hessian matrix obtained by double differencing, certainly
the hessian should be PROJCT'd or PURIFY'd. Note also that
EVIB must be chosen carefully, as decribed below.
On the other hand, identification of reactants and
products allows for much larger step sizes, and cruder
integration methods. In this type of IRC one might want to
be careful in leaving the saddle point (perhaps STRIDE
should be reduced to 0.10 or 0.05 for the first few steps
away from the transition state), but once a few points have
been taken, larger step sizes can be employed. In general,
the defaults in the $IRC group are set up for this latter,
cruder quality IRC. The STRIDE value for the GS2 method
can usually be safely larger than for other methods, no
matter what your interest in accuracy is.
The next few paragraphs describe the various
integrators, but note that GS2 is superior to the others.
The simplest method of determining an IRC is linear
gradient following, PACE=LINEAR. This method is also known
as Euler's method. If you are employing PACE=LINEAR, you
can select "stabilization" of the reaction path by the
Ishida, Morokuma, Komornicki method. This type of
corrector has no apparent mathematical basis, but works
rather well since the bisector usually intersects the
reaction path at right angles (for small step sizes). The
ELBOW variable allows for a method intermediate to LINEAR
and stabilized LINEAR, in that the stabilization will be
skipped if the gradients at the original IRC point, and at
the result of a linear prediction step form an angle
greater than ELBOW. Set ELBOW=180 to always perform the
stabilization.
A closely related method is PACE=QUAD, which fits a
quadratic polynomial to the gradient at the current and
immediately previous IRC point to predict the next point.
This pace has the same computational requirement as LINEAR,
and is slightly more accurate due to the reuse of the old
gradient. However, stabilization is not possible for this
pace, thus a stabilized LINEAR path is usually more
accurate than QUAD.
Two rather more sophisticated methods for integrating
the IRC equations are the fourth order Adams-Moulton
predictor-corrector (PACE=AMPC4) and fourth order Runge-
Kutta (PACE=RK4). AMPC4 takes a step towards the next IRC
point (prediction), and based on the gradient found at this
point (in the near vincinity of the next IRC point) obtains
a modified step to the desired IRC point (correction).
AMPC4 uses variable step sizes, based on the input STRIDE.
RK4 takes several steps part way toward the next IRC point,
and uses the gradient at these points to predict the next
IRC point. RK4 is one of the most accurate integration
method implemented in GAMESS, and is also the most time
consuming.
The Gonzalez-Schlegel 2nd order method (PACE=GS2) finds
the next IRC point by a constrained optimization on the
surface of a hypersphere, centered at a point 1/2 STRIDE
along the gradient vector leading from the previous IRC
point. By construction, the reaction path between two
successive IRC points is a circle tangent to the two
gradient vectors. The algorithm is much more robust for
large steps than the other methods, so it has been chosen
as the default method. Thus, the default for STRIDE is too
large for the other methods. The number of energy and
gradients need to find the next point varies with the
difficulty of the constrained optimization, but is
normally not very many points. Taking more than 2-3 steps
in this constrained optimization is indicative of reaction
path curvature, and thus it may help to reduce the step
size. Use a small GCUT (same value as OPTTOL) when trying
to integrate an IRC very accurately, to be sure the
hypersphere optimizations are well converged. Be sure to
provide the updated hessian from the previous run when
restarting PACE=GS2.
The number of wavefunction evaluations, and energy
gradients needed to jump from one point on the IRC to the
next point are summarized in the following table:
PACE # energies # gradients
---- ---------- -----------
LINEAR 1 1
stabilized
LINEAR 3 2
QUAD 1 1 (+ reuse of historical
gradient)
AMPC4 2 2 (see note)
RK4 4 4
GS2 2-4 2-4 (equal numbers)
Note that the AMPC4 method sometimes does more than one
correction step, with each such correction adding one more
energy and gradient to the calculation. You get what you
pay for in IRC calculations: the more energies and
gradients which are used, the more accurate the path found.
A description of these methods, as well as some others
that were found to be not as good is geven by Kim Baldridge
and Lisa Pederson, Pi Mu Epsilon J., 9, 513-521 (1993).
* * *
All methods are initiated by jumping from the saddle
point, parallel to the normal mode (CMODE) which has an
imaginary frequency. The jump taken is designed to lower
the energy by an amount EVIB. The actual distance taken is
thus a function of the imaginary frequency, as a smaller
FREQ will produce a larger initial jump. You can simply
provide a $HESS group instead of CMODE and FREQ, which
involves less typing. To find out the actual step taken
for a given EVIB, use EXETYP=CHECK. The direction of the
jump (towards reactants or products) is governed by FORWRD.
Note that if you have decided to use small step sizes, you
must employ a smaller EVIB to ensure a small first step.
The GS2 method begins by following the normal mode by one
half of STRIDE, and then performing a hypersphere
minimization about that point, so EVIB is irrelevant to
this PACE.
The only method which proves that a properly converged
IRC has been obtained is to regenerate the IRC with a
smaller step size, and check that the IRC is unchanged.
Again, note that the care with which an IRC must be
obtained is highly dependent on what use it is intended
for.
Some key IRC references are:
K.Ishida, K.Morokuma, A.Komornicki
J.Chem.Phys. 66, 2153-2156 (1977)
K.Muller
Angew.Chem., Int.Ed.Engl. 19, 1-13 (1980)
M.W.Schmidt, M.S.Gordon, M.Dupuis
J.Am.Chem.Soc. 107, 2585-2589 (1985)
B.C.Garrett, M.J.Redmon, R.Steckler, D.G.Truhlar,
K.K.Baldridge, D.Bartol, M.W.Schmidt, M.S.Gordon
J.Phys.Chem. 92, 1476-1488(1988)
K.K.Baldridge, M.S.Gordon, R.Steckler, D.G.Truhlar
J.Phys.Chem. 93, 5107-5119(1989)
C.Gonzalez, H.B.Schlegel
J.Chem.Phys. 90, 2154-2161(1989)
The IRC discussion closes with some practical tips:
The $IRC group has a confusing array of variables, but
fortunately very little thought need be given to most of
them. An IRC run is restarted by moving the coordinates of
the next predicted IRC point into $DATA, and inserting the
new $IRC group into your input file. You must select the
desired value for NPOINT. Thus, only the first job which
initiates the IRC requires much thought about $IRC.
The symmetry specified in the $DATA deck should be the
symmetry of the reaction path. If a saddle point happens
to have higher symmetry, use only the lower symmetry in the
$DATA deck when initiating the IRC. The reaction path will
have a lower symmetry than the saddle point whenever the
normal mode with imaginary frequency is not totally
symmetric. Be careful that the order and orientation of
the atoms corresponds to that used in the run which
generated the hessian matrix.
If you wish to follow an IRC for a different isotope,
use the $MASS group. If you wish to follow the IRC in
regular Cartesian coordinates, just enter unit masses for
each atom. Note that CMODE and FREQ are a function of the
atomic masses, so either regenerate FREQ and CMODE, or more
simply, provide the correct $HESS group.
Gradient Extremals
This section of the manual, as well as the source code
to trace gradient extremals was written by Frank Jensen of
the University of Aarhus.
A Gradient Extremal (GE) curve consists of points where
the gradient norm on a constant energy surface is
stationary. This is equivalent to the condition that
the gradient is an eigenvector of the Hessian. Such GE
curves radiate along all normal modes from a stationary
point, and the GE leaving along the lowest normal mode
from a minimum is the gentlest ascent curve. This is not
the same as the IRC curve connecting a minimum and a TS,
but may in some cases be close.
GEs may be divided into three groups: those leading
to dissociation, those leading to atoms colliding, and
those which connect stationary points. The latter class
allows a determination of many (all?) stationary points on
a PES by tracing out all the GEs. Following GEs is thus a
semi-systematic way of mapping out stationary points. The
disadvantages are:
i) There are many (but finitely many!) GEs for a
large molecule.
ii) Following GEs is computationally expensive.
iii) There is no control over what type of
stationary point (if any) a GE will lead to.
Normally one is only interested in minima and TSs, but
many higher order saddle points will also be found.
Furthermore, it appears that it is necessary to follow GEs
radiating also from TSs and second (and possibly also
higher) order saddle point to find all the TSs.
A rather complete map of the extremals for the H2CO
potential surface is available in a paper which explains
the points just raised in greater detail:
K.Bondensgaard, F.Jensen,
J.Chem.Phys. 104, 8025-8031(1996).
An earlier paper gives some of the properties of GEs:
D.K.Hoffman, R.S.Nord, K.Ruedenberg,
Theor. Chim. Acta 69, 265-279(1986).
There are two GE algorithms in GAMESS, one due to Sun
and Ruedenberg (METHOD=SR), which has been extended to
include the capability of locating bifurcation points and
turning points, and another due to Jorgensen, Jensen, and
Helgaker (METHOD=JJH):
J. Sun, K. Ruedenberg, J.Chem.Phys. 98, 9707-9714(1993)
P. Jorgensen, H. J. Aa. Jensen, T. Helgaker
Theor. Chim. Acta 73, 55 (1988).
The Sun and Ruedenberg method consist of a predictor
step taken along the tangent to the GE curve, followed by
one or more corrector steps to bring the geometry back to
the GE. Construction of the GE tangent and the corrector
step requires elements of the third derivative of the
energy, which is obtained by a numerical differentiation
of two Hessians. This puts some limitations on which
systems the GE algorithm can be used for. First, the
numerical differentiation of the Hessian to produce third
derivatives means that the Hessian should be calculated by
analytical methods, thus only those types of wavefunctions
where this is possible can be used. Second, each
predictor/corrector step requires at least two Hessians,
but often more. Maybe 20-50 such steps are necessary for
tracing a GE from one stationary point to the next. A
systematic study of all the GE radiating from a stationary
point increases the work by a factor of ~2*(3N-6). One
should thus be prepared to invest at least hundreds, and
more likely thousands, of Hessian calculations. In other
words, small systems, small basis sets, and simple wave-
functions.
The Jorgensen, Jensen, and Helgaker method consists of
taking a step in the direction of the chosen Hessian
eigenvector, and then a pure NR step in the perpendicular
modes. This requires (only) one Hessian calculation for
each step. It is not suitable for following GEs where the
GE tangent forms a large angle with the gradient, and it
is incapable of locating GE bifurcations.
Although experience is limited at present, the JJH
method does not appear to be suitable for following GEs in
general (at least not in the current implementation).
Experiment with it at your own risk!
The flow of the SR algorithm is as follows: A
predictor geometry is produced, either by jumping away
from a stationary point, or from a step in the tangent
direction from the previous point on the GE. At the
predictor geometry, we need the gradient, the Hessian, and
the third derivative in the gradient direction. Depending
on HSDFDB, this can be done in two ways. If .TRUE. the
gradient is calculated, 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. This corresponds to a double-
sided differentiation, and is the numerical most stable
method for getting the partial third derivative matrix.
If HSDFDB = .FALSE., the gradient and Hessian are
calculated at the current geometry, and one additional
Hessian is calculated at SNUMH distance in the gradient
direction. This corresponds to a single-sided differen-
tiation. In both cases, two full Hessian calculations are
necessary, but HSDFDB = .TRUE. require one additional
wavefunction and gradient calculation. This is usually
a fairly small price compared to two Hessians, and the
numerically better double-sided differentiation has
therefore been made the default.
Once the gradient, Hessian, and third derivative is
available, the corrector step and the new GE tangent are
constructed. If the corrector step is below a threshold,
a new predictor step is taken along the tangent vector.
If the corrector step is larger than the threshold, the
correction step is taken, and a new micro iteration is
performed. DELCOR thus determines how closely the GE will
be followed, and DPRED determine how closely the GE path
will be sampled.
The construction of the GE tangent and corrector step
involve solution of a set of linear equations, which in
matrix notation can be written as Ax=B. The A-matrix is
also the second derivative of the gradient norm on the
constant energy surface.
After each corrector step, various things are printed
to monitor the behavior: The projection of the gradient
along the Hessian eigenvalues (the gradient is parallel
to an eigenvector on the GE), the projection of the GE
tangent along the Hessian eigenvectors, and the overlap
of the Hessian eigenvectors with the mode being followed
from the previous (optimzed) geometry. The sign of these
overlaps are not significant, they just refer to an
arbitrary phase of the Hessian eigenvectors.
After the micro iterations has converged, the Hessian
eigenvector curvatures are also displayed, this is an
indication of the coupling between the normal modes. The
number of negative eigenvalues in the A-matrix is denoted
the GE index. If it changes, one of the eigenvalues must
have passed through zero. Such points may either be GE
bifurcations (where two GEs cross) or may just be "turning
points", normally when the GE switches from going uphill
in energy to downhill, or vice versa. The distinction is
made based on the B-element corresponding to the A-matrix
eigenvalue = 0. If the B-element = 0, it is a bifurcation,
otherwise it is a turning point.
If the GE index changes, a linear interpolation is
performed between the last two points to locate the point
where the A-matrix is singular, and the corresponding
B-element is determined. The linear interpolation points
will in general be off the GE, and thus the evaluation of
whether the B-element is 0 is not always easy. The
program additionally evaluates the two limiting vectors
which are solutions to the linear sets of equations, these
are also used for testing whether the singular point is a
bifurcation point or turning point.
Very close to a GE bifurcation, the corrector step
become numerically unstable, but this is rarely a problem
in practice. It is a priori expected that GE bifurcation
will occur only in symmetric systems, and the crossing GE
will break the symmetry. Equivalently, a crossing GE may
be encountered when a symmetry element is formed, however
such crossings are much harder to detect since the GE
index does not change, as one of the A-matrix eigenvalues
merely touches zero. The program prints an message if
the absolute value of an A-matrix eigenvalue reaches a
minimum near zero, as such points may indicate the
passage of a bifurcation where a higher symmetry GE
crosses. Run a movie of the geometries to see if a more
symmetric structure is passed during the run.
An estimate of the possible crossing GE direction is
made at all points where the A-matrix is singular, and two
perturbed geometries in the + and - direction are written
out. These may be used as predictor geometries for
following a crossing GE. If the singular geometry is a
turning point, the + and - geometries are just predictor
geometries on the GE being followed.
In any case, a new predictor step can be taken to trace
a different GE from the newly discovered singular point,
using the direction determined by interpolation from the
two end point tangents (the GE tangent cannot be uniquely
determined at a bifurcation point). It is not possible to
determine what the sign of IFOLOW should be when starting
off along a crossing GE at a bifurcation, one will have to
try a step to see if it returns to the bifurcation point
or not.
In order to determine whether the GE index change it
is necessary to keep track of the order of the A-matrix
eigenvalues. The overlap between successive eigenvectors
are shown as "Alpha mode overlaps".
Things to watch out for:
1) The numerical differentiation to get third derivatives
requires more accuracy than usual. The SCF convergence
should be at least 100 times smaller than SNUMH, and
preferably better. With the default SNUMH of 10**(-4)
the SCF convergence should be at least 10**(-6). Since
the last few SCF cycles are inexpensive, it is a good idea
to tighten the SCF convergence as much as possible, to
maybe 10**(-8) or better. You may also want to increase
the integral accuracy by reducing the cutoffs (ITOL and
ICUT) and possibly also try more accurate integrals
(INTTYP=HONDO). The CUTOFF in $TRNSFM may also be reduced
to produce more accurate Hessians. Don't attempt to use a
value for SNUMH below 10**(-6), as you simply can't get
enough accuracy. Since experience is limited at present,
it is recommended that some tests runs are made to learn
the sensitivity of these factors for your system.
2) GEs can be followed in both directions, uphill or
downhill. When stating from a stationary point, the
direction is implicitly given as away from the stationary
point. When starting from a non-stationary point, the "+"
and "-" directions (as chosen by the sign of IFOLOW)
refers to the gradient direction. The "+" direction is
along the gradient (energy increases) and "-" is opposite
to the gradient (energy decreases).
3) A switch from one GE to another may be seen when two
GE come close together. This is especially troublesome
near bifurcation points where two GEs actually cross. In
such cases a switch to a GE with -higher- symmetry may
occur without any indication that this has happened,
except possibly that a very large GE curvature suddenly
shows up. Avoid running the calculation with less
symmetry than the system actually has, as this increases
the likelihood that such switches occuring. Fix: alter
DPRED to avoid having the predictor step close to the
crossing GE.
4) "Off track" error message: The Hessian eigenvector
which is parallel to the gradient is not the same as
the one with the largest overlap to the previous
Hessian mode. This usually indicate that a GE switch
has occured (note that a switch may occur without this
error message), or a wrong value for IFOLOW when starting
from a non-stationary point. Fix: check IFOLOW, if it is
correct then reduce DPRED, and possibly also DELCOR.
5) Low overlaps of A-matrix eigenvectors. Small overlaps
may give wrong assignment, and wrong conclusions about GE
index change. Fix: reduce DPRED.
6) The interpolation for locating a point where one of the
A-matrix eigenvalues is zero fail to converge. Fix:
reduce DPRED (and possibly also DELCOR) to get a shorther
(and better) interpolation line.
7) The GE index changes by more than 1. A GE switch may
have occured, or more than one GE index change is located
between the last and current point. Fix: reduce DPRED to
sample the GE path more closely.
8) If SNRMAX is too large the algorithm may try to locate
stationary points which are not actually on the GE being
followed. Since GEs often pass quite near a stationary
point, SNRMAX should only be increased above the default
0.10 after some consideration.
Continuum Solvation Methods
In a very thorough 1994 review of continuum solvation
models, Tomasi and Persico divide the possible approaches
to the treatment of solvent effects into four categories:
a) virial equations of state, correlation functions
b) Monte Carlo or molecular dynamics simulations
c) continuum treatments
d) molecular treatments
The Effective Fragment Potential method, documented in the
following section of this chapter, falls into the latter
category, as each EFP solvent molecule is modeled as a
distinct object (discrete solvation). This section
describes the four continuum models which are implemented
in the standard version of GAMESS, and a fifth model which
can be interfaced.
Continuum models typically form a cavity of some sort
containing the solute molecule, while the solvent outside
the cavity is thought of as a continuous medium and is
categorized by a limited amount of physical data, such as
the dielectric constant. The electric field of the charged
particles comprising the solute interact with this
background medium, producing a polarization in it, which in
turn feeds back upon the solute's wavefunction.
Self Consistent Reaction Field (SCRF)
A simple continuum model is the Onsager cavity model,
often called the Self-Consistent Reaction Field, or SCRF
model. This represents the charge distribution of the
solute in terms of a multipole expansion. SCRF usually
uses an idealized cavity (spherical or ellipsoidal) to
allow an analytic solution to the interaction energy
between the solute multipole and the multipole which this
induces in the continuum. This method is implemented in
GAMESS in the simplest possible fashion:
i) a spherical cavity is used
ii) the molecular electrostatic potential of the
solute is represented as a dipole only, except
a monopole is also included for an ionic solute.
The input for this implementation of the Kirkwood-Onsager
model is provided in $SCRF.
Some references on the SCRF method are
1. J.G.Kirkwood J.Chem.Phys. 2, 351 (1934)
2. L.Onsager J.Am.Chem.Soc. 58, 1486 (1936)
3. O.Tapia, O.Goscinski Mol.Phys. 29, 1653 (1975)
4. M.M.Karelson, A.R.Katritzky, M.C.Zerner
Int.J.Quantum Chem., Symp. 20, 521-527 (1986)
5. K.V.Mikkelsen, H.Agren, H.J.Aa.Jensen, T.Helgaker
J.Chem.Phys. 89, 3086-3095 (1988)
6. M.W.Wong, M.J.Frisch, K.B.Wiberg
J.Am.Chem.Soc. 113, 4776-4782 (1991)
7. M.Szafran, M.M.Karelson, A.R.Katritzky, J.Koput,
M.C.Zerner J.Comput.Chem. 14, 371-377 (1993)
8. M.Karelson, T.Tamm, M.C.Zerner
J.Phys.Chem. 97, 11901-11907 (1993)
The method is very sensitive to the choice of the solute
RADIUS, but not very sensitive to the particular DIELEC of
polar solvents. The plots in reference 7 illustrate these
points very nicely. The SCRF implementation in GAMESS is
Zerner's Method A, described in the same reference. The
total solute energy includes the Born term, if the solute
is an ion. Another limitation is that a solute's electro-
static potential is not likely to be fit well as a dipole
moment only, for example see Table VI of reference 5 which
illustrates the importance of higher multipoles. Finally,
the restriction to a spherical cavity may not be very
representative of the solute's true shape. However, in the
special case of a roundish molecule, and a large dipole
which is geometry sensitive, the SCRF model may include
sufficient physics to be meaningful:
M.W.Schmidt, T.L.Windus, M.S.Gordon
J.Am.Chem.Soc. 117, 7480-7486(1995).
Most cases should choose PCM (next section) over SCRF!!!
Polarizable Continuum Model (PCM)
A much more sophisticated continuum method, named the
Polarizable Continuum Model, is also available. The PCM
method places a solute in a cavity formed by a union of
spheres centered on each atom. PCM includes a more exact
treatment of the electrostatic interaction of the solute
with the surrounding medium, on the cavity's surface. The
computational procedure divides this surface into many
small tesserae, each having a different "apparent surface
charge", reflecting the solute's and other tesserae's
electric field at each. These surface charges are the PCM
model's "solvation effect" and make contributions to the
energy and to the gradient of the solute.
Typically the cavity is defined as a union of atomic
spheres, which should be roughly 1.2 times the atomic van
der Waals radii. A technical difficulty caused by the
penetration of the solute's charge density outside this
cavity is dealt with by a renormalization. The solvent is
characterized by its dielectric constant, surface tension,
size, density, and so on. Procedures are provided not only
for the computation of the electrostatic interaction of the
solute with the apparent surface charges, but also for the
cavitation energy, and for the dispersion and repulsion
contributions to the solvation free energy.
Methodology for solving the Poisson equation to obtain
the "apparent surface charges" has progressed from D-PCM to
IEF-PCM to C-PCM over time, with the latter preferred.
Iterative solvers require far less computer resources than
direct solvers. Advancements have also been made in
schemes to divide the surface cavity into tiny tesserae.
As of fall 2008, the FIXPVA tessellation, which has smooth
switching functions for tesserae near sphere boundaries,
together with iterative C-PCM, gives very satisfactory
geometry optimizations for molecules of 100 atoms. The
FIXPVA tessellation was extended to work for cavitation
(ICAV), dispersion (IDP), and repulsion (IREP) options in
fall 2009, and dispersion/repulsion (IDISP) in spring 2010.
Other procedures remain, and make the input seem complex,
but their use is discouraged. Thus
$PCM SOLVNT=WATER $END
chooses iterative C-PCM (IEF=-10) and FIXPVA tessellation
(METHOD=4 in $TESCAV) to do basic electrostatics in an
accurate fashion.
The main input group is $PCM, with $PCMCAV providing
auxiliary cavity information. If any of the optional
energy computations are requested in $PCM, the additional
input groups $IEFPCM, $NEWCAV, $DISBS, or $DISREP may be
required.
It is useful to summarize the various cavities used by
PCM, since as many as three cavities may be used:
the basic cavity for electrostatics,
cavity for cavitation energy, if ICAV=1,
cavity for dispersion/repulsion, if IDISP=1.
The first and second share the same radii (see RADII in
$PCMCAV), which are scaled by ALPHA=1.2 for electrostatics,
but are used unscaled for cavitation. The dispersion
cavity is defined in $DISREP with a separate set of atomic
radii, and even solvent molecule radii! Only the
electrostatics cavity can use any of the GEPOL-GB, GEPOL-AS
or recommended FIXPVA tessellation, while the other two use
only the original GEPOL-GB.
Radii are an important part of the PCM parameterization.
Their values can have a significant impact on the quality
of the results. This is particularly true if the solute is
charged, and thus has a large electrostatic interaction
with the continuum. John Emsley's book "The Elements" is a
useful source of van der Waals and other radii.
PCM is at heart a means of treating the electrostatic
interactions between the solute's wavefunction and a
dielectric model for the bulk solvent. The former is
represented as an electron density from whatever quantum
mechanical treatment is used for the solute, and the latter
is a set of surface charges on the finite elements of the
cavity (tessellation). This leaves out other important
contributions to the solvation energy! These include the
energy needed to make a hole in the solvent (cavitation
energy), dispersion or repulsive interactions between the
solute and solvent, and in the SMD model (see below)
solvent structure changes such as would occur in the first
solvation shell. Some empirical formulae for such "CDS"
corrections are provided as keywords ICAV, IDISP, IREP/IDP,
which may not work with all wavefunctions, and may not be
compatible with gradients.
The SMD model gives an alternative set of such "CDS"
corrections, which are compatible with nuclear gradients:
see SMD=.TRUE. in $PCM. A more detailed description of SMD
is given in the paper cited below. The SMD solvent
parameters is described (as of 2010) in
http://comp.chem.umn.edu/solvation/mnsddb.pdf
This gives numerical parameters, all built into GAMESS, as
the various SOLX values, for SOLVNT=
ACETACID CLPROPAN PHOPH EGME E2PENTEN
ACETONE OCLTOLUE DPROAMIN MEACETAT PENTACET
ACETNTRL M-CRESOL DODECAN MEBNZATE PENTAMIN
ACETPHEN O-CRESOL MEG MEBUTATE PFB
ANILINE CYCHEXAN ETSH MEFORMAT BENZALCL
ANISOLE CYCHEXON ETHANOL MIBK PROPANAL
BENZALDH CYCPENTN ETOAC MEPROPYL PROPACID
BENZENE CYCPNTOL ETOME ISOBUTOL PROPANOL
BENZNTRL CYCPNTON EB TERBUTOL PROPNOL2
BENZYLCL DECLNCIS PHENETOL NMEANILN PROPNTRL
BRISOBUT DECLNTRA C6H5F MECYCHEX PROPENOL
BRBENZEN DECLNMIX FOCTANE NMFMIXTR PROPACET
BRETHANE DECANE FORMAMID ISOHEXAN PROPAMIN
BROMFORM DECANOL FORMACID MEPYRID2 PYRIDINE
BROCTANE EDB12 HEPTANE MEPYRID3 C2CL4
BRPENTAN DIBRMETN HEPTANOL MEPYRID4 THF
BRPROPA2 BUTYLETH HEPTNON2 C6H5NO2 SULFOLAN
BRPROPAN ODICLBNZ HEPTNON4 C2H5NO2 TETRALIN
BUTANAL EDC12 HEXADECN CH3NO2 THIOPHEN
BUTACID C12DCE HEXANE NTRPROP1 PHSH
BUTANOL T12DCE HEXNACID NTRPROP2 TOLUENE
BUTANOL2 DCM HEXANOL ONTRTOLU TBP
BUTANONE ETHER HEXANON2 NONANE TCA111
BUTANTRL ET2S HEXENE NONANOL TCA112
BUTILE DIETAMIN HEXYNE NONANONE TCE
NBA MI C6H5I OCTANE ET3N
NBUTBENZ DIPE IOBUTANE OCTANOL TFE222
SBUTBENZ DMDS C2H5I OCTANON2 TMBEN124
TBUTBENZ DMSO IOHEXDEC PENTDECN ISOCTANE
CS2 DMA CH3I PENTANAL UNDECANE
CARBNTET CISDMCHX IOPENTAN NPENTANE M-XYLENE
CLBENZEN DMF IOPROPAN PENTACID O-XYLENE
SECBUTCL DMEPEN24 CUMENE PENTANOL P-XYLENE
CHCL3 DMEPYR24 P-CYMENE PENTNON2 XYLENEMX
CLHEXANE DMEPYR26 MESITYLN PENTNON3
CLPENTAN DIOXANE METHANOL PENTENE
and provides a translation table to full chemical names, if
you can't guess from the input choices given above. The
translations can also be found in the source code. Two
important things to note about SMD are:
a) the atomic radii are changed, so although the
algorithms for electrostatics are those of standard PCM,
the numerical results for the electrostatics do change.
b) SMD's parameterization was developed for IEF-PCM
using GEPOL tessellation with a fine grid: IEF=-3 and
MTHALL=2, NTSALL=240. However it is considered acceptable
to use SMD's parameters, unchanged, with C-PCM and with the
FIXPVA tessellation, at default coarseness. Hence, input
such as
$pcm solvnt=dmso smd=.true. $end
is enough to carry out a SMD-style C-PCM treatment in DMSO.
c) The CDS correction involves cavitation, dispersion,
and as a collective "solvent structure contribution"
estimates for partial hydrogen bonding, repulsion, and
deviation of the dielectric constant from its bulk value.
d) See also SMVLE in the more sophisticated SS(V)PE
continuum model's description.
Solvation of course affects the non-linear optical
properties of molecules. The PCM implementation extends
RUNTYP=TDHF to include solvent effects. Both static and
frequency dependent hyperpolarizabilities can be found.
Besides the standard PCM electrostatic contribution, the
IREP and IDP keywords can be used to determine the effects
of repulsion and dispersion on the polarizabilities.
The implementation of the PCM model in GAMESS has
received considerable attention from Hui Li and Jan Jensen
at the University of Iowa, Iowa State University, and
University of Nebraska. This includes new iterative
techniques to solving the surface charge problem, new
tessellations that provide for numerically stable nuclear
gradients, the implementation of C-PCM equations, the
extension of PCM to all SCFTYPs and TDDFT, development of
an interface with the EFP model (quo vadis), and
heterogenous dielectric. Dmitri Fedorov at AIST has
interfaced PCM to the FMO method (quo vadis), and reduced
storage requirements.
Due to its sophistication, users of the PCM model are
strongly encouraged to read the primary literature:
Of particular relevance to PCM in GAMESS:
1) "Continuum solvation of large molecules described by
QM/MM: a semi-iterative implementation of the PCM/EFP
interface"
H.Li, C.S.Pomelli, J.H.Jensen
Theoret.Chim.Acta 109, 71-84(2003)
2) "Improving the efficiency and convergence of geometry
optimization with the polarizable continuum model: new
energy gradients and molecular surface tessellation"
H.Li, J.H.Jensen J.Comput.Chem. 25, 1449-1462(2004)
3) "The polarizable continuum model interfaced with the
Fragment Molecular Orbital method"
D.G.Fedorov, K.Kitaura, H.Li, J.H.Jensen, M.S.Gordon
J.Comput.Chem. 27, 976-985(2006)
4) "Energy gradients in combined Fragment Molecular Orbital
and Polarizable Continuum Model (FMO/PCM)"
H.Li, D.G.Fedorov, T.Nagata, K.Kitaura, J.H.Jensen,
M.S.Gordon J.Comput.Chem. 31, 778-790(2010)
5) "Continuous and smooth potential energy surface for
conductor-like screening solvation model using fixed points
with variable area"
P.Su, H.Li J.Chem.Phys. 130, 074109/1-13(2009)
6) "Heterogenous conductorlike solvation model"
D.Si, H.Li J.Chem.Phys. 131, 044123/1-8(2009)
7) "Quantum mechanical/molecular mechanical/continuum style
solvation model: linear response theory, variational
treatment, and nuclear gradients"
H.Li J.Chem.Phys. 131, 184103/1-8(2009)
8) "Smooth potential energy surface for cavitation,
dispersion, and repulsion free energies in polarizable
continuum model"
Y.Wang, H.Li J.Chem.Phys. 131, 206101/1-2(2009)
9) "Excited state geometry of photoactive yellow protein
chromophore: a combined conductorlike polarizable continuum
model and time-dependent density functional study"
Y.Wang, H.Li J.Chem.Phys. 133, 034108/1-11(2010)
Paper number 7 is about the treatment of QM systems with
the solvation models EFP and/or C-PCM.
SMD and its CDS cavitation/dispersion/solvent structure
corrections are described in
"Universal solution model based on solute electron density
and on a continuum model of the solvent defined by the bulk
dielectric constant and atomic surface tensions"
A.V.Marenich, C.J.Cramer, D.G.Truhlar
J.Phys.Chem.B 113, 6378-6396(2009)
General papers on PCM:
10) S.Miertus, E.Scrocco, J.Tomasi
Chem.Phys. 55, 117-129(1981)
11) J.Tomasi, M.Persico Chem.Rev. 94, 2027-2094(1994)
12) R.Cammi, J.Tomasi J.Comput.Chem. 16, 1449-1458(1995)
13) J.Tomasi, B.Mennucci, R.Cammi
Chem.Rev. 105, 2999-3093(2005)
The GEPOL-GB method for cavity construction:
14) J.L.Pascual-Ahuir, E.Silla, J.Tomasi, R.Bonaccorsi
J.Comput.Chem. 8, 778-787(1987)
Charge renormalization (see also ref. 12):
15) B.Mennucci, J.Tomasi J.Chem.Phys. 106, 5151-5158(1997)
Derivatives with respect to nuclear coordinates:
(energy gradient and hessian) See also paper 2 and 3.
16) R.Cammi, J.Tomasi J.Chem.Phys. 100, 7495-7502(1994)
17) R.Cammi, J.Tomasi J.Chem.Phys. 101, 3888-3897(1995)
18) M.Cossi, B.Mennucci, R.Cammi
J.Comput.Chem. 17, 57-73(1996)
Derivatives with respect to applied electric fields:
(polarizabilities and hyperpolarizabilities)
19) R.Cammi, J.Tomasi
Int.J.Quantum Chem. Symp. 29, 465-474(1995)
20) R.Cammi, M.Cossi, J.Tomasi
J.Chem.Phys. 104, 4611-4620(1996)
21) R.Cammi, M.Cossi, B.Mennucci, J.Tomasi
J.Chem.Phys. 105, 10556-10564(1996)
22) B. Mennucci, C. Amovilli, J. Tomasi
Chem.Phys.Lett. 286, 221-225(1998)
Cavitation energy:
23) R.A.Pierotti Chem.Rev. 76, 717-726(1976)
24) J.Langlet, P.Claverie, J.Caillet, A.Pullman
J.Phys.Chem. 92, 1617-1631(1988)
Dispersion and repulsion energies:
25) F.Floris, J.Tomasi J.Comput.Chem. 10, 616-627(1989)
26) C.Amovilli, B.Mennucci
J.Phys.Chem.B 101, 1051-1057(1997)
Integral Equation Formalism PCM. The first of these
deals with anisotropies, the last 2 with nuclear gradients.
27) E.Cances, B.Mennucci, J.Tomasi
J.Chem.Phys. 107, 3032-3041(1997)
28) B.Mennucci, E.Cances, J.Tomasi
J.Phys.Chem.B 101, 10506-17(1997)
29) B.Mennucci, R.Cammi, J.Tomasi
J.Chem.Phys. 109, 2798-2807(1998)
30) J.Tomasi, B.Mennucci, E.Cances
J.Mol.Struct.(THEOCHEM) 464, 211-226(1999)
31) E.Cances, B.Mennucci J.Chem.Phys. 109, 249-259(1998)
32) E.Cances, B.Mennucci, J.Tomasi
J.Chem.Phys. 109, 260-266(1998)
Conductor PCM (C-PCM):
33) V.Barone, M.Cossi J.Phys.Chem.A 102, 1995-2001(1998)
34) M.Cossi, N.Rega, G.Scalmani, V.Barone
J.Comput.Chem. 24, 669-681(2003)
C-PCM with TD-DFT:
35) M.Cossi, V.Barone J.Chem.Phys. 115, 4708-4717(2001)
See also paper #8 above for the coding in GAMESS.
At the present time, the PCM model in GAMESS has the
following limitations:
a) Any SCFTYP may be used (RHF to MCSCF). MP2 or DFT
may be used with any of the RHF, UHF, and ROHF
gradient programs. Closed shell TD-DFT excited
state gradients may also be used.
CI and Coupled Cluster programs are not available.
b) semi-empirical methods may not be used.
c) the only other solvent method that may be used at
used with PCM is the EFP model.
d) point group symmetry is switched off internally
during PCM.
e) The PCM model runs in parallel for IEF=3, -3, 10,
or -10 and for all 5 wavefunctions (energy or
gradient), but not for RUNTYP=TDHF jobs.
f) D-PCM stores electric field integrals at normals to
the surface elements on disk.
IEF-PCM and C-PCM using the explicit solver (+3 and
+10) store electric potential integrals at normals
to the surface on disk.
This is true even for direct AO integral runs, and
the file sizes may be considerable (basis set size
squared times the number of tesserae).
IEF-PCM and C-PCM with the iterative solvers do not
store the potential integrals, when IDIRCT=1 in the
$PCMITR group (this is the default)
g) nuclear derivatives are limited to gradients,
although theory for hessians is given in paper 17.
* * *
The only PCM method prior to Oct. 2000 was D-PCM, which
can be recovered by selecting IEF=0 and ICOMP=2 in $PCM.
The default PCM method between Oct. 2000 and May 2004 was
IEF-PCM, recoverable by IEF=-3 (but 3 for non-gradient
runs) and ICOMP=0. As of May 2004, the default PCM method
was changed to C-PCM (IEF=-10, ICOMP=0). The extension of
PCM to all SCFTYPs as of May 2004 involved a correction to
the MCSCF PCM operator, so that it would reproduce RHF
results when run on one determinant, meaning that it is
impossible to reproduce prior MCSCF PCM calculations.
The cavity definition was GEPOL-GB (MTHALL=1 in $TESCAV)
prior to May 2004, GEPOL-AS (MTHALL=2) from then until
September 2008, and FIXPVA (MTHALL=4) to the present time.
The option for generation of 'extra spheres' (RET in $PCM)
was changed from 0.2 to 100.0, to suppress these, in June
2003.
* * *
In general, use of PCM electrostatics is very simple, as
may be seen from exam31.inp supplied with the program.
The calculation shown next illustrates the use of some
of the older PCM options. Since methane is non-polar, its
internal energy change and the direct PCM electrostatic
interaction is smaller than the cavitation, repulsion, and
dispersion corrections. Note that the use of ICAV, IREP,
and IDP are currently incompatible with gradients, so a
reasonable calculation sequence might be to perform the
geometry optimization with PCM electrostatics turned on,
then perform an additional calculation to include the other
solvent effects, adding extra functions to improve the
dispersion correction.
! calculation of CH4 (metano), in PCM water.
!
! This input reproduces the data in Table 2, line 6, of
! C.Amovilli, B.Mennucci J.Phys.Chem.B 101, 1051-7(1997)
! To do this, we must use many original PCM options.
!
! The gas phase FINAL energy is -40.2075980292
! The FINAL energy in PCM water is -40.2048210283
! (lit.)
! FREE ENERGY IN SOLVENT = -25234.89 KCAL/MOL
! INTERNAL ENERGY IN SOLVENT = -25230.64 KCAL/MOL
! DELTA INTERNAL ENERGY = .01 KCAL/MOL ( 0.0)
! ELECTROSTATIC INTERACTION = -.22 KCAL/MOL (-0.2)
! PIEROTTI CAVITATION ENERGY = 5.98 KCAL/MOL ( 6.0)
! DISPERSION FREE ENERGY = -6.00 KCAL/MOL (-6.0)
! REPULSION FREE ENERGY = 1.98 KCAL/MOL ( 2.0)
! TOTAL INTERACTION = 1.73 KCAL/MOL ( 1.8)
! TOTAL FREE ENERGY IN SOLVENT= -25228.91 KCAL/MOL
!
$contrl scftyp=rhf runtyp=energy $end
$guess guess=huckel $end
$system mwords=2 $end
! the "W1 basis" input here exactly matches HONDO's DZP
$DATA
CH4...gas phase geometry...in PCM water
Td
Carbon 6.
DZV
D 1 ; 1 0.75 1.0
Hydrogen 1. 0.6258579976 0.6258579976 0.6258579976
DZV 0 1.20 1.15 ! inner and outer scale factors
P 1 ; 1 1.00 1.0
$END
! The reference cited used a value for H2O's solvent
! radius that differs from the built in value (RSOLV).
! The IEF, ICOMP, MTHALL, and RET keywords are set to
! duplicate the original code's published results,
! namely D-PCM and GEPOL-GB. This run doesn't put in
! any "extra spheres" but we try that option (RET)
! like it originally would have.
$PCM SOLVNT=WATER RSOLV=1.35 RET=0.2
IEF=0 ICOMP=2 IDISP=0 IREP=1 IDP=1 ICAV=1 $end
$TESCAV MTHALL=1 $END
$NEWCAV IPTYPE=2 ITSNUM=540 $END
! dispersion "W2 basis" uses exponents which are
! 1/3 of smallest exponent in "W1 basis" of $DATA.
$DISBS NADD=11 NKTYP(1)=0,1,2, 0,1, 0,1, 0,1, 0,1
XYZE(1)=0.0,0.0,0.0, 0.0511
0.0,0.0,0.0, 0.0382
0.0,0.0,0.0, 0.25
1.1817023, 1.1817023, 1.1817023, 0.05435467
1.1817023, 1.1817023, 1.1817023, 0.33333333
-1.1817023, 1.1817023,-1.1817023, 0.05435467
-1.1817023, 1.1817023,-1.1817023, 0.33333333
1.1817023,-1.1817023,-1.1817023, 0.05435467
1.1817023,-1.1817023,-1.1817023, 0.33333333
-1.1817023,-1.1817023, 1.1817023, 0.05435467
-1.1817023,-1.1817023, 1.1817023, 0.33333333 $end
SVPE and SS(V)PE.
The Surface Volume Polarization for Electrostatics
(SVPE), and an approximation to SVPE called the Surface and
Simulation of Volume Polarization for Electrostatics
(SS(V)PE) are continuum solvation models. Compared to
other continuum models, SVPE and SS(V)PE pay careful
attention to the problems of escaped charge, the shape of
the surface cavity, and to integration of the Poisson
equation for surface charges.
The original references for what is now called the
SVPE (surface and volume polarization for electrostatics)
method are the theory paper:
"Charge penetration in Dielectric Models of Solvation"
D.M.Chipman, J.Chem.Phys. 106, 10194-10206 (1997)
and two implementation papers:
"Volume Polarization in Reaction Field Theory"
C.-G.Zhan, J.Bentley, D.M.Chipman
J.Chem.Phys. 108, 177-192 (1998)
"New Formulation and Implementation for Volume
Polarization in Dielectric Continuum Theory"
D.M.Chipman, J.Chem.Phys. 124, 224111-1/10 (2006)
which should be cited in any publications that utilize the
SVPE code.
There are two options to include with SS(V)PE or SVPE
additional models that describe short-range solute-solvent
interactions to achieve a more complete description of
solvation energies. Both have their keywords merged into
the $SVP input group for convenience. One option to
include short-range interactions is CMIRS (Composite
Method for Implicit Representation of Solvent)
that combines the SS(V)PE dielectric continuum model with
the DEFESR (Dispersion, Exchange, and Field-Extremum Short-
Range) model.
A complete account of the original version
labeled CMIRS1.0 with application to hydration energies is
given in:
Hydration Energy from a Composite Method for Implicit
Representation of Solvent
A.Pomogaeva, D.M.Chipman
J.Chem.Theory Comput. 10, 211-219(2014)
which should be referenced in any publication that uses the
model. Applications to DMSO and acetonitrile solvents are
reported in
Composite Method for Implicit Representation of
Solvent in Diemthyl Sulfoxide and Acetonitrile
A.Pomogaeva, D.M. Chipman
J.Phys.Chem.A 119, 5173-5180(2015)
References to earlier works that develop the individual
components of the DEFESR parts of the full CMIRS recipe are:
Modeling short-range contributions to hydration
energies with minimal parameterization
A.Pomogaeva, D.W.Thompson, and D.M.Chipman
Chem.Phys.Lett. 511, 161-165(2011).
Field-Extremum Model for Short-Range Contributions
to Hydration Free Energy
A.Pomogaeva and D.M.Chipman
J.Chem.Theory Comput. 7, 3952-3960(2011).
New Implicit Solvation Models for Dispersion and
Exchange Energies
A.Pomogaeva and D.M.Chipman
J.Phys.Chem.A, 117, 5812-5820 (2013).
An error in the original CMIRS1.0 code for dispersion
energy has been reported in
Reparameterization of an Accurate, Few-Parameter
Implicit Solvation Model for Quantum Chemistry:
Composite Method for Implicit Representation of
Solvent, CMIRS v. 1.1
Zhi-Qiang You, John M. Herbert,
J.Chem.Theor.Comp. 12, 4338-4346 (2016)
This paper defines the CMIRS1.1 method giving revised
parameters for use with water, cyclohexane, benzene, DMSO,
and acetonitrile solvents obtained with the B3LYP/6-31+G*
method in conjunction with isodensity contours of 0.005
and 0.001 au. This error has been corrected in the current
Gamess code.
The other option to include short-range interactions is
the SMVLE (solvation model with volume and local
electrostatics) as described in
"Free energies of solvation with surface, volume, and
local electrostatic effects and atomic surface
tensions to represent the first solvation shell"
J.Liu, C.P.Kelly, A.C.Goren, A.V.Marenich,
C.J.Cramer, D.G.Truhlar, C.-G. Zhan
J.Chem.Theory Comput. 6, 1109-1117(2010).
Further information on the performance of SVPE and of
SS(V)PE can be found in:
"Comparison of Solvent Reaction Field Representations"
D.M.Chipman, Theor.Chem.Acc. 107, 80-89 (2002).
Details of the SS(V)PE convergence behavior and programming
strategy are in:
"Implementation of Solvent Reaction Fields for
Electronic Structure" D.M.Chipman, M.Dupuis,
Theor.Chem.Acc. 107, 90-102 (2002).
The SMVLE option (solvation model with volume and
local electrostatics) is described in
"Free energies of solvation with surface, volume, and
local electrostatic effects and atomic surface tensions to
represent the first solvation shell" J.Liu, C.P.Kelly,
A.C.Goren, A.V.Marenich, C.J.Cramer, D.G.Truhlar, C.-G.
Zhan J.Chem.Theory Comput. 6, 1109-1117(2010).
The SVPE and SS(V)PE models are like PCM and COSMO in
that they treat solvent as a continuum dielectric residing
outside a molecular-shaped cavity, determining the apparent
charges that represent the polarized dielectric by solving
Poisson's equation. The main difference between SVPE and
SS(V)PE is in treatment of volume polarization effects that
arise because of the tail of the electronic wave function
that penetrates outside the cavity, sometimes referred to
as the "escaped charge." SVPE treats volume polarization
effects explicitly by including apparent charges in the
volume outside the cavity as well as on the cavity surface.
With a sufficient number of grid points, SVPE can then
provide an exact treatment of charge penetration effects.
SS(V)PE, like PCM and COSMO, is an approximate treatment
that only uses apparent charges located on the cavity
surface. The SS(V)PE equation is particularly designed to
simulate as well as possible the influence of the missing
volume charges. For more information on the similarities
and differences of the SVPE and SS(V)PE models with other
continuum methods, see the paper "Comparison of Solvent
Reaction Field Representations" cited just above.
In addition, the cavity construction and Poisson
solver used in this implementation of SVPE and SS(V)PE also
receive careful numerical treatment. For example, the
cavity may be chosen to be an isodensity contour surface,
and the Lebedev grids for the Poisson solver can be chosen
very densely. The Lebedev grids used for surface
integration are taken from the Fortran translation by C.
van Wuellen of the original C programs developed by D.
Laikov. They were obtained from the CCL web site
www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-
Grids. A recent leading reference is V. I. Lebedev and D.
N. Laikov, Dokl. Math. 59, 477-481 (1999). All these grids
have octahedral symmetry and so are naturally adapted for
any solute having an Abelian point group. The larger and/or
the less spherical the solute may be, the more surface
points are needed to get satisfactory precision in the
results. Further experience will be required to develop
detailed recommendations for this parameter. Values as
small as 110 are usually sufficient for simple diatomics
and triatomics. The default value of 1202 has been found
adequate to obtain the energy to within 0.1 kcal/mol for
solutes the size of monosubstituted benzenes. The SVPE
method uses additional layers of points outside the cavity.
Typically just two layers are sufficient to converge the
direct volume polarization contribution to better than 0.1
kcal/mol.
The SVPE and SS(V)PE codes both report the amount of
solute charge penetrating outside the cavity as calculated
by Gauss' Law. The SVPE code additionally reports the same
quantity as alternatively calculated from the explicit
volume charges, and any substantial discrepancy between
these two determinations indicates that more volume
polarization layers should have been included for better
precision. The energy contribution from the outermost
volume polarization layer is also reported. If it is
significant then again more layers should have been
included. However, these tests are only diagnostic.
Passing them does not guarantee that enough layers are
included.
Analytic nuclear gradients are not yet available for
the SVPE or SS(V)PE energy, but numerical differentiation
will permit optimization of small solute molecules.
Wavefunctions may be any of the SCF type: RHF, UHF, ROHF,
GVB, and MCSCF, or the DFT analogs of some of these. In the
MCSCF implementation, no initial wavefunction is available
so the solvation code does not kick in until the second
iteration.
We close with a SVPE example. The gas phase energy,
obtained with no $SVP group, is -207.988975, and the run
just below gives the SVPE energy -208.006282. The free
energy of solvation, -10.860 kcal/mole, is the difference
of these, and is quoted at the right side of the 3rd line
from the bottom of Table 2 in the paper cited. The
"REACTION FIELD FREE ENERGY" for SVPE is -12.905 kcal/mole,
which is only part of the solvation free energy. There is
also a contribution due to the SCRF procedure polarizing
the wave function from its gas phase value, causing the
solute internal energy in dielectric to differ from that in
gas. Evaluating this latter contribution is what requires
the separate gas phase calculation. Changing the number of
layers (NVLPL) to zero produces the SS(V)PE approximation
to SVPE, E= -208.006208.
! SVPE solvation test...acetamide
! reproduce data in Table 2 of the paper on SVPE,
! D.M.Chipman J.Chem.Phys. 124, 224111/1-10(2006)
!
$contrl scftyp=rhf runtyp=energy $end
$system mwords=4 $end
$basis gbasis=n31 ngauss=6 ndfunc=1 npfunc=1 $end
$guess guess=moread norb=16 $end
$scf nconv=8 $end
$svp nvlpl=3 rhoiso=0.001 dielst=78.304 nptleb=1202 $end
$data
CH3CONH2 cgz geometry RHF/6-31G(d,p)
C1
C 6.0 1.361261 -0.309588 -0.000262
C 6.0 -0.079357 0.152773 -0.005665
H 1.0 1.602076 -0.751515 0.962042
H 1.0 1.537200 -1.056768 -0.767127
H 1.0 2.002415 0.542830 -0.168045
O 8.0 -0.387955 1.310027 0.002284
N 7.0 -1.002151 -0.840834 -0.011928
H 1.0 -1.961646 -0.589397 0.038911
H 1.0 -0.752774 -1.798630 0.035006
$end
gas phase vectors, E(RHF)= -207.9889751769
$VEC
1 1 1.18951670E-06 1.74015997E-05
...snipped...
$END
The example just above can be changed to a CMIRS run, by
changing NVLPL to zero and also requesting calculation of
the DEFESR contributions. The result should be a CMIRS1.1
total free energy of -208.013009. The input should
(1) in $CONTRL, specify DFTTYP=HFX, which requests a
Hartree-Fock calculation, but sets up DFT grids,
(2) in $DFT, request a very accurate grid
$dft method=grid nrad=96 nleb=1202 nrad0=96 nleb0=1202 $end
(3) in $SVP, invoke the keyword IDEF=1.
Adding the keyword EGAS=-207.9889751769 to any of
these examples allows reporting of the free energy of
solvation (-10.860 kcal/mol for SVPE, -10.814 for SS(V)PE,
-15.081 for CMIRS1.1). Note that this example invokes the
HF/6-31G** method for the wavefunction while inconsistently
utilizing CMIRS1.1 solvation parameters that were instead
optimized for use with the B3LYP/6-31+G* method. The latter
method gives much better agreement with experiment.
Conductor-like screening model (COSMO)
The COSMO (conductor-like screening model) represents a
different approach for carrying out polarized continuum
calculations. COSMO was originally developed by Andreas
Klamt, with extensions to ab initio computation in GAMESS
by Kim Baldridge.
In the COSMO method, the surrounding medium is modeled
as a conductor rather than as a dielectric in order to
establish the initial boundary conditions. The assumption
that the surrounding medium is well modeled as a conductor
simplifies the electrostatic computations and corrections
may be made a posteriori for dielectric behavior.
The original model of Klamt was introduced using a
molecular shaped cavity, which had open parts along the
crevices of intersecting atomic spheres. While having
considerable technical advantages, this approximation
causes artifacts in the context of the more generalized
theory, so the current method for cavity construction
includes a closure of the cavity to eliminate crevices or
pockets.
Two methodologies are implemented for treatment of the
outlying charge errors (OCE). The default is the well-
established double cavity procedure using a second, larger
cavity around the first one, and calculates OCE through the
difference between the potential on the inner and the outer
cavity. The second involves the calculation of distributed
multipoles up to hexadecapoles to represent the entire
charge distribution of the molecule within the cavity.
The COSMO model accounts only for the electrostatic
interactions between solvent and solute. Klamt has
proposed a novel statistical scheme to compute the full
solvation free energy for neutral solutes, COSMO-RS, which
is formulated for GAMESS by Peverati, Potier and Baldridge,
and is available as external plugin to the COSMOtherm
program by COSMOlogic GmbH&Co.
The iterative inclusion of non-electrostatic effects is
also possible right after a COSMO-RS calculation. The
DCOSMO-RS approach was implemented in GAMESS by Peverati,
Potier, and Baldridge, and more information is available on
Baldridge website at:
http://ocikbws.uzh.ch/gamess/
The simplicity of the COSMO model allows computation of
gradients, allowing optimization within the context of the
solvent. The method is programmed for RHF and UHF, all
corresponding kinds of DFT (including DFT-D), and the
corresponding MP2, energy and gradient.
Some references on the COSMO model are:
A.Klamt, G.Schuurman
J.Chem.Soc.Perkin Trans 2, 799-805(1993)
A.Klamt J.Phys.Chem. 99, 2224-2235(1995)
K.Baldridge, A.Klamt
J.Chem.Phys. 106, 6622-6633 (1997)
V.Jonas, K.Baldridge
J.Chem.Phys. 113, 7511-7518 (2000)
L.Gregerson, K.Baldridge
Helv.Chim.Acta 86, 4112-4132 (2003)
R.Peverati, Y.Potier, K.Baldridge
TO BE PUBLISHED SOON
Additional references on the COSMO-RS model, with
explanation of the methodology and program can be found:
A.Klamt, F.Eckert, W.Arlt
Annu.Rev.Chem.Biomol.Eng. 1, (2010)
The Effective Fragment Potential Method
The basic idea behind the effective fragment potential
(EFP) method is to replace the chemically inert part of a
system by EFPs, while performing a regular ab initio
calculation on the chemically active part. Here "inert"
means that no covalent bond breaking process occurs. This
"spectator region" consists of one or more "fragments",
which interact with the ab initio "active region" through
non-bonded interactions, and so of course these EFP
interactions affect the ab initio wavefunction. The EFP
particles can be closed shell or open shell (high spin
ROHF) based potentials. The "active region" can use nearly
every kind of wavefunction available in GAMESS.
A simple example of an active region might be a solute
molecule, with a surrounding spectator region of solvent
molecules represented by fragments. Each discrete solvent
molecule is represented by a single fragment potential, in
marked contrast to continuum models for solvation.
The quantum mechanical part of the system is entered in
the $DATA group, along with an appropriate basis. The EFPs
defining the fragments are input by means of a $EFRAG
group, and one or more $FRAGNAME groups describing each
fragment's EFP. These groups define non-bonded
interactions between the ab initio system and the
fragments, and also between the fragments. The former
interactions enter via one-electron operators in the ab
initio Hamiltonian, while the latter interactions are
treated by analytic functions. The only electrons
explicitly treated (with basis functions used to expand
occupied orbitals) are those in the active region, so there
are no new two electron terms. Thus the use of EFPs leads
to significant time savings, compared to full ab initio
calculations on the same system.
There are two types of EFP available in GAMESS, EFP1 and
EFP2. EFP1, the original method, employs a fitted
repulsive potential. EFP1 is primarily used to model water
molecules to study aqueous solvation effects, at the
RHF/DZP or DFT/DZP (specifically, B3LYP) levels, see
references 1-3 and 26, respectively. EFP2 is a more
general method that is applicable to any species, including
water, and its repulsive potential is obtained from first
principles. EFP2 has been extended to include other
effects as well, such as charge transfer and dispersion.
EFP2 forms the basis of the covalent EFP method described
below for modeling enzymes, see reference 14.
Parallelization of the EFP1 and EFP2 models is described
in reference 32.
MD simulations with EFP are described in reference 31.
The ab initio/EFP1, or pure EFP system can be wrapped in
a Polarizable Continuum Model, see references 23, 43, and
50.
terms in an EFP
The non-bonded interactions currently implemented are:
1) Coulomb interaction. The charge distribution of the
fragments is represented by an arbitrary number of charges,
dipoles, quadrupoles, and octopoles, which interact with
the ab initio hamiltonian as well as with multipoles on
other fragments (see reference 2 and 18). It is possible
to use a screening term that accounts for the charge
penetration (reference 17 and 42). This screening term is
automatically included for EFP1. Typically the multipole
expansion points are located on atomic nuclei and at bond
midpoints.
2) Dipole polarizability. An arbitrary number of dipole
polarizability tensors can be used to calculate the induced
dipole on a fragment due to the electric field of the ab
initio system as well as all the other fragments. These
induced dipoles interact with the ab initio system as well
as the other EFPs, in turn changing their electric fields.
All induced dipoles are therefore iterated to self-
consistency. Typically the polarizability tensors are
located at the centroid of charge of each localized orbital
of a fragment. See reference 41.
3) Repulsive potential. Two different forms are used in
EFP1: one for ab initio-EFP repulsion and one for EFP-EFP
repulsion. The form of the potentials is empirical, and
consists of distributed Gaussian or exponential functions,
respectively. The primary contribution to the repulsion is
the quantum mechanical exchange repulsion, but the fitting
technique used to develop this term also includes the
effects of charge transfer. Typically these fitted
potentials are located on each atomic nucleus within the
fragment (see reference 3). In EFP2, polarization energies
can also be augmented by screening terms, analogous to the
electrostatic screening, to prevent "polarization collapse"
(MS in preparation)
For EFP2, the third term is divided into separate analytic
formulae for different physical interactions:
a) exchange repulsion
b) dispersion
c) charge transfer
A summary of EFP2, and its contrast to EFP1 can be found in
reference 18 and 44. The repulsive potential for EFP2 is
based on an overlap expansion using localized molecular
orbitals, as described in references 5, 6, and 9.
Dispersion energy is described in reference 34, and charge
transfer in reference 39 (which supercedes reference 22's
formulae).
EFP2 potentials have no fitted parameters, and can be
automatically generated during a RUNTYP=MAKEFP job, as
described below.
constructing an EFP1
RUNTYP=MOROKUMA assists in the decomposition of inter-
molecular interaction energies into electrostatic,
polarization, charge transfer, and exchange repulsion
contributions. This is very useful in developing EFPs
since potential problems can be attributed to a particular
term by comparison to these energy components for a
particular system.
A molecular multipole expansion can be obtained using
$ELMOM. A distributed multipole expansion can be obtained
by either a Mulliken-like partitioning of the density
(using $STONE) or by using localized molecular orbitals
($LOCAL: DIPDCM and QADDCM). The dipole polarizability
tensor can be obtained during a Hessian run ($CPHF), and a
distributed LMO polarizability expression is also available
($LOCAL: POLDCM).
In EFP1, the repulsive potential is derived by fitting
the difference between ab initio computed intermolecular
interaction energies, and the form used for Coulomb and
polarizability interactions. This difference is obtained
at a large number of different interaction geometries, and
is then fitted. Thus, the repulsive term is implicitly a
function of the choices made in representing the Coulomb
and polarizability terms. Note that GAMESS currently does
not provide a way to obtain these EFP1 repulsive potential.
Since a user cannot generate all of the EFP1 terms
necessary to define a new $FRAGNAME group using GAMESS, in
practice the usage of EFP1 is limited to the internally
stored H2ORHF or H2ODFT potentials mentioned below.
constructing an EFP2
As noted above, the repulsive potential for EFP2 is
derived from a localized orbital overlap expansion. It is
generally recommended that one use at least a double zeta
plus diffuse plus polarization basis set, e.g. 6-31++G(d,p)
to generate the EFP2 repulsive potential. However, it has
been observed that 6-31G(d) works reasonably well due to a
fortuitous cancellation of errors. The EFP2 potential for
any molecule can be generated as follows:
(a) Choose a basis set and geometry for the molecule of
interest. The geometry is ordinarily optimized at your
choice of Hartree-Fock/MP2/CCSD(T), with your chosen basis
set, but this is not a requirement. It is good to recall,
however, that EFP internal geometries are fixed, so it is
important to give some thought to the chosen geometry.
(b) Perform a RUNTYP=MAKEFP run for the chosen molecule
using the chosen geometry in $DATA and the chosen basis set
in $BASIS. This will generate the entire EFP2 potential in
the run's .efp file. The only user-defined variable that
must be filled in is changing the FRAGNAME's group name, to
$C2H5OH or $DMSO, etc. This step can use RHF or ROHF to
describe the electronic structure of the system.
(c) Transfer the entire fragment potential for the molecule
to any input file in which this fragment is to be used.
Since the internal geometry of an EFP is fixed, one need
only specify the first three atoms of any fragment in order
to position them in $EFRAG. Coordinates of any other atoms
in the rigid fragment will be automatically determined by
the program.
If the EFP contains less than three atoms, you can still
generate a fragment potential. After a normal MAKEFP run,
add dummy atoms (e.g. in the X and/or Y directions) with
zero nuclear charges, and add corresponding dummy bond
midpoints too. Carefully insert zero entries in the
multipole sections, and in the electrostatic screening
sections, for each such dummy point, but don't add data to
any other kind of EFP term such as polarizability. This
trick gives the necessary 3 points for use in $EFRAG groups
to specify "rotational" positions of fragments.
current limitations
1. For EFP1, the energy and energy gradient are programmed,
which permits RUNTYP=ENERGY, GRADIENT, and numerical
HESSIAN. The necessary programing to use the EFP gradients
to move on the potential surface are programmed for
RUNTYP=OPTIMIZE, SADPOINT, IRC, and VSCF, but the other
gradient based potential surface explorations such as DRC
are not yet available. Finally, RUNTYP=PROP is also
permissible.
For EFP2, the gradient terms for ab initio-EFP interactions
have not yet been coded, so geometry optimizations are only
sensible for a COORD=FRAGONLY run; that is, a run in which
only EFP2 fragments are present.
2. The ab initio part of the system must be treated with
RHF, ROHF, UHF, the open shell SCF wavefunctions permitted
by the GVB code, or MCSCF. DFT analogs of RHF, ROHF, and
UHF may also be used. Correlated methods such as MP2 and
CI should not be used.
3. EFPs can move relative to the ab initio system and
relative to each other, but the internal structure of an
EFP is frozen.
4. The boundary between the ab initio system and EFP1's
must not be placed across a chemical bond. However, see
the discussion below regarding covalent bonds.
5. Calculations must be done in C1 symmetry at present.
6. Reorientation of the fragments and ab initio system is
not well coordinated. If you are giving Cartesian
coordinates for the fragments (COORD=CART in $EFRAG), be
sure to use $CONTRL's COORD=UNIQUE option so that the ab
initio molecule is not reoriented.
7. If you need IR intensities, you have to use NVIB=2. The
potential surface is usually very soft for EFP motions, and
double differenced Hessians should usually be obtained.
practical hints for using EFPs
At the present time, we have only two internally stored
EFP potentials suitable for general use. These model
water, using the fragment name H2ORHF or H2ODFT. The
H2ORHF numerical parameters are improved values over the
values which were presented and used in reference 2, and
they also include the improved EFP-EFP repulsive term
defined in reference 3. The H2ORHF water EFP was derived
from RHF/DH(d,p) computations on the water dimer system.
When you use it, therefore, the ab initio part of your
system should be treated at the SCF level, using a basis
set of the same quality (ideally DH(d,p), but probably
other DZP sets such as 6-31G(d,p) will give good results as
well). Use of better basis sets than DZP with this water
EFP has not been tested. Similarly, H2ODFT was developed
using B3LYP/DZP water wavefunctions, so this should be used
(rather than H2ORHF) if you are using DFT to treat the
solute. Since H2ODFT water parameters are obtained from a
correlated calculation, they can also be used when the
solute is treated by MP2.
As noted, effective fragments have frozen internal
geometries, and therefore only translate and rotate with
respect to the ab initio region. An EFP's frozen
coordinates are positioned to the desired location(s) in
$EFRAG as follows:
a) the corresponding points are found in $FRAGNAME.
b) Point -1- in $EFRAG and its FRAGNAME equivalent are
made to coincide.
c) The vector connecting -1- and -2- is aligned with
the corresponding vector connecting FRAGNAME points.
d) The plane defined by -1-, -2-, and -3- is made to
coincide with the corresponding FRAGNAME plane.
Therefore the 3 points in $EFRAG define only the relative
position of the EFP, and not its internal structure. So, if
the "internal structure" given by points in $EFRAG differs
from the true values in $FRAGNAME, then the order in which
the points are given in $EFRAG can affect the positioning
of the fragment. It may be easier to input water EFPs if
you use the Z-matrix style to define them, because then you
can ensure you use the actual frozen geometry in your
$EFRAG. Note that the H2ORHF EFP uses the frozen geometry
r(OH)=0.9438636, a(HOH)=106.70327, and the names of its 3
fragment points are ZO1, ZH2, ZH3.
* * *
Building a large cluster of EFP particles by hand can be
tedious. The RUNTYP=GLOBOP program described below has an
option for constructing dense clusters. The method tries
to place particles near the origin, but not colliding with
other EFP particles already placed there, so that the
clusters grow outwards from the center. Here are some
ideas:
a) place 100 water molecules, all with the same coords
in $EFRAG. This will build up a droplet of water
with particles close together, but not on top of
each other, with various orientations.
b) place 16 waters (same coords, all first) followed by
16 methanols (also sharing their same coords, after
all waters). A 50-50 mixture of 32 molecules will
be created, if you choose the default of picking the
particles randomly from the initial list of 32.
c) to solvate a solute, add the solute in the $DATA
group at or near the origin. Add the solvent
molecules near by (same coords is ok), and run
the globop run with RNDINI as demonstrated below.
(optional, add MCTYP=3 to $GLOBOP input)
Example, allowing the random cluster to have 20 geometry
optimization steps:
$contrl runtyp=globop coord=fragonly $end
$globop rndini=.true. riord=rand mcmin=.true.
mctyp=4 nblock=0 $end
$statpt nstep=20 $end
$efrag
coord=cart
FRAGNAME=WATER
O1 -2.8091763203009 -2.1942725073400 -0.2722207394107
H2 -2.3676165499399 -1.6856118830379 -0.9334073942601
H3 -2.1441965467625 -2.5006167998896 0.3234583094693
...repeat this 15 more times...
FRAGNAME=MeOH
O1 4.9515153249 .4286994611 1.3368662306
H2 5.3392575544 .1717424606 3.0555957053
C3 6.2191743799 2.5592349960 .4064662379
H4 5.7024200977 2.7548960076 -1.5604873643
H5 5.6658856694 4.2696553371 1.4008542042
H6 8.2588049857 2.3458272252 .5282762681
...repeat 15 more times...
$end
$water
...give a full EFP2 potential for water...
$end
$meoh
...give a full EFP2 potential for methanol...
$end
Note that the random cluster generation now proceeds into a
full Monte Carlo simulation.
* * *
The translations and rotations of EFPs with respect to
the ab initio system and one another are automatically
quite soft degrees of freedom. After all, the EFP model is
meant to handle weak interactions! Therefore the
satisfactory location of structures on these flat surfaces
will require use of a tight convergence on the gradient:
OPTTOL=0.00001 in the $STATPT group.
The effect of a bulk continuum surrounding the solute
plus EFP waters can be obtained by using the PCM model, see
reference 23 and 43. To do this, simply add a $PCM group
to your input, in addition to the $EFRAG. The simultaneous
use of EFP and PCM allows for gradients, so geometry
optimization can be performed.
global optimization
If there are a large number of particles to move (EFP
and/or FMO and/or atom groups), it is difficult to locate
the lowest energy structures by hand. Typically these are
numerous, and one would like to have a number of them, not
just the very lowest energy. The RUNTYP of GLOBOP contains
a Monte Carlo procedure to generate a random set of
starting structures to look for those with the lowest
energy at a single temperature. If desired, a simulated
annealing protocol to cool the temperature may be used.
These two procedures may be combined with a local minimum
search, at some or all of the randomly generated
structures. The local minimum search is controlled by the
usual geometry optimizer, namely $STATPT input, and thus
permits the optimization of any ab initio atoms.
The Monte Carlo procedure by default uses a Metropolis
algorithm to move just one of the fragments. The method of
Parks to move all fragments simultaneously is also allowed.
The present program was used to optimize the structure
of water clusters. Let us consider the case of the twelve
water cluster, for which the following ten structures were
published by Day, Pachter, Gordon, and Merrill:
1. (D2d)2 -0.170209 6. (D2d)(C2) -0.167796
2. (D2d)(S4) -0.169933 7. S6 -0.167761
3. (S4)2 -0.169724 8. cage b -0.167307
4. D3 -0.168289 9. cage a -0.167284
5. (C1c)(Cs) -0.167930 10. (C1c)(C1c) -0.167261
A test input using Metropolis style Monte Carlo to examine
300 geometries at each temperature value, using simulated
annealing cooling from 200 to 50 degrees, and with local
minimization every 10 structures was run ten times. Each
run sampled about 7000 geometries. One simulation found
structure 2, while two of the runs found structure 3. The
other seven runs located structures with energy values in
the range -0.163 to -0.164. In all cases the runs began
with the same initial geometry, but produced different
results due to the random number generation used in the
Monte Carlo. Clearly one must try a lot of simulations to
be confident about having found most of the low energy
structures. In particular, it is good to try more than one
initial structure, unlike what was done in this test.
Ab initio atoms can be addressed using FMO, either in
multiple fragments, or perhaps a single large fragment.
Alternatively, ab initio atoms can be put into groups and
used directly in globop, which for small systems has a
lower overhead than FMO. In the case of large molecules
separated into multiple fragments, the keywords NPRBND,
PRSEP, IBNDS, and INDEP are applicable. These specify the
atoms in each set of fragments or groups whose bond is cut
in the fragmentation process. The paired atoms are
constrained during the Monte Carlo procedure to ensure that
the bond is not spacially broken. In the case where a
fragment that is being translated or rotated is paired with
two or more fragments, the movement is repeated on all
attached fragments, after randomly choosing which pair is
the starting point. For example, given a molecule split
into five fragments such that:
A-B-C-D-E
where A,B,C,D,E are the fragments. If C is chosen for a
translation, either B or D will be randomly chosen to be
the starting pair. When B is chosen as the starting pair,
C, D, and E will all be translated by the same amount:
A-B--C-D-E
which maintains the relative position of C, D, and E.
Setting INDEP=1 will not propagate the translation:
A-B--CD-E
So that only C is moved.
The same approach is used for rotations. Since a small
translation or rotation can result in a significant change
in the total system, it is advised that case be taken when
using solvent molecules and to the size of boundary
conditions. If a propagated movement moves a fragment
outside the boundary, a warning will be printed and the
step will be discarded as a proximity alert. Also, the
pair binding is not implemented for RNDINI=.TRUE. To
initialize a set of solvent molecules around pair bonded
fragments, include the pair bonded fragments in IFXFMO.
The Metrpolis Monte Carlo procedure involves the movement
of groups that are internally rigid. To introduce some
internal flexibility for FMO and ab initio groups, a
secondary Monte Carlo search where the entire system is
held rigid while the atoms in one group are moved is
implemented. The secondary Monte Carlo occurs when a FMO
or ab initio group is translated and occurs for that group.
The lowest energy internal configuration for the secondary
Monte Carlo is used when evaluating the step of the primary
Monte Carlo search. The temperature at which the secondary
Monte Carlo is used in the case of simulated annealing is
set by SMTEMP and the number of steps in each secondary
search is given by NSMTP. To turn on this feature, set the
values of SMTEMP and NSMTP to non-zero values.
Monte Carlo references:
N.Metropolis, A.Rosenbluth, A.Teller
J.Chem.Phys. 21, 1087(1953).
G.T.Parks Nucl.Technol. 89, 233(1990).
Monte Carlo with local minimization:
Z.Li, H.A.Scheraga
Proc.Nat.Acad.Sci. USA 84, 6611(1987).
Simulated annealing reference:
S.Kirkpatrick, C.D.Gelatt, M.P.Vecci
Science 220, 671(1983).
The present program is described in reference 15. It is
patterned on the work of
D.J.Wales, M.P.Hodges Chem.Phys.Lett. 286, 65-72 (1998).
QM/MM across covalent bonds
Recent work by Visvaldas Kairys and Jan Jensen has made
it possible to extend the EFP methodology beyond the simple
solute/solvent case described above. When there is a
covalent bond between the portion of the system to be
modeled by quantum mechanics, and the portion which is to
be treated by EFP multipole and polarizability terms, an
additional layer is needed in the model. The covalent
linkage is not so simple as the interactions between closed
shell solute and solvent molecules. The "buffer zone"
between the quantum mechanics and the EFP consists of
frozen nuclei, and frozen localized orbitals, so that the
quantum mechanical region sees a orbital representation of
the closest particles, and multipoles etc. beyond that.
Since the orbitals in the buffer zone are frozen, it need
extend only over a few atoms in order to keep the orbitals
in the fully optimized quantum region within that region.
The general outline of this kind of computation is as
follows:
a) a full quantum mechanics computation on a system
containing the quantum region, the buffer region,
and a few atoms into the EFP region, to obtain the
frozen localized orbitals in the buffer zone.
This is called the "truncation run".
b) a full quantum mechanics computation on a system
with all quantum region atoms removed, and with
the frozen localized orbitals in the buffer zone.
The necessary multipole and polarizability data
to construct the EFP that will describes the EFP
region will be extracted from the wavefunction.
This is called the "MAKEFP run". It is possible
to use several such runs if the total EFP region
is quite large.
c) The intended QM/MM run(s), after combining the
information from these first two types of runs.
As an example, consider a protonated lysine residue
which one might want to consider quantum mechanically in a
protein whose larger parts are to be treated with an EFP.
The protonated lysine is
NH2
+ /
H3N(CH2)(CH2)(CH2)--(CH2)(CH)
\
COOH
The bonds which you see drawn show how the molecule is
partitioned between the quantum mechanical side chain, a
CH2CH group in the buffer zone, and eventually two
different EFPs may be substituted in the area of the NH2
and COOH groups to form the protein backbone.
The "truncation run" will be on the entire system as you
see it, with the 13 atoms in the side chain first in $DATA,
the 5 atoms in the buffer zone next in $DATA, and the
simplified EFP region at the end. This run will compute
the full quantum wavefunction by RUNTYP=ENERGY, followed by
the calculation of localized orbitals, and then truncation
of the localized orbitals that are found in the buffer zone
so that they contain no contribution from AOs outside the
buffer zone. The key input groups for this run are
$contrl
$truncn doproj=.true. plain=.true. natab=13 natbf=5 $end
This will generate a total of 6 localized molecular
orbitals in the buffer zone (one CC, three CH, two 1s inner
shells), expanded in terms of atomic orbitals located only
on those atoms.
The truncation run prepares template input files for
the next run, including adjustments of nuclear charges at
boundaries, etc.
The "MAKEFP" run drops all 13 atoms in the quantum
region, and uses the frozen orbitals just prepared to
obtain a wavefunction for the EFP region. The carbon atom
in the buffer zone that is connected to the now absent QM
region will have its nuclear charge changed from 6 to 5 to
account for a missing electron. The key input for this
RUNTYP=MAKEFP job is the six orbitals in $VEC, plus the
groups
$guess guess=huckel insorb=6 $end
$mofrz frz=.true. ifrz(1)=1,2,3,4,5,6 $end
$stone
QMMMbuf
$end
which will cause the wavefunction optimization for the
remaining atoms to optimize orbitals only in the NH2 and
COOH pieces. After this wavefunction is found, the run
extracts the EFP information needed for the QM/MM third
run(s). This means running the Stone analysis for
distributed multipoles, and obtaining a polarizability
tensor for each localized orbital in the EFP region.
The QM/MM run might be RUNTYP=OPTIMIZE, etc. depending
on what you want to do with the quantum atoms, and its
$DATA group will contain both the 13 fully optimized atoms,
and the 5 buffer atoms, and a basis set will exist on both
sets of atoms. The carbon atom in the buffer zone that
borders the EFP region will have its nuclear charge set to
4 since now two bonding electrons to the EFP region are
lost. $VEC input will provide the six frozen orbitals in
the buffer zone. The EFP atoms are defined in a fragment
potential group.
The QM/MM run could use RHF or ROHF wavefunctions, to
geometry optimize the locations of the quantum atoms (but
not of course the frozen buffer zone or the EFP piece). It
could remove the proton to compute the proton affinity at
that terminal nitrogen, hunt for transition states, and so
on. Presently the gradient for GVB and MCSCF is not quite
right, so their use is discouraged.
Input to control the QM/MM preparation is $TRUNCN and
$MOFRZ groups. There are a number of other parameters in
various groups, namely QMMMBUF in $STONE, MOIDON and POLNUM
in $LOCAL, NBUFFMO in $EFRAG, and INSORB in $GUESS that are
relevant to this kind of computation. For RUNTYP=MAKEFP,
the biggest choices are LOCAL=RUEDENBRG vs. BOYS, and
POLNUM in $LOCAL, otherwise this is pretty much a standard
RUNTYP=ENERGY input file.
Source code distributions of GAMESS contain a directory
named ~/gamess/tools/efp, which has various tools for EFP
manipulation in it, described in file readme.1st. A full
input file for the protonated lysine molecule is included,
with instructions about how to proceed to the next steps.
Tips on more specialized input possibilities are appended
to the file readme.1st.
Simpler potentials
Since the EFP model's electrostatics is a set of
distributed multipoles (monopole to octopole) and
distributed polarizabilities (dipole), it is possible to
generate some water potentials found in the literature by
setting many EFP terms to zero. It is also necessary to
provide a Lennard-Jones 6-12 repulsive potential, and then
make a choice to follow the EFP1 type formula for QM/EFP
repulsion. Accordingly, EFP1 type calculations can be made
with the following water potentials,
FRAGNAME=SPC, SPCE, TIP5P, TIP5PE, or POL5P
The Wikipedia page
http://en.wikipedia.org/wiki/Water_model
defines the first four of these, which are not polarizable
potentials. The same web site references the primary
literature, so that is not repeated here. POL5P is a
polarizable potential, with parameters given by
D.Si and H.Li J.Chem.Phys. 133, 144112/1-8(2010)
references
The first paper is more descriptive, while the second
presents a very detailed derivation of the EFP1 method.
Reference 18 is an overview article on EFP2. Reference 44
is the most recent review.
The model development papers are: 1, 2, 3, 5, 6, 9, 14,
17, 18, 22, 23, 26, 31, 32, 34, 39, 41, 42, 43, 44, 46, 50,
51, 55, 57, 58.
1. "Effective fragment method for modeling intermolecular
hydrogen bonding effects on quantum mechanical
calculations"
J.H.Jensen, P.N.Day, M.S.Gordon, H.Basch, D.Cohen,
D.R.Garmer, M.Krauss, W.J.Stevens in "Modeling the
Hydrogen Bond" (D.A. Smith, ed.) ACS Symposium Series
569, 1994, pp 139-151.
2. "An effective fragment method for modeling solvent
effects in quantum mechanical calculations".
P.N.Day, J.H.Jensen, M.S.Gordon, S.P.Webb,
W.J.Stevens, M.Krauss, D.Garmer, H.Basch, D.Cohen
J.Chem.Phys. 105, 1968-1986(1996).
3. "The effective fragment model for solvation: internal
rotation in formamide"
W.Chen, M.S.Gordon, J.Chem.Phys., 105, 11081-90(1996)
4. "Transphosphorylation catalyzed by ribonuclease A:
Computational study using ab initio EFPs"
B.D.Wladkowski, M. Krauss, W.J.Stevens
J.Am.Chem.Soc. 117, 10537-10545(1995)
5. "Modeling intermolecular exchange integrals between
nonorthogonal orbitals"
J.H.Jensen J.Chem.Phys. 104, 7795-7796(1996)
6. "An approximate formula for the intermolecular Pauli
repulsion between closed shell molecules"
J.H.Jensen, M.S.Gordon Mol.Phys. 89, 1313-1325(1996)
7. "A study of aqueous glutamic acid using the effective
fragment potential model"
P.N.Day, R.Pachter J.Chem.Phys. 107, 2990-9(1997)
8. "Solvation and the excited states of formamide"
M.Krauss, S.P.Webb J.Chem.Phys. 107, 5771-5(1997)
9. "An approximate formula for the intermolecular Pauli
repulsion between closed shell molecules. Application
to the effective fragment potential method"
J.H.Jensen, M.S.Gordon
J.Chem.Phys. 108, 4772-4782(1998)
10. "Study of small water clusters using the effective
fragment potential method"
G.N.Merrill, M.S.Gordon J.Phys.Chem.A 102, 2650-7(1998)
11. "Solvation of the Menshutkin Reaction: A Rigourous
test of the Effective Fragement Model"
S.P.Webb, M.S.Gordon J.Phys.Chem.A 103, 1265-73(1999)
12. "Evaluation of the charge penetration energy between
nonorthogonal molecular orbitals using the Spherical
Gaussian Overlap approximation"
V.Kairys, J.H.Jensen
Chem.Phys.Lett. 315, 140-144(1999)
13. "Solvation of Sodium Chloride: EFP study of NaCl(H2O)n"
C.P.Petersen, M.S.Gordon
J.Phys.Chem.A 103, 4162-6(1999)
14. "QM/MM boundaries across covalent bonds: frozen LMO
based approach for the Effective Fragment Potential
method"
V.Kairys, J.H.Jensen J.Phys.Chem.A 104, 6656-65(2000)
15. "A study of water clusters using the effective fragment
potential and Monte Carlo simulated annealing"
P.N.Day, R.Pachter, M.S.Gordon, G.N.Merrill
J.Chem.Phys. 112, 2063-73(2000)
16. "A combined discrete/continuum solvation model:
Application to glycine" P.Bandyopadhyay, M.S.Gordon
J.Chem.Phys. 113, 1104-9(2000)
17. "Evaluation of charge penetration between distributed
multipolar expansions"
M.A.Freitag, M.S.Gordon, J.H.Jensen, W.J.Stevens
J.Chem.Phys. 112, 7300-7306(2000)
18. "The Effective Fragment Potential Method: a QM-based MM
approach to modeling environmental effects in
chemistry"
M.S.Gordon, M.A.Freitag, P.Bandyopadhyay, J.H.Jensen,
V.Kairys, W.J.Stevens J.Phys.Chem.A 105, 293-307(2001)
19. "Accurate Intraprotein Electrostatics derived from
first principles: EFP study of proton affinities of
lysine 55 and tyrosine 20 in Turkey Ovomucoid"
R.M.Minikis, V.Kairys, J.H.Jensen
J.Phys.Chem.A 105, 3829-3837(2001)
20. "Active site structure & mechanism of Human Glyoxalase"
U.Richter, M.Krauss J.Am.Chem.Soc. 123, 6973-6982(2001)
21. "Solvent effect on the global and atomic DFT-based
reactivity descriptors using the EFP model. Solvation
of ammonia." R.Balawender, B.Safi, P.Geerlings
J.Phys.Chem.A 105, 6703-6710(2001)
22. "Intermolecular exchange-induction and charge transfer:
Derivation of approximate formulas using nonorthogonal
localized molecular orbitals."
J.H.Jensen J.Chem.Phys. 114, 8775-8783(2001)
23. "An integrated effective fragment-polarizable continuum
approach to solvation: Theory & application to glycine"
P.Bandyopadhyay, M.S.Gordon, B.Mennucci, J.Tomasi
J.Chem.Phys. 116, 5023-5032(2002)
24. "The prediction of protein pKa's using QM/MM: the pKa
of Lysine 55 in turkey ovomucoid third domain"
H.Li, A.W.Hains, J.E.Everts, A.D.Robertson, J.H.Jensen
J.Phys.Chem.B 106, 3486-3494(2002)
25. "Computational studies of aliphatic amine basicity"
D.C.Caskey, R.Damrauer, D.McGoff
J.Org.Chem. 67, 5098-5105(2002)
26. "Density Functional Theory based Effective Fragment
Potential" I.Adamovic, M.A.Freitag, M.S.Gordon
J.Chem.Phys. 118, 6725-6732(2003)
27. "Intraprotein electrostatics derived from first
principles: Divid-and-conquer approaches for QM/MM
calculations" P.A.Molina, H.Li, J.H.Jensen
J.Comput.Chem. 24, 1971-1979(2003)
28. "Formation of alkali metal/alkaline earth cation water
clusters, M(H2O)1-6, M=Li+, K+, Mg+2, Ca+2: an
effective fragment potential caase study"
G.N.Merrill, S.P.Webb, D.B.Bivin
J.Phys.Chem.A 107, 386-396(2003)
29. "Anion-water clusters A-(H2O)1-6, A=OH, F, SH, Cl, and
Br. An effective fragment potential test case"
G.N.Merrill, S.P.Webb
J.Phys.Chem.A 107,7852-7860(2003)
30. "The application of the Effective Fragment Potential to
molecular anion solvation: a study of ten oxyanion-
water clusters, A-(H2O)1-4"
G.N.Merrill, S.P.Webb J.Phys.Chem.A 108, 833-839(2004)
31. "The effective fragment potential: small clusters and
radial distribution functions"
H.M.Netzloff, M.S.Gordon J.Chem.Phys. 121, 2711-4(2004)
32. "Fast fragments: the development of a parallel
effective fragment potential method"
H.M.Netzloff, M.S.Gordon
J.Comput.Chem. 25, 1926-36(2004)
33. "Theoretical investigations of acetylcholine (Ach) and
acetylthiocholine (ATCh) using ab initio and effective
fragment potential methods"
J.Song, M.S.Gordon, C.A.Deakyne, W.Zheng
J.Phys.Chem.A 108, 11419-11432(2004)
34. "Dynamic polarizability, dispersion coefficient C6, and
dispersion energy in the effective fragment potential
method"
I.Adamovic, M.S.Gordon Mol.Phys. 103, 379-387(2005)
35. "Solvent effects on the SN2 reaction: Application of
the density functional theory-based effective fragment
potential method"
I.Adamovic, M.S.Gordon J.Phys.Chem.A 109, 1629-36(2005)
36. "Theoretical study of the solvation of fluorine and
chlorine anions by water"
D.D.Kemp, M.S.Gordon J.Phys.Chem.A 109, 7688-99(2005)
37. "Modeling styrene-styrene interactions"
I.Adamovic, H.Li, M.H.Lamm, M.S.Gordon
J.Phys.Chem.A 110, 519-525(2006)
38. "Methanol-water mixtures: a microsolvation study using
the Effective Fragment Potential method"
I.Adamovic, M.S.Gordon
J.Phys.Chem.A 110, 10267-10273(2006)
39. "Charge transfer interaction in the effective fragment
potential method" H.Li, M.S.Gordon, J.H.Jensen
J.Chem.Phys. 124, 214108/1-16(2006)
40. "Incremental solvation of nonionized and zwitterionic
glycine"
C.M.Aikens, M.S.Gordon
J.Am.Chem.Soc. 128, 12835-12850(2006)
41. "Gradients of the polarization energy in the Effective
Fragment Potential method"
H.Li, H.M.Netzloff, M.S.Gordon
J.Chem.Phys. 125, 194103/1-9(2006)
42. "Electrostatic energy in the Effective Fragment
Potential method: Theory and application to benzene
dimer"
L.V.Slipchenko, M.S.Gordon
J.Comput.Chem. 28, 276-291(2007)
43. "Polarization energy gradients in combined Quantum
Mechanics, Effective Fragment Potential, and
Polarizable Continuum Model Calculations"
H.Li, M.S.Gordon J.Chem.Phys. 126, 124112/1-10(2007)
44. "The Effective Fragment Potential: a general method
for predicting intermolecular interactions"
M.S.Gordon, L.V.Slipchenko, H.Li, J.H.Jensen
Annual Reports in Computational Chemistry, Volume 3,
pp 177-193 (2007).
45. "An Interpretation of the Enhancement of the Water
Dipole Moment Due to the Presence of Other Water
Molecules"
D.D.Kemp, M.S.Gordon
J.Phys.Chem.A 112, 4885-4894(2008)
46. "Solvent effects on optical properties of molecules: a
combined time-dependent density functional/effective
fragment potential approach"
S.Yoo, F.Zahariev, S.Sok, M.S.Gordon
J.Chem.Phys. 129, 144112/1-8(2008)
47. "Modeling pi-pi interactions with the effective
fragment potential method: The benzene dimer and
substituents"
T.Smith, L.V.Slipchenko, M.S.Gordon
J.Phys.Chem.A 112, 5286-5294(2008)
48. "Water-benzene interactions: An effective fragment
potential and correlated quantum chemistry study"
L.V.Slipchenko, M.S.Gordon
J.Phys.Chem.A 113, 2092-2102(2009)
49. "Ab initio QM/MM excited-state molecular dynamics study
of Coumarin 151 in water solution"
D.Kina, P.Arora, A.Nakayama, T.Noro, M.S.Gordon,
T.Taketsugu Int.J.Quantum Chem. 109, 2308-2318(2009)
50. "Damping functions in the effective fragment potential
method
L.V.Slipchenko, M.S.Gordon
Mol.Phys. 197, 999-1016 (2009)
51. "A combined effective fragment potential-fragment
molecular orbital method. 1. the energy expression"
T.Nagata, D.G.Fedorov, K.Kitaura, M.S.Gordon
J.Chem.Phys. 131, 024101/1-12(2009)
52. "Alanine: then there was water"
J.M.Mullin, M.S.Gordon
J.Phys.Chem.B 113, 8657-8669(2009)
53. "Water and Alanine: from puddles(32) to ponds(49)"
J.M.Mullin, M.S.Gordon
J.Phys.Chem.B 113, 14413-14420(2009)
54. "Structure of large nitrate-water clusters at ambient
temperatures: simulations with effective fragment
potentials and force fields with implications for
atmospheric chemistry"
Y.Miller, J.L.Thoman, D.D.Kemp, B.J.Finlayson-Pitts,
M.S.Gordon, D.J.Tobias, R.B.Gerber
J.Phys.Chem.A 113, 12805-12814(2009)
55. "Quantum mechanical/molecular mechanical/continuum
style solvation model: linear response theory,
variational treatment, and nuclear gradients"
H.Li J.Chem.Phys. 131, 184103/1-8(2009)
56. "Aqueous solvation of bihalide anions"
D.D.Kemp, M.S.Gordon
J.Phys.Chem.A 114, 1298-1303(2010)
57. "Exchange repulsion between effective fragment
potentials and ab initio molecules"
D.D.Kemp, J.M.Rintelman, M.S.Gordon, J.H.Jensen
Theoret.Chem.Acc. 125, 481-491(2010).
58. Modeling Solvent Effects on Electronic Excited States
A.DeFusco, N.Minezawa, L.V.Slipchenko, F.Zahariev,
M.S.Gordon J.Phys.Chem.Lett. 2, 2184-2192(2011)
The Fragment Molecular Orbital method
The method was proposed by Professor Kitaura and coworkers
in 1999, based on the Energy Decomposition Analysis (EDA,
sometimes called the Morokuma-Kitaura energy
decomposition). The FMO method is completely independent of
and bears no relation to:
1. Frontier molecular orbitals (FMO),
2. Fragment molecular orbitals (FMO).
The latter name is often used for the process of
construction of full molecular orbitals by combining MO
diagrams for parts of a molecule, ala Roald Hoffmann.
The effective fragment molecular orbital method (EFMO) is
closely related to but also bears significant difference to
FMO, and discussed below.
The FMO program was interfaced with GAMESS and follows
general GAMESS guidelines for code distribution and usage.
The users of the FMO program are requested to cite the
FMO3-RHF paper as the basic FMO reference,
D.G. Fedorov, K. Kitaura,
J. Chem. Phys. 120, 6832-6840(2004)
and other papers as appropriate (see below).
The basic idea of the method is to acknowledge the fact the
exchange and self-consistency are local in most molecules
(and clusters and molecular crystals), which permits
treating remote parts with Coulomb operators only, ignoring
the exchange. This idea further evolves into doing
molecular calculations, piecewise, with Coulomb fields due
to the remaining parts. In practice one divides the
molecule into fragments and performs n-mer calculations of
these in the Coulomb field of other fragments (n=1,2,3).
There are no empirical parameters, and the only departure
from ab initio rigor is the subjective fragmentation. It
has been observed that if performed physically reasonably,
the fragmentation scheme alters the results very little.
What changes the accuracy the most is the fragment size,
which also determines the computational efficiency of the
method.
The first question is how to get started. The easiest way
to prepare an FMO input file for GAMESS is to use the free
GUI software Facio, developed by M. Suenaga at Kyushu
University. It can do molecular modeling, automatic
fragmentation of peptides, nucleotides and saccharides and
create GAMESS/FMO input files:
http://www1.bbiq.jp/zzzfelis/Facio.html
If you have a cluster of identical molecules, you can
perform fragmentation with just one keyword ($FMO NACUT=).
Computationally, it is always better to partition in a
geometrical way (close parts together), so that the
distance-based approximations are more efficient. The
accuracy depends mainly upon the locality of the density
distribution, and the appropriateness of partitioning it
into fragments. There is no simple connexion between the
geometrical proximity of fragmentation and accuracy.
For basis sets you should use general guidelines and your
experience developed for ab initio methods. There is a file
provided (HMOs.txt) that contains hybrid molecular orbitals
(HMO) used to divide the MO space along fragmentation
points at covalent bonds. If your basis set is not there
you need to construct your own set of HMOs. See the example
file makeLMO.inp for this purpose.
Next you choose a wave function type. At present one can
use DFTB, RHF, DFT, TDDFT, MP2, CC, and MCSCF
(all except MCSCF support the 3-body expansion).
Geometry optimization can be
performed with all of these methods, except CC.
Note that presence of $FMO turns FMO on.
Surfaces and solids
Until 2008, for treating covalently connected fragments,
FMO had fully relaxed electron density of the detached
bonds. This method is now known as FMO/HOP (HOP=hybrid
orbital projection operator). It allows for a full
polarization of the system and is thus well suited to very
polar systems, such as proteins with charged residues. In
2008, an alternative fragmentation was suggested, based on
adaptive frozen orbitals (AFO), FMO/AFO. In it, the
electron density for each detached bond is first computed
in the automatically generated small model system (with the
bond intact), and in the FMO fragment calculations this
electron density is frozen. It was found that FMO/AFO works
quite well for surfaces and solids, where there is a dense
network of bonds to be detached in order to define
fragments (and the detached bonds interact quite strongly).
In addition, by restricting the polarization, FMO/AFO was
found to give a more balanced properties for large basis
sets (triple-zeta with polarization or larger), or in
comparing different isomers. However, for proteins with
charged residues the original FMO/HOP scheme has a better
accuracy (except large basis sets). At this point, FMO/AFO
was applied to zeolites only, and some more experience is
needed to give more practical advice to applications.
FMO/AFO is turned on by a nonzero rafo(1) parameter (rafo
array provides the thresholds to build model systems).
FMO variants
In 2007, Dahlke et al. introduced the Electrostatically
Embedded Many-Body Expansion method (see E. E. Dahlke and
D. G. Truhlar, J. Chem. Theory Comput. 4, 1-6 (2008) for
more recent work). This method is essentially FMO with the
RESPPC approximation (point charges for the electrostatic
field) applied to all fragments, with the further provision
that these charges may be defined at will (whereas RESPPC
uses Mulliken charges), and they are kept frozen (not
optimized, as in FMO). Next, Kamiya et al. suggested a fast
electron correlation method (M. Kamiya, S. Hirata, M.
Valiev, J. Chem. Phys. 128, 074103 (2008)), where again FMO
with the RESPPC approximation to all fragments is applied
with the further provision that the charges are derived
from the electrostatic potential (so called ESP charges),
and BSSE correction is added. The Dahlke's method was
generalized in GAMESS with the introduction of an arbitrary
hybrid approach, in which some fragments may have fixed and
some variationally optimized charges. This implementation
was employed in FMO-TDDFT calculations of solid state
quinacridone (see Ref. 16 below) by using DFT/PBC frozen
charges. The present energy only implementation is mostly
intended for such cases as that (i.e., TDDFT), and some
more work is needed to finish it for general calculations.
To turn this on, set RESPPC=-1 and define NOPFRG for frozen
charge fragments to 64, set frozen charges in ATCHRG.
Another FMO-like method is EFMO, see its own subsection
below. EFMO itself is related to several methods (PMISP: P.
Soederhjelm, U. Ryde, J. Phys. Chem. A 2009, 113, 617?627;
another is G. J. O. Beran, J. Chem. Phys. 2009, 130,
164115).
Effective fragment molecular orbital method (EFMO)
EFMO has been formulated by combining the physical models
in EFP and FMO, namely, in EFMO, fragments are computed
without the ESP (of FMO), and the polarization is estimated
using EFP models of fragment polarizabilities, which are
computed on the fly, so this can be thought of as
automatically generated potentials in EFP. Consequently,
close dimers are computed quantum-mechanically (without
ESP) and far dimers are computed using the electrostatic
multipole models of EFP. At present, only vacuum closed-
shell RHF and DFT are supported, for energy and gradient;
and only molecular clusters can be computed (no systems
with detached bonds). From the user point of view, EFMO
functionality is very intensively borrowed from FMO, and
the calculation setup is almost identical. Most additional
physical models such as PCM are not supported in EFMO. EFMO
should not be confused with FMO/EFP. The latter uses FMO
for some fragments and EFP for others. EFMO uses the same
model (EFMO), which is neither FMO nor EFP. For
approximations, EFMO at present has only RESDIM.
EFMO references are:
1. Effective Fragment Molecular Orbital Method: A Merger of
the Effective Fragment Potential and Fragment Molecular
Orbital Methods.
C. Steinmann, D. G. Fedorov, J. H. Jensen
J. Phys. Chem. A 114, 8705-8712 (2010).
2. The Effective Fragment Molecular Orbital Method for
Fragments Connected by Covalent Bonds.
C. Steinmann, D. G. Fedorov, J. H. Jensen
PLoS One, 7, e41117(2012).
3. Mapping enzymatic catalysis using the effective fragment
molecular orbital method: towards all ab initio
biochemistry.
C. Steinmann, D. G. Fedorov, J. H. Jensen
PLoS One 8, e60602 (2013).
Guidelines for approximations with FMO3
Three sets are suggested, for various accuracies:
low: resppc=2.5 resdim=2.5 ritrim(1)=0.001,-1,1.25
medium: resppc=2.5 resdim=3.25 ritrim(1)=1.25,-1,2.0
high: resppc=2.5 resdim=4.0 ritrim(1)=2,2,2
For correlated runs, add one more value to ritrim, equal to
the third element (i.e., 1.25 or 2.0). Note that gradient
runs do not support nonzero RESDIM and thus use RESDIM=0 if
gradient is to be computed. The "low" level of accuracy
for FMO3 has an error versus full ab initio similar to
FMO2, except for extended basis sets (6-311G** etc) where
it is substantially better than FMO2. Thus the low level is
only recommended for those large basis sets, and if a
better level cannot be afforded. The medium level is
recommended for production FMO3 runs; the high level is
mostly for accuracy evaluation in FMO development. The
cost is roughly: 3(low), 6(medium), 12(high). This means
that FMO3 with the medium level takes roughly six times
longer than FMO2.
Some of the default tolerances were changed as of January
2009, when FMO 3.2 was included in GAMESS. In general,
stricter parameters are now enforced when using FMO3, which
of course is intended to produce more accurate results. If
you wish to reproduce earlier results with the new code,
use the input to revert to the earlier values:
former -> FMO2 or FMO3 (as of 1/2009)
RESPPC: 2.0 2.0 2.50
RESDIM: 2.0 2.0 3.25
RCORSD: 2.0 2.0 3.25
RITRIM: 2.0,2.0,2.0,2.0 -> 1.25,-1.0,2.0,2.0 (FMO3 only)
MODESP: 1 0 1
MODGRD: 0 10 0
and two other settings which are not strictly speaking FMO
keywords may change FMO results:
MTHALL: 2 -> 4 (FMO/PCM only, see $TESCAV)
DFT grid: spherical -> Lebedev (FMO-DFT only, see $DFT)
Note that FMO2 energies printed during a FMO3 run will
differ from those in a FMO2 run, due to the different
tolerances used.
How to perform FMO-MCSCF calculations
Assuming that you are reasonably acquainted with ab initio
MCSCF, only FMO-specific points are highlighted. The active
space (the number of orbitals/electrons) is specified for
the MCSCF fragment. The number of core/virtual orbitals for
MCSCF dimers will be automatically computed. The most
important issue is the initial orbitals for the MCSCF
monomer. Just as for ab initio MCSCF, you should exercise
chemical knowledge and provide appropriate orbitals. There
are two basic ways to input MCSCF initial orbitals:
A) through the FMO monomer density binary file
B) by providing a text $VEC group.
The former way is briefly described in INPUT.DOC (see
orbital conversion). The latter way is really identical to
ab initio MCSCF, except the orbitals should be prepared for
the fragment (so in many cases you would have to get them
from an FMO calculation). Once you have the orbitals, put
them into $VEC1, and use the IJVEC option in $FMOPRP (e.g.,
if your MCSCF fragment is number 5, you would use $VEC1 and
ijvec(1)=5,0). For two-layer MCSCF the following
conditions apply. Usually one cannot simply use F40
restart, because its contents will be overwritten with RHF
orbitals and this will mess up your carefully chosen MCSCF
orbitals. Therefore, two ways exist. One is to modify A)
above by reordering the orbitals with something like
$guess guess=skip norder=1 iorder(28)=29,30,31,32,28 $end
Then the lower RHF layer will converge RHF orbitals that
you reorder with iorder in the same run (add 512 to nguess
in $FMO). This requires you know how to reorder before
running the job so it is not always convenient. Probably
the best way to run two-layer MCSCF is verbatim B) above,
so just provide MCSCF monomer orbitals in $VEC1. Finally,
it may happen that some MCSCF dimer will not converge.
Beside the usual MCSCF tricks to gain convergence as the
last resort you may be able to prepare good initial dimer
orbitals, put them into $VEC2 ($VEC3 etc) and read them
with ijvec. SOSCF is the preferred converger in FMO, and
the other one (FULLNR) has not been modified to eradicate
the artefacts of convergence (due to detached bonds). In
the bad cases you can try running one or two monomer SCF
iterations with FULLNR, stop the job and use its orbitals
in F40 to do a restart with SOSCF. We also found useful to
set CASDII=0.005 and nofo=10 in some cases running FOCAS
longer to get better orbitals for SOSCF.
How to perform multilayer runs
For some fragments you may like to specify a different
level of electron correlation and/or basis set. In a
typical case, you would use high level for the reaction
center and a lower level for the remaining part of the
system. The set up for multilayer runs is very similar to
the unilayer case. You only have to specify to what layer
each fragment belongs and for each layer define DFTTYP,
MPLEVL, SCFTYP as well as a basis set. If detached bonds
are present, appropriate HMOs should be defined. See the
paragraph above for multilayer MCSCF. Currently geometry
optimizations of multilayer runs require adding 128 to
NGUESS, if basis sets in layers differ from each other.
How to mix basis sets in FMO
You can mix basis sets in three ways.
a) single layer with ESP
You can specify different basis sets using the basis
set library in $DATA (such as C.1 and C.2 for CH3
and COO- carbons.
b) multilayer with ESP
Each layer can have its own basis set.
The difference between a 2-layer run with one basis set per
layer and a 1-layer run with 2-basis sets is significant:
in the former case the lower level densities are converged
with all fragments computed at the lower level. In the
latter case, the fragments are converged simultaneously,
each with its own basis set. In addition, dimer corrections
between layers will be computed differently: with the lower
basis set in the former case and with mixed basis set in
the latter. The latter approach may result in unphysical
polarization, so mixing basis sets is mainly intended to
add diffuse functions to anionic (e.g., carboxyl) groups,
not as a substitute for two-layer runs.
c) dual basis in an additive way
In this case, the whole system is calculated without
ESP and the smaller basis set, then with ESP and the smaller basis set,
and finally without ESP and the larger basis set.
This (c) way is often used with diffuse functions for
fragments with covalent boundaries.
How to perform FMO/PCM calculations
Solvent effects can be taken into account with PCM. PCM in
FMO is very similar to regular PCM. There is one basic
difference: in FMO/PCM the total electron density that
determines the electrostatic interaction is computed using
the FMO density expansion up to n-body terms. The cavity
is constructed surrounding the whole molecule, and the
whole cavity is used in each individual m-mer calculation.
There are several levels of accuracy (determined by the "n"
above), and the recommended level is FMO/PCM<1>,
specified by:
$pcm ief=-10 icomp=2 icav=1 idisp=1 ifmo=-1 $end
Many PCM options can be used as in the regular PCM. The
following restrictions apply:
IEF may be only -3 or -10, IDP must be 0.
Multilayer FMO runs are supported. Restarts are limited to
IREST=2, and in this case PCM charges (the ASCs) are not
recycled. However, the initial guess for the charges is
fairly reasonable, so IREST=2 may be useful although
reading the ASCs may be implemented in future.
Note for advanced users. IFMO < NBODY runs are permitted.
They are denoted by FMOm/PCM[n], where m=NBODY and n=IFMO.
In FMOm/PCM[n], the ASCs are computed with n-body level.
The difference between FMO2/PCM[1] and FMO2/PCM[1(2)] is
that in the former the ASCs are computed at the 1-body
level, whereas for the former at the 2-body level, but
without self-consistency (which would be FMO2/PCM[2]).
Probably, FMO3/PCM[2] should be regarded as the most
accurate and still affordable (with a few thousand nodes)
method. However, FMO3/PCM[1(2)] (specified with NBODY=3,
IFMO=2 and NPCMIT=2) is much cheaper and slightly less
accurate than FMO3/PCM[2]. FMO3/PCM[3] is the most
accurate and expensive level of all.
PCM<1> is very similar to PCM[1] except that the dimer
ES solvent contributions are halved.
How to perform FMO/EFP calculations
Solvent effects can also be taken into account with the
Effective Fragment Potential model. The presence of both
$FMO and $EFRAG groups selects FMO/EFP calculations. See
the $EFRAG group and the $FMO group for details.
In the FMO/EFP method, the Quantum Mechanical part of the
calculation in the usual EFP method is replaced by the FMO
method, which may save time for large molecules such as
proteins.
In the present version, only FMOn/EFP1 (water solvent only)
is available for RHF, DFT and MP2. One can use the MC
global optimization technique for FMO/EFP by RUNTYP=GLOBOP.
Of course, the group DDI (GDDI) parallelization technique
for the FMO method can be used.
Geometry optimization or saddle point search for FMO
The standard optimizers in GAMESS are now well
parallelized, and thus recommended to be used with FMO up
to the limit hardwired in GAMESS (2000 atoms). In practice,
if more than about 1000 atoms are present, numeric Hessian
updates often result in the improper curvature and
optimization stops. One can either do a restart, or use
RUNTYP=OPTFMO (which does not diagonalize the Hessian).
RUNTYP=OPTIMIZE applies to Cartesian coordinates or DLC.
RUNTYP=OPTFMO works only with Cartesian coordinates. If
your system has more than 2000 atoms you can consider
RUNTYP=OPTFMO, which can now use Hessian updates and
provides reasonable way to optimize although it is not as
good as the standard means in RUNTYP=OPTIMIZE.
A transition state search for FMO can be performed with
RUNTYP=SADPOINT using either Cartesian coordinates or DLC.
IRC calculations can be performed.
FMO hessian calculations
Analytic FMO Hessian with RUNTYP=HESSIAN may be computed
for RHF, ROHF, UHF and RDFT, in the gas phase (no
PCM, EFP etc).
Molecular dynamics with FMO
MD can be run for any FMO method, which has the gradient
implemented. However, in many cases the approximations in
the gradient for a particular method may lead to large
discrepancies in MD. The following methods have a fully
analytic gradient (which has to be turned on with $FMO
keyword MODGRD=42): FMO-RHF, FMO-DFT, FMO-MP2, FMO-RHF/EFP;
The RESPPC approximation adds a small error to the gradient,
so in MD larger values such as 3.0 should be used,
or 0.0 to turn RESPPC off.
Pair interaction energy decomposition analysis (PIEDA)
PIEDA can be performed for the PL0 and PL states. The PL0
state is the electronic state in which fragments are
polarised by the environment in its free (noninteracting)
state. The simplest example is that in a water cluster,
each molecule is computed in the electrostatic field
exerted by the electron densities of free water molecules.
The PL state is the FMO converged monomer state, that is,
the state in which fragments are polarised by the self-
consistently converged environment. Namely, following the
FMO prescription, fragments are recomputed in the external
field, until the latter converges. Using the PL0 state
requires a series of separate runs; and it also relies on a
"free state" which can be defined in many ways for
molecules with detached covalent bonds.
What should be done to do the PL0 state analysis?
1. run FMO0.
This computes the free state for each fragment, and those
electron densities are stored on file 30 (to be renamed
file 40 and reused in step 3).
2. compute BDA energies (if detached bonds are present),
using sample files in tools/fmo/pieda. This corrects for
artifacts of bond detaching, and involves running a model
system like H3C-CH3, to amend for C-C bond detaching.
3. Using results of (1) and (2), one can do the PL0
analysis. In addition to pasting the data from the two
punch files in steps 1,2 and the density file in step 1
should be provided.
What should be done to do the PL state analysis? The PL
state itself does not need either the free state or PL0
results. However, if the PL0 results are available,
coupling terms can be computed, and in this case IPIEDA is
set to 2; otherwise to 1.
So the easiest and frequently sufficient way to run PIEDA
is to set IPIEDA=1 and do not provide any data from
preceding calculations. The result of a PIEDA calculation
is a set of pair interaction energies (interfragment
interaction energies), decomposed into electrostatic,
exchange-repulsion, charge transfer and dispersion
contributions.
Finally, PIEDA (especially for the PL state) can be thought
of as FMO-EDA, EDA being the Kitaura-Morokuma decomposition
(RUNTYP=MOROKUMA). In fact, PIEDA (for the PL state) in
the case of just two fragments of standalone molecules is
entirely equivalent to EDA, which can be easily verified,
by running the full PIEDA analysis (ipieda=2). Note that
PIEDA can be run as direct SCF, whereas EDA cannot be, and
for large fragments PIEDA code can be used to perform EDA.
Also, EDA in GAMESS has no explicit dispersion.
In 2012, PIEDA/PCM was developed describing the solvent
screening. RO-PIEDA based on RO-(HF, MP2 or CC) may be
used for radicals. Grimme's dispersion models may be used
in PIEDA.
Excited states
At present, one can use CI, MCSCF, or TDDFT to compute
excited states in FMO. MCSCF is discussed separately
above, so here only TDDFT and CI are explained. They are
enabled by setting the IEXCIT option (EXCIT(1) defines the
excited state's fragment ID).
Two levels are implemented for TDDFT (FMO1-TDDFT and FMO2-
TDDFT). In the former, only monomer TDDFT calculations are
performed, whereas the latter adds pair corrections from
TDDFT dimers. PCM may be used for solvent effects with
TDDFT (PCM[1] is usually sufficient).
CI can only be done for CIS at the monomer level (nbody=1),
FMO1-CIS. The set-up for CI is similar to that for TDDFT.
Selective and subsystem FMO
Sometimes, one is interested only in some pair
interactions, for example, between ligand and protein, or
the opposite, only pair interactions within ligand. This
saves a lot of CPU time by omitting all other pair
calculations, but does not give the total properties. To
use this feature, define MOLFRG and MODMOL. RUNTYP=ENERGY
only is implemented.
In the subsystem analysis, one can divide fragments into
subsystems and obtain various properties of subsystems.
Frozen domain
To accelerate geometry optimisations, one can specify that
the electronic state of the first layer in a 2-layer FMO
can be computed at the initial geometry and consequently be
frozen. One can define the polarizable buffer (equal to
layer 2) and frozen domain (layer 1). Fragments in the
polarizable buffer which contain the atoms active in
geometry optimisation form the active domain. The
fragments in the active domain should have a nonzero
separation from the frozen domain. In FMO/FD all dimers in
the polarizable buffer are computed; in FMO/FDD only those
dimers which have at least one monomer in the active domain
are computed. FMO/FD and FNI/FDD are only implemented for
RUNTYP=OPTIMIZE. MODFD and IACTAT in $FMO specify
FMO/FD(D) atop of the usual multilayer FMO setup with some
atoms frozen in geometry optimization by the standard means
(i.e., IACTAT in $STATPT). Note that in FMO/FD(D) the
Hessian as used in RUNTYP=OPTIMIZE is formed only for the
atoms in the second layer, so this upper layer should not
have more than the GAMESS limit (currently, 2000 atoms).
IMOMM with FMO
IMOMM (namely, SIMOMM) calculations can be performed with
the "MO" in IMOMM treated using FMO, i.e., this is like
QM/MM but without electronic embedding of QM by MM.
This calculation uses Tinker, a plug-in source code,
available from the GAMESS web site. You should compile and
link in the Tinker plug-in by changing a single line in
comp/compall/lked,
set TINKER=false
into
set TINKER=true
In addition, you should change tinker_maxatm (the maximum
number of atoms in the whole system, as used by Tinker) in
mx_limits.src to match exactly the problem size selected in
Tinker's sizes.i include file. After changing this,
recompile and link GAMESS.
The input file style is in general like that of SIMOMM
(q.v.). Different from regular FMO, the atomic coordinates
are given in $TINXYZ, not in $FMOXYZ. The fragmentation in
FMO applies to QM atoms only, selected by IQMATM, and
numbered consequently in FMO, so that INDAT in $FMO applies
to the atoms renumbered from 1 (defined in IQMATM). Other
than $FMOXYZ being superceded by $TINXYZ, the rest of FMO
options is like in normal FMO. IMOMM based on FMO is
usually referred to as FMO/MM for short, although "FMO-
based SIMOMM" is probably easier to understand. The
somewhat tautological FMO-IMOMM has also been used by some.
Covalent boundaries between FMO and MM are supported (via
link atoms). FMO/MM can be used to run geometry
optimizations, whichis really what it is designed for.
Analyzing and visualizing the results
Annotated outputs provided in tools/fmo have matching
mathematical formulae added onto the outputs, for easier
reading.
Facio (http://www1.bbiq.jp/zzzfelis/Facio.html) can plot
various FMO properties such as interaction energies, using
interactive GUI viewers.
To plot orbitals for an n-mer, set NPUNCH=2 in $SCF and
PLTORB=.T. There are several ways to produce cube files
with electron densities.
Parallelization of FMO runs with GDDI
The FMO method has been developed within a 2-level
hierarchical parallelization scheme, group DDI (GDDI),
allowing massively parallel calculations. Different groups
of processors handle the various monomer, dimer, and maybe
trimer computations. The processor groups should be sized
so that GAMESS' innate scaling is fairly good, and the
fragments should be mapped onto the processor groups in an
intelligent fashion.
This is a very important and seemingly difficult issue. It
is very common to be able to speed up parallel runs at
least several times just by using GDDI better. First of
all, do not use plain DDI and always employ GDDI when
running FMO calculations. Next, learn that you can and
should divide nodes into groups to achieve better
performance. The very basic rule of thumb is to try to have
several times fewer groups than jobs. Since the number of
monomers and dimers is different, group division should
reflect that fact. Ergo, find a small parallel computer
with 8-32 nodes and experiment changing just two numbers:
ngrfmo(1)=N1,N2 and see how performance changes for your
particular system.
Limitations of the FMO method in GAMESS
1. Dimensions: in general none, except that the standard
GAMESS engines RUNTYP=OPTIMIZE and IRC are limited to 2000
atoms (for FD(D), domain B may not exceed this limit). The
limit can be increased by changing the source and
recompiling GAMESS (see elsewhere).
2. CHARMM may not be combined with FMO, and some other
extensions may not work. Not every illegal combination is
trapped, caveat emptor!
3. RUNTYP is limited to ENERGY, GRADIENT, OPTIMIZE, OPTFMO,
IRC, FMO0, MD, GLOBOP, SADPOINT, FMOHESS and RAMAN only:
Do not even try other ones!
4. Three-body FMO-MCSCF and FMO-TDDFT are not implemented.
5. No MOPAC semiempirical methods may be used, but DFTB was
interfaced with FMO..
What will work the same way as ab initio:
The various SCF convergers, all DFT functionals, in-core
integrals, direct SCF.
Restarts with the FMO method
RUNTYP=ENERGY can be restarted from anywhere before
trimers. To restart monomer SCF, copy file F40 with monomer
densities to the grandmaster node. To restart dimers,
provide file F40 and monomer energies ($FMOENM).
Optionally, some dimer energies can be supplied ($FMOEND)
to skip computation of corresponding dimers.
RUNTYP=GRADIENT can be easily restarted from monomer SCF
(which really means it is a restart of RUNTYP=ENERGY, since
gradient is computed at the end of this step). Provide
file F40. There is another restart option (1024 in $FMOPRP
irest=), supporting full gradient restart, requiring,
however, that you set this option in the original run
(whose results you use to restart). To use this option, you
would also need to keep (or save and restore) F38 files on
each node (they are different).
RUNTYP=OPTIMIZE can be restarted from anywhere within the
first RUNTYP=GRADIENT run (q.v.). In addition, by
replacing FMOXYZ group, one can restart at a different
geometry.
RUNTYP=OPTFMO can be restarted by providing a new set of
coordinates in $FMOXYZ and, optionally, by transferring
$OPTRST from the punch into the input file.
Note on accuracy
The FMO method is aimed at computation of large molecules.
This means that the total energy is large, for example, a
6646 atom molecule has a total energy of -165,676 Hartrees.
If one uses the standard accuracy of roughly 1e-9 (that
should be taken relatively), one gets an error as much as
0.001 hartree, in a single calculation. FMO involves many
ab initio single point calculations of fragments and their
n-mers, thus it can be expected that numeric accuracy is 1-
2 orders lower than that given by 1e-9. Therefore, it is
compulsory that accuracy should be raised, which is done by
default.
The following default parameters are reset by FMO:
ICUT/$CONTRL (9->12), ITOL/$CONTRL(20->24),
CONV/$SCF(1e-5 -> 1e-7),
CUTOFF/$MP2 (1e-9->1e-12), CUTTRF/$TRANS(1e-9->1e-10).
CVGTOL/$DET,$GUGDIA (1e-5 -> 1e-6)
This to some extent slows down the calculation (perhaps on
the order of 10-15%). It is suggested that you maintain
this accuracy for all final energetics. However, you may
be able to drop the accuracy a bit for the initial part of
geometry optimization if you are willing to do manual work
of adjusting accuracy in the input. It is recommended to
keep high accuracy at the flat surfaces (the final part of
optimizations) though. For DFT the numeric grid's accuracy
may be increased in accordance with the molecule size, e.g.
extending the default grid of 96*12*24 to 96*20*40.
However, some tests indicate that energy differences are
quite insensitive to this increase.
The principle FMO paper is:
Fragment molecular orbital method: an approximate
computational method for large molecules"
K. Kitaura, E. Ikeo, T. Asada, T. Nakano, M. Uebayasi
Chem. Phys. Lett., 313, 701(1999).
and the principal papers for FMO in GAMESS are:
1. D. G. Fedorov, K. Kitaura,
J. Chem. Phys. 120, 6832-6840(2004)
2. D. G. Fedorov, WIREs: Comp. Mol. Sc. 7 (2017) e1322.
Recent papers can be googled out. As the field is dynamically
changing, a list is not provided here.
EFMO references are given in its own subsection.
The Cluster-in-Molecules method
sequential and parallel execution
If the user is not interested in parallel CIM calculations,
MTDCIM must be set at 0, which is a default value, and no
additional steps have to be taken. If the user is
interested in a parallel execution, MTDCIM must initially
be set at 1 to prepare for individual subsystem GAMESS
runs, and then, after the individual subsystem runs are
completed, reset to 2 to complete the CIM calculation. When
MTDCIM is initially set at 1, multiple input files
$JOB.Sys-N.inp for individual subsystem calculations with
GAMESS, which can be run independent of one another, and
the $JOB.cim file, which contains the information about all
subsystems needed to complete the CIM calculation, are
automatically generated, and the program stops awaiting
further execution. Each subsystem N has then to be run as
an independent GAMESS calculation using the $JOB.Sys-N.inp
input file. This produces the $JOB.Sys-N.cim files which
contain the information about the correlation energy
contributions due to the occupied LMOs central in
subsystems. All $JOB.Sys-N.cim files resulting from the
individual subsystem calculations with GAMESS and the
$JOB.cim file are used to assemble the final results of a
CIM calculation for the entire system. In order to
accomplish this and complete the CIM run, one has to reset
MTDCIM in the main $JOB.inp input file to 2 and run GAMESS
again.
restarts
If any of the subsystem GAMESS calculations using the
$JOB.Sys-N.inp input files fails, the user can always rerun
it (editing the corresponding $JOB.Sys-N.inp file(s), if
need be), and then use MTDCIM=2 to complete the desired CIM
calculation for the entire system. This applies to
sequential and parallel CIM calculations. In the latter
case, this is a natural consequence of the way the parallel
execution is structured (see note 1). In the former case
(MTDCIM=0), if the entire calculation is completed, the
$JOB.Sys-N.* subsystem files are deleted, but if one of the
subsystem calculations fails, the program aborts, leaving
all $JOB.Sys-N.inp input files, the $JOB.Sys-N.cim output
files from the completed subsystem calculations, and the
$JOB.cim file on the disk. One can rerun the subsystem
GAMESS calculation that failed, editing the corresponding
$JOB.Sys-N.inp file if need be, and the remaining subsystem
calculations, and then, once all subsystem GAMESS runs are
completed, finish the calculation by using MTDCIM=2 in the
main $JOB.inp input. This has an advantage over the more
automated method of restarting the sequential CIM
calculations described below in that the Hartree-Fock and
orbital localization calculations for the entire system do
not have to be repeated.
If the sequential (MTDCIM=0) run does not complete due to
the failure of one of the subsystem GAMESS calculations,
one can also follow a simpler, more automated restart
strategy. At the time of failure, all $JOB.Sys-N.inp input
files, the $JOB.Sys-N.cim output files from the completed
subsystem calculations, and the $JOB.cim file are saved on
the disk. After inspecting the main output file, the user
can simply delete the $JOB.dat file and the $JOB.Sys-N.cim
file resulting from the failed subsystem calculation, edit
the corresponding $JOB.Sys-N.inp file, if necessary, and
rerun the GAMESS calculation with MTDCIM=0. The Hartree-
Fock and orbital localization calculations for the entire
system will be performed again, but the user will avoid the
need for running individual subsystem calculations one-by-
one, as described above.
the cimshell script
The Python script "cimshell" that automatically produces
typical $JOB.sh files for parallel OpenMP and MPI subsystem
calculations using the $JOB.Sys-N.inp files can be found in
$GMS_PATH/tools/cim/, where $GMS_PATH is the GAMESS main
directory. In order to run the subsystem calculations with
OpenMP or MPI, and with the help of the "cimshell" script,
the following steps should be performed:
1. Run GAMESS CIM calculation using MTDCIM=1 to produce
the $JOB.Sys-N.inp and $JOB.cim files.
2. After all subsystem $JOB.Sys-N.inp input files are
generated, use "cimshell" to automatically generate
the OpenMP or MPI script $JOB.sh for parallel
execution. By default, the "cimshell" program must
be run in the directory where the $JOB.inp and all
$JOB.Sys-N.inp files reside. For example, "cimshell
--np 4 $JOB" generates the $JOB.sh script for an
OpenMP parallel calculation on 4 cores, whereas
"cimshell --np 8 --para mpi --submit pbs --MPI_EXEC
mpiexec $JOB" generates the $JOB.sh script for an
MPI parallel calculation on 8 processors using the
PBS queue system for submitting the job. In these
two examples, we are assuming that the "cimshell"
has been copied to the directory where all
$JOB.Sys-N.inp files reside; this can be altered by
redefining the $PATH variable (adding
$GMS_PATH/tools/cim/ to it). Use "cimshell -h" for
more information about the "cimshell" options.
3. Run or submit the $JOB.sh script to have subsystem
calculations performed in parallel. In order to do
this, the user must compile the ompjob.for (the
OpenMP case) or mpijob.for (the MPI case) programs
that reside in $GMS_PATH/tools/cim/. The
corresponding executables used by $JOB.sh are called
ompjob and mpijob, respectively. The example of the
Makefile that can be used to install ompjob and
mpijob can be found in $GMS_PATH/tools/cim/.
4. After all subsystem calculations are completed, run
the final GAMESS CIM calculation using the original
input file $JOB.inp in which with MTDCIM=2. GAMESS
will automatically find the relevant $JOB.Sys-N.cim
and $JOB.cim files to complete the CIM calculation
for the entire system and print the final CIM
energies in the main output.
CIM references
THE FOLLOWING PAPERS SHOULD BE CITED WHEN USING
CLUSTER-IN-MOLECULE OPTIONS:
DUAL-ENVIRONMENT CIM (CIMTYP=DECIM)
W. LI, P. PIECUCH, J.R. GOUR, AND S. LI,
J. CHEM. PHYS. 131, 114109-1 - 114109-30 (2009).
SEE, ALSO, S. LI, J. SHEN, W. LI, AND Y. JIANG,
J. CHEM. PHYS. 125, 074109-1 - 074109-10 (2006)
SINGLE-ENVIRONMENT CIM (CIMTYP=SECIM,GSECIM)
W. LI, P. PIECUCH, J.R. GOUR, AND S. LI,
J. CHEM. PHYS. 131, 114109-1 - 114109-30 (2009);
W. LI AND P. PIECUCH, J. PHYS. CHEM. A 114, 8644-8657
(2010).
IN ADDITION, THE USE OF MULTI-LEVEL CIM SHOULD REFERENCE
W. LI AND P. PIECUCH, J. PHYS. CHEM. A 114, 6721-6727
(2010).
MOPAC Calculations within GAMESS
Parts of MOPAC 6.0 have been included in GAMESS giving
access to four semiempirical wavefunctions: MNDO, AM1,
PM3, and RM1. RM1 is the most recent parameterization,
replacing AM1 data for H, C-F, P-Cl, Br, and I. See G.
Bruno Rocha, R. Oliveira Freire, A. Mayall Simas, and
J.J.P.Stewart, J.Comput.Chem. 27, 1101-1111(2006).
These wavefunctions are quantum mechanical in nature
but neglect most two electron integrals, a deficiency that
is (hopefully) compensated for by introduction of empirical
parameters. The quantum mechanical nature of semiempirical
theory makes it quite compatible with the ab initio
methodology in GAMESS. As a result, very little of MOPAC
6.0 actually is incorporated into GAMESS. The part that
did survive is the code that evaluates
1) the one- and two-electron integrals,
2) the two-electron part of the Fock matrix,
3) the cartesian energy derivatives, and
4) the ZDO atomic charges and molecular dipole.
Everything else is actually GAMESS: coordinate input
(including point group symmetry), the SCF convergence
procedures, the matrix diagonalizer, the geometry searcher,
the numerical hessian driver, and so on. Most of the
output will look like an ab initio output.
It is extremely simple to perform these calculations.
All you need to do is specify GBASIS=MNDO, AM1, PM3, or RM1
in the $BASIS group. Note that this not only selects a
particular "Hamiltonian" (parameter set), it also picks a
Slater Type Orbital (STO) basis.
MOPAC parameters exist for the following elements. The
printout when you run will give you specific references for
each kind of atom. The quote on alkali's below means that
these elements are treated as "sparkles", rather than as
atoms with genuine basis functions.
For MNDO:
H
Li * B C N O F
Na' * Al Si P S Cl
K' * ... Zn * Ge * * Br
Rb' * ... * * Sn * * I
* * ... Hg * Pb *
For AM1: For PM3:
H H
* * B C N O F Li Be * C N O F
Na Mg Al Si P S Cl Na Mg Al Si P S Cl
K Ca ... Zn * Ge * * Br K Ca ... Zn Ga Ge As Se Br
Rb' * ... * * Sn * * I Rb' * ... Cd In Sn Sb Te I
* * ... Hg * * * * * ... Hg Tl Pb Bi
For RM1:
H
C N O F
P S Cl
Br
I
RM1 uses AM1 parameters for any element not listed above.
* * * * *
MOPAC will not work with every option in GAMESS: the
semiempirical wavefunctions must be RHF, UHF, and ROHF in
any combination with run types ENERGY, GRADIENT, OPTIMIZE,
SADPOINT, HESSIAN, and IRC. Note that nuclear hessian runs
use numerical finite differencing of analytic gradients.
MOPAC's CI and half electron methods are not supported.
Because the majority of the implementation is GAMESS
rather than MOPAC6, you will notice a few improvements.
Dynamic memory allocation is used, so GAMESS uses far less
memory for a given size of molecule. The starting orbitals
for SCF calculations are generated by a Huckel initial
guess routine. Spin restricted (high spin) ROHF can be
performed. Converged SCF orbitals will be labeled by their
symmetry type. Numerical hessians will make use of point
group symmetry, so that only the symmetry unique atoms need
to be displaced. Infrared intensities will be calculated
at the end of hessian runs. We have not at present used
the block diagonalizer during intermediate SCF iterations,
so that the run time for a single geometry point in GAMESS
is usually longer than in MOPAC. However, the geometry
optimizer in GAMESS can frequently optimize the structure
in fewer steps than the procedure in MOPAC. Orbitals and
hessians are punched out for convenient reuse in subsequent
calculations. Your molecular orbitals can be drawn with
the PLTORB graphics program, which has been taught about s
and p STO basis sets.
However, because of the STO basis set used in semi-
empirical runs, the various property calculations coded for
Gaussian (GTO) basis sets are unavailable. This means
$ELMOM, $ELPOT, etc. properties are unavailable. Note that
MOPAC6 did not include d STO integrals, so it is quite
impossible to run transition metals.
The PCM solvation model implemented in GAMESS can be
used with MOPAC runs by GAMESS.
To reduce CPU time, by default only the EXTRAP
convergence accelerator is used by the SCF procedures. For
difficult cases, the DIIS, RSTRCT, and/or SHIFT options
will work, but may add significantly to the run time. With
the Huckel guess procedure from GAMESS, most calculations
will converge acceptably without these special options.
The MOPAC implementation is able to run in parallel.
Semiempirical calculations are very fast. One of the
motives for the MOPAC implementation within GAMESS is to
take advantage of this speed. Semiempirical models can
rapidly provide reasonable starting geometries for ab
initio optimizations. Semiempirical hessian matrices are
obtained at virtually no computational cost, and may help
dramatically with an ab initio geometry optimization.
Simply use HESS=READ in $STATPT to use a MOPAC $HESS group
in an ab initio run.
It is important to exercise caution as semiempirical
methods can be dead wrong! The reasons for this are bad
parameters (in certain chemical situations), and the
underlying minimal basis set. A good question to ask
before using MOPAC is "how well is my system modeled by an
ab initio minimal basis set, such as STO-3G"? If the
answer is "not very well", there is a good chance that a
semiempirical description is equally poor.
Molecular Properties and Conversion Factors
These two papers are of general interest:
A.D.Buckingham, J.Chem.Phys. 30, 1580-1585(1959).
D.Neumann, J.W.Moskowitz J.Chem.Phys. 49, 2056-2070(1968).
The first deals with multipoles, and the second with other
properties such as electrostatic potentials.
All units are derived from the atomic units for distance
and the monopole electric charge, as given below.
distance 1 au = 5.291771E-09 cm
monopole 1 au = 4.803242E-10 esu
1 esu = sqrt(g-cm**3)/sec
dipole 1 au = 2.541766E-18 esu-cm
1 Debye = 1.0E-18 esu-cm
quadrupole 1 au = 1.345044E-26 esu-cm**2
1 Buckingham = 1.0E-26 esu-cm**2
octopole 1 au = 7.117668E-35 esu-cm**3
electric potential 1 au = 9.076814E-02 esu/cm
electric field 1 au = 1.715270E+07 esu/cm**2
1 esu/cm**2 = 1 dyne/esu
electric field gradient 1 au = 3.241390E+15 esu/cm**3
The atomic unit for the total electron density is
electron/bohr**3, but 1/bohr**3 for an orbital density.
The atomic unit for spin density is excess alpha spins per
unit volume, h/4*pi*bohr**3. Only the expectation value is
computed, with no constants premultiplying it.
IR intensities are printed in Debye**2/amu-Angstrom**2.
These can be converted into intensities as defined by
Wilson, Decius, and Cross's equation 7.9.25, in km/mole, by
multiplying by 42.255. If you prefer 1/atm-cm**2, use a
conversion factor of 171.65 instead. A good reference for
deciphering these units is A.Komornicki, R.L.Jaffe
J.Chem.Phys. 1979, 71, 2150-2155. A reference showing how
IR intensities change with basis improvements at the HF
level is Y.Yamaguchi, M.Frisch, J.Gaw, H.F.Schaefer,
J.S.Binkley, J.Chem.Phys. 1986, 84, 2262-2278.
Raman activities in A**4/amu multiply by 6.0220E-09 for
units of cm**4/g. One of the many sources explaining how
activity relates to intensity is D.Michalska, R.Wysokinski
Chem.Phys.Lett. 403, 211-217(2005)
Polarizabilities
Static polarizabilities are named alpha, beta, and gamma;
these are called the polarizability, hyperpolarizability,
and second hyperpolarizability. They are the 2nd, 3rd, and
4th derivatives of the energy with respect to uniform
applied electric fields, with the 1st derivative being the
dipole moment.
It is worth mentioning that a uniform (static) electric
field can be applied using $EFIELD, if you wish to develop
custom usages, but $EFIELD input must not be given for any
kind of run discussed below.
A general approach to computing static polarizabilities is
numerical differentiation, namely RUNTYP=FFIELD, which
should work for any energy method provided by GAMESS. A
sequence of computations with fields applied in the x, y,
and/or z directions will generate the three alpha, beta,
and gamma tensors. See $FFCALC for details. Analytic
computation of all three tensors is available for closed
shells only, see RUNTYP=TDHF and $TDHF input, or TDDFT=HPOL
and $TDDFT input. If you need to know just the static
alpha polarizability, see POLAR in $CPHF during any
analytic hessian job.
A break down of the static alpha polarizability in terms of
contributions from individual localized orbitals can be
obtained by setting POLDCM=.TRUE. in $LOCAL. Calculation
will be by analytic means, unless POLNUM in that group is
selected. This option is available only for SCFTYP=RHF.
The keyword LOCHYP in $FFCALC gives a similar analysis for
all three static polarizabilities, determined by numerical
differentiation.
Polarizabilities in a static electric field differ from
those in an oscillating field, such as a laser produces.
These are called frequency dependent alpha, beta, or gamma,
and in the limit of entering a zero frequency, become the
static quantities discussed just above.
For RHF cases, various frequency dependent alpha, beta, and
gamma polarizabilities can be generated, depending on the
experiment. A particularly easy one to understand is
'second harmonic generation', governed by a beta tensor
describing the absorption of two photons with the emission
of one photon at doubled frequency. See RUNTYP=TDHF, and
papers listed under $TDHF, for many other non-linear
optical experiments. A program for the computation of the
frequency dependent beta hyperpolarizability at the DFT
level is also available, for closed shell molecules: see
TDDFT=HPOL and keywords in $TDDFT input.
Nuclear derivatives of the dipole moment and the various
polarizabilities are also of interest. For example,
knowledge of the derivative of the dipole with respect to
nuclear coordinates yields the IR intensity. Similarly,
the nuclear derivative of the static alpha polarizability
gives Raman activities: see RUNTYP=RAMAN. Analytically
computed 1st or 2nd nuclear derivatives of static or
frequency dependent polarizabilities are available for
SCFTYP=RHF, see RUNTYP=TDHFX and $TDHFX, giving rise to
experimental observations such as resonance Raman and
hyper-Raman.
Finally, instead of considering polarizabilities to be a
function of real frequencies, they can be considered to be
dependent on the imaginary frequency. The imaginary
frequency dependent alpha polarizability can be computed
analytically for SCFTYP=RHF only, using POLDYN=.TRUE. in
$LOCAL. Integration of this quantity over the imaginary
frequency domain can be used to extract C6 dispersion
constants.
Polarizabilities are tensor quantities. There are a number
of different ways to define them, and various formulae to
extract "scalar" and "vector" quantites from the tensors.
A good reference for learning how to compare the output of
a theoretical program to experiment is
A.Willetts, J.E.Rice, D.M.Burland, D.P.Shelton
J.Chem.Phys. 97, 7590-7599(1992)
Localized Molecular Orbitals
Three different orbital localization methods are
implemented in GAMESS. The energy and dipole based
methods normally produce similar results, but see
M.W.Schmidt, S.Yabushita, M.S.Gordon in J.Chem.Phys.,
1984, 88, 382-389 for an interesting exception. You can
find references to the three methods at the beginning of
this chapter.
The method due to Edmiston and Ruedenberg works by
maximizing the sum of the orbitals' two electron self
repulsion integrals. Most people who think about the
different localization criteria end up concluding that
this one seems superior. The method requires the two
electron integrals, transformed into the molecular orbital
basis. Because only the integrals involving the orbitals
to be localized are needed, the integral transformation is
actually not very time consuming.
The Boys method maximizes the sum of the distances
between the orbital centroids, that is the difference in
the orbital dipole moments.
The population method due to Pipek and Mezey maximizes
a certain sum of gross atomic Mulliken populations. This
procedure will not mix sigma and pi bonds, so you will not
get localized banana bonds. Hence it is rather easy to
find cases where this method give different results than
the Ruedenberg or Boys approach.
GAMESS will localize orbitals for any kind of RHF, UHF,
ROHF, or MCSCF wavefunctions. The localizations will
automatically restrict any rotation that would cause the
energy of the wavefunction to be changed (the total
wavefunction is left invariant). As discussed below,
localizations for GVB or CI functions are not permitted.
The default is to freeze core orbitals. The localized
valence orbitals are scarcely changed if the core orbitals
are included, and it is usually convenient to leave them
out. Therefore, the default localizations are: RHF
functions localize all doubly occupied valence orbitals.
UHF functions localize all valence alpha, and then all
valence beta orbitals. ROHF functions localize all valence
doubly occupied orbitals, and all singly occupied orbitals,
but do not mix these two orbital spaces. MCSCF functions
localize all valence MCC type orbitals, and localize all
active orbitals, but do not mix these two orbital spaces.
To recover the invariant MCSCF function, you must be using
a FORS=.TRUE. wavefunction, and you must set GROUP=C1 in
$DRT, since the localized orbitals possess no symmetry.
In general, GVB functions are invariant only to
localizations of the NCO doubly occupied orbitals. Any
pairs must be written in natural form, so pair orbitals
cannot be localized. The open shells may be degenerate, so
in general these should not be mixed. If for some reason
you feel you must localize the doubly occupied space, do a
RUNTYP=PROP job. Feed in the GVB orbitals, but tell the
program it is SCFTYP=RHF, and enter a negative ICHARG so
that GAMESS thinks all orbitals occupied in the GVB are
occupied in this fictitous RHF. Use NINA or NOUTA to
localize the desired doubly occupied orbitals. Orbital
localization is not permitted for CI, because we cannot
imagine why you would want to do that anyway.
Boys localization of the core orbitals in molecules
having elements from the third or higher row almost never
succeeds. Boys localization including the core for second
row atoms will often work, since there is only one inner
shell on these. The Ruedenberg method should work for any
element, although including core orbitals in the integral
transformation is more expensive.
The easiest way to do localization is in the run which
generates the wavefunction, by selecting LOCAL=xxx in the
$CONTRL group. However, localization may be conveniently
done at any time after determination of the wavefunction,
by executing a RUNTYP=PROP job. This will require only
$CONTRL, $BASIS/$DATA, $GUESS (pick MOREAD), the converged
$VEC, possibly $SCF or $DRT to define your wavefunction,
and optionally some $LOCAL input.
There is an option to restrict all rotations that would
mix orbitals of different symmetries. SYMLOC=.TRUE. yields
only partially localized orbitals, but these still possess
symmetry. They are therefore very useful as starting
orbitals for MCSCF or GVB-PP calculations. Because they
still have symmetry, these partially localized orbitals run
as efficiently as the canonical orbitals. Because it is
much easier for a user to pick out the bonds which are to
be correlated, a significant number of iterations can be
saved, and convergence to false solutions is less likely.
* * *
The most important reason for localizing orbitals is to
analyze the wavefunction. A simple example is to look at
shapes of the orbitals with the MacMolPlt program. Or, you
might read the localized orbitals in during a RUNTYP=PROP
job to examine their Mulliken populations.
Localized orbitals are a particularly interesting way
to analyze MCSCF computations. The localized orbitals may
be oriented on each atom (see option ORIENT in $LOCAL) to
direct the orbitals on each atom towards their neighbors
for maximal bonding, and then print a bond order analysis.
The orientation procedure is newly programmed by J.Ivanic
and K.Ruedenberg, to deal with the situation of more than
one localized orbital occuring on any given atom. Some
examples of this type of analysis are
D.F.Feller, M.W.Schmidt, K.Ruedenberg
J.Am.Chem.Soc. 104, 960-967 (1982)
T.R.Cundari, M.S.Gordon
J.Am.Chem.Soc. 113, 5231-5243 (1991)
N.Matsunaga, T.R.Cundari, M.W.Schmidt, M.S.Gordon
Theoret.Chim.Acta 83, 57-68 (1992).
In addition, the energy of your molecule can be
partitioned over the localized orbitals so that you may
be able to understand the origin of barriers, etc. This
analysis can be made for the SCF energy, and also the MP2
correction to the SCF energy, which requires two separate
runs. An explanation of the method, and application to
hydrogen bonding may be found in J.H.Jensen, M.S.Gordon,
J.Phys.Chem. 99, 8091-8107(1995).
Analysis of the SCF energy is based on the localized
charge distribution (LCD) model: W.England and M.S.Gordon,
J.Am.Chem.Soc. 93, 4649-4657 (1971). This is implemented
for RHF and ROHF wavefunctions, and it requires use of
the Ruedenberg localization method, since it needs the
two electron integrals to correctly compute energy sums.
All orbitals must be included in the localization, even
the cores, so that the total energy is reproduced.
The LCD requires both electronic and nuclear charges
to be partitioned. The orbital localization automatically
accomplishes the former, but division of the nuclear
charge may require some assistance from you. The program
attempts to correctly partition the nuclear charge, if you
select the MOIDON option, according to the following: a
Mulliken type analysis of the localized orbitals is made.
This determines if an orbital is a core, lone pair, or
bonding MO. Two protons are assigned to the nucleus to
which any core or lone pair belongs. One proton is
assigned to each of the two nuclei in a bond. When all
localized orbitals have been assigned in this manner, the
total number of protons which have been assigned to each
nucleus should equal the true nuclear charge.
Many interesting systems (three center bonds, back-
bonding, aromatic delocalization, and all charged species)
may require you to assist the automatic assignment of
nuclear charge. First, note that MOIDON reorders the
localized orbitals into a consistent order: first comes
any core and lone pair orbitals on the 1st atom, then
any bonds from atom 1 to atoms 2, 3, ..., then any core
and lone pairs on atom 2, then any bonds from atom 2 to
3, 4, ..., and so on. Let us consider a simple case
where MOIDON fails, the ion NH4+. Assuming the nitrogen
is the 1st atom, MOIDON generates
NNUCMO=1,2,2,2,2
MOIJ=1,1,1,1,1
2,3,4,5
ZIJ=2.0,1.0,1.0,1.0,1.0,
1.0,1.0,1.0,1.0
The columns (which are LMOs) are allowed to span up to 5
rows (the nuclei), in situations with multicenter bonds.
MOIJ shows the Mulliken analysis thinks there are four
NH bonds following the nitrogen core. ZIJ shows that
since each such bond assigns one proton to nitrogen, the
total charge of N is +6. This is incorrect of course,
as indeed will always happen to some nucleus in a charged
molecule. In order for the energy analysis to correctly
reproduce the total energy, we must ensure that the
charge of nitrogen is +7. The least arbitrary way to
do this is to increase the nitrogen charge assigned to
each NH bond by 1/4. Since in our case NNUCMO and MOIJ
and much of ZIJ are correct, we need only override a
small part of them with $LOCAL input:
IJMO(1)=1,2, 1,3, 1,4, 1,5
ZIJ(1)=1.25, 1.25, 1.25, 1.25
which changes the charge of the first atom of orbitals
2 through 5 to 5/4, changing ZIJ to
ZIJ=2.0,1.25,1.25,1.25,1.25,
1.0, 1.0, 1.0, 1.0
The purpose of the IJMO sparse matrix pointer is to let
you give only the changed parts of ZIJ and/or MOIJ.
Another way to resolve the problem with NH4+ is to
change one of the 4 equivalent bond pairs into a "proton".
A "proton" orbital AH treats the LMO as if it were a
lone pair on A, and so assigns +2 to nucleus A. Use of
a "proton" also generates an imaginary orbital, with
zero electron occupancy. For example, if we make atom
2 in NH4+ a "proton", by
IPROT(1)=2
NNUCMO(2)=1
IJMO(1)=1,2,2,2 MOIJ(1)=1,0 ZIJ(1)=2.0,0.0
the automatic decomposition of the nuclear charges will be
NNUCMO=1,1,2,2,2,1
MOIJ=1,1,1,1,1,2
3,4,5
ZIJ=2.0,2.0,1.0,1.0,1.0,1.0
1.0,1.0,1.0
The 6th column is just a proton, and the decomposition
will not give any electronic energy associated with
this "orbital", since it is vacant. Note that the two ways
we have disected the nuclear charges for NH4+ will both
yield the correct total energy, but will give very
different individual orbital components. Most people
will feel that the first assignment is the least arbitrary,
since it treats all four NH bonds equivalently.
However you assign the nuclear charges, you must
ensure that the sum of all nuclear charges is correct.
This is most easily verified by checking that the energy
sum equals the total SCF energy of your system.
As another example, H3PO is studied in EXAM26.INP.
Here the MOIDON analysis decides the three equivalent
orbitals on oxygen are O lone pairs, assigning +2 to
the oxygen nucleus for each orbital. This gives Z(O)=9,
and Z(P)=14. The least arbitrary way to reduce Z(O)
and increase Z(P) is to recognize that there is some
backbonding in these "lone pairs" to P, and instead
assign the nuclear charge of these three orbitals by
1/3 to P, 5/3 to O.
Because you may have to make several runs, looking
carefully at the localized orbital output before the
correct nuclear assignments are made, there is an
option to skip directly to the decomposition when the
orbital localization has already been done. Use
$CONTRL RUNTYP=PROP
$GUESS GUESS=MOREAD NORB=
$VEC containing the localized orbitals!
$TWOEI
The latter group contains the necessary Coulomb and
exchange integrals, which are punched by the first
localization, and permits the decomposition to begin
immediately.
SCF level dipoles can also be analyzed using the
DIPDCM flag in $LOCAL. The theory of the dipole
analysis is given in the third paper of the LCD
sequence. The following list includes application of
the LCD analysis to many problems of chemical interest:
W.England, M.S.Gordon J.Am.Chem.Soc. 93, 4649-4657 (1971)
W.England, M.S.Gordon J.Am.Chem.Soc. 94, 4818-4823 (1972)
M.S.Gordon, W.England J.Am.Chem.Soc. 94, 5168-5178 (1972)
M.S.Gordon, W.England Chem.Phys.Lett. 15, 59-64 (1972)
M.S.Gordon, W.England J.Am.Chem.Soc. 95, 1753-1760 (1973)
M.S.Gordon J.Mol.Struct. 23, 399 (1974)
W.England, M.S.Gordon, K.Ruedenberg,
Theoret.Chim.Acta 37, 177-216 (1975)
J.H.Jensen, M.S.Gordon, J.Phys.Chem. 99, 8091-8107(1995)
J.H.Jensen, M.S.Gordon, J.Am.Chem.Soc. 117, 8159-8170(1995)
M.S.Gordon, J.H.Jensen, Acc.Chem.Res. 29, 536-543(1996)
* * *
It is also possible to analyze the MP2 correlation
correction in terms of localized orbitals, for the RHF
case. The method is that of G.Peterssen and M.L.Al-Laham,
J.Chem.Phys., 94, 6081-6090 (1991). Any type of localized
orbital may be used, and because the MP2 calculation
typically omits cores, the $LOCAL group will normally
include only valence orbitals in the localization. As
mentioned already, the analysis of the MP2 correction must
be done in a separate run from the SCF analysis, which must
include cores in order to sum up to the total SCF energy.
* * *
Typically, the results are most easily interpreted
by looking at "the bigger picture" than at "the details".
Plots of kinetic and potential energy, normally as a
function of some coordinate such as distance along an
IRC, are the most revealing. Once you determine, for
example, that the most significant contribution to the
total energy is the kinetic energy, you may wish to look
further into the minutia, such as the kinetic energies
of individual localized orbitals, or groups of LMOs
corresponding to an entire functional group.
Transition Moments and Spin-Orbit Coupling
A review of various ways of computing spin-orbit coupling:
D.G.Fedorov, S.Koseki, M.W.Schmidt, M.S.Gordon,
Int.Rev.Phys.Chem. 22, 551-592(2003)
GAMESS can compute transition moments and oscillator
strengths for the radiative transitions between states
written in terms of CI wavefunctions (GUGA only). The
transition moments are computed using both the "length
form" and the "velocity form". In a.u., where h-bar=m=1, we
start from
[A,q] = -i dA/dp
For A=H, dH/dp=p, and p= -i d/dq,
[H,q] = -i p = -d/dq.
For non-degenerate states,
=
(Ea-Eb) = -
This relates the dipole to the velocity form,
= -1/(Ea-Eb)
but the CI states will give different numbers for each
side, since the states aren't exact eigenfunctions.
Transition moment computation is OPERAT=DM in $TRANST. For
transition moments, the CI is necessarily performed on
states of the same multiplicity.
All other operators are various spin-orbit coupling
options. There are two kinds of calculations possible,
which we will call SO-CI and SO-MCQDPT. Note that there
is a hyphen in "spin-orbit CI" to avoid confusion with
"second order CI" in the sense of the SOCI keyword in $DRT
input. For SO-CI, the initial states may be any CI wave-
function that the GUGA package can generate. For SO-MCQDPT
the initial states for spin-orbit coupling are of CAS type,
and the operator mixing them corresponds to MCQDPT
generalised for spin-dependent operators (with certain
approximations).
GAMESS can compute the "microscopic Breit-Pauli
spin-orbit operator", which includes both a one and two
electron operator. Additional information is given in a
subsection below
states
For transition moments, the states are generated by
CI calculations using the GUGA package. These states are
the final states, and the results are just the transition
moments between these states. The states are defined by
$DRTx input groups.
For SO-CI, the energy of the CI states forms the
diagonal of a spin-orbit Hamiltonian, as in the state basis
the spin-free Hamiltonian is of course diagonal. Addition
of the Pauli-Breit operator does not change the diagonal,
but does add off-diagonal H-so elements. For SO-MCQDPT,
the spin-free MCQDPT matrix elements are expanded into
matrices corresponding to all Ms values for a pair of
multiplicities. These matrices are block-diagonal before
the addition of spin-orbit coupling terms, coupling Ms
values. The diagonalization of this spin-orbit Hamiltonian
gives new energy levels, and spin-mixed final states.
Optionally, the transition dipoles between the final states
can be computed. The input requirements are $DRTx or
$MCQDx groups which define the original pure spin states.
We will call the initial states CAS-CI, since most of
the time they will be MCSCF states. There may be cases
such as the Na example below where SCF orbitals are used,
or other cases where a FOCI or SOCI wavefunction might be
used for the initial states. Please keep in mind that the
term does not imply the states must be MCSCF states, just
that they commonly are.
In the above, x may vary from 1 to 64. The reason for
allowing such a large range is to permit the use of Abelian
point group symmetry during the generation of the initial
states. The best explanation will be an example, but the
number of these input groups depends on both the number of
orbital sets input, and how much symmetry is present. The
next two subsections discuss these points.
orbitals
The orbitals for transition moments or for SO-CI can be
one common set of orbitals used by all CI states. If one
set of orbitals is used, the transition moment or spin-
orbit coupling can be found for any type of GUGA CI wave-
function. Alternatively, two sets of orbitals (obtained by
separate MCSCF orbital optimizations) can be used. Two or
more separate CIs will be carried out. The two MO sets
must share a common set of frozen core orbitals, and the
CI must be of the complete active space type. These
restrictions are needed to leave the CI wavefunctions
invariant under the necessary rotation to corresponding
orbitals. The non-orthogonal procedure implemented is a
GUGA driven equivalent to the method of Lengsfield, et al.
Note that the FOCI and SOCI methods described by these
workers are not available in GAMESS.
If you would like to use separate orbitals during the
CI, they may be generated with the FCORE option in $MCSCF.
Typically you would optimize the ground state completely,
then use these MCSCF orbitals in an optimization of the
excited state, under the constraint of FCORE=.TRUE.
For SO-MCQDPT calculations, only one set of orbitals
may be input to describe all CAS-CI states. Typically that
orbital set will be obtained by state-averaged MCSCF, see
WSTATE in $DET/$DRT, and also in the $MCQDx input. Note
that although the RUNTYP=TRANSITN driver is tied to the
GUGA CI package, there is no reason the orbitals cannot be
obtained using the determinant CI package. In fact, for
the case of spin-orbit coupling, you might want to utilize
the ability to state average over several spins, see PURES
in $DET.
If there is no molecular symmetry present, transition
moment calculations will provide $DRT1 if there is one set
of orbitals, otherwise $DRT1 defines the CI based on $VEC1
and $DRT2 the CI based on $VEC2. Also for the case of no
symmetry, a spin-orbit job should enter one $DRTx or $MCQDx
for every spin multiplicity, and all states of the same
multiplicity have to be generated from $VEC1 or $VEC2,
according to IVEX input.
symmetry
The CAS-CI states are most efficiently generated using
symmetry, since states of different symmetry have zero
Hamiltonian matrix elements. It is probably more efficient
to do four CI calculations in the group C2v on A1, A2, B1,
and B2 symmetry, than one CI with a combined Hamiltonian
in C1 symmetry (unless the active space is very small), and
similar remarks apply to the SO-MCQDPT case. In order to
avoid repeatedly saying $DRTx or $MCQDx, the following few
paragraphs say $DRTx only.
Again supposing the group is C2v, and you are
interested in singlet-triplet coupling. After some
preliminary CI calculations, you might know that the lowest
8 states are two 1-a1, 1-b1, two 1-b2, one 3-a1, and two 3-
b2 states. In this case your input would consist of five
$DRTx, of which you can give the three singlets in any
order but these must preceed the two triplet input groups
to follow the rule of increasing multiplicity. Clearly it
is not possible to write a formula for how many $DRTx there
will be, this depends not only on the point group, but also
the chemistry of the situation.
If you are using two sets of orbitals, the generation
of the corresponding orbitals for the two sets will permute
the active orbitals in an unpredictable way. Use STSYM to
define the desired state symmetry, rather than relying on
the orbital order. It is easy and safer to be explicit
about the spatial orbital symmetry.
The users are encouraged to specify full symmetry in
their $DATA input even though they may choose to set the
symmetry in $DRTx to C1. The CI states will be labelled in
the group given in $DATA. The use of non-Abelian symmetry
is limited by the absence of non-Abelian CI or MCQDPT. In
this case the users can choose between setting full non-
Abelian symmetry in $DATA and C1 in $DRT or else an Abelian
subgroup in both $DATA and $DRT. The latter choice appears
to be most efficient at present.
An example of SO-MCQDPT illustrating how the carbon
atom of Kh symmetry (full rotation-reflection group) can be
entered in D2h, Kh's highest Abelian group. The run time is
considerably longer in C1 symmetry.
As another example, consider an organic molecule with a
singly excited state, where that state might be coupled to
low or high spin, and where these two states might be close
enough to have a strong spin-orbit coupling. If it happens
that the S1 and S0 states possess different symmetry, a
very reaasonable calculation would be to treat the S1 and
T1 state with the same $VEC2 orbitals, leaving the ground
state described by $VEC1. After doing an MCSCF on the S0
ground state for $VEC1, you could do a state-averaged MCSCF
for $VEC2 optimized for T1 and S1 simultaneously, using
PURES. The spin orbit job would obtain its initial states
from three GUGA CI computations, S0 from $VEC1 and $DRT1,
S1 from $VEC2 and $DRT2, and T1 from $VEC2 and $DRT3. Your
$TRANST would be NUMCI=3, IROOTS(1)=1,1,1, IVEX(1)=1,2,2.
Note that the second IROOTS value is 1 because S1 was
presumed to have a different symmetry than S0, so STSYM in
$DRT1 and $DRT2 will differ. The calculation just outlined
cannot be done if S0 and S1 have the same spatial symmetry,
as IROOTS(1)=1,2,1 to obtain S1 during the second CI will
bring in an additional S0 state (one expressed in terms of
the $VEC2, at slightly higher energy). This problem is the
origin of the statement several paragraphs above that a
system with no symmetry will have one $DRTx for every spin
multiplicity included.
For transition moments, which do not diagonalize a
matrix containing these duplicated states, it is OK to
proceed, provided you ignore all transition moments between
the same states obtained in the two different CIs.
spin orbit coupling
Spin-orbit coupling calculations are always performed
in a quasi-degenerate perturbative manner. Typically the
states close in energy are included into the spin-orbit
coupling matrix. "Close" has a easily understandable
meaning, since in the limit of small coupling the quasi-
degenerate treatment is reduced to a second order
perturbative treatment, that is, the affect of a state upon
the state of primary interest is given by the square of the
spin-orbit coupling matrix element divided by the
difference of the adiabatic energies. This is useful to
keep in mind when deciding how many CI states to include in
the matrix. The states that are included are treated in a
fashion that is equivalent to infinite order perturbation
theory (exact) whereas the states that are not included
make no contribution.
Spin-orbit runs can be done for even or odd numbers of
electrons (any spin), for more than two different spin
multiplicities at once, for general active spaces. At
times, when the spatial wavefunction is degenerate, a spin-
orbit run might involve only one spin multiplicity, e.g. a
triplet-pi state in a linear molecule. The most common
case is two different spins, with non-zero spin orbit
coupling possible only for delta-S=1: singlets spin-orbit
couple with triplets, but not with quintets. Use of three
spins, such as S=0,1,2, will generate couplings between
singlets and triplets, and between triplets and quintets,
which together engender an indirect singlet/quintet mixing.
As noted above, the treatment of spin-orbit involves
first obtaining a handful of spin-pure states, whose
energies form the diagonal of a model Hamiltonian. The
spin-orbit operator introduces off-diagonal couplings, and
the resulting small Hamiltonian is diagonalized. This
generates spin-mixed states in the model space. Since the
model states are not fully relaxed (internally contracted),
this is essentially a perturbative treatment: certainly
the spin-orbit effects have no influence on orbital
optimizations or potential energy surfaces when treated in
this manner, at the very last stage.
The Breit-Pauli spin-orbit operator contains a one
electron term arising from Pauli's reduction of the Dirac
hydrogenic equation to a single-component form, and a two
electron term due to Breit. Computation of the full Breit-
Pauli operator is OPERAT=HSO2 (or HSOFF). A close
approximation to the latter is HSO2P (P=partial), which
neglects all active-active two electron terms, which
usually do not contribute very much to the total coupling,
while saving substantial computer time.
HSO1 completely omits the two electron terms, so is
much faster than any of the two electron operators, but
represents a potentially much greater loss of accuracy.
HSO1's error can be remedied to some extent by regarding
the nuclear charge in the one electron term as an
adjustable parameter. In addition, these effective charges
are often used to compensate for missing nodes in valence
orbitals of ECP runs, in which case the ZEFF are typically
very far from the two nuclear charges. ZEFTYP selects some
built in values obtained by S.Koseki et al, but if you have
some favorite parameters, they can be read in as the ZEFF
input array. Effective charges may be used for any OPERAT,
but are most often used with HSO1.
Theoretical considerations indicate that the Breit-
Pauli operator is not variationally stable, if it were to
be used during the SCF iterations determining orbitals.
However, a first order Douglas-Kroll type correction to the
one electron part of the operator reduces its size by means
of certain kinematic factors, removing this problem. This
can produce substantially better results, even if the
operator is being treated pertubatively. This correction
(at first order) is automatically applied to the one
electron part of the Breit-Pauli operator by any run that
selects RELWFN=DK (at any order) in the scalar relativity
treatment during the variational steps prior to the spin-
orbit perturbation. See paper 32 below.
Because the diagonalization of the model spin-orbit
Hamiltonian leads to spin-mixed eigenvectors, approximate
wavefunctions including spin-orbit coupling are generated.
It is now possible to generate the density matrix for the
spin-mixed states, so that property values for spin-mixed
states can be found: see keyword ISTNO. The natural
orbitals of these spin-orbit density matrices turn out to
be good approximations to the two spinors of the large
components of full Dirac four component runs. The total
density of these spinors can be obtained for interpretation
purposes. Recognition that the spin-orbit coupling is a
rotational operator (L dot S, where L = R cross P) when it
acts on an orbital can lead to chemical interpretability of
the spin orbit results. See papers 40 and 41 below.
It is also possible to obtain the dipole transition
moments between the final spin-mixed wavefunctions, which
of course do not any longer have a rigourous S quantum no.
When the run is SO-MCQDPT, the transition moment are first
computed only between CAS states, and then combined with
the spin-mixed SO-MCQDPT coefficients.
technical matters:
The only practical limitation on the computation of
the Breit term is that HSO2FF is limited to 10 active
orbitals on 32 bit machines, and to about 26 active
orbitals on 64 bit machines. The spin-orbit matrix
elements vanish for |delta-S| > 1, but it is possible to
include three or more spins in the computation. Since
singlets interact with triplets, and triplets interact with
pentuplets, inclusion of S=0,1,2 simultaneously lets you
pick up the indirect interaction between singlets and
pentuplets that the intermediate triplets afford.
The choice between HSO2 and HSO2FF is very often in
favor of the former. HSO2 computes the matrix elements in
CSF basis and then contracts them with CI coefficients,
whereas HSO2FF uses a generalized density in AO basis
computed for each pair of states, thus HSO2 is much more
efficient in case of multiple states given in IROOTS.
HSO2FF takes less memory for integral storage, thus it can
be superior in case of small active spaces and large basis
sets, in part because it does not store 2e SOC integrals on
disk and secondly, it does not redundantly treat the same
pair of determinants if they appear in different CSFs. The
numerical results with HSO2 and HSO2FF should be identical
within machine and algorithmic accuracy.
Various symmetries are used to avoid computing zero
spin-orbit matrix elements. NOSYM in $TRANST allows some
control over this: NOSYM=1 gives up point group symmetry
completely, while 2 turns off additional symmetries such
as spin selection rules. HSO1,2,2P compute all matrix
elements in a group (i.e. between two $DRTx groups with
fixed Ms(ket)-Ms(bra)) if at least one of them does not
vanish by symmetry, and HSO2P actually avoids computation
for each pair of states if forbidden by symmetry. Setting
NOSYM=2 will cause HSO2FF to consider the elements between
two singlets, which are always calculated for HSO1,2,2P
when transition dipoles are requested as well.
SYMTOL has a dramatic effect on the run speed. This
cutoff is applied to CSF coefficcients, their products,
and these products times CSF orbital overlaps. The value
permits a tradeoff of accuracy for run time, and since the
error in the spin-orbit properties approaches SYMTOL mainly
for SOCI functions, it may be useful to increase SYMTOL to
save time for CAS or FOCI functions. Some experimenting
will tell you what you can get away with. SYMTOL is also
used during CI state symmetry assignment, for NOIRR=-1
in $DRT.
In case if you do not provide enough storage for the
form factors sorting then some extra disk space will be
used; the extra disk space can be eliminated if you set
SAVDSK=.TRUE. (the amount of savings depends on the active
space and memory provided, it some cases it can decrease
the disk space up to one order of magnitude). The form
factors are in binary format, and so can be transfered
between computers only if they have compatible binary
files. There is a built-in check for consistency of a
restart file DAFL30 with the current run parameters.
input nitty-gritty
The transition moment and spin-orbit coupling driver
is a rather restricted path through the GUGA CI part of
GAMESS. Note that $GUESS is not read, instead the MOs will
be MOREAD in a $VEC1 and perhaps a $VEC2 group. It is not
possible to reorder MOs. For SO-CI,
1) Give SCFTYP=NONE CITYP=GUGA MPLEVL=0.
2) $CIINP is not read. The CI is hardwired to consist
of CI DRT generation, integral transformation/sorting,
Hamiltonian generation, and diagonalization. This
means $DRT1 (and maybe $DRT2,...), $TRANS, $CISORT,
$GUGEM, and $GUGDIA input is read, and acted upon.
3) The density matrices are not generated, and so no
properties (other than the transition moment or the
spin-orbit coupling) are computed.
4) There is no restart capability provided, except for
saving some form-factor information.
5) $DRT1, $DRT2, $DRT3, ... must go from lowest to highest
multiplicity.
6) IROOTS will determine the number of CI states in each
CI for which the properties are calculated. Use
NSTATE to specify the number of CI states for the
CI Hamiltonian diagonalization. Sometimes the CI
convergence is assisted by requesting more roots
to be found in the diagonalization than you want to
include in the property calculation.
For SO-MCQDPT, the steps are
1) Give SCFTYP=NONE CITYP=NONE MPLEVL=2.
2) the number of roots in each MCQDPT is controlled by
$TRANST's IROOTS, and each such calculation is
defined by $MCQD1, $MCQD2, ... input. These must go
from lowest multiplicity to highest.
references
The review already mentioned:
"Spin-orbit coupling in molecules: chemistry beyond the
adiabatic approximation".
D.G.Fedorov, S.Koseki, M.W.Schmidt, M.S.Gordon,
Int.Rev.Phys.Chem. 22, 551-592(2003)
Reference for separate active orbital optimization:
1. B.H.Lengsfield, III, J.A.Jafri, D.H.Phillips,
C.W.Bauschlicher, Jr. J.Chem.Phys. 74,6849-6856(1981)
References for transition moments:
2a. H.C.Longuet-Higgins
Proc.Roy.Soc.(London) A235, 537-543(1956)
2b. F.Weinhold, J.Chem.Phys. 54,1874-1881(1970)
3. C.W.Bauschlicher, S.R.Langhoff
Theoret.Chim.Acta 79:93-103(1991)
4. "Intermediate Quantum Mechanics, 3rd Ed." Hans A.
Bethe, Roman Jackiw Benjamin/Cummings Publishing,
Menlo Park, CA (1986), chapters 10 and 11.
5. S.Koseki, M.S.Gordon
J.Mol.Spectrosc. 123, 392-404(1987)
References for HSO1 spin-orbit coupling, and Zeff values:
6. S.Koseki, M.W.Schmidt, M.S.Gordon
J.Phys.Chem. 96, 10768-10772 (1992)
7. S.Koseki, M.S.Gordon, M.W.Schmidt, N.Matsunaga
J.Phys.Chem. 99, 12764-12772 (1995)
8. N.Matsunaga, S.Koseki, M.S.Gordon
J.Chem.Phys. 104, 7988-7996 (1996)
9. S.Koseki, M.W.Schmidt, M.S.Gordon
J.Phys.Chem.A 102, 10430-10435 (1998)
10. S.Koseki, D.G.Fedorov, M.W.Schmidt, M.S.Gordon
J.Phys.Chem.A 105, 8262-8268 (2001)
11. S.Koseki, T.Matsushita, M.S.Gordon
J.Phys.Chem.A 110, 2560-2570(2006)
References for full Breit-Pauli spin-orbit coupling:
20. T.R.Furlani, H.F.King
J.Chem.Phys. 82, 5577-5583 (1985)
21. H.F.King, T.R.Furlani
J.Comput.Chem. 9, 771-778 (1988)
22. D.G.Fedorov, M.S.Gordon
J.Chem.Phys. 112, 5611-5623 (2000)
Paper 22 contains information on the HSO2P partial two
electron operator method.
Symmetry in spin-orbit coupling:
23. D.G.Fedorov, M.S.Gordon
ACS Symposium Series 828, pp 1-22(2002)
Reference for SO-MCQDPT:
25. D.G.Fedorov, J.P.Finley Phys.Rev.A 64, 042502 (2001)
Reference for Spin-Orbit with Model Core Potentials:
30. D.G.Fedorov, M.Klobukowski
Chem.Phys.Lett. 360, 223-228(2002)
31. T.Zeng, D.G.Fedorov, M.Klobukowski
J.Chem.Phys. 131, 124109/1-17(2009)
32. T.Zeng, D.G.Fedorov, M.Klobukowski
J.Chem.Phys. 132, 074102/1-15(2010)
The last two of these also discuss the 1st order Douglas-
Kroll transformation of the 1e- part of the spin orbit
operator.
Reference for properties, interpretations, and Spin-Orbit
Natural Spinors:
40. T. Zeng, D. G. Fedorov, M. W. Schmidt, M. Klobukowski
J. Chem. Phys. 134, 214107/1-9(2011)
41. T. Zeng, D. G. Fedorov, M. W. Schmidt, M. Klobukowski
J. Chem. Theor. Comput. 7, 2864-2875(2011)
with an application being
42. T.Zeng, D.G.Federov, M.W.Schmidt, M.Klobukowski
J.Chem.Theo.Comput. 8, 3061-3071(2012).
Recent applications (see also 32,40,41):
50. S.P.Webb, M.S.Gordon
J.Chem.Phys. 109, 919-927(1998)
51. D.G.Fedorov, M.Evans, Y.Song, M.S.Gordon, C.Y.Ng
J.Chem.Phys. 111, 6413-6421 (1999)
52. D.G.Fedorov, M.S.Gordon, Y.Song, C.Y.Ng
J.Chem.Phys. 115, 7393-7400 (2001)
53. B.J.Duke J.Comput.Chem. 22, 1552-1556 (2001)
54. C.M.Aikens, M.S.Gordon
J.Phys.Chem.A 107, 104-114(2003)
55. D.G.Fedorov, S.Koseki, M.W.Schmidt, M.S.Gordon
K.Hirao and Y.Ishikawa (eds.) Recent Advances in
Relativistic Molecular Theory, Vol. 5,
(World Scientific, Singapore), 2004, pp 107-136.
* * *
Special thanks to Bob Cave and Dave Feller for their
assistance in performing check spin-orbit coupling runs
with the MELDF programs. Special thanks to Tom Furlani for
contributing his 2e- spin-orbit code and answering many
questions about its interface. Special thanks to Haruyuki
Nakano for explaining the spin functions used in the MCQDPT
package.
examples
We end with 2 examples. Note that you must know what
you are doing with term symbols, J quantum numbers, point
group symmetry, and so on in order to make skillful use of
this part of the program. Seeing your final degeneracies
turn out like a text book says it should is beautiful!
! Compute the splitting of the famous sodium D line.
! Joseph von Fraunhofer (Denkschriften der Koeniglichen
! Akademie der Wissenschf. zu Muenchen, 5, 193(1814-1815))
! observed the sun through good prisms, finding 700 lines,
! and named the brightest ones A, B, C... just in order.
! He was able to resolve the D line into two lines, which
! occur at 5895.940 and 5889.973 Angstroms. It would take
! a century to understand the D line is Na's 3s <-> 3p
! transition, and that spin-orbit coupling is what splits
! the D line into two. Charlotte Moore's Atomic Energy
! Levels, volume 1, gives the experimental 2-P interval
! as 17.1963, since the three relevent levels are at
! 2-S-1/2= 0.0, 2-P-1/2= 16,956.183, 2-P-3/2= 16,973.379.
1. generate ground state 2-S orbitals by conventional ROHF.
the energy of the ground state is -161.8413919816
--- $contrl scftyp=rohf mult=2 $end
--- $system kdiag=3 memory=300000 $end
--- $guess guess=huckel $end
2. generate excited state 2-P orbitals, using a state-
averaged SCF wavefunction to ensure radial degeneracy of
the 3p shell is preserved. The open shell SCF energy is
-161.7682895801. The computation is both spin and space
restricted open shell SCF on the 2-P Russell-Saunders term.
Starting orbitals are reordered orbitals from step 1.
--- $contrl scftyp=gvb mult=2 $end
--- $system kdiag=3 memory=300000 $end
--- $guess guess=moread norb=13
--- norder=1 iorder(6)=7,8,9,6 $end
--- $scf nco=5 nseto=1 no(1)=3 rstrct=.t. couple=.true.
--- f(1)= 1.0 0.16666666666667
--- alpha(1)= 2.0 0.33333333333333 0.0
--- beta(1)= -1.0 -0.16666666666667 0.0 $end
3. compute spin-orbit coupling in the 2-P term. The use of
C1 symmetry in $DRT1 ensures that all three spatial CSFs
are kept in the CI function. In the preliminary CI, the
spin function is just the alpha spin doublet, and all three
roots should be degenerate, and furthermore equal to the
GVB energy at step 2. The spin-orbit coupling code uses
both doublet spin functions with each of the three spatial
wavefunctions, so the spin-orbit Hamiltonian is a 6x6
matrix. The two lowest roots of the full 6x6 spin-orbit
Hamiltonian are the doubly degenerate 2-P-1/2 level, while
the other four roots are the degenerate 2-P-3/2 level.
$contrl scftyp=none cityp=guga runtyp=transitn mult=2 $end
$system memory=2000000 $end
$basis gbasis=n31 ngauss=6 $end
$gugdia nstate=3 $end
$transt operat=hso1 numvec=1 numci=1 nfzc=5 nocc=8
iroots=3 zeff=10.04 $end
$drt1 group=c1 fors=.true. nfzc=5 nalp=1 nval=2 $end
$data
Na atom...2-P excited state...6-31G basis
Dnh 2
Na 11.0
$end
--- GVB ORBITALS --- GENERATED AT 7:46:08 CST 30-MAY-1996
Na atom...2-P excited state
E(GVB)= -161.7682895801, E(NUC)= .0000000000, 5
ITERS
$VEC1
1 1 9.97912679E-01 8.83038094E-03 0.00000000E+00...
... orbitals from step 2 go here ...
13 3-1.10674398E+00 0.00000000E+00 0.00000000E+00
$END
As an example of both SO-MCQDPT, and the use of as much
symmetry as possible, consider carbon. The CAS-CI uses
an active space of 2s,2p,3s,3p orbitals, and the spin-orbit
job includes all terms from the lowest configuration,
2s2,2p2. These terms are 3-P, 1-D, and 1-S. If you look
at table 58 in Herzberg's book on electronic spectra, you
will be able to see how the Kh spatial irreps P, D, S are
partitioned into the D2h irreps input below.
! C SO-MRMP on all levels in the s**2,p**2 configuration.
!
! levels CAS and MCQDPT
! 1 .0000 .0000 cm-1 3-P-0
! 2-4 12.6879-12.8469 13.2721-13.2722 3-P-1
! 5-9 37.8469-37.8470 39.5638-39.5639 3-P-2
! 10-14 12169.1275 10251.7910 1-D-2
! 15 19264.4221 21111.5130 1-S-0
!
! The active space consists of (2s,2p,3s,3p) with 4 e-.
! D2h symmetry speeds up the calculation considerably,
! on the same computer D2h = 78 and C1 = 424 seconds.
$contrl scftyp=none cityp=none mplevl=2
runtyp=transitn $end
$system memory=5000000 $end
!
! below is input to run in C1 subgroup
!
--- $transt operat=hso2 numvec=-2 numci=2 nfzc=1 nocc=9
--- iroots(1)=6,3 parmp=3
--- ivex(1)=1,1 $end
--- $mrmp mrpt=mcqdpt rdvecs=.t. $end
--- $MCQD1 nosym=1 nstate=6 mult=1 iforb=3
--- nmofzc=1 nmodoc=0 nmoact=8
--- wstate(1)=1,1,1,1,1,1 thrcon=1e-8 thrgen=1e-10 $END
--- $MCQD2 nosym=1 nstate=3 mult=3 iforb=3
--- nmofzc=1 nmodoc=0 nmoact=8
--- wstate(1)=1,1,1 thrcon=1e-8 thrgen=1e-10 $END
!
! below is input to run in D2h subgroup
!
$transt operat=hso2 numvec=-7 numci=7 nfzc=1 nocc=9
iroots(1)=3,1,1,1, 1,1,1 parmp=3
ivex(1)=1,1,1,1,1,1,1 $end
$mrmp mrpt=mcqdpt rdvecs=.t. $end
$MCQD1 nosym=-1 mult=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=Ag wstate(1)=1,1,1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD2 nosym=-1 mult=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=B1g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD3 nosym=-1 mult=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=B2g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD4 nosym=-1 mult=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=B3g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD5 nosym=-1 mult=3 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=B1g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD6 nosym=-1 mult=3 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=B2g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD7 nosym=-1 mult=3 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
stsym=B3g wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
!
! input to prepare the 3-P ground state orbitals
! great care is taken to create symmetry equivalent p's
!
--- $contrl scftyp=mcscf cityp=none mplevl=0
--- runtyp=energy mult=3 $end
--- $guess guess=moread norb=55 purify=.t. $end
--- $mcscf cistep=guga fullnr=.t. $end
--- $drt group=c1 fors=.true.
--- nmcc=1 ndoc=1 nalp=2 nval=5 $end
--- $gugdia nstate=9 maxdia=1000 $end
--- $gugdm2 wstate(1)=1,1,1 $end
!
$data
C...aug-cc-pvtz (10s,5p,2d,1f) -> [4s,3p,2d,1f]
(1s,1p,1d,1f)
Dnh 2
C 6.0
S 8
1 8236.000000 0.5310000000E-03
2 1235.000000 0.4108000000E-02
3 280.8000000 0.2108700000E-01
4 79.27000000 0.8185300000E-01
5 25.59000000 0.2348170000
6 8.997000000 0.4344010000
7 3.319000000 0.3461290000
8 0.3643000000 -0.8983000000E-02
S 8
1 8236.000000 -0.1130000000E-03
2 1235.000000 -0.8780000000E-03
3 280.8000000 -0.4540000000E-02
4 79.27000000 -0.1813300000E-01
5 25.59000000 -0.5576000000E-01
6 8.997000000 -0.1268950000
7 3.319000000 -0.1703520000
8 0.3643000000 0.5986840000
S 1
1 0.9059000000 1.000000000
S 1
1 0.1285000000 1.000000000
P 3
1 18.71000000 0.1403100000E-01
2 4.133000000 0.8686600000E-01
3 1.200000000 0.2902160000
P 1
1 0.3827000000 1.000000000
P 1
1 0.1209000000 1.000000000
D 1
1 1.097000000 1.000000000
D 1
1 0.3180000000 1.000000000
F 1
1 0.7610000000 1.000000000
S 1
1 0.440200000E-01 1.00000000
P 1
1 0.356900000E-01 1.00000000
D 1
1 0.100000000 1.00000000
F 1
1 0.268000000 1.00000000
$end
--- OPTIMIZED MCSCF MO-S --- GENERATED 22-AUG-2000
E(MCSCF)= -37.7282408589, 11 ITERS
$VEC1
1 1 9.75511467E-01 ...snipped...
$END
~~