CosmoSolvationParameters

class CosmoSolvationParameters(solvent_dielectric_constant=None, solvent_surface_tension=None, grid_spacing=None, atomic_radii=None, atom_grid_density=None, hydrogen_grid_density=None, minimum_grid_point_area=None, shifted_force_cutoff=None, shifted_force_damping=None)

Container class for the parameters used in COSMO solvent model calculations.

Parameters:

• solvent_dielectric_constant (float | Conductor) – The dielectric constant of the solvent. The flag Conductor assumes the solvent is a perfect conductor and gives an infinite dielectric constant
  Default: Conductor.

• solvent_surface_tension (PhysicalQuantity of type force per distance) – The surface tension of the surrounding solvent. Used to calculate the cavitation energy.

• grid_spacing (PhysicalQuantity of type length) – The spacing of the grid used to map the curvature of the cavity.
  Default: 0.3 Angstrom

• atomic_radii (dict) – A dictionary with custom atomic radii to use for different elements. Keys and values are PeriodicTableElements and PhysicalQuantity respectively.

• atom_grid_density (int) – The density of points used to generate the COSMO surface segments for non-hydrogen atoms. The total number of grid points is given as \(10n^2 + 2\).
  Default: 4, which corresponds to 162 grid points.

• hydrogen_grid_density (int) – The density of points used to generate the COSMO surface segments for hydrogen atoms. The total number of grid points is given as \(10n^2 + 2\).
  Default: 3, which corresponds to 92 grid points.

• minimum_grid_point_area (PhysicalQuantity of type area.) – The minimum allowed area associated with each grid point. Grid points with areas less than this limit are merged into surrounding grid points.
  Default: 0.1 Angstrom**2

• shifted_force_cutoff (PhysicalQuantity of type distance) – The cutoff used for estimating electrostatic interactions with the shifted-force potential. This is only used in surfaces in which shifted-force potentials are used.
  Default: 10 Angstrom.

• shifted_force_damping (PhysicalQuantity of type inverse distance) – The damping used for estimating electrostatic interactions with the shifted-force potential. This is only used in surfaces in which shifted-force potentials are used.
  Default: 0.1 Angstrom**-1.

atomGridDensity()

Returns:

The density of points used to generate the COSMO surface segments for non-hydrogen atoms.

Return type:

int

atomicRadii()

Returns:

A dictionary with custom atomic radii to use for different elements.

Return type:

dict | None

gridSpacing()

Returns:

The grid spacing of the tesserae grid.

Return type:

PhysicalQuantity of type length

hydrogenGridDensity()

Returns:

The density of points used to generate the COSMO surface segments for hydrogen atoms.

Return type:

int

minimumGridPointArea()

Returns:

The minimum allowed grid point area.

Return type:

PhysicalQuantity of type area

nlinfo()

Returns:

Calculator containing the info.

Return type:

dict

shiftedForceCutoff()

Returns:

The cutoff used in shifted-force surface electrostatic interactions.

Return type:

PhysicalQuantity of type length

shiftedForceDamping()

Returns:

The damping used in shifted-force surface electrostatic interactions.

Return type:

PhysicalQuantity of type length

solventDielectricConstant()

Returns:

The dielectric constant of the solvent.

Return type:

float

solventSurfaceTension()

Returns:

The surface tension of the solvent.

Return type:

float

uniqueString()

Return a unique string representing the state of the object.

\(\renewcommand\AA{\text{Å}}\)

Usage Examples

Calculate the optimized geometry and energy of a glycine zwitter ion in water.

# Create LCAOCalculator
solvation_parameters = CosmoSolvationParameters(
    solvent_dielectric_constant=78.3553,
    solvent_surface_tension=72.8*dyne/cm
)
calculator = LCAOCalculator(
    solvation_parameters=solvation_parameters
)

# Set Calculator
configuration.setCalculator(calculator)
configuration.update()
nlsave('glycine_zwitter_ion.hdf5', configuration)

# Optimize the geometry
optimized_configuration = OptimizeGeometry(
    configuration=configuration,
    trajectory_filename='glycine_zwitter_ion.hdf5'
)
nlsave('glycine_zwitter_ion.hdf5', optimized_configuration)

