Kappa minus one (#491)

* improve things

* refine potential before kappa transform

* project and store kappa minus one

* apply product rule

* rever product rule change

* use other derivative operator

* remove unused operator

* delete unused zora_kinetic operator

* remove old comments

* add some documentation

* rename kappa to chi

---------

Co-authored-by: moritzgubler <[email protected]>

Author: ilfreddy

src/qmoperators/one_electron/ZoraKineticOperator.h

@@ -47,22 +47,22 @@ ZoraOperator::ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inv
     if (k->hasReal()) {
         mrcpp::refine_grid(k->real(), 1);
         if (inverse) {
-            k->real().map([two_cc](double val) { return (two_cc - val) / two_cc; });
+            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); });
+            k->real().map([two_cc](double val) { return (val) / (two_cc - val); });
         }
         k->real().crop(proj_prec);
     }
 
-    RankZeroOperator &kappa = (*this);
-    kappa = k;
+    RankZeroOperator &chi = (*this);
+    chi = k;
     if (inverse) {
-        kappa.name() = "kappa_m1";
+        chi.name() = "chi_inv";
     } else {
-        kappa.name() = "kappa";
+        chi.name() = "chi";
     }
     auto plevel = Printer::getPrintLevel();
-    print_utils::qmfunction(2, "ZORA operator (" + kappa.name() + ")", *k, timer);
+    print_utils::qmfunction(2, "ZORA operator (" + chi.name() + ")", *k, timer);
 }
 
 } // namespace mrchem

src/qmoperators/one_electron/ZoraOperator.cpp

@@ -31,8 +31,22 @@ namespace mrchem {
 
 class QMPotential;
 
+/**
+ * @class ZoraOperator
+ * @brief Implements chi = kappa - 1 relativistic dampening function. This has to be done in order to
+ * avoid numerical instabilities in the ZORA operator. Whenever this operator is applied, the + 1 has to be added to chi i
+ * in order to get the kappa operator. kappa * phi = chi * phi + phi
+ * This has to be done manually.
+ */
 class ZoraOperator final : public RankZeroOperator {
 public:
+    /**
+     * @brief Constructor for the ZoraOperator that contains the chi = kappa - 1 function.
+     * @param vz The potential used to calculate the kappa function.
+     * @param c Speed of light.
+     * @param proj_prec The precision of the MW projection.
+     * @param inverse If true, the inverse of the chi function is calculated.
+     */
     ZoraOperator(QMPotential &vz, double c, double proj_prec, bool inverse = false);
 };
 

src/qmoperators/one_electron/ZoraOperator.h

@@ -96,11 +96,12 @@ void FockBuilder::setup(double prec) {
         mrcpp::print::value(3, "Light speed", c, "(au)", 5);
         mrcpp::print::separator(3, '-');
         auto vz = collectZoraBasePotential();
-        this->kappa = std::make_shared<ZoraOperator>(*vz, c, prec, false);
-        this->kappa_inv = std::make_shared<ZoraOperator>(*vz, c, prec, true);
+        // chi = kappa - 1. See ZoraOperator.h for more information.
+        this->chi = std::make_shared<ZoraOperator>(*vz, c, prec, false);
+        this->chi_inv = std::make_shared<ZoraOperator>(*vz, c, prec, true);
         this->zora_base = RankZeroOperator(vz);
-        this->kappa->setup(prec);
-        this->kappa_inv->setup(prec);
+        this->chi->setup(prec);
+        this->chi_inv->setup(prec);
         this->zora_base.setup(prec);
         mrcpp::print::footer(3, t_zora, 2);
     };
@@ -120,8 +121,8 @@ void FockBuilder::clear() {
     this->potential().clear();
     this->perturbation().clear();
     if (isZora()) {
-        this->kappa->clear();
-        this->kappa_inv->clear();
+        this->chi->clear();
+        this->chi_inv->clear();
         this->zora_base.clear();
     }
 }
@@ -180,7 +181,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->chi, Phi).real() + qmoperator::calc_kinetic_trace(momentum(), Phi);
     } else {
         E_kin = qmoperator::calc_kinetic_trace(momentum(), Phi);
     }
@@ -205,7 +206,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->chi, bra, ket) + qmoperator::calc_kinetic_matrix(momentum(), bra, ket);
     } else {
         T_mat = qmoperator::calc_kinetic_matrix(momentum(), bra, ket);
     }
@@ -247,12 +248,12 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita
     double two_cc = 2.0 * c * c;
     MomentumOperator &p = momentum();
     RankZeroOperator &V = potential();
-    RankZeroOperator &kappa = *this->kappa;
-    RankZeroOperator &kappa_m1 = *this->kappa_inv;
+    RankZeroOperator &chi = *this->chi;
+    RankZeroOperator &chi_inv = *this->chi_inv;
     RankZeroOperator &V_zora = this->zora_base;
 
-    RankZeroOperator operOne = 0.5 * tensor::dot(p(kappa), p);
-    RankZeroOperator operThree = kappa * V_zora;
+    RankZeroOperator operOne = 0.5 * tensor::dot(p(chi), p);
+    RankZeroOperator operThree = chi * V_zora + V_zora;
     operOne.setup(prec);
     operThree.setup(prec);
 
@@ -296,7 +297,11 @@ OrbitalVector FockBuilder::buildHelmholtzArgumentZORA(OrbitalVector &Phi, Orbita
     operOne.clear();
 
     Timer t_kappa;
-    auto out = kappa_m1(arg);
+    mrchem::OrbitalVector out = chi_inv(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);
     return out;
 }

src/qmoperators/two_electron/FockBuilder.cpp

@@ -105,8 +105,8 @@ class FockBuilder final {
     std::shared_ptr<XCOperator> xc{nullptr};
     std::shared_ptr<ReactionOperator> Ro{nullptr};       // Reaction field operator
     std::shared_ptr<ElectricFieldOperator> ext{nullptr}; // Total external potential
-    std::shared_ptr<ZoraOperator> kappa{nullptr};
-    std::shared_ptr<ZoraOperator> kappa_inv{nullptr};
+    std::shared_ptr<ZoraOperator> chi{nullptr};
+    std::shared_ptr<ZoraOperator> chi_inv{nullptr};
 
     std::shared_ptr<QMPotential> collectZoraBasePotential();
     OrbitalVector buildHelmholtzArgumentZORA(OrbitalVector &Phi, OrbitalVector &Psi, DoubleVector eps, double prec);