CosmoRSParameters

class CosmoRSParameters(alpha_prime=None, a_effective=None, f_correlation=None, lambda_ortho=None, c_hydrogen_bonding=None, sigma_hydrogen_bonding=None, hydrogen_bond_temperature=None, hydrogen_bond_donors=None, hydrogen_bond_acceptors=None, r_average=None, lambda_comb=None, lambda_reg=None, omega_ring=None, eta_gas=None, dispersion=None, pka_scale=None, pka_offset=None, sigma_potential_tolerance=None, sigma_potential_max_steps=None, bin_size=None)

Constructor for COSMO-RS parameters class. These parameters are used for calculating chemical potentials of different states. The default parameters are the generally recommended ones.

Parameters:

• alpha_prime (PhysicalQuantity of type energy per area) – Weak/misfit interaction coefficient.

• a_effective (PhysicalQuantity of type area) – The thermodynamically independent contact area.

• f_correlation (float) – Amount of correlation charge density to include in the misfit energy.

• lambda_ortho (float) – Amount of local charge density that contributes to the orthogonal charge density.

• c_hydrogen_bonding (PhysicalQuantity of type energy per area) – Hydrogen bonding coefficient for strong electrostatic interactions.

• sigma_hydrogen_bonding (PhysicalQuantity of type charge per area) – The charge density threshold above which the interaction is considered as including hydrogen bonding.

• hydrogen_bond_temperature (bool) – Whether or not the strength of hydrogen bonding is scaled according to the temperature.

• hydrogen_bond_donors (Sequence of PeriodicTableElement) – The elements that can act as donors for hydrogen bonding.

• hydrogen_bond_acceptors (Sequence of PeriodicTableElement) – The elements that can act as acceptors for hydrogen bonding.

• r_average (PhysicalQuantity of type length) – The sphere radius in which the surface charges are averaged.

• lambda_comb (float) – Factor used in calculating the combinatoral chemical potential.

• lambda_reg (float) – Factor for the energy difference between the exact and sigma-averaged screening charges.

• omega_ring (PhysicalQuantity of type energy) – Gas phase energy contribution for atoms in rings.

• eta_gas (float) – Scaling of the gas phase chemical potential with temperature.

• dispersion (dict) – The dispersion energy densities for elements.

• pka_scale (float) – The gradient of the pKa linear correction.

• pka_offset (float) – The offset of the pKa linear correction.

• sigma_potential_tolerance (PhysicalQuantity of type charge per area) – The convergence threshold for the sigma potential.

• sigma_potential_max_steps (int) – The maximum number of steps before the sigma potential is evaluated as non-converged.

• bin_size (Sequence of PhysicalQuantity of type charge per area) – The bin size for the sigma profile histogram.

aEffective()

Returns:

The thermodynamically independent contact area.

Return type:

PhysicalQuantity of type area

alphaPrime()

Returns:

Weak/misfit interaction coefficient.

Return type:

PhysicalQuantity of type energy per area

binSize()

Returns:

The histogram bin size of the charge densities.

Return type:

PhysicalQuantity of type charge per area

cHydrogenBonding()

Returns:

Hydrogen bonding coefficient for strong electrostatic interactions.

Return type:

PhysicalQuantity of type energy per area

createDatabaseEntry(file_path, name, reference=None, description=None)

Create a parameter database file.

Parameters:

• file_path (str) – The location of where the database entry file is to be saved.

• name (str) – The name of the parameter database entry.

• reference (str | None) – The optional reference for the source of the data.

• description (str | None) – Optional description of the database entry.

defaultDispersion()

Returns:

The average of the elemental dispersion parameters, used when the dispersion for the given element is undefined.

Return type:

PhysicalQuantity of type energy per area

dictionary()

Returns:

A dictionary of the parameters that defines the model.

Return type:

dict

difference(other, rel_tol=0.0001)

Return a dictionary of all of the parameters in the object that differ by the given factor from the parameters in another object.

dispersion()

Returns:

The dispersion energy densities for elements.

