Fix and test restart single hdf5 file

src/Control.h

@@ -371,6 +371,8 @@ class Control
     // 10 or larger means CG, otherwise MG V-cycles
     bool MGPoissonSolver() { return (diel_flag_ / 10 == 0); }
 
+    bool LangevinThermostat() { return (thermostat_type == 1); }
+
     //
     // data
     //

src/HDFrestart.cc

@@ -457,7 +457,8 @@ HDFrestart::HDFrestart(const std::string& filename, const pb::PEenv& pes,
     verbosity_ = 0;
     closed_    = false;
 
-    //(*MPIdata::sout)<<"HDFrestart::HDFrestart(), filename="<<filename<<endl;
+    //(*MPIdata::sout)<<"HDFrestart::HDFrestart(),
+    // filename="<<filename<<std::endl;
     setActivity();
 
     count_[0] = count_[1] = count_[2] = 1;
@@ -749,7 +750,7 @@ int writeListCentersAndRadii(
     }
 
     // if( onpe0 && ct.verbose>2 )
-    //    (*MPIdata::sout)<<"Write attribute "<<attname<<endl;
+    //    (*MPIdata::sout)<<"Write attribute "<<attname<<std::endl;
     herr_t status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &attr_data[0]);
     if (status < 0)
     {
@@ -1921,14 +1922,19 @@ int HDFrestart::readAtomicData(std::string datasetname, std::vector<int>& data)
     // send data to inactive PEs
     if (gather_data_x_) gatherDataXdir(data);
 
+    if (useHdf5p())
+    {
+        data.erase(std::remove(data.begin(), data.end(), -1), data.end());
+    }
+
     return 0;
 }
 
 int HDFrestart::readAtomicData(
     std::string datasetname, std::vector<double>& data)
 {
     if (onpe0)
-        (*MPIdata::sout) << "Read ionic positions from hdf5 file" << std::endl;
+        (*MPIdata::sout) << "Read atomic data from hdf5 file" << std::endl;
 
     if (active_)
     {
@@ -1968,6 +1974,10 @@ int HDFrestart::readAtomicData(
         }
     }
 
+    if (useHdf5p())
+    {
+        data.erase(std::remove(data.begin(), data.end(), 1e+32), data.end());
+    }
     if (gather_data_x_) gatherDataXdir(data);
 
     return 0;
@@ -2094,6 +2104,11 @@ int HDFrestart::readAtomicData(
         data.push_back(t);
     }
 
+    if (useHdf5p())
+    {
+        data.erase(std::remove(data.begin(), data.end(), ""), data.end());
+    }
+
     return 0;
 }
 
@@ -2128,7 +2143,6 @@ int HDFrestart::readRestartRandomStates(std::vector<unsigned short>& data)
                 dim = (int)H5Dget_storage_size(dataset_id)
                       / sizeof(unsigned short);
             }
-
             if (dim > 0)
             {
                 data.resize(dim);
@@ -2148,12 +2162,10 @@ int HDFrestart::readRestartRandomStates(std::vector<unsigned short>& data)
                 }
 
                 if (!data.empty())
-                    if (data[0] != data[0])
+                    if (std::isnan(data[0]))
                     {
                         MGMOL_HDFRESTART_FAIL(
-                            "ERROR: HDFrestart::readRestartRandomStates() "
-                            "--- data[0]="
-                            << data[0]);
+                            "readRestartRandomStates() is NaN");
                         return -2;
                     }
             }

src/Ions.cc

@@ -911,8 +911,10 @@ void Ions::initFromRestartFile(HDFrestart& h5_file)
     assert(at_numbers.size() == at_nlprojIds.size());
 
     num_ions_ = at_names.size();
-    mmpi.allreduce(&num_ions_, 1, MPI_SUM);
-
+    if (!h5_file.useHdf5p())
+    {
+        mmpi.allreduce(&num_ions_, 1, MPI_SUM);
+    }
     if (onpe0 && ct.verbose > 0)
     {
         (*MPIdata::sout) << "Ions::setFromRestartFile(), read " << num_ions_
@@ -947,9 +949,21 @@ void Ions::initFromRestartFile(HDFrestart& h5_file)
     }
     readRestartPositions(h5_file);
     readRestartVelocities(h5_file);
-    readRestartRandomStates(h5_file);
+    if (ct.LangevinThermostat()) readRestartRandomStates(h5_file);
     readLockedAtomNames(h5_file);
 
+    // remove atoms from local list if not local
+    for (std::vector<Ion*>::iterator it = local_ions_.begin();
+         it != local_ions_.end();)
+    {
+        double p[3];
+        (*it)->getPosition(p);
+        if (!inLocalIons(p[0], p[1], p[2]))
+            it = local_ions_.erase(it);
+        else
+            ++it;
+    }
+
     // rescale all velocities by factor specified in input
     rescaleVelocities(ct.VelocityScalingFactor());
 

