Environment

Classes for the solvent environment overlay

Cavity

class Cavity : public mrcpp::RepresentableFunction<3>

Interlocking spheres cavity centered on the nuclei of the molecule.

The Cavity class represents the following function Fosso-Tande2013

\[\begin{split} C(\mathbf{r}) = 1 - \prod^N_{i=1} (1-C_i(\mathbf{r})) \\ C_i(\mathbf{r}) = 1 - \frac{1}{2}\left( 1 + \textrm{erf}\left(\frac{|\mathbf{r} - \mathbf{r}_i| - R_i}{\sigma_i}\right) \right) \end{split}\]

where \(\mathbf{r}\) is the coordinate of a point in 3D space, \(\mathbf{r}_i\) is the coordinate of the i-th nucleus, \(R_i\) is the radius of the i-th sphere, and \(\sigma_i\) is the width of the transition between the inside and outside of the cavity. The transition has a sigmoidal shape, such that the boundary is a smooth function instead of sharp boundaries often seen in other continuum models. This function is \(1\) inside and \(0\) outside the cavity.

The radii are computed as:

\[ R_{i} = \alpha_{i} R_{0,i} + \beta_{i}\sigma_{i} \]

where:

  • \(R_{0,i}\) is the atomic radius. By default, the van der Waals radius.

  • \(\alpha_{i}\) is a scaling factor. By default, 1.1

  • \(\beta_{i}\) is a width scaling factor. By default, 0.5

  • \(\sigma_{i}\) is the width. By default, 0.2 bohr

Public Functions

Cavity(const std::vector<mrcpp::Coord<3>> &coords, const std::vector<double> &R, const std::vector<double> &alphas, const std::vector<double> &betas, const std::vector<double> &sigmas)

Initializes the members of the class and constructs the analytical gradient vector of the Cavity.

inline Cavity(const std::vector<mrcpp::Coord<3>> &coords, const std::vector<double> &R, double sigma)

Initializes the members of the class and constructs the analytical gradient vector of the Cavity.

This CTOR applies a single width factor to the cavity and does not modify the radii. That is, in the formula:

\[ R_{i} = \alpha_{i} R_{0,i} + \beta_{i}\sigma_{i} \]

for every atom \(i\), \(\alpha_{i} = 1.0\) and \(\beta_{i} = 0.0\).

double evalf(const mrcpp::Coord<3> &r) const override

Evaluates the value of the cavity at a 3D point \(\mathbf{r}\).

Parameters:

r – coordinate of 3D point at which the Cavity is to be evaluated at.

Returns:

double value of the Cavity at point \(\mathbf{r}\)

inline std::vector<mrcpp::Coord<3>> getCoordinates() const

Returns centers.

inline std::vector<double> getOriginalRadii() const

Returns radii_0.

inline std::vector<double> getRadii() const

Returns radii.

inline std::vector<double> getRadiiScalings() const

Returns alphas.

inline std::vector<double> getWidths() const

Returns sigmas.

inline std::vector<double> getWidthScalings() const

Returns betas.

void printParameters() const

Print parameters.

Protected Attributes

std::vector<double> radii_0

Contains the unscaled radius of each sphere in #Center.

std::vector<double> alphas

The radius scaling factor for each sphere.

std::vector<double> betas

The width scaling factor for each sphere.

std::vector<double> sigmas

The width for each sphere.

std::vector<double> radii

Contains the radius of each sphere in #Center. \(R_i = \alpha_{i} R_{0,i} + \beta_{i}\sigma_{i}\).

std::vector<mrcpp::Coord<3>> centers

Contains each of the spheres centered on the nuclei of the Molecule.

auto gradCavity(const mrcpp::Coord<3> &r, int index, const std::vector<mrcpp::Coord<3>> &centers, const std::vector<double> &radii, const std::vector<double> &widths) -> double

Constructs a single element of the gradient of the Cavity.

This constructs the analytical partial derivative of the Cavity \(C\) with respect to \(x\), \(y\) or \(z\) coordinates and evaluates it at a point \(\mathbf{r}\). This is given for \(x\) by

