The Polarizable Continuum Model (PCM)
The Polarizable Continuum Model (PCM) by Tomasi and coworkers is one of the most frequently used continuum solvation methods and has seen numerous variations over the years. The PCM model calculates the molecular free energy in solution as the sum over three terms:
Gsol = Ges + Gdr + Gcav (1)
These components represent the electrostatic (es) and the dispersion-repulsion (dr) contributions to the free energy, and the cavitation energy (cav). All three terms are calculated using a cavity defined through interlocking van der Waals-spheres centered at atomic positions. The reaction field is represented through point charges located on the surface of the molecular cavity (Apparent Surface Charge (ASC) model). The particular version of PCM that will be discussed here is the one using the United Atom for Hartree-Fock (UAHF) model to build the cavity. In this model the vdW-surface is constructed from spheres located on heavy (that is, non-hydrogen) elements only (United Atom approach). The vdW-radius of each atom is a function of atom type, connectivity, overall charge of the molecule, and the number of attached hydrogen atoms. In evaluating the three terms in equation (1) this cavity is used in slightly different ways.
While calculation of the cavitation energy Gcav uses the surface defined by the van der Waals-spheres, the solvent accessible surface is used to calculate the dispersion-repulsion contribution Gdr to the solution free energy. The latter differs from the former through additional consideration of the (idealized) solvent radius. The electrostatic contribution to the free energy in solution Ges uses an approximate version of the solvent excluding surface constructed through scaling all radii by a constant factor (e.g. 1.2 for water) and then adding some more spheres not centered on atoms in order to arrive at a somewhat smoother surface.
Localization and calculation of the surface charges is approached through systematic division of the spherical surface in small regions (tesserae) of known area and calculation of one point charge per surface element.
The implementation of the PCM/UAHF model in Gaussian 98 can be invoked using the SCRF keyword in combination with PCM-specific modifiers. The solvent can be specified using the Solvent= modifier to the SCRF keyword, acceptable solvent names being Water, DMSO, NitroMethane, Methanol, Ethanol, Acetone, DiChloroEthane, DiChloroMethane, TetraHydroFuran, Aniline, ChloroBenzene, Chloroform, Ether, Toluene, Benzene, CarbonTetrachloride, Cyclohexane, Heptane, and Acetonitrile. Additional options can be specified at the end of the input file and read in using the Read modifier of the SCRF keyword. The PCM solvation model is available for calculating energies and gradients at the Hartree-Fock and DFT levels. The output generated during PCM calculations can be dramatically extended using the DUMP option.
The following sample input illustrates the use of the corresponding keywords in the context of a single point calculation (without geometry optimization) on the aqueous solvation free energy of ethanol in its Cs symmetric conformation:
#P B3LYP/6-31G(d) scf=tight int=finegrid SCRF=(PCM,Read,Solvent=Water)
pcm/b3lyp/6-31G(d) sp ethanol in water (Cs)
0 1
O1
C2 1 r2
C3 2 r3 1 a3
H4 3 r4 2 a4 1 180.0
H5 3 r5 2 a5 4 d5
H6 3 r5 2 a5 4 -d5
H7 2 r7 3 a7 1 d7
H8 2 r7 3 a7 1 -d7
H9 1 r9 2 a9 3 180.0
r2=1.42492915
r3=1.51965095
r4=1.09569807
r5=1.09496362
r7=1.10264669
r9=0.96904984
a3=107.81130783
a4=110.63999342
a5=110.37205263
a7=109.90077195
a9=107.87777748
d5=-120.23659087
d7=-121.12750852
DUMP
|
 |
