1
(19 Sep 03)
***********************************
* *
* Section 4 - Further Information *
* *
***********************************
This section of the manual contains both references, and
hints on how to do things. The following is a list of
the topics covered:
o Computational References.
o Basis Set References, and descriptions.
o Spherical Harmonics.
o How to do RHF, ROHF, UHF, and GVB calculations.
General considerations Other open shell SCF cases
Direct SCF True GVB perfect pairing
Convergence accelerators The special case of TCSCF
High spin open shell SCF A caution about symmetry
o How to do MCSCF and CI calculations.
MCSCF implementation CSF CI
Orbital updates starting orbitals
CI coef optimization references
determinant CI
o Second order perturbation theory.
RHF and UHF MP2 GVB based MP2
high spin ROHF MP2 MCSCF based MP2
o Coupled-Cluster theory.
Available computations restarts
resource requirements references
o Density Functional Theory.
DFTTYP keywords DFT with grids
grid-free DFT references
o Geometry Searches and Internal Coordinates.
Quasi-Newton searches Practical matters
The nuclear hessian Saddle points
Coordinate choices Mode following
The role of symmetry
o Intrinsic Reaction Coordinate (IRC) methods.
o Gradient Extremals.
1
o Continuum solvation methods: SCRF, PCM, and COSMO.
o Effective Fragment Potential method.
Terms in an EFP Current Limitations
Constructing EFP1 or EFP2 Practical Hints
QM/MM across covalent bonds References
o MOPAC calculations within GAMESS.
o Molecular Properties, and conversion factors.
o Localization tips.
o Transition moments and spin-orbit coupling.
For people who are newcomers to computational chemistry,
it may be helpful to study an introductory book. The
volume by Frank Jensen is an outstanding survey of methods,
basis sets, properties, and other topics.
"Ab Initio Molecular Orbital Theory"
W.J.Hehre, L.Radom, J.A.Pople, P.v.R.Schleyer
Wiley and Sons, New York, 1986
"Modern Quantum Chemistry" (now a Dover paperback)
A.Szabo, N.S.Ostlund McGraw-Hill, 1989
"Quantum Chemistry, 5th Edition"
I.N.Levine Prentice Hall, 1999
"Introduction to Computational Chemistry"
F.Jensen Wiley and Sons, Chichester, 1999
"Introduction to Quantum Mechanics in Chemistry"
M.A.Ratner, G.C.Schatz Prentice Hall, 2000
1
Computational References
------------- ----------
GAMESS -
M.W.Schmidt, K.K.Baldridge, J.A.Boatz, S.T.Elbert,
M.S.Gordon, J.J.Jensen, S.Koseki, N.Matsunaga,
K.A.Nguyen, S.Su, T.L.Windus, M.Dupuis, J.A.Montgomery
J.Comput.Chem. 14, 1347-1363 (1993)
HONDO -
These papers describes many of the algorithms in detail,
and much of these applies also to GAMESS:
"The General Atomic and Molecular Electronic Structure
System: HONDO 7.0" M.Dupuis, J.D.Watts, H.O.Villar,
G.J.B.Hurst Comput.Phys.Comm. 52, 415-425(1989)
"HONDO: A General Atomic and Molecular Electronic
Structure System" M.Dupuis, P.Mougenot, J.D.Watts,
G.J.B.Hurst, H.O.Villar in "MOTECC: Modern Techniques
in Computational Chemistry" E.Clementi, Ed.
ESCOM, Leiden, the Netherlands, 1989, pp 307-361.
"HONDO: A General Atomic and Molecular Electronic
Structure System" M.Dupuis, A.Farazdel, S.P.Karna,
S.A.Maluendes in "MOTECC: Modern Techniques in
Computational Chemistry" E.Clementi, Ed.
ESCOM, Leiden, the Netherlands, 1990, pp 277-342.
M.Dupuis, S.Chin, A.Marquez in "Relativistic and Electron
Correlation Effects in Molecules", G.Malli, Ed. Plenum
Press, NY 1994, pp 315-338.
sp integrals and gradient integrals -
J.A.Pople, W.J.Hehre J.Comput.Phys. 27, 161-168(1978)
H.B.Schlegel, J.Chem.Phys. 77, 3676-3681(1982)
spdfg integrals -
"Numerical Integration Using Rys Polynomials"
H.F.King and M.Dupuis J.Comput.Phys. 21,144(1976)
"Evaluation of Molecular Integrals over Gaussian
Basis Functions"
M.Dupuis,J.Rys,H.F.King J.Chem.Phys. 65,111-116(1976)
"Molecular Symmetry and Closed Shell HF Calculations"
M.Dupuis and H.F.King Int.J.Quantum Chem. 11,613(1977)
"Computation of Electron Repulsion Integrals using
the Rys Quadrature Method"
J.Rys,M.Dupuis,H.F.King J.Comput.Chem. 4,154-157(1983)
spdfg gradient integrals -
"Molecular Symmetry. II. Gradient of Electronic Energy
with respect to Nuclear Coordinates"
M.Dupuis and H.F.King J.Chem.Phys. 68,3998(1978)
although the implementation is much newer than this paper.
1
spd hessian integrals -
"Molecular Symmetry. III. Second derivatives of Electronic
Energy with respect to Nuclear Coordinates"
T.Takada, M.Dupuis, H.F.King
J.Chem.Phys. 75, 332-336 (1981)
Q matrix, and integral transformation symmetry -
E.Hollauer, M.Dupuis J.Chem.Phys. 96, 5220 (1992)
spdfg effective core potentials (ECP) integrals -
C.F.Melius, W.A.Goddard Phys.Rev.A, 10,1528-1540(1974)
L.R.Kahn, P.Baybutt, D.G.Truhlar
J.Chem.Phys. 65, 3826-3853 (1976)
M.Krauss, W.J.Stevens Ann.Rev.Phys.Chem. 35, 357-385(1985)
J.Breidung, W.Thiel, A.Komornicki
Chem.Phys.Lett. 153, 76-81(1988)
B.M.Bode, M.S.Gordon J.Chem.Phys. 111, 8778-8784(1999)
See also the papers listed for SBKJC and HW basis sets.
Quantum fast multipole method -
E.O.Steinborn, K.Ruedenberg Adv.Quantum Chem. 7, 1-81(1973)
L.Greengard "The Rapid Evaluation of Potential Fields in
Particle Systems" (MIT, Cambridge, 1987)
C.H.Choi, J.Ivanic, M.S.Gordon, K.Ruedenberg
J.Chem.Phys. 111, 8825-8831(1999)
C.H.Choi, K.Ruedenberg, M.S.Gordon
J.Comput.Chem. 22, 1484-1501(2001)
RHF -
C.C.J.Roothaan Rev.Mod.Phys. 23, 69(1951)
UHF -
J.A.Pople, R.K.Nesbet J.Chem.Phys 22, 571 (1954)
high spin coupling ROHF -
R.McWeeny, G.Diercksen J.Chem.Phys. 49,4852-4856(1968)
M.F.Guest, V.R.Saunders, Mol.Phys. 28, 819-828(1974)
J.S.Binkley, J.A.Pople, P.A.Dobosh
Mol.Phys. 28, 1423-1429 (1974)
E.R.Davidson Chem.Phys.Lett. 21,565(1973)
K.Faegri, R.Manne Mol.Phys. 31,1037-1049(1976)
H.Hsu, E.R.Davidson, and R.M.Pitzer
J.Chem.Phys. 65,609(1976)
GVB and low spin coupling ROHF -
F.W.Bobrowicz and W.A.Goddard, in Modern Theoretical
Chemistry, Vol 3, H.F.Schaefer III, Ed., Chapter 4.
MCSCF - see reference list in the subsection below
1
determinant CI (full CI (ALDET) and general CI (GENCI)) -
J.Ivanic, K.Ruedenberg Theoret.Chem.Acc. 106, 339-351(2001)
determinant CI (occupationally restricted, ORMAS) -
J.Ivanic J.Chem.Phys. to appear November 2003
configuration state function CI (GUGA) -
B.Brooks and H.F.Schaefer J.Chem. Phys. 70,5092(1979)
B.Brooks, W.Laidig, P.Saxe, N.Handy, and H.F.Schaefer,
Physica Scripta 21, 312(1980).
CIS energy and gradient -
J.B.Foresman, M.Head-Gordon, J.A.Pople, M.J.Frisch
J.Phys.Chem. 96, 135-149(1992)
R.M.Shroll, W.D.Edwards
Int.J.Quantum Chem. 63, 1037-1049(1997)
closed and unrestricted open shell 2nd order Moller-Plesset -
J.A.Pople, J.S.Binkley, R.Seeger
Int. J. Quantum Chem. S10, 1-19(1976)
M.J.Frisch, M.Head-Gordon, J.A.Pople,
Chem.Phys.Lett. 166, 275-280(1990)
G.D.Fletcher, M.W.Schmidt, M.S.Gordon
Adv.Chem.Phys. 110, 267-294(1999)
C.M.Aikens, S.P.Webb, R.L.Bell, G.D.Fletcher, M.W.Schmidt,
M.S.Gordon Theoret.Chem.Acc., in press
spin restricted open shell MP2, so called ZAPT method -
T.J.Lee, D.Jayatilaka Chem.Phys.Lett. 201, 1-10(1993)
T.J.Lee, A.P.Rendell, K.G.Dyall, D.Jayatilaka
J.Chem.Phys. 100, 7400-7409(1994)
spin restricted open shell MP2, so called RMP method -
P.J.Knowles, J.S.Andrews, R.D.Amos, N.C.Handy, J.A.Pople
Chem.Phys.Lett. 186, 130-136 (1991)
W.J.Lauderdale, J.F.Stanton, J.Gauss, J.D.Watts,
R.J.Bartlett Chem.Phys.Lett. 187, 21-28(1991)
multiconfigurational quasidegenerate perturbation theory -
H.Nakano, J.Chem.Phys. 99, 7983-7992(1993)
Coupled-Cluster Program included in GAMESS -
P. Piecuch, S.A. Kucharski, K. Kowalski, and M. Musial,
Comput.Phys.Commun., 149, 71-96(2002)
Any publication describing the results of coupled-cluster
calculations obtained with GAMESS should give reference
this paper.
1
RHF/ROHF/TCSCF coupled perturbed Hartree Fock -
"Single Configuration SCF Second Derivatives on a Cray"
H.F.King, A.Komornicki in "Geometrical Derivatives of
Energy Surfaces and Molecular Properties" P.Jorgensen
J.Simons, Ed. D.Reidel, Dordrecht, 1986, pp 207-214.
Y.Osamura, Y.Yamaguchi, D.J.Fox, M.A.Vincent, H.F.Schaefer
J.Mol.Struct. 103, 183-186 (1983)
M.Duran, Y.Yamaguchi, H.F.Schaefer
J.Phys.Chem. 92, 3070-3075 (1988)
"A New Dimension to Quantum Chemistry" Y.Yamaguchi,
Y.Osamura, J.D.Goddard, H.F.Schaefer Oxford Press, NY 1994
harmonic vibrational analysis in Cartesian coordinates -
W.D.Gwinn J.Chem.Phys. 55,477-481(1971)
Normal coordinate decomposition analysis -
J.A.Boatz and M.S.Gordon,
J.Phys.Chem. 93, 1819-1826(1989).
Partial Hessian vibrational analysis -
H.Li, J.H.Jensen, Theoret.Chem.Acc. 107, 211-219(2002)
Raman intensity -
A.Komornicki, J.W.McIver J.Chem.Phys. 70, 2014-2016(1979)
G.B.Bacskay, S.Saebo, P.R.Taylor
Chem.Phys. 90, 215-224(1984)
anharmonic vibrational spectra (VSCF) -
G.M.Chaban, J.O.Jung, R.B.Gerber
J.Chem.Phys. 111, 1823-1829(1999)
N.Matsunaga, G.M.Chaban, R.B.Gerber
J.Chem.Phys. 117, 3541-3547(2002)
R.B.Gerber, J.O.Jung in "Computational Molecular
Spectroscopy" P.Jensen, P.R.Bunker, eds. Wiley and Sons,
Chichester, 2000, pp 365-390.
some applications of RUNTYP=VSCF:
G.M.Chaban, J.O.Jung, R.B.Gerber
J.Phys.Chem.A 104, 2772-2779(2000)
J.Lundell, G.M.Chaban, R.B.Gerber
Chem.Phys.Lett. 331, 308-316(2000)
G.M.Chaban, R.B.Gerber, K.C.Janda
J.Phys.Chem.A 105, 8323-8332(2001)
A.T.Kowal, Spectrochimica Acta A 58, 1055-1067(2002)
G.M.Chaban, S.S.Xantheas, R.B.Gerber
J.Phys.Chem.A 107, 4952-4956(2003)
1
Geometry optimization and saddle point location -
J.Baker J.Comput.Chem. 7, 385-395(1986).
T.Helgaker Chem.Phys.Lett. 182, 503-510(1991).
P.Culot, G.Dive, V.H.Nguyen, J.M.Ghuysen
Theoret.Chim.Acta 82, 189-205(1992).
Dynamic Reaction Coordinate (DRC) -
J.J.P.Stewart, L.P.Davis, L.W.Burggraf,
J.Comput.Chem. 8, 1117-1123 (1987)
S.A.Maluendes, M.Dupuis, J.Chem.Phys. 93, 5902-5911(1990)
T.Taketsugu, M.S.Gordon, J.Phys.Chem. 99, 8462-8471(1995)
T.Taketsugu, M.S.Gordon, J.Phys.Chem. 99, 14597-604(1995)
T.Taketsugu, M.S.Gordon, J.Chem.Phys. 103, 10042-9(1995)
M.S.Gordon, G.Chaban, T.Taketsugu
J.Phys.Chem. 100, 11512-11525(1996)
T.Takata, T.Taketsugu, K.Hirao, M.S.Gordon
J.Chem.Phys. 109, 4281-4289(1998)
Energy orbital localization -
C.Edmiston, K.Ruedenberg Rev.Mod.Phys. 35, 457-465(1963).
R.C.Raffenetti, K.Ruedenberg, C.L.Janssen, H.F.Schaefer,
Theoret.Chim.Acta 86, 149-165(1993)
Boys orbital localization -
S.F.Boys, "Quantum Science of Atoms, Molecules, and Solids"
P.O.Lowdin, Ed, Academic Press, NY, 1966, pp 253-262.
Population orbital localization -
J.Pipek, P.Z.Mezey J.Chem.Phys. 90, 4916(1989).
Mulliken Population Analysis -
R.S.Mulliken J.Chem.Phys. 23, 1833-1840, 1841-1846,
2338-2342, 2343-2346(1955)
so called "Lowdin Population Analysis" -
This should be described as "a Mulliken population analysis
(ref M1-M4 above) based on symmetrically orthogonalized
orbitals (ref L)", where reference L is
P.-O.Lowdin Adv.Chem.Phys. 5, 185-199(1970)
Bond orders and valences -
M.Giambiagi, M.Giambiagi, D.R.Grempel, C.D.Heymann
J.Chim.Phys. 72, 15-22(1975)
I.Mayer, Chem.Phys.Lett., 97,270-274(1983), 117,396(1985).
M.S.Giambiagi, M.Giambiagi, F.E.Jorge
Z.Naturforsch. 39a, 1259-73(1984)
I.Mayer, Theoret.Chim.Acta, 67,315-322(1985).
I.Mayer, Int.J.Quantum Chem., 29,73-84(1986).
I.Mayer, Int.J.Quantum Chem., 29,477-483(1986).
The same formula (apart from a factor of two) may also be
seen in equation 31 of the second of these papers (the bond
order formula in the 1st of these is not the same formula):
T.Okada, T.Fueno Bull.Chem.Soc.Japan 48, 2025-2032(1975)
T.Okada, T.Fueno Bull.Chem.Soc.Japan 49, 1524-1530(1976)
1
Direct SCF -
J.Almlof, K.Faegri, K.Korsell
J.Comput.Chem. 3, 385-399 (1982)
M.Haser, R.Ahlrichs
J.Comput.Chem. 10, 104-111 (1989)
DIIS (Direct Inversion in the Iterative Subspace) -
P.Pulay J.Comput.Chem. 3, 556-560(1982)
SOSCF -
G.Chaban, M.W.Schmidt, M.S.Gordon
Theor.Chem.Acc. 97, 88-95(1997)
T.H.Fischer, J.Almlof, J.Phys.Chem. 96,9768-74(1992)
Modified Virtual Orbitals (MVOs) -
C.W.Bauschlicher, Jr. J.Chem.Phys. 72,880-885(1980)
EVVRSP in memory diagonalization -
S.T.Elbert Theoret.Chim.Acta 71,169-186(1987)
Davidson eigenvector method -
E.R.Davidson J.Comput.Phys. 17,87(1975) and "Matrix
Eigenvector Methods" p. 95 in "Methods in Computational
Molecular Physics" ed. by G.H.F.Diercksen and S.Wilson
RESC -
T.Nakajima, K.Hirao Chem.Phys.Lett. 302, 383-391(1999)
T.Nakajima, T.Suzumura, K.Hirao
Chem.Phys.Lett. 304, 271(1999)
D.G.Fedorov, T.Nakajima, K.Hirao
Chem.Phys.Lett. 335, 183-187(2001)
NESC -
K.G.Dyall J.Comput.Chem. 23, 786-793(2002)
MOPAC 6 -
J.J.P.Stewart J.Computer-Aided Molecular Design
4, 1-105 (1990)
References for parameters for individual atoms may be
found on the printout from your runs.
MacMolPlt -
B.M.Bode, M.S.Gordon J.Mol.Graphics Mod. 16, 133-138(1998)
1
parallelization in GAMESS -
for SCF, the main GAMESS paper quoted above.
T.L.Windus, M.W.Schmidt, M.S.Gordon,
Chem.Phys.Lett., 216, 375-379(1993)
T.L.Windus, M.W.Schmidt, M.S.Gordon,
Theoret.Chim.Acta 89, 77-88 (1994)
T.L.Windus, M.W.Schmidt, M.S.Gordon, in "Parallel Computing
in Computational Chemistry", ACS Symposium Series 592,
Ed. by T.G.Mattson, ACS Washington, 1995, pp 16-28.
K.K.Baldridge, M.S.Gordon, J.H.Jensen, N.Matsunaga,
M.W.Schmidt, T.L.Windus, J.A.Boatz, T.R.Cundari
ibid, pp 29-46.
G.D.Fletcher, M.W.Schmidt, M.S.Gordon
Adv.Chem.Phys. 110, 267-294 (1999)
G.D.Fletcher, M.W.Schmidt, B.M.Bode, M.S.Gordon
Comput.Phys.Commun. 128, 190-200 (2000)
H.Umeda, S.Koseki, U.Nagashima, M.W.Schmidt
J.Comput.Chem. 22, 1243-1251 (2001)
1
Basis Set References
----- --- ----------
An excellent review of the relationship between
the atomic basis used, and the accuracy with which
various molecular properties will be computed is:
E.R.Davidson, D.Feller Chem.Rev. 86, 681-696(1986).
STO-NG H-Ne Ref. 1 and 2
Na-Ar, Ref. 2 and 3 **
K,Ca,Ga-Kr Ref. 4
Rb,Sr,In-Xe Ref. 5
Sc-Zn,Y-Cd Ref. 6
1) W.J.Hehre, R.F.Stewart, J.A.Pople
J.Chem.Phys. 51, 2657-2664(1969).
2) W.J.Hehre, R.Ditchfield, R.F.Stewart, J.A.Pople
J.Chem.Phys. 52, 2769-2773(1970).
3) M.S.Gordon, M.D.Bjorke, F.J.Marsh, M.S.Korth
J.Am.Chem.Soc. 100, 2670-2678(1978).
** the valence scale factors for Na-Cl are taken
from this paper, rather than the "official"
Pople values in Ref. 2.
4) W.J.Pietro, B.A.Levi, W.J.Hehre, R.F.Stewart,
Inorg.Chem. 19, 2225-2229(1980).
5) W.J.Pietro, E.S.Blurock, R.F.Hout,Jr., W.J.Hehre, D.J.
DeFrees, R.F.Stewart Inorg.Chem. 20, 3650-3654(1980).
6) W.J.Pietro, W.J.Hehre J.Comput.Chem. 4, 241-251(1983).
MINI/MIDI H-Xe Ref. 9
9) "Gaussian Basis Sets for Molecular Calculations"
S.Huzinaga, J.Andzelm, M.Klobukowski, E.Radzio-Andzelm,
Y.Sakai, H.Tatewaki Elsevier, Amsterdam, 1984.
The MINI bases are three gaussian expansions of each
atomic orbital. The exponents and contraction
coefficients are optimized for each element, and s and p
exponents are not constrained to be equal. As a result
these bases give much lower energies than does STO-3G.
The valence MINI orbitals of main group elements are
scaled by factors optimized by John Deisz at North Dakota
State University. Transition metal MINI bases are not
scaled. The MIDI bases are derived from the MINI sets by
floating the outermost primitive in each valence orbitals,
and renormalizing the remaining 2 gaussians. MIDI bases
are not scaled by GAMESS. The transition metal bases are
taken from the lowest SCF terms in the s**1,d**n
configurations.
1
3-21G H-Ne Ref. 10 (also 6-21G)
Na-Ar Ref. 11 (also 6-21G)
K,Ca,Ga-Kr,Rb,Sr,In-Xe Ref. 12
Sc-Zn Ref. 13
Y-Cd Ref. 14
10) J.S.Binkley, J.A.Pople, W.J.Hehre
J.Am.Chem.Soc. 102, 939-947(1980).
11) M.S.Gordon, J.S.Binkley, J.A.Pople, W.J.Pietro,
W.J.Hehre J.Am.Chem.Soc. 104, 2797-2803(1982).
12) K.D.Dobbs, W.J.Hehre J.Comput.Chem. 7,359-378(1986)
13) K.D.Dobbs, W.J.Hehre J.Comput.Chem. 8,861-879(1987)
14) K.D.Dobbs, W.J.Hehre J.Comput.Chem. 8,880-893(1987)
N-31G references for 4-31G 5-31G 6-31G
H 15 15 15
He 23 23 23
Li 19,24 19
Be 20,24 20
B 17 19
C-F 15 16 16
Ne 23 23
Na-Al 22
Si 21 **
P-Cl 18 22
Ar 22
K-Zn 25
15) R.Ditchfield, W.J.Hehre, J.A.Pople
J.Chem.Phys. 54, 724-728(1971).
16) W.J.Hehre, R.Ditchfield, J.A.Pople
J.Chem.Phys. 56, 2257-2261(1972).
17) W.J.Hehre, J.A.Pople J.Chem.Phys. 56, 4233-4234(1972).
18) W.J.Hehre, W.A.Lathan J.Chem.Phys. 56,5255-5257(1972).
19) J.D.Dill, J.A.Pople J.Chem.Phys. 62, 2921-2923(1975).
20) J.S.Binkley, J.A.Pople J.Chem.Phys. 66, 879-880(1977).
21) M.S.Gordon Chem.Phys.Lett. 76, 163-168(1980)
** - Note that the built in 6-31G basis for Si is
not that given by Pople in reference 22.
The Gordon basis gives a better wavefunction,
for a ROHF calculation in full atomic (Kh)
symmetry,
6-31G Energy virial
Gordon -288.828573 1.999978
Pople -288.828405 2.000280
See the input examples for how to run in Kh.
22) M.M.Francl, W.J.Pietro, W.J.Hehre, J.S.Binkley,
M.S.Gordon, D.J.DeFrees, J.A.Pople
J.Chem.Phys. 77, 3654-3665(1982).
23) Unpublished, copied out of GAUSSIAN82.
24) For Li and Be, 4-31G is actually a 5-21G expansion.
25) V.A.Rassolov, J.A.Pople, M.A.Ratner, T.L.Windus
J.Chem.Phys. 109, 1223-1229(1998)
1
Extended basis sets
--> 6-311G
28) R.Krishnan, J.S.Binkley, R.Seeger, J.A.Pople
J.Chem.Phys. 72, 650-654(1980).
--> valence double zeta "DZV" sets:
"DH" basis - DZV for H, Li-Ne, Al-Ar
30) T.H.Dunning, Jr., P.J.Hay Chapter 1 in "Methods of
Electronic Structure Theory", H.F.Shaefer III, Ed.
Plenum Press, N.Y. 1977, pp 1-27.
Note that GAMESS uses inner/outer scale factors of
1.2 and 1.15 for DH's hydrogen (since at least 1983).
To get Thom's usual basis, scaled 1.2 throughout:
HYDROGEN 1.0 x, y, z
DH 0 1.2 1.2
DZV for K,Ca
31) J.-P.Blaudeau, M.P.McGrath, L.A.Curtiss, L.Radom
J.Chem.Phys. 107, 5016-5021(1997)
"BC" basis - DZV for Ga-Kr
32) R.C.Binning, Jr., L.A.Curtiss
J.Comput.Chem. 11, 1206-1216(1990)
--> valence triple zeta "TZV" sets:
TZV for H,Li-Ne
40) T.H. Dunning, J.Chem.Phys. 55 (1971) 716-723.
TZV for Na-Ar - also known as the "MC" basis
41) A.D.McLean, G.S.Chandler
J.Chem.Phys. 72,5639-5648(1980).
TZV for K,Ca
42) A.J.H. Wachters, J.Chem.Phys. 52 (1970) 1033-1036.
(see Table VI, Contraction 3).
TZV for Sc-Zn (taken from HONDO 7)
This is Wachters' (14s9p5d) basis (ref 42) contracted
to (10s8p3d) with the following modifications
1. the most diffuse s removed;
2. additional s spanning 3s-4s region;
3. two additional p functions to describe the 4p;
4. (6d) contracted to (411) from ref 43,
except for Zn where Wachter's (5d)/[41]
and Hay's diffuse d are used.
43) A.K. Rappe, T.A. Smedley, and W.A. Goddard III,
J.Phys.Chem. 85 (1981) 2607-2611
1
Valence only basis sets (for use with corresponding ECPs)
SBKJC -31G splits, bigger for trans. metals (available Li-Rn)
50) W.J.Stevens, H.Basch, M.Krauss
J.Chem.Phys. 81, 6026-6033 (1984)
51) W.J.Stevens, M.Krauss, H.Basch, P.G.Jasien
Can.J.Chem. 70, 612-630 (1992)
52) T.R.Cundari, W.J.Stevens
J.Chem.Phys. 98, 5555-5565(1993)
HW -21 splits (sp exponents not shared)
transition metals (not built in at present, although
they will work if you type them in).
53) P.J.Hay, W.R.Wadt J.Chem.Phys. 82, 270-283 (1985)
main group (available Na-Xe)
54) W.R.Wadt, P.J.Hay J.Chem.Phys. 82, 284-298 (1985)
see also
55) P.J.Hay, W.R.Wadt J.Chem.Phys. 82, 299-310 (1985)
Polarization exponents
STO-NG*
60) J.B.Collins, P. von R. Schleyer, J.S.Binkley,
J.A.Pople J.Chem.Phys. 64, 5142-5151(1976).
3-21G*. See also reference 12.
61) W.J.Pietro, M.M.Francl, W.J.Hehre, D.J.DeFrees, J.A.
Pople, J.S.Binkley J.Am.Chem.Soc. 104,5039-5048(1982)
6-31G* and 6-31G**. See also reference 22 above.
62) P.C.Hariharan, J.A.Pople
Theoret.Chim.Acta 28, 213-222(1973)
multiple polarization, and f functions
63) M.J.Frisch, J.A.Pople, J.S.Binkley J.Chem.Phys.
80, 3265-3269 (1984).
1
STO-NG* means d orbitals are used on third row atoms only.
The original paper (ref 60) suggested z=0.09 for
Na and Mg, and z=0.39 for Al-Cl.
At NDSU we prefer to use the same exponents as
in 3-21G* and 6-31G*, so we know we're looking
at changes in the sp basis, not the d exponent.
3-21G* means d orbitals on main group elements in the
third and higher periods. Not defined for the
transition metals, where there are p's already
in the basis. Except for alkalis and alkali
earths, the 4th and 5th row zetas are from
Huzinaga, et al (ref 9). The exponents are
normally the same as for 6-31G*.
6-31G* means d orbitals on second and third row atoms.
We use Mark Gordon's z=0.395 for silicon, as well
as his fully optimized sp basis (ref 21).
This is often written 6-31G(d) today.
For the first row transition metals, the *
means an f function is added.
6-31G** means the same as 6-31G*, except that p functions
are added on hydrogens.
This is often written 6-31G(d,p) today.
6-311G** means p orbitals on H, and d orbitals elsewhere.
The exponents were derived from correlated atomic
states, and so are considerably tighter than the
polarizing functions used in 6-31G**, etc.
This is often written 6-311G(d,p) today.
The definitions for 6-31G* for C-F are disturbing in
that they treat these atoms the same. Dunning and Hay
(ref 30) have recommended a better set of exponents for
second row atoms and a slightly different value for H.
2p, 3p, 2d, 3p polarization sets are usually thought
of as arising from applying splitting factors to the
1p and 1d values. For example, SPLIT2=2.0, 0.5 means
to double and halve the single value. The default
values for SPLIT2 and SPLIT3 are taken from reference
63, and were derived with correlation in mind. The
SPLIT2 values often produce a higher (!) HF energy than
the singly polarized run, because the exponents are
split too widely. SPLIT2=0.4,1.4 will always lower the
SCF energy (the values are the unpublished personal
preference of MWS), and for SPLIT3 we might suggest
3.0,1.0,1/3.
With all this as background, we are ready to present
the table of polarization exponents built into GAMESS.
1
Built in polarization exponents, chosen by POLAR=
in the $BASIS group. The values are for d functions
unless otherwise indicated.
Please note that the names associated with each
column are only generally descriptive. For example, the
column marked "Pople" contains a value for Si with which
John Pople would not agree, and the Ga-Kr values in this
column are actually from the Huzinaga "green book". The
exponents for K-Kr under "Dunning" are from Curtiss, et
al., not Thom Dunning. And so on.
POPLE POPN311 DUNNING HUZINAGA HONDO7
------ ------- ------- -------- ------
H 1.1(p) 0.75(p) 1.0(p) 1.0(p) 1.0(p)
He 1.1(p) 0.75(p) 1.0(p) 1.0(p) 1.0(p)
Li 0.2 0.200 0.076(p)
Be 0.4 0.255 0.164(p) 0.32
B 0.6 0.401 0.70 0.388 0.50
C 0.8 0.626 0.75 0.600 0.72
N 0.8 0.913 0.80 0.864 0.98
O 0.8 1.292 0.85 1.154 1.28
F 0.8 1.750 0.90 1.496 1.62
Ne 0.8 2.304 1.00 1.888 2.00
Na 0.175 0.061(p) 0.157
Mg 0.175 0.101(p) 0.234
Al 0.325 0.198 0.311
Si 0.395 0.262 0.388
P 0.55 0.340 0.465
S 0.65 0.421 0.542
Cl 0.75 0.514 0.619
Ar 0.85 0.617 0.696
K 0.2 0.260 0.039(p)
Ca 0.2 0.229 0.059(p)
Sc-Zn 0.8(f) N/A N/A N/A N/A
Ga 0.207 0.141
Ge 0.246 0.202
As 0.293 0.273
Se 0.338 0.315
Br 0.389 0.338
Kr 0.443 0.318
Rb 0.11 0.034(p)
Sr 0.11 0.048(p)
A blank means the value equals the "Pople" column.
Common d polarization for all sets ("green book"):
In Sn Sb Te I Xe
0.160 0.183 0.211 0.237 0.266 0.297
Tl Pb Bi Po At Rn
0.146 0.164 0.185 0.204 0.225 0.247
1
f polarization functions, from reference 63:
Li Be B C N O F Ne
0.15 0.26 0.50 0.80 1.00 1.40 1.85 2.50
Na Mg Al Si P S Cl Ar
0.15 0.20 0.25 0.32 0.45 0.55 0.70 --
Anion diffuse functions
3-21+G, 3-21++G, etc.
70) T.Clark, J.Chandrasekhar, G.W.Spitznagel, P. von R.
Schleyer J.Comput.Chem. 4, 294-301(1983)
71) G.W.Spitznagel, Diplomarbeit, Erlangen, 1982.
Anions usually require diffuse basis functions to
properly represent their spatial diffuseness. The use of
diffuse sp shells on atoms in the second and third rows is
denoted by a + sign, also adding diffuse s functions on
hydrogen is symbolized by ++. These designations can be
applied to any of the Pople bases, e.g. 3-21+G, 3-21+G*,
6-31++G**. The following exponents are for L shells,
except for H. For H-F, they are taken from ref 70. For
Na-Cl, they are taken directly from reference 71. These
values may be found in footnote 13 of reference 63.
For Ga-Br, In-I, and Tl-At these were optimized for the
atomic ground state anion, using ROHF with a flexible ECP
basis set, by Ted Packwood at NDSU.
H
0.0360
Li Be B C N O F
0.0074 0.0207 0.0315 0.0438 0.0639 0.0845 0.1076
Na Mg Al Si P S Cl
0.0076 0.0146 0.0318 0.0331 0.0348 0.0405 0.0483
Ga Ge As Se Br
0.0205 0.0222 0.0287 0.0318 0.0376
In Sn Sb Te I
0.0223 0.0231 0.0259 0.0306 0.0368
Tl Pb Bi Po At
0.0170 0.0171 0.0215 0.0230 0.0294
Additional information about diffuse functions and also
Rydberg type exponents can be found in reference 30.
1
The following atomic energies are from UHF
calculations (RHF on 1-S states), with p orbitals not
symmetry equivalenced, and using the default molecular
scale factors. They should be useful in picking a basis
of the desired energy accuracy, and estimating the correct
molecular total energies.
Atom state STO-2G STO-3G 3-21G 6-31G
H 2-S -.454397 -.466582 -.496199 -.498233
He 1-S -2.702157 -2.807784 -2.835680 -2.855160
Li 2-S -7.070809 -7.315526 -7.381513 -7.431236
Be 1-S -13.890237 -14.351880 -14.486820 -14.566764
B 2-P -23.395284 -24.148989 -24.389762 -24.519492
C 3-P -36.060274 -37.198393 -37.481070 -37.677837
N 4-S -53.093007 -53.719010 -54.105390 -54.385008
O 3-P -71.572305 -73.804150 -74.393657 -74.780310
F 2-P -95.015084 -97.986505 -98.845009 -99.360860
Ne 1-S -122.360485 -126.132546 -127.803825 -128.473877
Na 2-S -155.170019 -159.797148 -160.854065 -161.841425
Mg 1-S -191.507082 -197.185978 -198.468103 -199.595219
Al 2-P -233.199965 -239.026471 -240.551046 -241.854186
Si 3-P -277.506857 -285.563052 -287.344431 -288.828598
P 4-S -327.564244 -336.944863 -339.000079 -340.689008
S 3-P -382.375012 -393.178951 -395.551336 -397.471414
Cl 2-P -442.206260 -454.546015 -457.276552 -459.442939
Ar 1-S -507.249273 -521.222881 -524.342962 -526.772151
SCF *
Atom state DH 6-311G MC limit
H 2-S -.498189 -.499810 -- -0.5
He 1-S -- -2.859895 -- -2.861680
Li 2-S -7.431736 -7.432026 -- -7.432727
Be 1-S -14.570907 -14.571874 -- -14.573023
B 2-P -24.526601 -24.527020 -- -24.529061
C 3-P -37.685571 -37.686024 -- -37.688619
N 4-S -54.397260 -54.397980 -- -54.400935
O 3-P -74.802707 -74.802496 -- -74.809400
F 2-P -99.395013 -99.394158 -- -99.409353
Ne 1-S -128.522354 -128.522553 -- -128.547104
Na 2-S -- -- -161.845587 -161.858917
Mg 1-S -- -- -199.606558 -199.614636
Al 2-P -241.855079 -- -241.870014 -241.876699
Si 3-P -288.829617 -- -288.847782 -288.854380
P 4-S -340.689043 -- -340.711346 -340.718798
S 3-P -397.468667 -- -397.498023 -397.504910
Cl 2-P -459.435938 -- -459.473412 -459.482088
Ar 1-S -- -- -526.806626 -526.817528
* M.W.Schmidt and K.Ruedenberg, J.Chem.Phys. 71,
3951-3962(1979). These are ROHF energies in Kh symmetry.
1
Spherical Harmonics
--------- ---------
The implementation of ISPHER in $CONTRL does not rely
on using a spherical harmonic basis set, in fact the atomic
basis remains the Cartesian Gaussians. Instead, certain
MOs formed from particular combinations of the Cartesian
Gaussians (for example, xx+yy+zz) are deleted from the MO
space. Thus a run with ISPHER=1 will have fewer MOs than
AOs. Since neither the occupied nor virtual MOs contain
any admixture of xx+yy+zz, the resulting energy and wave-
function is exactly equivalent to the use of a spherical
harmonic basis.
The log file output will contain expansions of each MO
in terms of 6 d's, 10 f's, and 15 g's, and the $VEC also
contains the same expansion over Cartesian Gaussians. Both
the matrix in your log file and in $VEC will contain fewer
MOs than AOs, the exact number of MOs used is printed in
the initial guess section of the log file. It should be
possible to read such $VEC groups into runs with different
settings of ISPHER, should you choose to do so.
The advantage of this approach is that intelligence in
the generation of symmetry orbitals combined with the
capability to drop linearly dependent MO combinations means
that the details of ISPHER are located only in the orbital
optimization code, where the variational spaces are simply
reduced in size to eliminate the undesired contaminant
functions. This means that none of the integral routines
need be modified, as the atomic basis remains the Cartesian
Gaussians. The disadvantage is that AO integral files run
over the Cartesian Gaussians, and thus are not reduced in
size. Of course transformed MO integrals and various
computations in correlated calculations are reduced in
size, since the number of MOs may be greatly reduced.
Computationally, the advantages of ISPHER=1 are not
limited to the reduced CPU time associated with fewer total
MOs. Questions about d orbital participation as measured
by Mulliken populations are cleanly addressed when the d's
usage in the MOs does not contain any contamination from
the s shape xx+yy+zz. Less obviously, the use of spherical
harmonics frequently greatly reduces problems with linear
dependency, that exhibit as poor SCF convergence.
1
How to do RHF, ROHF, UHF, and GVB calculations
--- -- -- --- ---- --- --- --- ------------
* * * General considerations * * *
These four SCF wavefunctions are all based on Fock
operator techniques, even though some GVB runs use more
than one determinant. Thus all of these have an intrinsic
N**4 time dependence, because they are all driven by
integrals in the AO basis. This similarity makes it
convenient to discuss them all together. In this section
we will use the term HF to refer generically to any of
these four wavefunctions, including the multi-determinate
GVB-PP functions. $SCF is the main input group for all
these HF wavefunctions.
As will be discussed below, in GAMESS the term ROHF
refers to high spin open shell SCF only, but other open
shell coupling cases are possible using the GVB code.
Analytic gradients are implemented for every possible
HF type calculation possible in GAMESS, and therefore
numerical hessians are available for each.
Analytic hessian calculation is implemented for RHF,
ROHF, and any GVB case with NPAIR=0 or NPAIR=1. Analytic
hessians are more accurate, and much more quickly computed
than numerical hessians, but require additional disk
storage to perform an integral transformation, and also
more physical memory.
The second order Moller-Plesset energy correction
(MP2) is implemented for RHF, UHF, ROHF, and MCSCF wave
functions. Analytic gradients may be obtained for MP2
with RHF or UHF reference wavefunctions, and MP2 level
properties are therefore available only for these two,
see MP2PRP in $MP2. All other cases give properties for
the SCF function.
Direct SCF is implemented for every possible HF type
calculation. The direct SCF method may not be used with
DEM convergence. Direct SCF may be used during energy,
gradient, numerical or analytic hessian, CI or MP2 energy
correction, or localized orbitals computations.
1
* * * direct SCF * * *
Normally, HF calculations proceed by evaluating a
large number of two electron repulsion integrals, and
storing these on a disk. This integral file is read back
in during each HF iteration to form the appropriate Fock
operators. In a direct HF, the integrals are not stored
on disk, but are instead reevaluated during each HF
iteration. Since the direct approach *always* requires
more CPU time, the default for DIRSCF in $SCF is .FALSE.
Even though direct SCF is slower, there are at least
three reasons why you may want to consider using it. The
first is that it may not be possible to store all of the
integrals on the disk drives attached to your computer.
Secondly, the index label packing scheme used by GAMESS
restricts the basis set size to no more than 361 if the
integrals are stored on disk, whereas for direct HF you
can (in principle) use up to 2047 basis functions.
Finally, what you are really interested in is reducing the
wall clock time to obtain your answer, not the CPU time.
Workstations have modest hardware (and sometimes software)
I/O capabilities. Other environments such as an IBM
mainframe shared by many users may also have very poor
CPU/wall clock performance for I/O bound jobs such as
conventional HF.
You can estimate the disk storage requirements for
conventional HF using a P or PK file by the following
formulae:
nint = 1/sigma * 1/8 * N**4
Mbytes = nint * x / 1024**2
Here N is the total number of basis functions in your
run, which you can learn from an EXETYP=CHECK run. The
1/8 accounts for permutational symmetry within the
integrals. Sigma accounts for the point group symmetry,
and is difficult to estimate accurately. Sigma cannot be
smaller than 1, in no symmetry (C1) calculations. For
benzene, sigma would be almost six, since you generate 6
C's and 6 H's by entering only 1 of each in $DATA. For
water sigma is not much larger than one, since most of the
basis set is on the unique oxygen, and the C2v symmetry
applies only to the H atoms. The factor x is 12 bytes per
integral for basis sets smaller than 255, and 16 otherwise.
Finally, since integrals that are very close to zero need
not be stored on disk, the actual power dependence is not
as bad as N**4, and in fact in the limit of very large
molecules can be as low as N**2. Thus plugging in sigma=1
should give you an upper bound to the actual disk space
needed. If the estimate exceeds your available disk
storage, your only recourse is direct HF.
1
What are the economics of direct HF? Naively, if we
assume the run takes 10 iterations to converge, we must
spend 10 times more CPU time doing the integrals on each
iteration. However, we do not have to waste any CPU time
reading blocks of integrals from disk, or in unpacking
their indices. We also do not have to waste any wall
clock time waiting for a relatively slow mechanical device
such as a disk to give us our data.
There are some less obvious savings too, as first
noted by Almlof. First, since the density matrix is known
while we are computing integrals, we can use the Schwarz
inequality to avoid doing some of the integrals. In a
conventional SCF this inequality is used to avoid doing
small integrals. In a direct SCF it can be used to avoid
doing integrals whose contribution to the Fock matrix is
small (density times integral=small). Secondly, we can
form the Fock matrix by calculating only its change since
the previous iteration. The contributions to the change
in the Fock matrix are equal to the change in the density
times the integrals. Since the change in the density goes
to zero as the run converges, we can use the Schwarz
screening to avoid more and more integrals as the
calculation progresses. The input option FDIFF in $SCF
selects formation of the Fock operator by computing only
its change from iteration to iteration. The FDIFF option
is not implemented for GVB since there are too many density
matrices from the previous iteration to store, but is the
default for direct RHF, ROHF, and UHF.
So, in our hypothetical 10 iteration case, we do not
spend as much as 10 times more time in integral
evaluation. Additionally, the run as a whole will not
slow down by whatever factor the integral time is
increased. A direct run spends no additional time summing
integrals into the Fock operators, and no additional time
in the Fock diagonalizations. So, generally speaking, a
RHF run with 10-15 iterations will slow down by a factor
of 2-4 times when run in direct mode. The energy gradient
time is unchanged by direct HF, and this is a large time
compared to HF energy, so geometry optimizations will be
slowed down even less. This is really the converse of
Amdahl's law: if you slow down only one portion of a
program by a large amount, the entire program slows down
by a much smaller factor.
To make this concrete, here are some times for GAMESS
for a job which is a RHF energy for a SbC4O2NH4. These
timings were obtained an extremely long time ago, on a
DECstation 3100 under Ultrix 3.1, which was running only
1
these tests, so that the wall clock times are meaningful.
This system is typical of Unix workstations in that it
uses SCSI disks, and the operating system is not terribly
good at disk I/O. By default GAMESS stores the integrals
on disk in the form of a P supermatrix, because this will
save time later in the SCF cycles. By choosing NOPK=1 in
$INTGRL, an ordinary integral file can be used, which
typically contains many fewer integrals, but takes more
CPU time in the SCF. Because the DECstation is not
terribly good at I/O, the wall clock time for the ordinary
integral file is actually less than when the supermatrix
is used, even though the CPU time is longer. The run
takes 13 iterations to converge, the times are in seconds.
P supermatrix ordinary file
# nonzero integrals 8,244,129 6,125,653
# blocks skipped 55,841 55,841
CPU time (ints) 709 636
CPU time (SCF) 1289 1472
CPU time (job total) 2123 2233
wall time (job total) 3468 3200
When the same calculation is run in direct mode
(integrals are processed like an ordinary integral disk
file when running direct),
iteration 1: FDIFF=.TRUE. FDIFF=.FALSE.
# nonzero integrals 6,117,416 6,117,416
# blocks skipped 60,208 60,208
iteration 13:
# nonzero integrals 3,709,733 6,122,912
# blocks skipped 105,278 59,415
CPU time (job total) 6719 7851
wall time (job total) 6764 7886
Note that elimination of the disk I/O dramatically
increases the CPU/wall efficiency. Here's the bottom line
on direct HF:
best direct CPU / best disk CPU = 6719/2123 = 3.2
best direct wall/ best disk wall= 6764/3200 = 2.1
Direct SCF is slower than conventional SCF, but not
outrageously so! From the data in the tables, we can see
that the best direct method spends about 6719-1472 = 5247
seconds doing integrals. This is an increase of about
5247/636 = 8.2 in the time spent doing integrals, in a run
which does 13 iterations (13 times evaluating integrals).
8.2 is less than 13 because the run avoids all CPU charges
related to I/O, and makes efficient use of the Schwarz
inequality to avoid doing many of the integrals in its
final iterations.
1
* * * convergence accelerators * * *
Generally speaking, the simpler the HF function, the
better its convergence. In our experience, the majority
of RHF and ROHF runs converge readily from GUESS=HUCKEL.
UHF often takes considerably more iterations than either
of these, due to the common occurence of heavy spin
contamination. GVB runs typically require GUESS=MOREAD,
although the Huckel guess usually works for NPAIR=0.
GVB cases with NPAIR greater than one are particularly
difficult.
Unfortunately, not all HF runs converge readily. The
best way to improve your convergence is to provide better
starting orbitals! In many cases, this means to MOREAD
orbitals from some simpler HF case. For example, if you
want to do a doublet ROHF, and the HUCKEL guess does not
seem to converge, do this: Do an RHF on the +1 cation.
RHF is typically more stable than ROHF, UHF, or GVB, and
cations are usually readily convergent. Then MOREAD the
cation's orbitals into the neutral calculation which you
wanted to do at first.
GUESS=HUCKEL does not always guess the correct
electronic configuration. It may be useful to use PRTMO
in $GUESS during a CHECK run to examine the starting
orbitals, and then reorder them with NORDER if that seems
appropriate.
Of course, by default GAMESS uses the convergence
procedures which are usually most effective. Still, there
are cases which are difficult, so the $SCF group permits
you to select several alternative methods for improving
convergence. Briefly, these are
EXTRAP. This extrapolates the three previous Fock
matrices, in an attempt to jump ahead a bit faster. This
is the most powerful of the old-fashioned accelerators,
and normally should be used at the beginning of any SCF
run. When an extrapolation occurs, the counter at the
left of the SCF printout is set to zero.
DAMP. This damps the oscillations between several
successive Fock matrices. It may help when the energy is
seen to oscillate wildly. Thinking about which orbitals
should be occupied initially may be an even better way to
avoid oscillatory behaviour.
1
SHIFT. This shifts the diagonal elements of the virtual
part of the Fock matrix up, in an attempt to uncouple the
unoccupied orbitals from the occupied ones. At convergence,
this has no effect on the orbitals, just their orbital
energies, but will produce different (and hopefully better)
orbitals during the iterations.
RSTRCT. This limits mixing of the occupied orbitals
with the empty ones, especially the flipping of the HOMO
and LUMO to produce undesired electronic configurations or
states. This should be used with caution, as it makes it
very easy to converge on incorrect electronic configurations,
especially if DIIS is also used. If you use this, be sure
to check your final orbital energies to see if they are
sensible. A lower energy for an unoccupied orbital than
for one of the occupied ones is a sure sign of problems.
DIIS. Direct Inversion in the Iterative Subspace is
a modern method, due to Pulay, using stored error and Fock
matrices from a large number of previous iterations to
interpolate an improved Fock matrix. This method was
developed to improve the convergence at the final stages
of the SCF process, but turns out to be quite powerful at
forcing convergence in the initial stages of SCF as well.
By giving ETHRSH as 10.0 in $SCF, you can practically
guarantee that DIIS will be in effect from the first
iteration. The default is set up to do a few iterations
with conventional methods (extrapolation) before engaging
DIIS. This is because DIIS can sometimes converge to
solutions of the SCF equations that do not have the lowest
possible energy. For example, the 3-A-2 small angle state
of SiLi2 (see M.S.Gordon and M.W.Schmidt, Chem.Phys.Lett.,
132, 294-8(1986)) will readily converge with DIIS to a
solution with a reasonable S**2, and an energy about 25
milliHartree above the correct answer. A SURE SIGN OF
TROUBLE WITH DIIS IS WHEN THE ENERGY RISES TO ITS FINAL
VALUE. However, if you obtain orbitals at one point on a
PES without DIIS, the subsequent use of DIIS with MOREAD
will probably not introduce any problems. Because DIIS
is quite powerful, EXTRAP, DAMP, and SHIFT are all turned
off once DIIS begins to work. DEM and RSTRCT will still
be in use, however.
SOSCF. Approximate second-order (quasi-Newton) SCF
orbital optimization. SOSCF will converge about as well as
DIIS at the initial geometry, and slightly better at
subsequent geometries. There's a bit less work solving the
SCF equations, too. The method kicks in after the orbital
gradient falls below SOGTOL. Some systems, particularly
transition metals with ECP basis sets, may have Huckel
orbitals for which the gradient is much larger than SOGTOL.
In this case it is probably better to use DIIS instead,
with a large ETHRSH, rather than increasing SOGTOL, since
you may well be outside the quadratic convergence region.
SOSCF does not exhibit true second order convergence since
it uses an approximation to the inverse hessian. SOSCF
1
will work for MOPAC runs, but is slower in this case.
SOSCF will work for UHF, but the convergence is slower than
DIIS. SOSCF will work for non-Abelian ROHF cases, but may
encounter problems if the open shell is degenerate.
DEM. Direct energy minimization should be your last
recourse. It explores the "line" between the current
orbitals and those generated by a conventional change in
the orbitals, looking for the minimum energy on that line.
DEM should always lower the energy on every iteration,
but is very time consuming, since each of the points
considered on the line search requires evaluation of a
Fock operator. DEM will be skipped once the density
change falls below DEMCUT, as the other methods should
then be able to affect final convergence. While DEM is
working, RSTRCT is held to be true, regardless of the
input choice for RSTRCT. Because of this, it behooves
you to be sure that the initial guess is occupying the
desired orbitals. DEM is available only for RHF. The
implementation in GAMESS resembles that of R.Seeger and
J.A.Pople, J.Chem.Phys. 65, 265-271(1976). Simultaneous
use of DEM and DIIS resembles the ADEM-DIOS method of
H.Sellers, Chem.Phys.Lett. 180, 461-465(1991). DEM does
not work with direct SCF.
* * * High spin open shell SCF (ROHF) * * *
Open shell SCF calculations are performed in GAMESS by
both the ROHF code and the GVB code. Note that when the
GVB code is executed with no pairs, the run is NOT a true
GVB run, and should be referred to in publications and
discussion as a ROHF calculation.
The ROHF module in GAMESS can handle any number of
open shell electrons, provided these have a high spin
coupling. Some commonly occurring cases are:
one open shell, doublet:
$CONTRL SCFTYP=ROHF MULT=2 $END
two open shells, triplet:
$CONTRL SCFTYP=ROHF MULT=3 $END
m open shells, high spin:
$CONTRL SCFTYP=ROHF MULT=m+1 $END
1
John Montgomery (then at United Technologies) is
responsible for the current ROHF implementation in GAMESS.
The following discussion is due to him:
The Fock matrix in the MO basis has the form
closed open virtual
closed F2 | Fb | (Fa+Fb)/2
-----------------------------------
open Fb | F1 | Fa
-----------------------------------
virtual (Fa+Fb)/2 | Fa | F0
where Fa and Fb are the usual alpha and beta Fock
matrices any UHF program produces. The Fock operators
for the doubly, singly, and zero occupied blocks can be
written as
F2 = Acc*Fa + Bcc*Fb
F1 = Aoo*Fa + Boo*Fb
F0 = Avv*Fa + Bvv*Fb
Some choices found in the literature for these
canonicalization coefficients are
Acc Bcc Aoo Boo Avv Bvv
Guest and Saunders 1/2 1/2 1/2 1/2 1/2 1/2
Roothaan single matrix -1/2 3/2 1/2 1/2 3/2 -1/2
Davidson 1/2 1/2 1 0 1 0
Binkley, Pople, Dobosh 1/2 1/2 1 0 0 1
McWeeny and Diercksen 1/3 2/3 1/3 1/3 2/3 1/3
Faegri and Manne 1/2 1/2 1 0 1/2 1/2
The choice of the diagonal blocks is arbitrary, as
ROHF is converged when the off diagonal blocks go to zero.
The exact choice for these blocks can however have an
effect on the convergence rate. This choice also affects
the MO coefficients, and orbital energies, as the
different choices produce different canonical orbitals
within the three subspaces. All methods, however, will
give identical total wavefunctions, and hence identical
properties such as gradients and hessians.
The default coupling case in GAMESS is the Roothaan
single matrix set. Note that pre-1988 versions of GAMESS
produced "Davidson" orbitals. If you would like to fool
around with any of these other canonicalizations, the Acc,
Aoo, Avv and Bcc, Boo, Bvv parameters can be input as the
first three elements of ALPHA and BETA in $SCF.
1
* * * Other open shell SCF cases (GVB) * * *
Genuine GVB-PP runs will be discussed later in this
section. First, we will consider how to do open shell SCF
with the GVB part of the program.
It is possible to do other open shell cases with the
GVB code, which can handle the following cases:
one open shell, doublet:
$CONTRL SCFTYP=GVB MULT=2 $END
$SCF NCO=xx NSETO=1 NO(1)=1 $END
two open shells, triplet:
$CONTRL SCFTYP=GVB MULT=3 $END
$SCF NCO=xx NSETO=2 NO(1)=1,1 $END
two open shells, singlet:
$CONTRL SCFTYP=GVB MULT=1 $END
$SCF NCO=xx NSETO=2 NO(1)=1,1 $END
Note that the first two cases duplicate runs which the
ROHF module can do better. Note that all of these cases
are really ROHF, since the default for NPAIR in $SCF is 0.
Many open shell states with degenerate open shells
(for example, in diatomic molecules) can be treated as
well. There is a sample of this in the 'Input Examples'
section of this manual.
If you would like to do any cases other than those
shown above, you must derive the coupling coefficients
ALPHA and BETA, and input them with the occupancies F in
the $SCF group.
Mariusz Klobukowski of the University of Alberta has
shown how to obtain coupling coefficients for the GVB open
shell program for many such open shell states. These can
be derived from the values in Appendix A of the book "A
General SCF Theory" by Ramon Carbo and Josep M. Riera,
Springer-Verlag (1978). The basic rule is
(1) F(i) = 1/2 * omega(i)
(2) ALPHA(i) = alpha(i)
(3) BETA(i) = - beta(i),
where omega, alpha, and beta are the names used by Ramon
in his Tables.
The variable NSETO should give the number of open
shells, and NO should give the degeneracy of each open
shell. Thus the 5-S state of carbon would have NSETO=2,
and NO(1)=1,3.
1
Some specific examples, for the lowest term in each
of the atomic P**N configurations are
! p**1 2-P state
$CONTRL SCFTYP=GVB MULT=2 $END
$SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE.
F(1)= 1.0 0.16666666666667
ALPHA(1)= 2.0 0.33333333333333 0.00000000000000
BETA(1)= -1.0 -0.16666666666667 -0.00000000000000 $END
! p**2 3-P state
$CONTRL SCFTYP=GVB MULT=3 $END
$SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE.
F(1)= 1.0 0.333333333333333
ALPHA(1)= 2.0 0.66666666666667 0.16666666666667
BETA(1)= -1.0 -0.33333333333333 -0.16666666666667 $END
! p**3 4-S state
$CONTRL SCFTYP=ROHF MULT=4 $END
! p**4 3-P state
$CONTRL SCFTYP=GVB MULT=3 $END
$SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE.
F(1)= 1.0 0.66666666666667
ALPHA(1)= 2.0 1.33333333333333 0.83333333333333
BETA(1)= -1.0 -0.66666666666667 -0.50000000000000 $END
! p**5 2-P state
$CONTRL SCFTYP=GVB MULT=2 $END
$SCF NCO=xx NSETO=1 NO=3 COUPLE=.TRUE.
F(1)= 1.0 0.83333333333333
ALPHA(1)= 2.0 1.66666666666667 1.33333333333333
BETA(1)= -1.0 -0.83333333333333 -0.66666666666667 $END
Be sure to give all the digits, as these are part of
a double precision energy formula.
Coupling constants for d**N configurations are
! d**1 2-D state
$CONTRL SCFTYP=GVB MULT=2 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.1
ALPHA(1)= 2.0, 0.20, 0.00
BETA(1)=-1.0,-0.10, 0.00 $END
! d**2 average of 3-F and 3-P states
$CONTRL SCFTYP=GVB MULT=3 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.2
ALPHA(1)= 2.0, 0.40, 0.05
BETA(1)=-1.0,-0.20,-0.05 $END
1
! d**3 average of 4-F and 4-P states
$CONTRL SCFTYP=GVB MULT=4 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.3
ALPHA(1)= 2.0, 0.60, 0.15
BETA(1)=-1.0,-0.30,-0.15 $END
! d**4 5-D state
$CONTRL SCFTYP=GVB MULT=5 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.4
ALPHA(1)= 2.0, 0.80, 0.30
BETA(1)=-1.0,-0.40,-0.30 $END
! d**5 6-S state
$CONTRL SCFTYP=ROHF MULT=6 $END
! d**6 5-D state
$CONTRL SCFTYP=GVB MULT=5 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.6
ALPHA(1)= 2.0, 1.20, 0.70
BETA(1)=-1.0,-0.60,-0.50 $END
! d**7 average of 4-F and 4-P states
$CONTRL SCFTYP=GVB MULT=4 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.7
ALPHA(1)= 2.0, 1.40, 0.95
BETA(1)=-1.0,-0.70,-0.55 $END
! d**8 average of 3-F and 3-P states
$CONTRL SCFTYP=GVB MULT=3 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.8
ALPHA(1)= 2.0, 1.60, 1.25
beta(1)=-1.0,-0.80,-0.65 $end
! d**9 2-D state
$CONTRL SCFTYP=GVB MULT=2 $END
$SCF NCO=xx NSETO=1 NO=5 COUPLE=.TRUE. F(1)=1.0,0.9
ALPHA(1)= 2.0, 1.80, 1.60
BETA(1)=-1.0,-0.90,-0.80 $END
The source for these values is R.Poirier, R.Kari, and
I.G.Csizmadia's book "Handbook of Gaussian Basis Sets",
Elsevier, Amsterdam, 1985.
Note that GAMESS can do a proper calculation on the ground
terms for the d**2, d**3, d**7, and d**8 configurations
only by means of state averaged MCSCF. For d**8, use
$CONTRL SCFTYP=MCSCF MULT=3 $END
$DRT GROUP=C1 FORS=.TRUE. NMCC=xx NDOC=3 NALP=2 $END
$GUGDIA NSTATE=10 $END
$GUGDM2 WSTATE(1)=1,1,1,1,1,1,1,0,0,0 $END
Open shell cases such as s**1,d**n are probably most easily
tackled with the state-averaged MCSCF program.
1
* * * True GVB perfect pairing runs * * *
True GVB runs are obtained by choosing NPAIR nonzero.
If you wish to have some open shell electrons in addition
to the geminal pairs, you may add the pairs to the end of
any of the GVB coupling cases shown above. The GVB module
assumes that you have reordered your MOs into the order:
NCO double occupied orbitals, NSETO sets of open shell
orbitals, and NPAIR sets of geminals (with NORDER=1 in the
$GUESS group).
Each geminal consists of two orbitals and contains two
singlet coupled electrons (perfect pairing). The first MO
of a geminal is probably heavily occupied (such as a
bonding MO u), and the second is probably weakly occupied
(such as an antibonding, correlating orbital v). If you
have more than one pair, you must be careful that the
initial MOs are ordered u1, v1, u2, v2..., which is -NOT-
the same order that RHF starting orbitals will be found
in. Use NORDER=1 to get the correct order.
These pair wavefunctions are actually a limited form
of MCSCF. GVB runs are much faster than MCSCF runs,
because the natural orbital u,v form of the wavefunction
permits a Fock operator based optimization. However,
convergence of the GVB run is by no means assured. The
same care in selecting the correlating orbitals that you
would apply to an MCSCF run must also be used for GVB
runs. In particular, look at the orbital expansions when
choosing the starting orbitals, and check them again after
the run converges.
GVB runs will be carried out entirely in orthonormal
natural u,v form, with strong orthogonality enforced on
the geminals. Orthogonal orbitals will pervade your
thinking in both initial orbital selection, and the entire
orbital optimization phase (the CICOEF values give the
weights of the u,v orbitals in each geminal). However,
once the calculation is converged, the program will
generate and print the nonorthogonal, generalized valence
bond orbitals. These GVB orbitals are an entirely
equivalent way of presenting the wavefunction, but are
generated only after the fact.
Convergence of true GVB runs is by no means as certain
as convergence of RHF, UHF, ROHF, or GVB with NPAIR=0.
You can assist convergence by doing a preliminary RHF or
ROHF calculation, and use these orbitals for GUESS=MOREAD.
Few, if any, GVB runs with NPAIR non-zero will converge
without using GUESS=MOREAD. Generation of MVOs during the
prelimnary SCF can also be advantageous. In fact, all the
advice outlined for MCSCF computations below is germane,
for GVB-PP is a type of MCSCF computation.
1
The total number of electrons in the GVB wavefunction
is given by the following formula:
NE = 2*NCO + sum 2*F(i)*NO(i) + 2*NPAIR
i
The charge is obtained by subtracting the total number of
protons given in $DATA. The multiplicity is implicit in
the choice of alpha and beta constants. Note that ICHARG
and MULT must be given correctly in $CONTRL anyway, as the
number of electrons from this formula is double checked
against the ICHARG value.
* * * the special case of TCSCF * * *
The wavefunction with NSETO=0 and NPAIR=1 is called
GVB-PP(1) by Goddard, two configuration SCF (TCSCF) by
Schaefer or Davidson, and CAS-SCF with two electrons in
two orbitals by others. Note that this is just semantics,
as these are all identical. This is a very important
type of wavefunction, as TCSCF is the minimum acceptable
treatment for singlet biradicals. The TCSCF wavefunction
can be obtained with SCFTYP=MCSCF, but it is usually much
faster to use the Fock based SCFTYP=GVB. Because of its
importance, the TCSCF function (if desired, with possible
open shells) permits analytic hessian computation.
* * * A caution about symmetry * * *
Caution! Some exotic calculations with the GVB
program do not permit the use of symmetry. The symmetry
algorithm in GAMESS was "derived assuming that the
electronic charge density transforms according to the
completely symmetric representation of the point group",
Dupuis/King, JCP, 68, 3998(1978). This may not be true
for certain open shell cases, and in fact during GVB runs,
it may not be true for closed shell singlet cases!
First, consider the following correct input for the
singlet-delta state of NH:
$CONTRL SCFTYP=GVB NOSYM=1 $END
$SCF NCO=3 NSETO=2 NO(1)=1,1 $END
for the x**1y**1 state, or for the x**2-y**2 state,
$CONTRL SCFTYP=GVB NOSYM=1 $END
$SCF NCO=3 NPAIR=1 CICOEF(1)=0.707,-0.707 $END
Neither gives correct results, unless you enter NOSYM=1.
The electronic term symbol is degenerate, a good tip off
that symmetry cannot be used. However, some degenerate
states can still use symmetry, because they use coupling
constants averaged over all degenerate states within a
single term, as is done in EXAM15 and EXAM16. Here the
"state averaged SCF" leads to a charge density which is
symmetric, and these runs can exploit symmetry.
1
Secondly, since GVB runs exploit symmetry for each
of the "shells", or type of orbitals, some calculations on
totally symmetric states may not be able to use symmetry.
An example is CO or N2, using a three pair GVB to treat
the sigma and pi bonds. Individual configurations such
as (sigma)**2,(pi-x)**2,(pi-y*)**2 do not have symmetric
charge densities since neither the pi nor pi* level is
completely filled. Correct answers for the sigma-plus
ground states result only if you input NOSYM=1.
Problems of the type mentioned should not arise if
the point group is Abelian, but will be fairly common in
linear molecules. Since GAMESS cannot detect that the GVB
electronic state is not totally symmetric (or averaged to
at least have a totally symmetric density), it is left up
to you to decide when to input NOSYM=1. If you have any
question about the use of symmetry, try it both ways. If
you get the same energy, both ways, it remains valid to
use symmetry to speed up your run.
And beware! Brain dead computations, such as RHF on
singlet O2, which actually is a half filled degenerate
shell, violate the symmetry assumptions, and also violate
nature. Use of partially filled degenerate shells always
leads to very wild oscillations in the RHF energy, which
is how the program tries to tell you to think first, and
compute second. Configurations such as pi**2, e**1, or
f2u**4 can be treated, but require GVB wavefunctions and
F, ALPHA, BETA values from the sources mentioned.
1
How to do MCSCF and CI calculations
--- -- -- ----- --- -- ------------
On the next pages, you will find older documentation for
the MCSCF program. In the summer of 2002, the changes to
lead to a truly scalable MCSCF program were begun. Since
they are not yet completed, the documentation has not been
brought up to date. In particular, the timing example cited
in PROG.DOC's section on parallel runs reflects the older
code version.
The partial integral transformation used to set up the
FULLNR converger has been changed to use distributed memory.
This will scale like the MP2 energy/gradient program, to
many nodes, and its EXETYP=CHECK should be performed in the
analogous style. The simpler transformation for FOCAS/SOSCF
does not use distributed memory, and in fact the scaling
of this MCSCF option is not very much changed yet. The
FULLNR step builds the orbital hessian in replicated memory
at present, so large MEMORY may be needed as well as MEMDDI.
The GUGA program used to be the only CI step enabled
for parallel execution, with modest scalability. The
determinant programs have been changed to permit replicated
computation on all nodes (zero speedup). Since the ALDET
program is much faster than GUGA, it is actually preferable
to use the non-scaling determinant codes instead of the
modestly scaling GUGA. Total run time will be less!
A chart of the MCSCF options in August 2003:
parallel run
transformation CI computation via CISTEP
converger memory GUGA ALDET GENCI ORMAS
--------- -------------- ---- ----- ----- -----
FOCAS replicated ok ok silly silly
SOSCF replicated ok ok ok ok
FULLNR distributed ok ok ok ok
QUAD serial ok xx xx xx
Jacobi serial ok ok ok ok
"ok" means you can run this computation.
"xx" means QUAD converger is coded only for CISTEP=GUGA.
"silly" means that this converger ignores active-active
rotations, and since most runs with CISTEP=GENCI do
not use a full CI space, these runs are likely to be
divergent, or perhaps converge to a false solution.
"serial" means this can only run sequentially at present.
ORMAS is a full CI within the subspaces, but rotations
between subspaces do affect the energy, so ORMAS is more
likely to require FULLNR over SOSCF than will ALDET.
A sample run, for no symmetry CI calculation with N
electrons in N orbitals should convince most people to use
the determinant CI code CISTEP=ALDET for must full active
space jobs:
N in N ALDET GENCI --- GUGA ---
8 0.8 1.4 0.7 0.5
10 7.9 38.0 19.1 32.6
12 227.5 3122.4 533.9 2208.7
14 7985.2 -- 15376.9 130855.2
The reason there are two numbers under GUGA is that the
first is for writing the loop info to disk (basically
computing H elements) and the second is for the actual
diagonalization. Note that the formation time alone is
greater than the entire ALDET computation, and that ALDET
also has no big disk file holding loops.
If you find something in the older documentation below
that conflicts with this chart, believe the chart.
1
Multi-configuration self consistent field (MCSCF)
wavefunctions are the most general SCF type, offering a
description of chemical processes involving the separation
of electrons (bond breaking, electronically excited states,
etc), which are often not well represented using the single
configuration SCF methods.
MCSCF wavefunctions, as the name implies, contain more
than one configuration, each of which is multiplied by a
"configuration interaction (CI) coefficient", determining
it weight. In addition, the orbitals which form each of
the configurations are optimized, just as in a simpler SCF,
to self consistency.
Typically each chemical problem requires that an MCSCF
wavefunction be designed to treat it, on a case by case
basis. For example, one may be interested in describing
the reactivity of a particular functional group, instead of
elsewhere in the molecule. This means some attention must
be paid in order to obtain correct results.
Procedures for the selection of configurations (which
amounts to choosing the number of active electrons and
active orbitals), for the two mathematical optimizations
just mentioned, ways to interpret the resulting MCSCF
wavefunction, and the treatment for dynamical correlation
not included in the MCSCF wavefunction are the focus of a
recent review article:
"The Construction and Interpretation
of MCSCF wavefunctions"
M.W.Schmidt and M.S.Gordon,
Ann.Rev.Phys.Chem. 49,233-266(1998)
One section of this is devoted to the problem of designing
the correct active space to treat your problem. Additional
reading is listed at the end of this section.
The most efficient technique implemented in GAMESS for
finding the dynamic correlation energy is second order
perturbation theory, in the variant known as MCQDPT.
MCQDPT is discussed in a different section of this chapter.
The use of CI, probably in the form of second order CI,
will be described below, en passant, during discussion of
the input defining the configurations for MCSCF. Selection
of a CI following some type of SCF (except UHF) is made
with CITYP in the $CONTRL group, and masterminded by the
$CIINP group.
1
--- MCSCF implementation ---
With the exception of the QUAD converger, the MCSCF
program is of the type termed "unfolded two step" by Roos.
This means the orbital and CI coefficient optimizations are
separated. The latter are obtained in a conventional CI
diagonalization, while the former are optimized by a
separate orbital improvement step.
Each MCSCF iteration (except for the JACOBI and QUAD
convergers) consists of the following steps:
1) transformation of AO integrals to the current MO basis,
2) generation of the Hamiltonian matrix and optimization
of the CI coefficients by a Davidson diagonalization,
3) generation of the first and second order density matrix,
4) improvement of the molecular orbitals.
The CI problem in steps two and three has three options
for the many electron basis, namely a full determinant or
a selected determinant or a full configuration state
function (CSF) list. The choice of these is determined
by CISTEP in $MCSCF. More will be said just below about
the differences between determinants and CSFs. The word
"configuration" is used in this section to refer to either
when a generic term is needed for the many-electron basis,
so please note there is a distinction between this and the
very similar term CSF.
The orbital problem in step four has four options,
namely FOCAS, SOSCF, FULLNR, and JACOBI, listed here in
order of their increasing mathematical sophistication,
convergence characteristics, and of course, their computer
resource requirements. Again, these are chosen by keywords
in the $MCSCF group. More will be said just below about the
relative merits of these.
Finally, we mention again the QUAD converger, which
works only for a CSF basis, in which the two optimization
problems are treated simultaneously, for modest numbers
of configuratations (50-100 is probably the limit). In
principle, this is the most robust method available, but
in practice, it has not received very much use compared
to the unfolded methods.
Depending on the converger chosen, the program will
select the appropriate kind of integral transformation.
There's seldom need to try to fine tune this, but note
that the $TRANS group does let you pick an AO integral
direct transformation with the DIRTRF flag.
On the first iteration at the first geometry, you will
receive the normal amount of output from each of these
stages, while each subsequent iterations will produce only
a single summarizing line.
1
--- orbital updates ---
There are presently five orbital improvement options,
namely FOCAS, SOSCF, FULLNR, JACOBI, and QUAD. All but the
JACOBI orbital update run in parallel. The convergers are
discussed briefly below, in order of increasing robustness.
FOCAS is a first order, complete active space MCSCF
optimization procedure. The FOCAS code was written by
Michel Dupuis and Antonio Marquez at IBM. It is based on a
novel approach due to Meier and Staemmler, using very fast
but numerous microiterations to improve the convergence of
what is intrinsically a first order method. Since FOCAS
requires only one virtual orbital index in the integral
transformation to compute the orbital gradient (aka the
Lagrangian), the total MCSCF job may take less time than
a second order method, even though it may require many
more iterations to converge. The use of microiterations is
crucial to FOCAS' ability to converge. It is important to
take a great deal of care choosing the starting orbitals.
SOSCF is a method built upon the FOCAS code, which
seeks to combine the speed of FOCAS with second order
convergence properties. Thus SOSCF is an approximate
Newton-Raphson, based on a diagonal guess at the orbital
hessian, and in fact has much in common with the SOSCF
option in $SCF. Its time requirements per iteration are
like FOCAS, with a convergence rate better than FOCAS but
not as good as true second order. Storage of only the
diagonal of the orbital hessian allows the SOSCF method
to be used with much larger basis sets than exact second
order methods. Because it usually requires the least CPU
time, disk space, and memory needs, SOSCF is the default.
Good convergence by the SOSCF method requires that you
prepare starting orbitals carefully, and read in all MOs
in $VEC, as providing canonicalized virtual orbitals
increases the diagonal dominance of the orbital hessian.
FULLNR means a full Newton-Raphson orbital improvement
step is taken, using the exact orbital hessian. FULLNR
is a quite powerful convergence method, and normally takes
the fewest iterations to converge. Computing the exact
orbital hessian requires two virtual orbital indices be
included in the transformation, making this step quite time
consuming, and of course memory for storage of the orbital
hessian must be available. Because both the transformation
and orbital improvement steps of FULLNR are time consuming,
FULLNR is not the default. You may want to try FULLNR when
convergence is difficult, assuming you have already tried
preparing good starting orbitals by the hints below.
1
The FULLNR MCSCF code in GAMESS is adapted from the
HONDO7 program, and was written by Michel Dupuis at IBM.
It uses the the augmented hessian matrix approach to solve
the Newton-Raphson equations. There are two suboptions
for computation of the orbital hessian. DM2 is the fastest
but takes more memory than TEI.
The JACOBI method was written by Joe Ivanic and Klaus
Ruedenberg. It uses a series of 2x2 orbital rotations by
by an angle predicted to lower the energy. This should
essentially ensure convergence after sweeping over all
possible orbital pairs enough times. The procedure was
created to converge selected (general) determinant MCSCF
functions, but of course it can be used will full lists
as well in difficult cases. The JACOBI calculation will
consist of a full four index transformation over all MOs
before the iterations begin. Each iteration consists of
1. a small 4 index transformation over active orbitals
2. optimization of the CI vector
3. generation of the 1e- and 2e- density matrices
4. sweeps over Jacobi rotations, using MO integrals in
memory to generate each rotation, with a subsequent
update after each pair is rotated.
5. when sufficient energy lowering has been achieved,
begin a new iteration.
This procedure never generates the orbital Lagrangian!
Unfortunately this means that at present it is not possible
to compute nuclear gradients. The energy should converge
to a value that is a function of ENGTOL (ACURCY is of
course irrelevant) but actually a bit better than ENGTOL.
QUAD uses a fully quadratic, or second order approach
and is thus the most powerful MCSCF converger. The QUAD
code is also adapted from Michel Dupuis's HONDO. QUAD runs
begin with unfolded FULLNR iterations, until the orbitals
approach convergence sufficiently. QUAD then begins the
simultaneous optimization of CI coefficients and orbitals,
and convergence should be obtained in 3-4 additional MCSCF
iterations. The QUAD method requires building the full
hessian, including orbital/orbital, orbital/CI, and CI/CI
blocks, which is a rather big matrix. QUAD may be helpful
in converging excited electronic states, but note that you
may not use state averaging with QUAD. QUAD is a memory
hog, and so may be used only for fairly small numbers of
configurations.
The input to control the orbital update step is the
$MCSCF group, where you should pick one of the convergence
procedures. Most of the input in this group is rather
specialized, but note in particular MAXIT and ACURCY which
control the convergence behaviour.
1
--- CI coefficient optimization ---
Determinants or configuration state functions (CSFs)
may be used to form the many electron basis set. It is
necessary to explain these in a bit of detail so that you
can understand the advantages of each.
A determinant is a simple object: a product of spin
orbitals with a given Sz quantum number, that is, the
number of alpha spins and number of beta spins are a
constant. Matrix elements involving determinants are
correspondingly simple, but unfortunately determinants
are not necessarily eigenfunctions of the S**2 operator.
To expand on this point, consider the four familiar
2e- functions which satisfy the Pauli principle. Here u,
v are space orbitals, and a, b are the alpha and beta spin
functions. As you know, the singlet and triplets are:
S1 = (uv + vu)/sqrt(2) * (ab - ba)/sqrt(2)
T1 = (uv - vu)/sqrt(2) * aa
T2 = (uv - vu)/sqrt(2) * (ab + ba)/sqrt(2)
T3 = (uv - vu)/sqrt(2) * bb
It is a simple matter to multiply out S1 and T2, and to
expand the two determinants which have Sz=0,
D1 = |ua vb| D2 = |va ub|
This reveals that
S1 = (D1+D2)/sqrt(2) or D1 = (S1 + T2)/sqrt(2)
T2 = (D1-D2)/sqrt(2) D2 = (S1 - T2)/sqrt(2)
Thus, one must take a linear combination of determinants in
order to have a wavefunction with the desired total spin.
There are two important points to note:
a) A two by two Hamiltonian matrix over D1 and D2 has
eigenfunctions with -different- spins, S=0 and S=1.
b) use of all determinants with Sz=0 does allow for the
construction of spin adapted states. D1+D2, or D1-D2,
are -not- spin contaminated.
By itself, a determinant such as D1 is said to be "spin
contaminated", being a fifty-fifty admixture of singlet and
triplet (it is curious that calculations with just one such
determinant are often called "singlet UHF"). Of course,
some determinants are spin adapted all by themselves, for
example the spin adapted functions T1 and T3 above are
single determinants, as are the closed shells
S2 = (uu) * (ab - ba)/sqrt(2).
S3 = (vv) * (ab - ba)/sqrt(2).
It is possible to perform a triplet calculation, with no
singlet states present, by choosing determinants with
Sz=1 such as T1, since then no state with Sz=0 as is
required when S=0 exists in the determinant basis set.
To summarize, the eigenfunctions of a Hamiltonian formed
by determinants with any particular Sz will be spin states
with S=Sz, S=Sz+1, S=Sz+2, ... but will not contain any S
values smaller than Sz.
1
CSFs are an antisymmetrized combination of a space
orbital product, and a spin adapted linear combination of
simple alpha-beta products. Namely, the following CSF
C1 = A (uv) * (ab-ba)/sqrt(2)
which has a singlet spin function is identical to S1 above
if you write out what the antisymmetrizer A does, and the
CSFs
C2 = A (uv) * aa
C3 = A (uv - vu)/sqrt(2) * (ab + ba)/sqrt(2)
C4 = A (uv) * bb
equal T1-T3. Since the three triplet CSFs have the same
energy, GAMESS works with the simpler form C2. Singlet
and triplet computations using CSFs are done in separate
runs, because when spin-orbit coupling is not considered,
the Hamiltonian is block diagonal in a CSF basis.
Technical information about the CSFs is that they use
Yamanouchi-Kotani spin couplings, and matrix elements are
obtained using a GUGA, or graphical unitary group approach.
Both determinant implementations and are primarily used
for MCSCF wavefunctions. The CSF code is capable of more
general CI computations, and so can be used for first or
second order CI computations. Other comparisons between
the determinant and CSF implementations, as they exist in
GAMESS today, are
determinants CSFs
parallel execution no yes
direct CI yes no
uses Abelian space symmetry yes yes
state average mixed spins yes no
first order density yes yes
state averaged densities yes yes
can form CI Lagrangian no yes
In nearly every circumstance the determinant CI will run
faster than GUGA, so it is the default. In addition, the
quality of the initial guess of the CI eigenvector in the
determinant code is much better than in the CSF code, so
the chances of it converging to an incorrect excited state
root is much less.
The next two sections describe in detail the input for
specification of the configurations, either determinants
or CSFs.
1
--- determinant CI ---
Three determinant CI codes are provided for MCSCF, one
for full CI spaces (ALDET), another for Occupation
Restricted Multiple Active Spaces (ORMAS), and one for
arbitrary (selected) determinant lists (GENCI). For
straight CI, but not MCSCF, there is a fourth program, the
full second order CI (FSOCI).
The simple $DET input group is basic to all determinant
CI codes. Keywords GROUP and ISTSYM specify the desired
spatial symmetry of the determinants. Most runs need give
only the orbital and electron counts: NCORE, NACT, and
NELS. The number of electrons is 2*NCORE+NELS, and will
be checked against the charge implied by ICHARG. The MULT
given in $CONTRL is used to determine the desired Sz value,
by extracting S from MULT=2S+1, then by default Sz=S. If
you wish to include lower spin multiplicities, which will
increase the CPU time of the run, but will let you know
what the energies of such states are, just input a smaller
value for SZ. The states whose orbitals will be MCSCF
optimized will be those having the requested MULT value,
unless you choose otherwise with the PURES flag.
The remaining parameters in the $DET group give extra
control over the diagonalization process. Most are not
given in normal circumstances, except NSTATE, which you
may need to adjust to produce enough roots of the desired
MULT value. The only important keyword which has not been
discussed is the WSTATE array, giving the weights for each
state in forming the first and second order density matrix
elements, which drive the orbital update methods. Note
that analytic gradients are available only when the WSTATE
array is a unit vector, corresponding to a pure state, such
as WSTATE(1)=0,1,0 which permits gradients of the first
excited state to be computed. When used for state averaged
MCSCF, WSTATE is normalized to a unit sum, thus input of
WSTATE(1)=1,1,1 really means a weight of 0.33333... for
the each of the states being averaged.
ORMAS is a program designed to limit the size of the
full CI problem, and may be useful when the number of
active orbitals is 10 or higher. By dividing your total
active space into multiple subspaces, and specifying a
range of electrons to occupy each subspace, most of the
full CI's effect can be included. ORMAS generates a full
CI within each orbital subspace, taking the product of
each small full CI to generate the determinant list.
Here are some ideas on how to use ORMAS, which is a very
flexible program:
1
a) single reference, arbitrary excition level CI-X
$det ncore=y nact=z nels=10 (y+z=entire basis)
$ormas nspace=2 mstart(1)=y+1,y+6 mine(1)=10-x,0
maxe(1)=10,x
excites the 5 doubly occupied orbitals, to the
desired excitation level of X.
b) simple product of active spaces
For example, consider furan, with two active
subspaces. Keeping the 5 true core and the 4 CH
bonds in the core space, the sigma subspace might
contains 5 ring sigma, one oxygen lone pair, and 5
ring sigma antibonds, with a total of 12 e-. The pi
active space contains 5 pi orbitals and 6 e-:
$det ncore=9 nact=16 nels=18
$ormas nspace=2 mstart(1)=10,21 mine(1)=12,6
maxe(1)=12,6
Having the minimum and maximum electron counts the
same is what makes this the simple product of two
separate active spaces. In other words, this is
similar to the QCAS procedure of Nakano and Hirao,
but ORMAS limits only the total electron counts,
not separately the numbers of alpha and beta e-,
in other words all spin couplings are used.
c) flexible occupancy between active subspaces
Imagine that you are interested in excited states of
formaldehyde, some of which will have Rydberg
character, dominated by single excitations into
diffuse orbitals. H2CO's valence states arise from 3
orbitals, the CO pi and pi* and one oxygen lone pair.
Placing the O sp lone pair and three sigma bonds into
the filled space, and centering diffuse s,p,d shells
on the carbon:
$det ncore=6 nact=12 nels=4
$ormas nspace=2 mstart(1)=7,10 mine(1)=3,0
maxe(1)=4,1
This is a 4e-, 3 orbital n,pi,pi* space to describe
valence states, and excites one electron into the 9
diffuse orbitals to describe Rydberg states. It is
many fewer determinants than a 4e- in 12 orbital FCI.
d) RAS-like CI
The previous example is reminiscent of Roos' RAS-SCF.
In fact ORMAS can do RAS-SCF, which is three spaces:
the lowest space is allowed to excite only a few
electrons, a middle space that is the rest, and a top
space into which only a few electrons can be excited.
Suppose there are 10 e-, 10 orbitals, that the bottom
and top spaces involve 3 orbitals, and that a "few"
means specifically 2 e-:
$det ncore=20 nact=10 nels=10 $end
$ormas nspace=3 mstart(1)=21,24,28 mine(1)=4,2,0
maxe(1)=6,6,2
However, ORMAS can have more than 3 orbital subspaces.
1
e) first or second order CI.
Consider C2H4, with a 4 orbital active space of CC
sigma, pi, pi*, and sigma*. In order to correlate
the four valence CH orbitals by double excitations,
an MCSCF based on $DET, followed by SOCI based on
$CIDET and $ORMAS, is:
$contrl scftyp=mcscf cityp=ormas
$mcscf cistep=aldet
$det ncore=6 nact=4 nels=4
$cidet ncore=2 nact=y nels=12 (y=rest of basis)
$ormas nspace=3 mstart(1)=3,7,11 mine(1)=6,2,0
maxe(1)=8,6,2
which permits singles and doubles out of the CH and
CC spaces, into the CC and external spaces.
ORMAS is a full CI (or several full CI's) within each
orbital subspace. However, ORMAS does not generate all
excitation levels between spaces (just those implied by the
minimum and maximum electron counts you give). This means
ORMAS MCSCF runs must optimize active-active rotations
between the subspaces, and therefore you should expect
better convergence from FULLNR than SOSCF.
ORMAS is sure to require orbital reordering. For the
furan example just mentioned, there is no reason to
expect that the RHF occupied orbitals will not have the
filled sigma and pi orbitals intermingled. You must use
the NORDER and IORDER keywords in $GUESS to carefully
partition starting orbitals into sigma and pi subspaces.
The selected (general) determinant list is used if
CISTEP=GENCI, and the list is controlled by two input
groups. The first is $GEN, which is identical to $DET
except for inclusion of an additional keyword GLIST=INPUT.
This reads the determinants (as space products) from an
additional input group $GCILST. Completely arbitrary
choices for the space products may be made, but peculiar
lists may lead to poor MCSCF convergence. The FOCAS
converger may not be used, as that assumes full CI spaces.
If you are doing straight CI calculations, the
required input for each determinant CITYP is:
ALDET needs $CIDET
ORMAS needs $CIDET and $ORMAS
GENCI needs $CIDET and $CIGEN and probably $GCILST
FSOCI needs $CIDET and $SODET
In other words, $CIDET replaces $DET, and $CIGEN replaces
$GEN, but the keywords in the group mean the same thing.
The reason for different names is to allow CI calculations
to follow MCSCF in the same run, without clashing input
group names.
1
--- CSF CI ---
The GUGA-based CSF package was originally a set of
different programs, so the input to control it is spread
over several input groups. The CSFs are specified by
a $CIDRT group in the case of CITYP=GUGA, and by a $DRT
group for MCSCF wavefunctions. Thus it is possible to
perform an MCSCF defined by a $DRT input (or perhaps using
$DET during the MCSCF), and follow this with a second order
CI defined by a $CIDRT group, in the same run.
The remaining input groups used by the GUGA CSFs are
$CISORT, $GUGEM, $GUGDIA, and $GUGDM2 for MCSCF runs, with
the latter two being the most important, and in the case
of CI computations, $GUGDM and possibly $LAGRAN groups are
relevant. Perhaps the most interesting variables outside
the $DRT/$CIDRT group are NSTATE in $GUGDIA to include
excited states in the CI computation, IROOT in $GUGDM to
select the state for properties, and WSTATE in $GUGDM2 to
control which (average) state's orbitals are optimized.
The $DRT and $CIDRT groups are almost the same, with
the only difference being orbitals restricted to double
occupancy are called MCC in $DRT, and FZC in $CIDET.
Therefore the rest of this section refers only to "$DRT".
The CSFs are specified by giving a reference CSF,
together with a maximum degree of electron excitation from
that single CSF. The MOs in the reference CSF are filled
in the order MCC or FZC first, followed by DOC, AOS, BOS,
ALP, VAL, and EXT (the Aufbau principle). AOS, BOS, and
ALP are singly occupied MOs. ALP means a high spin alpha
coupling, while AOS/BOS are an alpha/beta coupling to an
open shell singlet. This requires the value NAOS=NBOS,
and their MOs alternate. An example is
NFZC=1 NDOC=2 NAOS=2 NBOS=2 NALP=1 NVAL=3
which gives the reference CSF
FZC,DOC,DOC,AOS,BOS,AOS,BOS,ALP,VAL,VAL,VAL
This is a doublet state with five unpaired electrons. VAL
orbitals are unoccupied only in the reference CSF, they
will become occupied as the other CSFs are generated. This
is done by giving an excitation level, either explicitly by
the IEXCIT variable, or implicitly by the FORS, FOCI, or
SOCI flags. One of these four keywords must be chosen, and
during MCSCF runs, this is usually FORS.
Consider another simpler example, for an MCSCF run,
NMCC=3 NDOC=3 NVAL=2
which gives the reference CSF
MCC,MCC,MCC,DOC,DOC,DOC,VAL,VAL
having six electrons in five active orbitals. Usually,
MCSCF calculations are usually of the Full Optimized
Reaction Space (FORS) type. Some workers refer to FORS
as CASSCF, complete active space SCF. These are the same,
but the keyword is spelled FORS. In the present instance,
choosing FORS=.TRUE. gives an excitation level of 4, as
the 6 valence electrons have only 4 holes available for
1
excitation. MCSCF runs typically have only a small number
of VAL orbitals. It is common to summarize this example
as "six electrons in five orbitals".
The next example is a first or second order multi-
reference CI wavefunction, where
NFZC=3 NDOC=3 NVAL=2 NEXT=-1
leads to the reference CSF
FZC,FZC,FZC,DOC,DOC,DOC,VAL,VAL,EXT,EXT,...
FOCI or SOCI is chosen by selecting the appropriate flag,
the correct excitation level is automatically generated.
Note that the negative one for NEXT causes all remaining
MOs to be included in the external orbital space. One way
of viewing FOCI and SOCI wavefunctions is as all singles,
or all singles and doubles, from the entire MCSCF wave-
function as a reference. An equivalent way of saying this
is that all CSFs with N electrons (in this case N=6)
distributed in the valence orbitals in all ways (that is
the FORS MCSCF wavefunction) to make the base wavefunction.
To this, FOCI adds all CSFs with N-1 electrons in active
and 1 electron in external orbitals. SOCI adds all CSFs
with N-2 electrons in active orbitals and 2 in external
orbitals. SOCI is often prohibitively large, but is also
a very accurate wavefunction.
Sometimes people use the CI package for ordinary
single reference CI calculations, such as
NFZC=3 NDOC=5 NVAL=34
which means the reference RHF wavefunction is
FZC FZC FZC DOC DOC DOC VAL VAL ... VAL
and in this case NVAL is a large number conveying the
total number of -virtual- orbitals into which electrons
are excited. The excitation level would be given as
IEXCIT=2, perhaps, to perform a SD-CI. All excitations
smaller than the value of IEXCIT are automatically
included in the CI. Note that NVAL's spelling was chosen
to make the most sense for MCSCF calculations, and so it
is a bit of a misnomer here.
Before going on, there is a quirk related to single
reference CI that should be mentioned. Whenever the
single reference contains unpaired electrons, such as
NFZC=3 NDOC=4 NALP=2 NVAL=33
some "extra" CSFs will be generated. The reference here
can be abbreviated
2222 11 000 000 000 000 000 000 000 000 000 000 000
Supposing IEXCIT=2, the following CSF
2200 22 000 011 000 000 000 000 000 000 000 000 000
will be generated and used in the CI. Most people would
prefer to think of this as a quadruple excitation from
the reference, but acting solely on the reasoning that
no more than two electrons went into previously vacant
NVAL orbitals, the GUGA CSF package decides it is a double.
So, an open shell SD-CI calculation with GAMESS will not
give the same result as other programs, although the result
1
for any such calculation with these "extras" is correctly
computed. Note that if you also select the INTACT option,
the extra space products are eliminated, but that some of
the spin couplings for the truly IEXCIT'd space products
are also eliminated.
As was discussed above, the CSFs are automatically
spin-symmetry adapted, with S implicit in the reference
CSF. The spin quantum number you appear to be requesting
in $DRT (basically, S = NALP/2) will be checked against
the value of MULT in $CONTRL, and the total number of
electrons, 2*NMCC(or NFZC) + 2*NDOC + NAOS + NBOS + NALP
will be checked against the input given for ICHARG.
The CSF package is also able to exploit spatial
symmetry, which like the spin and charge, is implicitly
determined by the choice of the reference CSF. The keyword
GROUP in $DRT governs the use of spatial symmetry.
The CSF program works with Abelian point groups, which
are D2h and any of its subgroups. However, $DRT allows
the input of some (but not all) higher point groups. For
non-Abelian groups, the program automatically assigns the
orbitals to an irrep in the highest possible Abelian
subgroup. For the other non-Abelian groups, you must at
present select GROUP=C1. Note that when you are computing
a Hessian matrix, many of the displaced geometries are
asymmetric, hence you must choose C1 in $DRT (however, be
sure to use the highest symmetry possible in $DATA!).
The symmetry of the reference CSF given in your $DRT
determines the symmetry of the CSFs which are generated.
As an example, consider a molecule with Cs symmetry, and
these two reference CSFs
...MCC...DOC DOC VAL VAL
...MCC...DOC AOS BOS VAL
Suppose that the 2nd and 3rd active MOs have symmetries a'
and a". Both of these generate singlet wavefunctions,
with 4 electrons in 4 active orbitals, but the former
constructs 1-A' CSFs, while the latter generates 1-A"
CSFs. However, if the 2nd and 3rd orbitals have the same
symmetry type, an identical list of CSFs is generated.
In cases with high point group symmetry, it may be
possible to generate correct state degeneracies only by
using no symmetry (GROUP=C1) when generating CSFs. As
an example, consider the 2-pi ground state of NO. If you
use GROUP=C4V, which will be mapped into its highest
Abelian subgroup C2v, the two components of the pi state
will be seen as belonging to different irreps, B1 and B2.
The only way to ensure that both sets of CSFs are generated
is to enforce no symmetry at all, so that CSFs for both
components of the pi level are generated. This permits
state averaging (WSTATE(1)=0.5,0.5) to preserve cylindrical
symmetry. It is however perfectly feasible to use C4v or
D4h symmetry in $DRT when treating sigma states.
1
The use of spatial symmetry decreases the number of
CSFs, and thus the size of the Hamiltonian that must be
computed. In molecules with high symmetry, this may lead
to faster run times with the GUGA CSF code, compared to
the determinant code.
--- starting orbitals ---
The first step is to partition the orbital space into
core, active, and external sets, in a manner which is
sensible for your chemical problem. This is a bit of an
art, and the user is referred to the references quoted at
the end of this section. Having decided what MCSCF to
perform, you now must consider the more pedantic problem
of what orbitals to begin the MCSCF calculation with.
You should always start an MCSCF run with orbitals
from some other run, by means of GUESS=MOREAD. Do not
expect to be able to use HCORE or HUCKEL! Example 6 is a
poor example, converging only because methylene has so much
symmetry, and the basis is so small. If you are beginning
your MCSCF problem, use orbitals from some appropriate
converged SCF run. Thus, a realistic example of an MCSCF
calculation is examples 8 and 9. Once you get an MCSCF
to converge, you can and should use these MCSCF MOs (which
will be Schmidt orthogonalized) at other nearby geometries.
Starting from SCF orbitals can take a little bit of
care. Most of the time (but not always) the orbitals you
want to correlate will be the highest occupied orbitals in
the SCF. Fairly often, however, the correlating orbitals
you wish to use will not be the lowest unoccupied virtuals
of the SCF. You will soon become familiar with NORDER=1
in $GUESS, as reordering is needed in 50% or more cases.
The occupied and especially the virtual canonical SCF
MOs are often spread out over regions of the molecule other
than "where the action is". Orbitals which remedy this can
generated by two additional options at almost no CPU cost.
One way to improve upon the SCF orbitals as starting
MOs is to generate modified virtual orbitals (MVOs).
MVOs are obtained by diagonalizing the Fock operator of a
very positive ion, within the virtual orbital space only.
As implemented in GAMESS, MVOs can be obtained at the end
of any RHF, ROHF, or GVB run by setting MVOQ in $SCF
nonzero, at the cost of a single SCF cycle. Typically, we
use MVOQ=+6. Generating MVOs does not change any of the
occupied SCF orbitals of the original neutral, but gives
more valence-like LUMOs.
Another way to improve SCF starting orbitals is by
a partial localization of the occupied orbitals. Typically
MCSCF active orbitals are concentrated in the part of the
1
molecule where bonds are breaking, etc. Canonical SCF MOs
are normally more spread out. By choosing LOCAL=BOYS along
with SYMLOC=.TRUE. in $LOCAL, you can get orbitals which
are localized, but still retain orbital symmetry to help
speed the MCSCF along. In groups with an inversion center,
a SYMLOC Boys localization does not change the orbitals,
but you can instead use LOCAL=POP. Localization tends to
order the orbitals fairly randomly, so be prepared to
reorder them appropriately.
Pasting the virtuals from a MVOQ run onto the occupied
orbitals of a SYMLOC run (both can be done in the same SCF
computation) gives the best possible set of starting
orbitals. If you also take the time to design your active
space carefully, select the appropriate starting orbitals
from this combined $VEC, and inspect your converged results,
you will be able to carry out MCSCF computations correctly.
Convergence of MCSCF is by no means guaranteed. Poor
convergence can invariably be traced back to either a poor
initial selection of orbitals, or a poorly chosen number of
active orbitals. My advice is, before you even start:
"Look at the orbitals.
Then look at the orbitals again".
Later, if you have any trouble:
"Look at the orbitals some more".
Few people are able to see the orbital shapes in the LCAO
matrix in a log file, and so need a visualization program.
If you have a Macintosh, download a copy of MacMolPlt from
http://www.msg.ameslab.gov/GAMESS/GAMESS.html
for 2D or 3D plots, or use PLTORB under X-windows for 2D.
Even if you don't have any trouble, look at the
orbitals to see if they converged to what you expected,
and have reasonable occupation numbers. It is particularly
useful to check the oriented localized MCSCF orbitals (see
the discussion of this in the section on localized orbitals
in this section for more information). MCSCF is by no meand
the sort of "black box" that RHF is these days, so please
look very carefully at your final results.
--- miscellaneous hints ---
It is very helpful to execute a EXETYP=CHECK run
before doing any MCSCF or CI run. The CHECK run will tell
you the total number of configurations and check the charge
and multiplicity and electronic state symmetry, based on
your input. The CHECK run also lets the program feel out
the memory that will be required to actually do the run.
Thus the CHECK run can potentially prevent costly mistakes,
or tell you when a calculation is prohibitively large.
1
A very common MCSCF wavefunction has 2 electrons in 2
active MOs. This is the simplest possible wavefunction
describing a singlet diradical. While this function can be
obtained in an MCSCF run (using NACT=2 NELS=2 or NDOC=1
NVAL=1), it can be obtained much faster by use of the GVB
code, with one GVB pair. This GVB-PP(1) wavefunction is
also known in the literature as two configuration SCF, or
TCSCF. The two configurations of this GVB are equivalent
to the three configurations used in this MCSCF, as orbital
optimization in natural form (configurations 20 and 02)
causes the coefficient of the 11 configuration to vanish.
If you are using many GUGA CSFs (say 150,000 or more)
the main bottleneck in the MCSCF calculation is the
formation and diagonalization of the Hamiltonian, not the
integral transformation and orbital improvement steps.
In this case, you would be wise to switch to FULLNR, which
will minimize the total number of iterations. In addition,
each orbital improvement may contain some microiterations,
which consists of an integral transformation over the new
MOs, followed immediately by a orbital improvement, reusing
the current 2nd order density matrix. MICIT=2 in $MCSCF
may lead to better overall run times by doing two orbital
updates for every CI diagonalization step.
Since the determinant CI is a direct CI, it does not
have the bottleneck of storing a large disk file containing
Hamiltonian information. However, very large active spaces
containing 13 or 14 orbitals may result in more time being
spent in CI iterations than in the rest of the MCSCF steps.
The analogous trick to MICIT=2 is ITERMX in $DET or $GEN,
which may be set to a value like ITERMX=2 or ITERMX=3 to
improve the CI vectors only a bit. Since each iteration's
CI calculation starts with the previous iterations result,
the CI step will become fully converged during the MCSCF
cycles. The total run time may decrease, although there
may be a few additional MCSCF iterations required. For
small active spaces where the CI step takes trivial time,
you should use a bigger ITERMX to ensure fully converged
CI states are generated on every iteration.
If you choose to use ORMAS, a general determinant CI,
or if you select an excitation level IEXCIT smaller than
that needed to generate the FORS space, you must use the
SOSCF, JACOBI, or FULLNR method as these can optimize
active-active rotations. Be sure to set FORS=.FALSE. in
$MCSCF when for non-full CI cases, or else very poor
convergence will result. Actually, the convergence for
incomplete active spaces is likely to be poorer than for
full active spaces, anyway.
1
--- references ---
There are several review articles about MCSCF listed
below. Of these, the first two are a nice overview of the
subject, the final 3 are more technical.
1. "The Construction and Interpretation of MCSCF
wavefunctions"
M.W.Schmidt and M.S.Gordon,
Ann.Rev.Phys.Chem. 49,233-266(1998)
2a. "The Multiconfiguration SCF Method"
B.O.Roos, in "Methods in Computational Molecular
Physics", edited by G.H.F.Diercksen and S.Wilson
D.Reidel Publishing, Dordrecht, Netherlands,
1983, pp 161-187.
2b. "The Multiconfiguration SCF Method"
B.O.Roos, in "Lecture Notes in Quantum Chemistry",
edited by B.O.Roos, Lecture Notes in Chemistry v58,
Springer-Verlag, Berlin, 1994, pp 177-254.
3. "Optimization and Characterization of a MCSCF State"
J.Olsen, D.L.Yeager, P.Jorgensen
Adv.Chem.Phys. 54, 1-176(1983).
4. "Matrix Formulated Direct MCSCF and Multiconfiguration
Reference CI Methods"
H.-J.Werner, Adv.Chem.Phys. 69, 1-62(1987).
5. "The MCSCF Method"
R.Shepard, Adv.Chem.Phys. 69, 63-200(1987).
There is an entire section on the choice of active
spaces in Reference 1. As this is a matter of great
importance, here are two alternate presentations of the
design of active spaces:
6. "The CASSCF Method and its Application in Electronic
Structure Calculations"
B.O.Roos, in "Advances in Chemical Physics", vol.69,
edited by K.P.Lawley, Wiley Interscience, New York,
1987, pp 339-445.
7. "Are Atoms Intrinsic to Molecular Electronic
Wavefunctions?"
K.Ruedenberg, M.W.Schmidt, M.M.Gilbert, S.T.Elbert
Chem.Phys. 71, 41-49, 51-64, 65-78 (1982).
Two papers germane to the FOCAS implementation are
8. "An Efficient first-order CASSCF method based on
the renormalized Fock-operator technique."
U.Meier, V.Staemmler Theor.Chim.Acta 76, 95-111(1989)
9. "Modern tools for including electron correlation in
electronic structure studies"
M.Dupuis, S.Chen, A.Marquez, in "Relativistic and
Electron Correlation Effects in Molecules and
Solids", edited by G.L.Malli, Plenum, NY 1994
1
The paper germane to the the SOSCF method is
10. "Approximate second order method for orbital
optimization of SCF and MCSCF wavefunctions"
G.Chaban, M.W.Schmidt, M.S.Gordon
Theor.Chem.Acc. 97: 88-95(1997)
Two papers germane to the FULLNR implementation, and
one discussing the implementation details are
11. "General second order MCSCF theory: A Density Matrix
Directed Algorithm"
B.H.Lengsfield, III, J.Chem.Phys. 73,382-390(1980).
12. "The use of the Augmented Matrix in MCSCF Theory"
D.R.Yarkony, Chem.Phys.Lett. 77,634-635(1981).
13. M.Dupuis, P.Mougenot, J.D.Watts, in "Modern Techniques
in Theoretical Chemistry", E.Clementi, editor, ESCOM,
Leiden, 1989, chapter 7.
The paper describing the JACOBI converger is
14. "A MCSCF method for ground and excited states based on
full optimizatons of successive Jacobi rotations"
J.Ivanic, K.Ruedenberg J.Comput.Chem. 24, 1250-1262(2003)
For determinants and CSFs, respectively, see
15. "Identification of deadwood in configuration spaces
through general direct configuration interaction"
J.Ivanic, K.Ruedenberg
Theoret.Chem.Acc. 106, 339-351(2001)
16. "The GUGA approach to the electron correlation problem"
B.R.Brooks, H.F.Schaefer
J.Chem.Phys. 70, 5092-5106(1979)
The final references are simply some examples of FORS
MCSCF applications, the latter done with GAMESS.
16. D.F.Feller, M.W.Schmidt, K.Ruedenberg,
J.Am.Chem.Soc. 104, 960-967(1982).
17. M.W.Schmidt, P.N.Truong, M.S.Gordon,
J.Am.Chem.Soc. 109, 5217-5227(1987).
1
Second order perturbation theory
The perturbation theory techniques available in GAMESS
expand to the second order energy correction only, but
permit use of a broad range of zeroth order wavefunctions.
Since MP2 theory for systems well described as closed
shells recovers only about 80% of the correlation energy
(assuming the use of large basis sets), it is common to
extend the perturbative treatment to higher order, or to
use coupled cluster theory. While this is reasonable for
systems well described by RHF or UHF with small spin
contamination, this is probably a poor approach when the
system is not well described at zeroth order by these wave-
functions.
The input for second order pertubation calculations
based on SCFTYP=RHF, UHF, or ROHF is found in $MP2, while
for SCFTYP=MCSCF, see $MCQDPT.
--- RHF and UHF MP2
These methods are well defined, due to the uniqueness
of the Fock matrix definitions. These methods are also
well understood, and there is little need to say more.
One point which may not be commonly appreciated is that
the density matrix for the first order wavefunction for the
RHF and UHF case, which is generated during gradient runs
or if properties are requested in the $MP2 group, is of the
type known as "response density", which differs from the
more usual "expectation value density". The eigenvalues
of the response density matrix (which are the occupation
numbers of the MP2 natural orbitals) can therefore be
greater than 2 for frozen core orbitals, or even negative
values for the highest 'virtual' orbitals. The sum is
of course exactly the total number of electrons. We have
seen values outside the range 0-2 in several cases where
the single configuration RHF wavefunction was not an
appropriate description of the system, and thus these
occupancies may serve as a guide to the wisdom of using
a RHF reference. See
M.S.Gordon, M.W.Schmidt, G.M.Chaban, K.R.Glaesemann,
W.J.Stevens, C.Gonzalez J.Chem.Phys. 110,4199-4207(1999)
By default, frozen core MP2 calculations are performed.
1
--- high spin ROHF MP2
There are a number of open shell perturbation theories
described in the literature. It is important to note that
these methods give different results for the second order
energy correction, reflecting ambiguities in the selection
of the zeroth order Hamiltonian and in defining the ROHF
Fock matrices. Two of these are available in GAMESS.
One theory is known as RMP, which it should be pointed
out, is entirely equivalent to the ROHF-MBPT2 method. The
theory is as UHF-like as possible, and can be chosen in
GAMESS by selection of OSPT=RMP in $MP2. The second order
energy is defined by
1. P.J.Knowles, J.S.Andrews, R.D.Amos, N.C.Handy,
J.A.Pople Chem.Phys.Lett. 186, 130-136(1991)
2. W.J.Lauderdale, J.F.Stanton, J.Gauss, J.D.Watts,
R.J.Bartlett Chem.Phys.Lett. 187, 21-28(1991).
The submission dates are in inverse order of publication
dates, and -both- papers should be cited when using this
method. Here we will refer to the method as RMP in
keeping with much of the literature. The RMP method
diagonalizes the alpha and beta Fock matrices separately,
so that their occupied-occupied and virtual-virtual blocks
are canonicalized. This generates two distinct orbital
sets, whose double excitation contributions are processed
by the usual UHF MP2 program, but an additional energy
term from single excitations is required.
RMP's use of different orbitals for different spins adds
to the CPU time required for integral transformations, of
course. RMP is invariant under all of the orbital
transformations for which the ROHF itself is invariant.
Unlike UMP2, the second order RMP energy does not suffer
from spin contamination, since the reference ROHF wave-
function has no spin contamination. The RMP wavefunction,
however, is spin contaminated at 1st and higher order,
and therefore the 3rd and higher order RMP energies are
spin contaminated. Other workers have extended the RMP
theory to gradients and hessians at second order, and to
fourth order in the energy,
3. W.J.Lauderdale, J.F.Stanton, J.Gauss, J.D.Watts,
R.J.Bartlett J.Chem.Phys. 97, 6606-6620(1992)
4. J.Gauss, J.F.Stanton, R.J.Bartlett
J.Chem.Phys. 97, 7825-7828(1992)
5. D.J.Tozer, J.S.Andrews, R.D.Amos, N.C.Handy
Chem.Phys.Lett. 199, 229-236(1992)
6. D.J.Tozer, N.C.Handy, R.D.Amos, J.A.Pople, R.H.Nobes,
Y.Xie, H.F.Schaefer Mol.Phys. 79, 777-793(1993)
We deliberately omit references to the ROMP precurser to
the RMP formalism.
1
The ZAPT formalism is also implemented in GAMESS, as
OSPT=ZAPT in $MP2. Because this theory is not spin-
contaminated at any order, and has only a single set of
orbitals in the MO transformation, it is the default.
References for ZAPT (Z-averaged perturbation theory) are
7. T.J.Lee, D.Jayatilaka Chem.Phys.Lett. 201, 1-10(1993)
8. T.J.Lee, A.P.Rendell, K.G.Dyall, D.Jayatilaka
J.Chem.Phys. 100, 7400-7409(1994)
The formulae for the seven terms in the energy are most
clearly summarized in the paper
9. I.M.B.Nielsen, E.T.Seidl
J.Comput.Chem. 16, 1301-1313(1995)
We would like to thank Tim Lee for very gracious assistance
in the implementation of ZAPT.
There are a number of other open shell theories, with
names such as HC, OPT1, OPT2, and IOPT. The literature
for these is
10. I.Hubac, P.Carsky Phys.Rev.A 22, 2392-2399(1980)
11. C.Murray, E.R.Davidson
Chem.Phys.Lett. 187,451-454(1991)
12. C.Murray, E.R.Davidson
Int.J.Quantum Chem. 43, 755-768(1992)
13. P.M.Kozlowski, E.R.Davidson
Chem.Phys.Lett. 226, 440-446(1994)
14. C.W.Murray, N.C.Handy
J.Chem.Phys. 97, 6509-6516(1992)
15. T.D.Crawford, H.F.Schaefer, T.J.Lee
J.Chem.Phys. 105, 1060-1069(1996)
The latter two of these give comparisons of the various
high spin methods, and the numerical results in ref. 15
are the basis for the conventional wisdom that restricted
open shell theory is better convergent with order of the
perturbation level than unrestricted theory. Paper 8 has
some numerical comparisons of spin-restricted theories
as well. We are aware of one paper on low-spin coupled
open shell SCF perturbation theory
16. J.S.Andrews, C.W.Murray, N.C.Handy
Chem.Phys.Lett. 201, 458-464(1993)
but this is not implemented in GAMESS. See the MCQDPT
code for cases such as this.
1
--- GVB based MP2
This is not implemented in GAMESS. Note that the MCSCF
MP2 program discussed below should be able to develop the
perturbation correction for open shell singlets, by using
a $DRT input such as
NMCC=N/2-1 NDOC=0 NAOS=1 NBOS=1 NVAL=0
which generates a single CSF if the two open shells have
different symmetry, or for a one pair GVB function
NMCC=N/2-1 NDOC=1 NVAL=1
which generates a 3 CSF function entirely equivalent to
the two configuration TCSCF, a.k.a GVB-PP(1). For the
record, we note that if we attempt a triplet state with
the MCSCF program,
NMCC=N/2-1 NDOC=0 NALP=2 NVAL=0
we get a result equivalent to the OPT1 open shell method
described above, not the RMP result. It is possible to
generate the orbitals with a simpler SCF computation than
the MCSCF $DRT examples just given, and read them into the
MCSCF based MP2 program described below, by INORB=1.
--- MCSCF based MP2
Just as for the open shell case, there are several ways
to define a multireference perturbation theory. The most
noteworthy are the CASPT2 method of Roos' group, the MRMP2
method of Hirao, the MROPTn methods of Davidson, and the
MCQDPT2 method of Nakano. Although the results of each
method are different, energy differences should be rather
similar. In particular, the MCQDPT2 method implemented in
GAMESS gives results for the singlet-triplet splitting of
methylene in close agreement to CASPT2, MRMP2(Fav), and
MROPT1, and differs by 2 Kcal/mole from MRMP2(Fhs), and
the MROPT2 to MROPT4 methods.
The MCQDPT method implemented in GAMESS is a multistate
perturbation theory. If applied to 1 state, it is the same
as the MRMP model of Hirao. When applied to more than one
state, it is of the philosophy "perturb first, diagonalize
second". This means that perturbations are made to both
the diagonal and offdiagonal elements of an effective
Hamiltonian, whose dimension equals the number of states
being treated. The perturbed Hamiltonian is diagonalized
to give the corrected state energies. Diagonalization
after inclusion of the offdiagonal perturbation ensures
that avoided crossings of states of the same symmetry are
treated correctly. Such an avoided crossing is found in
the LiF molecule, as shown in the first of the two papers
on the MCQDPT method:
H.Nakano, J.Chem.Phys. 99, 7983-7992(1993)
H.Nakano, Chem.Phys.Lett. 207, 372-378(1993)
1
The closely related single state "diagonalize, then perturb"
MRMP model is discussed by
K.Hirao, Chem.Phys.Lett. 190, 374-380(1992)
K.Hirao, Chem.Phys.Lett. 196, 397-403(1992)
K.Hirao, Int.J.Quant.Chem. S26, 517-526(1992)
K.Hirao, Chem.Phys.Lett. 201, 59-66(1993)
Computation of reference weights and energy contributions
is illustrated by
H.Nakano, K.Nakayama, K.Hirao, M.Dupuis
J.Chem.Phys. 106, 4912-4917(1997)
T.Hashimoto, H.Nakano, K.Hirao
J.Mol.Struct.(THEOCHEM) 451, 25-33(1998)
Single state MCQDPT computations are very similiar to MRMP
computations. A beginning set of references to the other
multireference methods used includes:
P.M.Kozlowski, E.R.Davidson
J.Chem.Phys. 100, 3672-3682(1994)
K.G.Dyall J.Chem.Phys. 102, 4909-4918(1995)
B.O.Roos, K.Andersson, M.K.Fulscher, P.-A.Malmqvist,
L.Serrano-Andres, K.Pierloot, M.Merchan
Adv.Chem.Phys. 93, 219-331(1996).
and a review article is available comparing these methods,
E.R.Davidson, A.A.Jarzecki in "Recent Advances in Multi-
reference Methods" K.Hirao, Ed. World Scientific, 1999,
pp 31-63.
The MCQDPT code was written by Haruyuki Nakano, and was
interfaced to GAMESS by him in the summer of 1996. After
a few months experience, we can say that this code seems to
run in memory, disk, and CPU time comparable to the MCSCF
computation itself. It can be used for 150 to 250 AOs, for
example. A 2001 calculation with 351 AOs, 116 cores, and 8
active orbitals and 8 electrons was more heroic, requiring
an attached disk subsystem of 100 GBytes. Efficiency is
improved if you can add extra physical memory to reduce the
number of file reads.
We close the discussion with an input example which
illustrates RMP and MCQDPT computations on the ground state
of NH2 radical:
! 2nd order perturbation test on NH2, following
! T.J.Lee, A.P.Rendell, K.G.Dyall, D.Jayatilaka
! J.Chem.Phys. 100, 7400-7409(1994), Table III.
! State is 2-B-1, 69 AOs, 49 CSFs.
!
! For 1 CSF reference,
! E(ROHF) = -55.5836109825
! E(RMP) = -55.7772299929 (lit. RMP = -75.777230)
! E(MCQDPT) = -55.7830423024 (lit. OPT1= -75.783044)
! [E(MCQDPT) = -55.7830437413 at the lit's OPT1 geometry]
!
! For 49 CSF reference,
! E(MCSCF) = -55.6323324949
! E(MCQDPT) = -55.7857458575
!
1
$contrl scftyp=mcscf mplevl=2 runtyp=energy mult=2 $end
$system timlim=60 memory=1000000 $end
$guess guess=moread norb=69 $end
$mcscf fullnr=.true. $end
! Next two lines carry out a MCQDPT computation, after
! carrying out a full valence MCSCF orbital optimization.
$drt group=c2v fors=.t. nmcc=1 ndoc=3 nalp=1 nval=2 $end
$mcqdpt inorb=0 mult=2 nmofzc=1 nmodoc=0 nmoact=6
istsym=3 nstate=1 $end
! Next two lines carry out a single reference computation,
! using converged ROHF orbitals from the $VEC.
--- $drt group=c2v fors=.t. nmcc=4 ndoc=0 nalp=1 nval=0 $end
--- $mcqdpt inorb=1 nmofzc=1 nmodoc=3 nmoact=1
--- istsym=3 nstate=1 $end
$data
NH2...2-B-1...TZ2Pf basis, RMP geom. used by Lee, et al.
Cnv 2
Nitrogen 7.0
S 6
1 13520.0 0.000760
2 1999.0 0.006076
3 440.0 0.032847
4 120.9 0.132396
5 38.47 0.393261
6 13.46 0.546339
S 2
1 13.46 0.252036
2 4.993 0.779385
S 1 ; 1 1.569 1.0
S 1 ; 1 0.5800 1.0
S 1 ; 1 0.1923 1.0
P 3
1 35.91 0.040319
2 8.480 0.243602
3 2.706 0.805968
P 1 ; 1 0.9921 1.0
P 1 ; 1 0.3727 1.0
P 1 ; 1 0.1346 1.0
D 1 ; 1 1.654 1.0
D 1 ; 1 0.469 1.0
F 1 ; 1 1.093 1.0
Hydrogen 1.0 0.0 0.7993787 0.6359684
S 3 ! note that this is unscaled
1 33.64 0.025374
2 5.058 0.189684
3 1.147 0.852933
S 1 ; 1 0.3211 1.0
S 1 ; 1 0.1013 1.0
P 1 ; 1 1.407 1.0
P 1 ; 1 0.388 1.0
D 1 ; 1 1.057 1.0
$end
E(ROHF)= -55.5836109825, E(NUC)= 7.5835449477, 9 ITERS
$VEC ...omitted... $END
1
Coupled-Cluster theory
The coupled-cluster (CC) theory is widely recognized as
one of the most accurate methods for describing ground
electronic states of atoms and molecules. CC approaches
provide the best compromise between relatively low computer
costs and high accuracy. For example, the popular CCSD(T)
approach, which is a No**2 * Nu**4 procedure in the
iterative CCSD steps and a No**3 * Nu**4 procedure in the
non-iterative steps related to the calculation of triples
corrections, is capable of providing results of the CISDTQ
quality (which is a iterative No**4 * Nu**6 procedure).
Here, No and Nu are the numbers of correlated occupied and
unoccupied orbitals. Unlike CI methods, all standard CC
methods provide size extensive descriptions of molecular
systems.
Thanks to numerous advances in both the formal aspects
of the CC theory and the development of efficient computer
codes, the single-reference CC approach is nowadays
routinely used in calculations for closed-shell and simple
open-shell electronic ground states of atomic and molecular
systems. Extensions of the CC theory to quasi-degenerate
and excited states are possible, via the multi-reference,
equation-of-motion, and response formalisms, and it is
expected that at least some of these extensions will become
at least as popular as the existing multi-reference CI or
CASSCF methods. We should also add that the CC theory is
a fundamental many-body formalism, whose applicability
ranges from electronic structure of atoms and molecules and
nuclear physics to extended systems, phase transitions,
condensed matter theory, theories of homogeneous electron
gas, and relativistic quantum field theory, to mention a
few examples.
A number of review articles have been written over the
years and it is difficult to cite all of them here. We
recommend that users of GAMESS planning to use CC methods
read one or more reviews listed below:
"Coupled-cluster theory"
J. Paldus, in S. Wilson and G.H.F. Diercksen (Eds.),
Methods in Computational Molecular Physics, NATO Advanced
Study Institute, Series B: Physics, Vol. 293, Plenum, New
York, 1992, p. 99-194.
"Applications of post-Hartree-Fock methods: a tutorial."
R.J. Bartlett and J.F. Stanton, in K.B. Lipkowitz and
D.B.Boyd (Eds.), Reviews in Computational Chemistry,
Vol. 5, VCH Publishers, New York, 1994, p. 65-169.
"Coupled-Cluster Theory: an Overview of Recent Developments"
R.J. Bartlett, in D.R. Yarkony (Ed.), Modern Electronic
Structure Theory, Part I, World Scientific, Singapore,
1995, p. 1047-1131.
"Achieving chemical accuracy with coupled-cluster theory"
T.J. Lee and G.E. Scuseria, in S.R. Langhoff (Ed.),
Quantum Mechanical Electronic Structure Calculations with
1
Chemical Accuracy, Kluwer, Dordrecht, The Netherlands,
1995, p. 47-108.
"A Critical Assessment of Coupled Cluster Method in Quantum
Chemistry"
J. Paldus and X. Li, Adv. Chem. Phys. 110, 1 (1999),
T.D. Crawford and H.F. Schaefer III,
"An Introduction to Coupled Cluster Theory for Computational
Chemists"
T.D.Crawford, H.F.Schaefer in K.B. Lipkowitz and D.B.Boyd
(Eds.), Reviews in Computational Chemistry, Vol. 14, VCH
Publishers, New York, 2000, p. 33-136.
These reviews point to the other review articles and many
original papers. The list of original papers relevant to
CC methods implemented in GAMESS is provided below.
* * * available computations * * *
The CC program incorporated in GAMESS enables the user
to perform the standard LCCD, CCD, CCSD, CCSD[T] (also
known as CCSD+T(CCSD)) and CCSD(T) calculations and the
renormalized (R) and completely renormalized (CR) CCSD[T]
and CCSD(T) calculations for closed-shell RHF references.
Performance of the standard CC methods has been discussed
in a number of places (cf. the review articles mentioned
above). Methods such as CCSD(T) provide excellent results
for molecules in or near their equilibrium geometries.
Almost all standard CC methods are very good in describing
dynamical correlation, while being relatively inexpensive
and easy to use. We must remember, however, that they
should not be applied to bond breaking and quasi-degenerate
states when the RHF determinant is used as a reference.
The unique features of the CC code in GAMESS are the
renormalized (R) and completely renormalized (CR) CCSD[T]
and CCSD(T) methods, which are based on the more general
formalism of the method of moments of coupled-cluster
equations (MMCC). These new methods remove the pervasive
failing of the standard CCSD[T] and CCSD(T) approximations
at larger internuclear separations, while preserving the
simplicity and the relatively low cost of the standard
methods. In analogy to the CCSD[T] and CCSD(T) methods,
the R-CCSD[T], R-CCSD(T), CR-CCSD[T], and CR-CCSD(T)
approaches are based on an idea of improving the CCSD
results by adding a posteriori noniterative corrections to
CCSD energies. The CR-CCSD[T] and CR-CCSD(T) approaches
eliminate the unphysical humps on the potential energy
surfaces involving single bond breaking produced by the
standard CCSD[T] and CCSD(T) methods. The R-CCSD[T] and
R-CCSD(T) approaches may improve the standard CCSD[T] and
CCSD(T) results at intermediate internuclear separations,
but they usually fail at larger distances. The CR-CCSD[T]
and CR-CCSD(T) methods are much better in this regard,
since they provide a very good description of single bond
breaking at all internuclear separations. This includes
various cases of unimolecular dissociations and exchange
1
chemical reactions, in which single bonds break and form.
We do not recommend applying the CR-CCSD[T] and CR-CCSD(T)
approaches to multiple bond breaking (in this case, one
should resort to the completely renormalized CCSD(TQ) and
CCSDT(Q) approaches or to the so-called MMCC(2,6) method;
these approaches will be included in GAMESS in the future).
A detailed description of the R-CCSD[T], R-CCSD(T),
CR-CCSD[T], and CR-CCSD(T) and other MMCC methods can be
found in the papers at the very end of this subsection by
Piecuch, Kucharski, Kowalski, et al.
The R-CCSD[T], R-CCSD(T), CR-CCSD[T], and CR-CCSD(T)
methods are not strictly size extensive: there are unlinked
terms in the MBPT (many-body perturbation theory) expansions
of the renormalized and completely renormalized [T] and (T)
corrections to CCSD energies. This has no effect on bond
breaking (on the contrary, the CR-CCSD[T] and CR-CCSD(T)
potential surfaces are MUCH better than potential energy
surfaces obtained in the standard CCSD[T] and CCSD(T)
calculations), but lack of strict size extensivity may have
an effect on the results of calculations for very large and
extended systems. A lot depends on the values of the T2
amplitudes. If they are small, then the denominator
expressions which define the renormalized [T] and (T)
corrections are very close to 1, in which case there is no
major problem. If the T2 amplitudes are large, then the
denominators may become significantly greater than 1. This
behavior of the R-CCSD[T], R-CCSD(T), CR-CCSD[T], and
CR-CCSD(T) denominator expressions is useful for improving
the results for bond breaking, since the denominators
defining the renormalized [T] and (T) corrections damp the
unphysical values of the standard [T] and (T) corrections
at larger internuclear separations. However, for very
large systems, the denominators defining the renormalized
[T] and (T) corrections may "overdamp" the [T] and (T)
corrections. Numerical tests indicate that lack of strict
size extensivity has very little (fraction of a milli-
hartree) effect on the results of the CR-CCSD[T] and
CR-CCSD(T) calculations for smaller and medium size systems.
For large systems, such as the glycine dimer described by
the 6-31G basis set, the departure from size extensivity,
as measured by forming the difference of the sum of the
energies of isolated glycine molecules from the energy of
the dimer consisting of glycine molecules at very large
(200 bohr) distance, is ca. 3 millihartree (2 kcal/mol).
* * * resource requirements * * *
User can perform LCCD, CCD, and CCSD calculations, that
is without calculating the [T] and (T) corrections, or
calculate the entire set of the standard and renormalized
[T] and (T) corrections, in addition to the CCSD energies
(cf. the values of input variable CCTYP in $CONTRL). The
most expensive steps in CC calculations scale as follows:
1
LCCD, CCD or CCSD No**2 times Nu**4 (iterative)
CCSD[T], CCSD(T),
R-CCSD[T], R-CCSD(T),
CR-CCSD[T], CR-CCSD(T). No**3 times Nu**4 (non-iterative)
The cost of calculating the standard CCSD[T] and CCSD(T)
energies and the cost of calculating the R-CCSD[T] and
R-CCSD(T) energies are essentially the same. The cost of
calculating the triples corrections of the CR-CCSD[T] and
CR-CCSD(T) approaches is essentially twice the cost of
calculating the standard CCSD[T] and CCSD(T) corrections.
Although the triples corections may be seen to grow as the
seventh power of the system size, they often require less
time than the sixth power iterations of the CCSD step,
while providing a great increase in accuracy.
Rough estimates of the memory required are:
CCSD 4 No**2 times Nu**2 + No times Nu**3
CCSD[T], CCSD(T),
R-CCSD[T], R-CCSD(T) 4 No**2 times Nu**2 + No times Nu**3
CR-CCSD[T], CR-CCSD(T) No**2 times Nu**2 + 2 * No times Nu**3
(a faster algorithm), or
4 No**2 times Nu**2 + No times Nu**3
(a somewhat slower but less memory
consuming algorithm)
The program automatically selects the algorithm for the
CR-CCSD[T] and CR-CCSD(T) calculations, depending on the
amount of available memory. The above estimates are rough.
The time required for calculating the CR-CCSD[T] and
CR-CCSD(T) triples corrections is only twice the cost of
calculating the standard CCSD[T] and CCSD(T) corrections.
Thus, by just doubling the CPU time for the noniterative
triples corrections and by selecting CCTYP = CR-CC, we gain
access to all six noniterative triples corrections (the
CCSD[T], CCSD(T), R-CCSD[T], R-CCSD(T), CR-CCSD[T], and
CR-CCSD(T) energies) plus, of course, to the CCSD energy.
At the same time, the CR-CCSD[T] and CR-CCSD(T) results for
stretched nuclear geometries are better than the results of
the standard CCSD[T] and CCSD(T) calculations. In some
cases, choosing CCTYP = R-CC might be reasonable, too.
The choice CCTYP=R-CC gives five different energies (CCSD,
CCSD[T], CCSD(T), R-CCSD[T], and R-CCSD(T)) for the price
of three (CCSD, CCSD[T], and CCSD(T)) as the there is no
extra time needed for the R- theories compared to the
standard ones.
We encourage the user to read the paper
P.Piecuch, S.A.Kucharski, K.Kowalski, M.Musial
Comput.Phys.Comm., 149, 71-96(2002)
where time and memory requirements for various types of CC
calculations are described in considerable detail.
1
* * * restarts * * *
The CC code incorporated in GAMESS is quite good in
converging the CCSD equations with the default guess for
cluster amplitudes. The code is designed to converge in
relatively few iterations for significantly stretched
nuclear geometries, where it is not unusual to obtain large
cluster amplitudes whose absolute values are close to 1.
This is accomplished by combining the standard Jacobi
algorithm with the DIIS extrapolation method of Pulay. The
maximum number of amplitude vectors used in the DIIS
extrapolation procedure is defined by the input variable
MXDIIS. The default for MXDIIS is as follows:
MXDIIS = 5, for 5 < No*Nu,
MXDIIS = 3, for 2 < No*Nu < 6,
MXDIIS = 0, for No*Nu < 3.
Thus, in the vast majority of cases, the default value of
MXDIIS is 5. However, for very small problems, when the
DIIS expansion subspace leads to singular systems of linear
equations, it is necessary to reduce the value of MXDIIS to
2-4 (we chose 3) or switch off DIIS altogether (which is
the case when MXDIIS = 0).
It may, of course, happen that the solver for the CCSD
equations does not converge, in spite of increasing the
maximum number of iterations (input variable MAXCC; the
default value is 30) and in spite of changing the default
value of MXDIIS. In order to facilitate the calculations
in all such cases, we included the restart option in the CC
codes incorporated in GAMESS. Thus, the user can restart a
CCSD (or (L)CCD) calculation from the restart file created
by an earlier CC calculation. In order to use the restart
option, the user must save the disk file CCREST from the
previous CC run (cf. the GAMESS script rungms). A restart
is invoked by entering a nonzero value for IREST, which
should be the number of the last iteration completed, and
must be some value greater than or equal 3. Examples of
using the restart option include the following situations:
o The CCSD program did not converge in MAXCC iterations,
but there is a chance to converge it if the value of
MAXCC is increased. User restarts the calculation with
the increased value of MAXCC.
o User ran a CCSD calculation, obtaining the converged CCSD
energy, but later decided to run a CR-CCSD(T) calculation.
Instead of running the entire CCSD --> CR-CCSD(T) task
again, user restarts the calculation after changing the
value of input variable CCTYP, and entering IREST to
reuse the previous CCSD amplitudes, proceeding at once
to the non-iterative triples.
1
o The CCSD program diverged for some geometry with a
significantly stretched bond. User performs an extra
calculation for a different nuclear geometry, for which
it is easier to converge the CCSD equations, and restarts
the calculation from the restart file generated by an
extra calculation. This technique of restarting the CC
calculations from the cluster amplitudes obtained for a
neighboring nuclear geometry is particularly useful for
scanning PESs and for calculating energy derivatives by
numerical differentiation.
* * * references * * *
The rest of this section is references to the original
formulation of various areas in Coupled-Cluster Theory.
Electronic structure:
J. Cizek, J. Chem. Phys. 45, 4256 (1966)
J. Cizek, Adv. Chem. Phys. 14, 35 (1969),
J. Cizek and J. Paldus, Int. J. Quantum Chem. 5, 359 (1971),
Nuclear theory:
F. Coester, Nucl. Phys. 7, 421 (1958),
F. Coester and H. Kuemmel, Nucl. Phys. 17, 477 (1960).
Coupled-Cluster Method with Doubles (CCD) -
J. Cizek, J. Chem. Phys. 45, 4256 (1966),
J. Cizek, Adv. Chem. Phys. 14, 35 (1969),
J. Cizek and J. Paldus, Int. J. Quantum Chem. 5, 359 (1971),
J.A. Pople, R. Krishnan, H.B. Schlegel, and J.S. Binkley,
Int. J. Quantum Chem. Symp. 14, 545 (1978),
R.J. Bartlett and G.D. Purvis,
Int. J. Quantum Chem. Symp. 14, 561 (1978).
J. Paldus, J. Chem. Phys. 67, 303 (1977)
[orthogonally spin-adapted formulation].
Linearized Coupled-Cluster Method with Doubles (LCCD;
cf., also, D-MBPT(infinity), CEPA(0))
J. Cizek, J. Chem. Phys. 45, 4256 (1966)
J. Cizek, Adv. Chem. Phys. 14, 35 (1969)
R.J. Bartlett, I. Shavitt, Chem.Phys.Lett.50, 190 (1977)
57, 157 (1978) [Erratum]
R. Ahlrichs, Comp. Phys. Commun. 17, 31 (1979).
Coupled-Cluster Method with Singles and Doubles (CCSD) -
G.D.Purvis III, R.J.Bartlett, J.Chem.Phys. 76, 1910 (1982)
[spin-orbital formulation],
P.Piecuch and J.Paldus, Int.J.Quantum Chem. 36, 429 (1989)
[orthogonally spin-adapted formulation],
G.E.Scuseria, A.C.Scheiner, T.J.Lee, J.E.Rice,
H.F.Schaefer III, J. Chem. Phys. 86, 2881 (1987)
[non-orthogonally spin-adapted formulation],
G.E. Scuseria, C.L. Janssen, H.F.Schaefer III
J. Chem. Phys. 89, 7382 (1988)
[non-orthogonally spin-adapted formulation]
T.J. Lee and J.E. Rice, Chem. Phys. Lett. 150, 406 (1988)
[non-orthogonally spin-adapted formulation].
1
Coupled-Cluster Method with Singles and Doubles and
Noniterative Triples, CCSD[T] = CCSD+T(CCSD) -
M. Urban, J. Noga, S. J. Cole, and R. J. Bartlett,
J. Chem. Phys. 83, 4041 (1985),
P. Piecuch and J. Paldus, Theor. Chim. Acta 78, 65 (1990)
[orthogonally spin-adapted formulation].
Coupled-Cluster Method with Singles and Doubles and
Noniterative Triples, CCSD(T) -
K. Raghavachari, G. W. Trucks, J. A. Pople, M. Head-Gordon,
Chem. Phys. Lett. 157, 479 (1989).
Renormalized and Completely Renormalized Coupled-Cluster
Methods, Method of Moments of Coupled-Cluster Equations
(Ground States) -
P.Piecuch, K.Kowalski,
in J.Leszczynski (Ed.), Computational Chemistry:
Reviews of Current Trends, Vol. 5, World Scientific,
Singapore, 2000, p.1
K.Kowalski, P.Piecuch, J.Chem.Phys. 113, 18-35(2000)
K.Kowalski, P.Piecuch, J.Chem.Phys. 113, 5644-5652(2000)
K.Kowalski, P.Piecuch, Chem.Phys.Lett. 344, 165 (2001)
P.Piecuch, S.A.Kucharski, K.Kowalski,
Chem.Phys.Lett. 344, 176(2001)
P.Piecuch, S.A.Kucharski, V.Spirko, K.Kowalski,
J.Chem.Phys. 115, 5796 (2001)
P.Piecuch, K.Kowalski, I.S.O.Pimienta, S.A.Kucharski,
in M.R.Hoffmann, K.G.Dyall (Eds.), Low-Lying Potential
Energy Surfaces, ACS Symposium Series, Vol. 828, Am.
Chem. Society, Washington, D.C., 2002, p. XXX,
M.J.McGuire, K.Kowalski, P.Piecuch,
J.Chem.Phys., in press (2002)
P.Piecuch, S.A.Kucharski, K.Kowalski, M.Musial,
Comput.Phys.Comm., 149, 71-96(2002).
T1 diagnostic:
T.J.Lee, P.R.Taylor Int.J.Quantum Chem., S23, 199-207(1989)
T1>0.02 indicates CCSD may not correct a system which is
not very single reference in nature. (T) corrections
tolerate somewhat greater singles amplitudes.
1
Density Functional Theory
------- ---------- ------
There are actually two DFT programs in GAMESS, one
using the typical grid quadrature for integration of
functionals, and one using resolution of the identity to
avoid the need for grids. The default program is
METHOD=GRID is discussed below, but we will begin with a
description of METHOD=GRIDFREE. The final section is
references to various functionals, and other topics of
interest.
* * * DFTTYP keywords * * *
For convenience in comparing to other DFT programs, the
following table matches DFTTYP to their input keywords:
GAMESS NWChem MOLPRO2000 Gaussian94
Slater slater S HFS
Gill gill96
PBE xpbe96
SVWN slater vwn_5 S,VWN SVWN5
SLYP slater lyp S,LYP SLYP
Becke becke88 B HFB
BVWN becke88 vwn_5 B,VWN BVWN5
BLYP becke88 lyp B,LYP BLYP
B3LYP -see note- B3LYP
This table applies to both grid and grid-free DFTTYPs.
Note that B3LYP in GAMESS is based in part on the VWN5
correlation functional. Since there are five formulae
in the VWN paper for local correlation, other programs
may use other formulae, and therefore generate different
B3LYP energies. For example, NWChem's manual says it uses
the "VWN 1 functional with RPA parameters as opposed to the
prescribed Monte Carlo parameters", but NWChem can be made
to use the VWN5 formula by
xc HFexch 0.20 slater 0.80 becke88 nonlocal 0.72
lyp 0.81 vwn_5 0.19
If you use some other program, its B3LYP energy will be
different from GAMESS if it does not employ VWN5.
* * * Grid-free DFT * * *
The grid-free code is a research tool into the use of
the resolution of the identity to simplify evaluation of
integrals over functionals, rather than quadrature grids.
This trades the use of finite grids and their associated
errors for the use of a finite basis set used to expand
the identity, with an associated truncation error. The
present choice of auxiliary basis sets was obtained by
tests on small 2nd row molecules like NH3 and N2, and
hence the built in bases for the 3rd row are not as well
developed. Auxiliary bases for the remaining elements do
not exist at the present time.
1
The grid-free Becke/6-31G(d) energy at a C1 AM1 geometry
for ethanol is -154.084592, while the result from a run
using the "army grade grid" is -154.105052. So, the error
using the AUX3 RI basis is about 5 milliHartree per 2nd row
atom (the H's must account for some of the error too). The
energy values are probably OK, the differences noted should
by and large cancel when comparing two different geometries.
The grid-free gradient code contains some numerical
inaccuracies, possibly due to the manner in which the RI
is implemented for the gradient. Computed gradients
consequently may not be very reliable. For example, a
Becke/6-31G(d) geometry optimization of water started from
the EXAM08 geometry behaves as:
FINAL E= -76.0439853638, RMS GRADIENT = .0200293
FINAL E= -76.0413274662, RMS GRADIENT = .0231574
FINAL E= -76.0455283912, RMS GRADIENT = .0045887
FINAL E= -76.0457360477, RMS GRADIENT = .0009356
FINAL E= -76.0457239113, RMS GRADIENT = .0001222
FINAL E= -76.0457216186, RMS GRADIENT = .0000577
FINAL E= -76.0457202264, RMS GRADIENT = .0000018
FINAL E= -76.0457202253, RMS GRADIENT = .0000001
Examination shows that the point on the PES where the
gradient is zero is not where the energy is lowest, in
fact the 4th geometry is the lowest encountered.
The behavior for Becke/6-31G(d) ethanol is as follows:
FINAL E= -154.0845920132, RMS GRADIENT = .0135540
FINAL E= -154.0933138447, RMS GRADIENT = .0052778
FINAL E= -154.0885472996, RMS GRADIENT = .0009306
FINAL E= -154.0886268185, RMS GRADIENT = .0002043
FINAL E= -154.0886352947, RMS GRADIENT = .0000795
FINAL E= -154.0885599794, RMS GRADIENT = .0000342
FINAL E= -154.0885514829, RMS GRADIENT = .0000679
FINAL E= -154.0884955093, RMS GRADIENT = .0000205
FINAL E= -154.0886438244, RMS GRADIENT = .0000330
FINAL E= -154.0886596883, RMS GRADIENT = .0000325
FINAL E= -154.0886094081, RMS GRADIENT = .0000120
FINAL E= -154.0886054003, RMS GRADIENT = .0000109
FINAL E= -154.0885939751, RMS GRADIENT = .0000152
FINAL E= -154.0886711482, RMS GRADIENT = .0000439
FINAL E= -154.0886972557, RMS GRADIENT = .0000230
with similar fluctuations through a total of 50 steps
without locating a zero gradient. Note that the second
energy above is substantially below all later points, so
geometry optimizations with the grid-free DFT gradient
code are at this time unsatisfactory.
* * * DFT with grids * * *
The default METHOD=GRID produces good energy and
gradient quantities. For example, when running the same
Becke/6-31G(d) ethanol test case, the default grid will
produce an ethanol energy -154.104863 to be compared to
-154.105052 using the "army grade" grid. Thus energy
errors should be less than a tenth of a milliHartree per
1
atom with the standard grid. Note that the energies are
nonetheless a function of the grid size, just as they are
a function of the basis used, so you must only compare
runs that use the same grid size. The grid code will
give gradient vectors that are accurate, and lead to
satisfactory geometry optimizations. This means that
DFT frequencies obtained by numerical differentiation
should also be OK. RUNTYP=ENERGY, GRADIENT, HESSIAN,
and their chemical combinations for OPTIMIZE, SADPOINT,
IRC, DRC, VSCF, RAMAN, and FFIELD should all work, but
TRANSITN, MOROKUMA, and TDHF should not be used with DFT.
The grid DFT uses symmetry during the quadrature in two
ways. First, the integration runs only over grid points
placed around the symmetry unique atoms. Also, "octant
symmetry" is implemented using an appropriate Abelian
subgroup of the full group. Your run should be done in the
full non-Abelian group, so that the usual integrals and the
SCF steps can exploit full symmetry. The grid evaluation
automatically uses an appropriate subgroup to reduce the
number of grid points for atoms that lie on symmetry axes
or planes. For example, in Cs, atoms lying in the xy plane
will be integrated only over the upper hemisphere of their
grid points. Octant symmetry is not used for any of these:
a) if a non-standard axis orientation is input in $DATA
b) if the angular grid size (NTHE,NTHE0,NPHI,NPI0) is not
a multiple of the octant symmetry factors, such as
NTHE=15 in C2v. The permissible values depend on the
group, but NTHE a multiple of 2 and NPHI a multiple of
4 is generally safe.
* * * references * * *
An excellent overview can be found in Chapter 6 of Frank
Jensen's book. Two other monographs are
"Density Functional Theory of Atoms and Molecules"
R.G.Parr, W.Yang Oxford Scientific, 1989
"A Chemist's Guide to Density Functional Theory"
W.Koch, M.C.Holthausen Wiley-VCH 2001
A delightful thought provoking paper on the relationship of
DFT to conventional quantum mechanics using wavefunctions:
P.M.W.Gill Aust.J.Chem. 54, 661-662(2001)
Some reading on the grid-free approach to density
functional theory is:
Y.C.Zheng, J.Almlof
Chem.Phys.Lett. 214, 397-401(1996)
Y.C.Zheng, J.Almlof
J.Mol.Struct.(Theochem) 288, 277(1996)
K.Glaesemann, M.S.Gordon
J.Chem.Phys. 108, 9959-9969(1998)
K.Glaesemann, M.S.Gordon
J.Chem.Phys. 110, 6580-6582(1999)
K.Glaesemann, M.S.Gordon
J.Chem.Phys. 112, 10738-10745(2000)
1
Three papers on gridding:
C.W.Murray, N.C.Handy, G.L.Laming
Mol.Phys. 78, 997-1014(1993)
P.M.W.Gill, B.G.Johnson, J.A.Pople
Chem.Phys.Lett. 209, 506-512(1993)
R.Lindh, P.-A.Malmqvist, L.Gagliardi
Theoret.Chem.Acc. 106, 178-187(2001)
Slater exchange:
J.C.Slater Phys.Rev. 81, 385-390(1951)
XALPHA is Slater with alpha=0.70
BECKE (often called B88) exchange:
A.D.Becke Phys.Rev. A38, 3098-3100(1988)
GILL (often called G96) exchange:
P.M.W.Gill Mol.Phys. 89, 433-445(1996)
PBE exchange/correlation functional:
J.P.Perdew, K.Burke, M.Ernzerhof
Phys.Rev.Lett. 77, 3865-8(1996); Err. 78,1396(1997)
M.Ernzerhof, G.E.Scuseria
J.Chem.Phys. 110, 5029-5036(1999)
Note that only the exchange functional is in GAMESS.
PBE has not been tested to the same level as the others.
Depristo/Kress exchange:
A.E.DePristo, J.E.Kress J.Chem.Phys. 86, 1425-1428(1987)
VWN (meaning specifically VWN5) correlation:
S.H.Vosko, L.Wilk, M.Nusair
Can.J.Phys. 58, 1200-1211(1980)
This paper has five formulae in it, and since the 5th is
a good quality fit, it states "since formula 5 is easiest
to implement in LSDA calculations, we recommend its use".
PWLOC:
J.D.Perdew, Y.Wang Phys.Rev. B45, 13244-13249(1992)
LYP correlation:
C.Lee, W.Yang, R.G.Parr Phys.Rev. B37, 785-789(1988)
For practical purposes this is used in a transformed way,
involving the square of the density gradient:
B.Miehlich, A.Savin, H.Stoll, H.Preuss
Chem.Phys.Lett. 157, 200-206(1989)
PW91 correlation: (not presently implemented in GAMESS)
J.P.Perdew, J.A.Chevray, S.H.Vosko, K.A.Jackson,
M.R.Pederson, D.J.Singh, C.Fiolhais
Phys.Rev. B46, 6671-6687(1992)
OP (One-parameter Progressive) correlation:
T.Tsuneda, K.Hirao Chem.Phys.Lett. 268, 510-520(1997)
T.Tsuneda, T.Suzumura, K.Hirao
J.Chem.Phys. 110, 10664-10678(1999)
various WIGNER exchange/correlation functionals:
Q.Zhao, R.G.Parr Phys.Rev. A46, 5320-5323(1992)
1
CAMA/CAMB exchange/correlation functionals:
G.J.Laming, V.Termath, N.C.Handy
J.Chem.Phys. 99. 8765-8773(1993)
B3LYP hybrid:
A.D.Becke J.Chem.Phys. 98, 5648-5642(1993)
P.J.Stephens, F.J.Devlin, C.F.Chablowski, M.J.Frisch
J.Phys.Chem. 98, 11623-11627(1994)
R.H.Hertwig, W.Koch Chem.Phys.Lett. 268, 345-351(1997)
The first paper is actually on B3PW91, never mentioning the
B3LYP hybrid, and optimizes the mixing of five functionals
with PW91 as the correlation GGA. The second paper then
proposed use of LYP in place of PW91, without reoptimizing
the mixing ratios of the hybrid. The final paper discusses
the controversy surrounding which VWN functional is used
in the hybrid. GAMESS uses VWN5 in its B3LYP hybrid.
HALF exchange:
This is programmed as 50% HF plus 50% B88 exchange.
BHHLYP exchange/correlation:
This is 50% HF plus 50% B88, with LYP correlation.
Note: neither is the HALF-AND-HALF exchange/correlation:
A.D.Becke J.Chem.Phys. 98, 1372-1377(1993)
which he defined as 50% HF + 50% SVWN.
See http://www.dl.ac.uk/DFTlib/contents.html for other
information about functionals.
The paper of Johnson, Gill, and Pople listed below has a
useful summary of formulae, and details about a gradient
implementation. A paper on 1st and 2nd derivatives of DFT
with respect to nuclear coordinates and applied fields is
A.Komornicki, G.Fitzgerald
J.Chem.Phys. 98, 1398-1421(1993)
A few of the many papers assessing the accuracy of DFT:
B.Miehlich, A.Savin, H.Stoll, H.Preuss
Chem.Phys.Lett. 157, 200-206(1989)
B.G.Johnson, P.M.W.Gill, J.A.Pople
J.Chem.Phys. 98, 5612-5626(1993)
N.Oliphant, R.J.Bartlett
J.Chem.Phys. 100, 6550-6561(1994)
L.A.Curtiss, K.Raghavachari, P.C.Redfern, J.A.Pople
J.Chem.Phys. 106, 1063-1079(1997)
E.R.Davidson Int.J.Quantum Chem. 69, 241-245(1998)
B.J.Lynch, D.G.Truhlar
J.Phys.Chem.A 105, 2936-2941(2001)
R.A.Pascal J.Phys.Chem.A 105, 9040-9048(2001)
Of course there are assessments in many of the functional
papers as well!
On the accuracy of DFT for large molecule thermochemistry:
L.A.Curtiss, K.Ragavachari, P.C.Redfern, J.A.Pople
J.Chem.Phys. 112, 7374-7383(2000)
P.C.Redfern, P.Zapol, L.A.Curtiss, K.Ragavachari
J.Phys.Chem.A 104, 5850-5854(2000)
1
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.
1
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.
1
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 somewhat dated summary of various 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 Odense University, 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, 305-310(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 Hessian * * *
Although quasi-Newton methods require only an
approximation to the true hessian, the choice of this
matrix has a great affect on convergence of the geometry
search.
There is a procedure contained within GAMESS for
guessing a diagonal, positive definite hessian matrix,
HESS=GUESS. If you are using Cartesian coordinates, the
guess hessian is 1/3 times the unit matrix. The guess is
more sophisticated when internal coordinates are defined,
as empirical rules will be used to estimate stretching
and bending force constants. Other force constants are set
to 1/4. The diagonal 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.
1
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 eigen-
values 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.
1
* * * 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.
1
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.
1
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
1
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."
1
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.
1
* * * 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.
1
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.
1
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.
1
* * * 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.
1
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).
1
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!
1
IRC 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 AMPC4 or
RK4), 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), and in the initiation
of the IRC run. The latter might require a hessian matrix
obtained by double differencing, certainly the hessian
should be 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.
1
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 the most accurate integration method
implemented in GAMESS, and is also the most time consuming.
The Gonzalez-Schlegel 2nd order method finds the next
IRC point by a constrained optimization on the surface of
a hypersphere, centered at 1/2 STRIDE along the gradient
vector leading from the previous IRC point. By construction,
the reaction path between two successive IRC points is
thus 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. Be sure
to provide the updated hessian from the previous run when
restarting PACE=GS2.
1
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 corection 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 Journal, 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.
1
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.Gonzales, 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.
1
Gradient Extremals
-------- ---------
This section of the manual, as well as the source code
to trace gradient extremals was written by Frank Jensen of
Odense University.
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).
1
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
1
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.
1
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.
1
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.
1
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. This section describes the three continuum
models which are implemented in the standard version of
GAMESS, and a fourth 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.
* * *
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)
1
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).
* * *
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 also includes a more
exact treatment of the electrostatic interaction with the
surrounding medium, as the electrostatic potential of the
solute generates an 'apparent surface charge' on the
cavity's surface. The computational procedure divides
this surface into small tesserae, on which the charge (and
contributions to the gradient) are evaluated. Typically
the spheres defining the cavity are taken to be 1.2 times
the van der Waals radii. A technical difficulty caused by
the penetration of the solute 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 dispersion and repulsion
contributions to the solvation free energy.
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.
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.
1
Due to its sophistication, users of the PCM model are
strongly encouraged to read the primary literature. The
first references use the boundary element method for solving
the apparent surface charge problem. Recently the integral
equation formalism has been developed as a more numerically
suitable method.
General papers on the PCM method:
1) S.Miertus, E.Scrocco, J.Tomasi
Chem.Phys. 55, 117-129(1981)
2) J.Tomasi, M.Persico Chem.Rev. 94, 2027-2094(1994)
3) R.Cammi, J.Tomasi J.Comput.Chem. 16, 1449-1458(1995)
The GEPOL method for cavity construction:
4) J.L.Pascual-Ahuir, E.Silla, J.Tomasi, R.Bonaccorsi
J.Comput.Chem. 8, 778-787(1987)
Charge renormalization (see also ref. 3):
5) B.Mennucci, J.Tomasi J.Chem.Phys. 106, 5151-5158(1997)
Derivatives with respect to nuclear coordinates:
(energy gradient and hessian)
6) R.Cammi, J.Tomasi J.Chem.Phys. 100, 7495-7502(1994)
7) R.Cammi, J.Tomasi J.Chem.Phys. 101, 3888-3897(1995)
8) M.Cossi, B.Mennucci, R.Cammi
J.Comput.Chem. 17, 57-73(1996)
Derivatives with respect to applied electric fields:
(polarizabilities and hyperpolarizabilities)
9) R.Cammi, J.Tomasi
Int.J.Quantum Chem. Symp. 29, 465-474(1995)
10) R.Cammi, M.Cossi, J.Tomasi
J.Chem.Phys. 104, 4611-4620(1996)
11) R.Cammi, M.Cossi, B.Mennucci, J.Tomasi
J.Chem.Phys. 105, 10556-10564(1996)
12) B. Mennucci, C. Amovilli, J. Tomasi
J.Chem.Phys. submitted.
Cavitation energy:
13) R.A.Pierotti Chem.Rev. 76, 717-726(1976)
14) J.Langlet, P.Claverie, J.Caillet, A.Pullman
J.Phys.Chem. 92, 1617-1631(1988)
Dispersion and repulsion energies:
15) F.Floris, J.Tomasi J.Comput.Chem. 10, 616-627(1989)
16) C.Amovilli, B.Mennucci
J.Phys.Chem.B 101, 1051-1057(1997)
1
Integral Equation Formalism papers. The first of these
deals with anisotropies, the last 2 with nuclear gradients.
17) E.Cances, B.Mennucci, J.Tomasi
J.Chem.Phys. 107, 3032-3041(1997)
18) B.Mennucci, E.Cances, J.Tomasi
J.Phys.Chem.B 101, 10506-17(1997)
19) B.Mennucci, R.Cammi, J.Tomasi
J.Chem.Phys. 109, 2798-2807(1998)
20) J.Tomasi, B.Mennucci, E.Cances
J.Mol.Struct.(THEOCHEM) 464, 211-226(1999)
21) E.Cances, B.Mennucci J.Chem.Phys. 109, 249-259(1998)
22) E.Cances, B.Mennucci, J.Tomasi
J.Chem.Phys. 109, 260-266(1998)
At the present time, the PCM model in GAMESS has the
following limitations:
a) SCFTYP=RHF and MCSCF, only.
b) point group symmetry is switched off internally
during PCM.
c) The PCM model runs in parallel only for IEF=3 and
SCFTYP=RHF, but not for TDHF jobs.
d) electric field integrals at normals to the surface
elements are stored on disk, even during DIRSCF
runs. The file size may be considerable.
e) To minimize common block storage, the maximum
number of spheres forming the cavity is 100, with
an upper limit on the number of surface tesserae
set to 2500. These may be raised by the 'mung'
script listed in the Programming chapter.
f) nuclear derivatives are limited to gradients,
although theory for hessians is given in Ref. 7.
The calculation shown on the next page illustrates the
use of most 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.
1
! 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. B101, 1051-7(1997)
!
! The gas phase FINAL energy is -40.2075980280
! The FINAL energy in PCM water= -40.2143590161
! (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 memory=300000 $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
! reference cited used value for H2O's solvent radius
! which differs from the built in constants.
$PCM IEF=0 ICOMP=2 IREP=1 IDP=1 ICAV=1
SOLVNT=WATER RSOLV=1.35 $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
1
COSMO (conductor-like screening model) represents a
different approach for carrying out polarized continuum
calculations. The model 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 modelled as a conductor
simplifies the electrostatic computations and corrections
may be made a posteriori for dielectric behavior.
The current implementation of COSMO involves the
computation of distributed multipoles up to hexadecapoles
to represent the charge distribution of the molecule
within the cavity. The multipole moments induce the
formation of charges on the surface of the cavity that
contains the molecule. These charges are then fed back
into the SCF, and both the molecular wavefunction and
the surface charges are iterated to self-consistency.
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.
At present, 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, which
will be formulated for GAMESS by Baldridge et al.
The simplicity of the COSMO model allows computation of
gradients, allowing optimization within the context of the
solvent. The method is programmed for closed shell RHF
energy and gradient, and the MP2 energy correction may be
obtained.
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)
1
A final possible continuum treatment is the "solution
model 5" approach. Ab initio SM5 is described in
J.Li, G.D.Hawkins, C.J.Cramer, D.G.Truhlar
Chem.Phys.Lett. 288, 293-298(1998)
SM5 represents the molecule's electrostatics as a set of
atomic point charges. These are chosen by a procedure
based on correcting Lowdin atomic charges according to a
quadratic function of the computed Mayer bond orders,
which is designed to better reproduce experimental dipole
moments. These charges are called "charge model 2", and
CM2 is described in
J.Li, T.Zhu, C.J.Cramer, D.G.Truhlar
J.Phys.Chem.A 102, 1820-1831(1998)
In addition to a self-consistent reaction field treatment
of the CM2 electrostatics, SM5 includes a term accounting
for the following first solvation shell effects: cavity
creation, dispersion, and changes in solvent structure,
which are modeled by atomic surface tension parameters.
It is possible to use this code simply to extract gas phase
CM2 charges. The implementation is termed GAMESOL (one S),
by
J.Li, G.D.Hawkins, D.A.Liotard, C.J.Cramer, D.G.Truhlar
and is available at
http://comp.chem.umn.edu/gamesol
After signing a license not much more stringent than the
license for GAMESS itself, you can obtain the new source
code from U. Minnesota. The interface is not clean, as
considerable code is inserted directly into RHFUHF and
other GAMESS modules, so you must be very careful to obtain
code that matches the dates on the top of your original
GAMESS source files.
1
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. 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 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 RHF/DZP water
molecules to study aqueous solvation effects, see references
1-3. Co-workers at NIST have also used EFP1 to model parts
of enzymes, see reference 4. EFP2 is a more general method
that is applicable to any species, including water, and its
repulsive potential is obtained from first principles. EFP2
forms the basis of the covalent EFP method described below.
*** 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. It is possible to input a screening term
that accounts for the charge penetration. This screening
term is automatically included for EFP1. Typically the
multipole expansion points are located on atomic nuclei
and at bond midpoints.
1
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.
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 atomic nuclei within the fragment.
The repulsive potential for EFP2 was derived based on an
overlap expansion using localized molecular orbitals, as
described in references 4-6. The EFP2 repulsive potential
has no fitted parameters, and it can be automatically
generated during a RUNTYP=MAKEFP job, as described below.
*** Constructing an EFP1 using GAMESS ***
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 molecular 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,
or the charge penetration screening parameters.
1
Since for EFP1 a user cannot generate all terms necessary
to define a new $FRAGNAME group using GAMESS, in practice
the usage of EFP1 is limited to using the internally stored
H2OEF2 potential mentioned below.
*** Constructing an EFP2 using GAMESS ***
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 moleulce can be generated as follows:
(a) Choose a basis set and geometry for the molecule of
interest. The geometry is ordinarily optimized at the
Hartree-Fock level of theory with the 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 .dat file. The only user-defined variable
that must be filled in is the FRAGNAME.
(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. The coordinates of any other
atoms are automatically fixed by the program.
*** 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, and IRC (see reference 3), 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 fragments are present.
2. The ab initio system must be treated with RHF, ROHF,
UHF, the open shell SCF wavefunctions permitted by the GVB
code, or MCSCF. The correlated methods such as MP2 and CI
should not be used, since the available H2OEF2 potential
1
was derived at the RHF level, and therefore does not
contain dispersion terms. A correlated computation on
the ab initio system without these terms in the EFP will
probably lead to unphysical results.
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.
Enter NOSYM=1 in $CONTRL to accomplish this.
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.
8. For EFP2, the charge penetration (screening) term is
not currently implemented for automatic generation by
RUNTYP=MAKEFP, so generation of an EFP2 as described above
will not include screening. This is considered to be a
minor omission that will be corrected in a future release.
*** Global optimization ***
If there are a large number of effective fragments, 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 the Metropolis
algorithm to move just one of the effective fragments. If
desired, the method of Parks to move all fragments at once
may be tried, by changing ALPHA from zero and setting
BOLTWT=AVESTEP instead of STANDARD.
1
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.
If there is an ab initio portion present in your
system, it is probably impractical to carry out a simulated
annealing protocol. However, a single temperature Monte
Carlo calculation may be feasible. In particular, you may
wish to avoid the local minimization steps, and instead
manually examine the structures from the Monte Carlo steps
in order to choose a few for full geometry optimization.
Note that SMODIF input can allow the ab initio part of the
system to participate in the Monte Carlo jumps. However,
this should be done with caution.
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 the paper by Paul Day,
Ruth Pachter, Mark Gordon, and Grant Merrill listed in the
EFP references at the end of this section. It resembles
the work of
D.J.Wales, M.P.Hodges Chem.Phys.Lett. 286, 65-72 (1998).
1
*** Practical hints for using EFPs ***
At the present time, we have only one EFP suitable for
general use. This EFP models water, and its numerical
parameters are internally stored, using the fragment name
H2OEF2. These numerical parameters are improved values
over the H2OEF1 set which were presented and used in
reference 2, and they also include the improved EFP-EFP
repulsive term defined in reference 3. The H2OEF2 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.
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 H2OEF2 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.
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.
To do this, simply add a $PCM group to your input, in
addition to the $EFRAG. The simultaneous use of EFP and
PCM is presently limited to energy calculations, so any
geometry optimization must be done with only $EFRAG input.
1
*** 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.
1
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, followed by 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.
1
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.
*** References ***
The first of these is more descriptive, and the second
has a very detailed derivation of the method. The 13th
paper in the list is an overview article.
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. "A study of aqueous glutamic acid using the effective
fragment potential model"
P.N.Day, R.Pachter J.Chem.Phys. 107, 2990-9(1997)
6. "Solvation and the excited states of formamide"
M.Krauss, S.P.Webb J.Chem.Phys. 107, 5771-5(1997)
1
7. "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)
8. "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)
9. "Solvation of Sodium Chloride: EFP study of NaCl(H2O)n"
C.P.Petersen, M.S.Gordon J.Phys.Chem.A 103, 4162-6(1999)
10. "QM/MM boundaries across covalent bonds: a frozen LMO based
approach for the Effective Fragment Potential method"
V.Kairys, J.H.Jensen J.Phys.Chem.A 104, 6656-65(2000)
11. "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)
12. "A combined discrete/continuum solvation model:
Application to glycine" P.Bandyopadhyay, M.S.Gordon
J.Chem.Phys. 113, 1104-9(2000)
13. "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)
14. "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)
15. "Active site structure and mechanism of Human Glyoxalase I"
U.RIchter, M.Krauss J.Am.Chem.Soc. 123, 6973-6982(2001)
16. "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)
17. "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)
18. "An integrated effective fragment-polarizable continuum
approach to solvation: Theory and application to glycine"
P.Banyopadhyay, M.S.Gordon, B.Mennucci, J.Tomasi
J.Chem.Phys. 116, 5023-5032(2002)
1
MOPAC calculations within GAMESS
----- ------------ ------ ------
Parts of MOPAC 6.0 have been included in GAMESS so
that the GAMESS user has access to three semiempirical
wavefunctions: MNDO, AM1 and PM3. 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 make it in
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 one of these
calculations. All you need to do is specify GBASIS=MNDO,
AM1, or PM3 in the $BASIS group. Note that this not only
picks a particular Slater orbital basis, it also selects a
particular "hamiltonian", namely a particular parameter
set.
MNDO, AM1, and PM3 will not work with every option in
GAMESS. Currently the semiempirical wavefunctions support
SCFTYP=RHF, UHF, and ROHF in any combination with
RUNTYP=ENERGY, GRADIENT, OPTIMIZE, SADPOINT, HESSIAN, and
IRC. Note that all hessian runs are numerical finite
differencing. The MOPAC CI and half electron methods are
not supported.
Because the majority of the implementation is GAMESS
rather than MOPAC you will notice a few improvments.
Dynamic memory allocation is used, so that 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.
1
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 basis sets are unavailable. This means $ELMOM,
$ELPOT, etc. properties are unavailable. Likewise the
solvation models do not work with semi-empirical runs.
Note that MOPAC6 did not include d STO functions, and it
is therefore quite impossible to run transition metals.
To reduce CPU time, only the EXTRAP convergence
accelerator is used by the SCF procdures. For difficult
cases, the DIIS, RSTRCT, and/or SHIFT options will work,
but may add significantly to the run time. With the
Huckel guess, most calculations will converge acceptably
without these special options.
MOPAC parameters exist for the following elements.
The quote 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 * Be * C N O F
Na' * Al Si P S Cl Na' Mg Al Si P S Cl
K' * ... Zn * Ge * * Br K' * ... Zn Ga Ge As Se Br
Rb' * ... * * Sn * * I Rb' * ... Cd In Sn Sb Te I
* * ... Hg * * * * * ... Hg Tl Pb Bi
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 MNDO, AM1 or PM3 is "how well is my system
modeled with an ab initio minimal basis sets, such as
STO-3G?" If the answer is "not very well" there is a good
chance that a semiempirical description is equally poor.
1
Molecular Properties
--------- ----------
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).
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 electron density is electron/bohr**3
for the total density, and 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 intensities in A**4/amu multiply by 6.0220E-09 for
units of cm**4/g.
1
Localization tips
------------ ----
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.
1
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 make
plots of the orbitals with either the MacMolPlt or PLTORB
graphics codes, or perhaps to read the localized orbitals
in during a RUNTYP=PROP job to examine their Mulliken
populations.
1
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
1
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.
1
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)
1
* * *
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.
1
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
moments are computed using both the "length (dipole) form"
and "velocity form". The two values will be slightly
different as the CI wavefunction does not exactly satisfy
the Hellmann-Feynman theorem. This basic 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. The full Breit-Pauli operator can be
computed exactly (OPERAT=HSO2), or approximated in two
ways: complete elimination of the 2e- term, whose absence
can be approximately accounted for by means of effective
nuclear charges (HSO1), or by inclusion of only the core-
active 2e- terms which typically account for 90% or more
of the two electron term, while saving most of the 2e-
terms' CPU cost (HSO2P).
Spin-orbit runs can be done for general spins, for
more than two different spin multiplicities at once, for
general active spaces. At times, when the spatial wave-
function is degenerate, a spin-orbit run may involve only
one spin multiplicity, e.g. a triplet-pi state in a linear
molecule. The most common case is two different spins.
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. Compared to older versions,
the basis set has been fully generalized to allow any s, p,
d, f, g, or L functions.
1
* * * * 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 H-so elements off diagonal. 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.
1
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 ISTSYM 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.
1
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 which is actually Kh symmetry (full rotation-reflection
group) can be entered in D2h, the 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 ISTSYM 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 details * * * *
Spin-orbit coupling is 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
1
given by the square of the spin-orbit coupling matrix
element divided by the difference of the adibatic 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.
The choice between HSO2 and HSO2FF is very often in
favour of the former. HSO2 computes the matrix elements in
CSF basis and then contracts them with CI coefficients,
whereas HSO2FF uses a generalised 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.
The spin-orbit operator contains a one electron term
arising from Pauli's reduction of the hydrogenic Dirac
equation to one-component form, and a two electron term
added by Breit. 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.
As an approximation, the nuclear charge appearing in
the one electron term can be regarded as an empirical scale
factor, compensating for the omission of the two electron
operator. 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.
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
1
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 HSO2PP 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.
1
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 diagonalisation. 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.
1
* * * * 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:
2. F.Weinhold, J.Chem.Phys. 54,1874-1881(1970)
3. C.W.Bauschlicher, S.R.Langhoff
Theoret.Chim.Acta 79:93-103(1991)
4. "Intramediate 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 Zeff spin-orbit coupling, and ZEFTYP 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)
References for full Breit-Pauli spin-orbit coupling:
11. T.R.Furlani, H.F.King
J.Chem.Phys. 82, 5577-5583 (1985)
12. H.F.King, T.R.Furlani
J.Comput.Chem. 9, 771-778 (1988)
13. D.G.Fedorov, M.S.Gordon
J.Chem.Phys. 112, 5611-5623 (2000)
with the latter including information on the partial
two electron operator method.
Reference for SO-MCQDPT:
14. D.G.Fedorov, J.P.Finley Phys.Rev.A 64, 042502 (2001)
Recent applications:
15. D.G.Fedorov, M.Evans, Y.Song, M.S.Gordon, C.Y.Ng
J.Chem.Phys. 111, 6413-6421 (1999)
16. D.G.Fedorov, M.S.Gordon, Y.Song, C.Y.Ng
J.Chem.Phys. 115, 7393-7400 (2001)
17. B.J.Duke J.Comput.Chem. 22, 1552-1556 (2001)
1
* * *
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.
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.
!
! The two SCF energies below give an excitation energy
! of 16,044 cm-1 to the 2-P term. The computed spin-orbit
! levels are at RELATIVE E=-10.269 and 5.135 cm-1, which
! means the 2-P level interval is 15.404 cm-1.
!
! Charlotte Moore's Atomic Energy Levels, volume 1, gives
! the experimental 2-P interval as 17.1963, the 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=.true. 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
1
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
!
1
! 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
--- $MCQD1 nosym=1 nstate=6 mult=1 INORB=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 INORB=1 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
$MCQD1 nosym=-1 mult=1 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=1 wstate(1)=1,1,1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD2 nosym=-1 mult=1 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=2 wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD3 nosym=-1 mult=1 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=3 wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD4 nosym=-1 mult=1 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=4 wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD5 nosym=-1 mult=3 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=2 wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD6 nosym=-1 mult=3 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=3 wstate(1)=1 thrcon=1e-8 thrgen=1e-10 $END
$MCQD7 nosym=-1 mult=3 INORB=1 iforb=3
nmofzc=1 nmodoc=0 nmoact=8
istsym=4 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
!
1
$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 ...
$END