diff --git a/CMakeLists.txt b/CMakeLists.txt
index 75752a3a..6c48a2e9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -178,7 +178,7 @@ if(${MGMOL_WITH_CLANG_FORMAT})
   find_package(CLANG_FORMAT)
   if(${CLANG_FORMAT_FOUND})
     message(STATUS "Indent with clang-format")
-    file(GLOB_RECURSE FORMAT_SOURCES src/*.cc src/*.h tests/*.cc tests/*.h)
+    file(GLOB_RECURSE FORMAT_SOURCES src/*.cc src/*.h tests/*.cc tests/*.h drivers/*.cc)
     add_custom_target(format
       COMMAND ${CLANG_FORMAT_EXECUTABLE} -i -style=file ${FORMAT_SOURCES}
       DEPENDS ${FORMAT_SOURCES})
@@ -271,5 +271,7 @@ link_libraries(${LIBROM_LIB})
 # add subdirectories for source files, tests
 add_subdirectory(src)
 
+add_subdirectory(drivers)
+
 add_subdirectory(tests)
 
diff --git a/drivers/CMakeLists.txt b/drivers/CMakeLists.txt
new file mode 100644
index 00000000..4eb1c741
--- /dev/null
+++ b/drivers/CMakeLists.txt
@@ -0,0 +1,11 @@
+include_directories( ${CMAKE_SOURCE_DIR}/src )
+
+add_executable(check_input check_input.cc)
+
+add_executable(example1 example1.cc)
+
+target_include_directories(check_input PRIVATE ${Boost_INCLUDE_DIRS})
+target_include_directories(example1 PRIVATE ${Boost_INCLUDE_DIRS})
+
+target_link_libraries(check_input mgmol_src)
+target_link_libraries(example1 mgmol_src)
diff --git a/drivers/check_input.cc b/drivers/check_input.cc
new file mode 100644
index 00000000..14dbaa09
--- /dev/null
+++ b/drivers/check_input.cc
@@ -0,0 +1,93 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+
+#include "Control.h"
+#include "ExtendedGridOrbitals.h"
+#include "LocGridOrbitals.h"
+#include "MGmol.h"
+#include "MGmol_MPI.h"
+#include "MPIdata.h"
+#include "mgmol_run.h"
+
+#include <cassert>
+#include <iostream>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+int main(int argc, char** argv)
+{
+    int mpirc = MPI_Init(&argc, &argv);
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Initialization failed!!!" << std::endl;
+        MPI_Abort(MPI_COMM_WORLD, 0);
+    }
+
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    mgmol_init(comm);
+
+    // read runtime parameters
+    std::string input_filename("");
+    std::string lrs_filename;
+    std::string constraints_filename("");
+
+    float total_spin = 0.;
+    bool with_spin   = false;
+
+    po::variables_map vm;
+
+    // read options from PE0 only
+    if (MPIdata::onpe0)
+    {
+        read_config(argc, argv, vm, input_filename, lrs_filename,
+            constraints_filename, total_spin, with_spin);
+    }
+
+    MGmol_MPI::setup(comm, std::cout, with_spin);
+    MGmol_MPI& mmpi      = *(MGmol_MPI::instance());
+    MPI_Comm global_comm = mmpi.commGlobal();
+
+    Control::setup(global_comm, with_spin, total_spin);
+    Control& ct = *(Control::instance());
+
+    ct.setOptions(vm);
+
+    int ret = ct.checkOptions();
+    if (ret < 0) return ret;
+
+    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(lrs_filename);
+
+    // Enter main scope
+    {
+        MGmolInterface* mgmol;
+        if (ct.isLocMode())
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+        else
+            mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+
+        *MPIdata::sout << " Input parameters OK\n";
+
+        delete mgmol;
+    } // close main scope
+
+    mgmol_finalize();
+
+    mpirc = MPI_Finalize();
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
+    }
+
+    return 0;
+}
diff --git a/drivers/example1.cc b/drivers/example1.cc
new file mode 100644
index 00000000..561f5e08
--- /dev/null
+++ b/drivers/example1.cc
@@ -0,0 +1,167 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+
+#include "Control.h"
+#include "ExtendedGridOrbitals.h"
+#include "LocGridOrbitals.h"
+#include "MGmol.h"
+#include "MGmol_MPI.h"
+#include "MPIdata.h"
+#include "mgmol_run.h"
+
+#include <cassert>
+#include <iostream>
+#include <time.h>
+#include <vector>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+int main(int argc, char** argv)
+{
+    int mpirc = MPI_Init(&argc, &argv);
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Initialization failed!!!" << std::endl;
+        MPI_Abort(MPI_COMM_WORLD, 0);
+    }
+
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    /*
+     * Initialize general things, like magma, openmp, IO, ...
+     */
+    mgmol_init(comm);
+
+    /*
+     * read runtime parameters
+     */
+    std::string input_filename("");
+    std::string lrs_filename;
+    std::string constraints_filename("");
+
+    float total_spin = 0.;
+    bool with_spin   = false;
+
+    po::variables_map vm;
+
+    // read from PE0 only
+    if (MPIdata::onpe0)
+    {
+        read_config(argc, argv, vm, input_filename, lrs_filename,
+            constraints_filename, total_spin, with_spin);
+    }
+
+    MGmol_MPI::setup(comm, std::cout, with_spin);
+    MGmol_MPI& mmpi      = *(MGmol_MPI::instance());
+    MPI_Comm global_comm = mmpi.commGlobal();
+
+    /*
+     * Setup control struct with run time parameters
+     */
+    Control::setup(global_comm, with_spin, total_spin);
+    Control& ct = *(Control::instance());
+
+    ct.setOptions(vm);
+
+    int ret = ct.checkOptions();
+    if (ret < 0) return ret;
+
+    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(lrs_filename);
+
+    // Enter main scope
+    {
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "Construct MGmol object..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+
+        MGmolInterface* mgmol;
+        if (ct.isLocMode())
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+        else
+            mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "MGmol setup..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+        mgmol->setup();
+
+        if (MPIdata::onpe0)
+        {
+            std::cout << "-------------------------" << std::endl;
+            std::cout << "Setup done..." << std::endl;
+            std::cout << "-------------------------" << std::endl;
+        }
+
+        // here we just use the atomic positions read in and used
+        // to initialize MGmol
+        std::vector<double> positions;
+        mgmol->getAtomicPositions(positions);
+        std::vector<short> anumbers;
+        mgmol->getAtomicNumbers(anumbers);
+        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++;
+            }
+        }
+
+        // compute energy and forces using all MPI tasks
+        // expect positions to be replicated on all MPI tasks
+        std::vector<double> forces;
+        mgmol->evaluateEnergyAndForces(positions, anumbers, forces);
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+            }
+        }
+
+        delete mgmol;
+
+    } // close main scope
+
+    mgmol_finalize();
+
+    mpirc = MPI_Finalize();
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
+    }
+
+    time_t tt;
+    time(&tt);
+    if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl;
+
+    return 0;
+}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index e33d03d2..e602101f 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -191,13 +191,15 @@ install(TARGETS mgmol_src DESTINATION lib)
 add_executable(mgmol-opt main.cc)
 target_include_directories (mgmol-opt PRIVATE ${Boost_INCLUDE_DIRS})
 
-target_link_libraries(mgmol-opt mgmol_src ${link_libs})
-target_link_libraries(mgmol-opt ${SCALAPACK_LIBRARIES})
-target_link_libraries(mgmol-opt ${HDF5_LIBRARIES})
-target_link_libraries(mgmol-opt ${HDF5_HL_LIBRARIES})
-target_link_libraries(mgmol-opt ${BLAS_LIBRARIES})
-target_link_libraries(mgmol-opt ${LAPACK_LIBRARIES})
-target_link_libraries(mgmol-opt ${Boost_LIBRARIES})
+target_link_libraries(mgmol_src ${link_libs})
+target_link_libraries(mgmol_src ${SCALAPACK_LIBRARIES})
+target_link_libraries(mgmol_src ${LAPACK_LIBRARIES})
+target_link_libraries(mgmol_src ${BLAS_LIBRARIES})
+target_link_libraries(mgmol_src ${HDF5_LIBRARIES})
+target_link_libraries(mgmol_src ${HDF5_HL_LIBRARIES})
+target_link_libraries(mgmol_src ${Boost_LIBRARIES})
+
+target_link_libraries(mgmol-opt mgmol_src)
 if (${OPENMP_CXX_FOUND})
   target_link_libraries(mgmol-opt OpenMP::OpenMP_CXX)
 endif()
diff --git a/src/DFTsolver.cc b/src/DFTsolver.cc
index c1484160..91c4a0f6 100644
--- a/src/DFTsolver.cc
+++ b/src/DFTsolver.cc
@@ -339,13 +339,15 @@ int DFTsolver<OrbitalsType>::solve(OrbitalsType& orbitals,
             else
             {
                 bool updateDM = false;
+                bool updatedS = false;
                 if (!ct.fullyOccupied())
                 {
                     orbitals.computeGramAndInvS();
                     dm_strategy_->dressDM();
                     updateDM = true;
+                    updatedS = true;
                 }
-                orbitals.orthonormalizeLoewdin(true, nullptr, updateDM);
+                orbitals.orthonormalizeLoewdin(updatedS, nullptr, updateDM);
 
                 orbitals_stepper_->restartMixing();
             }
diff --git a/src/Electrostatic.cc b/src/Electrostatic.cc
index b0afeb76..80cb754b 100644
--- a/src/Electrostatic.cc
+++ b/src/Electrostatic.cc
@@ -36,7 +36,7 @@ Timer Electrostatic::solve_tm_("Electrostatic::solve");
 
 Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
     const double screening_const)
