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
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 

   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
This gives numerical parameters, all built into GAMESS, as 
the various SOLX values, for SOLVNT=
  NBA       MI        C6H5I     OCTANE    ET3N    
  CS2       DMA       CH3I      PENTANAL  UNDECANE
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 
   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 

                        * * *

   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.)
!  DELTA INTERNAL ENERGY       =       .01 KCAL/MOL ( 0.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)
 $contrl scftyp=rhf runtyp=energy $end
 $guess  guess=huckel $end
 $system mwords=2 $end
!    the "W1 basis" input here exactly matches HONDO's DZP
CH4...gas phase PCM water

Carbon      6.
   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

!    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.
         IEF=0 ICOMP=2 IDISP=0 IREP=1 IDP=1 ICAV=1 $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


     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
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 

      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 

      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 

     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
CH3CONH2 cgz geometry RHF/6-31G(d,p)
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
gas phase vectors, E(RHF)=     -207.9889751769
 1  1 1.18951670E-06 1.74015997E-05

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 

    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:


    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
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