glycine_zwitter_ion.py

Notes

This object contains the parameters used in a COnductor-like Screening MOdel (COSMO) calculation[1][2]. COSMO is a continuum solvation method that can be added to an LCAO-DFT calculation to model the effect of a solvent surrounding the molecule or interacting with a surface. This is done by creating a discretized surface mesh around the molecule or surface. This mesh is then polarized according to the electrostatic potential from the molecule and the dielectric constant of the solvent. During the SCF optimization of the electron density, the potential from the cavity mesh is included. Similarly, in each SCF step the new electron density is also used to update the mesh charges. In this way the molecule and surrounding solvent are mutually polarized to a self consistent solution.

Mathematically, the total COSMO dielectric energy is given as:

\[E = \frac{1}{2}\mathbf{q}A\mathbf{q} + f(\epsilon)\left(\mathbf{z}B\mathbf{q} + \mathbf{c}\mathbf{q}\right)\]

Here the vector \(\mathbf{q}\) represents the charges on the cavity mesh. The matrix \(A\) is the Coulomb matrix giving the interaction energy between the mesh charges. Likewise \(B\) is a Coulomb matrix giving the interaction energy between the atomic nuclear charges and the mesh charges. The vector \(\mathbf{c}\) completes the electrostatic interactions giving the interaction between the electron density within the cavity and the mesh charges. This is given as:

\[c_i = \int_V \frac{\rho(\mathbf{r})}{|\mathbf{r}_i-\mathbf{r}|} dr\]

Finally \(f(\epsilon)\) is a scaling function of the dielectric constant \(\epsilon\). This is given as:

\[f(\epsilon) = \frac{\epsilon - 1}{\epsilon + 0.5}\]

This allows the mesh charges can be calculated by minimizing the dielectric energy. This gives the equation for the charges

\[\mathbf{q} = f(\epsilon) A^{-1}\left(\mathbf{z}B + \mathbf{c}\right)\]

An additional correction is added to the charges to account for the electron density outside the solvent cavity. The correction assumes that the total charge of the cavity surface is the negative of the system charge multiplied by the scaling function. This is implemented by adding a Lagrange multiplier to solving the surface charges that constrains the total charge.

Analytic atomic forces for COSMO are also available. These are calculated by taking the derivatives of the dielectric energy. This enables COSMO to be included in geometry optimizations which can have a significant effect in highly charged configurations, such as molecular zwitter ions. Here the inclusion of the COSMO can stabilize the zwitter ion form, making this more energetically favorable then the neutral molecule.

In the CosmoSolvationParameters the dielectric constant \(\epsilon\) is given using the solvent_dielectric_constant argument. This controls the strength of the polarization of the solvent around the molecule, and can be used to approximate different solvents. In the COSMO-RS model, the surrounding solvent is assumed to be a perfect conductor. This can be specified with the flag Conductor, which gives an infinite dielectric constant for the solvent corresponding to \(f(\epsilon) = 1\).

In addition to the dielectric energy, an estimate of the non-electrostatic energy of solvation, including the cavitation and dispersion-repulsion energies can be included. This is done with a linear regression model that includes the solvent cavity surface area and surface tension as parameters. The solvent surface tension is specified with the argument solvent_surface_tension. If a solvent surface tension is not given, the non-electrostatic energy is not calculated, and simply given as zero. This is appropriate when performing calculations the the COSMO-RS model.

One of the important parts of the COSMO method is the solvent cavity construction. In QuantumATK the cavity is constructed by first calculating an pseudo-density isosurface around the solvent facing atoms[3]. This gives a smooth surface that approximates the solvent-accessible surface. Grid points from each atom are then projected out to this pseudo-density isosurface to form the mesh used to solve the COSMO surface charges. Surface area is also assigned to each grid point based on their position on the isosurface.

