Use components

external/upstream/fetch_mrcpp.cmake

@@ -39,7 +39,7 @@ else()
     GIT_REPOSITORY
       https://github.com/MRChemSoft/mrcpp.git
     GIT_TAG
-      720133372c9717134c5a01e963cb9804a1e8c36e
+      d9b07c274b0d8d87cc57ecf23fc9f4139467fc6d
   )
 
   FetchContent_GetProperties(mrcpp_sources)

src/analyticfunctions/NuclearGradientFunction.h

@@ -31,7 +31,7 @@
 
 namespace mrchem {
 
-class NuclearGradientFunction : public mrcpp::RepresentableFunction<3> {
+    class NuclearGradientFunction : public mrcpp::RepresentableFunction<3, double> {
 public:
     /*!
      *  @brief NuclearGradientFunction represents the function: Z * [x,y,z]/|r - o|^3

src/driver.cpp

@@ -390,7 +390,7 @@ bool driver::scf::guess_orbitals(const json &json_guess, Molecule &mol) {
         success = false;
     }
     for (const auto &phi_i : Phi) {
-        double err = (mrcpp::mpi::my_orb(phi_i)) ? std::abs(phi_i.norm() - 1.0) : 0.0;
+        double err = (mrcpp::mpi::my_func(phi_i)) ? std::abs(phi_i.norm() - 1.0) : 0.0;
         if (err > 0.01) MSG_WARN("MO not normalized!");
     }
 
@@ -622,7 +622,7 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) {
         if (line) plt.linePlot(npts, rho, fname);
         if (surf) plt.surfPlot(npts, rho, fname);
         if (cube) plt.cubePlot(npts, rho, fname);
-        rho.free(NUMBER::Total);
+        rho.free();
         mrcpp::print::time(1, fname, t_lap);
 
         if (orbital::size_singly(Phi) > 0) {
@@ -633,7 +633,7 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) {
             if (surf) plt.surfPlot(npts, rho, fname);
             if (cube) plt.cubePlot(npts, rho, fname);
             mrcpp::print::time(1, fname, t_lap);
-            rho.free(NUMBER::Total);
+            rho.free();
 
             t_lap.start();
             fname = path + "/rho_a";
@@ -642,15 +642,15 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) {
             if (surf) plt.surfPlot(npts, rho, fname);
             if (cube) plt.cubePlot(npts, rho, fname);
             mrcpp::print::time(1, fname, t_lap);
-            rho.free(NUMBER::Total);
+            rho.free();
 
             t_lap.start();
             fname = path + "/rho_b";
             density::compute(-1.0, rho, Phi, DensityType::Beta);
             if (line) plt.linePlot(npts, rho, fname);
             if (surf) plt.surfPlot(npts, rho, fname);
             if (cube) plt.cubePlot(npts, rho, fname);
-            rho.free(NUMBER::Total);
+            rho.free();
             mrcpp::print::time(1, fname, t_lap);
         }
     }
@@ -660,10 +660,10 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) {
         if (orb_idx[0] < 0) {
             // Plotting ALL orbitals
             for (auto i = 0; i < Phi.size(); i++) {
-                if (not mrcpp::mpi::my_orb(Phi[i])) continue;
+                if (not mrcpp::mpi::my_func(Phi[i])) continue;
                 t_lap.start();
                 std::stringstream name;
-                name << path << "/phi_" << Phi[i].printSpin() << "_scf_idx_" << i;
+                name << path << "/phi_" << Orbital(Phi[i]).printSpin() << "_scf_idx_" << i;
                 if (line) plt.linePlot(npts, Phi[i], name.str());
                 if (surf) plt.surfPlot(npts, Phi[i], name.str());
                 if (cube) plt.cubePlot(npts, Phi[i], name.str());
@@ -672,7 +672,7 @@ void driver::scf::plot_quantities(const json &json_plot, Molecule &mol) {
         } else {
             // Plotting some orbitals
             for (auto &i : orb_idx) {
-                if (not mrcpp::mpi::my_orb(Phi[i])) continue;
+                if (not mrcpp::mpi::my_func(Phi[i])) continue;
                 t_lap.start();
                 std::stringstream name;
                 auto sp = 'u';

src/environment/GPESolver.cpp

@@ -60,7 +60,7 @@ GPESolver::GPESolver(const Permittivity &e, const Density &rho_nuc, PoissonOpera
         , poisson(P) {}
 
 GPESolver::~GPESolver() {
-    this->rho_nuc.free(NUMBER::Real);
+    this->rho_nuc.free();
     clear();
 }
 
@@ -80,35 +80,34 @@ void GPESolver::computeDensities(const Density &rho_el, Density &rho_out) {
 
     switch (this->density_type) {
         case SCRFDensityType::TOTAL:
-            mrcpp::cplxfunc::add(rho_out, 1.0, rho_el, 1.0, this->rho_nuc, -1.0);
+            mrcpp::add(rho_out, 1.0, rho_el, 1.0, this->rho_nuc, -1.0);
             break;
         case SCRFDensityType::ELECTRONIC:
-            // using const_cast is absolutely EVIL here and in principle shouldn't be needed!
-            // however MRCPP is not everywhere const-correct, so here we go!
-            mrcpp::cplxfunc::deep_copy(rho_out, const_cast<Density &>(rho_el));
+            mrcpp::deep_copy(rho_out, rho_el);
             break;
         case SCRFDensityType::NUCLEAR:
-            mrcpp::cplxfunc::deep_copy(rho_out, this->rho_nuc);
+            mrcpp::deep_copy(rho_out, this->rho_nuc);
             break;
         default:
             MSG_ABORT("Invalid density type");
     }
     print_utils::qmfunction(3, "Vacuum density", rho_out, timer);
 }
 
-void GPESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFunction &out_gamma) {
-
+void GPESolver::computeGamma(mrcpp::CompFunction<3> &potential, mrcpp::CompFunction<3> &out_gamma) {
     auto d_V = mrcpp::gradient(*derivative, potential.real()); // FunctionTreeVector
     resetComplexFunction(out_gamma);
 
     for (int d = 0; d < 3; d++) {
         auto C_pin = this->epsilon.getCavity_p();
         mrcpp::AnalyticFunction<3> d_cav(C_pin->getGradVector()[d]);
-        mrcpp::ComplexFunction cplxfunc_prod;
-        mrcpp::cplxfunc::multiply(cplxfunc_prod, get_func(d_V, d), d_cav, this->apply_prec, 1);
+        mrcpp::CompFunction<3> cplxfunc_prod;
+
+        mrcpp::FunctionTree<3, double>& Tree =  get_func(d_V, d);
+        mrcpp::multiply(cplxfunc_prod, Tree, d_cav, this->apply_prec, 1);
         // add result into out_gamma
         if (d == 0) {
-            mrcpp::cplxfunc::deep_copy(out_gamma, cplxfunc_prod);
+            mrcpp::deep_copy(out_gamma, cplxfunc_prod);
         } else {
             out_gamma.add(1.0, cplxfunc_prod);
         }
@@ -118,32 +117,31 @@ void GPESolver::computeGamma(mrcpp::ComplexFunction &potential, mrcpp::ComplexFu
     mrcpp::clear(d_V, true);
 }
 
-mrcpp::ComplexFunction GPESolver::solvePoissonEquation(const mrcpp::ComplexFunction &in_gamma, const Density &rho_el) {
-    mrcpp::ComplexFunction Poisson_func;
-    mrcpp::ComplexFunction rho_eff;
-    mrcpp::ComplexFunction first_term;
-    mrcpp::ComplexFunction Vr_np1;
-    Vr_np1.alloc(NUMBER::Real);
+mrcpp::CompFunction<3> GPESolver::solvePoissonEquation(const mrcpp::CompFunction<3> &in_gamma, const Density &rho_el) {
+    mrcpp::CompFunction<3> Poisson_func;
+    mrcpp::CompFunction<3> rho_eff;
+    mrcpp::CompFunction<3> first_term;
+    mrcpp::CompFunction<3> Vr_np1;
+    Vr_np1.func_ptr->isreal = 1;
+    Vr_np1.alloc(1);
 
     auto eps_inv_func = mrcpp::AnalyticFunction<3>([this](const mrcpp::Coord<3> &r) { return 1.0 / this->epsilon.evalf(r); });
     Density rho_tot(false);
     computeDensities(rho_el, rho_tot);
 
-    mrcpp::cplxfunc::multiply(first_term, rho_tot, eps_inv_func, this->apply_prec);
+    mrcpp::multiply(first_term, rho_tot, eps_inv_func, this->apply_prec);
 
-    mrcpp::cplxfunc::add(rho_eff, 1.0, first_term, -1.0, rho_tot, -1.0);
-    rho_tot.free(NUMBER::Real);
+    mrcpp::add(rho_eff, 1.0, first_term, -1.0, rho_tot, -1.0);
+    rho_tot.free();
 
-    mrcpp::cplxfunc::add(Poisson_func, 1.0, in_gamma, 1.0, rho_eff, -1.0);
+    mrcpp::add(Poisson_func, 1.0, in_gamma, 1.0, rho_eff, -1.0);
     mrcpp::apply(this->apply_prec, Vr_np1.real(), *poisson, Poisson_func.real());
     return Vr_np1;
 }
 
-void GPESolver::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::ComplexFunction &func, KAIN &kain) {
-    OrbitalVector phi_n(0);
-    OrbitalVector dPhi_n(0);
-    phi_n.push_back(Orbital(SPIN::Paired));
-    dPhi_n.push_back(Orbital(SPIN::Paired));
+void GPESolver::accelerateConvergence(mrcpp::CompFunction<3> &dfunc, mrcpp::CompFunction<3> &func, KAIN &kain) {
+    OrbitalVector phi_n(1);
+    OrbitalVector dPhi_n(1);
 
     phi_n[0] = func;
     dPhi_n[0] = dfunc;
@@ -160,7 +158,7 @@ void GPESolver::accelerateConvergence(mrcpp::ComplexFunction &dfunc, mrcpp::Comp
     dPhi_n.clear();
 }
 
-void GPESolver::runMicroIterations(const mrcpp::ComplexFunction &V_vac, const Density &rho_el) {
+void GPESolver::runMicroIterations(const mrcpp::CompFunction<3> &V_vac, const Density &rho_el) {
     KAIN kain(this->history);
     kain.setLocalPrintLevel(10);
 
@@ -171,30 +169,30 @@ void GPESolver::runMicroIterations(const mrcpp::ComplexFunction &V_vac, const De
     while (update >= this->conv_thrs && iter <= max_iter) {
         Timer t_iter;
         // solve the poisson equation
-        mrcpp::ComplexFunction V_tot;
-        mrcpp::ComplexFunction gamma_n;
-        mrcpp::ComplexFunction dVr_n;
+        mrcpp::CompFunction<3> V_tot;
+        mrcpp::CompFunction<3> gamma_n;
+        mrcpp::CompFunction<3> dVr_n;
 
-        mrcpp::cplxfunc::add(V_tot, 1.0, this->Vr_n, 1.0, V_vac, -1.0);
+        mrcpp::add(V_tot, 1.0, this->Vr_n, 1.0, V_vac, -1.0);
         computeGamma(V_tot, gamma_n);
         auto Vr_np1 = solvePoissonEquation(gamma_n, rho_el);
         norm = Vr_np1.norm();
 
         // use a convergence accelerator
 
-        mrcpp::cplxfunc::add(dVr_n, 1.0, Vr_np1, -1.0, this->Vr_n, -1.0);
+        mrcpp::add(dVr_n, 1.0, Vr_np1, -1.0, this->Vr_n, -1.0);
         update = dVr_n.norm();
 
         if (iter > 1 and this->history > 0) {
             accelerateConvergence(dVr_n, Vr_n, kain);
-            Vr_np1.free(NUMBER::Real);
-            mrcpp::cplxfunc::add(Vr_np1, 1.0, Vr_n, 1.0, dVr_n, -1.0);
+            Vr_np1.free();
+            mrcpp::add(Vr_np1, 1.0, Vr_n, 1.0, dVr_n, -1.0);
         }
 
         // set up for next iteration
         resetComplexFunction(this->Vr_n);
-        mrcpp::cplxfunc::deep_copy(this->Vr_n, Vr_np1);
-        Vr_np1.free(NUMBER::Real);
+        mrcpp::deep_copy(this->Vr_n, Vr_np1);
+        Vr_np1.free();
 
         printConvergenceRow(iter, norm, update, t_iter.elapsed());
 
@@ -233,25 +231,29 @@ void GPESolver::printConvergenceRow(int i, double norm, double update, double ti
     println(3, o_txt.str());
 }
 
-mrcpp::ComplexFunction &GPESolver::solveEquation(double prec, const Density &rho_el) {
+mrcpp::CompFunction<3> &GPESolver::solveEquation(double prec, const Density &rho_el) {
     this->apply_prec = prec;
     Density rho_tot(false);
     computeDensities(rho_el, rho_tot);
     Timer t_vac;
-    mrcpp::ComplexFunction V_vac;
-    V_vac.alloc(NUMBER::Real);
+    mrcpp::CompFunction<3> V_vac;
+    V_vac.func_ptr->isreal = 1;
+    V_vac.alloc(1);
     mrcpp::apply(this->apply_prec, V_vac.real(), *poisson, rho_tot.real());
-    rho_tot.free(NUMBER::Real);
+    rho_tot.free();
     print_utils::qmfunction(3, "Vacuum potential", V_vac, t_vac);
 
     // set up the zero-th iteration potential and gamma, so the first iteration gamma and potentials can be made
 
     Timer t_gamma;
-    if (not this->Vr_n.hasReal()) {
-        mrcpp::ComplexFunction gamma_0;
-        mrcpp::ComplexFunction V_tot;
-        computeGamma(V_vac, gamma_0);
-        this->Vr_n = solvePoissonEquation(gamma_0, rho_el);
+    if (Vr_n.Ncomp() == 0) {
+         mrcpp::CompFunction<3> gamma_0;
+         mrcpp::CompFunction<3> V_tot;
+         gamma_0.func_ptr->isreal = 1;
+         gamma_0.alloc(1);
+
+         computeGamma(V_vac, gamma_0);
+         this->Vr_n = solvePoissonEquation(gamma_0, rho_el);
     }
 
     // update the potential/gamma before doing anything with them
@@ -263,15 +265,15 @@ mrcpp::ComplexFunction &GPESolver::solveEquation(double prec, const Density &rho
 }
 
 auto GPESolver::computeEnergies(const Density &rho_el) -> std::tuple<double, double> {
-    auto Er_nuc = 0.5 * mrcpp::cplxfunc::dot(this->rho_nuc, this->Vr_n).real();
-    auto Er_el = 0.5 * mrcpp::cplxfunc::dot(rho_el, this->Vr_n).real();
+    auto Er_nuc = 0.5 * mrcpp::dot(this->rho_nuc, this->Vr_n).real();
+    auto Er_el = 0.5 * mrcpp::dot(rho_el, this->Vr_n).real();
     return {Er_el, Er_nuc};
 }
 
-void GPESolver::resetComplexFunction(mrcpp::ComplexFunction &function) {
-    if (function.hasReal()) function.free(NUMBER::Real);
-    if (function.hasImag()) function.free(NUMBER::Imag);
-    function.alloc(NUMBER::Real);
+void GPESolver::resetComplexFunction(mrcpp::CompFunction<3> &function) {
+    function.func_ptr->isreal = 1;
+    function.func_ptr->iscomplex = 0;
+    function.alloc(1);
 }
 
 void GPESolver::printParameters() const {

src/environment/GPESolver.h

@@ -129,7 +129,7 @@ class GPESolver {
     // 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;
+    mrcpp::CompFunction<3> Vr_n;
 
     std::shared_ptr<mrcpp::DerivativeOperator<3>> derivative;
     std::shared_ptr<mrcpp::PoissonOperator> poisson;
@@ -152,14 +152,14 @@ class GPESolver {
 
     /** @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.
+     * @param out_gamma CompFunction<3> 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);
+    virtual void computeGamma(mrcpp::CompFunction<3> &potential, mrcpp::CompFunction<3> &out_gamma);
 
     /** @brief Iterates once through the Generalized Poisson equation to compute the reaction potential
      * @param ingamma the surface charge distribution
@@ -172,14 +172,14 @@ class GPESolver {
      * 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, const Density &rho_el);
+    mrcpp::CompFunction<3> solvePoissonEquation(const mrcpp::CompFunction<3> &ingamma, const Density &rho_el);
 
     /** @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);
+    void accelerateConvergence(mrcpp::CompFunction<3> &dfunc, mrcpp::CompFunction<3> &func, KAIN &kain);
 
     /** @brief Iterates through the application of the Poisson operator to Solve the Generalized Poisson equation
      *  @param V_vac the vacuum potential
@@ -195,7 +195,7 @@ class GPESolver {
      *  -# 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, const Density &rho_el);
+    void runMicroIterations(const mrcpp::CompFunction<3> &V_vac, const Density &rho_el);
 
     /** @brief Setups and computes the reaction potential through the microiterations
      * @param V_vac the vacuum potential
@@ -209,13 +209,13 @@ class GPESolver {
      * 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 Density &rho_el);
+    mrcpp::CompFunction<3> &solveEquation(double prec, const Density &rho_el);
 
     /** @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.
+     * @param function the CompFunction<3> to reset
+     * @details This is done to avoid memory leaks when the CompFunction<3> is used in the micro-iterations.
      */
-    void resetComplexFunction(mrcpp::ComplexFunction &function);
+    void resetComplexFunction(mrcpp::CompFunction<3> &function);
 
     virtual void printParameters() const;
     void printConvergenceRow(int i, double norm, double update, double time) const;

src/environment/LPBESolver.cpp

@@ -56,9 +56,9 @@ LPBESolver::LPBESolver(const Permittivity &e,
                        SCRFDensityType density_type)
         : PBESolver(e, k, rho_nuc, P, D, kain_hist, max_iter, dyn_thrs, density_type) {}
 // TODO separate this for the linear and non-linear solver
-void LPBESolver::computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) {
+void LPBESolver::computePBTerm(mrcpp::CompFunction<3> &V_tot, const double salt_factor, mrcpp::CompFunction<3> &pb_term) {
     resetComplexFunction(pb_term);
-    mrcpp::cplxfunc::multiply(pb_term, V_tot, this->kappa, this->apply_prec);
+    mrcpp::multiply(pb_term, V_tot, this->kappa, this->apply_prec);
     pb_term.rescale(salt_factor / (4.0 * mrcpp::pi));
 }
 

src/environment/LPBESolver.h

@@ -61,10 +61,10 @@ class LPBESolver final : public PBESolver {
     /** @brief Computes the PB term
      * @param[in] V_tot the total potential
      * @param[in] salt_factor the salt factor deciding how much of the total concentration to include in the PB term
-     * @param[out] pb_term the ComplexFunction in which to store the result
+     * @param[out] pb_term the CompFunction<3> in which to store the result
      * @details The PB term is computed as \f$ \kappa^2 V_{tot} \f$ and returned.
      */
-    void computePBTerm(mrcpp::ComplexFunction &V_tot, const double salt_factor, mrcpp::ComplexFunction &pb_term) override;
+    void computePBTerm(mrcpp::CompFunction<3> &V_tot, const double salt_factor, mrcpp::CompFunction<3> &pb_term) override;
     std::string solver_name{"Linearized Poisson-Boltzmann"};
 };
 } // namespace mrchem