Extract number empty orbitals from restart file

src/Control.h

@@ -242,6 +242,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); }
 
@@ -263,6 +264,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); }

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/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

tests/RestartEnergyAndForces/test.py

@@ -0,0 +1,88 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test RestartEnergyAndForces...")
+
+nargs=len(sys.argv)
+
+mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3]
+for i in range(4,nargs-7):
+  mpicmd = mpicmd + " "+sys.argv[i]
+print("MPI run command: {}".format(mpicmd))
+
+mgmol_exe = sys.argv[nargs-6]
+test_exe = sys.argv[nargs-5]
+input1 = sys.argv[nargs-4]
+input2 = sys.argv[nargs-3]
+coords = sys.argv[nargs-2]
+print("coordinates file: %s"%coords)
+
+#create links to potentials files
+dst1 = 'pseudo.H_ONCV_PBE_SG15'
+src1 = sys.argv[-1] + '/' + dst1
+
+dst2 = 'pseudo.O_ONCV_PBE_SG15'
+src2 = sys.argv[-1] + '/' + dst2
+
+if not os.path.exists(dst1):
+  print("Create link to %s"%dst1)
+  os.symlink(src1, dst1)
+
+if not os.path.exists(dst2):
+  print("Create link to %s"%dst2)
+  os.symlink(src2, dst2)
+
+#run mgmol
+command = "{} {} -c {} -i {}".format(mpicmd,mgmol_exe,input1,coords)
+print("Run command: {}".format(command))
+
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+
+#analyse output
+ref_energy=1.e18
+for line in lines:
+  if line.count(b'%%'):
+    print(line)
+    words=line.split()
+    words=words[5].split(b',')[0]
+    energy = words.decode()
+  if line.count(b'achieved'):
+    ref_energy=energy
+    break
+
+#sys.exit(0)
+
+#run test
+command = "{} {} -c {} -i {}".format(mpicmd,test_exe,input2,coords)
+print("Run command: {}".format(command))
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+
+test_energy=1.e18
+for line in lines:
+  if line.count(b'%%'):
+    print(line)
+    words=line.split()
+    words=words[5].split(b',')[0]
+    energy = words.decode()
+  if line.count(b'Eks'):
+    print(line)
+    words=line.split()
+    print(words)
+    test_energy = words[2]
+    break
+
+
+tol = 1.e-6
+diff=eval(test_energy)-eval(ref_energy)
+print(diff)
+if abs(diff)>tol:
+  print("Energies differ: {} vs {} !!!".format(ref_energy,test_energy))
+  sys.exit(1)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)