Update usage of Hamiltonian class

src/DMStrategyFactory.cc

@@ -7,6 +7,7 @@ DMStrategy<LocGridOrbitals>* DMStrategyFactory<LocGridOrbitals,
                                                                           comm,
     std::ostream& os, Ions& ions, Rho<LocGridOrbitals>* rho,
     Energy<LocGridOrbitals>* energy, Electrostatic* electrostat,
+    Hamiltonian<LocGridOrbitals>* hamiltonian,
     MGmol<LocGridOrbitals>* mgmol_strategy,
     ProjectedMatricesInterface* /*proj_matrices*/, LocGridOrbitals* orbitals,
     const bool short_sighted)
@@ -16,7 +17,7 @@ DMStrategy<LocGridOrbitals>* DMStrategyFactory<LocGridOrbitals,
         DMStrategy<LocGridOrbitals>* dm_strategy
             = new HamiltonianMVP_DMStrategy<VariableSizeMatrix<sparserow>,
                 ProjectedMatricesSparse, LocGridOrbitals>(comm, os, ions, rho,
-                energy, electrostat, mgmol_strategy, orbitals);
+                energy, electrostat, hamiltonian, mgmol_strategy, orbitals);
 
         return dm_strategy;
     }
@@ -27,7 +28,7 @@ DMStrategy<LocGridOrbitals>* DMStrategyFactory<LocGridOrbitals,
                 dist_matrix::DistMatrix<DISTMATDTYPE>,
                 ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>,
                 LocGridOrbitals>(comm, os, ions, rho, energy, electrostat,
-                mgmol_strategy, orbitals);
+                hamiltonian, mgmol_strategy, orbitals);
 
         return dm_strategy;
     }
@@ -38,6 +39,7 @@ DMStrategy<LocGridOrbitals>* DMStrategyFactory<LocGridOrbitals,
     ReplicatedMatrix>::createHamiltonianMVP_DMStrategy(MPI_Comm comm,
     std::ostream& /*os*/, Ions& /*ions*/, Rho<LocGridOrbitals>* /*rho*/,
     Energy<LocGridOrbitals>* /*energy*/, Electrostatic* /*electrostat*/,
+    Hamiltonian<LocGridOrbitals>* /*hamiltonian*/,
     MGmol<LocGridOrbitals>* /*mgmol_strategy*/,
     ProjectedMatricesInterface* /*proj_matrices*/,
     LocGridOrbitals* /*orbitals*/, const bool /*short_sighted*/)
@@ -55,6 +57,7 @@ DMStrategy<ExtendedGridOrbitals>* DMStrategyFactory<ExtendedGridOrbitals,
                                                                           comm,
     std::ostream& os, Ions& ions, Rho<ExtendedGridOrbitals>* rho,
     Energy<ExtendedGridOrbitals>* energy, Electrostatic* electrostat,
+    Hamiltonian<ExtendedGridOrbitals>* hamiltonian,
     MGmol<ExtendedGridOrbitals>* mgmol_strategy,
     ProjectedMatricesInterface* /*proj_matrices*/,
     ExtendedGridOrbitals* orbitals, const bool short_sighted)
@@ -64,8 +67,8 @@ DMStrategy<ExtendedGridOrbitals>* DMStrategyFactory<ExtendedGridOrbitals,
     DMStrategy<ExtendedGridOrbitals>* dm_strategy
         = new HamiltonianMVP_DMStrategy<dist_matrix::DistMatrix<DISTMATDTYPE>,
             ProjectedMatrices<dist_matrix::DistMatrix<DISTMATDTYPE>>,
-            ExtendedGridOrbitals>(
-            comm, os, ions, rho, energy, electrostat, mgmol_strategy, orbitals);
+            ExtendedGridOrbitals>(comm, os, ions, rho, energy, electrostat,
+            hamiltonian, mgmol_strategy, orbitals);
 
     return dm_strategy;
 }