-    : laptype_(lap_type)
+    : laptype_(lap_type), poisson_solver_(nullptr)
 {
     assert(bcPoisson[0] >= 0);
     assert(bcPoisson[1] >= 0);
@@ -166,6 +166,8 @@ Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
 
 Electrostatic::~Electrostatic()
 {
+    assert(poisson_solver_ != nullptr);
+
     delete poisson_solver_;
     if (grhod_ != nullptr) delete grhod_;
     if (grhoc_ != nullptr) delete grhoc_;
diff --git a/src/Ions.cc b/src/Ions.cc
index 53ee0cf4..940e12e8 100644
--- a/src/Ions.cc
+++ b/src/Ions.cc
@@ -10,6 +10,7 @@
 #include "Ions.h"
 #include "Control.h"
 #include "HDFrestart.h"
+#include "MGmol_MPI.h"
 #include "MGmol_blas1.h"
 #include "MPIdata.h"
 #include "Mesh.h"
@@ -169,12 +170,9 @@ void Ions::computeMaxNumProjs()
                          << " initialized" << std::endl;
 #endif
 
-    std::vector<Ion*>::iterator ion = list_ions_.begin();
-    while (ion != list_ions_.end())
+    for (auto& ion : list_ions_)
     {
-        max_num_proj_ = std::max(max_num_proj_, (*ion)->nProjectors());
-
-        ion++;
+        max_num_proj_ = std::max(max_num_proj_, ion->nProjectors());
     }
 
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
@@ -210,11 +208,9 @@ void Ions::setup()
     if (ct.verbose > 0)
         printWithTimeStamp("Ions::setup()... individual ions...", std::cout);
 
-    std::vector<Ion*>::iterator ion = list_ions_.begin();
-    while (ion != list_ions_.end())
+    for (auto& ion : list_ions_)
     {
-        (*ion)->setup();
-        ion++;
+        ion->setup();
     }
 
     setMapVL();
@@ -243,7 +239,7 @@ Ions::~Ions()
 void Ions::setupListOverlappingIons()
 {
     Control& ct = *(Control::instance());
-    if (ct.verbose > 0)
+    if (ct.verbose > 1)
         printWithTimeStamp("Ions::setupListOverlappingIons()...", std::cout);
 
     overlappingNL_ions_.clear();
@@ -269,7 +265,7 @@ void Ions::setupListOverlappingIons()
 void Ions::setupInteractingIons()
 {
     Control& ct = *(Control::instance());
-    if (ct.verbose > 0)
+    if (ct.verbose > 1)
         printWithTimeStamp("Ions::setupInteractingIons()...", std::cout);
 
     ions_setupInteractingIons_tm.start();
@@ -1569,7 +1565,7 @@ Ion* Ions::findLocalIon(const int index) const
     return nullptr;
 }
 
-void Ions::setPositions(const std::vector<double>& tau)
+void Ions::setLocalPositions(const std::vector<double>& tau)
 {
     assert(tau.size() == 3 * local_ions_.size());
 
@@ -1586,69 +1582,6 @@ void Ions::setPositions(const std::vector<double>& tau)
     setup_ = false;
 }
 
-void Ions::get_positions(std::vector<std::vector<double>>& rr) const
-{
-    assert(rr.size() == local_ions_.size());
-    if (local_ions_.empty()) return;
-    std::vector<double> tau(3);
-    int i                                 = 0;
-    std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-    while (ion != local_ions_.end())
-    {
-        tau[0] = (*ion)->position(0);
-        tau[1] = (*ion)->position(1);
-        tau[2] = (*ion)->position(2);
-        rr[i]  = tau;
-        ion++;
-        i++;
-    }
-}
-void Ions::set_positions(const std::vector<std::vector<double>>& rr)
-{
-    assert(rr.size() == local_ions_.size());
-
-    if (local_ions_.empty()) return;
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    int i                           = 0;
-    while (ion != local_ions_.end())
-    {
-        assert(rr[i].size() == 3);
-        (*ion)->setPosition(rr[i][0], rr[i][1], rr[i][2]);
-        ion++;
-        i++;
-    }
-}
-void Ions::get_forces(std::vector<std::vector<double>>& ff) const
-{
-    assert(ff.size() == local_ions_.size());
-
-    if (local_ions_.empty()) return;
-    int i                                 = 0;
-    std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-    while (ion != local_ions_.end())
-    {
-        ff[i][0] = (*ion)->force(0);
-        ff[i][1] = (*ion)->force(1);
-        ff[i][2] = (*ion)->force(2);
-        ion++;
-        i++;
-    }
-}
-void Ions::set_forces(const std::vector<std::vector<double>>& ff)
-{
-    assert(ff.size() == local_ions_.size());
-
-    if (local_ions_.empty()) return;
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    int i                           = 0;
-    while (ion != local_ions_.end())
-    {
-        (*ion)->setForce(ff[i][0], ff[i][1], ff[i][2]);
-        ion++;
-        i++;
-    }
-}
-
 int Ions::readAtoms(const std::string& filename, const bool cell_relative)
 {
 
@@ -1785,6 +1718,7 @@ int Ions::readAtomsFromXYZ(
             }
             ++count;
         }
+        delete tfile;
     }
 
     mmpi.bcastGlobal(&count, 1);
@@ -1793,6 +1727,17 @@ int Ions::readAtomsFromXYZ(
     mmpi.bcastGlobal(&crds[0], 3 * natoms);
     mmpi.bcastGlobal(&spec[0], natoms);
 
+    return setAtoms(crds, spec);
+}
+
+int Ions::setAtoms(
+    const std::vector<double>& crds, const std::vector<short>& spec)
+{
+    MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
+    Control& ct(*(Control::instance()));
+
+    const int natoms = crds.size() / 3;
+
     double velocity[3] = { 0., 0., 0. };
     bool locked        = false;
     for (int ia = 0; ia < natoms; ++ia)
@@ -1820,7 +1765,7 @@ int Ions::readAtomsFromXYZ(
         }
         if (spname.compare("") == 0)
         {
-            (*MPIdata::serr) << "Ions::readAtomsFromXYZ() --- ERROR: unknown "
+            (*MPIdata::serr) << "Ions::setAtoms() --- ERROR: unknown "
                                 "species for atomic number "
                              << spec[ia] << std::endl;
             return -1;
@@ -1841,7 +1786,7 @@ int Ions::readAtomsFromXYZ(
 
         // Populate list_ions_ list
         // std::cout<<"crds: "<<crds[3*ia+0]<<", "<<crds[3*ia+1]<<",
-        // "<<crds[3*ia+2]<<endl;
+        // "<<crds[3*ia+2]<<std::endl;
         if (inListIons(crds[3 * ia + 0], crds[3 * ia + 1], crds[3 * ia + 2]))
         {
             list_ions_.push_back(new_ion);
@@ -1878,8 +1823,6 @@ int Ions::readAtomsFromXYZ(
     //    std::cout<<mmpi.mype()<<"...list size = "<<list_ions_.size()<<" local
     //    ions size = "<<local_ions_.size()<<endl;
 
-    if (mmpi.PE0()) delete tfile;
-
     return natoms;
 }
 
@@ -2118,8 +2061,6 @@ int Ions::read1atom(std::ifstream* tfile, const bool cell_relative)
     // Ion needs to be created by each MPI task to set global ids
     Ion* new_ion = new Ion(species_[isp], name_read, crds, velocity, locked);
 
-    // assert( (new_ion->locked() == false) || (new_ion->locked() == true) );
-
 #ifdef DEBUG
     if (onpe0) (*MPIdata::sout) << "Ion read..." << std::endl;
 #endif
@@ -2131,11 +2072,11 @@ int Ions::read1atom(std::ifstream* tfile, const bool cell_relative)
         // populate local_ions_ list
         if (inLocalIons(crds[0], crds[1], crds[2]))
         {
-            (new_ion)->set_here(true);
+            new_ion->set_here(true);
             local_ions_.push_back(new_ion);
         }
         else
-            (new_ion)->set_here(false);
+            new_ion->set_here(false);
     }
     else
     {
@@ -2193,7 +2134,7 @@ void Ions::setVelocities(const std::vector<double>& tau0,
     }
 }
 
-void Ions::getPositions(std::vector<double>& tau) const
+void Ions::getLocalPositions(std::vector<double>& tau) const
 {
     assert(tau.size() == 3 * local_ions_.size());
 
@@ -2207,6 +2148,42 @@ void Ions::getPositions(std::vector<double>& tau) const
     }
 }
 
+void Ions::getPositions(std::vector<double>& tau)
+{
+    std::vector<double> tau_local(3 * local_ions_.size());
+
+    getLocalPositions(tau_local);
+
+    MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+    mmpi.allGatherV(tau_local, tau);
+}
+
+void Ions::getAtomicNumbers(std::vector<short>& atnumbers)
+{
+    std::vector<short> local_atnumbers;
+
+    for (auto& ion : local_ions_)
+    {
+        local_atnumbers.push_back(ion->atomic_number());
+    }
+
+    MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+    mmpi.allGatherV(local_atnumbers, atnumbers);
+}
+
+void Ions::getForces(std::vector<double>& forces)
+{
+    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);
+}
+
 void Ions::setTau0()
 {
     assert(tau0_.size() == 3 * local_ions_.size());
@@ -2236,6 +2213,20 @@ void Ions::setPositionsToTau0()
     }
 }
 
+void Ions::setPositions(
+    const std::vector<double>& tau, const std::vector<short>& anumbers)
+{
+    assert(tau.size() == anumbers.size() * 3);
+
+    // clear previous data
+    clearLists();
+
+    num_ions_ = setAtoms(tau, anumbers);
+
+    // setup required after updating local ions positions
+    setup();
+}
+
 void Ions::setVelocitiesToVel()
 {
     assert(velocity_.size() == 3 * local_ions_.size());
@@ -2249,7 +2240,7 @@ void Ions::setVelocitiesToVel()
     }
 }
 
-void Ions::getForces(std::vector<double>& tau) const
+void Ions::getLocalForces(std::vector<double>& tau) const
 {
     assert(tau.size() == 3 * local_ions_.size());
 
@@ -2380,12 +2371,10 @@ double Ions::computeMaxNLprojRadius() const
 {
     double radius = 0.;
 
-    std::vector<Ion*>::const_iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    for (auto& iion : local_ions_)
     {
-        double r = (*iion)->radiusNLproj();
+        double r = iion->radiusNLproj();
         radius   = r > radius ? r : radius;
-        iion++;
     }
 
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
@@ -3061,7 +3050,7 @@ void Ions::updateForcesInteractingIons()
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
 
     // get computed forces into fion_
-    getForces(fion_);
+    getLocalForces(fion_);
 
     // initialize with local names and forces
     DistributedIonicData forces_data(local_names_, fion_);
@@ -3136,6 +3125,18 @@ void Ions::updateTaupInteractingIons()
     }
 }
 
+void Ions::clearLists()
+{
+    local_ions_.clear();
+    std::vector<Ion*>::iterator ion = list_ions_.begin();
+    while (ion != list_ions_.end())
+    {
+        delete *ion;
+        ion++;
+    }
+    list_ions_.clear();
+}
+
 // update list of local ions
 void Ions::updateListIons()
 {
@@ -3160,15 +3161,7 @@ void Ions::updateListIons()
     // Note: this is based on data from MD std::vectors
 
     // First cleanup list_ions_
-    local_ions_.clear();
-    // delete current ions from list
-    std::vector<Ion*>::iterator ion = list_ions_.begin();
-    while (ion != list_ions_.end())
-    {
-        delete *ion;
-        ion++;
-    }
-    list_ions_.clear();
+    clearLists();
 
     // Update list starting with local ion data.
     // This enables overlapping data accumulation with communication.
diff --git a/src/Ions.h b/src/Ions.h
index 2351d678..772c51d1 100644
--- a/src/Ions.h
+++ b/src/Ions.h
@@ -16,7 +16,6 @@
 #include <map>
 #include <vector>
 
-#include "DataDistribution.h"
 #include "DistributedIonicData.h"
 #include "Ion.h"
 #include "hdf5.h"
@@ -37,7 +36,12 @@ class Ions
     const std::vector<Species>& species_;
 
     std::vector<Ion*> list_ions_;
-    std::vector<Ion*> local_ions_; // centered in local sub-domain
+
+    /*
+     * ions located in local sub-domain
+     */
+    std::vector<Ion*> local_ions_;
+
     std::vector<Ion*> interacting_ions_; // for ion-ion interactions
     std::vector<Ion*>
         overlappingNL_ions_; // with projectors overlapping local sub-domain
@@ -62,7 +66,6 @@ class Ions
     void readRestartPositions(HDFrestart& h5_file);
     int read1atom(std::ifstream* tfile, const bool cell_relative);
 
-    // void associate2PE();
     void setupInteractingIons();
     void setupListOverlappingIons();
     void setMapVL();
@@ -158,6 +161,7 @@ class Ions
     void gatherForces(
         std::vector<double>& forces, const int root, const MPI_Comm comm) const;
     bool hasLockedAtoms() const;
+    void clearLists();
 
 public:
     Ions(const double lat[3], const std::vector<Species>& sp);
@@ -247,11 +251,7 @@ class Ions
     // check if ion is in list of local ions
     bool isLocal(const std::string& ion_name) const;
 
-    void setPositions(const std::vector<double>& tau);
-    void get_positions(std::vector<std::vector<double>>& r) const;
-    void set_positions(const std::vector<std::vector<double>>& r);
-    void get_forces(std::vector<std::vector<double>>& f) const;
-    void set_forces(const std::vector<std::vector<double>>& f);
+    void setLocalPositions(const std::vector<double>& tau);
 
     void lockAtom(const std::string& name);
 
@@ -260,6 +260,8 @@ class Ions
     int readAtoms(const std::string& filename, const bool cell_relative);
     int readAtoms(std::ifstream* tfile, const bool cell_relative);
     void initFromRestartFile(HDFrestart& h5_file);
+    int setAtoms(
+        const std::vector<double>& crds, const std::vector<short>& spec);
 
     int getNValenceElectrons() const;
     void syncForces();
@@ -268,9 +270,15 @@ class Ions
     void setTau0();
     void setPositionsToTau0();
     void setVelocitiesToVel();
+    void setPositions(
+        const std::vector<double>& tau, const std::vector<short>& anumbers);
+
+    void getLocalPositions(std::vector<double>& tau) const;
+    void getPositions(std::vector<double>& tau);
+    void getAtomicNumbers(std::vector<short>& atnumbers);
 
-    void getPositions(std::vector<double>& tau) const;
-    void getForces(std::vector<double>& tau) const;
+    void getForces(std::vector<double>& forces);
+    void getLocalForces(std::vector<double>& tau) const;
     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 39f30378..c91120e9 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -113,73 +113,29 @@ template <class OrbitalsType>
 MGmol<OrbitalsType>::MGmol(MPI_Comm comm, std::ostream& os,
     std::string input_filename, std::string lrs_filename,
     std::string constraints_filename)
-    : os_(os)
+    : os_(os), comm_(comm)
 {
-    comm_ = comm;
-
-    constraints_ = new ConstraintSet();
+    constraints_.reset(new ConstraintSet());
 
     mgmol::out = &os;
 
-    geom_optimizer_ = nullptr;
-    local_cluster_  = nullptr;
-    proj_matrices_  = nullptr;
-    dm_strategy_    = nullptr;
-
-    h5f_file_ = nullptr;
-
-    aomm_ = nullptr;
-
-    spreadf_ = nullptr;
-
-    spread_penalty_ = nullptr;
-
-    orbitals_precond_ = nullptr;
-
-    forces_ = nullptr;
-
-    energy_ = nullptr;
+    current_orbitals_ = nullptr;
 
     setupFromInput(input_filename);
 
+    /*
+     * Extra setup if using localization regions
+     */
     setupLRs(lrs_filename);
 
-    setupConstraintsFromInput(constraints_filename);
-
-    setup();
+    if (!constraints_filename.empty())
+        setupConstraintsFromInput(constraints_filename);
 }
 
 template <class OrbitalsType>
 MGmol<OrbitalsType>::~MGmol()
 {
-    delete electrostat_;
-    delete rho_;
-    delete constraints_;
-    delete xcongrid_;
-    delete energy_;
-    if (hamiltonian_ != nullptr) delete hamiltonian_;
-    if (geom_optimizer_ != nullptr) delete geom_optimizer_;
-
-    if (currentMasks_ != nullptr) delete currentMasks_;
-    if (corrMasks_ != nullptr) delete corrMasks_;
-
-    if (aomm_ != nullptr) delete aomm_;
-
-    delete current_orbitals_;
-    delete ions_;
-    delete g_kbpsi_;
-
-    delete proj_matrices_;
-    if (local_cluster_ != nullptr) delete local_cluster_;
-
-    if (h5f_file_ != nullptr) delete h5f_file_;
-
-    if (spreadf_ != nullptr) delete spreadf_;
-
-    if (spread_penalty_ != nullptr) delete spread_penalty_;
-
-    delete forces_;
-    if (dm_strategy_ != nullptr) delete dm_strategy_;
+    if (current_orbitals_) delete current_orbitals_;
 }
 
 template <>
@@ -192,18 +148,17 @@ void MGmol<LocGridOrbitals>::initialMasks()
     if (ct.verbose > 0)
         printWithTimeStamp("MGmol<OrbitalsType>::initialMasks()...", os_);
 
-    currentMasks_ = new MasksSet(false, ct.getMGlevels());
+    currentMasks_
+        = std::shared_ptr<MasksSet>(new MasksSet(false, ct.getMGlevels()));
     currentMasks_->setup(lrs_);
 
-    corrMasks_ = new MasksSet(true, 0);
+    corrMasks_ = std::shared_ptr<MasksSet>(new MasksSet(true, 0));
     corrMasks_->setup(lrs_);
 }
 
 template <>
 void MGmol<ExtendedGridOrbitals>::initialMasks()
 {
-    currentMasks_ = nullptr;
-    corrMasks_    = nullptr;
 }
 
 template <class OrbitalsType>
@@ -227,7 +182,7 @@ int MGmol<OrbitalsType>::initial()
 
     pb::Lap<ORBDTYPE>* lapop
         = ct.Mehrstellen() ? hamiltonian_->lapOper() : nullptr;
-    g_kbpsi_ = new KBPsiMatrixSparse(lapop);
+    g_kbpsi_.reset(new KBPsiMatrixSparse(lapop));
 
     check_anisotropy();
 
@@ -269,7 +224,7 @@ int MGmol<OrbitalsType>::initial()
         // set number of iterations to 10.
         if (ct.load_balancing_alpha > 0.0)
         {
-            local_cluster_ = new ClusterOrbitals(lrs_);
+            local_cluster_.reset(new ClusterOrbitals(lrs_));
             local_cluster_->setup();
             local_cluster_->computeClusters(ct.load_balancing_max_iterations);
         }
@@ -289,29 +244,31 @@ int MGmol<OrbitalsType>::initial()
     {
 #ifdef HAVE_MAGMA
         if (use_replicated_matrix)
-            proj_matrices_ = new ProjectedMatricesMehrstellen<ReplicatedMatrix>(
-                ct.numst, with_spin, ct.occ_width);
+            proj_matrices_.reset(
+                new ProjectedMatricesMehrstellen<ReplicatedMatrix>(
+                    ct.numst, with_spin, ct.occ_width));
         else
 #endif
-            proj_matrices_ = new ProjectedMatricesMehrstellen<
+            proj_matrices_.reset(new ProjectedMatricesMehrstellen<
                 dist_matrix::DistMatrix<DISTMATDTYPE>>(
-                ct.numst, with_spin, ct.occ_width);
+                ct.numst, with_spin, ct.occ_width));
     }
     else if (ct.short_sighted)
