From 8b84de49791659a4d874c617f71a920fb08e9138 Mon Sep 17 00:00:00 2001 From: Jean-Luc Fattebert Date: Mon, 24 Feb 2025 15:03:24 -0500 Subject: [PATCH] Rho and VH restart (#311) * enable restart with consistent rho and VHartree --------- Co-authored-by: Seung Whan Chung --- src/Electrostatic.h | 3 + src/MGmol.h | 4 +- src/Potentials.cc | 6 + src/Potentials.h | 2 + src/md.cc | 13 ++ tests/CMakeLists.txt | 15 ++ tests/RhoVhRestart/h2o.xyz | 6 + tests/RhoVhRestart/md.cfg | 30 ++++ tests/RhoVhRestart/mgmol.cfg | 28 ++++ tests/RhoVhRestart/restart.cfg | 25 +++ tests/RhoVhRestart/test.py | 79 +++++++++ tests/RhoVhRestart/testRhoVhRestart.cc | 220 +++++++++++++++++++++++++ 12 files changed, 430 insertions(+), 1 deletion(-) create mode 100644 tests/RhoVhRestart/h2o.xyz create mode 100644 tests/RhoVhRestart/md.cfg create mode 100644 tests/RhoVhRestart/mgmol.cfg create mode 100644 tests/RhoVhRestart/restart.cfg create mode 100755 tests/RhoVhRestart/test.py create mode 100644 tests/RhoVhRestart/testRhoVhRestart.cc diff --git a/src/Electrostatic.h b/src/Electrostatic.h index 9beee899..06a9e17c 100644 --- a/src/Electrostatic.h +++ b/src/Electrostatic.h @@ -47,6 +47,9 @@ class Electrostatic ~Electrostatic(); static Timer solve_tm() { return solve_tm_; } + pb::GridFunc* 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 index bccbc7af..b18465a1 100644 --- 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> getHamiltonian() { return hamiltonian_; } + std::shared_ptr> getRho() { return rho_; } void run() override; diff --git a/src/Potentials.cc b/src/Potentials.cc index ebe42320..96375c6d 100644 --- 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 index cb0f8e05..8e9b763c 100644 --- 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 index 8ec20074..b1ffb619 100644 --- a/src/md.cc +++ b/src/md.cc @@ -702,6 +702,19 @@ void MGmol::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 index baf1c451..838ce6a4 100644 --- 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 new file mode 100644 index 00000000..d5171c8b --- /dev/null +++ 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 new file mode 100644 index 00000000..2b8a378b --- /dev/null +++ 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 new file mode 100644 index 00000000..eee7f11c --- /dev/null +++ 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 new file mode 100644 index 00000000..20f0293a --- /dev/null +++ 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 new file mode 100755 index 00000000..a34b962f --- /dev/null +++ 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) diff --git a/tests/RhoVhRestart/testRhoVhRestart.cc b/tests/RhoVhRestart/testRhoVhRestart.cc new file mode 100644 index 00000000..09da3320 --- /dev/null +++ b/tests/RhoVhRestart/testRhoVhRestart.cc @@ -0,0 +1,220 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge +// National Laboratory. +// LLNL-CODE-743438 +// All rights reserved. +// This file is part of MGmol. For details, see https://github.com/llnl/mgmol. +// Please also read this link https://github.com/llnl/mgmol/LICENSE + +#include "Control.h" +#include "Electrostatic.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "Poisson.h" +#include "Potentials.h" +#include "mgmol_run.h" + +#include +#include +#include +#include +#include + +#include +namespace po = boost::program_options; + +template +int testRhoRestart(MGmolInterface* mgmol_) +{ + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + + MGmol* mgmol = static_cast*>(mgmol_); + std::shared_ptr> rho = mgmol->getRho(); + + /* save density from the restart file to elsewhere */ + std::vector rho0(rho->rho_[0].size()); + rho0 = rho->rho_[0]; + + /* recompute rho from the orbital */ + rho->update(*mgmol->getOrbitals()); + + /* check if the recomputed density is the same */ + for (int d = 0; d < (int)rho0.size(); d++) + { + double error = abs(rho0[d] - rho->rho_[0][d]) / abs(rho0[d]); + if (error > 1e-10) + { + printf("rank %d, rho[%d]=%.15e, rho0[%d]=%.15e\n", rank, d, + rho->rho_[0][d], d, rho0[d]); + std::cerr << "Density is inconsistent!!!" << std::endl; + return -1; + } + } + if (rank == 0) std::cout << "Density is consistent..." << std::endl; + + return 0; +} + +template +int testPotRestart(MGmolInterface* mgmol_) +{ + Control& ct = *(Control::instance()); + + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + + MGmol* mgmol = static_cast*>(mgmol_); + Potentials& pot = mgmol->getHamiltonian()->potential(); + Poisson* poisson = mgmol->electrostat_->getPoissonSolver(); + std::shared_ptr> rho = mgmol->getRho(); + + /* GridFunc initialization inputs */ + short bc[3]; + for (int d = 0; d < 3; d++) + bc[d] = ct.bcPoisson[d]; + + /* save potential from the restart file to elsewhere */ + pb::GridFunc vh0_gf(mygrid, bc[0], bc[1], bc[2]); + vh0_gf.assign((pot.vh_rho()).data(), 'd'); + double n = vh0_gf.norm2(); + std::cout << "Norm2 of vh = " << n << std::endl; + + std::vector vh0(pot.size()); + const std::vector& d_vhrho(pot.vh_rho()); + for (int d = 0; d < (int)vh0.size(); d++) + vh0[d] = d_vhrho[d]; + + /* recompute potential */ + pb::GridFunc grho(mygrid, bc[0], bc[1], bc[2]); + grho.assign(&rho->rho_[0][0]); + pb::GridFunc* grhoc = mgmol->electrostat_->getRhoc(); + + poisson->solve(grho, *grhoc); + const pb::GridFunc& vh(poisson->vh()); + + pb::GridFunc error_gf(vh0_gf); + error_gf -= vh; + + double rel_error = error_gf.norm2() / vh0_gf.norm2(); + if (rank == 0) + { + printf("FOM potential relative error: %.3e\n", rel_error); + } + if (rel_error > 1e-9) + { + if (rank == 0) + { + std::cerr << "Potential is inconsistent!!!" << std::endl; + } + return -1; + } + if (rank == 0) std::cout << "Potential is consistent..." << std::endl; + + return 0; +} + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + /* + * Initialize general things, like magma, openmp, IO, ... + */ + mgmol_init(comm); + + /* + * read runtime parameters + */ + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + /* + * Setup control struct with run time parameters + */ + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + int status = 0; + + // Enter main scope + { + MGmolInterface* mgmol = new MGmol(global_comm, + *MPIdata::sout, input_filename, lrs_filename, constraints_filename); + + mgmol->setup(); + + /* load a restart file */ + MGmol* mgmol_ext + = dynamic_cast*>(mgmol); + mgmol_ext->loadRestartFile(ct.restart_file); + + if (MPIdata::onpe0) + std::cout << "=============================" << std::endl; + if (MPIdata::onpe0) std::cout << "testRhoRestart..." << std::endl; + status = testRhoRestart(mgmol); + if (status < 0) return status; + + if (MPIdata::onpe0) + std::cout << "=============================" << std::endl; + if (MPIdata::onpe0) std::cout << "testPotRestart..." << std::endl; + status = testPotRestart(mgmol); + if (status < 0) return status; + + delete mgmol; + + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + time_t tt; + time(&tt); + if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; + + return 0; +}