@@ -75,6 +78,7 @@ DMStrategy<ExtendedGridOrbitals>* DMStrategyFactory<ExtendedGridOrbitals,
     ReplicatedMatrix>::createHamiltonianMVP_DMStrategy(MPI_Comm comm,
     std::ostream& os, Ions& ions, Rho<ExtendedGridOrbitals>* rho,
     Energy<ExtendedGridOrbitals>* energy, Electrostatic* electrostat,
+    Hamiltonian<ExtendedGridOrbitals>* hamiltonian,
     MGmol<ExtendedGridOrbitals>* mgmol_strategy,
     ProjectedMatricesInterface* /*proj_matrices*/,
     ExtendedGridOrbitals* orbitals, const bool short_sighted)
@@ -83,8 +87,9 @@ DMStrategy<ExtendedGridOrbitals>* DMStrategyFactory<ExtendedGridOrbitals,
 
     DMStrategy<ExtendedGridOrbitals>* dm_strategy
         = new HamiltonianMVP_DMStrategy<ReplicatedMatrix,
-            ProjectedMatrices<ReplicatedMatrix>, ExtendedGridOrbitals>(
-            comm, os, ions, rho, energy, electrostat, mgmol_strategy, orbitals);
+            ProjectedMatrices<ReplicatedMatrix>, ExtendedGridOrbitals>(comm, os,
+            ions, rho, energy, electrostat, hamiltonian, mgmol_strategy,
+            orbitals);
 
     return dm_strategy;
 }

src/DMStrategyFactory.h

