diff --git a/src/ABPG.cc b/src/ABPG.cc
index 62cb6ba1..107edc63 100644
--- a/src/ABPG.cc
+++ b/src/ABPG.cc
@@ -38,7 +38,7 @@ void ABPG<T>::setup(T& orbitals)
 //
 // orthof=true: wants orthonormalized updated wave functions
 template <class T>
-int ABPG<T>::updateWF(T& orbitals, Ions& /*ions*/, const double precond_factor,
+int ABPG<T>::updateWF(T& orbitals, Ions& ions, const double precond_factor,
     const bool /*orthof*/, T& work_orbitals, const bool accelerate,
     const bool print_res, const double atol)
 {
@@ -51,8 +51,8 @@ int ABPG<T>::updateWF(T& orbitals, Ions& /*ions*/, const double precond_factor,
     T res("Residual", orbitals, false);
 
     const bool check_res = (atol > 0.);
-    double normRes       = mgmol_strategy_->computeResidual(
-        orbitals, work_orbitals, res, (print_res || check_res), check_res);
+    double normRes = mgmol_strategy_->computeResidual(orbitals, work_orbitals,
+        ions, res, (print_res || check_res), check_res);
     if (normRes < atol && check_res)
     {
         abpg_nl_update_tm_.stop();
diff --git a/src/DFTsolver.cc b/src/DFTsolver.cc
index 91c4a0f6..206fba91 100644
--- a/src/DFTsolver.cc
+++ b/src/DFTsolver.cc
@@ -196,8 +196,8 @@ double DFTsolver<OrbitalsType>::evaluateEnergy(
 
     // Get the new total energy
     const double ts = 0.5 * proj_matrices_->computeEntropy(); // in [Ha]
-    eks_history_[0]
-        = energy_->evaluateTotal(ts, proj_matrices_, orbitals, print_flag, os_);
+    eks_history_[0] = energy_->evaluateTotal(
+        ts, proj_matrices_, ions_, orbitals, print_flag, os_);
 
     sum_eig_[1] = sum_eig_[0];
     sum_eig_[0] = 2. * proj_matrices_->getEigSum(); // 2.*sum in [Ry]
diff --git a/src/DavidsonSolver.cc b/src/DavidsonSolver.cc
index 0dc88ef4..a049fde7 100644
--- a/src/DavidsonSolver.cc
+++ b/src/DavidsonSolver.cc
@@ -500,7 +500,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 ts0 = evalEntropy(projmatrices, (ct.verbose > 1), os_);
                 e0  = energy_->evaluateTotal(
-                    ts0, projmatrices, orbitals, printE, os_);
+                    ts0, projmatrices, ions_, orbitals, printE, os_);
 
                 retval = checkConvergence(e0, outer_it, ct.conv_tol);
                 if (retval == 0 || (outer_it == ct.max_electronic_steps))
@@ -549,7 +549,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 ts0 = evalEntropy(proj_mat2N_.get(), (ct.verbose > 1), os_);
                 e0  = energy_->evaluateTotal(
-                    ts0, proj_mat2N_.get(), orbitals, printE, os_);
+                    ts0, proj_mat2N_.get(), ions_, orbitals, printE, os_);
             }
 
             // 2N x 2N target...
@@ -621,8 +621,8 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 const double ts1
                     = evalEntropy(proj_mat2N_.get(), (ct.verbose > 2), os_);
-                const double e1 = energy_->evaluateTotal(
-                    ts1, proj_mat2N_.get(), orbitals, ct.verbose - 1, os_);
+                const double e1 = energy_->evaluateTotal(ts1, proj_mat2N_.get(),
+                    ions_, orbitals, ct.verbose - 1, os_);
 
                 // line minimization
                 beta = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);
diff --git a/src/Energy.cc b/src/Energy.cc
index 1cee6d67..f53dbc93 100644
--- a/src/Energy.cc
+++ b/src/Energy.cc
@@ -25,11 +25,10 @@
 #define RY2HA 0.5
 
 template <class T>
-Energy<T>::Energy(const pb::Grid& mygrid, const Ions& ions,
-    const Potentials& pot, const Electrostatic& es, const Rho<T>& rho,
-    const XConGrid& xc, SpreadPenaltyInterface<T>* spread_penalty)
+Energy<T>::Energy(const pb::Grid& mygrid, const Potentials& pot,
+    const Electrostatic& es, const Rho<T>& rho, const XConGrid& xc,
+    SpreadPenaltyInterface<T>* spread_penalty)
     : mygrid_(mygrid),
-      ions_(ions),
       pot_(pot),
       es_(es),
       rho_(rho),
@@ -59,7 +58,7 @@ double Energy<T>::getEVrhoRho() const
 }
 
 template <class T>
