Merge release to ROMFPMD (#298)

src/Control.cc

@@ -1119,78 +1119,6 @@ void Control::printThermostatInfo(std::ostream& os) const
     }
 }
 
-int Control::readOccupations(std::ifstream* tfile)
-{
-    int count = 0;
-    float nel = 0.;
-    do
-    {
-        float t1 = 0.;
-        int nst  = 0;
-        if (mype_ == 0)
-        {
-#ifdef DEBUG
-            (*MPIdata::sout) << " Occupations of states..." << std::endl;
-#endif
-            (*tfile) >> nst;
-            if (nst <= 0)
-            {
-                (*MPIdata::sout)
-                    << "Control::readOccupations: numst=" << numst << std::endl;
-                (*MPIdata::sout) << "Control::readOccupations: nst=" << nst
-                                 << ", count=" << count << std::endl;
-                (*MPIdata::sout) << "Control::readOccupations: Bad repeat "
-                                    "count for state occupations"
-                                 << std::endl;
-                return -1;
-            }
-            if ((count + nst) > numst)
-            {
-                (*MPIdata::sout) << "Control::readOccupations: Occupations "
-                                    "specified for too many states"
-                                 << std::endl;
-                return -1;
-            }
-
-            (*tfile) >> t1;
-            if (t1 < 0.)
-            {
-                (*MPIdata::sout)
-                    << "Control::readOccupations: occupation=" << t1
-                    << std::endl;
-                (*MPIdata::sout) << "Control::readOccupations: occupation "
-                                    "should be a positive number"
-                                 << std::endl;
-                return -1;
-            }
-            finishRead(*tfile);
-        }
-        int mpirc = MPI_Bcast(&nst, 1, MPI_INT, 0, comm_global_);
-        if (mpirc != MPI_SUCCESS)
-        {
-            (*MPIdata::sout)
-                << "MPI Bcast of occupation numbers failed!!!" << std::endl;
-            return -1;
-        }
-        mpirc = MPI_Bcast(&t1, 1, MPI_FLOAT, 0, comm_global_);
-        if (mpirc != MPI_SUCCESS)
-        {
-            (*MPIdata::sout)
-                << "MPI Bcast of occupation failed!!!" << std::endl;
-            return -1;
-        }
-        nel += nst * t1;
-        count += nst;
-
-    } while (count < numst);
-
-    nel_ = (int)nel;
-
-    nempty_ = (2 * numst - (int)nel) / 2;
-
-    return count;
-}
-
 void Control::setLocMode(const float radius, const float lx, const float ly,
     const float lz, const float mind_centers)
 {

src/Control.h

@@ -249,6 +249,7 @@ class Control
         const float total_spin, std::string run_directory = ".");
 
     void setDefaultValues();
+    bool withSpin() { return with_spin_; }
 
     bool globalColoring() const { return (coloring_algo_ / 10 == 0); }
 
@@ -270,6 +271,8 @@ class Control
 
     float getSpin() const { return total_spin_; }
 
+    void setNempty(const int nempty) { nempty_ = nempty; }
+
     short getMGlevels() { return mg_levels_; }
 
     bool withPreconditioner() const { return (mg_levels_ >= 0); }
@@ -301,7 +304,6 @@ class Control
         const float dtol, const short kim, const short itmax, const short lfil,
         const short maxfill, const short ilutype);
     void setSpreadRadius();
-    int readOccupations(std::ifstream* tfile);
     bool checkTimeout() { return timeout_.check(); }
 
     bool occupationWidthIsZero() { return occ_width < 1.e-12; }

src/HDFrestart.cc

@@ -2725,6 +2725,21 @@ void HDFrestart::gatherDataXdir(std::vector<T>& data)
     }
 }
 
+int HDFrestart::countFunctionObjects(std::string& name) const
+{
+    int count = 0;
+    int found = 0;
+    do
+    {
+        std::string datasetname(getDatasetName(name, count));
+        // check if dataset exists...
+        found = checkDataExists(datasetname);
+        if (found) count++;
+    } while (found); // dataset exists
+
+    return count;
+}
+
 template int HDFrestart::read_1func_hdf5(float*, const std::string&);
 template int HDFrestart::read_1func_hdf5(double*, const std::string&);
 

src/HDFrestart.h

@@ -173,6 +173,7 @@ class HDFrestart
         return checkDataExistsLocal(datasetname.c_str());
     }
 