\[ \frac{\partial C\left(\mathbf{r}\right)}{\partial x} = \left(1 - C{\left(\mathbf{r} \right)}\right) \sum_{i=1}^{N} - \frac{\left(x-{x}_{i}\right)e^{- \frac{\operatorname{s_{i}}^{2}{\left(\mathbf{r} \right)}}{\sigma^{2}}}} {\sqrt{\pi}\sigma\left(0.5 \operatorname{erf}{\left(\frac{\operatorname{s_{i}}{\left(\mathbf{r} \right)}}{\sigma} \right)} + 0.5\right) \left| \mathbf{r} - \mathbf{r}_{i} \right|} \]
where the subscript \(i\) is the index related to each sphere in the cavity, and \(\operatorname{s}\) is the signed normal distance from the surface of each sphere.

Parameters:
  • r – The coordinates of a test point in 3D space.

  • index – An integer that defines the variable of differentiation (0->x, 1->y and 2->z).

  • centers – A vector containing the coordinates of the centers of the spheres in the cavity.

  • radii – A vector containing the radii of the spheres.

  • width – A double value describing the width of the transition at the boundary of the spheres.

Returns:

A double number which represents the value of the differential (w.r.t. x, y or z) at point r.

Permittivity

class Permittivity : public mrchem::StepFunction

Permittivity function related to a substrate molecule and a solvent continuum. The Permittivity class represents the following function Fosso-Tande2013.

\[ \epsilon(\mathbf{r}) = \epsilon_{in}\exp\left(\left(\log\frac{\epsilon_{out}}{\epsilon_{in}} \right) \left(1 - C(\mathbf{r})\right)\right) \]
where \(\mathbf{r}\) is the coordinate of a point in 3D space, \( C \) is the #cavity function of the substrate, and \(\epsilon_{in}\) and \( \epsilon_{out} \) are the dielectric constants describing, respectively, the permittivity inside and outside the #cavity of the substrate.

Public Functions

Permittivity(std::shared_ptr<Cavity> cavity, double epsilon_in, double epsilon_out, std::string formulation)

Standard constructor. Initializes the #cavity, #epsilon_in and #epsilon_out with the input parameters.

Parameters:
  • cavity – interlocking spheres of Cavity class.

  • epsilon_in – permittivity inside the #cavity.

  • epsilon_out – permittivity outside the #cavity.

  • formulation – Decides which formulation of the Permittivity function to implement, only exponential available as of now.

double evalf(const mrcpp::Coord<3> &r) const override

Evaluates Permittivity at a point in 3D space with respect to the state of #inverse.

Parameters:

r – coordinates of a 3D point in space.

Returns:

\(\frac{1}{\epsilon(\mathbf{r})}\) if #inverse is true, and \( \epsilon(\mathbf{r})\) if #inverse is false.

DHScreening

class DHScreening : public mrchem::StepFunction

Square of the Debye-Huckel Screening parameter.

This is used for the Poisson-Boltzmann solver. The DHScreening function is defined as

\[\begin{split} \kappa^2(\mathbf{r}) = \begin{cases} \kappa^2_{out} & \text{if } \mathbf{r} \notin \Omega_{ion} \\ 0.0 & \text{if } \mathbf{r} \in \Omega_{ion} \end{cases} \end{split}\]
This can be parametrized a number of ways. The one used here is
\[ \kappa^2(\mathbf{r}) = (1 - C_{ion}(\mathbf{r})) \kappa^2_{out} \]
Where \(C_{ion}(\mathbf{r})\) is the ion accessible Cavity function.

Public Functions

DHScreening(std::shared_ptr<Cavity> cavity_ion, double kappa_out, const std::string &formulation)

Standard constructor. Initializes the #cavity_ion and #kappa_out with the input parameters.

#kappa_out is given by

\[ \kappa = \sqrt{\frac{2000 I_{0} e^2 N_a I_0}{\epsilon_{out} \epsilon_{in} k_B T}} \]
where \(N_a\) is the Avogadro constant, e is the elementary charge, \(I_0\) is the concentration of the ions, \(k_B\) is the Boltzmann constant, \(T\) is the temperature, \(\epsilon_{out}\) is the permittivity of the solvent and \(\epsilon_{in}\) is the permittivity of free space.

