Commit a11607a1 authored by aheld's avatar aheld

solvation: adds documentation for continuum solvent model

* adds a continuum solvent model tutorial
* updates projects page in devel section
* adds continuum solvent feature to index page
* adds citation information to FAQ
* adds reference to doc/GPAW.bib
parent 91e1cae0
......@@ -138,3 +138,15 @@ grid and basis representations.}
doi = {10.1103/PhysRevB.87.235132},
url = {}
author = {Held, Alexander and Walter, Michael},
title = {Simplified continuum solvent model with a smooth cavity based on volumetric data},
journal = {J. Chem. Phys.},
year = {2014},
volume = {141},
number = {17},
pages = {174108},
url = {},
doi = {10.1063/1.4900838}
Continuum Solvent Model (CSM)
Alexander Held
See :ref:`continuum_solvent_model`.
Polarizable continuum model (PCM)
Alexander Held
A polarizable continuum model according to [1]_ is being implemented
(still to be checked into the repository).
Recent features:
* extend bmgs solver to deal with global coefficients
* new poisson solver for poisson equation with dielectric
* cavity from van der Waals radii
* calculate surface energy from cavity
Planned features:
* forces
* charged systems
* cavity from density cut-off
.. [1] V. M. Sanchez, M. Sued and D. A. Scherlis,
J. Chem. Phys. 131, 174108 (2009).
......@@ -24,6 +24,6 @@ Ongoing Projects
......@@ -62,6 +62,14 @@ If you use the :ref:`gw_tutorial`, please cite also:
If you use the :ref:`continuum_solvent_model`, please cite also:
| A. Held and M. Walter
| `Simplified continuum solvent model with a smooth cavity based on volumetric data`__
| The Journal of Chemical Physics Vol. **141**, 174108, 2014
BibTex (:svn:`doc/GPAW.bib`):
.. literalinclude:: GPAW.bib
......@@ -38,7 +38,7 @@ Quick links to all features:
- :ref:`GW <gw_theory>`
- :ref:`BSE <bse>`
* - :ref:`Parallelization <parallel_runs>`
- :ref:`Continuum Solvent Model <continuum_solvent_model>`
Watch us on GPAWTV
.. _continuum_solvent_model:
Continuum Solvent Model (CSM)
Theoretical Background
A continuum solvent model (CSM) has been added to the GPAW
implementation as described in Reference [#HW14]_. The model is based
on a smooth cavity `g(\br)` represented on a real space grid. The
electrostatic screening effect of the solvent is modeled as a
spatially varying permittivity (relative dielectric constant)
`\epsilon(\br)`, which is related to the cavity by
.. math:: \epsilon(\br) = 1 + (\epsilon_{\infty} - 1) g(\br).
The bulk static dielectric constant `\epsilon_{\infty}` is taken from
experiment for the solvent in use (e.g. `\epsilon_{\infty} \approx 80`
for water at room temperature). The electrostatic energy for the total
charge distribution `\rho(\br)` of the solute (electron density +
atomic nuclei) is calculated as
.. math:: W_{\mathrm{el}} = \frac{1}{2} \int \rho(\br) \Phi(\br) \mathrm{d}\br
where `\Phi(\br)` is the electrostatic potential in the presence of
the dielectric. It is obtained by solving the Poisson equation
including the spatially varying permittivity:
.. math:: \nabla \Big(\epsilon(\br) \nabla \Phi(\br)\Big) = -4\pi \rho(\br)
This is done numerically with a finite difference multi-grid solver as
outlined in Reference [#SSS09]_.
As not only the electrostatic interaction but also cavity formation
and short-range interactions between the solute and the solvent
contribute to the solvation Gibbs energy, additional short-range
interacitons can be included in the calculation. The ansatz chosen in
Reference [#HW14]_ is to describe the cavity formation energy and all
short-range interactions as a term proportional to the surface area of
the cavity. The surface area is not well defined for a smooth cavity
`g(\br)`. The approach of Im. *et al.* [#IBR98]_ is taken, where the
surface area `A` is calculated from the gradient of the cavity:
.. math:: A = \int \| \nabla g(\br) \| \mathrm{d}\br
The cavity formation energy and short-range interaction is then given
as `G_{\mathrm{cav+sr}} = \gamma A` with a proportionality constant
`\gamma` that has the dimensions of a surface tension.
The cavity is constructed from an effective repulsive potential
`u(\br)` via the Boltzmann equation
.. math:: g(\br) = \exp \Big(-\frac{u(\br)}{k_{\mathrm{B}} T}\Big).
The effective potential is taken as the repulsive part of the
Lennard-Jones potential and constructed from the positions
`\mathbf{R}_a` of the atomic nuclei and a set of atomic radii
`R_a^{\mathrm{vdW}}` as
.. math:: u(\br) = u_0 \sum_a \Big(\frac{R_a^{\mathrm{vdW}}}{\|\br - \mathbf{R}_a \|} \Big)^{12}.
The single free parameter `u_0` describes the strength of the repulsion at
the atomic radii. It controls the size of the cavity and can be fitted
to experimental partial molar volumes (per atom) `\bar{v}_M^{\infty}`
at infinite dilution via the compressibility equation
.. math:: \int (1 - g(\br)) \mathrm{d}\br = \bar{v}_M^{\infty} - \kappa_T k_{\mathrm{B}} T,
where `\kappa_T` is the isothermal compressibility of the bulk solvent
(usually negligible compared to `\bar{v}_M^{\infty}` except for very
small molecules). For water as a solvent, Bondi van der Waals radii
with a modified value for hydrogen (1.09 Å) are a good
choice for the atomic radii of the effective repulsive potential
Alltogether, the model has three parameters. The static dielectric
constant `\epsilon_{\infty}` of the solvent is taken directly from
experimental data. The parameter `u_0` is fitted to experimental
volumes. Afterwards, the parameter `\gamma` can be fitted to
experimental solvation Gibbs energies. Parameters for water (at room
temperature) as a solvent are given in Reference [#HW14]_. A
preliminary parametrization for dimethylsulfoxide (DMSO) is also given
there. The accuracy for aqueous solvation Gibbs energies compared to
experimental values is about 5 kJ per mole for neutral molecules and
13 kJ per mole for cations. It is not recommended to apply the model
to anions in water without further modification as short-range
interactions between anions and water are not well represented using
the parameters optimized for neutral molecules [#HW14]_.
A lot of the parts that make up the model (cavity, dielectric
function, Poisson solver, non-electrostatic interactions) can be
replaced easily by alternative implementations as they are represented
as python classes in the source code. Some alternative models already
exist in the GPAW source code, yet they are not well tested and
therefore not recommended for production use. Nevertheless they can
serve as an example on how to add your own solvation
interactions/models to GPAW. Also refer to Section III of Reference
[#HW14]_ for a more precise description of the general framework of
continuum solvent models with a smooth cavity employed in the GPAW
implementation. In the following section, the usage of the model
described above is documented in terms of a use case.
Usage Example: Ethanol in Water
As a usage example, the calculation of the solvation Gibbs energy of
ethanol in water is demonstrated. For simplicity, the solvation Gibbs
energy is calculated as the total energy difference of a calculation
with the continuum solvent model and a gas phase calculation using a
fixed geometry of the ethanol molecule. In principle, a relaxation of
the solute can be done also with the continuum solvent model. The
following annotated script demonstrates how to perform the calculation
using the continuum solvent model:
.. literalinclude::
The calculated value for the solvation Gibbs energy should be about
-4.3 kcal per mole. The experimental value is -5.0 kcal per mole
[#KCT06]_. Please refer to the Epydoc documentation of the
:epydoc:`gpaw.solvation` module or the source code for further reading
about the usage of the :class:`~gpaw.solvation.SolvationGPAW`
calculator class or model specific parts.
.. [#HW14] A. Held and M. Walter,
`Simplified continuum solvent model with a smooth cavity based on volumetric data <>`_,
*J. Chem. Phys.* **141**, 174108 (2014).
.. [#SSS09] V. M. Sanchez, M. Sued and D. A. Scherlis,
`First-principles molecular dynamics simulations at solid-liquid interfaces with a continuum solvent <>`_,
*J. Chem. Phys.* **131** 174108 (2009).
.. [#IBR98] W. Im, D. Beglov and B. Roux,
`Continuum solvation model: Computation of electrostatic forces from numerical solutions to the Poisson-Boltzmann equation <>`_,
*Comput. Phys. Commun.* **111** 59 (1998).
.. [#KCT06] C. P. Kelly, C. J. Cramer and D. G. Truhlar,
`Aqueous Solvation Free Energies of Ions and Ion-Water Clusters Based on an Accurate Value for the Absolute Aqueous Solvation Free Energy of the Proton <>`_,
*J. Phys. Chem. B* **110** 16066 (2006)
from gpaw import GPAW
from ase.structure import molecule
from ase.units import mol, kJ, kcal, Pascal, m
from import vdw_radii
from ase.parallel import parprint
from gpaw.solvation import (
SolvationGPAW, # the solvation calculator
EffectivePotentialCavity, # cavity using an effective potential
Power12Potential, # a specific effective potential
LinearDielectric, # rule to construct permittivity function from the cavity
GradientSurface, # rule to calculate the surface area from the cavity
SurfaceInteraction # rule to calculate non-electrostatic interactions
# all parameters on the user side of the solvation API follow the ASE unit conventions (eV, Angstrom, ...)
# non-solvent related DFT parameters
h = 0.2
vac = 5.0
# solvent parameters for water from J. Chem. Phys. 141, 174108 (2014)
u0 = 0.180 # eV
epsinf = 78.36 # dimensionless
gamma = 18.4 * 1e-3 * Pascal * m # convert from dyne / cm to eV / Angstrom ** 2
T = 298.15 # Kelvin
vdw_radii = vdw_radii[:]
vdw_radii[1] = 1.09
# atomic_radii expected by gpaw.solvation.Power12Potential have to be a callable
# mapping an Atoms object to an iterable of floats representing the atomic radii of
# the atoms in the same order as in the Atoms object (in Angstrom)
atomic_radii = lambda atoms: [vdw_radii[n] for n in atoms.numbers]
# create Atoms object for ethanol and add vacuum
atoms = molecule('CH3CH2OH')
# perform gas phase calculation
atoms.calc = GPAW(xc='PBE', h=h)
Egasphase = atoms.get_potential_energy()
# perform calculation with continuum solvent model from J. Chem. Phys. 141, 174108 (2014)
atoms.calc = SolvationGPAW(
xc='PBE', h=h,
effective_potential=Power12Potential(atomic_radii, u0),
Ewater = atoms.get_potential_energy()
# calculate solvation Gibbs energy in various units
DGSol_eV = Ewater - Egasphase
DGSol_kJ_per_mol = DGSol_eV / (kJ / mol)
DGSol_kcal_per_mol = DGSol_eV / (kcal / mol)
parprint('calculated Delta Gsol = %.0f meV = %.1f kJ / mol = %.1f kcal / mol' % (DGSol_eV * 1000., DGSol_kJ_per_mol, DGSol_kcal_per_mol))
......@@ -51,4 +51,4 @@ Specialized tutorials
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment