diff --git a/src/Ions.cc b/src/Ions.cc
index a372593d..cb6dddf5 100644
--- a/src/Ions.cc
+++ b/src/Ions.cc
@@ -1324,6 +1324,39 @@ void Ions::setLocalForces(
     }
 }
 
+void Ions::setLocalForces(
+    const std::vector<double>& forces, const std::vector<double>& coords)
+{
+    assert(forces.size() == coords.size());
+
+    // tolerance can be pretty loose, as long as it does not allow mix up
+    // with coordinates of other atoms
+    const double tol = 1.e-2;
+
+    // loop over global list of forces and coordinates
+    std::vector<double>::const_iterator cit = coords.begin();
+    for (auto fit = forces.begin(); fit != forces.end(); fit += 3)
+    {
+        // find possible matching ion
+        for (auto& ion : local_ions_)
+        {
+            double p[3];
+            ion->getPosition(&p[0]);
+            double d2 = (p[0] - (*cit)) * (p[0] - (*cit))
+                        + (p[1] - (*(cit + 1))) * (p[0] - (*(cit + 1)))
+                        + (p[2] - (*(cit + 2))) * (p[0] - (*(cit + 2)));
+            double d = std::sqrt(d2);
+            if (d < tol)
+            {
+                ion->set_force(0, *fit);
+                ion->set_force(1, *(fit + 1));
+                ion->set_force(2, *(fit + 2));
+            }
+        }
+        cit += 3;
+    }
+}
+
 // Writes out the postions of the ions and the current forces on them by root
 void Ions::printForcesGlobal(std::ostream& os, const int root) const
 {
diff --git a/src/Ions.h b/src/Ions.h
index f39e1fed..aaf48e79 100644
--- a/src/Ions.h
+++ b/src/Ions.h
@@ -308,6 +308,12 @@ class Ions
     void setLocalForces(const std::vector<double>& forces,
         const std::vector<std::string>& names);
 
+    /*!
+     * set forces for ions in local_ions_ based on coordinates matching
+     */
+    void setLocalForces(
+        const std::vector<double>& forces, const std::vector<double>& coords);
+
     void syncData(const std::vector<Species>& sp);
     // void syncNames(const int nions, std::vector<std::string>& local_names,
     // std::vector<std::string>& names);
diff --git a/src/MGmol.cc b/src/MGmol.cc
index 97e1ee6f..4ccb8e90 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -1503,8 +1503,6 @@ double MGmol<OrbitalsType>::evaluateDMandEnergyAndForces(Orbitals* orbitals,
 
     ions.getForces(forces);
 
-    ions_->getForces(forces);
-
     return eks;
 }
 
diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h
index dc54a8e7..3c2384fd 100644
--- a/src/MGmolInterface.h
+++ b/src/MGmolInterface.h
@@ -25,6 +25,7 @@ class MGmolInterface
     virtual int setupConstraintsFromInput(const std::string input_file) = 0;
     virtual void setup()                                                = 0;
     virtual void run()                                                  = 0;
+    virtual void cleanup()                                              = 0;
 
     virtual double evaluateEnergyAndForces(const std::vector<double>& tau,
         const std::vector<short>& atnumbers, std::vector<double>& forces)
diff --git a/tests/WFEnergyAndForces/testWFEnergyAndForces.cc b/tests/WFEnergyAndForces/testWFEnergyAndForces.cc
index c9e28d48..64765c7a 100644
--- a/tests/WFEnergyAndForces/testWFEnergyAndForces.cc
+++ b/tests/WFEnergyAndForces/testWFEnergyAndForces.cc
@@ -190,6 +190,7 @@ int main(int argc, char** argv)
             }
         }
 
+        mgmol->cleanup();
         delete mgmol;
 
     } // close main scope
diff --git a/tests/testIons.cc b/tests/testIons.cc
index 6bea5d64..33c3cd24 100644
--- a/tests/testIons.cc
+++ b/tests/testIons.cc
@@ -6,6 +6,29 @@
 
 #include <random>
 
+// check that all forces components have integer values larger than 0
+// and differ from each other
+int checkForces(std::vector<double>& forces)
+{
+    const double tol = 1.e-14;
+
+    for (auto f0 = forces.begin(); f0 != forces.end(); f0++)
+    {
+        std::cout << "f0 = " << *f0 << std::endl;
+        for (auto f1 = f0 + 1; f1 != forces.end(); f1++)
+        {
+            // make sure each force component is different
+            if (std::abs(*f0 - *f1) < tol || *f1 < tol || *f0 < tol)
+            {
+                std::cerr << "f0 = " << *f0 << ", f1 = " << *f1 << std::endl;
+                return 1;
+            }
+        }
+    }
+
+    return 0;
+}
+
 int main(int argc, char** argv)
 {
     int status = 0;
@@ -164,7 +187,7 @@ int main(int argc, char** argv)
 
     std::vector<double> forces(3 * na);
     // set forces to a different arbitrary value for each component
-    int i = 0;
+    int i = 1;
     for (auto& f : forces)
     {
         f = (double)i;
@@ -189,20 +212,29 @@ int main(int argc, char** argv)
 
     ions.getForces(forces);
     if (myrank == 0)
-        for (auto f0 = forces.begin(); f0 != forces.end(); f0++)
+    {
+        int status = checkForces(forces);
+        if (status > 0) MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
+    }
+
+    // test Ions::setLocalForces based on coordinates matching
+    {
+        std::vector<double> positions;
+        std::vector<short> anumbers;
+        ions.getPositions(positions);
+
+        ions.setLocalForces(forces, positions);
+
+        ions.printForcesGlobal(std::cout);
+
+        ions.getForces(forces);
+        if (myrank == 0)
         {
-            std::cout << "f0 = " << *f0 << std::endl;
-            for (auto f1 = f0 + 1; f1 != forces.end(); f1++)
-            {
-                // make sure each force component is different
-                if (std::abs(*f0 - *f1) < 1.e-14)
-                {
-                    std::cerr << "f0 = " << *f0 << ", f1 = " << *f1
-                              << std::endl;
-                    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
-                }
-            }
+            int status = checkForces(forces);
+            if (status > 0) MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
         }
+    }
+
     mpirc = MPI_Finalize();
     if (mpirc != MPI_SUCCESS)
     {