Return type:

dict

etaGas()

Returns:

Scaling of the gas phase chemical potential with temperature.

Return type:

float

fCorrelation()

Returns:

Amount of correlation charge density to include in the misfit energy.

Return type:

float

hydrogenBondAcceptors()

Returns:

The elements that can act as acceptors of hydrogen bonding

Return type:

Sequence of PeriodicTableElement

hydrogenBondDonors()

Returns:

The elements that can act as donors of hydrogen bonding

Return type:

Sequence of PeriodicTableElement

hydrogenBondTemperature()

Returns:

Whether or not hydrogen bonding is scaled with temperature.

Return type:

bool

lambdaComb()

Returns:

Factor used in calculating the combinatoral chemical potential.

Return type:

float

lambdaOrtho()

Returns:

Amount of local charge density that contributes to the orthogonal charge density.

Return type:

float

lambdaReg()

Returns:

Factor for the energy difference between the exact and averaged screening charges.

Return type:

float

nlinfo()

Returns:

The nlinfo.

Return type:

dict

omegaRing()

Returns:

Gas phase energy contribution for atoms in rings.

Return type:

PhysicalQuantity of type energy

pkaOffset()

Returns:

The intercept of the linear correction for the pKa.

Return type:

float

pkaScale()

Returns:

The gradient of the linear correction for the pKa

Return type:

float

rAverage()

Returns:

The sphere radius for which the surface charges are averaged.

Return type:

PhysicalQuantity of type length

sigmaHydrogenBonding()

Returns:

The charge density hydrogen bonding threshold.

Return type:

PhysicalQuantity of type charge per area

sigmaPotentialMaxSteps()

Returns:

The number of steps before the sigma potential is evaluated as non-converged.

Return type:

int

sigmaPotentialTolerance()

Returns:

The threshold under which the sigma potential is evaluated as converged.

Return type:

PhysicalQuantity of type charge per area

uniqueString()

Return a unique string representing the state of the object.

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

Usage Examples

Set some custom parameters for the COSMO-RS model and use them to calculate the vapor pressure of water at 300K.

# Load the water COSMO-RS species
database = CosmoRSSpeciesDatabase()
water = database.exportSpecies('water')

# Create the parameter set with custom parameters for the gas chemical potential
parameters = CosmoRSParameters(
    r_average=0.5*Ang,
    omega_ring=-0.21*kiloCaloriePerMol,
    eta_gas=-9.15,
    lambda_reg=0.8,
    lambda_ortho=0.816,
    dispersion={
        Hydrogen: -4.10*kiloCaloriePerMol/nm**2,
        Oxygen: -4.20*kiloCaloriePerMol/nm**2,
    }
)

# Calculate the vapor pressure using the different parameters
liquid = CosmoRealSolvent(300*Kelvin, water, parameters)
gas = CosmoRealGas(300*Kelvin, parameters)
mu_liquid = liquid.chemicalPotential(water)
mu_gas = gas.chemicalPotential(water)
vapor_pressure = numpy.exp((mu_liquid - mu_gas) / (300 * Kelvin * boltzmann_constant)) * bar

nlprint(f'The calculated vapor pressure is {vapor_pressure}')

cosmors_parameters_example.py

Notes

The CosmoRSParameters class is a container class for the parameters used in the COSMO-RS method. These parameters are used in the calculation of the chemical potential of molecules in different states. They are mostly empirical parameters that are fitted to large sets of experimental data. The equations used in COSMO-RS are published in a number of places[1][2][3][4]. Care must be taken when reading these sources however, as there are some mistakes and oversights in the original literature.

A number of the parameters are used for calculating either the liquid or solvent chemical potentials. When taking surface charges from a COSMO-DFT calculation the surface charge densities \(\sigma^*\) may fluctuate in an nonphysical manner. This is due to the linear fitting algorithm used to calculate the charge densities. The charge densities used in the COSMO-RS method are smoothed using a sigma averaging method. Here the smoothed charge densities \(\sigma\) are calculated as:

\[\sigma_i = \frac{\sum_j \sigma_j^* I_{ij}}{\sum_j I_{ij}}\]

\[I_{ij} = \frac{r^2_j}{r^2_j + r^2_{av}} \exp\left(-\frac{d^2_{ij}}{r^2_{av}}\right)\]

Here \(d_{ij}\) is the distance between two surface segments, \(r_j\) is the average radius of the surface segment and \(r_{av}\) is the averaging radius. The averaging radius is one of the parameters in the COSMO-RS model, and can be set using the argument r_average.

In addition to the sigma averaged charge densities \(\sigma\), COSMO-RS also uses orthogonal charges \(\sigma^\perp\). These are created by first creating a set of charge densities \(\sigma^\circ\) that averaged from the \(\sigma^*\) charge densities using double the averaging radius. The orthogonal charges are then defined as:

\[\sigma^\perp = \sigma^\circ - \lambda_{ortho}\sigma\]

Here \(\lambda_{ortho}\) is another COSMO-RS parameter. This can be set in the CosmoRSParameters object using the argument lambda_ortho.

With these two charge densities, the \(\sigma\)-potential \(\bar{\mu}(\sigma, \sigma^\perp)\) is then calculated by iteratively summing over all of the surface segments of each component of the solvent. The equation for this can be given as:

\[\bar{\mu}(\sigma, \sigma^\perp) = -\ln\left[\frac{1}{\sum_i x_i A_i} \sum_i x_i \sum_{\nu \in X_i} s_v \exp\left(\bar{\mu}(\sigma_\nu, \sigma^\perp_\nu) - \frac{a_{eff}}{RT}E(\sigma, \sigma^\perp, \sigma_\nu, \sigma_\nu^\perp) \right) \right]\]

Here \(x_i\) is the mole fraction of each solvent component, \(A_i\) is the total molecule surface area, \(E(\sigma, \sigma^\perp, \sigma_\nu, \sigma_\nu^\perp)\) is the surface segment interaction energy. The effective contact area \(a_{eff}\) is another free parameter in the COSMO-RS method. This can be set using the a_effective argument. Note also that since the \(\sigma\)-potential \(\bar{\mu}\) appears on both the right and left side of the equation, the sigma potential for each surface segment is calculated iteratively.

The surface segment interaction energy is the sum of two additional terms, the misfit energy and the hydrogen bonding energy. The misfit energy, \(E_{misfit}\) which measures the energy of interaction between two surface segments can be given as:

