Rho and VH restart (#311)

* enable restart with consistent rho and VHartree --------- Co-authored-by: Seung Whan Chung <[email protected]>
LLNL · Feb 24, 2025 · 8b84de4 · 8b84de4
1 parent 77bd87e
commit 8b84de4
Show file tree

Hide file tree

Showing 12 changed files with 430 additions and 1 deletion.
diff --git a/src/Electrostatic.h b/src/Electrostatic.h
@@ -47,6 +47,9 @@ class Electrostatic
     ~Electrostatic();
     static Timer solve_tm() { return solve_tm_; }
 
+    pb::GridFunc<RHODTYPE>* getRhoc() { return grhoc_; }
+    Poisson* getPoissonSolver() { return poisson_solver_; }
+
     void setup(const short max_sweeps);
     void setupPB(const double e0, const double rho0, const double drho0,
         Potentials& pot);

diff --git a/src/MGmol.h b/src/MGmol.h
@@ -150,7 +150,6 @@ class MGmol : public MGmolInterface
     void initialMasks();
     int setupLRsFromInput(const std::string filename);
 
-    void setup();
     int setupLRs(const std::string input_file) override;
     int setupFromInput(const std::string input_file) override;
     int setupConstraintsFromInput(const std::string input_file) override;
@@ -177,12 +176,15 @@ class MGmol : public MGmolInterface
 
     ~MGmol() override;
 
+    void setup();
+
     /* access functions */
     OrbitalsType* getOrbitals() { return current_orbitals_; }
     std::shared_ptr<Hamiltonian<OrbitalsType>> getHamiltonian()
     {
         return hamiltonian_;
     }
+    std::shared_ptr<Rho<OrbitalsType>> getRho() { return rho_; }
 
     void run() override;
 

diff --git a/src/Potentials.cc b/src/Potentials.cc
@@ -939,6 +939,12 @@ int Potentials::read(HDFrestart& h5f_file)
         h5f_file.read_1func_hdf5(vepsilon_.data(), "VDielectric");
     }
 
+    std::string datasetname("Preceding_Hartree");
+    if (h5f_file.checkDataExists(datasetname))
+    {
+        h5f_file.read_1func_hdf5(vh_rho_backup_.data(), datasetname);
+    }
+
     return 0;
 }
 

diff --git a/src/Potentials.h b/src/Potentials.h
@@ -206,6 +206,8 @@ class Potentials
      */
     void backupVh();
 
+    void resetVhRho2Backup() { vh_rho_ = vh_rho_backup_; }
+
 #ifdef HAVE_TRICUBIC
     void readExternalPot(const string filename, const char type);
     void setupVextTricubic();

diff --git a/src/md.cc b/src/md.cc
@@ -702,6 +702,19 @@ void MGmol<OrbitalsType>::loadRestartFile(const std::string filename)
             ierr = proj_matrices_->readWFDM(h5file);
     }
 
+    if (h5file.checkDataExists("Preceding_Hartree"))
+    {
+        ions_->readRestartPreviousPositions(h5file);
+        ions_->resetPositionsToPrevious();
+        ions_->setup();
+
+        Potentials& pot = hamiltonian_->potential();
+        pot.initialize(*ions_);
+        if (onpe0) std::cout << "Reset VhRho to backup..." << std::endl;
+        pot.resetVhRho2Backup();
+        electrostat_->setupRhoc(pot.rho_comp());
+    }
+
     ierr = h5file.close();
     mmpi.allreduce(&ierr, 1, MPI_MIN);
     if (ierr < 0)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
@@ -250,6 +250,8 @@ add_executable(testDensityMatrix
                ${CMAKE_SOURCE_DIR}/src/ReplicatedWorkSpace.cc
                ${CMAKE_SOURCE_DIR}/src/hdf_tools.cc
                ${CMAKE_SOURCE_DIR}/tests/ut_magma_main.cc)
+add_executable(testRhoVhRestart
+               ${CMAKE_SOURCE_DIR}/tests/RhoVhRestart/testRhoVhRestart.cc)
 add_executable(testEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc)
 add_executable(testWFEnergyAndForces
@@ -382,6 +384,17 @@ add_test(NAME testRestartEnergyAndForces
          ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/restart.cfg
          ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/h2o.xyz
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+add_test(NAME testRhoVhRestart
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/RhoVhRestart/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
+         ${CMAKE_CURRENT_BINARY_DIR}/testRhoVhRestart
+         ${CMAKE_CURRENT_SOURCE_DIR}/RhoVhRestart/mgmol.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/RhoVhRestart/md.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/RhoVhRestart/restart.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/RhoVhRestart/h2o.xyz
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+
 
 if(${MAGMA_FOUND}) 
   add_test(NAME testOpenmpOffload
@@ -577,6 +590,7 @@ target_include_directories(testDensityMatrix PRIVATE ${Boost_INCLUDE_DIRS} ${HDF
 target_include_directories(testGramMatrix PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testAndersonMix PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testIons PRIVATE ${Boost_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS})
+target_include_directories(testRhoVhRestart PRIVATE ${Boost_INCLUDE_DIRS})
 
 target_link_libraries(testMPI PRIVATE MPI::MPI_CXX)
 target_link_libraries(testBlacsContext PRIVATE ${SCALAPACK_LIBRARIES}
@@ -589,6 +603,7 @@ 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})
+target_link_libraries(testRhoVhRestart mgmol_src)
 
 if(${MAGMA_FOUND})
   target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES}

diff --git a/tests/RhoVhRestart/h2o.xyz b/tests/RhoVhRestart/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
+
diff --git a/tests/RhoVhRestart/md.cfg b/tests/RhoVhRestart/md.cfg
@@ -0,0 +1,30 @@
+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=MD
+[MD]
+num_steps=5
+dt=40.
+[Quench]
+max_steps=24
+atol=1.e-8
+[Restart]
+input_level=4
+input_filename=WF
+output_level=4
+output_filename=WF_MD
diff --git a/tests/RhoVhRestart/mgmol.cfg b/tests/RhoVhRestart/mgmol.cfg
@@ -0,0 +1,28 @@
+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
diff --git a/tests/RhoVhRestart/restart.cfg b/tests/RhoVhRestart/restart.cfg
@@ -0,0 +1,25 @@
+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=24
+atol=1.e-8
+[Restart]
+input_level=4
+input_filename=WF_MD
diff --git a/tests/RhoVhRestart/test.py b/tests/RhoVhRestart/test.py
@@ -0,0 +1,79 @@
+#!/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-7]
+test_exe = sys.argv[nargs-6]
+input1 = sys.argv[nargs-5]
+input2 = sys.argv[nargs-4]
+input3 = 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')
+
+flag=0
+for line in lines:
+  if line.count(b'Run ended'):
+    flag=1
+
+if flag==0:
+  print("Initial quench failed to complete!")
+  sys.exit(1)
+
+#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')
+
+flag=0
+for line in lines:
+  if line.count(b'Run ended'):
+    flag=1
+
+if flag==0:
+  print("MD failed to complete!")
+  sys.exit(1)
+
+#run test
+command = "{} {} -c {} -i {}".format(mpicmd,test_exe,input3,coords)
+print("Run command: {}".format(command))
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+for line in lines:
+  print(line)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)