@@ -26,7 +26,8 @@ class DMStrategyFactory
 public:
     static DMStrategy<OrbitalsType>* create(MPI_Comm comm, std::ostream& os,
         Ions& ions, Rho<OrbitalsType>* rho, Energy<OrbitalsType>* energy,
-        Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
+        Electrostatic* electrostat, Hamiltonian<OrbitalsType>* hamiltonian,
+        MGmol<OrbitalsType>* mgmol_strategy,
         ProjectedMatricesInterface* proj_matrices, OrbitalsType* orbitals)
     {
         Control& ct     = *(Control::instance());
@@ -36,14 +37,14 @@ class DMStrategyFactory
         if (ct.DM_solver() == DMNonLinearSolverType::MVP)
         {
             dm_strategy = new MVP_DMStrategy<OrbitalsType, MatrixType>(comm, os,
-                ions, rho, energy, electrostat, mgmol_strategy,
+                ions, rho, energy, electrostat, hamiltonian, mgmol_strategy,
                 orbitals->getOverlappingGids(), proj_matrices, ct.use_old_dm());
         }
         else if (ct.DM_solver() == DMNonLinearSolverType::HMVP)
         {
             dm_strategy = createHamiltonianMVP_DMStrategy(comm, os, ions, rho,
-                energy, electrostat, mgmol_strategy, proj_matrices, orbitals,
-                ct.short_sighted);
+                energy, electrostat, hamiltonian, mgmol_strategy, proj_matrices,
+                orbitals, ct.short_sighted);
         }
         else
         {
@@ -84,6 +85,7 @@ class DMStrategyFactory
     static DMStrategy<OrbitalsType>* createHamiltonianMVP_DMStrategy(
         MPI_Comm comm, std::ostream& os, Ions& ions, Rho<OrbitalsType>* rho,
         Energy<OrbitalsType>* energy, Electrostatic* electrostat,
+        Hamiltonian<OrbitalsType>* hamiltonian,
         MGmol<OrbitalsType>* mgmol_strategy,
         ProjectedMatricesInterface* proj_matrices, OrbitalsType* orbitals,
         const bool);

src/Hamiltonian.cc

@@ -216,6 +216,18 @@ void Hamiltonian<ExtendedGridOrbitals>::addHlocal2matrix(
     phi1.addDotWithNcol2Matrix(*hlphi_, hij);
 }
 
+template <>
+template <>
+void Hamiltonian<LocGridOrbitals>::addHlocal2matrix(LocGridOrbitals& phi1,
+    LocGridOrbitals& phi2, ReplicatedMatrix& hij, const bool force)
+{
+    applyLocal(phi2, force);
+
+    // phi1.addDotWithNcol2Matrix(*hlphi_, hij);
+    std::cerr << "Not implemented!" << std::endl;
+    MPI_Abort(MPI_COMM_WORLD, EXIT_FAILURE);
+}
+
 template <class T>
 void Hamiltonian<T>::addHlocalij(
     T& phi1, T& phi2, ProjectedMatricesInterface* proj_matrices)
@@ -293,3 +305,8 @@ template void Hamiltonian<ExtendedGridOrbitals>::addHlocalij(
     ExtendedGridOrbitals&, ProjectedMatricesInterface* proj_matrices);
 template void Hamiltonian<LocGridOrbitals>::addHlocal2matrix(LocGridOrbitals&,
     LocGridOrbitals&, VariableSizeMatrix<sparserow>& mat, const bool force);
+template void Hamiltonian<LocGridOrbitals>::addHlocal2matrix(LocGridOrbitals&,
+    LocGridOrbitals&, dist_matrix::DistMatrix<DISTMATDTYPE>& hij,
+    const bool force);
+template void Hamiltonian<LocGridOrbitals>::addHlocal2matrix(
+    LocGridOrbitals&, LocGridOrbitals&, ReplicatedMatrix&, const bool force);

src/HamiltonianMVPSolver.cc

@@ -37,8 +37,9 @@ template <class MatrixType, class ProjMatrixType, class OrbitalsType>
 HamiltonianMVPSolver<MatrixType, ProjMatrixType,
     OrbitalsType>::HamiltonianMVPSolver(std::ostream& os, Ions& ions,
     Rho<OrbitalsType>* rho, Energy<OrbitalsType>* energy,
-    Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
-    const int numst, const short n_inner_steps, const MatrixType& hinit,
+    Electrostatic* electrostat, Hamiltonian<OrbitalsType>* hamiltonian,
+    MGmol<OrbitalsType>* mgmol_strategy, const int numst,
+    const short n_inner_steps, const MatrixType& hinit,
     const bool try_shorter_intervals)
     : os_(os),
       n_inner_steps_(n_inner_steps),
@@ -50,6 +51,7 @@ HamiltonianMVPSolver<MatrixType, ProjMatrixType,
     rho_            = rho;
     energy_         = energy;
     electrostat_    = electrostat;
+    hamiltonian_    = hamiltonian;
     mgmol_strategy_ = mgmol_strategy;
 
     numst_ = numst;
@@ -149,7 +151,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
         // compute new h11 for the current potential by adding local part to
         // nonlocal components
         h11 = h11nl;
-        mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
+        hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
 
         projmatrices->assignH(h11);
         projmatrices->setHB2H();
@@ -177,7 +179,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
         // update H and compute energy at midpoint
         h11 = h11nl;
-        mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
+        hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
 
         projmatrices->assignH(h11);
         projmatrices->setHB2H();
@@ -212,7 +214,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
         // update H with new potential
         h11 = h11nl;
-        mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
+        hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
 
         projmatrices->assignH(h11);
         projmatrices->setHB2H();
@@ -268,7 +270,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
                 // update H
                 h11 = h11nl;
-                mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
+                hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
 
                 projmatrices->assignH(h11);
                 projmatrices->setHB2H();

src/HamiltonianMVPSolver.h

@@ -11,6 +11,7 @@
 #define MGMOL_HAMILTONIANMVP_SOLVER_H_
 
 #include "Energy.h"
+#include "Hamiltonian.h"
 #include "MGmol.h"
 #include "Rho.h"
 #include "Timer.h"
@@ -38,6 +39,7 @@ class HamiltonianMVPSolver
     Rho<OrbitalsType>* rho_;
     Energy<OrbitalsType>* energy_;
     Electrostatic* electrostat_;
+    Hamiltonian<OrbitalsType>* hamiltonian_;
     MGmol<OrbitalsType>* mgmol_strategy_;
 
     int numst_;
@@ -66,6 +68,7 @@ class HamiltonianMVPSolver
 public:
     HamiltonianMVPSolver(std::ostream& os, Ions& ions, Rho<OrbitalsType>* rho,
         Energy<OrbitalsType>* energy, Electrostatic* electrostat,
+        Hamiltonian<OrbitalsType>* hamiltonian,
         MGmol<OrbitalsType>* mgmol_strategy, const int numst,
         const short n_inner_steps, const MatrixType& hinit,
         const bool try_shorter_intervals = false);

src/HamiltonianMVP_DMStrategy.cc

@@ -22,14 +22,15 @@ template <class MatrixType, class ProjMatrixType, class OrbitalsType>
 HamiltonianMVP_DMStrategy<MatrixType, ProjMatrixType,
     OrbitalsType>::HamiltonianMVP_DMStrategy(MPI_Comm comm, std::ostream& os,
     Ions& ions, Rho<OrbitalsType>* rho, Energy<OrbitalsType>* energy,
-    Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
-    OrbitalsType* orbitals)
+    Electrostatic* electrostat, Hamiltonian<OrbitalsType>* hamiltonian,
+    MGmol<OrbitalsType>* mgmol_strategy, OrbitalsType* orbitals)
     : comm_(comm),
       os_(os),
       ions_(ions),
       rho_(rho),
       energy_(energy),
       electrostat_(electrostat),
