diff --git a/src/Hamiltonian.cc b/src/Hamiltonian.cc
index 77a27093..aed9dbad 100644
--- a/src/Hamiltonian.cc
+++ b/src/Hamiltonian.cc
@@ -222,6 +222,8 @@ void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocal2matrix(
     LocGridOrbitals<ORBDTYPE>& phi1, LocGridOrbitals<ORBDTYPE>& phi2,
     ReplicatedMatrix& hij, const bool force)
 {
+    (void)phi1;
+    (void)phi2;
     (void)hij;
 
     applyLocal(phi2, force);
@@ -284,39 +286,5 @@ void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocal2matrix(
     mat.insertMatrixElements(ss, phi1.getOverlappingGids(), ct.numst);
 }
 
-template Hamiltonian<LocGridOrbitals<ORBDTYPE>>::Hamiltonian();
-template Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>::Hamiltonian();
-
-template Hamiltonian<LocGridOrbitals<ORBDTYPE>>::~Hamiltonian();
-template Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>::~Hamiltonian();
-
-template void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::setup(
-    pb::Grid const&, int);
-template void Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>::setup(
-    pb::Grid const&, int);
-
-template const LocGridOrbitals<ORBDTYPE>&
-Hamiltonian<LocGridOrbitals<ORBDTYPE>>::applyLocal(
-    LocGridOrbitals<ORBDTYPE>&, const bool);
-template const ExtendedGridOrbitals<ORBDTYPE>&
-Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>::applyLocal(
-    ExtendedGridOrbitals<ORBDTYPE>&, const bool);
-template void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocalij(
-    LocGridOrbitals<ORBDTYPE>&, LocGridOrbitals<ORBDTYPE>&,
-    ProjectedMatricesInterface* proj_matrices);
-template void Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>::addHlocalij(
-    ExtendedGridOrbitals<ORBDTYPE>&, ExtendedGridOrbitals<ORBDTYPE>&,
-    ProjectedMatricesInterface* proj_matrices);
-template void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocalij(
-    LocGridOrbitals<ORBDTYPE>&, ProjectedMatricesInterface* proj_matrices);
-template void Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>::addHlocalij(
-    ExtendedGridOrbitals<ORBDTYPE>&, ProjectedMatricesInterface* proj_matrices);
-template void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocal2matrix(
-    LocGridOrbitals<ORBDTYPE>&, LocGridOrbitals<ORBDTYPE>&,
-    VariableSizeMatrix<sparserow>& mat, const bool force);
-template void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocal2matrix(
-    LocGridOrbitals<ORBDTYPE>&, LocGridOrbitals<ORBDTYPE>&,
-    dist_matrix::DistMatrix<DISTMATDTYPE>& hij, const bool force);
-template void Hamiltonian<LocGridOrbitals<ORBDTYPE>>::addHlocal2matrix(
-    LocGridOrbitals<ORBDTYPE>&, LocGridOrbitals<ORBDTYPE>&, ReplicatedMatrix&,
-    const bool force);
+template class Hamiltonian<LocGridOrbitals<ORBDTYPE>>;
+template class Hamiltonian<ExtendedGridOrbitals<ORBDTYPE>>;
diff --git a/src/LocGridOrbitals.h b/src/LocGridOrbitals.h
index a19f42a6..a2cd525d 100644
--- a/src/LocGridOrbitals.h
+++ b/src/LocGridOrbitals.h
@@ -35,6 +35,7 @@ class ProjectedMatricesInterface;
 class LocalizationRegions;
 class MasksSet;
 class Masks4Orbitals;
+class ReplicatedMatrix;
 
 template <typename ScalarType>
 class LocGridOrbitals : public Orbitals
@@ -355,6 +356,12 @@ class LocGridOrbitals : public Orbitals
 
     void addDotWithNcol2Matrix(
         LocGridOrbitals&, dist_matrix::DistMatrix<DISTMATDTYPE>&) const;
+    void addDotWithNcol2Matrix(LocGridOrbitals&, ReplicatedMatrix&) const
+    {
+        std::cerr << "LocGridOrbitals::addDotWithNcol2Matrix not implemented "
+                     "for ReplicatedMatrix"
+                  << std::endl;
+    }
 
     void scal(const double alpha)
     {
diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc
index 3ab04651..3be7fa10 100644
--- a/src/MVPSolver.cc
+++ b/src/MVPSolver.cc
@@ -211,6 +211,8 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 
         kbpsi.computeHvnlMatrix(&kbpsi, ions_, h11_nl);
 
+        OrbitalsType hphi("MVP_hphi", orbitals);
+
         for (int inner_it = 0; inner_it < n_inner_steps_; inner_it++)
         {
             if (onpe0 && ct.verbose > 1)
@@ -239,7 +241,8 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
             // compute h11 for the current potential by adding local part to
             // nonlocal components
             MatrixType h11(h11_nl);
-            hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11, false);
+            hamiltonian_->applyLocal(numst_, orbitals, hphi);
+            orbitals.addDotWithNcol2Matrix(hphi, h11);
 
             current_proj_mat->assignH(h11);
             current_proj_mat->setHB2H();
@@ -318,8 +321,8 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 
                     // update h11
                     h11 = h11_nl;
-                    hamiltonian_->addHlocal2matrix(
-                        orbitals, orbitals, h11, false);
+                    hamiltonian_->applyLocal(numst_, orbitals, hphi);
+                    orbitals.addDotWithNcol2Matrix(hphi, h11);
 
                     proj_mat_work_->assignH(h11);
                     proj_mat_work_->setHB2H();