The construction of COSMO solvent cavity can be modified in a number of ways. The spacing of the grid used to calculate the initial pseudo-density is given with the argument grid_spacing. Lowering the distance between grid points places the calculated isodensity surface closer to the true isosurface, at the cost of more computation time and memory. The default value is sufficient for most COSMO calculations. The atomic radii, which determines the size of the COSMO cavity, can also be specified with the atomic_radii argument. This takes a dictionary of elements and corresponding atomic radii for the elements whose radii are to be modified. Here again the default radii are optimized for the COSMO method and are suitable for most calculations [4].

Once the pseudo-density surface is created, the COSMO grid points are created by projecting points from each atom out to the surface. The number of these projected points can be set, altering their density on the surface. Different densities are allowed for hydrogen and non-hydrogen atoms, as hydrogen is significantly smaller than other atoms. The arguments to set the number of points are hydrogen_grid_density and atom_grid_density. These each take an integer that specifies how many times the atomic grid is divided. For a given number \(n\), the number of grid points is \(10n^2 + 2\). For hydrogen and non-hydrogen atoms the defaults are 3 and 4 respectively, which corresponds to 92 and 162 points. Increasing the grid density increases the accuracy of the surface charge distribution at the cost of more memory and computational time. As with the other parameters, the defaults here should be sufficient for most calculations. When the points are projected from the atoms, they are assigned a surface area based on the segments of the pseudo-density surface that are closest to the point. Points with a surface area below a set limit are removed, and their surface area is assigned to surrounding grid points. This limit can be set with the argument minimum_grid_point_araea. By default this is set to 0.1\(\AA^2\) . Lowering this limit can cause numerical problems with large surface charges resulting in the COSMO calculation.

For surfaces, to account for the 2D periodicity of the structure a damped shifted-force potential is used in place of the normal Coulomb potential[5]. The equation for this potential is:

\[\DeclareMathOperator\erfc{erfc} V(r) = q_iq_j \left[ \frac{\erfc(\alpha r)}{r} - \frac{\erfc(\alpha R_c)}{R_c} + \left( \frac{\erfc(\alpha R_c)}{R_c^2} + \frac{2\alpha \exp(-\alpha^2 R_c^2)}{\sqrt{\pi} R_c} \right) (r-R_c) \right]\]

Here \(r\) is the distance between the points, \(R_c\) is a cutoff, beyond which the potential is zero and \(\alpha\) is a damping factor. The cutoff and damping factor can be set using the arguments shifted_force_cutoff and shifted_force_damping, respectively.

When performing COSMO calculations the boundary conditions of the simulation can be important. For neutral molecules it is generally accurate enough to use 3D-FFT boundary conditions. For charged molecules it is recommended to use Multipole boundary conditions with the conjugate gradient Poisson solver ParallelConjugateGradientSolver. These boundary conditions correctly model the overall charge of the molecule.

In bulk slab systems, the surface needs to be orientated in the C direction of the unit cell. It is then assumed that the system is only periodic in the A and B direcations. A COSMO surface is added only to the top layer of the surface, as it is assumed that the lower face continues into the bulk material. It is recommended that the bulk slabs 2D-FFT boundary conditions are used for A and B, and Dirichlet and Neumann are used for the bottom and top, respectively, for the C direction. With these boundary conditions the slab should be in the middle of the cell, so that the potential value and gradient have space to reach zero at each end of the cell. The same boundary conditions should also be used with COSMO calculations using a SurfaceConfiguration.

[1]

A. Klamt and G. Schüürmann. COSMO: A new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2, 99(7):799–805, feb 1993. doi:10.1039/P29930000799.

[2]

Andreas Klamt. The COSMO and COSMO-RS solvation models. WIREs Comput Mol Sci, 1(7):699, sep 2011. doi:10.1002/wcms.56.

[3]

Andreas Klamt and Michael Diedenhofen. A refined cavity construction algorithm for the conductor-like screening model. J. Comp. Chem., 39(21):1648–1655, 2018. doi:10.1002/jcc.25342.

[4]

Andreas Klamt. Refinements, parameterization, and the complete COSMO-RS, pages 109–125. Elsevier, 2005.

[5]

C. J. Fennell and J. D. Gezelter. Is the Ewald summation still necessary? Pairwise alternatives to the accepted standard for long-range electrostatics. J. Chem. Phys., 2006. doi:http://dx.doi.org/10.1063/1.2206581.