+      hamiltonian_(hamiltonian),
       global_indexes_(orbitals->getOverlappingGids()),
       mgmol_strategy_(mgmol_strategy)
 {
@@ -44,8 +45,9 @@ HamiltonianMVP_DMStrategy<MatrixType, ProjMatrixType,
 
     solver_
         = new HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>(
-            os_, ions_, rho_, energy_, electrostat_, mgmol_strategy_, ct.numst,
-            ct.dm_inner_steps, projmatrices->getH(), true);
+            os_, ions_, rho_, energy_, electrostat_, hamiltonian_,
+            mgmol_strategy_, ct.numst, ct.dm_inner_steps, projmatrices->getH(),
+            true);
 }
 
 template <class MatrixType, class ProjMatrixType, class OrbitalsType>

src/HamiltonianMVP_DMStrategy.h

@@ -11,13 +11,13 @@
 #define MGMOL_HamiltonianMVP_DMStrategy_H
 
 #include "DMStrategy.h"
+#include "Electrostatic.h"
 #include "Energy.h"
 #include "HamiltonianMVPSolver.h"
+#include "Ions.h"
 #include "MGmol.h"
 #include "Rho.h"
 
-class Ions;
-class Electrostatic;
 template <class T>
 class MGmol;
 
@@ -32,6 +32,7 @@ class HamiltonianMVP_DMStrategy : public DMStrategy<OrbitalsType>
     Rho<OrbitalsType>* rho_;
     Energy<OrbitalsType>* energy_;
     Electrostatic* electrostat_;
+    Hamiltonian<OrbitalsType>* hamiltonian_;
     const std::vector<std::vector<int>>& global_indexes_;
     MGmol<OrbitalsType>* mgmol_strategy_;
 
@@ -40,8 +41,8 @@ class HamiltonianMVP_DMStrategy : public DMStrategy<OrbitalsType>
 public:
     HamiltonianMVP_DMStrategy(MPI_Comm comm, std::ostream& os, Ions& ions,
         Rho<OrbitalsType>* rho, Energy<OrbitalsType>* energy,
-        Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
-        OrbitalsType* orbitals);
+        Electrostatic* electrostat, Hamiltonian<OrbitalsType>* hamiltonian,
+        MGmol<OrbitalsType>* mgmol_strategy, OrbitalsType* orbitals);
 
     ~HamiltonianMVP_DMStrategy() override;
 

src/MGmol.cc

@@ -461,12 +461,13 @@ int MGmol<OrbitalsType>::initial()
         dm_strategy_.reset(
             DMStrategyFactory<OrbitalsType, ReplicatedMatrix>::create(comm_,
                 os_, *ions_, rho_.get(), energy_.get(), electrostat_.get(),
-                this, proj_matrices_.get(), current_orbitals_));
+                hamiltonian_.get(), this, proj_matrices_.get(),
+                current_orbitals_));
     else
         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_));
+            rho_.get(), energy_.get(), electrostat_.get(), hamiltonian_.get(),
+            this, proj_matrices_.get(), current_orbitals_));
 
     // theta = invB * Hij
     proj_matrices_->updateThetaAndHB();