-        proj_matrices_ = new ProjectedMatricesSparse(
-            ct.numst, ct.occ_width, lrs_, local_cluster_);
+        proj_matrices_.reset(new ProjectedMatricesSparse(
+            ct.numst, ct.occ_width, lrs_, local_cluster_.get()));
     else
 #ifdef HAVE_MAGMA
         if (use_replicated_matrix)
-        proj_matrices_ = new ProjectedMatrices<ReplicatedMatrix>(
-            ct.numst, with_spin, ct.occ_width);
+        proj_matrices_.reset(new ProjectedMatrices<ReplicatedMatrix>(
+            ct.numst, with_spin, ct.occ_width));
     else
 #endif
-        proj_matrices_
-            = new ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>(
-                ct.numst, with_spin, ct.occ_width);
+        proj_matrices_.reset(
+            new ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>(
+                ct.numst, with_spin, ct.occ_width));
 
-    forces_ = new Forces<OrbitalsType>(hamiltonian_, rho_, proj_matrices_);
+    forces_.reset(new Forces<OrbitalsType>(
+        hamiltonian_.get(), rho_.get(), proj_matrices_.get()));
 
     if (ct.verbose > 0)
         printWithTimeStamp(
@@ -326,8 +283,8 @@ int MGmol<OrbitalsType>::initial()
         printWithTimeStamp("MGmol<OrbitalsType>::initial(), create T...", os_);
 
     current_orbitals_ = new OrbitalsType("Primary", mygrid, mymesh->subdivx(),
-        ct.numst, ct.bcWF, proj_matrices_, lrs_, currentMasks_, corrMasks_,
-        local_cluster_, true);
+        ct.numst, ct.bcWF, proj_matrices_.get(), lrs_, currentMasks_.get(),
+        corrMasks_.get(), local_cluster_.get(), true);
 
     increaseMemorySlotsForOrbitals<MemorySpaceType>();
 
@@ -372,7 +329,7 @@ int MGmol<OrbitalsType>::initial()
                 "MGmol<OrbitalsType>::initial(), init wf and masks...", os_);
 
         // Make temp mask for initial random wave functions
