Add documentation to the solvent classes

Author: Gabrielgerez

doc/index.rst

@@ -42,6 +42,8 @@ Features in MRChem-1.1:
         - Electric field
     + Solvent effects
         - Cavity-free PCM
+        - Poisson-Boltzmann PCM
+        - Linearized Poisson-Boltzmann PCM
 * Properties:
     + Ground state energy
     + Dipole moment

doc/programmers/code_reference/environment.rst

@@ -21,11 +21,38 @@ Permittivity
    :protected-members:
    :private-members: 
 
-SCRF
+DHScreening
 ------------
 
-.. doxygenclass:: mrchem::SCRF
+.. doxygenclass:: mrchem::DHScreening
    :project: MRChem
    :members:  
    :protected-members:
    :private-members: 
+
+GPESolver
+------------
+
+.. doxygenclass:: mrchem::GPESolver
+   :project: MRChem
+   :members:  
+   :protected-members:
+   :private-members: 
+
+PBESolver
+------------
+
+.. doxygenclass:: mrchem::PBESolver
+   :project: MRChem
+   :members:  
+   :protected-members:
+   :private-members: 
+
+LPBESolver
+------------
+
+.. doxygenclass:: mrchem::LPBESolver
+   :project: MRChem
+   :members:  
+   :protected-members:
+   :private-members: 

doc/programming.bib

@@ -333,7 +333,6 @@ @article{Fosso-Tande2013
 abstract = {We describe and present results of the implementation of the surface and volume polarization for electrostatics (SVPE) solvation model. Unlike most other implementations of the solvation model where the solute and the solvent are described with multiple numerical representations, our implementation uses a multiresolution, adaptive multiwavelet basis to describe both the solute and the solvent. This requires reformulation to use integral equations throughout as well as a conscious management of numerical properties of the basis. {\textcopyright} 2013 Elsevier B.V. All rights reserved.},
 author = {Fosso-Tande, Jacob and Harrison, Robert J.},
 doi = {10.1016/j.cplett.2013.01.065},
-file = {:home/ggerez/.local/share/data/Mendeley Ltd./Mendeley Desktop/Downloaded/Fosso-Tande, Harrison - 2013 - Implicit solvation models in a multiresolution multiwavelet basis.pdf:pdf},
 issn = {00092614},
 journal = {Chem. Phys. Lett.},
 pages = {179--184},
@@ -343,3 +342,20 @@ @article{Fosso-Tande2013
 year = {2013}
 }
 
+@article{gerez2023,
+author = {Gerez S, Gabriel A. and Di Remigio Eikås, Roberto and Jensen, Stig Rune and Bjørgve, Magnar and Frediani, Luca},
+title = {Cavity-Free Continuum Solvation: Implementation and Parametrization in a Multiwavelet Framework},
+journal = {Journal of Chemical Theory and Computation},
+volume = {19},
+number = {7},
+pages = {1986-1997},
+year = {2023},
+doi = {10.1021/acs.jctc.2c01098},
+    note ={PMID: 36933225},
+URL = { 
+        https://doi.org/10.1021/acs.jctc.2c01098
+},
+eprint = { 
+        https://doi.org/10.1021/acs.jctc.2c01098
+}
+}

doc/users/user_inp.rst

@@ -194,7 +194,7 @@ as the ``world_origin`` is the true origin).
 WaveFunction
 ------------
 
-Here we give the wavefunction method and whether we run spin restricted (alpha
+Here we give the wavefunction method, environment used (for solvent models) and whether we run spin restricted (alpha
 and beta spins are forced to occupy the same spatial orbitals) or not (method
 must be specified, otherwise defaults are shown):
 
@@ -203,6 +203,7 @@ must be specified, otherwise defaults are shown):
     WaveFunction {
       method = <wavefunction_method>        # Core, Hartree, HF or DFT
       restricted = true                     # Spin restricted/unrestricted
+      environment = pcm                     # Environment (pcm, pcm-pb, pcm-lpb) defaults to none
     }
 
 There are currently four methods available: Core Hamiltonian, Hartree,
@@ -212,6 +213,11 @@ B3LYP``), *or* you can set ``method = DFT`` and specify a "non-standard"
 functional in the separate DFT section (see below). See
 :ref:`User input reference` for a list of available default functionals.
 
+The solvent model implemented is a cavity free PCM, described in :cite:`gerez2023`. 
+In this model we have implemented the Generalized Poisson equation solver, keyword ``pcm``, a 
+Poisson-Boltzmann solver, keyword ``pcm-pb`` and a Linearized Poisson-Boltzmann solver, keyword ``pcm-lpb``. 
+Further details for the calculation have to be included in the ``PCM`` section, see :ref: `User input reference` for details.
+
 .. note::
 
     Restricted open-shell wavefunctions are not supported.

src/environment/Cavity.cpp

@@ -37,8 +37,7 @@ namespace mrchem {
 namespace detail {
 /**  @relates mrchem::Cavity
  *   @brief Constructs a single element of the gradient of the Cavity.
- *
- *   This constructs the analytical partial derivative of the Cavity \f$C\f$ with respect to \f$x\f$, \f$y\f$ or \f$z\f$
+ *   @details This constructs the analytical partial derivative of the Cavity \f$C\f$ with respect to \f$x\f$, \f$y\f$ or \f$z\f$
  * coordinates and evaluates it at a point \f$\mathbf{r}\f$. This is given for \f$x\f$ by
  * \f[
  *    \frac{\partial C\left(\mathbf{r}\right)}{\partial x} = \left(1 - C{\left(\mathbf{r} \right)}\right)

src/environment/Cavity.h

@@ -34,7 +34,7 @@
 namespace mrchem {
 /** @class Cavity
  * @brief Interlocking spheres cavity centered on the nuclei of the molecule.
- * The Cavity class represents the following function \cite Fosso-Tande2013
+ * @details The Cavity class represents the following function \cite Fosso-Tande2013
  *
  * \f[
  *    C(\mathbf{r})   = 1 - \prod^N_{i=1} (1-C_i(\mathbf{r})) \\
@@ -59,7 +59,7 @@ namespace mrchem {
  *  - \f$R_{0,i}\f$ is the atomic radius. By default, the van der Waals radius.
  *  - \f$\alpha_{i}\f$ is a scaling factor. By default, 1.1
  *  - \f$\beta_{i}\f$ is a width scaling factor. By default, 0.5
- *  - \f$\sigma_{i}\f$ is the width. By default, 0.2
+ *  - \f$\sigma_{i}\f$ is the width. By default, 0.2 bohr
  */
 
 class Cavity final : public mrcpp::RepresentableFunction<3> {

src/environment/DHScreening.h

@@ -34,10 +34,19 @@ namespace mrchem {
 /** @class DHScreening
  *
  * @brief Square of the Debye-Huckel Screening parameter.
- * TODO write proper docs for this function
- *
- * Here the \f$\bar{\kappa}\f$ is dependent on a separate cavity function with its own radius.
- *
+ * @details
+ * This is used for the Poisson-Boltzmann solver. The DHScreening function is defined as
+ * \f[
+ * \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}
+ * \f]
+ * This can be parametrized a number of ways. The one used here is
+ * \f[
+ * \kappa^2(\mathbf{r}) = (1 - C_{ion}(\mathbf{r})) \kappa^2_{out}
+ * \f]
+ * Where \f$C_{ion}(\mathbf{r})\f$ is the ion accessible Cavity function.
  */
 
 class Cavity;
@@ -47,40 +56,38 @@ class DHScreening final : public mrcpp::RepresentableFunction<3> {
     /** @brief Standard constructor. Initializes the #cavity_ion and #kappa_out with the input parameters.
      *  @param cavity_ion interlocking spheres of Cavity class.
      *  @param kappa_out value of the screening function outside the #cavity_ion.
-     *  @param formulation Decides which formulation of the #DHScreening function to implement, only continuous screening function available
-     * available as of now.
+     *  @param formulation Decides which formulation of the #DHScreening function to implement, only continuous screening function available.
+     * @details #kappa_out is given by
+     * \f[
+     * \kappa = \sqrt{\frac{2000 I_{0} e^2 N_a I_0}{\epsilon_{out} \epsilon_{in} k_B T}}
+     * \f]
+     * where \f$N_a\f$ is the Avogadro constant, e is the elementary charge, \f$I_0\f$ is the concentration of the ions,
+     * \f$k_B\f$ is the Boltzmann constant, \f$T\f$ is the temperature, \f$\epsilon_{out}\f$ is the permittivity of the solvent and \f$\epsilon_{in}\f$ is the permittivity of free space.
      */
-    DHScreening(const Cavity & cavity_ion, double kappa_out, const std::string & formulation);
+    DHScreening(const Cavity &cavity_ion, double kappa_out, const std::string &formulation);
 
-    /** @brief Evaluates DHScreening at a point in 3D space with respect to the state of #inverse.
+    /** @brief Evaluates DHScreening at a point in 3D space.
      *  @param r coordinates of a 3D point in space.
-     *  @return \f$\frac{1}{\epsilon(\mathbf{r})}\f$ if #inverse is true, and \f$ \epsilon(\mathbf{r})\f$ if #inverse is
-     *  false.
+     *  @return Value at point r.
      */
     double evalf(const mrcpp::Coord<3> &r) const override;
 
-    /** @brief Calls the Cavity::getCoordinates() method of the #cavity instance. */
     auto getCoordinates() const { return this->cavity_ion.getCoordinates(); }
 
-    /** @brief Calls the Cavity::getRadii() method of the #cavity instance. */
     auto getRadii() const { return this->cavity_ion.getRadii(); }
 
-    /** @brief Returns the value of #kappa_out. */
     auto getKOut() const { return this->kappa_out; }
 
-    /** @brief Returns the cavity */
     Cavity getCavity() const { return this->cavity_ion; }
 
-    /** @brief Returns the formulation */
     std::string getFormulation() const { return this->formulation; }
 
-    /** @brief Print parameters */
     void printParameters() const;
 
 private:
-    double kappa_out;        //!< Dielectric constant describing the permittivity of the solvent.
+    double kappa_out;        //!< value of the Debye-Huckel screening constant outside the ion accessible cavity.
     std::string formulation; //!< Formulation of the DHScreening function. Only linear variable is used now.
-    Cavity cavity_ion;       //!< A Cavity class instance.
+    Cavity cavity_ion;       //!< A Cavity class instance which describes the ion accessible cavity.
 };
 
 } // namespace mrchem

src/environment/GPESolver.h

@@ -32,7 +32,27 @@
 
 namespace mrchem {
 /** @class GPESolver
- *  @brief class that performs the computation of the  ReactionPotential, named Self Consistent Reaction Field.
+ *  @brief Solves the Generalized Poisson equation iteratively
+ *  @details The Generalized Poisson equation is given by
+ * \f[
+ * \nabla \cdot \left( \epsilon(\mathbf{r}) \nabla V(\mathbf{r}) \right) = -4\pi \rho(\mathbf{r})
+ * \f]
+ * where \f$\epsilon(\mathbf{r})\f$ is the permittivity, \f$V(\mathbf{r})\f$ is the total electrostatic potential and \f$\rho(\mathbf{r})\f$ is the molecular charge density defined as:
+ * \f[
+ * \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r})
+ * \f]
+ * where \f$\rho_{el}\f$ is the electronic charge density and \f$\rho_{nuc}\f$ 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 \f$\mathcal{P}\f$  :cite:`Fosso-Tande2013`
+ * \f[
+ * V_R(\mathbf{r}) = \mathcal{P} \star \left[ \rho_{eff}(\mathbf{r}) - \rho(\mathbf{r}) + \gamma_s(\mathbf{r}) \right]
+ * \f]
+ * where \f$\gamma_s(\mathbf{r})\f$ is the surface charge distribution describing the polarization at the surface, \f$\rho_{eff}(\mathbf{r})\f$ is the effective charge density given by
+ * \f$\frac{\rho(\mathbf{r})}{\epsilon(\mathbf{r})}\f$ and \f$V_R(\mathbf{r})\f$ 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`
  */
 class Nuclei;
 class KAIN;
@@ -49,6 +69,12 @@ class GPESolver {
 
     ~GPESolver();
 
+    /** @brief Sets the convergence threshold for the micro-iterations, used with dynamic thresholding.
+     *  @param prec value to set the convergence threshold to
+     *  @return the current convergence threshold.
+     *  @details 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.
+     */
     double setConvergenceThreshold(double prec);
 
     mrcpp::ComplexFunction &getCurrentReactionPotential() { return this->Vr_n; }
@@ -58,12 +84,23 @@ class GPESolver {
     void updateMOResidual(double const err_t) { this->mo_residual = err_t; }
 
     friend class ReactionPotential;
+
+    /** @brief Computes the energy contributions from the reaction potential
+     * @param rho_el the electronic charge density
+     * @return a tuple containing the electronic and nuclear energy contributions
+     * @details We compute the reaction energy through the following integral:
+     * \f[
+     * 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}
+     * \f]
+     * Each term represents the electronic and nuclear contributions to the reaction energy, respectively.
+     * We compute each term separately, and return a tuple containing both.
+     */
     auto computeEnergies(const Density &rho_el) -> std::tuple<double, double>;
 
 protected:
     void clear();
     bool dynamic_thrs;
-    std::string density_type;
+    std::string density_type; //!< Decides which density we will use for computing the reaction potential, options are ``total``, ``electronic`` and ``nuclear``.
     std::string solver_name = "Generalized Poisson";
 
     int max_iter;
@@ -79,24 +116,93 @@ class GPESolver {
     // another one could be to define a representable function which only has the exact analytical form of the nuclear contribution.
     // Since we already have \nabla^2 V_nuc = -4*pi*rho_nuc (nuclear potential) we could use this as a way to bypass computing rho_nuc at all
     // Same with the coulomb potential, which is basically what is left from V_vac after subtracting V_nuc. in one way we could just precompute both and
-    // just iterate through V_R only. Only issue here is (1 -1\varepsilon)/\varepsilon * \rho_nuc as I am not sure how to represent this as an analytitcal function,
+    // just iterate through V_R only. Only issue here is (1 -1\epsilon)/\epsilon * \rho_nuc as I am not sure how to represent this as an analytitcal function,
     // maybe after applying the Poisson operator?
 
     mrcpp::ComplexFunction Vr_n;
 
     std::shared_ptr<mrcpp::DerivativeOperator<3>> derivative;
     std::shared_ptr<mrcpp::PoissonOperator> poisson;
 
+    /** @brief computes density wrt. the density_type variable
+     * @param Phi the molecular orbitals
+     * @param rho_out Density function in which the density will be computed.
+     * @details The total charge density is given by the sum of the electronic and nuclear charge densities:
+     * \f[
+     * \rho(\mathbf{r}) = \rho_{el}(\mathbf{r}) + \rho_{nuc}(\mathbf{r})
+     * \f]
+     * where \f$\rho_{el}\f$ is the electronic charge density and \f$\rho_{nuc}\f$ 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``.
+     */
     void computeDensities(OrbitalVector &Phi, Density &rho_out);
+
+    /** @brief Computes the surface charge distibution due to polarization at the solute-solvent boundary
+     * @param potential Potential used to compute \f$\nabla V(\mathbf{r})\f$
+     * @param out_gamma ComplexFunction in which the surface charge distribution will be computed.
+     * @details The surface charge distribution is given by
+     * \f[
+     * \gamma_s(\mathbf{r}) = \frac{\log \frac{\epsilon_{in}}{\epsilon_{out}}}{4 \pi } \left( \nabla C(\mathbf{r}) \cdot \nabla V(\mathbf{r})\right)
+     * \f]
+     * where \f$\epsilon_{in}\f$ is the permittivity inside the cavity and \f$\epsilon_{out}\f$ is the permittivity outside the cavity.
+     */
     virtual void computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma);
 
+    /** @brief Iterates once through the Generalized Poisson equation to compute the reaction potential
+     * @param ingamma the surface charge distribution
+     * @param Phi the molecular orbitals
+     * @return the reaction potential
+     * @details Constructs the effective charge density \f$\rho_{eff}(\mathbf{r})\f$ and the Poisson operator \f$\mathcal{P}\f$ as:
+     * \f[
+     * V_R(\mathbf{r}) = \mathcal{P} \star \left[ \rho_{eff}(\mathbf{r}) - \rho(\mathbf{r}) + \gamma_s(\mathbf{r}) \right]
+     * \f]
+     * where \f$\gamma_s(\mathbf{r})\f$ is the surface charge distribution describing the polarization at the surface, \f$\rho_{eff}(\mathbf{r})\f$ is the effective charge density given by
+     * \f$\frac{\rho(\mathbf{r})}{\epsilon(\mathbf{r})}\f$ and \f$V_R(\mathbf{r})\f$ is the reaction potential.
+     */
     mrcpp::ComplexFunction solvePoissonEquation(const mrcpp::ComplexFunction &ingamma, mrchem::OrbitalVector &Phi);
 
+    /** @brief Uses KAIN to accelerate convergece of the reaction potential
+     * @param dfunc the current update of the reaction potential
+     * @param func the current reaction potential
+     * @param kain the KAIN object
+     */
     void accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain);
 
+    /** @brief Iterates through the application of the Poisson operator to Solve the Generalized Poisson equation
+     *  @param V_vac the vacuum potential
+     *  @param Phi_p the molecular orbitals
+     * @details 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 \f$V(\mathbf{r}) = V_{vac}(\mathbf{r}) + V_R(\mathbf{r})\f$
+     *  -# Compute the surface charge distribution \f$\gamma_s(\mathbf{r})\f$ with #computeGamma
+     *  -# Compute a new reaction potential \f$V_R(\mathbf{r})\f$ by application of the Poisson operator with #solvePoissonEquation
+     *  -# Calculate the update of the reaction potential as \f$\Delta V_R(\mathbf{r}) = V_R(\mathbf{r}) - V_R^{old}(\mathbf{r})\f$
+     *  -# Accelerate convergence of the reaction potential through KAIN
+     *  -# Update the reaction potential as \f$V_R(\mathbf{r}) = V_R^{old}(\mathbf{r}) + \Delta V_R(\mathbf{r})\f$
+     *  -# Check if the reaction potential has converged, if not, repeat from step 1.
+     */
     void runMicroIterations(const mrcpp::ComplexFunction &V_vac, std::shared_ptr<mrchem::OrbitalVector> Phi_p);
+
+    /** @brief Setups and computes the reaction potential through the microiterations
+     * @param V_vac the vacuum potential
+     * @param Phi_p the molecular orbitals
+     * @return The converged reaction potential for the current SCF iteration
+     * @details An initial guess of the reaction potential is computed with the following steps:
+     * -# Set the total potential as \f$V(\mathbf{r}) = V_{vac}(\mathbf{r})\f$
+     * -# Compute the surface charge distribution \f$\gamma_s(\mathbf{r})\f$ from this potential
+     * -# Iterate once through the application of the Poisson operator to return the initial guess of the reaction potential \f$V_R(\mathbf{r})\f$
+     *
+     * 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.
+     */
     mrcpp::ComplexFunction &solveEquation(double prec, const std::shared_ptr<mrchem::OrbitalVector> &Phi_p);
 
+    /** @brief Frees the memory used by the FunctionTrees of the input Complexfunction and reallocates them.
+     * @param function the ComplexFunction to reset
+     * @details This is done to avoid memory leaks when the ComplexFunction is used in the micro-iterations.
+     */
     void resetComplexFunction(mrcpp::ComplexFunction &function);
 
     virtual void printParameters() const;