+    // Returns zero (false), a positive (true) or a negative (failure) value.
     herr_t checkDataExistsLocal(const char* const datasetname) const
     {
         herr_t err_id = 0;
@@ -282,6 +283,8 @@ class HDFrestart
     int getMDstepFromFile() const;
     int getFromFile(const std::string& attname) const;
 
+    int countFunctionObjects(std::string& name) const;
+
     hid_t createPlist()
     {
         hid_t plist_id = H5P_DEFAULT;

src/Hamiltonian.cc

@@ -59,9 +59,9 @@ const T& Hamiltonian<T>::applyLocal(T& phi, const bool force)
     if (onpe0)
     {
         (*MPIdata::sout) << "Hamiltonian<T>::applyLocal(), new_index ="
-                         << new_index << endl;
+                         << new_index << std::endl;
         (*MPIdata::sout) << "Hamiltonian<T>::applyLocal(), itindex_  ="
-                         << itindex_ << endl;
+                         << itindex_ << std::endl;
     }
 #endif
     if (force || new_index != itindex_)
@@ -76,7 +76,8 @@ const T& Hamiltonian<T>::applyLocal(T& phi, const bool force)
         if (onpe0)
             (*MPIdata::sout)
                 << "Hamiltonian<T>::hlphi up to date, itindex_=" << itindex_
-                << endl;
+                << ", Potential index=" << pot_->getIterativeIndex()
+                << std::endl;
 #endif
     }
     return *hlphi_;

src/MVPSolver.cc

@@ -143,7 +143,6 @@ void MVPSolver<OrbitalsType, MatrixType>::buildTarget_MVP(
     proj_mat_work_->setHB2H();
 
     proj_mat_work_->updateDM(orbitals_index);
-
     target = proj_mat_work_->dm();
 
     if (ct.verbose > 2)
@@ -286,7 +285,7 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                     std::cout << "MVP --- Target energy..." << std::endl;
                 proj_mat_work_->setDM(target, orbitals.getIterativeIndex());
                 proj_mat_work_->computeOccupationsFromDM();
-                if (ct.verbose > 2) current_proj_mat->printOccupations(os_);
+                if (ct.verbose > 2) proj_mat_work_->printOccupations(os_);
                 const double nel = proj_mat_work_->getNel();
                 if (onpe0 && ct.verbose > 1)
                     os_ << "MVP --- Number of electrons at beta=1 : " << nel
@@ -346,7 +345,6 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                         os_ << "Number of electrons for interpolated DM = "
                             << pnel << std::endl;
                 }
-                // if( onpe0 )os_<<"Rho..."<<endl;
                 rho_->computeRho(orbitals, *work_);
             }
 

src/setup.cc

@@ -50,6 +50,7 @@ int MGmol<OrbitalsType>::setupFromInput(const std::string filename)
     if (ct.isLocMode()) mymesh->subdivGridx(ct.getMGlevels());
 
     const pb::PEenv& myPEenv = mymesh->peenv();
+
     if (ct.restart_info > 0)
         h5f_file_.reset(
             new HDFrestart(ct.restart_file, myPEenv, ct.restart_file_type));
@@ -58,8 +59,20 @@ int MGmol<OrbitalsType>::setupFromInput(const std::string filename)
     if (status == -1) return -1;
 
     const short myspin = mmpi.myspin();
-    const int nval     = ions_->getNValenceElectrons();
-    ct.setNumst(myspin, nval);
+    const int nel      = ions_->getNValenceElectrons();
+    // for the case of extended wavefunctions, we can determine the number
+    // of empty states from the number of wavefunctions in restart file
+    if (ct.restart_info > 2 && !ct.short_sighted)
+    {
+        std::string name = "Function";
+        int count        = h5f_file_->countFunctionObjects(name);
+        std::cout << "found " << count << " functions in restart file..."
+                  << std::endl;
+        int nempty = ct.withSpin() ? count - nel : count - int(0.5 * nel);
+        ct.setNempty(nempty);
+    }
+    ct.setNumst(myspin, nel);
+
     ct.setTolEnergy();
     ct.setSpreadRadius();
 

tests/CMakeLists.txt

@@ -256,6 +256,8 @@ add_executable(testWFEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/WFEnergyAndForces/testWFEnergyAndForces.cc)
 add_executable(testDMandEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc)