@@ -1502,8 +1503,8 @@ double MGmol<OrbitalsType>::evaluateDMandEnergyAndForces(Orbitals* orbitals,
     {
         std::shared_ptr<DMStrategy<OrbitalsType>> dm_strategy(
             DMStrategyFactory<OrbitalsType, ReplicatedMatrix>::create(comm_,
-                os_, ions, rho_.get(), energy_.get(), electrostat_.get(), this,
-                proj_matrices_.get(), dorbitals));
+                os_, ions, rho_.get(), energy_.get(), electrostat_.get(),
+                hamiltonian_.get(), this, proj_matrices_.get(), dorbitals));
 
         dm_strategy->update(*dorbitals);
     }
@@ -1512,8 +1513,8 @@ double MGmol<OrbitalsType>::evaluateDMandEnergyAndForces(Orbitals* orbitals,
         std::shared_ptr<DMStrategy<OrbitalsType>> dm_strategy(
             DMStrategyFactory<OrbitalsType,
                 dist_matrix::DistMatrix<double>>::create(comm_, os_, ions,
-                rho_.get(), energy_.get(), electrostat_.get(), this,
-                proj_matrices_.get(), dorbitals));
+                rho_.get(), energy_.get(), electrostat_.get(),
+                hamiltonian_.get(), this, proj_matrices_.get(), dorbitals));
 
         dm_strategy->update(*dorbitals);
     }

src/MGmol.h

@@ -254,12 +254,6 @@ class MGmol : public MGmolInterface
         const Ions& ions, const KBPsiMatrixSparse* const kbpsi,
         ProjectedMatricesInterface*);
 
-    template <class MatrixType>
-    void addHlocal2matrix(
-        OrbitalsType& orbitalsi, OrbitalsType& orbitalsj, MatrixType& mat);
-    void addHlocal2matrix(OrbitalsType& orbitalsi, OrbitalsType& orbitalsj,
-        VariableSizeMatrix<SparseRow>& mat);
-
     void update_pot(const pb::GridFunc<POTDTYPE>& vh_init, const Ions& ions);
     void update_pot(const Ions& ions);
     int quench(OrbitalsType& orbitals, Ions& ions, const int max_steps,

src/MVPSolver.cc

@@ -44,8 +44,8 @@ double evalEntropyMVP(ProjectedMatricesInterface* projmatrices,
 template <class OrbitalsType, class MatrixType>
 MVPSolver<OrbitalsType, MatrixType>::MVPSolver(MPI_Comm comm, std::ostream& os,
     Ions& ions, Rho<OrbitalsType>* rho, Energy<OrbitalsType>* energy,
-    Electrostatic* electrostat, MGmol<OrbitalsType>* mgmol_strategy,
-    const int numst, const double kbT,
+    Electrostatic* electrostat, Hamiltonian<OrbitalsType>* hamiltonian,
+    MGmol<OrbitalsType>* mgmol_strategy, const int numst, const double kbT,
     const std::vector<std::vector<int>>& global_indexes,
     const short n_inner_steps, const double mixing, const double tol_de0,
     const bool use_old_dm)
@@ -69,6 +69,7 @@ MVPSolver<OrbitalsType, MatrixType>::MVPSolver(MPI_Comm comm, std::ostream& os,
     rho_            = rho;
     energy_         = energy;
     electrostat_    = electrostat;
+    hamiltonian_    = hamiltonian;
     mgmol_strategy_ = mgmol_strategy;
 
     work_ = new MatrixType("workMVP", numst_, numst_);
@@ -238,7 +239,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
             // compute h11 for the current potential by adding local part to
             // nonlocal components
             MatrixType h11(h11_nl);
-            mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
+            hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
 
             current_proj_mat->assignH(h11);
             current_proj_mat->setHB2H();
@@ -318,8 +319,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                     // update h11
                     {
                         h11 = h11_nl;
-                        mgmol_strategy_->addHlocal2matrix(
-                            orbitals, orbitals, h11);
+                        hamiltonian_->addHlocal2matrix(orbitals, orbitals, h11);
                     }
 
                     proj_mat_work_->assignH(h11);