mpi working

Author: moritzgubler

src/qmoperators/two_electron/XCPotential.cpp

@@ -64,7 +64,7 @@ void XCPotential::setup(double prec) {
 
     auto &grid = this->mrdft->grid().get();
     mrcpp::FunctionTreeVector<3> xc_inp = setupDensities(prec, grid);
-    xc_out = this->mrdft->evaluate(xc_inp);
+    mrcpp::FunctionTreeVector<3> xc_out = this->mrdft->evaluate(xc_inp);
 
     // Fetch energy
     this->energy = this->mrdft->functional().XCenergy;
@@ -97,6 +97,7 @@ void XCPotential::setup(double prec) {
         auto t = timer.elapsed();
         mrcpp::print::tree(2, "XC operator", totNodes, totSize, t);
     }
+    mrcpp::clear(xc_out, true);
     mrcpp::print::footer(3, timer, 2);
 }
 
@@ -105,7 +106,6 @@ void XCPotential::clear() {
     this->energy = 0.0;
     for (auto &rho : this->densities) rho.free(NUMBER::Total);
     mrcpp::clear(this->potentials, true);
-    mrcpp::clear(xc_out, true);
     clearApplyPrec();
 }
 

src/qmoperators/two_electron/XCPotential.h

@@ -63,8 +63,11 @@ class XCPotential : public QMPotential {
             , mrdft(std::move(F)) {}
     ~XCPotential() override = default;
 
+    /**
+     * @brief Get the XC potential. For unrestricted calculations, the potential is a vector of two functions.
+     */
     std::shared_ptr<mrcpp::FunctionTreeVector<3>> getPotentialVector() { 
-        return std::make_shared<mrcpp::FunctionTreeVector<3>>(xc_out); 
+        return std::make_shared<mrcpp::FunctionTreeVector<3>>(potentials); 
     }
 
     friend class XCOperator;
@@ -88,7 +91,6 @@ class XCPotential : public QMPotential {
     Orbital apply(Orbital phi) override;
     Orbital dagger(Orbital phi) override;
     QMOperatorVector apply(std::shared_ptr<QMOperator> &O) override;
-    mrcpp::FunctionTreeVector<3> xc_out;
 };
 
 } // namespace mrchem

src/surface_forces/xcStress.cpp

@@ -81,7 +81,7 @@ std::vector<Matrix3d> xcGGAStress(unique_ptr<mrdft::MRDFT> &mrdft_p, mrcpp::Func
         pos[1] = gridPos(i, 1);
         pos[2] = gridPos(i, 2);
         for (int j = 0; j < 3; j++) {
-            out[i](j, j) = xcOUT(i, 0) - rhoGrid(i) * std::get<1>(xc_pots[1])->evalf(pos);
+            out[i](j, j) = xcOUT(i, 0) - rhoGrid(i) * std::get<1>(xc_pots[0])->evalf(pos);
         }
         for (int j1 = 0; j1 < 3; j1++) {
             for (int j2 = 0; j2 < 3; j2++) {
@@ -122,7 +122,7 @@ std::vector<Matrix3d> xcGGASpinStress(unique_ptr<mrdft::MRDFT> &mrdft_p, mrcpp::
         pos[1] = gridPos(i, 1);
         pos[2] = gridPos(i, 2);
         for (int j = 0; j < 3; j++) {
-            out[i](j, j) = xc(i, 0) - std::get<1>(xc_pots[1])->evalf(pos) * rhoGridAlpha(i) - std::get<1>(xc_pots[2])->evalf(pos) * rhoGridBeta(i);
+            out[i](j, j) = xc(i, 0) - std::get<1>(xc_pots[0])->evalf(pos) * rhoGridAlpha(i) - std::get<1>(xc_pots[1])->evalf(pos) * rhoGridBeta(i);
         }
         for (int j1 = 0; j1 < 3; j1++) {
             for (int j2 = 0; j2 < 3; j2++) {