+add_executable(testRestartEnergyAndForces
+               ${CMAKE_SOURCE_DIR}/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc)
 
 if(${MAGMA_FOUND})
   add_executable(testOpenmpOffload
@@ -371,6 +373,15 @@ add_test(NAME testDMandEnergyAndForces
          ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/lrs.in
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+add_test(NAME testRestartEnergyAndForces
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
+         ${CMAKE_CURRENT_BINARY_DIR}/testRestartEnergyAndForces
+         ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/mgmol.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/restart.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/h2o.xyz
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 
 if(${MAGMA_FOUND}) 
   add_test(NAME testOpenmpOffload
@@ -567,6 +578,7 @@ target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX)
 target_link_libraries(testEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src)
+target_link_libraries(testRestartEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testIons PRIVATE mgmol_src)
 target_link_libraries(testDensityMatrix PRIVATE ${HDF5_LIBRARIES})
 

tests/FIRE/test.py

@@ -49,7 +49,11 @@
 restart_file = subprocess.check_output(command,shell=True)
 restart_file=str(restart_file[:-1],'utf-8')
 print(restart_file)
-os.symlink(restart_file, 'wave.out')
+try:
+  os.symlink(restart_file, 'wave.out')
+except FileExistsError:
+  os.remove('wave.out')
+  os.symlink(restart_file, 'wave.out')
 
 command = "{} {} -c {}".format(mpicmd,exe,inp2)
 print(command)

tests/LBFGS/test.py

@@ -49,7 +49,12 @@
 restart_file = subprocess.check_output(command,shell=True)
 restart_file=str(restart_file[:-1],'utf-8')
 print(restart_file)
-os.symlink(restart_file, 'wave.out')
+
+try:
+  os.symlink(restart_file, 'wave.out')
+except FileExistsError:
+  os.remove('wave.out')
+  os.symlink(restart_file, 'wave.out')
 
 command = "{} {} -c {} -i {} -l {}".format(mpicmd,exe,inp2,coords,lrs)
 print(command)
@@ -75,6 +80,7 @@
 os.remove('wave.out')
 
 if force>tol:
+  print("Force larger than tol {}".format(tol))
   sys.exit(1)
 
 sys.exit(0)

tests/MD_D72/test.py

@@ -48,7 +48,11 @@
 restart_file=str(restart_file[:-1],'utf-8')
 print(restart_file)
 
-os.symlink(restart_file, 'wave.out')
+try:
+  os.symlink(restart_file, 'wave.out')
+except FileExistsError:
+  os.remove('wave.out')
+  os.symlink(restart_file, 'wave.out')
 
 command = "{} {} -c {} -i {}".format(mpicmd,exe,inp2,coords)
 output2 = subprocess.check_output(command,shell=True)

tests/RestartEnergyAndForces/h2o.xyz

@@ -0,0 +1,6 @@
+3
+https://pubchem.ncbi.nlm.nih.gov/compound/Water
+O  2.5369   -0.1550    0.0
+H  3.0739    0.1550    0.0
+H  2.0000    0.1550    0.0
+

tests/RestartEnergyAndForces/mgmol.cfg

@@ -0,0 +1,33 @@
+verbosity=2
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx=64
+ny=64
+nz=64
+[Domain]
+ox=-3.4
+oy=-6.4
+oz=-6.4
+lx=12.8
+ly=12.8
+lz=12.8
+[Potentials]
+pseudopotential=pseudo.O_ONCV_PBE_SG15
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+[Run]
+type=QUENCH
+[Quench]
+solver=PSD
+max_steps=120
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=1.5
+nempty=2
+[Restart]
+output_level=3
+output_filename=WF
+[DensityMatrix]
+solver=MVP
+nb_inner_it=1

tests/RestartEnergyAndForces/restart.cfg

@@ -0,0 +1,32 @@
+verbosity=2
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx=64
+ny=64
+nz=64
+[Domain]
+ox=-3.4
+oy=-6.4
+oz=-6.4
+lx=12.8
+ly=12.8
+lz=12.8
+[Potentials]
+pseudopotential=pseudo.O_ONCV_PBE_SG15
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+[Run]
+type=QUENCH
+[Quench]
+solver=PSD
+max_steps=120
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=1.5
+[Restart]
+input_level=3
+input_filename=WF
+[DensityMatrix]
+solver=MVP
+nb_inner_it=20