Merge release (#332)

* Fix EnergyAndForces tests (#277)

* have them work in debug mode too

* Move factor 4pi out og linear solvers (#278)

* Move some code into PoissonSolverFactory (#279)

* Clean up class Potentials (#280)

* Clean up class Ions, add test for it (#281)

* Add test MD_MVP (#290)

* Clean up code related to DM restart data (#292)

* Write dm (#291)

* Update use of DM in restart

* Remove unused function in Control (#294)

* Change symlink to restart in tests (#295)

* enable tests when old link present

* Extract number empty orbitals from restart file (#296)

* Clean up MD_IonicStepper restart data write (#297)

* Clean up code based on compiler warnings (#299)

* Add getForces for evaluateDMandEnergyAndForces (#300)

* Add check for compatibility MVP/Mehrstellen (#301)

* Add support for Br atom (#302)

* Clean up some HDFrestart functions (#303)

* Fix and test restart single hdf5 file (#305)

* Save Hartree potential for write in restart file (#306)

* Speed-up recently added tests (#309)

* simply use a smaller domain/mesh

* Encapsulate some functions in Potentials (#310)

* some functions were called at random places and confusing

* Remove confusing 0 in naming restart files (#308)

* use added integer only in case of fail/retry

* Add functionalities for extra info in restart file (#312)

* MPI abort (#313)

* Fix error code used by MPI_Abort

* Rho and VH restart  (#311)

* enable restart with consistent rho and VHartree

---------

Co-authored-by: Seung Whan Chung <[email protected]>

* Strenghten testIons (#315)

* More clean up in class Potentials (#316)

* remove some unused functions
* make use of class MGmol_MPI
* change some function names for clarity
* fix some minor bug in iterative indexes

* Misc code fixes (#317)

* exposed by porting on new environment

* Use unique restart filenames in test HDF5single (#318)

* Swap ions in existing test (#319)

* strengthen testing by swapping ions in testRestartEnergyAndForces

* Code clean up (#320)

* use more const
* use initNuc() instead of moveVnuc()
* remove incorrect assert
* use MGmol_MPI in more cases
* setup g_kbpsi_ inside initNuc()
* rename initNuc() into setupPotentials()

* Added functionalities to set local forces (#321)

* Update mixed precision code (#322)

* enable cmake build with mixed precision
* enable use of BLIS
* isolate blas3 loop rewrite in separate files
* fix misc issues with mixed precision code

* Clean up and fixes Ions (#326)

* reset static index counter ion Ion when needed
* simplify some loops
* make Ions copy constructor private

* Fixes for build without HDF5P (#324)

* Fix testRhoVhRestart (#325)

* Fix a few more issues with class Ions (#328)

* Introduce new constructor for class Ions (#330)

* modify various functions to use argument of type Ions instead of MGmol::ions_

---------

Co-authored-by: Jean-Luc Fattebert <[email protected]>
Co-authored-by: Seung Whan Chung <[email protected]>

Author: 3 people

src/ABPG.cc

@@ -38,7 +38,7 @@ void ABPG<T>::setup(T& orbitals)
 //
 // orthof=true: wants orthonormalized updated wave functions
 template <class T>
-int ABPG<T>::updateWF(T& orbitals, Ions& /*ions*/, const double precond_factor,
+int ABPG<T>::updateWF(T& orbitals, Ions& ions, const double precond_factor,
     const bool /*orthof*/, T& work_orbitals, const bool accelerate,
     const bool print_res, const double atol)
 {
@@ -51,8 +51,8 @@ int ABPG<T>::updateWF(T& orbitals, Ions& /*ions*/, const double precond_factor,
     T res("Residual", orbitals, false);
 
     const bool check_res = (atol > 0.);
-    double normRes       = mgmol_strategy_->computeResidual(
-        orbitals, work_orbitals, res, (print_res || check_res), check_res);
+    double normRes = mgmol_strategy_->computeResidual(orbitals, work_orbitals,
+        ions, res, (print_res || check_res), check_res);
     if (normRes < atol && check_res)
     {
         abpg_nl_update_tm_.stop();

src/DFTsolver.cc

@@ -196,8 +196,8 @@ double DFTsolver<OrbitalsType>::evaluateEnergy(
 
     // Get the new total energy
     const double ts = 0.5 * proj_matrices_->computeEntropy(); // in [Ha]
-    eks_history_[0]
-        = energy_->evaluateTotal(ts, proj_matrices_, orbitals, print_flag, os_);
+    eks_history_[0] = energy_->evaluateTotal(
+        ts, proj_matrices_, ions_, orbitals, print_flag, os_);
 
     sum_eig_[1] = sum_eig_[0];
     sum_eig_[0] = 2. * proj_matrices_->getEigSum(); // 2.*sum in [Ry]

src/DavidsonSolver.cc

@@ -500,7 +500,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 ts0 = evalEntropy(projmatrices, (ct.verbose > 1), os_);
                 e0  = energy_->evaluateTotal(
-                    ts0, projmatrices, orbitals, printE, os_);
+                    ts0, projmatrices, ions_, orbitals, printE, os_);
 
                 retval = checkConvergence(e0, outer_it, ct.conv_tol);
                 if (retval == 0 || (outer_it == ct.max_electronic_steps))
@@ -549,7 +549,7 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 ts0 = evalEntropy(proj_mat2N_.get(), (ct.verbose > 1), os_);
                 e0  = energy_->evaluateTotal(
-                    ts0, proj_mat2N_.get(), orbitals, printE, os_);
+                    ts0, proj_mat2N_.get(), ions_, orbitals, printE, os_);
             }
 
             // 2N x 2N target...
@@ -621,8 +621,8 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 const double ts1
                     = evalEntropy(proj_mat2N_.get(), (ct.verbose > 2), os_);
-                const double e1 = energy_->evaluateTotal(
-                    ts1, proj_mat2N_.get(), orbitals, ct.verbose - 1, os_);
+                const double e1 = energy_->evaluateTotal(ts1, proj_mat2N_.get(),
+                    ions_, orbitals, ct.verbose - 1, os_);
 
                 // line minimization
                 beta = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);

src/Energy.cc

@@ -25,11 +25,10 @@
 #define RY2HA 0.5
 
 template <class T>
-Energy<T>::Energy(const pb::Grid& mygrid, const Ions& ions,
-    const Potentials& pot, const Electrostatic& es, const Rho<T>& rho,
-    const XConGrid& xc, SpreadPenaltyInterface<T>* spread_penalty)
+Energy<T>::Energy(const pb::Grid& mygrid, const Potentials& pot,
+    const Electrostatic& es, const Rho<T>& rho, const XConGrid& xc,
+    SpreadPenaltyInterface<T>* spread_penalty)
     : mygrid_(mygrid),
-      ions_(ions),
       pot_(pot),
       es_(es),
       rho_(rho),
@@ -59,7 +58,7 @@ double Energy<T>::getEVrhoRho() const
 }
 
 template <class T>
-double Energy<T>::evaluateEnergyIonsInVext()
+double Energy<T>::evaluateEnergyIonsInVext(Ions& ions)
 {
     double energy = 0.;
 
@@ -69,12 +68,12 @@ double Energy<T>::evaluateEnergyIonsInVext()
     //(*MPIdata::sout)<<"Energy<T>::evaluateEnergyIonsInVext()"<<std::endl;
     double position[3];
     std::vector<double> positions;
-    positions.reserve(3 * ions_.local_ions().size());
+    positions.reserve(3 * ions.local_ions().size());
 
     // loop over ions
     int nions                             = 0;
-    std::vector<Ion*>::const_iterator ion = ions_.local_ions().begin();
-    while (ion != ions_.local_ions().end())
+    std::vector<Ion*>::const_iterator ion = ions.local_ions().begin();
+    while (ion != ions.local_ions().end())
     {
         (*ion)->getPosition(position);
         positions.push_back(position[0]);
@@ -88,9 +87,9 @@ double Energy<T>::evaluateEnergyIonsInVext()
     pot_.getValVext(positions, val);
 
     // loop over ions again
-    ion           = ions_.local_ions().begin();
+    ion           = ions.local_ions().begin();
     int ion_index = 0;
-    while (ion != ions_.local_ions().end())
+    while (ion != ions.local_ions().end())
     {
         const double z = (*ion)->getZion();
         // int ion_index=(*ion)->index();
@@ -112,16 +111,16 @@ double Energy<T>::evaluateEnergyIonsInVext()
 
 template <class T>
 double Energy<T>::evaluateTotal(const double ts, // in [Ha]
-    ProjectedMatricesInterface* projmatrices, const T& phi, const int verbosity,
-    std::ostream& os)
+    ProjectedMatricesInterface* projmatrices, Ions& ions, const T& phi,
+    const int verbosity, std::ostream& os)
 {
     eval_te_tm_.start();
 
     Control& ct = *(Control::instance());
 
-    const double eself = ions_.energySelf();
-    const double ediff = ions_.energyDiff(ct.bcPoisson);
-    const double eipot = evaluateEnergyIonsInVext();
+    const double eself = ions.energySelf();
+    const double ediff = ions.energyDiff(ct.bcPoisson);
+    const double eipot = evaluateEnergyIonsInVext(ions);
 
     const double eigsum = 0.5 * projmatrices->getExpectationH();
 

src/Energy.h

@@ -11,6 +11,7 @@
 #define MGMOL_ENERGY_H
 
 #include "Grid.h"
+#include "Ions.h"
 #include "Rho.h"
 #include "SpreadPenaltyInterface.h"
 #include "Timer.h"
@@ -20,7 +21,6 @@
 #include <vector>
 
 class Potentials;
-class Ions;
 class Electrostatic;
 class ProjectedMatricesInterface;
 class XConGrid;
@@ -29,7 +29,6 @@ template <class T>
 class Energy
 {
     const pb::Grid& mygrid_;
-    const Ions& ions_;
     const Potentials& pot_;
     const Electrostatic& es_;
     const Rho<T>& rho_;
@@ -45,16 +44,15 @@ class Energy
     double getEVrhoRho() const;
 
 public:
-    Energy(const pb::Grid&, const Ions&, const Potentials&,
-        const Electrostatic&, const Rho<T>&, const XConGrid&,
-        SpreadPenaltyInterface<T>*);
+    Energy(const pb::Grid&, const Potentials&, const Electrostatic&,
+        const Rho<T>&, const XConGrid&, SpreadPenaltyInterface<T>*);
 
     static Timer eval_te_tm() { return eval_te_tm_; }
 
-    double evaluateTotal(const double ts, ProjectedMatricesInterface*,
+    double evaluateTotal(const double ts, ProjectedMatricesInterface*, Ions&,
         const T& phi, const int, std::ostream&);
 
-    double evaluateEnergyIonsInVext();
+    double evaluateEnergyIonsInVext(Ions&);
 
     void saveVofRho();
 };

src/GrassmanLineMinimization.cc

@@ -29,7 +29,7 @@ Timer GrassmanLineMinimization<T>::update_states_tm_("Grassman_update_states");
 //
 // orthof=true: wants orthonormalized updated wave functions
 template <class T>
-int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& /*ions*/,
+int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& ions,
     const double precond_factor, const bool orthof, T& work_orbitals,
     const bool accelerate, const bool print_res, const double atol)
 {
@@ -61,7 +61,7 @@ int GrassmanLineMinimization<T>::updateWF(T& orbitals, Ions& /*ions*/,
     // Update wavefunctions
     const bool check_res = (atol > 0.);
     double normRes = mgmol_strategy_->computeResidual(orbitals, work_orbitals,
-        *new_grad_, (print_res || check_res), check_res);
+        ions, *new_grad_, (print_res || check_res), check_res);
     if (normRes < atol && check_res)
     {
         nl_update_tm_.stop();

src/HamiltonianMVPSolver.cc

@@ -160,8 +160,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
         // compute energy at origin
         const int printE = (ct.verbose > 1) ? 1 : 0;
-        double e0
-            = energy_->evaluateTotal(ts0, projmatrices, orbitals, printE, os_);
+        double e0        = energy_->evaluateTotal(
+            ts0, projmatrices, ions_, orbitals, printE, os_);
 
         //
         // compute energy at end for new H
@@ -189,8 +189,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
         projmatrices->setHB2H();
 
         // compute energy at end (beta=1.)
-        double e1
-            = energy_->evaluateTotal(ts1, projmatrices, orbitals, printE, os_);
+        double e1 = energy_->evaluateTotal(
+            ts1, projmatrices, ions_, orbitals, printE, os_);
 
         //
         // evaluate energy at mid-point
@@ -226,8 +226,8 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
         projmatrices->setHB2H();
 
         // compute energy at midpoint
-        double ei
-            = energy_->evaluateTotal(tsi, projmatrices, orbitals, printE, os_);
+        double ei = energy_->evaluateTotal(
+            tsi, projmatrices, ions_, orbitals, printE, os_);
 
         // line minimization
         double beta
@@ -285,7 +285,7 @@ int HamiltonianMVPSolver<MatrixType, ProjMatrixType, OrbitalsType>::solve(
 
                 // compute energy at end (beta=1.)
                 ei = energy_->evaluateTotal(
-                    tsi, projmatrices, orbitals, printE, os_);
+                    tsi, projmatrices, ions_, orbitals, printE, os_);
 
                 // line minimization
                 beta = minQuadPolynomialFrom3values(