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}\).
-
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>> ¢ers, 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
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
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 aretotal
,electronic
andnuclear
.- 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:
Compute the total potential as \(V(\mathbf{r}) = V_{vac}(\mathbf{r}) + V_R(\mathbf{r})\)
Compute the surface charge distribution \(\gamma_s(\mathbf{r})\) with computeGamma
Compute a new reaction potential \(V_R(\mathbf{r})\) by application of the Poisson operator with solvePoissonEquation
Calculate the update of the reaction potential as \(\Delta V_R(\mathbf{r}) = V_R(\mathbf{r}) - V_R^{old}(\mathbf{r})\)
Accelerate convergence of the reaction potential through KAIN
Update the reaction potential as \(V_R(\mathbf{r}) = V_R^{old}(\mathbf{r}) + \Delta V_R(\mathbf{r})\)
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:
Set the total potential as \(V(\mathbf{r}) = V_{vac}(\mathbf{r})\)
Compute the surface charge distribution \(\gamma_s(\mathbf{r})\) from this potential
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
andnuclear
.
-
double setConvergenceThreshold(double prec)¶
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\)
-
virtual void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) override¶
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
-
virtual void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) override¶