project and store kappa minus one

MRChemSoft · Aug 9, 2024 · 6c4c3d5 · 6c4c3d5
1 parent 1f3200f
commit 6c4c3d5
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 12 deletions.
diff --git a/src/qmoperators/one_electron/ZoraKineticOperator.h b/src/qmoperators/one_electron/ZoraKineticOperator.h
@@ -40,7 +40,7 @@ class ZoraKineticOperator final : public RankZeroOperator {
     ZoraKineticOperator(MomentumOperator p, ZoraOperator kappa) {
         // Invoke operator= to assign *this operator
         RankZeroOperator &t = (*this);
-        t = 0.5 * (p[0] * kappa * p[0] + p[1] * kappa * p[1] + p[2] * kappa * p[2]);
+        t = 0.5 * (p[0] * kappa * p[0] + p[1] * kappa * p[1] + p[2] * kappa * p[2]) + 0.5 * (p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
         t.name() = "T_zora";
     }
 };

diff --git a/src/qmoperators/one_electron/ZoraOperator.cpp b/src/qmoperators/one_electron/ZoraOperator.cpp
@@ -49,15 +49,16 @@ ZoraOperator::ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inv
         if (inverse) {
             k->real().map([two_cc](double val) { return (two_cc - val) / two_cc - 1.0; });
         } else {
-            k->real().map([two_cc](double val) { return two_cc / (two_cc - val) - 1.0; });
+            // k->real().map([two_cc](double val) { return two_cc / (two_cc - val) - 1.0; });
+            k->real().map([two_cc](double val) { return (val) / (two_cc - val); });
         }
         k->real().crop(proj_prec);
-        auto const_analytic = [](const mrcpp::Coord<3>& r) {
-            return 1.0;
-        };
-        mrcpp::ComplexFunction const_func;
-        mrcpp::cplxfunc::project(const_func, const_analytic, mrcpp::NUMBER::Real, proj_prec);
-        k->add(1.0, const_func);
+        // auto const_analytic = [](const mrcpp::Coord<3>& r) {
+        //     return 1.0;
+        // };
+        // mrcpp::ComplexFunction const_func;
+        // mrcpp::cplxfunc::project(const_func, const_analytic, mrcpp::NUMBER::Real, proj_prec);
+        // k->add(1.0, const_func);
     }
 
     RankZeroOperator &kappa = (*this);

diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp
@@ -182,7 +182,7 @@ SCFEnergy FockBuilder::trace(OrbitalVector &Phi, const Nuclei &nucs) {
 
     // Kinetic part
     if (isZora()) {
-        E_kin = qmoperator::calc_kinetic_trace(momentum(), *this->kappa, Phi).real();
+        E_kin = qmoperator::calc_kinetic_trace(momentum(), *this->kappa, Phi).real() + qmoperator::calc_kinetic_trace(momentum(), Phi);
     } else {
         E_kin = qmoperator::calc_kinetic_trace(momentum(), Phi);
     }
@@ -207,7 +207,7 @@ ComplexMatrix FockBuilder::operator()(OrbitalVector &bra, OrbitalVector &ket) {
 
     ComplexMatrix T_mat = ComplexMatrix::Zero(bra.size(), ket.size());
     if (isZora()) {
-        T_mat = qmoperator::calc_kinetic_matrix(momentum(), *this->kappa, bra, ket);
+        T_mat = qmoperator::calc_kinetic_matrix(momentum(), *this->kappa, bra, ket) + qmoperator::calc_kinetic_matrix(momentum(), bra, ket);
     } else {
         T_mat = qmoperator::calc_kinetic_matrix(momentum(), bra, ket);
     }
@@ -258,7 +258,7 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita
     nabla.setup(prec);
 
     RankZeroOperator operOne = 0.5 * tensor::dot(nabla(kappa), p);
-    RankZeroOperator operThree = kappa * V_zora;
+    RankZeroOperator operThree = kappa * V_zora + V_zora;
     operOne.setup(prec);
     operThree.setup(prec);
 
@@ -302,7 +302,11 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita
     operOne.clear();
 
     Timer t_kappa;
-    auto out = kappa_m1(arg);
+    mrchem::OrbitalVector out = kappa_m1(arg);
+    for (int i = 0; i < arg.size(); i++) {
+        if (not mrcpp::mpi::my_orb(out[i])) continue;
+        out[i].add(1.0, arg[i]);
+    }
     mrcpp::print::time(2, "Applying kappa inverse", t_kappa);
     nabla.clear();
     return out;