Continuum Solvation Methods
In a very thorough 1994 review of continuum solvation
models, Tomasi and Persico divide the possible approaches
to the treatment of solvent effects into four categories:
a) virial equations of state, correlation functions
b) Monte Carlo or molecular dynamics simulations
c) continuum treatments
d) molecular treatments
The Effective Fragment Potential method, documented in the
following section of this chapter, falls into the latter
category, as each EFP solvent molecule is modeled as a
distinct object (discrete solvation). This section
describes the four continuum models which are implemented
in the standard version of GAMESS, and a fifth model which
can be interfaced.
Continuum models typically form a cavity of some sort
containing the solute molecule, while the solvent outside
the cavity is thought of as a continuous medium and is
categorized by a limited amount of physical data, such as
the dielectric constant. The electric field of the charged
particles comprising the solute interact with this
background medium, producing a polarization in it, which in
turn feeds back upon the solute's wavefunction.
Self Consistent Reaction Field (SCRF)
A simple continuum model is the Onsager cavity model,
often called the Self-Consistent Reaction Field, or SCRF
model. This represents the charge distribution of the
solute in terms of a multipole expansion. SCRF usually
uses an idealized cavity (spherical or ellipsoidal) to
allow an analytic solution to the interaction energy
between the solute multipole and the multipole which this
induces in the continuum. This method is implemented in
GAMESS in the simplest possible fashion:
i) a spherical cavity is used
ii) the molecular electrostatic potential of the
solute is represented as a dipole only, except
a monopole is also included for an ionic solute.
The input for this implementation of the Kirkwood-Onsager
model is provided in $SCRF.
Some references on the SCRF method are
1. J.G.Kirkwood J.Chem.Phys. 2, 351 (1934)
2. L.Onsager J.Am.Chem.Soc. 58, 1486 (1936)
3. O.Tapia, O.Goscinski Mol.Phys. 29, 1653 (1975)
4. M.M.Karelson, A.R.Katritzky, M.C.Zerner
Int.J.Quantum Chem., Symp. 20, 521-527 (1986)
5. K.V.Mikkelsen, H.Agren, H.J.Aa.Jensen, T.Helgaker
J.Chem.Phys. 89, 3086-3095 (1988)
6. M.W.Wong, M.J.Frisch, K.B.Wiberg
J.Am.Chem.Soc. 113, 4776-4782 (1991)
7. M.Szafran, M.M.Karelson, A.R.Katritzky, J.Koput,
M.C.Zerner J.Comput.Chem. 14, 371-377 (1993)
8. M.Karelson, T.Tamm, M.C.Zerner
J.Phys.Chem. 97, 11901-11907 (1993)
The method is very sensitive to the choice of the solute
RADIUS, but not very sensitive to the particular DIELEC of
polar solvents. The plots in reference 7 illustrate these
points very nicely. The SCRF implementation in GAMESS is
Zerner's Method A, described in the same reference. The
total solute energy includes the Born term, if the solute
is an ion. Another limitation is that a solute's electro-
static potential is not likely to be fit well as a dipole
moment only, for example see Table VI of reference 5 which
illustrates the importance of higher multipoles. Finally,
the restriction to a spherical cavity may not be very
representative of the solute's true shape. However, in the
special case of a roundish molecule, and a large dipole
which is geometry sensitive, the SCRF model may include
sufficient physics to be meaningful:
M.W.Schmidt, T.L.Windus, M.S.Gordon
J.Am.Chem.Soc. 117, 7480-7486(1995).
Most cases should choose PCM (next section) over SCRF!!!
Polarizable Continuum Model (PCM)
A much more sophisticated continuum method, named the
Polarizable Continuum Model, is also available. The PCM
method places a solute in a cavity formed by a union of
spheres centered on each atom. PCM includes a more exact
treatment of the electrostatic interaction of the solute
with the surrounding medium, on the cavity's surface. The
computational procedure divides this surface into many
small tesserae, each having a different "apparent surface
charge", reflecting the solute's and other tesserae's
electric field at each. These surface charges are the PCM
model's "solvation effect" and make contributions to the
energy and to the gradient of the solute.
Typically the cavity is defined as a union of atomic
spheres, which should be roughly 1.2 times the atomic van
der Waals radii. A technical difficulty caused by the
penetration of the solute's charge density outside this
cavity is dealt with by a renormalization. The solvent is
characterized by its dielectric constant, surface tension,
size, density, and so on. Procedures are provided not only
for the computation of the electrostatic interaction of the
solute with the apparent surface charges, but also for the
cavitation energy, and for the dispersion and repulsion
contributions to the solvation free energy.
Methodology for solving the Poisson equation to obtain
the "apparent surface charges" has progressed from D-PCM to
IEF-PCM to C-PCM over time, with the latter preferred.
Iterative solvers require far less computer resources than
direct solvers. Advancements have also been made in
schemes to divide the surface cavity into tiny tesserae.
As of fall 2008, the FIXPVA tessellation, which has smooth
switching functions for tesserae near sphere boundaries,
together with iterative C-PCM, gives very satisfactory
geometry optimizations for molecules of 100 atoms. The
FIXPVA tessellation was extended to work for cavitation
(ICAV), dispersion (IDP), and repulsion (IREP) options in
fall 2009, and dispersion/repulsion (IDISP) in spring 2010.
Other procedures remain, and make the input seem complex,
but their use is discouraged. Thus
$PCM SOLVNT=WATER $END
chooses iterative C-PCM (IEF=-10) and FIXPVA tessellation
(METHOD=4 in $TESCAV) to do basic electrostatics in an
accurate fashion.
The main input group is $PCM, with $PCMCAV providing
auxiliary cavity information. If any of the optional
energy computations are requested in $PCM, the additional
input groups $IEFPCM, $NEWCAV, $DISBS, or $DISREP may be
required.
It is useful to summarize the various cavities used by
PCM, since as many as three cavities may be used:
the basic cavity for electrostatics,
cavity for cavitation energy, if ICAV=1,
cavity for dispersion/repulsion, if IDISP=1.
The first and second share the same radii (see RADII in
$PCMCAV), which are scaled by ALPHA=1.2 for electrostatics,
but are used unscaled for cavitation. The dispersion
cavity is defined in $DISREP with a separate set of atomic
radii, and even solvent molecule radii! Only the
electrostatics cavity can use any of the GEPOL-GB, GEPOL-AS
or recommended FIXPVA tessellation, while the other two use
only the original GEPOL-GB.
Radii are an important part of the PCM parameterization.
Their values can have a significant impact on the quality
of the results. This is particularly true if the solute is
charged, and thus has a large electrostatic interaction
with the continuum. John Emsley's book "The Elements" is a
useful source of van der Waals and other radii.
PCM is at heart a means of treating the electrostatic
interactions between the solute's wavefunction and a
dielectric model for the bulk solvent. The former is
represented as an electron density from whatever quantum
mechanical treatment is used for the solute, and the latter
is a set of surface charges on the finite elements of the
cavity (tessellation). This leaves out other important
contributions to the solvation energy! These include the
energy needed to make a hole in the solvent (cavitation
energy), dispersion or repulsive interactions between the
solute and solvent, and in the SMD model (see below)
solvent structure changes such as would occur in the first
solvation shell. Some empirical formulae for such "CDS"
corrections are provided as keywords ICAV, IDISP, IREP/IDP,
which may not work with all wavefunctions, and may not be
compatible with gradients.
The SMD model gives an alternative set of such "CDS"
corrections, which are compatible with nuclear gradients:
see SMD=.TRUE. in $PCM. A more detailed description of SMD
is given in the paper cited below. The SMD solvent
parameters is described (as of 2010) in
http://comp.chem.umn.edu/solvation/mnsddb.pdf
This gives numerical parameters, all built into GAMESS, as
the various SOLX values, for SOLVNT=
ACETACID CLPROPAN PHOPH EGME E2PENTEN
ACETONE OCLTOLUE DPROAMIN MEACETAT PENTACET
ACETNTRL M-CRESOL DODECAN MEBNZATE PENTAMIN
ACETPHEN O-CRESOL MEG MEBUTATE PFB
ANILINE CYCHEXAN ETSH MEFORMAT BENZALCL
ANISOLE CYCHEXON ETHANOL MIBK PROPANAL
BENZALDH CYCPENTN ETOAC MEPROPYL PROPACID
BENZENE CYCPNTOL ETOME ISOBUTOL PROPANOL
BENZNTRL CYCPNTON EB TERBUTOL PROPNOL2
BENZYLCL DECLNCIS PHENETOL NMEANILN PROPNTRL
BRISOBUT DECLNTRA C6H5F MECYCHEX PROPENOL
BRBENZEN DECLNMIX FOCTANE NMFMIXTR PROPACET
BRETHANE DECANE FORMAMID ISOHEXAN PROPAMIN
BROMFORM DECANOL FORMACID MEPYRID2 PYRIDINE
BROCTANE EDB12 HEPTANE MEPYRID3 C2CL4
BRPENTAN DIBRMETN HEPTANOL MEPYRID4 THF
BRPROPA2 BUTYLETH HEPTNON2 C6H5NO2 SULFOLAN
BRPROPAN ODICLBNZ HEPTNON4 C2H5NO2 TETRALIN
BUTANAL EDC12 HEXADECN CH3NO2 THIOPHEN
BUTACID C12DCE HEXANE NTRPROP1 PHSH
BUTANOL T12DCE HEXNACID NTRPROP2 TOLUENE
BUTANOL2 DCM HEXANOL ONTRTOLU TBP
BUTANONE ETHER HEXANON2 NONANE TCA111
BUTANTRL ET2S HEXENE NONANOL TCA112
BUTILE DIETAMIN HEXYNE NONANONE TCE
NBA MI C6H5I OCTANE ET3N
NBUTBENZ DIPE IOBUTANE OCTANOL TFE222
SBUTBENZ DMDS C2H5I OCTANON2 TMBEN124
TBUTBENZ DMSO IOHEXDEC PENTDECN ISOCTANE
CS2 DMA CH3I PENTANAL UNDECANE
CARBNTET CISDMCHX IOPENTAN NPENTANE M-XYLENE
CLBENZEN DMF IOPROPAN PENTACID O-XYLENE
SECBUTCL DMEPEN24 CUMENE PENTANOL P-XYLENE
CHCL3 DMEPYR24 P-CYMENE PENTNON2 XYLENEMX
CLHEXANE DMEPYR26 MESITYLN PENTNON3
CLPENTAN DIOXANE METHANOL PENTENE
and provides a translation table to full chemical names, if
you can't guess from the input choices given above. The
translations can also be found in the source code. Two
important things to note about SMD are:
a) the atomic radii are changed, so although the
algorithms for electrostatics are those of standard PCM,
the numerical results for the electrostatics do change.
b) SMD's parameterization was developed for IEF-PCM
using GEPOL tessellation with a fine grid: IEF=-3 and
MTHALL=2, NTSALL=240. However it is considered acceptable
to use SMD's parameters, unchanged, with C-PCM and with the
FIXPVA tessellation, at default coarseness. Hence, input
such as
$pcm solvnt=dmso smd=.true. $end
is enough to carry out a SMD-style C-PCM treatment in DMSO.
c) The CDS correction involves cavitation, dispersion,
and as a collective "solvent structure contribution"
estimates for partial hydrogen bonding, repulsion, and
deviation of the dielectric constant from its bulk value.
d) See also SMVLE in the more sophisticated SS(V)PE
continuum model's description.
Solvation of course affects the non-linear optical
properties of molecules. The PCM implementation extends
RUNTYP=TDHF to include solvent effects. Both static and
frequency dependent hyperpolarizabilities can be found.
Besides the standard PCM electrostatic contribution, the
IREP and IDP keywords can be used to determine the effects
of repulsion and dispersion on the polarizabilities.
The implementation of the PCM model in GAMESS has
received considerable attention from Hui Li and Jan Jensen
at the University of Iowa, Iowa State University, and
University of Nebraska. This includes new iterative
techniques to solving the surface charge problem, new
tessellations that provide for numerically stable nuclear
gradients, the implementation of C-PCM equations, the
extension of PCM to all SCFTYPs and TDDFT, development of
an interface with the EFP model (quo vadis), and
heterogenous dielectric. Dmitri Fedorov at AIST has
interfaced PCM to the FMO method (quo vadis), and reduced
storage requirements.
Due to its sophistication, users of the PCM model are
strongly encouraged to read the primary literature:
Of particular relevance to PCM in GAMESS:
1) "Continuum solvation of large molecules described by
QM/MM: a semi-iterative implementation of the PCM/EFP
interface"
H.Li, C.S.Pomelli, J.H.Jensen
Theoret.Chim.Acta 109, 71-84(2003)
2) "Improving the efficiency and convergence of geometry
optimization with the polarizable continuum model: new
energy gradients and molecular surface tessellation"
H.Li, J.H.Jensen J.Comput.Chem. 25, 1449-1462(2004)
3) "The polarizable continuum model interfaced with the
Fragment Molecular Orbital method"
D.G.Fedorov, K.Kitaura, H.Li, J.H.Jensen, M.S.Gordon
J.Comput.Chem. 27, 976-985(2006)
4) "Energy gradients in combined Fragment Molecular Orbital
and Polarizable Continuum Model (FMO/PCM)"
H.Li, D.G.Fedorov, T.Nagata, K.Kitaura, J.H.Jensen,
M.S.Gordon J.Comput.Chem. 31, 778-790(2010)
5) "Continuous and smooth potential energy surface for
conductor-like screening solvation model using fixed points
with variable area"
P.Su, H.Li J.Chem.Phys. 130, 074109/1-13(2009)
6) "Heterogenous conductorlike solvation model"
D.Si, H.Li J.Chem.Phys. 131, 044123/1-8(2009)
7) "Quantum mechanical/molecular mechanical/continuum style
solvation model: linear response theory, variational
treatment, and nuclear gradients"
H.Li J.Chem.Phys. 131, 184103/1-8(2009)
8) "Smooth potential energy surface for cavitation,
dispersion, and repulsion free energies in polarizable
continuum model"
Y.Wang, H.Li J.Chem.Phys. 131, 206101/1-2(2009)
9) "Excited state geometry of photoactive yellow protein
chromophore: a combined conductorlike polarizable continuum
model and time-dependent density functional study"
Y.Wang, H.Li J.Chem.Phys. 133, 034108/1-11(2010)
Paper number 7 is about the treatment of QM systems with
the solvation models EFP and/or C-PCM.
SMD and its CDS cavitation/dispersion/solvent structure
corrections are described in
"Universal solution model based on solute electron density
and on a continuum model of the solvent defined by the bulk
dielectric constant and atomic surface tensions"
A.V.Marenich, C.J.Cramer, D.G.Truhlar
J.Phys.Chem.B 113, 6378-6396(2009)
General papers on PCM:
10) S.Miertus, E.Scrocco, J.Tomasi
Chem.Phys. 55, 117-129(1981)
11) J.Tomasi, M.Persico Chem.Rev. 94, 2027-2094(1994)
12) R.Cammi, J.Tomasi J.Comput.Chem. 16, 1449-1458(1995)
13) J.Tomasi, B.Mennucci, R.Cammi
Chem.Rev. 105, 2999-3093(2005)
The GEPOL-GB method for cavity construction:
14) J.L.Pascual-Ahuir, E.Silla, J.Tomasi, R.Bonaccorsi
J.Comput.Chem. 8, 778-787(1987)
Charge renormalization (see also ref. 12):
15) B.Mennucci, J.Tomasi J.Chem.Phys. 106, 5151-5158(1997)
Derivatives with respect to nuclear coordinates:
(energy gradient and hessian) See also paper 2 and 3.
16) R.Cammi, J.Tomasi J.Chem.Phys. 100, 7495-7502(1994)
17) R.Cammi, J.Tomasi J.Chem.Phys. 101, 3888-3897(1995)
18) M.Cossi, B.Mennucci, R.Cammi
J.Comput.Chem. 17, 57-73(1996)
Derivatives with respect to applied electric fields:
(polarizabilities and hyperpolarizabilities)
19) R.Cammi, J.Tomasi
Int.J.Quantum Chem. Symp. 29, 465-474(1995)
20) R.Cammi, M.Cossi, J.Tomasi
J.Chem.Phys. 104, 4611-4620(1996)
21) R.Cammi, M.Cossi, B.Mennucci, J.Tomasi
J.Chem.Phys. 105, 10556-10564(1996)
22) B. Mennucci, C. Amovilli, J. Tomasi
Chem.Phys.Lett. 286, 221-225(1998)
Cavitation energy:
23) R.A.Pierotti Chem.Rev. 76, 717-726(1976)
24) J.Langlet, P.Claverie, J.Caillet, A.Pullman
J.Phys.Chem. 92, 1617-1631(1988)
Dispersion and repulsion energies:
25) F.Floris, J.Tomasi J.Comput.Chem. 10, 616-627(1989)
26) C.Amovilli, B.Mennucci
J.Phys.Chem.B 101, 1051-1057(1997)
Integral Equation Formalism PCM. The first of these
deals with anisotropies, the last 2 with nuclear gradients.
27) E.Cances, B.Mennucci, J.Tomasi
J.Chem.Phys. 107, 3032-3041(1997)
28) B.Mennucci, E.Cances, J.Tomasi
J.Phys.Chem.B 101, 10506-17(1997)
29) B.Mennucci, R.Cammi, J.Tomasi
J.Chem.Phys. 109, 2798-2807(1998)
30) J.Tomasi, B.Mennucci, E.Cances
J.Mol.Struct.(THEOCHEM) 464, 211-226(1999)
31) E.Cances, B.Mennucci J.Chem.Phys. 109, 249-259(1998)
32) E.Cances, B.Mennucci, J.Tomasi
J.Chem.Phys. 109, 260-266(1998)
Conductor PCM (C-PCM):
33) V.Barone, M.Cossi J.Phys.Chem.A 102, 1995-2001(1998)
34) M.Cossi, N.Rega, G.Scalmani, V.Barone
J.Comput.Chem. 24, 669-681(2003)
C-PCM with TD-DFT:
35) M.Cossi, V.Barone J.Chem.Phys. 115, 4708-4717(2001)
See also paper #8 above for the coding in GAMESS.
At the present time, the PCM model in GAMESS has the
following limitations:
a) Any SCFTYP may be used (RHF to MCSCF). MP2 or DFT
may be used with any of the RHF, UHF, and ROHF
gradient programs. Closed shell TD-DFT excited
state gradients may also be used.
CI and Coupled Cluster programs are not available.
b) semi-empirical methods may not be used.
c) the only other solvent method that may be used at
used with PCM is the EFP model.
d) point group symmetry is switched off internally
during PCM.
e) The PCM model runs in parallel for IEF=3, -3, 10,
or -10 and for all 5 wavefunctions (energy or
gradient), but not for RUNTYP=TDHF jobs.
f) D-PCM stores electric field integrals at normals to
the surface elements on disk.
IEF-PCM and C-PCM using the explicit solver (+3 and
+10) store electric potential integrals at normals
to the surface on disk.
This is true even for direct AO integral runs, and
the file sizes may be considerable (basis set size
squared times the number of tesserae).
IEF-PCM and C-PCM with the iterative solvers do not
store the potential integrals, when IDIRCT=1 in the
$PCMITR group (this is the default)
g) nuclear derivatives are limited to gradients,
although theory for hessians is given in paper 17.
* * *
The only PCM method prior to Oct. 2000 was D-PCM, which
can be recovered by selecting IEF=0 and ICOMP=2 in $PCM.
The default PCM method between Oct. 2000 and May 2004 was
IEF-PCM, recoverable by IEF=-3 (but 3 for non-gradient
runs) and ICOMP=0. As of May 2004, the default PCM method
was changed to C-PCM (IEF=-10, ICOMP=0). The extension of
PCM to all SCFTYPs as of May 2004 involved a correction to
the MCSCF PCM operator, so that it would reproduce RHF
results when run on one determinant, meaning that it is
impossible to reproduce prior MCSCF PCM calculations.
The cavity definition was GEPOL-GB (MTHALL=1 in $TESCAV)
prior to May 2004, GEPOL-AS (MTHALL=2) from then until
September 2008, and FIXPVA (MTHALL=4) to the present time.
The option for generation of 'extra spheres' (RET in $PCM)
was changed from 0.2 to 100.0, to suppress these, in June
2003.
* * *
In general, use of PCM electrostatics is very simple, as
may be seen from exam31.inp supplied with the program.
The calculation shown next illustrates the use of some
of the older PCM options. Since methane is non-polar, its
internal energy change and the direct PCM electrostatic
interaction is smaller than the cavitation, repulsion, and
dispersion corrections. Note that the use of ICAV, IREP,
and IDP are currently incompatible with gradients, so a
reasonable calculation sequence might be to perform the
geometry optimization with PCM electrostatics turned on,
then perform an additional calculation to include the other
solvent effects, adding extra functions to improve the
dispersion correction.
! calculation of CH4 (metano), in PCM water.
!
! This input reproduces the data in Table 2, line 6, of
! C.Amovilli, B.Mennucci J.Phys.Chem.B 101, 1051-7(1997)
! To do this, we must use many original PCM options.
!
! The gas phase FINAL energy is -40.2075980292
! The FINAL energy in PCM water is -40.2048210283
! (lit.)
! FREE ENERGY IN SOLVENT = -25234.89 KCAL/MOL
! INTERNAL ENERGY IN SOLVENT = -25230.64 KCAL/MOL
! DELTA INTERNAL ENERGY = .01 KCAL/MOL ( 0.0)
! ELECTROSTATIC INTERACTION = -.22 KCAL/MOL (-0.2)
! PIEROTTI CAVITATION ENERGY = 5.98 KCAL/MOL ( 6.0)
! DISPERSION FREE ENERGY = -6.00 KCAL/MOL (-6.0)
! REPULSION FREE ENERGY = 1.98 KCAL/MOL ( 2.0)
! TOTAL INTERACTION = 1.73 KCAL/MOL ( 1.8)
! TOTAL FREE ENERGY IN SOLVENT= -25228.91 KCAL/MOL
!
$contrl scftyp=rhf runtyp=energy $end
$guess guess=huckel $end
$system mwords=2 $end
! the "W1 basis" input here exactly matches HONDO's DZP
$DATA
CH4...gas phase geometry...in PCM water
Td
Carbon 6.
DZV
D 1 ; 1 0.75 1.0
Hydrogen 1. 0.6258579976 0.6258579976 0.6258579976
DZV 0 1.20 1.15 ! inner and outer scale factors
P 1 ; 1 1.00 1.0
$END
! The reference cited used a value for H2O's solvent
! radius that differs from the built in value (RSOLV).
! The IEF, ICOMP, MTHALL, and RET keywords are set to
! duplicate the original code's published results,
! namely D-PCM and GEPOL-GB. This run doesn't put in
! any "extra spheres" but we try that option (RET)
! like it originally would have.
$PCM SOLVNT=WATER RSOLV=1.35 RET=0.2
IEF=0 ICOMP=2 IDISP=0 IREP=1 IDP=1 ICAV=1 $end
$TESCAV MTHALL=1 $END
$NEWCAV IPTYPE=2 ITSNUM=540 $END
! dispersion "W2 basis" uses exponents which are
! 1/3 of smallest exponent in "W1 basis" of $DATA.
$DISBS NADD=11 NKTYP(1)=0,1,2, 0,1, 0,1, 0,1, 0,1
XYZE(1)=0.0,0.0,0.0, 0.0511
0.0,0.0,0.0, 0.0382
0.0,0.0,0.0, 0.25
1.1817023, 1.1817023, 1.1817023, 0.05435467
1.1817023, 1.1817023, 1.1817023, 0.33333333
-1.1817023, 1.1817023,-1.1817023, 0.05435467
-1.1817023, 1.1817023,-1.1817023, 0.33333333
1.1817023,-1.1817023,-1.1817023, 0.05435467
1.1817023,-1.1817023,-1.1817023, 0.33333333
-1.1817023,-1.1817023, 1.1817023, 0.05435467
-1.1817023,-1.1817023, 1.1817023, 0.33333333 $end
SVPE and SS(V)PE.
The Surface Volume Polarization for Electrostatics
(SVPE), and an approximation to SVPE called the Surface and
Simulation of Volume Polarization for Electrostatics
(SS(V)PE) are continuum solvation models. Compared to
other continuum models, SVPE and SS(V)PE pay careful
attention to the problems of escaped charge, the shape of
the surface cavity, and to integration of the Poisson
equation for surface charges.
The original references for what is now called the
SVPE (surface and volume polarization for electrostatics)
method are the theory paper:
"Charge penetration in Dielectric Models of Solvation"
D.M.Chipman, J.Chem.Phys. 106, 10194-10206 (1997)
and two implementation papers:
"Volume Polarization in Reaction Field Theory"
C.-G.Zhan, J.Bentley, D.M.Chipman
J.Chem.Phys. 108, 177-192 (1998)
"New Formulation and Implementation for Volume
Polarization in Dielectric Continuum Theory"
D.M.Chipman, J.Chem.Phys. 124, 224111-1/10 (2006)
which should be cited in any publications that utilize the
SVPE code.
There are two options to include with SS(V)PE or SVPE
additional models that describe short-range solute-solvent
interactions to achieve a more complete description of
solvation energies. Both have their keywords merged into
the $SVP input group for convenience. One option to
include short-range interactions is CMIRS1.0 (Composite
Method for Implicit Representation of Solvent, Version 1.0)
that combines the SS(V)PE dielectric continuum model with
the DEFESR (Dispersion, Exchange, and Field-Extremum Short-
Range) model. A complete account of CMIRS1.0 with
application to hydration energies is given in:
?Hydration Energy from a Composite Method for Implicit
Representation of Solvent?
A.Pomogaeva, D.M.Chipman
J.Chem.Theory Comput. 10, 211-219(2014)
which should be referenced in any publication that uses the
model. References to earlier works that develop the
individual components of the DEFESR parts of the full CMIRS
recipe are:
?Modeling short-range contributions to hydration
energies with minimal parameterization?
A.Pomogaeva, D.W.Thompson, and D.M.Chipman
Chem.Phys.Lett. 511, 161-165(2011).
?Field-Extremum Model for Short-Range Contributions
to Hydration Free Energy?
A.Pomogaeva and D.M.Chipman
J.Chem.Theory Comput. 7, 3952-3960(2011).
?New Implicit Solvation Models for Dispersion and
Exchange Energies?
A.Pomogaeva and D.M.Chipman
J.Phys.Chem.A, 117, 5812-5820 (2013).
The other option to include short-range interactions is the
SMVLE (solvation model with volume and local
electrostatics) as described in
"Free energies of solvation with surface, volume, and
local electrostatic effects and atomic surface
tensions to represent the first solvation shell"
J.Liu, C.P.Kelly, A.C.Goren, A.V.Marenich,
C.J.Cramer, D.G.Truhlar, C.-G. Zhan
J.Chem.Theory Comput. 6, 1109-1117(2010).
Further information on the performance of SVPE and of
SS(V)PE can be found in:
"Comparison of Solvent Reaction Field Representations"
D.M.Chipman, Theor.Chem.Acc. 107, 80-89 (2002).
Details of the SS(V)PE convergence behavior and programming
strategy are in:
"Implementation of Solvent Reaction Fields for
Electronic Structure" D.M.Chipman, M.Dupuis,
Theor.Chem.Acc. 107, 90-102 (2002).
The SMVLE option (solvation model with volume and
local electrostatics) is described in
"Free energies of solvation with surface, volume, and
local electrostatic effects and atomic surface tensions to
represent the first solvation shell" J.Liu, C.P.Kelly,
A.C.Goren, A.V.Marenich, C.J.Cramer, D.G.Truhlar, C.-G.
Zhan J.Chem.Theory Comput. 6, 1109-1117(2010).
The SVPE and SS(V)PE models are like PCM and COSMO in
that they treat solvent as a continuum dielectric residing
outside a molecular-shaped cavity, determining the apparent
charges that represent the polarized dielectric by solving
Poisson's equation. The main difference between SVPE and
SS(V)PE is in treatment of volume polarization effects that
arise because of the tail of the electronic wave function
that penetrates outside the cavity, sometimes referred to
as the "escaped charge." SVPE treats volume polarization
effects explicitly by including apparent charges in the
volume outside the cavity as well as on the cavity surface.
With a sufficient number of grid points, SVPE can then
provide an exact treatment of charge penetration effects.
SS(V)PE, like PCM and COSMO, is an approximate treatment
that only uses apparent charges located on the cavity
surface. The SS(V)PE equation is particularly designed to
simulate as well as possible the influence of the missing
volume charges. For more information on the similarities
and differences of the SVPE and SS(V)PE models with other
continuum methods, see the paper "Comparison of Solvent
Reaction Field Representations" cited just above.
In addition, the cavity construction and Poisson
solver used in this implementation of SVPE and SS(V)PE also
receive careful numerical treatment. For example, the
cavity may be chosen to be an isodensity contour surface,
and the Lebedev grids for the Poisson solver can be chosen
very densely. The Lebedev grids used for surface
integration are taken from the Fortran translation by C.
van Wuellen of the original C programs developed by D.
Laikov. They were obtained from the CCL web site
www.ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-
Grids. A recent leading reference is V. I. Lebedev and D.
N. Laikov, Dokl. Math. 59, 477-481 (1999). All these grids
have octahedral symmetry and so are naturally adapted for
any solute having an Abelian point group. The larger and/or
the less spherical the solute may be, the more surface
points are needed to get satisfactory precision in the
results. Further experience will be required to develop
detailed recommendations for this parameter. Values as
small as 110 are usually sufficient for simple diatomics
and triatomics. The default value of 1202 has been found
adequate to obtain the energy to within 0.1 kcal/mol for
solutes the size of monosubstituted benzenes. The SVPE
method uses additional layers of points outside the cavity.
Typically just two layers are sufficient to converge the
direct volume polarization contribution to better than 0.1
kcal/mol.
The SVPE and SS(V)PE codes both report the amount of
solute charge penetrating outside the cavity as calculated
by Gauss' Law. The SVPE code additionally reports the same
quantity as alternatively calculated from the explicit
volume charges, and any substantial discrepancy between
these two determinations indicates that more volume
polarization layers should have been included for better
precision. The energy contribution from the outermost
volume polarization layer is also reported. If it is
significant then again more layers should have been
included. However, these tests are only diagnostic.
Passing them does not guarantee that enough layers are
included.
The SVPE and SS(V)PE models treat the electrostatic
interaction between a quantum solute and a classical
dielectric continuum solvent. No treatment is yet
implemented for cavitation, dispersion, or any of a variety
of other specific solvation effects. Note that corrections
for these latter effects that might be reported by other
programs are generally not transferable. The reason is
that they are usually parameterized to improve the ultimate
agreement with experiment. In addition to providing
corrections for the physical effects advertised, they
therefore also inherently include contributions that help
to make up for any deficiencies in the electrostatic
description. Consequently, they are appropriate only for
use with the particular electrostatic model in which they
were originally developed.
Analytic nuclear gradients are not yet available for
the SVPE or SS(V)PE energy, but numerical differentiation
will permit optimization of small solute molecules.
Wavefunctions may be any of the SCF type: RHF, UHF, ROHF,
GVB, and MCSCF, or the DFT analogs of some of these. In the
MCSCF implementation, no initial wavefunction is available
so the solvation code does not kick in until the second
iteration.
We close with a SVPE example. The gas phase energy,
obtained with no $SVP group, is -207.988975, and the run
just below gives the SVPE energy -208.006282. The free
energy of solvation, -10.860 kcal/mole, is the difference
of these, and is quoted at the right side of the 3rd line
from the bottom of Table 2 in the paper cited. The
"REACTION FIELD FREE ENERGY" for SVPE is -12.905 kcal/mole,
which is only part of the solvation free energy. There is
also a contribution due to the SCRF procedure polarizing
the wave function from its gas phase value, causing the
solute internal energy in dielectric to differ from that in
gas. Evaluating this latter contribution is what requires
the separate gas phase calculation. Changing the number of
layers (NVLPL) to zero produces the SS(V)PE approximation
to SVPE, E= -208.006208.
! SVPE solvation test...acetamide
! reproduce data in Table 2 of the paper on SVPE,
! D.M.Chipman J.Chem.Phys. 124, 224111/1-10(2006)
!
$contrl scftyp=rhf runtyp=energy $end
$system mwords=4 $end
$basis gbasis=n31 ngauss=6 ndfunc=1 npfunc=1 $end
$guess guess=moread norb=16 $end
$scf nconv=8 $end
$svp nvlpl=3 rhoiso=0.001 dielst=78.304 nptleb=1202 $end
$data
CH3CONH2 cgz geometry RHF/6-31G(d,p)
C1
C 6.0 1.361261 -0.309588 -0.000262
C 6.0 -0.079357 0.152773 -0.005665
H 1.0 1.602076 -0.751515 0.962042
H 1.0 1.537200 -1.056768 -0.767127
H 1.0 2.002415 0.542830 -0.168045
O 8.0 -0.387955 1.310027 0.002284
N 7.0 -1.002151 -0.840834 -0.011928
H 1.0 -1.961646 -0.589397 0.038911
H 1.0 -0.752774 -1.798630 0.035006
$end
gas phase vectors, E(RHF)= -207.9889751769
$VEC
1 1 1.18951670E-06 1.74015997E-05
...snipped...
$END
The example just above can be changed to a CMIRS run, by
changing NVLPL to zero and also requesting calculation of
the DEFESR contributions. The result should be a CMIRS1.0
total free energy of -208.013517. The input should
(1) in $CONTRL, specify DFTTYP=HFX, which requests a
Hartree-Fock calculation, but sets up DFT grids,
(2) in $DFT, request a very accurate grid
$dft method=grid nrad=96 nleb=1202 nrad0=96 nleb0=1202 $end
(3) in $SVP, invoke the keyword IDEF=1.
Adding the keyword EGAS=-207.9889751769 to any of
these examples allows reporting of the free energy of
solvation (-10.860 kcal/mol for SVPE, -10.814 for SS(V)PE,
-15.400 for CMIRS1.0).
Conductor-like screening model (COSMO)
The COSMO (conductor-like screening model) represents a
different approach for carrying out polarized continuum
calculations. COSMO was originally developed by Andreas
Klamt, with extensions to ab initio computation in GAMESS
by Kim Baldridge.
In the COSMO method, the surrounding medium is modeled
as a conductor rather than as a dielectric in order to
establish the initial boundary conditions. The assumption
that the surrounding medium is well modeled as a conductor
simplifies the electrostatic computations and corrections
may be made a posteriori for dielectric behavior.
The original model of Klamt was introduced using a
molecular shaped cavity, which had open parts along the
crevices of intersecting atomic spheres. While having
considerable technical advantages, this approximation
causes artifacts in the context of the more generalized
theory, so the current method for cavity construction
includes a closure of the cavity to eliminate crevices or
pockets.
Two methodologies are implemented for treatment of the
outlying charge errors (OCE). The default is the well-
established double cavity procedure using a second, larger
cavity around the first one, and calculates OCE through the
difference between the potential on the inner and the outer
cavity. The second involves the calculation of distributed
multipoles up to hexadecapoles to represent the entire
charge distribution of the molecule within the cavity.
The COSMO model accounts only for the electrostatic
interactions between solvent and solute. Klamt has
proposed a novel statistical scheme to compute the full
solvation free energy for neutral solutes, COSMO-RS, which
is formulated for GAMESS by Peverati, Potier and Baldridge,
and is available as external plugin to the COSMOtherm
program by COSMOlogic GmbH&Co.
The iterative inclusion of non-electrostatic effects is
also possible right after a COSMO-RS calculation. The
DCOSMO-RS approach was implemented in GAMESS by Peverati,
Potier, and Baldridge, and more information is available on
Baldridge website at:
http://ocikbws.uzh.ch/gamess/
The simplicity of the COSMO model allows computation of
gradients, allowing optimization within the context of the
solvent. The method is programmed for RHF and UHF, all
corresponding kinds of DFT (including DFT-D), and the
corresponding MP2, energy and gradient.
Some references on the COSMO model are:
A.Klamt, G.Schuurman
J.Chem.Soc.Perkin Trans 2, 799-805(1993)
A.Klamt J.Phys.Chem. 99, 2224-2235(1995)
K.Baldridge, A.Klamt
J.Chem.Phys. 106, 6622-6633 (1997)
V.Jonas, K.Baldridge
J.Chem.Phys. 113, 7511-7518 (2000)
L.Gregerson, K.Baldridge
Helv.Chim.Acta 86, 4112-4132 (2003)
R.Peverati, Y.Potier, K.Baldridge
TO BE PUBLISHED SOON
Additional references on the COSMO-RS model, with
explanation of the methodology and program can be found:
A.Klamt, F.Eckert, W.Arlt
Annu.Rev.Chem.Biomol.Eng. 1, (2010)
created on 7/7/2017