-        if (ct.init_loc == 1 && currentMasks_ != nullptr)
+        if (ct.init_loc == 1 && currentMasks_)
         {
             float cut_init = ct.initRadius();
             assert(cut_init > 0.);
@@ -383,7 +340,7 @@ int MGmol<OrbitalsType>::initial()
         current_orbitals_->initWF(lrs_);
 
         // initialize masks again
-        if (ct.init_loc == 1 && currentMasks_ != nullptr)
+        if (ct.init_loc == 1 && currentMasks_)
         {
             currentMasks_->initialize(lrs_, 0);
         }
@@ -427,9 +384,10 @@ int MGmol<OrbitalsType>::initial()
     }
 
     if (ct.verbose > 0) printWithTimeStamp("Initialize XC functional...", os_);
-    xcongrid_ = XCfunctionalFactory<OrbitalsType>::create(
-        ct.xctype, mmpi.nspin(), *rho_, pot);
-    assert(xcongrid_ != nullptr);
+    xcongrid_
+        = std::shared_ptr<XConGrid>(XCfunctionalFactory<OrbitalsType>::create(
+            ct.xctype, mmpi.nspin(), *rho_, pot));
+    assert(xcongrid_);
 
     // initialize nl potentials with restart values if possible
     //    if( ct.restart_info>1 )
@@ -473,7 +431,7 @@ int MGmol<OrbitalsType>::initial()
     {
         Vector3D origin(mygrid.origin(0), mygrid.origin(1), mygrid.origin(2));
         Vector3D ll(mygrid.ll(0), mygrid.ll(1), mygrid.ll(2));
-        spreadf_ = new SpreadsAndCenters<OrbitalsType>(origin, ll);
+        spreadf_.reset(new SpreadsAndCenters<OrbitalsType>(origin, ll));
     }
 
     bool energy_with_spread_penalty = false;
@@ -481,26 +439,30 @@ int MGmol<OrbitalsType>::initial()
     {
         if (ct.isSpreadFunctionalVolume())
         {
-            spread_penalty_ = new SpreadPenaltyVolume<OrbitalsType>(spreadf_,
-                ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(),
-                ct.spreadPenaltyDampingFactor());
+            spread_penalty_.reset(
+                new SpreadPenaltyVolume<OrbitalsType>(spreadf_.get(),
+                    ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(),
+                    ct.spreadPenaltyDampingFactor()));
         }
         else if (ct.isSpreadFunctionalEnergy())
         {
             energy_with_spread_penalty = true;
-            spread_penalty_ = new EnergySpreadPenalty<OrbitalsType>(spreadf_,
-                ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor());
+            spread_penalty_.reset(
+                new EnergySpreadPenalty<OrbitalsType>(spreadf_.get(),
+                    ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor()));
         }
         else
-            spread_penalty_ = new SpreadPenalty<OrbitalsType>(spreadf_,
-                ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(),
-                ct.spreadPenaltyDampingFactor());
+            spread_penalty_.reset(
+                new SpreadPenalty<OrbitalsType>(spreadf_.get(),
+                    ct.spreadPenaltyTarget(), ct.spreadPenaltyAlphaFactor(),
+                    ct.spreadPenaltyDampingFactor()));
     }
 
-    SpreadPenaltyInterface<OrbitalsType>* spread_penalty
+    std::shared_ptr<SpreadPenaltyInterface<OrbitalsType>> spread_penalty
         = energy_with_spread_penalty ? spread_penalty_ : nullptr;
-    energy_ = new Energy<OrbitalsType>(
-        mygrid, *ions_, pot, *electrostat_, *rho_, *xcongrid_, spread_penalty);
+    energy_ = std::shared_ptr<Energy<OrbitalsType>>(
+        new Energy<OrbitalsType>(mygrid, *ions_, pot, *electrostat_, *rho_,
+            *xcongrid_, spread_penalty.get()));
 
     if (ct.verbose > 0) printWithTimeStamp("Setup matrices...", os_);
 
@@ -509,15 +471,16 @@ int MGmol<OrbitalsType>::initial()
     // HMVP algorithm requires that H is initialized
 #ifdef HAVE_MAGMA
     if (use_replicated_matrix)
-        dm_strategy_
-            = DMStrategyFactory<OrbitalsType, ReplicatedMatrix>::create(comm_,
-                os_, *ions_, rho_, energy_, electrostat_, this, proj_matrices_,
-                current_orbitals_);
+        dm_strategy_.reset(
+            DMStrategyFactory<OrbitalsType, ReplicatedMatrix>::create(comm_,
+                os_, *ions_, rho_.get(), energy_.get(), electrostat_.get(),
+                this, proj_matrices_.get(), current_orbitals_));
     else
 #endif
-        dm_strategy_ = DMStrategyFactory<OrbitalsType,
-            dist_matrix::DistMatrix<double>>::create(comm_, os_, *ions_, rho_,
-            energy_, electrostat_, this, proj_matrices_, current_orbitals_);
+        dm_strategy_.reset(DMStrategyFactory<OrbitalsType,
+            dist_matrix::DistMatrix<double>>::create(comm_, os_, *ions_,
+            rho_.get(), energy_.get(), electrostat_.get(), this,
+            proj_matrices_.get(), current_orbitals_));
 
     // theta = invB * Hij
     proj_matrices_->updateThetaAndHB();
@@ -609,7 +572,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_, *current_orbitals_, 2, os_);
+        ts, proj_matrices_.get(), *current_orbitals_, 2, os_);
 }
 
 template <class OrbitalsType>
@@ -624,7 +587,7 @@ void MGmol<OrbitalsType>::printMM()
         ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>* projmatrices
             = dynamic_cast<
                 ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>*>(
-                proj_matrices_);
+                proj_matrices_.get());
         assert(projmatrices != nullptr);
         projmatrices->printHamiltonianMM(tfileh);
     }
@@ -697,7 +660,7 @@ void MGmol<OrbitalsType>::write_header()
         ct.printPoissonOptions(os_);
     } // onpe0
 
