diff --git a/src/Electrostatic.h b/src/Electrostatic.h
index f49ad12f..619da18a 100644
--- a/src/Electrostatic.h
+++ b/src/Electrostatic.h
@@ -49,7 +49,6 @@ class Electrostatic
 
     const bool isDielectric() { return diel_flag_; }
     pb::GridFunc<RHODTYPE>* getRhoc() { return grhoc_; }
-
     Poisson* getPoissonSolver() { return poisson_solver_; }
 
     void setup(const short max_sweeps);
diff --git a/src/Ions.cc b/src/Ions.cc
index f10fb834..3b128b0e 100644
--- a/src/Ions.cc
+++ b/src/Ions.cc
@@ -109,14 +109,12 @@ Ions::Ions(const double lat[3], const std::vector<Species>& sp) : species_(sp)
 
 Ions::Ions(const Ions& ions, const double shift[3]) : species_(ions.species_)
 {
-    std::vector<Ion*>::const_iterator ion = ions.list_ions_.begin();
-    while (ion != ions.list_ions_.end())
+    for (const auto& ion : ions.list_ions_)
     {
-        Ion* newion = new Ion(**ion);
+        Ion* newion = new Ion(*ion);
         newion->shiftPosition(shift);
         newion->setup();
         list_ions_.push_back(newion);
-        ion++;
     }
     for (short i = 0; i < 3; ++i)
         lattice_[i] = ions.lattice_[i];
@@ -947,10 +945,11 @@ void Ions::initFromRestartFile(HDFrestart& h5_file)
     num_ions_ = at_names.size();
 #ifdef MGMOL_USE_HDF5P
     if (!h5_file.useHdf5p())
+#endif
     {
         mmpi.allreduce(&num_ions_, 1, MPI_SUM);
     }
-#endif
+
     if (onpe0 && ct.verbose > 0)
     {
         (*MPIdata::sout) << "Ions::setFromRestartFile(), read " << num_ions_
@@ -2211,7 +2210,7 @@ void Ions::getLocalPositions(std::vector<double>& tau) const
 void Ions::getLocalNames(std::vector<std::string>& names) const
 {
     names.clear();
-    for (auto& ion : local_ions_)
+    for (const auto& ion : local_ions_)
     {
         names.push_back(ion->name());
     }
@@ -2219,14 +2218,14 @@ void Ions::getLocalNames(std::vector<std::string>& names) const
 
 void Ions::getNames(std::vector<std::string>& names) const
 {
-    names.clear();
-    for (auto& ion : list_ions_)
-    {
-        names.push_back(ion->name());
-    }
+    std::vector<std::string> local_names;
+    getLocalNames(local_names);
+
+    MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+    mmpi.allGatherV(local_names, names);
 }
 
-void Ions::getPositions(std::vector<double>& tau)
+void Ions::getPositions(std::vector<double>& tau) const
 {
     std::vector<double> tau_local(3 * local_ions_.size());
 
@@ -2236,10 +2235,9 @@ void Ions::getPositions(std::vector<double>& tau)
     mmpi.allGatherV(tau_local, tau);
 }
 
-void Ions::getAtomicNumbers(std::vector<short>& atnumbers)
+void Ions::getAtomicNumbers(std::vector<short>& atnumbers) const
 {
     std::vector<short> local_atnumbers;
-
     for (auto& ion : local_ions_)
     {
         local_atnumbers.push_back(ion->atomic_number());
@@ -2249,15 +2247,11 @@ void Ions::getAtomicNumbers(std::vector<short>& atnumbers)
     mmpi.allGatherV(local_atnumbers, atnumbers);
 }
 
-void Ions::getForces(std::vector<double>& forces)
+void Ions::getForces(std::vector<double>& forces) const
 {
     std::vector<double> forces_local(3 * local_ions_.size());
-
     getLocalForces(forces_local);
 
-    int n = getNumIons();
-    forces.resize(3 * n);
-
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
     mmpi.allGatherV(forces_local, forces);
 }
diff --git a/src/Ions.h b/src/Ions.h
index 3420d725..3d382eb7 100644
--- a/src/Ions.h
+++ b/src/Ions.h
@@ -292,10 +292,10 @@ class Ions
     void getLocalPositions(std::vector<double>& tau) const;
     void getLocalNames(std::vector<std::string>& names) const;
     void getNames(std::vector<std::string>& names) const;
-    void getPositions(std::vector<double>& tau);
-    void getAtomicNumbers(std::vector<short>& atnumbers);
+    void getPositions(std::vector<double>& tau) const;
+    void getAtomicNumbers(std::vector<short>& atnumbers) const;
 
-    void getForces(std::vector<double>& forces);
+    void getForces(std::vector<double>& forces) const;
     void getLocalForces(std::vector<double>& tau) const;
 
     /*!
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 0b4dfc53..02edbe8b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -557,14 +557,16 @@ add_test(NAME testMD_D72
          ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/lrs.in
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
-add_test(NAME testHDF5single
-         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/test.py
-         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
-         ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
-         ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/mgmol.cfg
-         ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/md.cfg
-         ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/h2o.xyz
-         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+if(MGMOL_USE_HDF5P)
+  add_test(NAME testHDF5single
+           COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/test.py
+           ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+           ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
+           ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/mgmol.cfg
+           ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/md.cfg
+           ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/h2o.xyz
+           ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+endif()
 add_test(NAME testMD_MVP
          COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/test.py
          ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
diff --git a/tests/MVP/mvp.cfg b/tests/MVP/mvp.cfg
index bb3ab52a..868e5703 100644
--- a/tests/MVP/mvp.cfg
+++ b/tests/MVP/mvp.cfg
@@ -32,4 +32,3 @@ solver=MVP
 nb_inner_it=2
 [Restart]
 output_level=2
-output_type=single_file
diff --git a/tests/RhoVhRestart/testRhoVhRestart.cc b/tests/RhoVhRestart/testRhoVhRestart.cc
index 09da3320..ce40272d 100644
--- a/tests/RhoVhRestart/testRhoVhRestart.cc
+++ b/tests/RhoVhRestart/testRhoVhRestart.cc
@@ -46,8 +46,8 @@ int testRhoRestart(MGmolInterface* mgmol_)
     /* check if the recomputed density is the same */
     for (int d = 0; d < (int)rho0.size(); d++)
     {
-        double error = abs(rho0[d] - rho->rho_[0][d]) / abs(rho0[d]);
-        if (error > 1e-10)
+        double error = abs(rho0[d] - rho->rho_[0][d]);
+        if (error > 1e-10 * abs(rho0[d]))
         {
             printf("rank %d, rho[%d]=%.15e, rho0[%d]=%.15e\n", rank, d,
                 rho->rho_[0][d], d, rho0[d]);
diff --git a/tests/testIons.cc b/tests/testIons.cc
index 39c96f4a..1741e0f7 100644
--- a/tests/testIons.cc
+++ b/tests/testIons.cc
@@ -107,6 +107,7 @@ int main(int argc, char** argv)
         ions.getAtomicNumbers(anumbers);
         if (myrank == 0)
         {
+            std::cout << "Positions:" << std::endl;
             int i = 0;
             for (auto& position : positions)
             {
@@ -121,7 +122,7 @@ int main(int argc, char** argv)
         MPI_Barrier(MPI_COMM_WORLD);
 
         // swap x and z
-        for (size_t i = 0; i < positions.size() - 2; i++)
+        for (size_t i = 0; i < positions.size() - 2; i += 3)
         {
             double x         = positions[i];
             double z         = positions[i + 2];
@@ -162,10 +163,13 @@ int main(int argc, char** argv)
     MPI_Barrier(MPI_COMM_WORLD);
 
     std::vector<double> forces(3 * na);
-    // arbitrary value
-    const double fval = 1.12;
+    // set forces to a different arbitrary value for each component
+    int i = 0;
     for (auto& f : forces)
-        f = fval;
+    {
+        f = (double)i;
+        i++;
+    }
 
     ions.getNames(names);
     ions.setLocalForces(forces, names);
@@ -177,13 +181,29 @@ int main(int argc, char** argv)
     ions.getLocalForces(lforces);
     for (auto& f : lforces)
     {
-        if (std::abs(f - fval) > 1.e-14)
+        if (std::fmod(f, 1.) > 1.e-14)
         {
             std::cerr << "f = " << f << std::endl;
             MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
         }
     }
 
+    ions.getForces(forces);
+    if (myrank == 0)
+        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) < 1.e-14)
+                {
+                    std::cerr << "f0 = " << *f0 << ", f1 = " << *f1
+                              << std::endl;
+                    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
+                }
+            }
+        }
     mpirc = MPI_Finalize();
     if (mpirc != MPI_SUCCESS)
     {