-double Energy<T>::evaluateEnergyIonsInVext()
+double Energy<T>::evaluateEnergyIonsInVext(Ions& ions)
 {
     double energy = 0.;
 
@@ -69,12 +68,12 @@ double Energy<T>::evaluateEnergyIonsInVext()
     //(*MPIdata::sout)<<"Energy<T>::evaluateEnergyIonsInVext()"<<std::endl;
     double position[3];
     std::vector<double> positions;
-    positions.reserve(3 * ions_.local_ions().size());
+    positions.reserve(3 * ions.local_ions().size());
 
     // loop over ions
     int nions                             = 0;
-    std::vector<Ion*>::const_iterator ion = ions_.local_ions().begin();
-    while (ion != ions_.local_ions().end())
+    std::vector<Ion*>::const_iterator ion = ions.local_ions().begin();
+    while (ion != ions.local_ions().end())
     {
         (*ion)->getPosition(position);
         positions.push_back(position[0]);
@@ -88,9 +87,9 @@ double Energy<T>::evaluateEnergyIonsInVext()
     pot_.getValVext(positions, val);
 
     // loop over ions again
-    ion           = ions_.local_ions().begin();
+    ion           = ions.local_ions().begin();
     int ion_index = 0;
-    while (ion != ions_.local_ions().end())
+    while (ion != ions.local_ions().end())
     {
         const double z = (*ion)->getZion();
         // int ion_index=(*ion)->index();
@@ -112,16 +111,16 @@ double Energy<T>::evaluateEnergyIonsInVext()
 
 template <class T>
 double Energy<T>::evaluateTotal(const double ts, // in [Ha]
-    ProjectedMatricesInterface* projmatrices, const T& phi, const int verbosity,
-    std::ostream& os)
+    ProjectedMatricesInterface* projmatrices, Ions& ions, const T& phi,
+    const int verbosity, std::ostream& os)
 {
     eval_te_tm_.start();
 
     Control& ct = *(Control::instance());
 
-    const double eself = ions_.energySelf();
-    const double ediff = ions_.energyDiff(ct.bcPoisson);
-    const double eipot = evaluateEnergyIonsInVext();
+    const double eself = ions.energySelf();
+    const double ediff = ions.energyDiff(ct.bcPoisson);
+    const double eipot = evaluateEnergyIonsInVext(ions);
 
     const double eigsum = 0.5 * projmatrices->getExpectationH();
 
diff --git a/src/Energy.h b/src/Energy.h
index e6cafd98..d55d4040 100644
--- a/src/Energy.h
+++ b/src/Energy.h
@@ -11,6 +11,7 @@
 #define MGMOL_ENERGY_H
 
 #include "Grid.h"
+#include "Ions.h"
 #include "Rho.h"
 #include "SpreadPenaltyInterface.h"
 #include "Timer.h"
@@ -20,7 +21,6 @@
 #include <vector>
 
 class Potentials;
-class Ions;
 class Electrostatic;
 class ProjectedMatricesInterface;
 class XConGrid;
@@ -29,7 +29,6 @@ template <class T>
 class Energy
 {
     const pb::Grid& mygrid_;
-    const Ions& ions_;
     const Potentials& pot_;
     const Electrostatic& es_;
     const Rho<T>& rho_;
@@ -45,16 +44,15 @@ class Energy
     double getEVrhoRho() const;
 
 public:
-    Energy(const pb::Grid&, const Ions&, const Potentials&,
-        const Electrostatic&, const Rho<T>&, const XConGrid&,
-        SpreadPenaltyInterface<T>*);
+    Energy(const pb::Grid&, const Potentials&, const Electrostatic&,
+        const Rho<T>&, const XConGrid&, SpreadPenaltyInterface<T>*);
 
     static Timer eval_te_tm() { return eval_te_tm_; }
 
-    double evaluateTotal(const double ts, ProjectedMatricesInterface*,
+    double evaluateTotal(const double ts, ProjectedMatricesInterface*, Ions&,
         const T& phi, const int, std::ostream&);
 
-    double evaluateEnergyIonsInVext();
+    double evaluateEnergyIonsInVext(Ions&);
 
     void saveVofRho();
 };
diff --git a/src/GrassmanLineMinimization.cc b/src/GrassmanLineMinimization.cc
index a83d6b08..f642e0f4 100644
--- a/src/GrassmanLineMinimization.cc
+++ b/src/GrassmanLineMinimization.cc
@@ -29,7 +29,7 @@ Timer GrassmanLineMinimization<T>::update_states_tm_("Grassman_update_states");
 //
 // orthof=true: wants orthonormalized updated wave functions
 template <class T>
-int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& /*ions*/,
+int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& ions,
     const double precond_factor, const bool orthof, T& work_orbitals,
     const bool accelerate, const bool print_res, const double atol)
 {
@@ -61,7 +61,7 @@ int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& /*ions*/,
     // Update wavefunctions
     const bool check_res = (atol > 0.);
     double normRes = mgmol_strategy_->computeResidual(orbitals, work_orbitals,
-        *new_grad_, (print_res || check_res), check_res);
+        ions, *new_grad_, (print_res || check_res), check_res);
     if (normRes < atol && check_res)
     {
         nl_update_tm_.stop();
diff --git a/src/HamiltonianMVPSolver.cc b/src/HamiltonianMVPSolver.cc
index 4bc017dc..f90dd228 100644
--- a/src/HamiltonianMVPSolver.cc
+++ b/src/HamiltonianMVPSolver.cc
@@ -160,8 +160,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
         // compute energy at origin
         const int printE = (ct.verbose > 1) ? 1 : 0;
-        double e0
-            = energy_->evaluateTotal(ts0, projmatrices, orbitals, printE, os_);
+        double e0        = energy_->evaluateTotal(
+            ts0, projmatrices, ions_, orbitals, printE, os_);
 
         //
         // compute energy at end for new H
@@ -189,8 +189,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
         projmatrices->setHB2H();
 
         // compute energy at end (beta=1.)
-        double e1
-            = energy_->evaluateTotal(ts1, projmatrices, orbitals, printE, os_);
+        double e1 = energy_->evaluateTotal(
+            ts1, projmatrices, ions_, orbitals, printE, os_);
 
         //
         // evaluate energy at mid-point
@@ -226,8 +226,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
         projmatrices->setHB2H();
 
         // compute energy at midpoint
-        double ei
-            = energy_->evaluateTotal(tsi, projmatrices, orbitals, printE, os_);
+        double ei = energy_->evaluateTotal(
+            tsi, projmatrices, ions_, orbitals, printE, os_);
 
         // line minimization
         double beta
@@ -285,7 +285,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
                 // compute energy at end (beta=1.)
                 ei = energy_->evaluateTotal(
-                    tsi, projmatrices, orbitals, printE, os_);
+                    tsi, projmatrices, ions_, orbitals, printE, os_);
 
                 // line minimization
                 beta = minQuadPolynomialFrom3values(
diff --git a/src/Ions.cc b/src/Ions.cc
index 3b128b0e..ef9d687c 100644
--- a/src/Ions.cc
+++ b/src/Ions.cc
@@ -70,8 +70,18 @@ void writeData2d(HDFrestart& h5f_file, std::string datasetname,
     }
 }
 
-Ions::Ions(const double lat[3], const std::vector<Species>& sp) : species_(sp)
+Ions::Ions(const double lat[3], const std::vector<Species>& sp)
+    : species_(sp), setup_(false), has_locked_atoms_(false)
 {
+    setupSubdomains(lat);
+}
+
+void Ions::setupSubdomains(const double lat[3])
+{
+    MGmol_MPI& mmpi          = *(MGmol_MPI::instance());
+    Mesh* mymesh             = Mesh::instance();
+    const pb::PEenv& myPEenv = mymesh->peenv();
+
     for (short i = 0; i < 3; i++)
     {
         assert(lat[i] > 0.);
@@ -79,15 +89,10 @@ Ions::Ions(const double lat[3], const std::vector<Species>& sp) : species_(sp)
         {
             (*MPIdata::serr) << "Ions constructor: lattice[" << i
                              << "]=" << lat[i] << "!!!" << std::endl;
-            exit(2);
+            mmpi.abort();
         }
         lattice_[i] = lat[i];
     }
-    setup_            = false;
-    has_locked_atoms_ = false;
-
-    Mesh* mymesh             = Mesh::instance();
-    const pb::PEenv& myPEenv = mymesh->peenv();
 
     for (short i = 0; i < 3; i++)
         div_lattice_[i] = lattice_[i] / (double)(myPEenv.n_mpi_task(i));
@@ -107,6 +112,20 @@ Ions::Ions(const double lat[3], const std::vector<Species>& sp) : species_(sp)
         MPI_Cart_shift(cart_comm_, dir, disp, &source_[dir], &dest_[dir]);
 }
 
+Ions::Ions(const std::vector<double>& p, const std::vector<short>& anum,
+    const double lat[3], const std::vector<Species>& sp)
+    : species_(sp)
+{
+    setupSubdomains(lat);
+
+    double rmax = getMaxListRadius();
+    setupListIonsBoundaries(rmax);
+
+    num_ions_ = setAtoms(p, anum);
+
+    setup();
+}
+
 Ions::Ions(const Ions& ions, const double shift[3]) : species_(ions.species_)
 {
     for (const auto& ion : ions.list_ions_)
@@ -181,8 +200,6 @@ void Ions::setup()
 
     ions_setup_tm.start();
 
-    Ion::resetIndexCount();
-
     updateListIons();
 
     //#ifndef NDEBUG
@@ -197,7 +214,7 @@ void Ions::setup()
     has_locked_atoms_ = hasLockedAtoms();
 
     // initialize data for constraints
-    setupContraintsData(interacting_ions_);
+    setupContraintsData();
 
     computeMaxNumProjs();
 
@@ -297,13 +314,13 @@ void Ions::setupInteractingIons()
 
 // setup arrays to be used in constraints enforcement
 // using references to local_ions and extra "dummy" data
-void Ions::setupContraintsData(std::vector<Ion*>& ions_for_constraints)
+void Ions::setupContraintsData()
 {
     Control& ct = *(Control::instance());
 
     if (ct.verbose > 0)
         printWithTimeStamp("Ions::setupContraintsData()...", std::cout);
-    const int nnloc = ions_for_constraints.size() - local_ions_.size();
+    const int nnloc = interacting_ions_.size() - local_ions_.size();
     // std::cout<<"interacting_ions_.size()="<<interacting_ions_.size()<<endl;
     // std::cout<<"local_ions_.size()="<<local_ions_.size()<<endl;
     assert(nnloc >= 0);
@@ -322,12 +339,11 @@ void Ions::setupContraintsData(std::vector<Ion*>& ions_for_constraints)
     interacting_atmove_.clear();
     interacting_pmass_.clear();
 
-    int ia                                = 0; // count local ions
-    int ib                                = 0; // count non-local ions
-    std::vector<Ion*>::const_iterator ion = ions_for_constraints.begin();
-    while (ion != ions_for_constraints.end())
+    int ia = 0; // count local ions
+    int ib = 0; // count non-local ions
+    for (auto& ion : interacting_ions_)
     {
-        if (isLocal((*ion)->name()))
+        if (isLocal(ion->name()))
         {
             assert(3 * ia < static_cast<int>(tau0_.size()));
             assert(3 * ia < static_cast<int>(taup_.size()));
@@ -355,12 +371,12 @@ void Ions::setupContraintsData(std::vector<Ion*>& ions_for_constraints)
             assert(3 * ib < static_cast<int>(fion_dummy_.size()));
 
             // fill up dummy data
-            (*ion)->getPosition(&tau0_dummy_[3 * ib]);
-            (*ion)->getPosition(&taup_dummy_[3 * ib]);
-            (*ion)->getForce(&fion_dummy_[3 * ib]);
-            names_dummy_.push_back((*ion)->name());
-            pmass_dummy_.push_back((*ion)->getMass());
-            atmove_dummy_.push_back(!(*ion)->locked());
+            ion->getPosition(&tau0_dummy_[3 * ib]);
+            ion->getPosition(&taup_dummy_[3 * ib]);
+            ion->getForce(&fion_dummy_[3 * ib]);
+            names_dummy_.push_back(ion->name());
+            pmass_dummy_.push_back(ion->getMass());
+            atmove_dummy_.push_back(!ion->locked());
 
             // set references to dummy data
             interacting_tau0_.push_back(&tau0_dummy_[3 * ib]);
@@ -378,8 +394,6 @@ void Ions::setupContraintsData(std::vector<Ion*>& ions_for_constraints)
 
             ++ib;
         }
-
-        ++ion;
     }
 }
 
@@ -1635,14 +1649,15 @@ Ion* Ions::findLocalIon(const int index) const
     return nullptr;
 }
 
-void Ions::setLocalPositions(const std::vector<double>& tau)
+void Ions::setLocalPositions(const std::vector<double>& positions)
 {
-    assert(tau.size() == 3 * local_ions_.size());
+    assert(positions.size() == 3 * local_ions_.size());
 
     int ia = 0;
     for (auto& ion : local_ions_)
     {
-        ion->setPosition(tau[3 * ia + 0], tau[3 * ia + 1], tau[3 * ia + 2]);
+        ion->setPosition(positions[3 * ia + 0], positions[3 * ia + 1],
+            positions[3 * ia + 2]);
         ia++;
     }
 
@@ -1799,6 +1814,8 @@ int Ions::readAtomsFromXYZ(
 int Ions::setAtoms(
     const std::vector<double>& crds, const std::vector<short>& anum)
 {
+    // we are setting up a new list of atoms and need an index starting at 0
+    Ion::resetIndexCount();
     const int natoms = crds.size() / 3;
 
     double velocity[3] = { 0., 0., 0. };
@@ -1858,6 +1875,7 @@ void Ions::addIonToList(const Species& sp, const std::string& name,
 
     // create a new Ion
     Ion* new_ion = new Ion(sp, name, &crds[0], velocity, locked);
+    // std::cout << "New Ion index = " << new_ion->index() << std::endl;
     new_ion->bcast(mmpi.commGlobal());
 
     if (inListIons(crds[0], crds[1], crds[2]))
@@ -2195,14 +2213,14 @@ void Ions::setVelocities(const std::vector<double>& tau0,
     }
 }
 
-void Ions::getLocalPositions(std::vector<double>& tau) const
+void Ions::getLocalPositions(std::vector<double>& pos) const
 {
-    assert(tau.size() == 3 * local_ions_.size());
+    assert(pos.size() == 3 * local_ions_.size());
 
     int ia = 0;
     for (auto& ion : local_ions_)
     {
-        ion->getPosition(&tau[3 * ia]);
+        ion->getPosition(&pos[3 * ia]);
         ia++;
     }
 }
@@ -2225,14 +2243,14 @@ void Ions::getNames(std::vector<std::string>& names) const
     mmpi.allGatherV(local_names, names);
 }
 
-void Ions::getPositions(std::vector<double>& tau) const
+void Ions::getPositions(std::vector<double>& pos) const
 {
-    std::vector<double> tau_local(3 * local_ions_.size());
+    std::vector<double> pos_local(3 * local_ions_.size());
 
-    getLocalPositions(tau_local);
+    getLocalPositions(pos_local);
 
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
-    mmpi.allGatherV(tau_local, tau);
+    mmpi.allGatherV(pos_local, pos);
 }
 
 void Ions::getAtomicNumbers(std::vector<short>& atnumbers) const
@@ -2282,14 +2300,14 @@ void Ions::setPositionsToTau0()
 }
 
 void Ions::setPositions(
-    const std::vector<double>& tau, const std::vector<short>& anumbers)
+    const std::vector<double>& p, const std::vector<short>& anumbers)
 {
-    assert(tau.size() == anumbers.size() * 3);
+    assert(p.size() == anumbers.size() * 3);
 
     // clear previous data
     clearLists();
 
-    num_ions_ = setAtoms(tau, anumbers);
+    num_ions_ = setAtoms(p, anumbers);
 
     // setup required after updating local ions positions
     setup();
@@ -2308,15 +2326,15 @@ void Ions::setVelocitiesToVel()
     }
 }
 
-void Ions::getLocalForces(std::vector<double>& tau) const
+void Ions::getLocalForces(std::vector<double>& f) const
 {
-    assert(tau.size() == 3 * local_ions_.size());
+    assert(f.size() == 3 * local_ions_.size());
 
     int ia = 0;
-    for (auto& ion : local_ions_)
+    for (const auto& ion : local_ions_)
     {
-        assert(3 * ia + 2 < (int)tau.size());
-        ion->getForce(&tau[3 * ia]);
+        assert(3 * ia + 2 < (int)f.size());
+        ion->getForce(&f[3 * ia]);
         ia++;
     }
 }
@@ -2713,6 +2731,8 @@ void Ions::gatherForces(std::vector<double>& forces, const int root) const
     for (auto& ion : local_ions_)
     {
         const int index = ion->index();
+        // std::cout << "index = " << index << std::endl;
+        assert(index < num_ions_);
         assert(forces.size() >= 3 * index);
         assert(index < num_ions_);
         ion->getForce(&forces[3 * index]);
@@ -3107,80 +3127,48 @@ void Ions::augmentIonsData(const int nsteps, const int dir, const int disp,
 
 void Ions::updateForcesInteractingIons()
 {
-    MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
-
     // get computed forces into fion_
     getLocalForces(fion_);
 
-    // initialize with local names and forces
-    DistributedIonicData forces_data(local_names_, fion_);
-
-    for (short dir = 0; dir < 3; dir++)
-    {
-        // initial data
-        const int lsize = forces_data.size();
-        int maxsize     = lsize;
-        mmpi.allreduce(&maxsize, 1, MPI_MAX);
-
-        // send local data
-        DistributedIonicData data2send(forces_data);
-
-        // send right to left
-        int disp = -1;
-        forces_data.augmentData(
-            lstep_[dir], dir, disp, lsize, maxsize, data2send);
-
-        // send left to right
-        disp = 1;
-        forces_data.augmentData(
-            rstep_[dir], dir, disp, lsize, maxsize, data2send);
-    }
-
-    int ia = 0;
-    for (std::vector<std::string>::const_iterator it
-         = interacting_names_.begin();
-         it != interacting_names_.end(); ++it)
-    {
-        double* force = interacting_fion_[ia];
-        forces_data.getData(*it, force);
-        ++ia;
-    }
+    updateDataInteractingIons(fion_, interacting_fion_);
 }
 
 void Ions::updateTaupInteractingIons()
+{
+    updateDataInteractingIons(taup_, interacting_taup_);
+}
+
+void Ions::updateDataInteractingIons(
+    std::vector<double>& data, std::vector<double*>& interacting_data)
 {
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
 
     // initialize with local names and forces
-    DistributedIonicData taup_data(local_names_, taup_);
+    DistributedIonicData ddata(local_names_, data);
 
     for (short dir = 0; dir < 3; dir++)
     {
         // initial data
-        const int lsize = taup_data.size();
+        const int lsize = ddata.size();
         int maxsize     = lsize;
         mmpi.allreduce(&maxsize, 1, MPI_MAX);
 
-        DistributedIonicData data2send(taup_data);
+        DistributedIonicData data2send(ddata);
 
         // send right to left
         int disp = -1;
-        taup_data.augmentData(
-            lstep_[dir], dir, disp, lsize, maxsize, data2send);
+        ddata.augmentData(lstep_[dir], dir, disp, lsize, maxsize, data2send);
 
         // send left to right
         disp = 1;
-        taup_data.augmentData(
-            rstep_[dir], dir, disp, lsize, maxsize, data2send);
+        ddata.augmentData(rstep_[dir], dir, disp, lsize, maxsize, data2send);
     }
 
     int ia = 0;
-    for (std::vector<std::string>::const_iterator it
-         = interacting_names_.begin();
-         it != interacting_names_.end(); ++it)
+    for (const auto& name : interacting_names_)
     {
-        double* taup = interacting_taup_[ia];
-        taup_data.getData(*it, taup);
+        double* ptr = interacting_data[ia];
+        ddata.getData(name, ptr);
         ++ia;
     }
 }
@@ -3203,12 +3191,11 @@ void Ions::updateListIons()
     // collect local_ions data
     // assume local_ions size is same as size of ions_names std::vector
     ions_data_.clear();
-
     {
         IonData idata;
-        for (auto& lion : local_ions_)
+        for (const auto& ion : local_ions_)
         {
-            lion->getIonData(idata);
+            ion->getIonData(idata);
 
             // populate ions_data_ list
             ions_data_.push_back(idata);
@@ -3307,20 +3294,20 @@ void Ions::initStepperData()
 {
     clearStepperData();
 
-    for (auto& lion : local_ions_)
+    for (const auto& ion : local_ions_)
     {
-        local_names_.push_back(lion->name());
-        atmove_.push_back(!lion->locked());
-        pmass_.push_back(lion->getMass());
-        gids_.push_back(lion->index());
+        local_names_.push_back(ion->name());
+        atmove_.push_back(!ion->locked());
+        pmass_.push_back(ion->getMass());
+        gids_.push_back(ion->index());
 
         for (short i = 0; i < 3; i++)
         {
-            taum_.push_back(lion->getPreviousPosition(i));
-            tau0_.push_back(lion->position(i));
-            fion_.push_back(lion->force(i));
-            velocity_.push_back(lion->velocity(i));
-            rand_states_.push_back(lion->randomState(i));
+            taum_.push_back(ion->getPreviousPosition(i));
+            tau0_.push_back(ion->position(i));
+            fion_.push_back(ion->force(i));
+            velocity_.push_back(ion->velocity(i));
+            rand_states_.push_back(ion->randomState(i));
         }
     }
     // initialize taup to enable computing velocities
diff --git a/src/Ions.h b/src/Ions.h
index 3d382eb7..f39e1fed 100644
--- a/src/Ions.h
+++ b/src/Ions.h
@@ -35,6 +35,10 @@ class Ions
 
     const std::vector<Species>& species_;
 
+    bool setup_;
+
+    bool has_locked_atoms_;
+
     std::vector<Ion*> list_ions_;
 
     /*!
@@ -59,8 +63,6 @@ class Ions
     int dest_[3];
     MPI_Comm cart_comm_; // MPI cartesian communicator for data distribution
 
-    bool has_locked_atoms_;
-
     /*!
      * Prevent usage of copy constructor by making it private and
      * non-implemented
@@ -73,6 +75,7 @@ class Ions
     void readRestartPositions(HDFrestart& h5_file);
     int read1atom(std::ifstream* tfile, const bool cell_relative);
 
+    void setupSubdomains(const double lat[3]);
     void setupInteractingIons();
     void setupListOverlappingIons();
     void setMapVL();
@@ -129,14 +132,12 @@ class Ions
     int rstep_[3]; // number of steps to the right for each dimension to gather
                    // data
 
-    bool setup_;
-
     void gatherLockedData(std::vector<int>& locked_data, const int root) const;
     void computeNumIons(void);
     int readNatoms(const std::string& input_file, const bool cell_relative);
     int readNatoms(std::ifstream* tfile, const bool cell_relative);
     int readAtomsFromXYZ(const std::string& filename, const bool cell_relative);
-    void setupContraintsData(std::vector<Ion*>&);
+    void setupContraintsData();
     void clearStepperData();
     void initStepperData();
     void computeMaxNumProjs();
@@ -180,6 +181,9 @@ class Ions
 
     Ions(const Ions&, const double shift[3]);
 
+    Ions(const std::vector<double>& p, const std::vector<short>& anum,
+        const double lat[3], const std::vector<Species>& sp);
+
     ~Ions();
 
     void setup();
@@ -357,9 +361,14 @@ class Ions
         return *spi;
     }
 
+    const std::vector<Species>& getSpecies() const { return species_; }
+
     void updateForcesInteractingIons();
     void updateTaupInteractingIons();
 
+    void updateDataInteractingIons(
+        std::vector<double>& data, std::vector<double*>& interacting_data);
+
     /*!
      * Calculate minimum distance between local pairs
      */
diff --git a/src/MGmol.cc b/src/MGmol.cc
index 4d69387d..018a07a2 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -455,9 +455,8 @@ int MGmol<OrbitalsType>::initial()
 
     std::shared_ptr<SpreadPenaltyInterface<OrbitalsType>> spread_penalty
         = energy_with_spread_penalty ? spread_penalty_ : nullptr;
-    energy_ = std::shared_ptr<Energy<OrbitalsType>>(
-        new Energy<OrbitalsType>(mygrid, *ions_, pot, *electrostat_, *rho_,
-            *xcongrid_, spread_penalty.get()));
+    energy_ = std::shared_ptr<Energy<OrbitalsType>>(new Energy<OrbitalsType>(
+        mygrid, pot, *electrostat_, *rho_, *xcongrid_, spread_penalty.get()));
 
     if (ct.verbose > 0) printWithTimeStamp("Setup matrices...", os_);
 
@@ -568,7 +567,7 @@ void MGmol<OrbitalsType>::finalEnergy()
     // Get the total energy
     const double ts = 0.5 * proj_matrices_->computeEntropy(); // in [Ha]
     total_energy_   = energy_->evaluateTotal(
-        ts, proj_matrices_.get(), *current_orbitals_, 2, os_);
+        ts, proj_matrices_.get(), *ions_, *current_orbitals_, 2, os_);
 }
 
 template <class OrbitalsType>
@@ -828,7 +827,7 @@ void MGmol<OrbitalsType>::setupPotentials(Ions& ions)
     pot.initialize(ions);
 
     if (ct.verbose > 0) printWithTimeStamp("Setup kbpsi...", os_);
-    g_kbpsi_->setup(*ions_);
+    g_kbpsi_->setup(ions);
 
     electrostat_->setupRhoc(pot.rho_comp());
 
@@ -1170,8 +1169,8 @@ void MGmol<OrbitalsType>::precond_mg(OrbitalsType& phi)
 
 template <class OrbitalsType>
 double MGmol<OrbitalsType>::computeResidual(OrbitalsType& orbitals,
-    OrbitalsType& work_orbitals, OrbitalsType& res, const bool print_residual,
-    const bool norm_res)
+    OrbitalsType& work_orbitals, Ions& ions, OrbitalsType& res,
+    const bool print_residual, const bool norm_res)
 
 {
     assert(orbitals.getIterativeIndex() >= 0);
@@ -1190,7 +1189,7 @@ double MGmol<OrbitalsType>::computeResidual(OrbitalsType& orbitals,
 
     // get H*psi stored in work_orbitals.psi
     // and psi^T H psi in Hij
-    getHpsiAndTheta(*ions_, orbitals, work_orbitals);
+    getHpsiAndTheta(ions, orbitals, work_orbitals);
 
     double norm2Res = computeConstraintResidual(
         orbitals, work_orbitals, res, print_residual, norm_res);
@@ -1428,17 +1427,22 @@ double MGmol<OrbitalsType>::evaluateEnergyAndForces(Orbitals* orbitals,
 
     Control& ct = *(Control::instance());
 
-    ions_->setPositions(tau, atnumbers);
+    // create a new temporary Ions object to be used for
+    // energy end forces calculation
+    Mesh* mymesh            = Mesh::instance();
+    const pb::Grid& mygrid  = mymesh->grid();
+    const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) };
+    Ions ions(tau, atnumbers, lattice, ions_->getSpecies());
 
-    setupPotentials(*ions_);
+    setupPotentials(ions);
 
     double eks              = 0.;
     OrbitalsType* dorbitals = dynamic_cast<OrbitalsType*>(orbitals);
-    quench(*dorbitals, *ions_, ct.max_electronic_steps, 20, eks);
+    quench(*dorbitals, ions, ct.max_electronic_steps, 20, eks);
 
-    force(*dorbitals, *ions_);
+    force(*dorbitals, ions);
 
-    ions_->getForces(forces);
+    ions.getForces(forces);
 
     return eks;
 }
@@ -1450,37 +1454,42 @@ double MGmol<OrbitalsType>::evaluateDMandEnergyAndForces(Orbitals* orbitals,
 {
     OrbitalsType* dorbitals = dynamic_cast<OrbitalsType*>(orbitals);
 
-    ions_->setPositions(tau, atnumbers);
+    // create a new temporary Ions object to be used for
+    // energy end forces calculation
+    Mesh* mymesh            = Mesh::instance();
+    const pb::Grid& mygrid  = mymesh->grid();
+    const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) };
+    Ions ions(tau, atnumbers, lattice, ions_->getSpecies());
 
-    setupPotentials(*ions_);
+    setupPotentials(ions);
 
     // initialize electronic density
     rho_->update(*dorbitals);
 
     // initialize potential
-    update_pot(*ions_);
+    update_pot(ions);
 
     // initialize projected matrices
-    updateHmatrix(*dorbitals, *ions_);
+    updateHmatrix(*dorbitals, ions);
     proj_matrices_->updateThetaAndHB();
 
     // compute DM
     std::shared_ptr<DMStrategy<OrbitalsType>> dm_strategy(
         DMStrategyFactory<OrbitalsType,
-            dist_matrix::DistMatrix<double>>::create(comm_, os_, *ions_,
+            dist_matrix::DistMatrix<double>>::create(comm_, os_, ions,
             rho_.get(), energy_.get(), electrostat_.get(), this,
             proj_matrices_.get(), dorbitals));
 
     dm_strategy->update(*dorbitals);
 
     // evaluate energy and forces
-    double ts = 0.;
-    double eks
-        = energy_->evaluateTotal(ts, proj_matrices_.get(), *dorbitals, 2, os_);
+    double ts  = 0.;
+    double eks = energy_->evaluateTotal(
+        ts, proj_matrices_.get(), ions, *dorbitals, 2, os_);
 
-    force(*dorbitals, *ions_);
+    force(*dorbitals, ions);
 
-    ions_->getForces(forces);
+    ions.getForces(forces);
 
     return eks;
 }
diff --git a/src/MGmol.h b/src/MGmol.h
index 12197dfd..c85e7b59 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -335,17 +335,13 @@ class MGmol : public MGmolInterface
     void precond_mg(OrbitalsType& orbitals);
     void setGamma(const pb::Lap<ORBDTYPE>& lapOper, const Potentials& pot);
     double computeResidual(OrbitalsType& orbitals, OrbitalsType& work_orbitals,
-        OrbitalsType& res, const bool print_residual, const bool norm_res);
+        Ions& ions, OrbitalsType& res, const bool print_residual,
+        const bool norm_res);
     void applyAOMMprojection(OrbitalsType&);
     void force(OrbitalsType& orbitals, Ions& ions)
     {
         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 dc54a8e7..9a9bf8a6 100644
--- a/src/MGmolInterface.h
+++ b/src/MGmolInterface.h
@@ -40,9 +40,6 @@ 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/src/MVPSolver.cc b/src/MVPSolver.cc
index 7ab9f665..6990e29d 100644
--- a/src/MVPSolver.cc
+++ b/src/MVPSolver.cc
@@ -256,7 +256,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
             const double ts0
                 = evalEntropyMVP(current_proj_mat, (ct.verbose > 1), os_);
             const double e0 = energy_->evaluateTotal(
-                ts0, current_proj_mat, orbitals, printE, os_);
+                ts0, current_proj_mat, ions_, orbitals, printE, os_);
 
             MatrixType target("target", numst_, numst_);
 
@@ -309,7 +309,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                 const double ts1
                     = evalEntropyMVP(proj_mat_work_, (ct.verbose > 2), os_);
                 const double e1 = energy_->evaluateTotal(
-                    ts1, proj_mat_work_, orbitals, ct.verbose - 1, os_);
+                    ts1, proj_mat_work_, ions_, orbitals, ct.verbose - 1, os_);
 
                 // line minimization
                 const double beta
diff --git a/src/PolakRibiereSolver.cc b/src/PolakRibiereSolver.cc
index 6b5f88de..3f46a00e 100644
--- a/src/PolakRibiereSolver.cc
+++ b/src/PolakRibiereSolver.cc
@@ -243,9 +243,9 @@ double PolakRibiereSolver<OrbitalsType>::evaluateEnergy(
     const OrbitalsType& orbitals, const bool print_flag)
 {
     // Get the new total energy
-    const double ts = 0.5 * proj_matrices_->computeEntropy(); // in [Ha]
-    const double energy
-        = energy_->evaluateTotal(ts, proj_matrices_, orbitals, print_flag, os_);
+    const double ts     = 0.5 * proj_matrices_->computeEntropy(); // in [Ha]
+    const double energy = energy_->evaluateTotal(
+        ts, proj_matrices_, ions_, orbitals, print_flag, os_);
 
     // Control& ct(*(Control::instance()));
     // if( ct.verbose>2 && onpe0 )os_<<"energy="<<energy<<endl;
@@ -369,7 +369,7 @@ int PolakRibiereSolver<OrbitalsType>::solve(OrbitalsType& orbitals,
 
             // evaluate residuals, preconditioned residuals for current orbitals
             double normRes = mgmol_strategy_->computeResidual(orbitals,
-                work_orbitals, *r_k_, (print_res || ct.checkResidual()),
+                work_orbitals, ions_, *r_k_, (print_res || ct.checkResidual()),
                 ct.checkResidual());
             if (normRes < ct.conv_tol && ct.checkResidual())
             {
diff --git a/src/lbfgsrlx.cc b/src/lbfgsrlx.cc
index 6698eef0..ce05c1a3 100644
--- a/src/lbfgsrlx.cc
+++ b/src/lbfgsrlx.cc
@@ -14,7 +14,6 @@
 #include "Energy.h"
 #include "Ions.h"
 #include "LBFGS.h"
-#include "LBFGS_IonicStepper.h"
 #include "LocalizationRegions.h"
 #include "MGmol.h"
 #include "MGmol_blas1.h"
diff --git a/src/quench.cc b/src/quench.cc
index 073798f6..6eed74bb 100644
--- a/src/quench.cc
+++ b/src/quench.cc
@@ -600,8 +600,8 @@ int MGmol<OrbitalsType>::quench(OrbitalsType& orbitals, Ions& ions,
                 << " TS             [Ha] = " << ts << std::endl;
         }
     }
-    last_eks
-        = energy_->evaluateTotal(ts, proj_matrices_.get(), orbitals, 2, os_);
+    last_eks = energy_->evaluateTotal(
+        ts, proj_matrices_.get(), ions, orbitals, 2, os_);
 
     if (ct.computeCondGramMD())
     {
diff --git a/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc b/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc
index 3027ab55..b1fb0094 100644
--- a/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc
+++ b/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc
@@ -124,8 +124,6 @@ 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();
@@ -195,8 +193,6 @@ int main(int argc, char** argv)
             }
         }
 
-        mgmol->setPositions(positions, anumbers);
-
         //
         // evaluate energy and forces with wavefunctions just read
         //