-    if (current_orbitals_ != nullptr && ct.verbose > 0)
+    if (current_orbitals_ && ct.verbose > 0)
     {
         current_orbitals_->printNumst(os_);
         current_orbitals_->printChromaticNumber(os_);
@@ -802,9 +765,9 @@ void MGmol<OrbitalsType>::printEigAndOcc()
 #ifdef HAVE_MAGMA
         // try with ReplicatedMatrix first
         {
-            ProjectedMatrices<ReplicatedMatrix>* projmatrices
-                = dynamic_cast<ProjectedMatrices<ReplicatedMatrix>*>(
-                    proj_matrices_);
+            std::shared_ptr<ProjectedMatrices<ReplicatedMatrix>> projmatrices
+                = std::dynamic_pointer_cast<
+                    ProjectedMatrices<ReplicatedMatrix>>(proj_matrices_);
             if (projmatrices)
             {
                 projmatrices->printEigenvalues(os_);
@@ -815,10 +778,10 @@ void MGmol<OrbitalsType>::printEigAndOcc()
 #endif
         if (!printflag)
         {
-            ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>*
-                projmatrices
-                = dynamic_cast<
-                    ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>*>(
+            std::shared_ptr<
+                ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>
+                projmatrices = std::dynamic_pointer_cast<
+                    ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>(
                     proj_matrices_);
             assert(projmatrices);
 
@@ -1062,21 +1025,23 @@ double MGmol<OrbitalsType>::get_evnl(const Ions& ions)
     double val;
     if (ct.short_sighted)
     {
-        ProjectedMatricesSparse* projmatrices
-            = dynamic_cast<ProjectedMatricesSparse*>(proj_matrices_);
+        std::shared_ptr<ProjectedMatricesSparse> projmatrices
+            = std::dynamic_pointer_cast<ProjectedMatricesSparse>(
+                proj_matrices_);
         assert(projmatrices);
 
-        val = g_kbpsi_->getEvnl(ions, projmatrices);
+        val = g_kbpsi_->getEvnl(ions, projmatrices.get());
     }
     else
     {
-        ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>* projmatrices
-            = dynamic_cast<
-                ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>*>(
+        std::shared_ptr<
+            ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>
+            projmatrices = std::dynamic_pointer_cast<
+                ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>(
                 proj_matrices_);
         assert(projmatrices);
 
-        val = g_kbpsi_->getEvnl(ions, projmatrices);
+        val = g_kbpsi_->getEvnl(ions, projmatrices.get());
     }
 
     evnl_tm_.stop();
@@ -1102,11 +1067,11 @@ void MGmol<OrbitalsType>::setup()
         printWithTimeStamp("MGmol<OrbitalsType>::setup()...", os_);
 
     if (ct.verbose > 0) printWithTimeStamp("Setup VH...", os_);
-    electrostat_ = new Electrostatic(
-        ct.getPoissonFDtype(), ct.bcPoisson, ct.screening_const);
+    electrostat_ = std::shared_ptr<Electrostatic>(new Electrostatic(
+        ct.getPoissonFDtype(), ct.bcPoisson, ct.screening_const));
     electrostat_->setup(ct.vh_init);
 
-    rho_ = new Rho<OrbitalsType>();
+    rho_ = std::shared_ptr<Rho<OrbitalsType>>(new Rho<OrbitalsType>());
     rho_->setVerbosityLevel(ct.verbose);
 
 #ifdef HAVE_MAGMA
@@ -1205,17 +1170,18 @@ template <class OrbitalsType>
 void MGmol<OrbitalsType>::setGamma(
     const pb::Lap<ORBDTYPE>& lapOper, const Potentials& pot)
 {
-    assert(orbitals_precond_ != nullptr);
+    assert(orbitals_precond_);
 
     Control& ct = *(Control::instance());
 
-    orbitals_precond_->setGamma(lapOper, pot, ct.getMGlevels(), proj_matrices_);
+    orbitals_precond_->setGamma(
+        lapOper, pot, ct.getMGlevels(), proj_matrices_.get());
 }
 
 template <class OrbitalsType>
 void MGmol<OrbitalsType>::precond_mg(OrbitalsType& phi)
 {
-    assert(orbitals_precond_ != nullptr);
+    assert(orbitals_precond_);
 
     orbitals_precond_->precond_mg(phi);
 }
@@ -1446,11 +1412,44 @@ template <class OrbitalsType>
 void MGmol<OrbitalsType>::addResidualSpreadPenalty(
     OrbitalsType& phi, OrbitalsType& res)
 {
-    assert(spread_penalty_ != nullptr);
+    assert(spread_penalty_);
 
     spread_penalty_->addResidual(phi, res);
 }
 
+template <class OrbitalsType>
+void MGmol<OrbitalsType>::getAtomicPositions(std::vector<double>& tau)
+{
+    ions_->getPositions(tau);
+}
+
+template <class OrbitalsType>
+void MGmol<OrbitalsType>::getAtomicNumbers(std::vector<short>& an)
+{
+    ions_->getAtomicNumbers(an);
+}
+
+template <class OrbitalsType>
+double MGmol<OrbitalsType>::evaluateEnergyAndForces(
+    const std::vector<double>& tau, std::vector<short>& atnumbers,
+    std::vector<double>& forces)
+{
+    assert(tau.size() == 3 * atnumbers.size());
+
+    Control& ct = *(Control::instance());
+
+    ions_->setPositions(tau, atnumbers);
+
+    double eks = 0.;
+    quench(current_orbitals_, *ions_, ct.max_electronic_steps, 20, eks);
+
+    force(*current_orbitals_, *ions_);
+
+    ions_->getForces(forces);
+
+    return eks;
+}
+
 template class MGmol<LocGridOrbitals>;
 template class MGmol<ExtendedGridOrbitals>;
 template int MGmol<LocGridOrbitals>::initial<MemorySpace::Host>();
diff --git a/src/MGmol.h b/src/MGmol.h
index 1ed04041..ef7a020b 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -44,6 +44,7 @@ class IonicAlgorithm;
 #include "AOMMprojector.h"
 #include "ClusterOrbitals.h"
 #include "DMStrategy.h"
+#include "Energy.h"
 #include "ExtendedGridOrbitals.h"
 #include "Forces.h"
 #include "Ions.h"
@@ -55,6 +56,8 @@ class IonicAlgorithm;
 #include "SpreadPenaltyInterface.h"
 #include "SpreadsAndCenters.h"
 
+#include <memory>
+
 template <class OrbitalsType>
 class MGmol : public MGmolInterface
 {
@@ -63,48 +66,49 @@ class MGmol : public MGmolInterface
 
     MPI_Comm comm_;
 
-    XConGrid* xcongrid_;
+    std::shared_ptr<XConGrid> xcongrid_;
 
     OrbitalsType* current_orbitals_;
 
-    AOMMprojector* aomm_;
+    std::shared_ptr<AOMMprojector> aomm_;
 
-    Ions* ions_;
+    std::shared_ptr<Ions> ions_;
 
-    Rho<OrbitalsType>* rho_;
+    std::shared_ptr<Rho<OrbitalsType>> rho_;
 
-    Energy<OrbitalsType>* energy_;
+    std::shared_ptr<Energy<OrbitalsType>> energy_;
 
-    Hamiltonian<OrbitalsType>* hamiltonian_;
+    std::shared_ptr<Hamiltonian<OrbitalsType>> hamiltonian_;
 
-    Forces<OrbitalsType>* forces_;
+    std::shared_ptr<Forces<OrbitalsType>> forces_;
 
-    MasksSet* currentMasks_;
-    MasksSet* corrMasks_;
+    std::shared_ptr<MasksSet> currentMasks_;
+    std::shared_ptr<MasksSet> corrMasks_;
 
-    // ProjectedMatrices* proj_matrices_;
-    ProjectedMatricesInterface* proj_matrices_;
+    std::shared_ptr<ProjectedMatricesInterface> proj_matrices_;
 
-    IonicAlgorithm<OrbitalsType>* geom_optimizer_;
+    std::shared_ptr<IonicAlgorithm<OrbitalsType>> geom_optimizer_;
 
     std::shared_ptr<LocalizationRegions> lrs_;
 
-    ClusterOrbitals* local_cluster_;
+    std::shared_ptr<ClusterOrbitals> local_cluster_;
+
+    std::shared_ptr<KBPsiMatrixSparse> g_kbpsi_;
 
-    KBPsiMatrixSparse* g_kbpsi_;
+    std::shared_ptr<SpreadsAndCenters<OrbitalsType>> spreadf_;
 
-    SpreadsAndCenters<OrbitalsType>* spreadf_;
+    std::shared_ptr<SpreadPenaltyInterface<OrbitalsType>> spread_penalty_;
 
-    SpreadPenaltyInterface<OrbitalsType>* spread_penalty_;
+    std::shared_ptr<DMStrategy<OrbitalsType>> dm_strategy_;
 
-    DMStrategy<OrbitalsType>* dm_strategy_;
+    std::shared_ptr<HDFrestart> h5f_file_;
 
-    HDFrestart* h5f_file_;
+    std::shared_ptr<OrbitalsPreconditioning<OrbitalsType>> orbitals_precond_;
 
     double total_energy_;
-    ConstraintSet* constraints_;
+    std::shared_ptr<ConstraintSet> constraints_;
 
-    OrbitalsExtrapolation<OrbitalsType>* orbitals_extrapol_ = nullptr;
+    std::shared_ptr<OrbitalsExtrapolation<OrbitalsType>> orbitals_extrapol_;
 
     float md_time_;
     int md_iteration_;
@@ -168,17 +172,30 @@ class MGmol : public MGmolInterface
     static Timer comp_res_tm_;
     static Timer init_nuc_tm_;
 
-    OrbitalsPreconditioning<OrbitalsType>* orbitals_precond_;
-
 public:
-    Electrostatic* electrostat_;
+    std::shared_ptr<Electrostatic> electrostat_;
 
     MGmol(MPI_Comm comm, std::ostream& os, std::string input_filename,
         std::string lrs_filename, std::string constraints_filename);
 
     ~MGmol() override;
 
+    /* access functions */
+    OrbitalsType* getOrbitals() { return current_orbitals_; }
+    std::shared_ptr<Hamiltonian<OrbitalsType>> getHamiltonian() { return hamiltonian_; }
+
     void run() override;
+
+    double evaluateEnergyAndForces(const std::vector<double>& tau,
+        std::vector<short>& atnumbers, std::vector<double>& forces);
+
+    /*
+     * get internal atomic positions
+     */
+    void getAtomicPositions(std::vector<double>& tau);
+
+    void getAtomicNumbers(std::vector<short>& an);
+
     void initNuc(Ions& ions);
     void initKBR();
 
@@ -276,10 +293,6 @@ class MGmol : public MGmolInterface
     double get_evnl(const Ions& ions);
     void sebprintPositions();
     void sebprintForces();
-    void get_positions(std::vector<std::vector<double>>& r);
-    void set_positions(std::vector<std::vector<double>>& r);
-    void get_forces(std::vector<std::vector<double>>& f);
-    void set_forces(std::vector<std::vector<double>>& f);
     int nions() { return ions_->getNumIons(); }
     double getTotalEnergy();
     void cleanup();
@@ -306,7 +319,7 @@ class MGmol : public MGmolInterface
         forces_->force(orbitals, ions);
     }
 
-    OrbitalsType* loadOrbitalFromRestartFile(const std::string filename);
+    void loadRestartFile(const std::string filename);
 
 #ifdef MGMOL_HAS_LIBROM
     int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals);
diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h
index 9932dd9f..046459d6 100644
--- a/src/MGmolInterface.h
+++ b/src/MGmolInterface.h
@@ -11,6 +11,7 @@
 #define MGMOLINTERFACE_H
 
 #include <cstring>
+#include <vector>
 
 class MGmolInterface
 {
@@ -24,6 +25,11 @@ class MGmolInterface
     virtual int setupConstraintsFromInput(const std::string input_file) = 0;
     virtual void setup()                                                = 0;
     virtual void run()                                                  = 0;
+    virtual double evaluateEnergyAndForces(const std::vector<double>& tau,
+        std::vector<short>& atnumbers, std::vector<double>& forces)
+        = 0;
+    virtual void getAtomicPositions(std::vector<double>& tau) = 0;
+    virtual void getAtomicNumbers(std::vector<short>& an)     = 0;
 };
 
 #endif
diff --git a/src/MGmol_prototypes.h b/src/MGmol_prototypes.h
index 207201d7..02e37890 100644
--- a/src/MGmol_prototypes.h
+++ b/src/MGmol_prototypes.h
@@ -25,7 +25,7 @@ double getLAeigen(const double tol, const int maxit, Ions& ions);
 int read_config(int argc, char** argv,
     boost::program_options::variables_map& vm, std::string& input_file,
     std::string& lrs_filename, std::string& constraints_filename,
-    float& total_spin, bool& with_spin, bool& tcheck);
+    float& total_spin, bool& with_spin);
 #ifdef MGMOL_HAS_LIBROM
 void setupROMConfigOption(po::options_description &rom_cfg);
 #endif
diff --git a/src/OrbitalsExtrapolation.h b/src/OrbitalsExtrapolation.h
index 47c1a21b..e2d0ff0c 100644
--- a/src/OrbitalsExtrapolation.h
+++ b/src/OrbitalsExtrapolation.h
@@ -26,8 +26,6 @@ class OrbitalsExtrapolation
 
     virtual ~OrbitalsExtrapolation();
 
-    virtual bool extrapolatedH() const { return false; }
-
     virtual void clearOldOrbitals();
     bool getRestartData(OrbitalsType& orbitals);
     virtual void setupPreviousOrbitals(OrbitalsType** orbitals,
diff --git a/src/Rho.cc b/src/Rho.cc
index b456f6c1..ca0f20ea 100644
--- a/src/Rho.cc
+++ b/src/Rho.cc
@@ -175,6 +175,12 @@ void Rho<OrbitalsType>::rescaleTotalCharge()
             (*MPIdata::sout) << " Rescaling factor: " << t1 << std::endl;
             (*MPIdata::sout) << " Num. electrons: " << nel << std::endl;
         }
+        if (fabs(t1 - 1.) > 0.02)
+        {
+            std::cerr << "Error on density too large to continue. Abort."
+                      << std::endl;
+            mmpi.abort();
+        }
 
         for (int ispin = 0; ispin < nspin; ispin++)
             LinearAlgebraUtils<MemorySpace::Host>::MPscal(
diff --git a/src/computeHij.cc b/src/computeHij.cc
index 4942959a..75b8ef86 100644
--- a/src/computeHij.cc
+++ b/src/computeHij.cc
@@ -280,13 +280,14 @@ template <class OrbitalsType>
 void MGmol<OrbitalsType>::getKBPsiAndHij(OrbitalsType& orbitals, Ions& ions,
     KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix<DISTMATDTYPE>& hij)
 {
-    getKBPsiAndHij(orbitals, orbitals, ions, kbpsi, proj_matrices_, hij);
+    getKBPsiAndHij(orbitals, orbitals, ions, kbpsi, proj_matrices_.get(), hij);
 }
 
 template <class OrbitalsType>
 void MGmol<OrbitalsType>::getKBPsiAndHij(OrbitalsType& orbitals, Ions& ions)
 {
-    getKBPsiAndHij(orbitals, orbitals, ions, g_kbpsi_, proj_matrices_);
+    getKBPsiAndHij(
+        orbitals, orbitals, ions, g_kbpsi_.get(), proj_matrices_.get());
 }
 
 template <class OrbitalsType>
@@ -396,7 +397,7 @@ template <class OrbitalsType>
 void MGmol<OrbitalsType>::getHpsiAndTheta(
     Ions& ions, OrbitalsType& phi, OrbitalsType& hphi)
 {
-    getHpsiAndTheta(ions, phi, hphi, g_kbpsi_);
+    getHpsiAndTheta(ions, phi, hphi, g_kbpsi_.get());
 }
 
 template <class OrbitalsType>
@@ -434,10 +435,10 @@ void MGmol<OrbitalsType>::getHpsiAndTheta(Ions& ions, OrbitalsType& phi,
         proj_matrices_->clearSparseH();
 
         // Compute the contribution of the non-local potential into sh
-        kbpsi->computeHvnlMatrix(ions, proj_matrices_);
+        kbpsi->computeHvnlMatrix(ions, proj_matrices_.get());
 
         // add local part of H to sh
-        hamiltonian_->addHlocalij(phi, proj_matrices_);
+        hamiltonian_->addHlocalij(phi, proj_matrices_.get());
 
         energy_->saveVofRho();
 
diff --git a/src/lbfgsrlx.cc b/src/lbfgsrlx.cc
index 3881abab..b2aff795 100644
--- a/src/lbfgsrlx.cc
+++ b/src/lbfgsrlx.cc
@@ -36,15 +36,14 @@ void MGmol<OrbitalsType>::lbfgsrlx(OrbitalsType** orbitals, Ions& ions)
     Control& ct = *(Control::instance());
 
     LBFGS<OrbitalsType> lbfgs(orbitals, ions, *rho_, *constraints_, lrs_,
-        local_cluster_, *currentMasks_, *corrMasks_, *electrostat_, ct.dt,
+        local_cluster_.get(), *currentMasks_, *corrMasks_, *electrostat_, ct.dt,
         *this);
 
     DFTsolver<OrbitalsType>::resetItCount();
 
-    lbfgs.init(h5f_file_);
+    lbfgs.init(h5f_file_.get());
 
-    delete h5f_file_;
-    h5f_file_ = nullptr;
+    h5f_file_.reset();
 
     // additional quench to compensate random start
     if (ct.restart_info < 3)
diff --git a/src/main.cc b/src/main.cc
index f6f396d2..3fcbf791 100644
--- a/src/main.cc
+++ b/src/main.cc
@@ -29,39 +29,16 @@
 #include "MGmol.h"
 #include "MGmol_MPI.h"
 #include "MPIdata.h"
-#include "Mesh.h"
 #include "mgmol_run.h"
-#include "tools.h"
 
 #include <cassert>
-#include <fenv.h>
 #include <iostream>
-#include <iterator>
-#include <sys/cdefs.h>
 #include <time.h>
 #include <vector>
 
-#include <mpi.h>
-
 #include <boost/program_options.hpp>
 namespace po = boost::program_options;
 
-//#include "MemTrack.h"
-
-/*
-void trapfpe () {
-  feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
-}
-*/
-
-// A helper function
-template <class T>
-std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
-{
-    copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " "));
-    return os;
-}
-
 int main(int argc, char** argv)
 {
     int mpirc = MPI_Init(&argc, &argv);
@@ -73,33 +50,37 @@ int main(int argc, char** argv)
 
     MPI_Comm comm = MPI_COMM_WORLD;
 
+    /*
+     * Initialize general things, like magma, openmp, IO, ...
+     */
     mgmol_init(comm);
 
-    // read runtime parameters
+    /*
+     * read runtime parameters
+     */
     std::string input_filename("");
     std::string lrs_filename;
     std::string constraints_filename("");
-    bool tcheck = false;
 
     float total_spin = 0.;
     bool with_spin   = false;
 
     po::variables_map vm;
 
-    // use configure file if it can be found
-    // std::string config_filename("mgmol.cfg");
-
-    // read options from PE0 only
+    // read from PE0 only
     if (MPIdata::onpe0)
     {
         read_config(argc, argv, vm, input_filename, lrs_filename,
-            constraints_filename, total_spin, with_spin, tcheck);
+            constraints_filename, total_spin, with_spin);
     }
 
     MGmol_MPI::setup(comm, std::cout, with_spin);
     MGmol_MPI& mmpi      = *(MGmol_MPI::instance());
     MPI_Comm global_comm = mmpi.commGlobal();
 
+    /*
+     * Setup control struct with run time parameters
+     */
     Control::setup(global_comm, with_spin, total_spin);
     Control& ct = *(Control::instance());
 
@@ -108,11 +89,6 @@ int main(int argc, char** argv)
     int ret = ct.checkOptions();
     if (ret < 0) return ret;
 
-    unsigned ngpts[3]    = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] };
-    double origin[3]     = { ct.ox_, ct.oy_, ct.oz_ };
-    const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ };
-    Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type);
-
     mmpi.bcastGlobal(input_filename);
     mmpi.bcastGlobal(lrs_filename);
 