Parameters:
  • cavity_ion – interlocking spheres of Cavity class.

  • kappa_out – value of the screening function outside the #cavity_ion.

  • formulation – Decides which formulation of the DHScreening function to implement, only continuous screening function available.

double evalf(const mrcpp::Coord<3> &r) const override

Evaluates DHScreening at a point in 3D space.

Parameters:

r – coordinates of a 3D point in space.

Returns:

Value at point r.

Private Members

std::string formulation = {"Continuous Screening Function"}

Formulation of the DHScreening function. Only linear variable is used now.

GPESolver

class GPESolver

Solves the Generalized Poisson equation iteratively.

The Generalized Poisson equation is given by

\[ \nabla \cdot \left( \epsilon(\mathbf{r}) \nabla V(\mathbf{r}) \right) = -4\pi \rho(\mathbf{r}) \]
where \(\epsilon(\mathbf{r})\) is the permittivity, \(V(\mathbf{r})\) is the total electrostatic potential and \(\rho(\mathbf{r})\) is the molecular charge density defined as:
\[ \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r}) \]
where \(\rho_{el}\) is the electronic charge density and \(\rho_{nuc}\) is the nuclear charge density. The Generalized Poisson equation is solved iteratively through a set of micro-iteration on each SCF-iteration by appliation of the Poisson operator \(\mathcal{P}\) :cite:Fosso-Tande2013
\[ V_R(\mathbf{r}) = \mathcal{P} \star \left[ \rho_{eff}(\mathbf{r}) - \rho(\mathbf{r}) + \gamma_s(\mathbf{r}) \right] \]
where \(\gamma_s(\mathbf{r})\) is the surface charge distribution describing the polarization at the surface, \(\rho_{eff}(\mathbf{r})\) is the effective charge density given by \(\frac{\rho(\mathbf{r})}{\epsilon(\mathbf{r})}\) and \(V_R(\mathbf{r})\) is the reaction potential.

We utilize a so-called dynamic threshold to more easily converge the reaction potential. This is done by setting the convergence threshold of the micro-iterations to the MO update of the previous SCF iteration, unless the MO update is small enough (once the quality of the MOs is good enough, we use the default convergence threshold). Another optimization used is that we utilize the previous SCF converged Reaction potential as an initial guess for the next micro-iterations. These procedures are investigated and explained in :cite:gerez2023

Subclassed by mrchem::PBESolver

Public Functions

double setConvergenceThreshold(double prec)

Sets the convergence threshold for the micro-iterations, used with dynamic thresholding.

will check if the MO update is small enough (ten times as big) wrt. to the scf convergence threshold, if so, it will use the default convergence threshold. If not, it will use the MO update as the convergence threshold.

Parameters:

prec – value to set the convergence threshold to

Returns:

the current convergence threshold.

auto computeEnergies(const Density &rho_el) -> std::tuple<double, double>

Computes the energy contributions from the reaction potential.

We compute the reaction energy through the following integral:

\[ E_{R} = \frac{1}{2}\int \rho_{el}(\mathbf{r}) V_R(\mathbf{r}) d\mathbf{r} + \frac{1}{2} \int \rho_{nuc}(\mathbf{r}) V_R(\mathbf{r}) d\mathbf{r} \]
Each term represents the electronic and nuclear contributions to the reaction energy, respectively. We compute each term separately, and return a tuple containing both.

Parameters:

rho_el – the electronic charge density

Returns:

a tuple containing the electronic and nuclear energy contributions

Protected Functions

void computeDensities(const Density &rho_el, Density &rho_out)

computes density wrt. the density_type variable

The total charge density is given by the sum of the electronic and nuclear charge densities:

\[ \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r}) \]
where \(\rho_{el}\) is the electronic charge density and \(\rho_{nuc}\) is the nuclear charge density. The nuclear charge density is stored in the class variable rho_nuc, while we compute the electronic charge density from the molecular orbitals. The class variable density_type decides the density which will be computed in rho_out, options are total, electronic and nuclear.

Parameters:
  • Phi – the molecular orbitals

  • rho_out – Density function in which the density will be computed.

virtual void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma)

Computes the surface charge distibution due to polarization at the solute-solvent boundary.

The surface charge distribution is given by