\[E_{misfit}(\sigma_i, \sigma_j, \sigma_i^\perp, \sigma_j^\perp) = \frac{\alpha'}{2}(\sigma_i + \sigma_j)[(\sigma_i + \sigma_j) + f_{corr}(\sigma_i^\perp + \sigma_j^\perp)]\]

Here the charge densities are taken from two different surface elements \(i\) and \(j\). This equation also contains two additional free parameters. \(\alpha'\) represents the strength of the misfit energy interaction, and can be set with the input argument alpha_prime. Likewise \(f_{corr}\) represents the strength of the correlation charge density interaction, and can be set with the argument f_correlation.

The hydrogen bonding energy \(E_{hb}\) utilizes only the local charge densities \(\sigma\). This energy can be given as:

\[E_{hb}(\sigma_i, \sigma_j) = c_{hb}f_{hb}(T) \max(0, \max[\sigma_i, \sigma_j] - \sigma_{hb}) \min(0, \min[\sigma_i, \sigma_j] + \sigma_{hb})\]

In this equation \(c_{hb}\) represents the strength of the hydrogen bonding interaction, while \(\sigma_{hb}\) represents the charge density limit above which interactions are deemed to include hydrogen bonding. These are both free parameters, and can be set with the c_hydrogen_bonding and sigma_hydrogen_bonding respectively. Scaling of the hydrogen bonding interaction can also be included through the function \(f_{hb}(T)\). This can be given as:

\[f_{hb}(T) = \frac{T \ln(1+\exp(\bar{E}/RT) / 200)}{\bar{T} \ln(1+\exp(\bar{E}/R\bar{T}) / 200)}\]

Here \(\bar{E}\) is a reference energy, set to 20kJ/mol and \(\bar{T}\) is a reference temperature, set to 298.15K. This temperature scaling lowers the hydrogen bonding energy at temperatures above 298.15K. As this scaling has no free parameters, it can only be either enabled or disabled. This is done with the argument hydrogen_bond_temperature.

Since the \(\sigma\)-potential is calculated by summing over specific surface elements, it is possible to also use the attached atom of each surface segment to determine whether or not to include hydrogen bonding. This can be useful in cases, such as when modeling metal cations, where molecules may have large charge densities from atoms that do not have any hydrogen bonding character. The elements to include as hydrogen bonding donors or acceptors can be specified with the arguments hydrogen_bond_donors and hydrogen_bond_acceptors, respectively. These both take a list of elements or the flag All. The flag All indicates that hydrogen bonding is to be included for all surface segments withing the hydrogen bonding range of charge density regardless of attached atom.

Once the \(\sigma\)-potential is calculated, the residual chemical potential \(\mu^{res}_X\) of a molecule \(X\) can simply be calculated by summing the interactions of the solute surface segments with those of the solvent. This is done through the equation:

\[\mu^{res}_X = \frac{RT}{a_{eff}} \sum_{i \in X} s_\nu \bar{\mu}(\sigma_i, \sigma_i^\perp)\]

The total chemical potential of a solvated molecule is the sum of both the residual and combinatorial chemical potential. This latter chemical potential includes cavitation and entropic terms. The combinatatorial chemical potential \(\mu^{comb}_X\) can be given as:

\[\mu^{comb}_X = \lambda_{comb} RT\left[\ln(V^X) + \left(1 - \frac{V^X}{V_{av}} - \ln(V^{av}) \right) + \left(1 - \frac{A^X}{A_{av}} - \ln(A^{av}) \right) \right]\]

Here \(A^X\) and \(V^X\) are the surface area and volume of the solute species. Likewise \(A^{av}\) and \(V^{av}\) are the average surface area and volume of the components of the solvent. All areas and volumes are also divided by the reference of \(1\AA^2\) and \(1\AA^3\), respectively, to render them dimensionless. The free parameter \(\lambda_{comb}\) sets the relative size of the combinatorial chemical potential. This can be specified with the argument lambda_comb.

A number of free parameters are also used in the estimation of the gas phase chemical potential. In COSMO-RS the gas chemical potential of a substance is defined so that the vapor pressure \(p\) can be given as:

\[p = p_0 \exp([\mu^{liq} - \mu^{gas}] / RT)\]

Here \(p_0\) is a reference pressure, defined to be 1 bar. The overall equation of the gas phase chemical potential of a molecule \(\mu^{gas}_X\) can be given as:

\[\mu^{gas}_X = \Delta^X - \sum_k \gamma_k A_k^X -\omega n^X_{ring} + \eta RT\]

This expression contains a number of terms. The first, \(\Delta^X\) represents the difference in solvation energy between the gas and solvent states. This can be given as:

\[\Delta^X = E^{gas} - E^{COSMO} + \frac{\lambda_{reg}}{2} \sum_i \Phi_i s_\nu (\sigma_i - \sigma_i^*)\]

Here \(E^{gas}\) and \(E^{COSMO}\) are the total energy of the molecule in the gas and COSMO solvent model, respectively. The COSMO energy includes both the dielectric energy and the polarization of the molecule resulting from the charge of the solvent cavity. The final sum is a consequence of the COSMO dielectric energy being calculated with the raw charges. In COSMO-RS smoothed charges are always used, and so the final sum accounts for the energy associated with smoothing the charge density distribution. Here \(\Phi\) is the molecular electrostatic potential at each surface element. The free parameter \(\lambda_{reg}\) estimates how much the total COSMO energy is raised by the molecule being polarized by the solvent cavity surface. This can be set with the argument lambda_reg.

Going back to the equation for the total gas chemical potential, the second term accounts for the loss of dispersion interactions from the solvent in transferring the molecule to the gas phase. The dispersion energies are simply calculated by scaling the total surface area associated with each element, \(A_k^X\) with a scaling factor \(\gamma_k\). These scaling factors are different for each element. Since the dispersion energy depends only on the solvent cavity, it is the same for any solvent. For that reason it is not included in the solvent chemical potential, but rather in the gas chemical potential. The scaling factors for each element can be set with the argument dispersion. This takes a dictionary, where the keys are each element to be set and the value is the scaling factor.

The third term in the gas chemical potential is a special correction for molecules with rings. Here \(n_{ring}\) represents the number of atoms in cyclic structures and \(\eta\) is an energy scaling factor. This can be set with the argument omega_ring. The final term in the gas chemical potential accounts for the temperature of the gas. Here \(\eta\) is a dimensionless scaling factor. This can be set with the argument eta_gas.

When calculating pK_a values, a linear correction can be applied to the value calculated by COSMO-RS to give an estimate of the actual pK_a value[5]. This corrects for both systematic deficiencies in the underlying DFT calculation and the COSMO-RS methodology. Using this linear correction the pK_a can be given as:

\[pK_a = A \times pK_a^{raw} + B\]

Here \(pK_a^{raw}\) is the pK_a value directly calculated by COSMO-RS. The factors \(A\) and \(B\) can be set with the arguments pka_scale and pka_offset, respectively.

The final parameters in the CosmoRSParameters object relate to numerical details in how the COSMO-RS calculation is performed. Since the \(\sigma\)-potential is calculated iteratively, this is done until either the \(\sigma\)-potential reaches a set convergence limit or a maximum number of steps has been reached. The convergence limit can be set with the argument sigma_potential_tolerance. Likewise the maximum number of steps can be set with the argument sigma_potential_max_steps. In some cases, such as at low temperature, the \(\sigma\)-potential may not converge using the default settings. In these cases either the convergence limit can be raised or the number of steps increased until convergence is reached. When calculating plots of either the \(\sigma\)-profile or the \(\sigma\)-potential the charge density is divided into a number of steps or bins. The argument bin_size sets the spacing of these steps. Since the COSMO-RS calculations are done by summing over surface segments, this argument is not used in calculating chemical potentials.

When creating a CosmoRSParameters object, any number or combination of parameters may be specified. The default values for all parameters is taken from the quantumatk_u2022-12 parameter set. The values contained in this parameter set can be shown using the COSMO-RS Parameters dialog in the COSMO-RS Analyzer.

Once a parameter set is created, it can be saved for further use in QuantumATK. Parameter sets can be added to QuantumATK through the CosmoRSParameterSets class, which interfaces to a local database of parameter sets. Parameter sets can also be added manually. The method createDatabaseEntry creates a representation of the parameter set in a file. The parameter set is then added to the database by placing it in the .quantumatk/databases/cosmors_parameters/ in your home directory.

[1]

Andreas Klamt, Volker Jonas, Thorsten Bürger, and John C. W. Lohrenz. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A., 102:5074–5085, 1998. doi:10.1021/jp980017s.

[2]

Andreas Klamt. The basic COSMO-RS, pages 83–107. Elsevier, 2005.

[3]

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

[4]

Cory C. Pye, Tom Ziegler, Erik van Lenthe, and Jaap N. Louwen. An implementation of the conductor-like screening model of solvation withing the Amsterdam density functional package - Part II. COSMO for real solvents. Can. J. Chem., 87:790–797, 2009. doi:10.1139/V09-008.

[5]

Andreas Klamt, Frank Eckert, Michael Diedenhofen, and Michael E. Beck. First principles calculations of aqueous pKa values for organic and inorganic acids using COSMO-RS reveal an inconsistency in the slope of the pKa scale. J. Phys. Chem. A., 107:9380–9386, 2003. doi:10.1021/jp034688o.