Continuum Solvent Model (CSM)

Theoretical Background

A continuum solvent model (CSM) has been added to the GPAW implementation as described in Reference [1]. 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

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

\[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:

\[\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 [2].

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 [1] 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. [3] is taken, where the surface area \(A\) is calculated from the gradient of the cavity:

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

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

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

\[\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 [1].

Altogether, 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 [1]. 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 [1].

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 [1] 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:

from gpaw import GPAW
from ase.build import molecule
from ase.units import mol, kJ, kcal, Pascal, m
from ase.data.vdw 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 func 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.copy()
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)


def atomic_radii(atoms):
    return [vdw_radii[n] for n in atoms.numbers]


# create Atoms object for ethanol and add vacuum
atoms = molecule('CH3CH2OH')
atoms.center(vacuum=vac)

# perform gas phase calculation
atoms.calc = GPAW(mode='fd', xc='PBE', h=h, txt='gasphase.txt')
Egasphase = atoms.get_potential_energy()

# perform calculation with continuum solvent model from
# J. Chem. Phys. 141, 174108 (2014)
atoms.calc = SolvationGPAW(
    mode='fd', xc='PBE', h=h, txt='water.txt',
    cavity=EffectivePotentialCavity(
        effective_potential=Power12Potential(atomic_radii, u0),
        temperature=T,
        surface_calculator=GradientSurface()),
    dielectric=LinearDielectric(epsinf=epsinf),
    interactions=[SurfaceInteraction(surface_tension=gamma)])
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))

The calculated value for the solvation Gibbs energy should be about -4.5 kcal per mole. The experimental value is -5.0 kcal per mole [4]. Please refer to the solvation module source code for further reading about the usage of the SolvationGPAW calculator class or model specific parts.

There is also a helper function to use the solvation parameters for water as in the above example:

import gpaw.solvation as solv

# ...

# convenient way to use HW14 water parameters:
calc = solv.SolvationGPAW(
    mode='fd', xc='PBE', h=0.2,  # non-solvent DFT parameters
    **solv.get_HW14_water_kwargs()
)

# ...

class gpaw.solvation.calculator.SolvationGPAW(restart=None, cavity=None, dielectric=None, interactions=None, **gpaw_kwargs)

Subclass of gpaw.GPAW calculator with continuum solvent model.

See also Section III of A. Held and M. Walter, J. Chem. Phys. 141, 174108 (2014).

Constructor for SolvationGPAW class.

Additional arguments not present in GPAW class: cavity – A Cavity instance. dielectric – A Dielectric instance. interactions – A list of Interaction instances.

References

[1] (1,2,3,4,5,6)

A. Held and M. Walter, Simplified continuum solvent model with a smooth cavity based on volumetric data, J. Chem. Phys. 141, 174108 (2014).

[2]

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

[3]

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

[4]

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)