The additional output caused by the PCM solvation model is produced by link 502 responsible for the SCF calculation:
------------------------------------------------------------------------------
Solvent: WATER
Model : PCM/UAHF, Icomp = 4
Version: MATRIX INVERSION
Cavity : PENTAKISDODECAHEDRA with 60 initial tesserae
------------------------------------------------------------------------------
Nord Group Hybr Charge Alpha Radius Bonded to
1 OH sp3 0.00 1.20 1.590 C2 [s]
2 CH2 sp3 0.00 1.20 1.860 O1 [s] C3 [s]
3 CH3 sp3 0.00 1.20 1.950 C2 [s]
------------------------------------------------------------------------------
------------------------------------------------------
Dielectric Const = 78.39000
High.Fr.D.Const = 1.77600
d(Diel.Const.)/dT = -0.35620
Molar Volume = 18.07000
Therm.Exp.Coeff. = 0.00026
Radius = 1.38500
Absolute temper. = 298.00000
Number of spheres = 3
OMEGA = 40.00000
RET = 0.20000
FRO = 0.70000
Accuracy = 0.1D-05
------------------------------------------------------
The first four lines repeat settings specific for water or default settings for the PCM model in general. The solute cavity is constructed from vdW-spheres represented through regular pentakisdodecahedra, dividing each sphere's surface in 60 elements of equal size.
The following four lines list the results of the UAHF analysis, identifying only three (united atom) centers. For each center the assumed hybridization is listed together with its formal charge, the final radius, and the solvent-specific scaling parameter Alpha. The latter usually defaults to 1.2, but can be specified directly using the option
ALPHA=x.x
The final part of the output lists solvent-specific parameters such as the dielectric constant and the effective solvent radius, together with some more default PCM settings such as the number of initial spheres and the current values of parameters OMEGA, RET, and FRO. These latter three parameters govern the process of adding more spheres (not located at atomic positions) in order to smooth the surface. New values for OMEGA can be set using the option:
OMEGA=n.n
meaningful values ranging from 40.0 to 90.0 (higher values giving less added spheres). New values for FRO can be specified with:
FRO=m.m
meaningful values ranging from 0.7 to 0.2 (smaller values giving less added spheres). RET specifies the minimum radius of added new spheres and new values can be specified with:
RET=l.l
increasing values giving less added spheres and very large values avoid additional spheres completely.
The PCM algorithm then first runs through a gas phase energy calculation in order to obtain a reference point for the subsequent solvation free energy calculation. After completion of the gas phase SCF cycle, details of the iterative process of cavity generation are listed:
------------------------------------------------------
---------- CAVITY for ELECTROSTATIC term -----------
------------------------------------------------------
------- The SOLUTE is enclosed in ONE CAVITY -------
Total N.of Tesserae = 132
Surface Area (Ang**2) = 97.71529
Volume (Ang**3) = 84.60281
------------------------------------------------------
Original Sphere On Atom Re0 Alpha Surface
1 O1 1.590 1.200 24.08112
2 C2 1.860 1.200 27.26919
3 C3 1.950 1.200 46.36498
------------------------------------------------------
------------------------------------------------------------------------------
AT CONVERGENCE
132 Tesserae over a maximum of 1500
Surface Area (Ang**2) = 97.71529
Volume (Ang**3) = 84.60281
Escaped Charge= 0.13334
Error on NUCLEAR pol.charges = 0.21898 Error on ELECTR. pol.charges =-0.33812
------------------------------------------------------------------------------
------------------------------------------------------------------------------
dG(solv)/dEps (kcal/mol) = 0.00000
------------------------------------------------------
IN VACUO Dipole moment (Debye):
X= 0.0176 Y= 1.5625 Z= 0.0000 Tot= 1.5626
IN SOLUTION Dipole moment (Debye):
X= 0.1053 Y= 1.9429 Z= -0.0019 Tot= 1.9457
------------------------------------------------------
Tessera X Y Z QTot QSN QSE
1 2.84136 0.54870 3.42913 0.00354 -0.17977 0.18331
.
.
.
132 -2.70116 -2.26783 -4.13447 -0.00393 -0.25391 0.24997
In this (well behaved) case there is no need for additional spheres. The overall surface is represented by 132 surface elements ("Tesserae"). The molecular volume as defined through the current surface does not contain all the electron density of the system due to the long (in fact: infinite) tails of the electronic wavefunction, giving rise to some "Escaped Charge". At the center of each surface element sits one surface charge "QTot" containing one component from the nuclear charges of the solute and one from the electronic charges of the solute.
The results obtained during solution of the electronic Schroedinger equation (including the additional effects of the reaction field) are then given as:
SCF Done: E(RB+HF-LYP) = -155.041616090 A.U. after 10 cycles
Convg = 0.1222D-08 -V/T = 2.0093
S**2 = 0.0000
KE= 1.536131750668D+02 PE=-5.252140797608D+02 EE= 1.349508863800D+02
------------------------------------------------------
-------------- VARIATIONAL PCM RESULTS -------------
------------------------------------------------------
<Psi(0)| H |Psi(0)> (a.u.) = -155.033805
<Psi(0)|H+V(0)/2|Psi(0)> (a.u.) = -155.040760
<Psi(0)|H+V(f)/2|Psi(0)> (a.u.) = -155.041609
<Psi(f)| H |Psi(f)> (a.u.) = -155.032883
<Psi(f)|H+V(f)/2|Psi(f)> (a.u.) = -155.041616
Total free energy in sol.
(with non electrost.terms) (a.u.) = -155.041162
------------------------------------------------------
(Unpol.Solute)-Solvent (kcal/mol) = -4.36
(Polar.Solute)-Solvent (kcal/mol) = -5.48
Solute Polarization (kcal/mol) = 0.58
Total Electrostatic (kcal/mol) = -4.90
------------------------------------------------------
Cavitation energy (kcal/mol) = 8.92
Dispersion energy (kcal/mol) = -11.39
Repulsion energy (kcal/mol) = 2.75
Total non electr. (kcal/mol) = 0.29
------------------------------------------------------
DeltaG (solv) (kcal/mol) = -4.62
------------------------------------------------------
What is listed here as:
<Psi(0)| H |Psi(0)> (a.u.) = -155.033805
is the unperturbed gas phase SCF solution used as the reference for all subsequent steps. The following energy described as:
<Psi(0)|H+V(0)/2|Psi(0)> (a.u.) = -155.040760
includes the interaction of the unpolarized solute with the unpolarized solvent. Comparison with the gas phase reference energy yields the corresponding interaction energy:
(Unpol.Solute)-Solvent (kcal/mol) = -4.36
After reporting the total energy corresponding to the interaction of the unpolarized solute with the polarized solvent, the next important information is that on the total energy of the polarized solute:
<Psi(f)| H |Psi(f)> (a.u.) = -155.032883
The energy difference to the umpolarized gas phase total energy is listed as:
Solute Polarization (kcal/mol) = 0.58
and should always be positive. The fully interacting system of polarized solvent with polarized solute gives rise to total energy:
<Psi(f)|H+V(f)/2|Psi(f)> (a.u.) = -155.041616
and the corresponding electrostatic interaction energy is listed as:
Total Electrostatic (kcal/mol) = -4.90
The non-electrostatic parts of the solvation free energy are given together in one block terminated by the sum of the cavitation and the dispersion-repulsion energies:
Total non electr. (kcal/mol) = 0.29
The sum of the electrostatic and non-electrostatic contributions gives the overall free energy of solvation:
DeltaG (solv) (kcal/mol) = -4.62
It should be recognized, however, that the solvation free energy calculated here refers to a motionless system in the gas phase at 0 Kelvin. In order to arrive at thermochemically meaningful free energies at finite temperatures these solvation free energies have to be augmented with a standard treatment of gas phase thermochemistry. For the calculation of the solvation free energy itself as the difference in gas and solution free energies this is often neglected and the value produced by PCM is used directly in comparison to experimental values. For ethanol, the experimental free energy of solvation has been measured as -5.0 kcal/mol. Compared to this value the PCM prediction of -4.6 kcal/mol can be considered to be rather accurate.
When calculating solvation free energies as the difference of gas phase and solution phase free energies attention must be paid to the definition of the respective standard states. In case both the gas and solution phase concentrations are given in molar values (mol/l), gas and solution phase data can be compared directly. Gas phase values refer, however, often to a partial pressure of 1 atm. Assuming ideal gas behaviour, this corresponds to 1/24.46 mol/l at 298.15K.
An improved prediction of solvation free energies should also cover the effects of structural relaxation in solution. Geometry optimizations using the PCM model are possible, but much more time consuming than gas phase optimizations. This is not only due to higher CPU times for each of the energy and gradient calculations, but also due to a much slower convergence of the optimization process and frequent oscillations. Two options that are helpful in alleviating some of the convergence problems are TSNUM and TSARE:
TSNUM specifies the number of surface elements (tesserea) for each sphere. The PCM algorithm selects regular polyhedra whose number of surface elements are as close as possible to TSNUM. Aside from the default value of 60, some larger values such as 64, 80, or 100 may be helpful in reducing the oscillatory behaviour in some geometry optimizations.
TSARE specifies the area of the surface elements in units of (Angstroms)2. Meaningful values range from 0.4 to 0.2, smaller values leading to a larger number of surface elements. Setting the size of the surface elements to a particular value leads to surface elements of equal size, regardless of the radius of the sphere (this is not the case when TSNUM settings are changed).
In both cases, however, the total energies depend on the actual choice of surface elements and a comparison of results for different systems or different conformers is only meaningful with the same choice of options. For the ethanol example used here the geometry optimization does not converge within 23 steps using the default settings, but can be brought to convergence within 7 steps using TSARE=0.3, yielding a final total energy in solution of -155.043097702 au. Compared to a PCM single point calculation using the gas phase geometry (with total energy of -155.041371 au) this implies a structural relaxation energy of -1.1 kcal/mol and thus an "improved" prediction of the solvation free energy of -5.7 kcal/mol.
One problem when applying the current implementation of the PCM/UAHF model to reaction pathways in solution is a direct consequence of using hybridization and connectivity data for the derivation of the vdW radii. As both hybridization and connectivity are not going to change smoothly but suddenly from one point to the next along a reaction pathway, the UAHF approach is bound to produce sudden breaks in the free energy of solvation profile. These problems can, in principle, be circumvented by smoothly scaling from one set of radii to another, or by avoiding the UAHF approach altogether and choosing radii that don't depend on connectivity. One common choice for the latter case are Pauling radii available with the option RADII=Pauling.
More recent versions of Gaussian contain a strongly revised version of the PCM model, in particular with respect to the definition of the solute cavity. This has also brought some changes to the list of acceptable keywords. The following example illustrates the keywords PCMDOC (formerly DUMP), RADII, and SCFVAC.
#P B3LYP/6-31G(d) scf=tight int=finegrid SCRF=(PCM,Read,Solvent=Water)
pcm/b3lyp/6-31G(d) sp ethanol in water (Cs)
0 1
O1
C2 1 r2
C3 2 r3 1 a3
H4 3 r4 2 a4 1 180.0
H5 3 r5 2 a5 4 d5
H6 3 r5 2 a5 4 -d5
H7 2 r7 3 a7 1 d7
H8 2 r7 3 a7 1 -d7
H9 1 r9 2 a9 3 180.0
r2=1.42492915
r3=1.51965095
r4=1.09569807
r5=1.09496362
r7=1.10264669
r9=0.96904984
a3=107.81130783
a4=110.63999342
a5=110.37205263
a7=109.90077195
a9=107.87777748
d5=-120.23659087
d7=-121.12750852
PCMDOC
RADII=UAHF
SCFVAC
|
|
The last of the keywords forces the program to first calculate the unperturbed gas phase wavefunction and afterwards the one in solution. While this is the default in earlier versions of Gaussian, the current default is to only perform the solution calculation. A second major change concerns the choice of default radii. While UAHF radii have been the default choice in earlier versions of Gaussian, the newer version use United Atom Topological Model (UA0) radii by default. This is not necessarily the best choice and some experimentation is required as to which electronic structure method fits which type of radii best. United Atom radii optimized to fit density functional methods particularly well have been optimized for the PBE/6-31G(d) level of theory and can be used with RADII=UAKS. The TSARE keyword has retained its previous meaning, but defaults now to a value of 0.2 (Angstroms)2.