diff --git a/src/qmoperators/one_electron/ZoraKineticOperator.h b/src/qmoperators/one_electron/ZoraKineticOperator.h
deleted file mode 100644
index 1f763761a..000000000
--- a/src/qmoperators/one_electron/ZoraKineticOperator.h
+++ /dev/null
@@ -1,48 +0,0 @@
-/*
- * MRChem, a numerical real-space code for molecular electronic structure
- * calculations within the self-consistent field (SCF) approximations of quantum
- * chemistry (Hartree-Fock and Density Functional Theory).
- * Copyright (C) 2023 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
- *
- * This file is part of MRChem.
- *
- * MRChem is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * MRChem is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with MRChem. If not, see .
- *
- * For information on the complete list of contributors to MRChem, see:
- *
- */
-
-#pragma once
-
-#include "tensor/RankZeroOperator.h"
-
-#include "MomentumOperator.h"
-#include "ZoraOperator.h"
-
-namespace mrchem {
-
-class ZoraKineticOperator final : public RankZeroOperator {
-public:
- ZoraKineticOperator(std::shared_ptr> D, ZoraOperator kappa)
- : ZoraKineticOperator(MomentumOperator(D), kappa) {}
-
- 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.name() = "T_zora";
- }
-};
-
-} // namespace mrchem
diff --git a/src/qmoperators/one_electron/ZoraOperator.cpp b/src/qmoperators/one_electron/ZoraOperator.cpp
index e2b38f505..d4ed96a81 100644
--- a/src/qmoperators/one_electron/ZoraOperator.cpp
+++ b/src/qmoperators/one_electron/ZoraOperator.cpp
@@ -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
diff --git a/src/qmoperators/one_electron/ZoraOperator.h b/src/qmoperators/one_electron/ZoraOperator.h
index 8dd2ee7ce..2cadfb83c 100644
--- a/src/qmoperators/one_electron/ZoraOperator.h
+++ b/src/qmoperators/one_electron/ZoraOperator.h
@@ -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);
};
diff --git a/src/qmoperators/two_electron/FockBuilder.cpp b/src/qmoperators/two_electron/FockBuilder.cpp
index ad9aaa5d9..e36df3d1f 100644
--- a/src/qmoperators/two_electron/FockBuilder.cpp
+++ b/src/qmoperators/two_electron/FockBuilder.cpp
@@ -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(*vz, c, prec, false);
- this->kappa_inv = std::make_shared(*vz, c, prec, true);
+ // chi = kappa - 1. See ZoraOperator.h for more information.
+ this->chi = std::make_shared(*vz, c, prec, false);
+ this->chi_inv = std::make_shared(*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;
}
diff --git a/src/qmoperators/two_electron/FockBuilder.h b/src/qmoperators/two_electron/FockBuilder.h
index 53db4958e..34cb0ca41 100644
--- a/src/qmoperators/two_electron/FockBuilder.h
+++ b/src/qmoperators/two_electron/FockBuilder.h
@@ -105,8 +105,8 @@ class FockBuilder final {
std::shared_ptr xc{nullptr};
std::shared_ptr Ro{nullptr}; // Reaction field operator
std::shared_ptr ext{nullptr}; // Total external potential
- std::shared_ptr kappa{nullptr};
- std::shared_ptr kappa_inv{nullptr};
+ std::shared_ptr chi{nullptr};
+ std::shared_ptr chi_inv{nullptr};
std::shared_ptr collectZoraBasePotential();
OrbitalVector buildHelmholtzArgumentZORA(OrbitalVector &Phi, OrbitalVector &Psi, DoubleVector eps, double prec);