src/restart.cc

@@ -117,7 +117,7 @@ int MGmol<OrbitalsType>::write_hdf5(HDFrestart& h5f_file,
         ions.writeAtomicIDs(h5f_file);
         ions.writeAtomicNLprojIDs(h5f_file);
         ions.writePositions(h5f_file);
-        ions.writeRandomStates(h5f_file);
+        if (ct.LangevinThermostat()) ions.writeRandomStates(h5f_file);
         ions.writeVelocities(h5f_file);
         ions.writeForces(h5f_file);
 

tests/CMakeLists.txt

@@ -504,6 +504,14 @@ add_test(NAME testMD_D72
          ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/lrs.in
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+add_test(NAME testHDF5single
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
+         ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/mgmol.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/md.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/HDF5single/h2o.xyz
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 add_test(NAME testMD_MVP
          COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/test.py
          ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}

tests/DMandEnergyAndForces/test.py

@@ -3,6 +3,7 @@
 import os
 import subprocess
 import string
+import shutil
 
 print("Test DMandEnergyAndForces...")
 
@@ -33,6 +34,8 @@
 output = subprocess.check_output(command,shell=True)
 lines=output.split(b'\n')
 
+shutil.rmtree('WF')
+
 #analyse output
 energies=[]
 for line in lines:

tests/HDF5single/h2o.xyz

@@ -0,0 +1,6 @@
+3
+
+O  0.00   0.00  0.00
+H  -0.76  0.59  0.00
+H  0.76   0.59 0.00
+

tests/HDF5single/md.cfg

@@ -0,0 +1,34 @@
+verbosity=3
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx=48
+ny=48
+nz=48
+[Domain]
+ox=-4.5
+oy=-4.5
+oz=-4.5
+lx=9.
+ly=9.
+lz=9.
+[Potentials]
+pseudopotential=pseudo.O_ONCV_PBE_SG15
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+[Run]
+type=MD
+[MD]
+num_steps=5
+dt=40.
+[Quench]
+max_steps=24
+atol=1.e-8
+[Restart]
+input_level=4
+input_filename=WF
+input_type=single_file
+output_level=4
+output_filename=WF_MD
+output_type=single_file
+[Coloring]
+scope=global

tests/HDF5single/mgmol.cfg

@@ -0,0 +1,31 @@
+verbosity=2
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx=48
+ny=48
+nz=48
+[Domain]
+ox=-4.5
+oy=-4.5
+oz=-4.5
+lx=9.
+ly=9.
+lz=9.
+[Potentials]
+pseudopotential=pseudo.O_ONCV_PBE_SG15
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+[Run]
+type=QUENCH
+[Quench]
+max_steps=120
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=1.5
+[Restart]
+output_level=4
+output_filename=WF
+output_type=single_file
+[Coloring]
+scope=global

tests/HDF5single/test.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test test_rho_restart...")
+
+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-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 to generate initial ground state
+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')
+
+#run MD
+command = "{} {} -c {} -i {}".format(mpicmd,mgmol_exe,input2,coords)
+print("Run command: {}".format(command))
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+
+os.remove('WF')
+
+print("Check energy conservation...")
+tol = 1.e-4
+energy = 0.
+count = 0
+for line in lines:
+  if line.count(b'Total') and line.count(b'Energy'):
+    print(line)
+    count=count+1
+    words=line.split()
+
+    energy=eval(words[2])
+    if count==1:
+      first_energy=energy
+
+    if count>1 and abs(energy-first_energy)>tol:
+      print("ERROR Energy = {} != {}".format(energy,first_energy))
+      sys.exit(1)
+
+if count<4:
+  print("ERROR needs 4 energy values for checking conservation!")
+  sys.exit(1)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)

tests/RestartEnergyAndForces/test.py

@@ -3,6 +3,7 @@
 import os
 import subprocess
 import string
+import shutil
 
 print("Test RestartEnergyAndForces...")
 
@@ -54,14 +55,14 @@
     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')
 
+shutil.rmtree('WF')
+
 test_energy=1.e18
 for line in lines:
   if line.count(b'%%'):

tests/WFEnergyAndForces/test.py

@@ -3,6 +3,7 @@
 import os
 import subprocess
 import string
+import shutil
 
 print("Test WFEnergyAndForces...")
 
@@ -37,6 +38,8 @@
 output = subprocess.check_output(command,shell=True)
 lines=output.split(b'\n')
 
+shutil.rmtree('WF')
+
 #analyse output
 energies=[]
 for line in lines: