some cleaning

external/upstream/fetch_mrcpp.cmake

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

src/qmfunctions/Orbital.h

@@ -69,7 +69,11 @@ class Orbital : public mrcpp::CompFunction<3> {
     void loadOrbital(const std::string &file);
 };
 
-//using OrbitalVector = mrcpp::CompFunctionVector<3>;
+// All MPI processes have a vector of full length, but
+// only "my_func" are fully defined.
+// The others orbitals (not my_func) have only basic data (spin etc) but no trees defined.
+// Vector of orbitals where all orbitals are defined for all MPI, should not have typoe
+// OrbitalVectors, but directly vector<Orbital>.
 class OrbitalVector : public mrcpp::CompFunctionVector {
 public:
     OrbitalVector(int N = 0) : mrcpp::CompFunctionVector(N) {}

src/qmfunctions/density_utils.cpp

@@ -157,14 +157,12 @@ void density::compute_local_X(double prec, Density &rho, OrbitalVector &Phi, Orb
     if (rho.Ncomp() == 0) rho.alloc(0);
 
     // Compute local density from own orbitals
-    rho.real().setZero();
     for (int i = 0; i < Phi.size(); i++) {
         if (mrcpp::mpi::my_func(Phi[i])) {
             if (not mrcpp::mpi::my_func(X[i])) MSG_ABORT("Inconsistent MPI distribution");
             Orbital phi_i = Phi[i];
             double occ = density::compute_occupation(phi_i, spin);
             if (std::abs(occ) < mrcpp::MachineZero) continue; // next orbital if this one is not occupied!
-
             Density rho_i(false);
             mrcpp::multiply(rho_i, phi_i, X[i], mult_prec);
             rho.add(2.0 * occ, rho_i);
@@ -177,8 +175,7 @@ void density::compute_local_XY(double prec, Density &rho, OrbitalVector &Phi, Or
     int N_el = orbital::get_electron_number(Phi);
     double mult_prec = prec;       // prec for rho_i = |x_i><phi_i| + |phi_i><y_i|
     double add_prec = prec / N_el; // prec for rho = sum_i rho_i
-    std::cout<<"ERROR"<<std::endl;
-    MSG_ERROR(" NOT implemented");
+
     if (Phi.size() != X.size()) MSG_ERROR("Size mismatch");
     if (Phi.size() != Y.size()) MSG_ERROR("Size mismatch");
 
@@ -198,10 +195,8 @@ void density::compute_local_XY(double prec, Density &rho, OrbitalVector &Phi, Or
             Density rho_x(false);
             Density rho_y(false);
             MSG_ERROR("Complex trees not yet included");
-            //            mrcpp::multiply(rho_x, X[i], Phi[i].dagger(), mult_prec);
-            //            mrcpp::multiply(rho_y, Phi[i], Y[i].dagger(), mult_prec);
-            mrcpp::multiply(rho_x, X[i], Phi[i], mult_prec);
-            mrcpp::multiply(rho_y, Phi[i], Y[i], mult_prec);
+            mrcpp::multiply(rho_x, Phi[i], X[i], mult_prec, false, false, true); //last true means complex conjugate of Phi[i]
+            mrcpp::multiply(rho_y, Y[i], Phi[i], mult_prec, false, false, true);//last true means complex conjugate of Y[i]
             rho.add(occ, rho_x);
             rho.add(occ, rho_y);
             rho.crop(add_prec);

src/qmfunctions/orbital_utils.cpp

@@ -68,80 +68,6 @@ OrbitalData getOrbitalData(const Orbital &orb) {
  * Orbital related standalone functions *
  ****************************************/
 
-/** @brief Compute <bra|ket> = int bra^\dag(r) * ket(r) dr.
- *
- *  Notice that the <bra| position is already complex conjugated.
- *  Alpha spin is orthogonal to beta spin, but paired orbitals are
- *  not necessarily orthogonal to alpha/beta orbitals.
- *
- */
-ComplexDouble orbital::dot(Orbital bra, Orbital ket) {
-    if ((bra.spin() == SPIN::Alpha) and (ket.spin() == SPIN::Beta)) return 0.0;
-    if ((bra.spin() == SPIN::Beta) and (ket.spin() == SPIN::Alpha)) return 0.0;
-    return mrcpp::dot(bra, ket);
-}
-
-
-/** @brief Compute <bra|ket> = int bra^\dag(r) * ket(r) dr.
- *
- *  Notice that the <bra| position is already complex conjugated.
- *  Alpha spin is orthogonal to beta spin, but paired orbitals are
- *  not necessarily orthogonal to alpha/beta orbitals.
- *
- */
-ComplexDouble orbital::dot(mrcpp::CompFunction<3> bra, mrcpp::CompFunction<3>ket) {
-    if ((bra.data().n1[0] == SPIN::Alpha) and (ket.data().n1[0] == SPIN::Beta)) return 0.0;
-    if ((bra.data().n1[0] == SPIN::Beta) and (ket.data().n1[0] == SPIN::Alpha)) return 0.0;
-    return mrcpp::dot(bra, ket);
-}
-
-/** @brief Compute the diagonal dot products <bra_i|ket_i>
- *
- * MPI: dot product is computed by the ket owner and the corresponding
- *      bra is communicated. The resulting vector is allreduced, and
- *      the foreign bra's are cleared.
- *
- */
-ComplexVector orbital::dot(OrbitalVector &Bra, OrbitalVector &Ket) {
-    if (Bra.size() != Ket.size()) MSG_ABORT("Size mismatch");
-
-    int N = Bra.size();
-    ComplexVector result = ComplexVector::Zero(N);
-    for (int i = 0; i < N; i++) {
-        // The bra is sent to the owner of the ket
-        if (mrcpp::mpi::my_func(Bra[i]) != mrcpp::mpi::my_func(Ket[i])) {
-            int tag = 8765 + i;
-            int src = (Bra[i].getRank()) % mrcpp::mpi::wrk_size;
-            int dst = (Ket[i].getRank()) % mrcpp::mpi::wrk_size;
-            if (mrcpp::mpi::my_func(Bra[i])) mrcpp::mpi::send_function(Bra[i], dst, tag, mrcpp::mpi::comm_wrk);
-            if (mrcpp::mpi::my_func(Ket[i])) mrcpp::mpi::recv_function(Bra[i], src, tag, mrcpp::mpi::comm_wrk);
-        }
-        result[i] = orbital::dot(Bra[i], Ket[i]);
-        if (not mrcpp::mpi::my_func(Bra[i])) Bra[i].free();
-    }
-    mrcpp::mpi::allreduce_vector(result, mrcpp::mpi::comm_wrk);
-    return result;
-}
-
-/** @brief Compute <bra|ket> = int |bra^\dag(r)| * |ket(r)| dr.
- *
- */
-double orbital::node_norm_dot(Orbital bra, Orbital ket) {
-    if ((bra.spin() == SPIN::Alpha) and (ket.spin() == SPIN::Beta)) return 0.0;
-    if ((bra.spin() == SPIN::Beta) and (ket.spin() == SPIN::Alpha)) return 0.0;
-    return mrcpp::node_norm_dot(bra, ket);
-}
-
-
-/** @brief Compute <bra|ket> = int |bra^\dag(r)| * |ket(r)| dr.
- *
-
-double orbital::node_norm_dot(mrcpp::CompFunction<3> bra, mrcpp::CompFunction<3> ket, bool exact) {
-    if ((bra.spin() == SPIN::Alpha) and (ket.spin() == SPIN::Beta)) return 0.0;
-    if ((bra.spin() == SPIN::Beta) and (ket.spin() == SPIN::Alpha)) return 0.0;
-    return mrcpp::node_norm_dot(bra, ket, exact);
-}*/
-
 /** @brief Compare spin and occupation of two orbitals
  *
  *  Returns true if orbital parameters are the same.
@@ -499,7 +425,7 @@ void orbital::normalize(OrbitalVector &Phi) {
 
 /** @brief In place orthogonalize against inp. Private function. */
 void orbital::orthogonalize(double prec, Orbital &&phi, Orbital psi) {
-    ComplexDouble overlap = orbital::dot(psi, phi);
+    ComplexDouble overlap = mrcpp::dot(psi, phi);
     double sq_norm = psi.getSquareNorm();
     if (std::abs(overlap) > prec) phi.add(-1.0 * overlap / sq_norm, psi);
 }

src/qmfunctions/orbital_utils.h

@@ -36,12 +36,6 @@ bool compare(const Orbital &phi_a, const Orbital &phi_b);
 int compare_spin(const Orbital &phi_a, const Orbital &phi_b);
 int compare_occupation(const Orbital &phi_a, const Orbital &phi_b);
 
-ComplexDouble dot(Orbital bra, Orbital ket);
-ComplexDouble dot(mrcpp::CompFunction<3> bra, mrcpp::CompFunction<3>ket);
-ComplexVector dot(OrbitalVector &Bra, OrbitalVector &Ket);
-double node_norm_dot(Orbital bra, Orbital ket);
-double node_norm_dot(mrcpp::CompFunction<3> bra, mrcpp::CompFunction<3>  ket, bool exact);
-
 void normalize(Orbital phi);
 OrbitalChunk get_my_chunk(OrbitalVector &Phi);
 void orthogonalize(double prec, Orbital &&phi, Orbital psi);

src/qmoperators/QMDerivative.cpp

@@ -63,7 +63,7 @@ Orbital QMDerivative::apply(Orbital inp) {
     auto dir = this->apply_dir;
     auto &D = *this->derivative;
 
-    Orbital out = inp.paramCopy();
+    Orbital out = inp.paramCopy(true);
     ComplexDouble metric[4][4];
     ComplexDouble diago = 1.0;
     if (isImag()) diago = {0.0, 1.0};

src/qmoperators/QMIdentity.cpp

@@ -37,7 +37,7 @@ namespace mrchem {
 /** Identity operator is a deep copy */
 Orbital QMIdentity::apply(Orbital inp) {
     if (this->apply_prec < 0.0) MSG_ERROR("Uninitialized operator");
-    Orbital out = inp.paramCopy();
+    Orbital out;
     mrcpp::deep_copy(out, inp);
 
     return out;

src/qmoperators/QMPotential.cpp

@@ -89,7 +89,7 @@ Orbital QMPotential::apply(Orbital inp) {
 Orbital QMPotential::dagger(Orbital inp) {
     if (this->apply_prec < 0.0) MSG_ERROR("Uninitialized operator");
 
-    Orbital out = inp.paramCopy();
+    Orbital out = inp.paramCopy(true);
     calc(out, inp, true);
     return out;
 }

src/qmoperators/QMSpin.cpp

@@ -64,7 +64,7 @@ Orbital QMSpin::apply(Orbital inp) {
             MSG_ABORT("Cannot apply spin operator on paired orbital");
     }
 
-    Orbital out = inp.paramCopy();
+    Orbital out;
     mrcpp::deep_copy(out, inp);
     out.rescale(coef);
 

src/qmoperators/two_electron/ExchangePotential.cpp

@@ -170,7 +170,7 @@ void ExchangePotential::calcExchange_kij(double prec, Orbital phi_k, Orbital phi
     // if the product is smaller than the target precision,
     // the result is expected to be negligible
     Timer timer_ij;
-    Orbital rho_ij = phi_i.paramCopy();
+    Orbital rho_ij = phi_i.paramCopy(true);
     mrcpp::multiply(rho_ij, phi_i, phi_j, prec_m1, true, true, true);
     timer_ij.stop();
     if (rho_ij.norm() < prec) return;
@@ -197,8 +197,7 @@ void ExchangePotential::calcExchange_kij(double prec, Orbital phi_k, Orbital phi
     }
     // compute V_ij = P[rho_ij]
     Timer timer_p;
-    Orbital V_ij = rho_ij.paramCopy();
-    V_ij.alloc(0);
+    Orbital V_ij = rho_ij.paramCopy(true);
     if (RealOrbitals) {
         mrcpp::apply(prec_p, *V_ij.CompD[0], P, *rho_ij.CompD[0], phi_opt_vec_real, -1, true);
     } else {

src/qmoperators/two_electron/ExchangePotentialD1.cpp

@@ -102,7 +102,7 @@ int ExchangePotentialD1::testInternal(Orbital phi_p) const {
  * particular exchange contribution has been precomputed.
  */
 Orbital ExchangePotentialD1::apply(Orbital phi_p) {
-    Orbital out_p = phi_p.paramCopy();
+    Orbital out_p = phi_p.paramCopy(true);
     if (this->apply_prec < 0.0) {
         MSG_ERROR("Uninitialized operator");
         return out_p;
@@ -272,8 +272,8 @@ void ExchangePotentialD1::setupInternal(double prec) {
             for (int i = 0; i < iorb_vec.size(); i++) {
                 int iorb = itasks[task][i];
                 Orbital &phi_i = iorb_vec[i];
-                Orbital ex_jji = phi_i.paramCopy();
-                Orbital ex_iij = phi_j.paramCopy();
+                Orbital ex_jji = phi_i.paramCopy(true);
+                Orbital ex_iij = phi_j.paramCopy(true);
 
                 // compute K_iij and K_jji in one operation
                 double j_fac = getSpinFactor(phi_i, phi_j);
@@ -311,7 +311,7 @@ void ExchangePotentialD1::setupInternal(double prec) {
             }
             // add all contributions to ex_j,
             if (mrcpp::mpi::bank_size > 0 and iijfunc_vec.size() > 0) {
-                Orbital ex_j = phi_j.paramCopy();
+                Orbital ex_j = phi_j.paramCopy(true);
                 t_add.resume();
                 mrcpp::linear_combination(ex_j, coef_vec, iijfunc_vec, prec);
                 t_add.stop();
@@ -361,7 +361,7 @@ void ExchangePotentialD1::setupInternal(double prec) {
             if (tot >= totmax) {
                 // we sum the contributions so far before fetching new ones
                 t_add.resume();
-                auto tmp_j = Ex[j].paramCopy();
+                auto tmp_j = Ex[j].paramCopy(true);
                 mrcpp::linear_combination(tmp_j, coef_vec, iijfunc_vec, prec);
                 Ex[j].add(1.0, tmp_j);
                 tmp_j.free();
@@ -374,7 +374,7 @@ void ExchangePotentialD1::setupInternal(double prec) {
         }
         if (iijfunc_vec.size() > 0) {
             t_add.resume();
-            auto tmp_j = Ex[j].paramCopy();
+            auto tmp_j = Ex[j].paramCopy(true);
             mrcpp::linear_combination(tmp_j, coef_vec, iijfunc_vec, prec);
             Ex[j].add(1.0, tmp_j);
             tmp_j.free();
@@ -426,7 +426,7 @@ Orbital ExchangePotentialD1::calcExchange(Orbital phi_p) {
 
         double spin_fac = getSpinFactor(phi_i, phi_p);
         if (std::abs(spin_fac) >= mrcpp::MachineZero) {
-            Orbital ex_iip = phi_p.paramCopy();
+            Orbital ex_iip = phi_p.paramCopy(true);
             calcExchange_kij(precf, phi_i, phi_i, phi_p, ex_iip);
             coef_vec.push_back(spin_fac / phi_i.getSquareNorm());
             func_vec.push_back(ex_iip);
@@ -436,7 +436,7 @@ Orbital ExchangePotentialD1::calcExchange(Orbital phi_p) {
     }
 
     // compute ex_p = sum_i c_i*ex_iip
-    Orbital ex_p = phi_p.paramCopy();
+    Orbital ex_p = phi_p.paramCopy(true);
     //Eigen::Map<ComplexVector> coefs(coef_vec.data(), coef_vec.size());
     mrcpp::linear_combination(ex_p, coef_vec, func_vec, prec);
     print_utils::qmfunction(4, "Applied exchange", ex_p, timer);

src/qmoperators/two_electron/ExchangePotentialD2.cpp

@@ -96,7 +96,7 @@ void ExchangePotentialD2::clearBank() {
 Orbital ExchangePotentialD2::apply(Orbital phi_p) {
     if (this->apply_prec < 0.0) {
         MSG_ERROR("Uninitialized operator");
-        return phi_p.paramCopy();
+        return phi_p.paramCopy(true);
     }
 
     Timer timer;
@@ -123,8 +123,8 @@ Orbital ExchangePotentialD2::apply(Orbital phi_p) {
 
         double spin_fac = getSpinFactor(phi_i, phi_p);
         if (std::abs(spin_fac) >= mrcpp::MachineZero) {
-            Orbital ex_xip = phi_p.paramCopy();
-            Orbital ex_iyp = phi_p.paramCopy();
+            Orbital ex_xip = phi_p.paramCopy(true);
+            Orbital ex_iyp = phi_p.paramCopy(true);
             calcExchange_kij(precf, x_i, phi_i, phi_p, ex_xip);
             calcExchange_kij(precf, phi_i, y_i, phi_p, ex_iyp);
             func_vec.push_back(ex_xip);
@@ -138,7 +138,7 @@ Orbital ExchangePotentialD2::apply(Orbital phi_p) {
     }
 
     // compute out_p = sum_i c_i*(ex_xip + ex_iyp)
-    Orbital out_p = phi_p.paramCopy();
+    Orbital out_p = phi_p.paramCopy(true);
     //Eigen::Map<ComplexVector> coefs(coef_vec.data(), coef_vec.size());
     mrcpp::linear_combination(out_p, coef_vec, func_vec, prec);
     print_utils::qmfunction(4, "Applied exchange", out_p, timer);
@@ -155,7 +155,7 @@ Orbital ExchangePotentialD2::apply(Orbital phi_p) {
 Orbital ExchangePotentialD2::dagger(Orbital phi_p) {
     if (this->apply_prec < 0.0) {
         MSG_ERROR("Uninitialized operator");
-        return phi_p.paramCopy();
+        return phi_p.paramCopy(true);
     }
 
     Timer timer;
@@ -182,8 +182,8 @@ Orbital ExchangePotentialD2::dagger(Orbital phi_p) {
 
         double spin_fac = getSpinFactor(phi_i, phi_p);
         if (std::abs(spin_fac) >= mrcpp::MachineZero) {
-            Orbital ex_ixp = phi_p.paramCopy();
-            Orbital ex_yip = phi_p.paramCopy();
+            Orbital ex_ixp = phi_p.paramCopy(true);
+            Orbital ex_yip = phi_p.paramCopy(true);
             calcExchange_kij(precf, phi_i, x_i, phi_p, ex_ixp);
             calcExchange_kij(precf, y_i, phi_i, phi_p, ex_yip);
             func_vec.push_back(ex_ixp);
@@ -197,7 +197,7 @@ Orbital ExchangePotentialD2::dagger(Orbital phi_p) {
     }
 
     // compute ex_p = sum_i c_i*(ex_ixp + ex_yip)
-    Orbital ex_p = phi_p.paramCopy();
+    Orbital ex_p = phi_p.paramCopy(true);
     //Eigen::Map<ComplexVector> coefs(coef_vec.data(), coef_vec.size());
     mrcpp::linear_combination(ex_p, coef_vec, func_vec, prec);
     print_utils::qmfunction(4, "Applied exchange", ex_p, timer);

src/scf_solver/Accelerator.cpp

@@ -166,7 +166,7 @@ bool Accelerator::verifyOverlap(OrbitalVector &Phi) {
                 auto &last_i = this->orbitals[nHistory][i];
                 if (not mrcpp::mpi::my_func(last_i)) MSG_ABORT("MPI rank mismatch");
                 auto sqNorm = phi_i.getSquareNorm();
-                auto overlap = orbital::dot(phi_i, last_i);
+                auto overlap = mrcpp::dot(phi_i, last_i);
                 if (std::abs(overlap) < 0.5 * sqNorm) {
                     mrcpp::print::value(this->pl + 2, "Overlap not verified ", std::abs(overlap));
                     out(i) = 1;

src/scf_solver/HelmholtzVector.cpp

@@ -143,8 +143,7 @@ Orbital HelmholtzVector::apply(int i, const Orbital &phi) const {
     if (std::abs(mu_i.imag()) > mrcpp::MachineZero) MSG_ABORT("Mu cannot be complex");
     mrcpp::HelmholtzOperator H(*MRA, mu_i.real(), this->prec);
 
-    Orbital out = phi.paramCopy();
-    out.alloc(0);
+    Orbital out = phi.paramCopy(true);
     ComplexDouble metric[4][4];
     for (int i=0; i<4; i++){
         for (int j=0; j<4; j++){