@@ -126,14 +102,10 @@ int main(int argc, char** argv)
             mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
                 input_filename, lrs_filename, constraints_filename);
 
-        if (!tcheck)
-        {
-            mgmol->run();
-        }
-        else
-        {
-            *MPIdata::sout << " Input parameters OK\n";
-        }
+        mgmol->setup();
+
+        mgmol->run();
+
         delete mgmol;
 
     } // close main scope
@@ -150,9 +122,5 @@ int main(int argc, char** argv)
     time(&tt);
     if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl;
 
-    // MemTrack::TrackDumpBlocks();
-
-    //    MemTrack::TrackListMemoryUsage();
-
     return 0;
 }
diff --git a/src/md.cc b/src/md.cc
index 089deb82..939c6e62 100644
--- a/src/md.cc
+++ b/src/md.cc
@@ -91,11 +91,6 @@ void MGmol<OrbitalsType>::postWFextrapolation(OrbitalsType* orbitals)
         else
             proj_matrices_->setDMto2InvS();
     }
-    else
-    {
-        if (orbitals_extrapol_->extrapolatedH())
-            dm_strategy_->update(*orbitals);
-    }
 }
 
 template <class OrbitalsType>
@@ -137,8 +132,8 @@ OrbitalsType* MGmol<OrbitalsType>::new_orbitals_with_current_LRs(bool setup)
 
     // need to build new orbitals as masks have changed
     OrbitalsType* new_orbitals = new OrbitalsType("NewMasks", mygrid,
-        mymesh->subdivx(), ct.numst, ct.bcWF, proj_matrices_, lrs_,
-        currentMasks_, corrMasks_, local_cluster_, setup);
+        mymesh->subdivx(), ct.numst, ct.bcWF, proj_matrices_.get(), lrs_,
+        currentMasks_.get(), corrMasks_.get(), local_cluster_.get(), setup);
 
     return new_orbitals;
 }
@@ -317,8 +312,8 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
     int size_tau = (int)tau0.size();
     DFTsolver<OrbitalsType>::resetItCount();
 
-    orbitals_extrapol_ = OrbitalsExtrapolationFactory<OrbitalsType>::create(
-        ct.WFExtrapolation());
+    orbitals_extrapol_.reset(OrbitalsExtrapolationFactory<OrbitalsType>::create(
+        ct.WFExtrapolation()));
 
     MD_IonicStepper* stepper = new MD_IonicStepper(
         ct.dt, atmove, tau0, taup, taum, fion, pmass, rand_states);
@@ -349,8 +344,9 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         constraints_->enforceConstraints(20);
 
         stepper->updateTau();
-        ions.setPositions(tau0);
+        ions.setLocalPositions(tau0);
 
+        // setup required after updating local ions positions
         ions.setup();
     }
 
@@ -376,8 +372,8 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
                 if (onpe0) os_ << "Create new orbitals_minus1..." << std::endl;
 
                 orbitals_extrapol_->setupPreviousOrbitals(&current_orbitals_,
-                    proj_matrices_, lrs_, local_cluster_, currentMasks_,
-                    corrMasks_, *h5f_file_);
+                    proj_matrices_.get(), lrs_, local_cluster_.get(),
+                    currentMasks_.get(), corrMasks_.get(), *h5f_file_);
 
                 // need to reset a few things as we just read new orbitals
                 (*orbitals)->computeGramAndInvS();
@@ -395,8 +391,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
             moveVnuc(ions);
         }
 
-        delete h5f_file_;
-        h5f_file_ = nullptr;
+        h5f_file_.reset();
     }
 
     // additional SC steps to compensate random start
@@ -437,7 +432,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
             if (ct.adaptiveLRs())
             {
                 assert(lrs_);
-                adaptLR(spreadf_, nullptr);
+                adaptLR(spreadf_.get(), nullptr);
 
                 last_move_is_small = lrs_->moveIsSmall();
 
@@ -495,7 +490,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         force(**orbitals, ions);
 
         // set fion
-        ions.getForces(fion);
+        ions.getLocalForces(fion);
 
         // constraints need to be added again and setup as atoms
         // may have moved and local atoms are not the same anymore
@@ -682,11 +677,11 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
     }
 
     delete stepper;
-    delete orbitals_extrapol_;
+    orbitals_extrapol_.reset();
 }
 
 template <class OrbitalsType>
