diff --git a/src/MGmol.h b/src/MGmol.h
index b18465a1..350e8511 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -341,6 +341,11 @@ class MGmol : public MGmolInterface
     {
         forces_->force(orbitals, ions);
     }
+    void setPositions(const std::vector<double>& positions,
+        const std::vector<short>& atnumbers)
+    {
+        ions_->setPositions(positions, atnumbers);
+    }
 
     /*
      * simply dump current state
diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h
index 9a9bf8a6..dc54a8e7 100644
--- a/src/MGmolInterface.h
+++ b/src/MGmolInterface.h
@@ -40,6 +40,9 @@ class MGmolInterface
 
     virtual void getAtomicPositions(std::vector<double>& tau) = 0;
     virtual void getAtomicNumbers(std::vector<short>& an)     = 0;
+    virtual void setPositions(const std::vector<double>& positions,
+        const std::vector<short>& atnumbers)
+        = 0;
     virtual std::shared_ptr<ProjectedMatricesInterface> getProjectedMatrices()
         = 0;
     virtual void dumpRestart() = 0;
diff --git a/tests/RestartEnergyAndForces/test.py b/tests/RestartEnergyAndForces/test.py
index b62d39f8..585722f0 100755
--- a/tests/RestartEnergyAndForces/test.py
+++ b/tests/RestartEnergyAndForces/test.py
@@ -64,7 +64,13 @@
 shutil.rmtree('WF')
 
 test_energy=1.e18
+l=-1
 for line in lines:
+  if line.count(b'Positions'):
+    l=0
+  if l>=0 and l<4:
+    print(line)
+    l=l+1
   if line.count(b'%%'):
     print(line)
     words=line.split()
diff --git a/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc b/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc
index e323afbc..3027ab55 100644
--- a/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc
+++ b/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc
@@ -124,6 +124,8 @@ int main(int argc, char** argv)
             }
         }
 
+        mgmol->setPositions(positions, anumbers);
+
         Mesh* mymesh             = Mesh::instance();
         const pb::Grid& mygrid   = mymesh->grid();
         const pb::PEenv& myPEenv = mymesh->peenv();
@@ -165,6 +167,36 @@ int main(int argc, char** argv)
         projmatrices->setDMuniform(ct.getNelSpin(), 0);
         projmatrices->printDM(std::cout);
 
+        // swap H and O to make sure order of atoms in list does not matter
+        double x     = positions[0];
+        double y     = positions[1];
+        double z     = positions[2];
+        positions[0] = positions[3];
+        positions[1] = positions[4];
+        positions[2] = positions[5];
+        positions[3] = x;
+        positions[4] = y;
+        positions[5] = z;
+        short tmp    = anumbers[0];
+        anumbers[0]  = anumbers[1];
+        anumbers[1]  = tmp;
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Positions:" << std::endl;
+            std::vector<short>::iterator ita = anumbers.begin();
+            for (std::vector<double>::iterator it = positions.begin();
+                 it != positions.end(); it += 3)
+            {
+                std::cout << *ita;
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+                ita++;
+            }
+        }
+
+        mgmol->setPositions(positions, anumbers);
+
         //
         // evaluate energy and forces with wavefunctions just read
         //