Bugfix: critical sign error in magnetic properties

src/qmoperators/QMDerivative.cpp

@@ -79,12 +79,11 @@ Orbital QMDerivative::apply(Orbital inp) {
         if (inp.hasImag()) {
             out.alloc(NUMBER::Real);
             mrcpp::apply(out.real(), D, inp.imag(), dir);
-            if (inp.conjugate()) out.real().rescale(-1.0);
+            if (!inp.conjugate()) out.real().rescale(-1.0);
         }
         if (inp.hasReal()) {
             out.alloc(NUMBER::Imag);
             mrcpp::apply(out.imag(), D, inp.real(), dir);
-            out.imag().rescale(-1.0);
         }
     }
     return out;
@@ -121,12 +120,11 @@ QMOperatorVector QMDerivative::apply(QMOperator_p &O) {
             if (V_inp->hasImag()) {
                 V_out->alloc(NUMBER::Real);
                 mrcpp::apply(V_out->real(), D, V_inp->imag(), d);
-                if (V_inp->conjugate()) V_out->real().rescale(-1.0);
+                if (!V_inp->conjugate()) V_out->real().rescale(-1.0);
             }
             if (V_inp->hasReal()) {
                 V_out->alloc(NUMBER::Imag);
                 mrcpp::apply(V_out->imag(), D, V_inp->real(), d);
-                V_out->imag().rescale(-1.0);
             }
         }
         out.push_back(V_out);

src/qmoperators/one_electron/MomentumOperator.h

@@ -27,22 +27,21 @@
 
 #include "tensor/RankOneOperator.h"
 
-#include "qmoperators/QMDerivative.h"
+#include "qmoperators/one_electron/NablaOperator.h"
 
 namespace mrchem {
 
 class MomentumOperator final : public RankOneOperator<3> {
 public:
-    MomentumOperator(std::shared_ptr<mrcpp::DerivativeOperator<3>> D) {
-        auto p_x = std::make_shared<QMDerivative>(0, D, true);
-        auto p_y = std::make_shared<QMDerivative>(1, D, true);
-        auto p_z = std::make_shared<QMDerivative>(2, D, true);
+    MomentumOperator(std::shared_ptr<mrcpp::DerivativeOperator<3>> D)
+            : MomentumOperator(NablaOperator(D, true)) {}
 
+    MomentumOperator(NablaOperator D) {
         // Invoke operator= to assign *this operator
         RankOneOperator<3> &p = (*this);
-        p[0] = p_x;
-        p[1] = p_y;
-        p[2] = p_z;
+        p[0] = -1.0*D[0];
+        p[1] = -1.0*D[1];
+        p[2] = -1.0*D[2];
         p[0].name() = "p[x]";
         p[1].name() = "p[y]";
         p[2].name() = "p[z]";

src/qmoperators/one_electron/NablaOperator.h

@@ -33,10 +33,10 @@ namespace mrchem {
 
 class NablaOperator final : public RankOneOperator<3> {
 public:
-    NablaOperator(std::shared_ptr<mrcpp::DerivativeOperator<3>> D) {
-        auto d_x = std::make_shared<QMDerivative>(0, D);
-        auto d_y = std::make_shared<QMDerivative>(1, D);
-        auto d_z = std::make_shared<QMDerivative>(2, D);
+    NablaOperator(std::shared_ptr<mrcpp::DerivativeOperator<3>> D, bool imag = false) {
+        auto d_x = std::make_shared<QMDerivative>(0, D, imag);
+        auto d_y = std::make_shared<QMDerivative>(1, D, imag);
+        auto d_z = std::make_shared<QMDerivative>(2, D, imag);
 
         // Invoke operator= to assign *this operator
         RankOneOperator<3> &d = (*this);

src/scf_solver/HelmholtzVector.cpp

@@ -152,7 +152,7 @@ Orbital HelmholtzVector::apply(int i, Orbital &phi) const {
     if (phi.hasImag()) {
         out.alloc(NUMBER::Imag);
         mrcpp::apply(this->prec, out.imag(), H, phi.imag(), -1, true); // Absolute prec
-        double sign = (phi.conjugate()) ? -1.0 : 1.0;
+        double sign = (phi.conjugate()) ? 1.0 : -1.0;
         out.imag().rescale(sign / (2.0 * mrcpp::pi));
     }
     return out;