\[ \gamma_s(\mathbf{r}) = \frac{\log \frac{\epsilon_{in}}{\epsilon_{out}}}{4 \pi } \left( \nabla C(\mathbf{r}) \cdot \nabla V(\mathbf{r})\right) \]
where \(\epsilon_{in}\) is the permittivity inside the cavity and \(\epsilon_{out}\) is the permittivity outside the cavity.

Parameters:
  • potential – Potential used to compute \(\nabla V(\mathbf{r})\)

  • out_gamma – ComplexFunction in which the surface charge distribution will be computed.

mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma, const Density &rho_el)

Iterates once through the Generalized Poisson equation to compute the reaction potential.

Constructs the effective charge density \(\rho_{eff}(\mathbf{r})\) and the Poisson operator \(\mathcal{P}\) as:

\[ V_R(\mathbf{r}) = \mathcal{P} \star \left[ \rho_{eff}(\mathbf{r}) - \rho(\mathbf{r}) + \gamma_s(\mathbf{r}) \right] \]
where \(\gamma_s(\mathbf{r})\) is the surface charge distribution describing the polarization at the surface, \(\rho_{eff}(\mathbf{r})\) is the effective charge density given by \(\frac{\rho(\mathbf{r})}{\epsilon(\mathbf{r})}\) and \(V_R(\mathbf{r})\) is the reaction potential.

Parameters:
  • ingamma – the surface charge distribution

  • Phi – the molecular orbitals

Returns:

the reaction potential

void accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain)

Uses KAIN to accelerate convergece of the reaction potential.

Parameters:
  • dfunc – the current update of the reaction potential

  • func – the current reaction potential

  • kain – the KAIN object

void runMicroIterations(const mrcpp::ComplexFunction &V_vac, const Density &rho_el)

Iterates through the application of the Poisson operator to Solve the Generalized Poisson equation.

Iterating through the application of the Poisson operator is done through a set of micro-iterations, where the convergence threshold is set to the MO update of the previous SCF iteration. The micro-iterations are done through the following steps:

  1. Compute the total potential as \(V(\mathbf{r}) = V_{vac}(\mathbf{r}) + V_R(\mathbf{r})\)

  2. Compute the surface charge distribution \(\gamma_s(\mathbf{r})\) with computeGamma

  3. Compute a new reaction potential \(V_R(\mathbf{r})\) by application of the Poisson operator with solvePoissonEquation

  4. Calculate the update of the reaction potential as \(\Delta V_R(\mathbf{r}) = V_R(\mathbf{r}) - V_R^{old}(\mathbf{r})\)

  5. Accelerate convergence of the reaction potential through KAIN

  6. Update the reaction potential as \(V_R(\mathbf{r}) = V_R^{old}(\mathbf{r}) + \Delta V_R(\mathbf{r})\)

  7. Check if the reaction potential has converged, if not, repeat from step 1.

Parameters:
  • V_vac – the vacuum potential

  • Phi_p – the molecular orbitals

mrcpp::ComplexFunction &solveEquation(double prec, const Density &rho_el)

Setups and computes the reaction potential through the microiterations.

An initial guess of the reaction potential is computed with the following steps:

  1. Set the total potential as \(V(\mathbf{r}) = V_{vac}(\mathbf{r})\)

  2. Compute the surface charge distribution \(\gamma_s(\mathbf{r})\) from this potential

  3. Iterate once through the application of the Poisson operator to return the initial guess of the reaction potential \(V_R(\mathbf{r})\)

the method then runs the micro-iterations through runMicroIterations and returns the converged reaction potential. If this is not the first SCF iteration, the previous converged reaction potential is used as an initial guess for the micro-iterations.

Parameters:
  • V_vac – the vacuum potential

  • Phi_p – the molecular orbitals

Returns:

The converged reaction potential for the current SCF iteration

void resetComplexFunction(mrcpp::ComplexFunction &function)

Frees the memory used by the FunctionTrees of the input Complexfunction and reallocates them.

This is done to avoid memory leaks when the ComplexFunction is used in the micro-iterations.

Parameters:

function – the ComplexFunction to reset

Protected Attributes

SCRFDensityType density_type

Decides which density we will use for computing the reaction potential, options are total, electronic and nuclear.

PBESolver