-OrbitalsType* MGmol<OrbitalsType>::loadOrbitalFromRestartFile(
+void MGmol<OrbitalsType>::loadRestartFile(
     const std::string filename)
 {
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
@@ -700,53 +695,15 @@ OrbitalsType* MGmol<OrbitalsType>::loadOrbitalFromRestartFile(
            || (ct.AtomsDynamic() == AtomsDynamicType::Quench));
 
     HDFrestart h5file(filename, myPEenv, ct.out_restart_file_type);
-    int ierr;
-
-    OrbitalsType* restart_orbitals = nullptr;
-
-    /* This corresponds to MGmol<OrbitalsType>::initial */
-    {
-        // Copy from current orbital, instead of constructing brand-new one
-        restart_orbitals
-            = new OrbitalsType("ForLoading", *current_orbitals_, false);
 
-        /* This corresponds to MGmol<OrbitalsType>::read_restart_data */
-        {
-            ierr = restart_orbitals->read_func_hdf5(h5file);
-        } // read_restart_data
-    } // initial()
-
-    /* This corresponds to MGmol<OrbitalsType>::md */
+    int ierr = read_restart_data(h5file, *rho_, *current_orbitals_);
+    if (ierr < 0)
     {
-        int flag_extrapolated_data = 0;
         if (onpe0)
-        {
-            flag_extrapolated_data
-                = h5file.dset_exists("ExtrapolatedFunction0000");
-            if (flag_extrapolated_data == 0)
-                flag_extrapolated_data
-                    = h5file.dset_exists("ExtrapolatedFunction0");
-        }
-        MPI_Bcast(&flag_extrapolated_data, 1, MPI_INT, 0, comm_);
-
-        /*
-            If extrapolated function exists,
-            then function is set as previous orbitals,
-            while the extrapolated function is set as the current orbitals.
-            This is how the restart file is saved via dumprestartFile.
+            (*MPIdata::serr) << "loadRestartFile: failed to read the restart file." << std::endl;
 
-            For now, we just enforce to not use the restart files with extrapolation.
-        */
-        if (flag_extrapolated_data)
-        {
-            if (onpe0)
-                (*MPIdata::serr) << "loadRestartFile: does not support restart files with extrapolation." << std::endl;
-
-            global_exit(0);
-        }
-
-        /* main workflow delete h5f_file_ here, meaning the loading is over. */
-    } // md()
+        global_exit(0);
+    }
 
     ierr = h5file.close();
     mmpi.allreduce(&ierr, 1, MPI_MIN);
@@ -755,20 +712,10 @@ OrbitalsType* MGmol<OrbitalsType>::loadOrbitalFromRestartFile(
         if (onpe0)
             (*MPIdata::serr)
                 << "loadRestartFile: cannot close file..." << std::endl;
-        return nullptr;
+        return;
     }
 
-    /*
-        In returning restart_orbitals,
-        we hope that the wavefunctions in restart_orbitals are all set.
-        At least the following functions should return proper data loaded from
-       the file:
-
-            restart_orbitals->getLocNumpt()
-            restart_orbitals->chromatic_number()
-            restart_orbitals->getPsi(idx)   (for int idx)
-    */
-    return restart_orbitals;
+    return;
 }
 
 template class MGmol<LocGridOrbitals>;
diff --git a/src/mlwf.cc b/src/mlwf.cc
index cc4fdeb8..3efdfeb2 100644
--- a/src/mlwf.cc
+++ b/src/mlwf.cc
@@ -279,9 +279,9 @@ int MGmol<OrbitalsType>::get_NOLMO(NOLMOTransform& noot, OrbitalsType& orbitals,
     }
     sincos.clear();
 
-    ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>* projmatrices
-        = dynamic_cast<
-            ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>*>(
+    std::shared_ptr<ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>
+        projmatrices = std::dynamic_pointer_cast<
+            ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>>(
             proj_matrices_);
     assert(projmatrices);
 
diff --git a/src/quench.cc b/src/quench.cc
index de01dad3..614aa636 100644
--- a/src/quench.cc
+++ b/src/quench.cc
@@ -326,8 +326,9 @@ bool MGmol<OrbitalsType>::rotateStatesPairsOverlap(
         spreadf2st.print(os_);
         orbitals.computeGramAndInvS();
 
-        ProjectedMatricesSparse* projmatrices
-            = dynamic_cast<ProjectedMatricesSparse*>(proj_matrices_);
+        std::shared_ptr<ProjectedMatricesSparse> projmatrices
+            = std::dynamic_pointer_cast<ProjectedMatricesSparse>(
+                proj_matrices_);
         if (onpe0)
             std::cout << "Gram Matrix for pair of states after transformation:"
                       << std::endl;
@@ -398,7 +399,7 @@ void MGmol<OrbitalsType>::disentangleOrbitals(OrbitalsType& orbitals,
 template <>
 void MGmol<LocGridOrbitals>::applyAOMMprojection(LocGridOrbitals& orbitals)
 {
-    aomm_ = new AOMMprojector(orbitals, lrs_);
+    aomm_.reset(new AOMMprojector(orbitals, lrs_));
     aomm_->projectOut(orbitals);
 }
 
@@ -422,8 +423,9 @@ int MGmol<LocGridOrbitals>::outerSolve(LocGridOrbitals& orbitals,
         case OuterSolverType::ABPG:
         case OuterSolverType::NLCG:
         {
-            DFTsolver<LocGridOrbitals> solver(hamiltonian_, proj_matrices_,
-                energy_, electrostat_, this, ions, rho_, dm_strategy_, os_);
+            DFTsolver<LocGridOrbitals> solver(hamiltonian_.get(),
+                proj_matrices_.get(), energy_.get(), electrostat_.get(), this,
+                ions, rho_.get(), dm_strategy_.get(), os_);
 
             retval = solver.solve(
                 orbitals, work_orbitals, ions, max_steps, iprint, last_eks);
@@ -433,9 +435,9 @@ int MGmol<LocGridOrbitals>::outerSolve(LocGridOrbitals& orbitals,
 
         case OuterSolverType::PolakRibiere:
         {
-            PolakRibiereSolver<LocGridOrbitals> solver(hamiltonian_,
-                proj_matrices_, energy_, electrostat_, this, ions, rho_,
-                dm_strategy_, os_);
+            PolakRibiereSolver<LocGridOrbitals> solver(hamiltonian_.get(),
+                proj_matrices_.get(), energy_.get(), electrostat_.get(), this,
+                ions, rho_.get(), dm_strategy_.get(), os_);
 
             retval = solver.solve(
                 orbitals, work_orbitals, ions, max_steps, iprint, last_eks);
@@ -466,8 +468,9 @@ int MGmol<OrbitalsType>::outerSolve(OrbitalsType& orbitals,
         case OuterSolverType::ABPG:
         case OuterSolverType::NLCG:
         {
-            DFTsolver<OrbitalsType> solver(hamiltonian_, proj_matrices_,
-                energy_, electrostat_, this, ions, rho_, dm_strategy_, os_);
+            DFTsolver<OrbitalsType> solver(hamiltonian_.get(),
+                proj_matrices_.get(), energy_.get(), electrostat_.get(), this,
+                ions, rho_.get(), dm_strategy_.get(), os_);
 
             retval = solver.solve(
                 orbitals, work_orbitals, ions, max_steps, iprint, last_eks);
@@ -477,9 +480,9 @@ int MGmol<OrbitalsType>::outerSolve(OrbitalsType& orbitals,
 
         case OuterSolverType::PolakRibiere:
         {
-            PolakRibiereSolver<OrbitalsType> solver(hamiltonian_,
-                proj_matrices_, energy_, electrostat_, this, ions, rho_,
-                dm_strategy_, os_);
+            PolakRibiereSolver<OrbitalsType> solver(hamiltonian_.get(),
+                proj_matrices_.get(), energy_.get(), electrostat_.get(), this,
+                ions, rho_.get(), dm_strategy_.get(), os_);
 
             retval = solver.solve(
                 orbitals, work_orbitals, ions, max_steps, iprint, last_eks);
@@ -499,8 +502,9 @@ int MGmol<OrbitalsType>::outerSolve(OrbitalsType& orbitals,
 #else
             DavidsonSolver<OrbitalsType, dist_matrix::DistMatrix<DISTMATDTYPE>>
 #endif
-                solver(os_, *ions_, hamiltonian_, rho_, energy_, electrostat_,
-                    this, gids, ct.dm_mix, with_spin);
+                solver(os_, *ions_, hamiltonian_.get(), rho_.get(),
+                    energy_.get(), electrostat_.get(), this, gids, ct.dm_mix,
+                    with_spin);
 
             retval = solver.solve(orbitals, work_orbitals);
             break;
@@ -558,9 +562,9 @@ int MGmol<OrbitalsType>::quench(OrbitalsType* orbitals, Ions& ions,
         applyAOMMprojection(*orbitals);
     }
 
-    orbitals_precond_ = new OrbitalsPreconditioning<OrbitalsType>();
+    orbitals_precond_.reset(new OrbitalsPreconditioning<OrbitalsType>());
     orbitals_precond_->setup(
-        *orbitals, ct.getMGlevels(), ct.lap_type, currentMasks_, lrs_);
+        *orbitals, ct.getMGlevels(), ct.lap_type, currentMasks_.get(), lrs_);
 
     // solve electronic structure problem
     // (inner iterations)
@@ -570,11 +574,9 @@ int MGmol<OrbitalsType>::quench(OrbitalsType* orbitals, Ions& ions,
 
     if (ct.use_kernel_functions)
     {
-        delete aomm_;
-        aomm_ = nullptr;
+        aomm_.reset();
     }
-    delete orbitals_precond_;
-    orbitals_precond_ = nullptr;
+    orbitals_precond_.reset();
 
     // Get the n.l. energy
     // TODO: Fix bug where energy vs. time output is incorrect if get_evnl is
@@ -599,7 +601,8 @@ int MGmol<OrbitalsType>::quench(OrbitalsType* orbitals, Ions& ions,
                 << " TS             [Ha] = " << ts << std::endl;
         }
     }
-    last_eks = energy_->evaluateTotal(ts, proj_matrices_, *orbitals, 2, os_);
+    last_eks
+        = energy_->evaluateTotal(ts, proj_matrices_.get(), *orbitals, 2, os_);
 
     if (ct.computeCondGramMD())
     {
diff --git a/src/readInput.cc b/src/readInput.cc
index b0a47411..79fbecf1 100644
--- a/src/readInput.cc
+++ b/src/readInput.cc
@@ -175,7 +175,7 @@ int MGmol<OrbitalsType>::readCoordinates(
 
     // setup ions
     const std::vector<Species>& sp(ct.getSpecies());
-    ions_ = new Ions(lattice, sp);
+    ions_.reset(new Ions(lattice, sp));
 
     if (ct.restart_info > 0
         && ct.override_restart == 0) // read restart ionic positions
@@ -217,7 +217,7 @@ int MGmol<OrbitalsType>::readCoordinates(
 
     // setup ions
     const std::vector<Species>& sp(ct.getSpecies());
-    ions_ = new Ions(lattice, sp);
+    ions_.reset(new Ions(lattice, sp));
 
     if (ct.restart_info > 0
         && ct.override_restart == 0) // read restart ionic positions
diff --git a/src/read_config.cc b/src/read_config.cc
index 6ce97648..74533dfb 100644
--- a/src/read_config.cc
+++ b/src/read_config.cc
@@ -21,8 +21,7 @@ namespace po = boost::program_options;
 
 int read_config(int argc, char** argv, po::variables_map& vm,
     std::string& input_file, std::string& lrs_filename,
-    std::string& constraints_filename, float& total_spin, bool& with_spin,
-    bool& tcheck)
+    std::string& constraints_filename, float& total_spin, bool& with_spin)
 {
     // use configure file if it can be found
     // std::string config_filename("mgmol.cfg");
@@ -378,10 +377,6 @@ int read_config(int argc, char** argv, po::variables_map& vm,
             return 0;
         }
 
-        if (vm.count("check"))
-        {
-            tcheck = true;
-        }
         if (vm.count("spin"))
         {
             total_spin = vm["spin"].as<float>();
diff --git a/src/rom_main.cc b/src/rom_main.cc
index dc877a58..3718078f 100644
--- a/src/rom_main.cc
+++ b/src/rom_main.cc
@@ -50,7 +50,6 @@ int main(int argc, char** argv)
     std::string input_filename("");
     std::string lrs_filename;
     std::string constraints_filename("");
-    bool tcheck = false;
 
     float total_spin = 0.;
     bool with_spin   = false;
@@ -64,7 +63,7 @@ int main(int argc, char** argv)
     if (MPIdata::onpe0)
     {
         read_config(argc, argv, vm, input_filename, lrs_filename,
-            constraints_filename, total_spin, with_spin, tcheck);
+            constraints_filename, total_spin, with_spin);
     }
 
     MGmol_MPI::setup(comm, std::cout, with_spin);
@@ -79,6 +78,11 @@ int main(int argc, char** argv)
     int ret = ct.checkOptions();
     if (ret < 0) return ret;
 
+    unsigned ngpts[3]    = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] };
+    double origin[3]     = { ct.ox_, ct.oy_, ct.oz_ };
+    const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ };
+    Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type);
+
     mmpi.bcastGlobal(input_filename);
     mmpi.bcastGlobal(lrs_filename);
 
@@ -86,37 +90,19 @@ int main(int argc, char** argv)
     {
         MGmolInterface* mgmol;
         if (ct.isLocMode())
-            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout);
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
         else
-            mgmol
-                = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout);
-
-        unsigned ngpts[3]    = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] };
-        double origin[3]     = { ct.ox_, ct.oy_, ct.oz_ };
-        const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ };
-        Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type);
-
-        mgmol->setupFromInput(input_filename);
-
-        if (ct.isLocMode() || ct.init_loc == 1) mgmol->setupLRs(lrs_filename);
+            mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
 
-        mgmol->setupConstraintsFromInput(constraints_filename);
+        mgmol->setup();
 
-        mgmol_setup();
-
-        if (!tcheck)
-        {
-            mgmol->setup();
-
-            if (ct.isLocMode())
-                readRestartFiles<LocGridOrbitals>(mgmol);
-            else
-                readRestartFiles<ExtendedGridOrbitals>(mgmol);
-        }
+        if (ct.isLocMode())
+            readRestartFiles<LocGridOrbitals>(mgmol);
         else
-        {
-            *MPIdata::sout << " Input parameters OK\n";
-        }
+            readRestartFiles<ExtendedGridOrbitals>(mgmol);
+
         delete mgmol;
 
     } // close main scope
diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc
index 4594dd3a..1788281d 100644
--- a/src/rom_workflows.cc
+++ b/src/rom_workflows.cc
@@ -40,7 +40,7 @@ void readRestartFiles(MGmolInterface *mgmol_)
 
     /* Read the first snapshot to determin dimension and number of snapshots */
     filename = string_format(rom_options.restart_file_fmt, minidx);
-    orbitals = mgmol->loadOrbitalFromRestartFile(filename);
+    mgmol->loadRestartFile(filename);
     const int dim = orbitals->getLocNumpt();
     const int chrom_num = orbitals->chromatic_number();
     const int totalSamples = orbitals->chromatic_number() * num_restart;
@@ -54,7 +54,7 @@ void readRestartFiles(MGmolInterface *mgmol_)
     for (int k = minidx; k <= maxidx; k++)
     {
         filename = string_format(rom_options.restart_file_fmt, k);
-        orbitals = mgmol->loadOrbitalFromRestartFile(filename);
+        mgmol->loadRestartFile(filename);
         assert(dim == orbitals->getLocNumpt());
         assert(chrom_num == orbitals->chromatic_number());
 
diff --git a/src/runfire.cc b/src/runfire.cc
index d71c52a8..3239e1b1 100644
--- a/src/runfire.cc
+++ b/src/runfire.cc
@@ -33,13 +33,12 @@ void MGmol<OrbitalsType>::runfire(OrbitalsType** orbitals, Ions& ions)
 
     DFTsolver<OrbitalsType>::resetItCount();
 
-    orbitals_extrapol_ = OrbitalsExtrapolationFactory<OrbitalsType>::create(
-        ct.WFExtrapolation());
+    orbitals_extrapol_.reset(OrbitalsExtrapolationFactory<OrbitalsType>::create(
+        ct.WFExtrapolation()));
 
-    fire.init(h5f_file_);
+    fire.init(h5f_file_.get());
 
-    delete h5f_file_;
-    h5f_file_ = nullptr;
+    h5f_file_.reset();
 
     // additional quench to compensate random start
     if (ct.restart_info < 3)
@@ -129,7 +128,7 @@ void MGmol<OrbitalsType>::runfire(OrbitalsType** orbitals, Ions& ions)
 
     } // end for steps
 
-    delete orbitals_extrapol_;
+    orbitals_extrapol_.reset();
 
     // final dump
     if (ct.out_restart_info > 0)
diff --git a/src/setup.cc b/src/setup.cc
index 62d9f835..f022678c 100644
--- a/src/setup.cc
+++ b/src/setup.cc
@@ -28,7 +28,15 @@ int MGmol<OrbitalsType>::setupFromInput(const std::string filename)
 
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
 
-    hamiltonian_    = new Hamiltonian<OrbitalsType>();
+    /*
+     * Setup global mesh for calculations
+     */
+    unsigned ngpts[3]    = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] };
+    double origin[3]     = { ct.ox_, ct.oy_, ct.oz_ };
+    const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ };
+    Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type);
+
+    hamiltonian_.reset(new Hamiltonian<OrbitalsType>());
     Potentials& pot = hamiltonian_->potential();
 
     ct.registerPotentials(pot);
