Merge release (#354)

* 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 <seung.chung@austin.utexas.edu>

* 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_

* Add cleanup and MGmolInterface (#331)

* Add function to set local forces (#333)

* based on matching coordinates to local ions

* Fix setLocalForces (#335)

* Fix test WFEnergyAndForces (#336)

* writing restart twice in same directory was failing

* Add mixing option in MVP (#337)

* Update MVP with mixing (#339)

* use tolerance to terminate iterations
* tune verbosity

* Update DensityMatrix class (#342)

* add extrapolation function
* remove unused function
* simplify use of iterative index

* Fixing DM tolerance (#341)

* Fix Density Matrix tolerance read

---------

Co-authored-by: Jean-Luc Fattebert <fattebertj@ornl.gov>

* Rm oversubscribe from CMakeLists.txt (#346)

* already set in build script

* Adapt to new environment on condo @ ORNL (#348)

* Remove unused function

* Clean up use of coords file when restarting

* use coordinates file to override atomic coordinates when using
  restart wavefunctions

* Add test restart with override atomic coordinates

* Add some missing include files

* required for gcc12

* Print MVP timers

---------

Co-authored-by: Jean-Luc Fattebert <fattebertj@ornl.gov>
Co-authored-by: Seung Whan Chung <seung.chung@austin.utexas.edu>

Author: 3 people

src/Control.cc

@@ -112,7 +112,6 @@ Control::Control()
     out_restart_file_type             = -1;
     spread_radius                     = -1.;
     iprint_residual                   = -1;
-    override_restart                  = -1;
     dot_product_type                  = -1;
     spread_penalty_damping_           = -1;
     spread_penalty_alpha_             = -1;
@@ -386,7 +385,7 @@ void Control::sync(void)
         short_buffer[47] = out_restart_file_naming_strategy;
         short_buffer[48] = enforceVmass0;
         short_buffer[49] = dm_inner_steps;
-        short_buffer[50] = override_restart;
+        short_buffer[50] = -1;
         short_buffer[51] = fgmres_kim;
         short_buffer[52] = fgmres_maxits;
         short_buffer[53] = ilu_type;
@@ -599,41 +598,41 @@ void Control::sync(void)
     out_restart_file_naming_strategy = short_buffer[47];
     enforceVmass0                    = short_buffer[48];
     dm_inner_steps                   = short_buffer[49];
-    override_restart                 = short_buffer[50];
-    fgmres_kim                       = short_buffer[51];
-    fgmres_maxits                    = short_buffer[52];
-    ilu_type                         = short_buffer[53];
-    ilu_lof                          = short_buffer[54];
-    ilu_maxfil                       = short_buffer[55];
-    coloring_algo_                   = short_buffer[56];
-    diel_flag_                       = short_buffer[57];
-    poisson_pc_nu1                   = short_buffer[58];
-    poisson_pc_nu2                   = short_buffer[59];
-    poisson_pc_nlev                  = short_buffer[60];
-    system_charge_                   = short_buffer[61];
-    md_print_freq                    = short_buffer[62];
-    use_kernel_functions             = short_buffer[63];
-    ngpts_[0]                        = short_buffer[64];
-    ngpts_[1]                        = short_buffer[65];
-    ngpts_[2]                        = short_buffer[66];
-    computeCondGram_                 = short_buffer[67];
-    lrs_extrapolation                = short_buffer[68];
-    parallel_transport               = (bool)short_buffer[69];
-    with_spin_                       = (bool)short_buffer[70];
-    conv_criterion_                  = short_buffer[71];
-    load_balancing_max_iterations    = short_buffer[72];
-    load_balancing_modulo            = short_buffer[73];
-    write_clusters                   = short_buffer[74];
-    DM_solver_                       = short_buffer[75];
-    dm_algo_                         = short_buffer[80];
-    dm_approx_order                  = short_buffer[81];
-    dm_approx_ndigits                = short_buffer[82];
-    dm_approx_power_maxits           = short_buffer[83];
-    spread_penalty_type_             = short_buffer[84];
-    dm_use_old_                      = short_buffer[85];
-    max_electronic_steps_tight_      = short_buffer[86];
-    hartree_reset_                   = short_buffer[88];
-    MD_last_step_                    = short_buffer[89];
+    //...                 = short_buffer[50];
+    fgmres_kim                    = short_buffer[51];
+    fgmres_maxits                 = short_buffer[52];
+    ilu_type                      = short_buffer[53];
+    ilu_lof                       = short_buffer[54];
+    ilu_maxfil                    = short_buffer[55];
+    coloring_algo_                = short_buffer[56];
+    diel_flag_                    = short_buffer[57];
+    poisson_pc_nu1                = short_buffer[58];
+    poisson_pc_nu2                = short_buffer[59];
+    poisson_pc_nlev               = short_buffer[60];
+    system_charge_                = short_buffer[61];
+    md_print_freq                 = short_buffer[62];
+    use_kernel_functions          = short_buffer[63];
+    ngpts_[0]                     = short_buffer[64];
+    ngpts_[1]                     = short_buffer[65];
+    ngpts_[2]                     = short_buffer[66];
+    computeCondGram_              = short_buffer[67];
+    lrs_extrapolation             = short_buffer[68];
+    parallel_transport            = (bool)short_buffer[69];
+    with_spin_                    = (bool)short_buffer[70];
+    conv_criterion_               = short_buffer[71];
+    load_balancing_max_iterations = short_buffer[72];
+    load_balancing_modulo         = short_buffer[73];
+    write_clusters                = short_buffer[74];
+    DM_solver_                    = short_buffer[75];
+    dm_algo_                      = short_buffer[80];
+    dm_approx_order               = short_buffer[81];
+    dm_approx_ndigits             = short_buffer[82];
+    dm_approx_power_maxits        = short_buffer[83];
+    spread_penalty_type_          = short_buffer[84];
+    dm_use_old_                   = short_buffer[85];
+    max_electronic_steps_tight_   = short_buffer[86];
+    hartree_reset_                = short_buffer[88];
+    MD_last_step_                 = short_buffer[89];
     poisson_lap_type_ = static_cast<PoissonFDtype>(short_buffer[90]);
 
     numst    = int_buffer[0];
@@ -1787,7 +1786,6 @@ void Control::setOptions(const boost::program_options::variables_map& vm)
         // options not available in configure file
         lr_updates_type         = 0;
         precond_factor_computed = false;
-        override_restart        = 0;
         mix_pot                 = 1.;
         project_out_psd         = 0;
         multipole_order         = 1;

src/Control.h

@@ -562,7 +562,6 @@ class Control
     short out_restart_file_naming_strategy;
     short restart_file_type;
     short out_restart_file_type;
-    short override_restart;
 
     short verbose;
 

src/DensityMatrix.cc

@@ -22,12 +22,6 @@
 
 const double factor_kernel4dot = 10.;
 
-#define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
-    {                                                                          \
-        std::cerr << "DensityMatrix failure:" << std::endl;                    \
-        std::cerr << "Error Message: " << X << std::endl;                      \
-    }
-
 #define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
     {                                                                          \
         std::cerr << "DensityMatrix failure:" << std::endl;                    \

src/MGmol.cc

@@ -45,6 +45,7 @@
 #include "MGmol.h"
 #include "MLWFTransform.h"
 #include "MPIdata.h"
+#include "MVPSolver.h"
 #include "MasksSet.h"
 #include "Mesh.h"
 #include "OrbitalsPreconditioning.h"
@@ -901,6 +902,9 @@ void MGmol<OrbitalsType>::printTimers()
     AndersonMix<OrbitalsType>::update_tm().print(os_);
     proj_matrices_->printTimers(os_);
     ShortSightedInverse::printTimers(os_);
+    if (std::is_same<OrbitalsType, ExtendedGridOrbitals>::value)
+        MVPSolver<ExtendedGridOrbitals,
+            dist_matrix::DistMatrix<DISTMATDTYPE>>::printTimers(os_);
     VariableSizeMatrixInterface::printTimers(os_);
     DataDistribution::printTimers(os_);
     PackedCommunicationBuffer::printTimers(os_);

src/MGmol.h

@@ -229,7 +229,6 @@ class MGmol : public MGmolInterface
 
     void printEigAndOcc();
 
-    int readCoordinates(std::ifstream* tfile, const bool cell_relative);
     int readCoordinates(const std::string& filename, const bool cell_relative);
     double computeConstraintResidual(OrbitalsType& orbitals,
         const OrbitalsType& hphi, OrbitalsType& res, const bool print_residual,

src/MVPSolver.cc

@@ -398,12 +398,8 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
 template <class OrbitalsType, class MatrixType>
 void MVPSolver<OrbitalsType, MatrixType>::printTimers(std::ostream& os)
 {
-    if (onpe0)
-    {
-        os << std::setprecision(2) << std::fixed << std::endl;
-        solve_tm_.print(os);
-        target_tm_.print(os);
-    }
+    solve_tm_.print(os);
+    target_tm_.print(os);
 }
 
 template class MVPSolver<LocGridOrbitals,

src/MVPSolver.h

@@ -16,7 +16,6 @@
 
 class Ions;
 class Electrostatic;
-class ProjectedMatrices2N;
 template <class OrbitalsType>
 class ProjectedMatrices;
 
@@ -67,7 +66,7 @@ class MVPSolver
     ~MVPSolver();
 
     int solve(OrbitalsType& orbitals);
-    void printTimers(std::ostream& os);
+    static void printTimers(std::ostream& os);
 };
 
 #endif

src/main.cc

@@ -58,7 +58,7 @@ int main(int argc, char** argv)
     /*
      * read runtime parameters
      */
-    std::string input_filename("");
+    std::string coords_filename("");
     std::string lrs_filename;
     std::string constraints_filename("");
 
@@ -70,7 +70,7 @@ int main(int argc, char** argv)
     // read from PE0 only
     if (MPIdata::onpe0)
     {
-        read_config(argc, argv, vm, input_filename, lrs_filename,
+        read_config(argc, argv, vm, coords_filename, lrs_filename,
             constraints_filename, total_spin, with_spin);
     }
 
@@ -89,18 +89,18 @@ int main(int argc, char** argv)
     int ret = ct.checkOptions();
     if (ret < 0) return ret;
 
-    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(coords_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);
+                coords_filename, lrs_filename, constraints_filename);
         else
             mgmol = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout,
-                input_filename, lrs_filename, constraints_filename);
+                coords_filename, lrs_filename, constraints_filename);
 
         mgmol->setup();
 

src/md.cc

@@ -340,7 +340,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 
     constraints_->printConstraints(os_);
 
-    if (ct.restart_info > 0 && !ct.override_restart)
+    if (ct.restart_info > 0)
     {
         if (onpe0) os_ << "Use restart file to initialize MD..." << std::endl;
         stepper->init(*h5f_file_);

src/readInput.cc

@@ -164,7 +164,7 @@ int MGmol<OrbitalsType>::readLRsFromInput(std::ifstream* tfile)
 
 template <class OrbitalsType>
 int MGmol<OrbitalsType>::readCoordinates(
-    std::ifstream* tfile, const bool cell_relative)
+    const std::string& coords_filename, const bool cell_relative)
 {
     Control& ct = *(Control::instance());
     if (ct.verbose > 0) printWithTimeStamp("Read atomic coordinates...", os_);
@@ -177,71 +177,26 @@ int MGmol<OrbitalsType>::readCoordinates(
     const std::vector<Species>& sp(ct.getSpecies());
     ions_.reset(new Ions(lattice, sp));
 
-    if (ct.restart_info > 0
-        && ct.override_restart == 0) // read restart ionic positions
+    if (ct.restart_info > 0 && coords_filename.empty())
     {
-        if (ct.restart_info > 0)
-        {
-            if (onpe0 && ct.verbose > 0)
-            {
-                os_ << "Initialize ionic positions from restart file "
-                    << ct.restart_file << std::endl;
-            }
-            ions_->initFromRestartFile(*h5f_file_);
-        }
-    }
-    else
-    {
-        // Coordinates and species type for each ion.
-        int info = ions_->readAtoms(tfile, cell_relative);
-
-        return info;
-    }
-
-    const int num_ions = ions_->getNumIons();
-    if (onpe0) os_ << num_ions << " ions in simulation" << std::endl;
-
-    return 0;
-}
-
-template <class OrbitalsType>
-int MGmol<OrbitalsType>::readCoordinates(
-    const std::string& filename, const bool cell_relative)
-{
-    Control& ct = *(Control::instance());
-    if (ct.verbose > 0) printWithTimeStamp("Read atomic coordinates...", os_);
-    Mesh* mymesh           = Mesh::instance();
-    const pb::Grid& mygrid = mymesh->grid();
-
-    const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) };
-
-    // setup ions
-    const std::vector<Species>& sp(ct.getSpecies());
-    ions_.reset(new Ions(lattice, sp));
-
-    if (ct.restart_info > 0
-        && ct.override_restart == 0) // read restart ionic positions
-    {
-        if (ct.restart_info > 0)
+        // read restart atomic positions
+        if (onpe0 && ct.verbose > 0)
         {
-            if (onpe0 && ct.verbose > 0)
-            {
-                os_ << "Initialize ionic positions from restart file "
-                    << ct.restart_file << std::endl;
-            }
-            ions_->initFromRestartFile(*h5f_file_);
+            os_ << "Initialize atomic positions from restart file "
+                << ct.restart_file << std::endl;
         }
+        ions_->initFromRestartFile(*h5f_file_);
     }
     else
     {
         // Coordinates and species type for each ion.
-        int info = ions_->readAtoms(filename, cell_relative);
+        int info = ions_->readAtoms(coords_filename, cell_relative);
 
         return info;
     }
 
     const int num_ions = ions_->getNumIons();
-    if (onpe0) os_ << num_ions << " ions in simulation" << std::endl;
+    if (onpe0) os_ << num_ions << " atoms in simulation" << std::endl;
 
     return 0;
 }