class PBESolver : public mrchem::GPESolver

Solves the Poisson-Boltzmann equation iteratively.

The Poisson-Boltzmann equation is solved iteratively using the SCRF procedure outlined in GPESolver. The Poisson-Boltzmann equation models the electrostatic potential in a solvent with electrolytes. The general equation for electrolyte solutions is given by

\[ \nabla \cdot \epsilon \nabla V_{tot} = -4\pi\left(\rho_{el} + \rho_{nuc} + \rho_{ext}\right) \]
where \( V_{tot} \) is the total electrostatic potential, \( \epsilon\) is the permittivity function of the solvent, \(\rho_{el}\) is the electronic charge density, \(\rho_{nuc}\) is the nuclear charge density and \(\rho_{ext}\) is the external charge density. In the general form for the Poisson-Boltzmann equation, the external charge density is approximated by assuming a boltzmann distribution of the ions.
\[ \rho_{ext} = \sum_{i}^{N_{ion}} q_i e I_{0, i}\exp\left(-\frac{q_i e V_{tot}}{k_B T}\right) \]
where \(I_{0, i}\) is the concentration of the i-th ion species, \(q_i\) is the charge of the i-th ion species, \(k_B\) is the Boltzmann constant and \(T\) is the temperature. In this implementation we assume a 1:1 ( \(I_{0, 0} = I_{0, 1}\)) electrolyte soluttion of ions of same opposite charges ( \(z_i = +1, -1\)). This simplifies the external density to
\[ \rho_{ext} = -2 e I_{0} \sinh\left(\frac{e V_{tot}}{2 k_B T}\right) \]
where \(I_{0}\) is the concentration of the ions. We can plug this into the first equation (and massage terms a bit) to arrive at the Poisson-Boltzmann equation for 1:1 electrolyte solution
\[ \nabla^2 V_{R} = -4\pi\frac{1-\epsilon}{\epsilon}\left(\rho_{el} + \rho_{nuc}\right) + \gamma_s - \kappa^2 \sinh\left(V_{tot}\right) \]
where \(\gamma_s\) is the surface charge density, \(\kappa\) is obtained from the DHScreening class and \(V_{R}\) is the reaction potential.

Subclassed by mrchem::LPBESolver

Protected Functions

virtual void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) override

constructs the surface chage distribution and adds it to the PB term

Method follows the implementation in GPESolver::computeGamma, but adds the PB term to the surface charge distribution.

Parameters:
  • potential[in] the potential to compute \(\nabla V\) from

  • out_gamma[out] the ComplexFunction in which to store the result

virtual void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term)

Computes the PB term.

The PB term is computed as \( \kappa^2 \sinh(V_{tot}) \) and returned.

Parameters:
  • V_tot[in] the total potential

  • salt_factor[in] the salt factor deciding how much of the total concentration to include in the PB term

  • pb_term[out] the ComplexFunction in which to store the result

Protected Attributes

DHScreening kappa

the DHScreening object used to compute the PB term \(\kappa\)

LPBESolver

class LPBESolver : public mrchem::PBESolver

Solves the Linearized Poisson-Boltzmann equation iteratively.

The Linearized Poisson-Boltzmann equation is solved iteratively using the SCRF procedure outlined in GPESolver and PBESolver. The linearized Poisson-Boltzmann equation is a further simplification of the Poisson-Boltzmann equation, outlined in PBESolver, where the PB term is expanded and only the linear term is included. This is a good approximation for low ionic strength solutions. The linearized Poisson-Boltzmann equation is given by

\[ \nabla^2 V_{R} = -4\pi\frac{1-\epsilon}{\epsilon}\left(\rho_{el} + \rho_{nuc}\right) + \gamma_s - \kappa^2 V_{tot} \]
where \(\gamma_s\) is the surface charge density, \(\kappa\) is obtained from the DHScreening class and \(V_{R}\) is the reaction potential.

Protected Functions

virtual void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) override

Computes the PB term.

The PB term is computed as \( \kappa^2 V_{tot} \) and returned.

Parameters:
  • V_tot[in] the total potential

  • salt_factor[in] the salt factor deciding how much of the total concentration to include in the PB term

  • pb_term[out] the ComplexFunction in which to store the result