@@ -43,8 +51,8 @@ int MGmol<OrbitalsType>::setupFromInput(const std::string filename)
 
     const pb::PEenv& myPEenv = mymesh->peenv();
     if (ct.restart_info > 0)
-        h5f_file_
-            = new HDFrestart(ct.restart_file, myPEenv, ct.restart_file_type);
+        h5f_file_.reset(
+            new HDFrestart(ct.restart_file, myPEenv, ct.restart_file_type));
 
     int status = readCoordinates(filename, false);
     if (status == -1) return -1;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 4c68adb7..fe7b571b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -482,7 +482,7 @@ add_test(NAME ChebyshevMVP
          ${CMAKE_CURRENT_SOURCE_DIR}/Chebyshev/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 
-if(NOT ${MAGMA_FOUND})
+if(NOT ${MGMOL_WITH_MAGMA})
   add_test(NAME testShortSighted
            COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/ShortSighted/test.py
            ${MPIEXEC} --oversubscribe ${MPIEXEC_NUMPROC_FLAG} 5 ${MPIEXEC_PREFLAGS}
diff --git a/tests/SiH4/testSiH4.py b/tests/SiH4/testSiH4.py
index 1a0965ac..c7f81bdc 100755
--- a/tests/SiH4/testSiH4.py
+++ b/tests/SiH4/testSiH4.py
@@ -44,16 +44,21 @@
 lines=output.split(b'\n')
 
 tol = 5.e-3
+found_forces = False
 for line in lines:
-  num_matches = line.count(b'%%')
-  if num_matches:
+  if line.count(b'%%'):
     print(line)
-  num_matches = line.count(b'##')
-  if num_matches:
+  if line.count(b'##'):
     words=line.split()
     if len(words)==8:
       print(line)
+      found_forces = True
       for i in range(5,8):
         if abs(eval(words[i]))>tol:
           sys.exit(1)
 
+if (not found_forces):
+  print("no forces found")
+  sys.exit(1)
+
+sys.exit(0)
diff --git a/util/xyz2in.py b/util/xyz2in.py
index 3555e138..bb2e022d 100644
--- a/util/xyz2in.py
+++ b/util/xyz2in.py
@@ -5,9 +5,11 @@
 # This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 # Please also read this link https://github.com/llnl/mgmol/LICENSE
 #
-# Python program to convert .xyz file into mgmol input
+# Python program to convert .xyz file into mgmol input coordinates file
+# Optional arguments: [lx ly lz] to define box size and map all atoms into
+#                     box (0,0,0)-(lx,ly,lz) using periodic boundary conditions
 #
-# use: python coords.xyz > coords.in
+# use: python coords.xyz [lx ly lz] > coords.in
 #-------------------------------------------------------------------------------
 import sys, string
 
@@ -15,6 +17,16 @@
 
 #read file
 ifile=open(sys.argv[1],'r')
+
+lx = 0.
+ly = 0.
+lz = 0.
+if( len(sys.argv) > 2 ):
+  lx = eval(sys.argv[2])
+  ly = eval(sys.argv[3])
+  lz = eval(sys.argv[4])
+
+
 lines=ifile.readlines()
 
 count=0
@@ -22,13 +34,31 @@
 dummy=0   #unused flag set to 0
 
 for line in lines: ## loop over lines of file
-  words=line.split()
-  if len(words)>1:
-    if words[0][0:1]!='#':
-      name=words[0]+str(count)
+  if count>1:
+    words=line.split()
+    if len(words)>1:
+      name=words[0]+str(count-2)
       x=eval(words[1])*ang2bohr
       y=eval(words[2])*ang2bohr
       z=eval(words[3])*ang2bohr
-      
+
+      if lx > 0.:
+        if x<0:
+          x = x +lx
+        if x>lx:
+          x = x -lx
+
+      if ly > 0.:
+        if y<0:
+          y = y +ly
+        if y>ly:
+          y = y -ly
+
+      if lz > 0.:
+        if z<0:
+          z = z +lz
+        if z>lz:
+          z = z -lz
+
       print(name,'\t',dummy,'\t',x,'\t',y,'\t',z,'\t',movable)
-      count=count+1
+  count=count+1