diff --git a/source/source_esolver/esolver_double_xc.cpp b/source/source_esolver/esolver_double_xc.cpp index b6f4ceaa19..4fa764622a 100644 --- a/source/source_esolver/esolver_double_xc.cpp +++ b/source/source_esolver/esolver_double_xc.cpp @@ -1,7 +1,6 @@ #include "esolver_double_xc.h" #include "source_hamilt/module_ewald/H_Ewald_pw.h" -#include "source_hamilt/module_vdw/vdw.h" #include "source_hamilt/module_xc/xc_functional.h" #ifdef __MLALGO #include "source_lcao/module_deepks/LCAO_deepks.h" @@ -121,13 +120,9 @@ void ESolver_DoubleXC::before_scf(UnitCell& ucell, const int istep) ESolver_KS_LCAO::before_scf(ucell, istep); //---------------------------------------------------------- - //! calculate D2 or D3 vdW + //! Reuse the vdW correction prepared by ESolver_FP::before_scf. //---------------------------------------------------------- - auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); - if (vdw_solver != nullptr) - { - this->pelec_base->f_en.evdw = vdw_solver->get_energy(); - } + this->pelec_base->f_en.evdw = this->pelec->f_en.evdw; //---------------------------------------------------------- //! calculate ewald energy @@ -398,6 +393,7 @@ void ESolver_DoubleXC::cal_force(BaseCell& basecell, ModuleBase::matrix& this->deepks.dpks_out_type = "base"; // for deepks method fsl.getForceStress(ucell, + this->get_vdw_result(), PARAM.inp.cal_force, PARAM.inp.cal_stress, PARAM.inp.test_force, diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp index 76dc665d6b..754471313d 100644 --- a/source/source_esolver/esolver_fp.cpp +++ b/source/source_esolver/esolver_fp.cpp @@ -190,11 +190,18 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) GlobalV::ofs_running, GlobalV::ofs_warning); } - //! calculate D2 or D3 vdW + //! Evaluate the vdW correction once for this ionic configuration. + this->vdw_result_.reset(); auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) { - this->pelec->f_en.evdw = vdw_solver->get_energy(); + const vdw::VdwRequest request(PARAM.inp.cal_force, PARAM.inp.cal_stress); + this->vdw_result_.reset(new vdw::VdwResult(vdw_solver->evaluate(request))); + this->pelec->f_en.evdw = this->vdw_result_->energy; + } + else + { + this->pelec->f_en.evdw = 0.0; } //! calculate ewald energy diff --git a/source/source_esolver/esolver_fp.h b/source/source_esolver/esolver_fp.h index afeb1248e1..13e2bd3bd4 100644 --- a/source/source_esolver/esolver_fp.h +++ b/source/source_esolver/esolver_fp.h @@ -12,6 +12,12 @@ #include "source_pw/module_pwdft/vl_pw.h" // local pseudopotential #include +#include + +namespace vdw +{ +struct VdwResult; +} //! The First-Principles (FP) Energy Solver Class /** @@ -42,6 +48,11 @@ class ESolver_FP : public ESolver virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver); + const vdw::VdwResult* get_vdw_result() const { return this->vdw_result_.get(); } + + //! vdW correction evaluated once for the current ionic configuration. + std::unique_ptr vdw_result_; + //! These pointers will be deleted in the free_pointers() function every ion step. elecstate::ElecState* pelec = nullptr; ///< Electronic states diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp index 666ea1278d..741fa3fc1f 100644 --- a/source/source_esolver/esolver_ks_lcao.cpp +++ b/source/source_esolver/esolver_ks_lcao.cpp @@ -248,7 +248,7 @@ void ESolver_KS_LCAO::cal_force(BaseCell& basecell, ModuleBase::matrix& deepks.dpks_out_type = "tot"; // for deepks method - fsl.getForceStress(ucell, PARAM.inp.cal_force, PARAM.inp.cal_stress, + fsl.getForceStress(ucell, this->get_vdw_result(), PARAM.inp.cal_force, PARAM.inp.cal_stress, PARAM.inp.test_force, PARAM.inp.test_stress, this->gd, this->pv, this->pelec, this->dmat, this->psi, two_center_bundle_, orb_, force, this->scs, diff --git a/source/source_esolver/esolver_ks_pw.cpp b/source/source_esolver/esolver_ks_pw.cpp index 387137e349..d5cedd6123 100644 --- a/source/source_esolver/esolver_ks_pw.cpp +++ b/source/source_esolver/esolver_ks_pw.cpp @@ -369,6 +369,7 @@ void ESolver_KS_PW::cal_force(BaseCell& basecell, ModuleBase::matrix& // Calculate forces ff.cal_force(ucell, force, + this->get_vdw_result(), *this->pelec, this->pw_rhod, &ucell.symm, @@ -395,6 +396,7 @@ void ESolver_KS_PW::cal_stress(BaseCell& basecell, ModuleBase::matrix ss.cal_stress(stress, ucell, + this->get_vdw_result(), this->dftu, this->locpp, this->ppcell, diff --git a/source/source_esolver/esolver_of.cpp b/source/source_esolver/esolver_of.cpp index 0539e4a81b..6437ace50c 100644 --- a/source/source_esolver/esolver_of.cpp +++ b/source/source_esolver/esolver_of.cpp @@ -558,7 +558,8 @@ void ESolver_OF::cal_force(BaseCell& basecell, ModuleBase::matrix& force) // here nullptr is for DFT+U, which may cause bugs, mohan note 2025-11-07 // solvent can be used? mohan ask 2025-11-07 - ff.cal_force(ucell, force, *pelec, this->pw_rho, &ucell.symm, &sf, this->solvent, nullptr, &this->locpp); + ff.cal_force(ucell, force, this->get_vdw_result(), *pelec, this->pw_rho, &ucell.symm, &sf, + this->solvent, nullptr, &this->locpp); } /** @@ -577,6 +578,6 @@ void ESolver_OF::cal_stress(BaseCell& basecell, ModuleBase::matrix& stress) this->pphi_, this->pw_rho, kinetic_stress_); // kinetic stress OF_Stress_PW ss(this->pelec, this->pw_rho); - ss.cal_stress(stress, kinetic_stress_, ucell, &ucell.symm, this->locpp, &sf, &kv); + ss.cal_stress(stress, kinetic_stress_, ucell, this->get_vdw_result(), &ucell.symm, this->locpp, &sf, &kv); } } // namespace ModuleESolver diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp index a5873c80d9..01aa74ec27 100644 --- a/source/source_hamilt/module_vdw/test/vdw_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp @@ -25,8 +25,8 @@ * - vdw::make_vdw(): * Based on the value of INPUT.vdw_method, construct * Vdwd2 or Vdwd3 class, and do the initialization. -* - vdw::get_energy()/vdw::get_force()/vdw::get_stress(): -* Calculate the VDW (d2, d3_0 and d3_bj types) enerygy, force, stress. +* - vdw::Vdw::evaluate(): +* Calculate the requested vdW energy, force and stress in one evaluation. * - Vdwd2Parameters::initial_parameters() * - Vdwd3Parameters::initial_parameters() */ @@ -313,21 +313,26 @@ TEST_F(vdwd2Test, D2R0ZeroQuit) vdwd2_test.parameter().R0_["Si"] = 0.0; testing::internal::CaptureStdout(); - EXPECT_EXIT(vdwd2_test.get_energy(), ::testing::ExitedWithCode(1), ""); + EXPECT_EXIT(vdwd2_test.evaluate(vdw::VdwRequest(false, false)), ::testing::ExitedWithCode(1), ""); std::string output = testing::internal::GetCapturedStdout(); } TEST_F(vdwd2Test, D2GetEnergy) { auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene,-0.034526673470525196,1E-10); } TEST_F(vdwd2Test, D2GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.034526673470525196, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, -0.00078824525563651242,1e-12); EXPECT_NEAR(force[0].y, 2.6299822052061785e-08,1e-12); EXPECT_NEAR(force[0].z, 2.6299822050796364e-08,1e-12); @@ -339,7 +344,11 @@ TEST_F(vdwd2Test, D2GetForce) TEST_F(vdwd2Test, D2GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.034526673470525196, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, -0.00020532319044269705,1e-12); EXPECT_NEAR(stress.e12, -3.5642821939401251e-08,1e-12); EXPECT_NEAR(stress.e13, -3.5642821939437223e-08,1e-12); @@ -433,14 +442,19 @@ TEST_F(vdwd3Test, D30Period) TEST_F(vdwd3Test, D30GetEnergy) { auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene,-0.20932367230529664,1E-10); } TEST_F(vdwd3Test, D30GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.20932367230529664, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, -0.032450975169023302,1e-12); EXPECT_NEAR(force[0].y, 0.0,1e-12); EXPECT_NEAR(force[0].z, 0.0,1e-12); @@ -452,7 +466,11 @@ TEST_F(vdwd3Test, D30GetForce) TEST_F(vdwd3Test, D30GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.20932367230529664, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, -0.0011141545452036336,1e-12); EXPECT_NEAR(stress.e12, 0.0,1e-12); EXPECT_NEAR(stress.e13, 0.0,1e-12); @@ -468,7 +486,8 @@ TEST_F(vdwd3Test, D3bjGetEnergy) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene,-0.047458675421836918,1E-10); } @@ -476,7 +495,11 @@ TEST_F(vdwd3Test, D3bjGetForce) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.047458675421836918, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, -0.0026006968781200602,1e-12); EXPECT_NEAR(force[0].y, 0.0,1e-12); EXPECT_NEAR(force[0].z, 0.0,1e-12); @@ -489,7 +512,11 @@ TEST_F(vdwd3Test, D3bjGetStress) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.047458675421836918, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, -0.00014376286737216365,1e-12); EXPECT_NEAR(stress.e12, 0.0,1e-12); EXPECT_NEAR(stress.e13, 0.0,1e-12); @@ -538,14 +565,19 @@ class vdwd3abcTest: public testing::Test TEST_F(vdwd3abcTest, D30GetEnergy) { auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene,-0.11487062308916372,1E-10); } TEST_F(vdwd3abcTest, D30GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.11487062308916372, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, 0.030320738678429094,1e-12); EXPECT_NEAR(force[0].y, 0.025570534655235538,1e-12); EXPECT_NEAR(force[0].z, 0.025570534655235538,1e-12); @@ -557,7 +589,11 @@ TEST_F(vdwd3abcTest, D30GetForce) TEST_F(vdwd3abcTest, D30GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.11487062308916372, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, -0.00023421562840819491,1e-12); EXPECT_NEAR(stress.e12, -0.00015112406243413323,1e-12); EXPECT_NEAR(stress.e13, -0.00015112406243413302,1e-12); @@ -573,7 +609,8 @@ TEST_F(vdwd3abcTest, D3bjGetEnergy) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene,-0.030667806197006021,1E-10); } @@ -581,7 +618,11 @@ TEST_F(vdwd3abcTest, D3bjGetForce) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.030667806197006021, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, -0.0010630099217696475,1e-12); EXPECT_NEAR(force[0].y, -0.0010031953309458587,1e-12); EXPECT_NEAR(force[0].z, -0.0010031953309458642,1e-12); @@ -594,7 +635,11 @@ TEST_F(vdwd3abcTest, D3bjGetStress) { input.vdw_method = "d3_bj"; auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.030667806197006021, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, -3.3803329202372578e-05,1e-12); EXPECT_NEAR(stress.e12, 5.1291622417145846e-06,1e-12); EXPECT_NEAR(stress.e13, 5.1291622417145889e-06,1e-12); @@ -643,7 +688,8 @@ class vdwd4Test: public testing::Test TEST_F(vdwd4Test, D4GetEnergy) { auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene, -0.04998837990336073, 1E-10); } @@ -652,14 +698,19 @@ TEST_F(vdwd4Test, D4GetEnergyForChargedSystem) input.nelec = 7.0; auto vdw_solver = vdw::make_vdw(ucell, input); - const double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene, -0.04359451765256733, 1E-10); } TEST_F(vdwd4Test, D4GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.04998837990336073, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, -0.0023357259921368717, 1e-12); EXPECT_NEAR(force[0].y, 0.0, 1e-12); EXPECT_NEAR(force[0].z, 0.0, 1e-12); @@ -671,7 +722,11 @@ TEST_F(vdwd4Test, D4GetForce) TEST_F(vdwd4Test, D4GetStress) { auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.04998837990336073, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, 0.00015830384474877792, 1e-12); EXPECT_NEAR(stress.e12, 0.0, 1e-12); EXPECT_NEAR(stress.e13, 0.0, 1e-12); @@ -687,7 +742,8 @@ TEST_F(vdwd4Test, D4SGetEnergy) { input.vdw_d4_model = "d4s"; auto vdw_solver = vdw::make_vdw(ucell, input); - double ene = vdw_solver->get_energy(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(false, false)); + const double ene = result.energy; EXPECT_NEAR(ene, -0.05638517144755526, 1E-10); } @@ -695,7 +751,11 @@ TEST_F(vdwd4Test, D4SGetForce) { input.vdw_d4_model = "d4s"; auto vdw_solver = vdw::make_vdw(ucell, input); - std::vector> force = vdw_solver->get_force(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, false)); + EXPECT_NEAR(result.energy, -0.05638517144755526, 1E-10); + ASSERT_TRUE(result.has_force); + EXPECT_FALSE(result.has_stress); + const std::vector>& force = result.force; EXPECT_NEAR(force[0].x, -0.005448661796788402, 1e-12); EXPECT_NEAR(force[0].y, 0.0, 1e-12); EXPECT_NEAR(force[0].z, 0.0, 1e-12); @@ -708,7 +768,11 @@ TEST_F(vdwd4Test, D4SGetStress) { input.vdw_d4_model = "d4s"; auto vdw_solver = vdw::make_vdw(ucell, input); - ModuleBase::Matrix3 stress = vdw_solver->get_stress(); + const vdw::VdwResult result = vdw_solver->evaluate(vdw::VdwRequest(true, true)); + EXPECT_NEAR(result.energy, -0.05638517144755526, 1E-10); + ASSERT_TRUE(result.has_force); + ASSERT_TRUE(result.has_stress); + const ModuleBase::Matrix3& stress = result.stress; EXPECT_NEAR(stress.e11, 0.00013831119855416262, 1e-12); EXPECT_NEAR(stress.e12, 0.0, 1e-12); EXPECT_NEAR(stress.e13, 0.0, 1e-12); diff --git a/source/source_hamilt/module_vdw/vdw.h b/source/source_hamilt/module_vdw/vdw.h index 1aeb4618d1..747aa5d3f1 100644 --- a/source/source_hamilt/module_vdw/vdw.h +++ b/source/source_hamilt/module_vdw/vdw.h @@ -2,7 +2,9 @@ #define VDW_H #include +#include #include + #include "source_cell/unitcell.h" #include "vdw_parameters.h" #include "vdwd2_parameters.h" @@ -11,54 +13,64 @@ namespace vdw { - template - std::unique_ptr make_unique(Args &&... args) { - return std::unique_ptr(new T(std::forward(args)...)); - +template +std::unique_ptr make_unique(Args&&... args) +{ + return std::unique_ptr(new T(std::forward(args)...)); } +struct VdwRequest +{ + VdwRequest(const bool force_in, const bool stress_in) : force(force_in), stress(stress_in) {} + + bool force; + bool stress; +}; + +struct VdwResult +{ + VdwResult() : energy(0.0), has_force(false), has_stress(false) + { + stress.Zero(); + } + + double energy; + std::vector> force; + ModuleBase::Matrix3 stress; + bool has_force; + bool has_stress; +}; + class Vdw { public: - Vdw(const UnitCell &unit_in) : ucell_(unit_in) {}; + Vdw(const UnitCell& unit_in) : ucell_(unit_in) {} virtual ~Vdw() = default; - inline double get_energy(bool cal=true) { - if (cal) { cal_energy(); } - return energy_; - } - inline const std::vector> &get_force(bool cal=true) { - if (cal) { cal_force(); } - return force_; - } - inline const ModuleBase::Matrix3 &get_stress(bool cal=true) { - if (cal) { cal_stress(); } - return stress_; + VdwResult evaluate(const VdwRequest& request) + { + VdwResult result; + evaluate_impl(request, result); + return result; } protected: - const UnitCell &ucell_; - - double energy_ = 0; - std::vector> force_; - ModuleBase::Matrix3 stress_; + const UnitCell& ucell_; - virtual void cal_energy() = 0; - virtual void cal_force() = 0; - virtual void cal_stress() = 0; + virtual void evaluate_impl(const VdwRequest& request, VdwResult& result) = 0; }; /** * @brief make vdw correction object - * + * * @param ucell UnitCell instance * @param input Parameter instance * @param plog optional, for logging the parameter setting process - * @return std::unique_ptr + * @return std::unique_ptr */ -std::unique_ptr make_vdw(const UnitCell &ucell, - const Input_para &input, +std::unique_ptr make_vdw(const UnitCell& ucell, + const Input_para& input, std::ofstream* plog = nullptr); } // namespace vdw diff --git a/source/source_hamilt/module_vdw/vdwd2.cpp b/source/source_hamilt/module_vdw/vdwd2.cpp index 74a48633f6..9489ea62a7 100644 --- a/source/source_hamilt/module_vdw/vdwd2.cpp +++ b/source/source_hamilt/module_vdw/vdwd2.cpp @@ -7,89 +7,91 @@ #include "vdwd2.h" #include "source_base/timer.h" +#include +#include + namespace vdw { -void Vdwd2::cal_energy() +void Vdwd2::evaluate_impl(const VdwRequest& request, VdwResult& result) { - ModuleBase::TITLE("Vdwd2", "energy"); - ModuleBase::timer::start("Vdwd2", "energy"); - para_.initset(ucell_); - energy_ = 0; - - auto energy = [&](double r, - double R0_sum, - double C6_product, - double r_sqr, - int, - int, - const ModuleBase::Vector3 &, - const ModuleBase::Vector3 &) { - const double tmp_damp_recip = 1 + exp(-para_.damping() * (r / R0_sum - 1)); - energy_ -= C6_product / pow(r_sqr, 3) / tmp_damp_recip / 2; - }; - index_loops(energy); - energy_ *= para_.scaling(); - ModuleBase::timer::end("Vdwd2", "energy"); -} + ModuleBase::TITLE("Vdwd2", "evaluate"); + ModuleBase::timer::start("Vdwd2", "evaluate"); -void Vdwd2::cal_force() -{ - ModuleBase::TITLE("Vdwd2", "force"); - ModuleBase::timer::start("Vdwd2", "force"); para_.initset(ucell_); - force_.clear(); - force_.resize(ucell_.nat); - - auto force = [&](double r, - double R0_sum, - double C6_product, - double r_sqr, - int it1, - int ia1, - const ModuleBase::Vector3 &tau1, - const ModuleBase::Vector3 &tau2) { - const double tmp_exp = exp(-para_.damping() * (r / R0_sum - 1)); - const double tmp_factor = C6_product / pow(r_sqr, 3) / r / (1 + tmp_exp) - * (-6 / r + tmp_exp / (1 + tmp_exp) * para_.damping() / R0_sum); - force_[ucell_.itia2iat(it1, ia1)] += tmp_factor * (tau1 - tau2); - }; - index_loops(force); - std::for_each(force_.begin(), force_.end(), [&](ModuleBase::Vector3 &f) { - f *= para_.scaling() / ucell_.lat0; - }); - ModuleBase::timer::end("Vdwd2", "force"); -} + if (request.force) + { + result.force.resize(ucell_.nat); + } + if (request.stress) + { + result.stress.Zero(); + } -void Vdwd2::cal_stress() -{ - ModuleBase::TITLE("Vdwd2", "stress"); - ModuleBase::timer::start("Vdwd2", "stress"); - para_.initset(ucell_); - stress_.Zero(); - - auto stress = [&](double r, - double R0_sum, - double C6_product, - double r_sqr, - int it1, - int ia1, - const ModuleBase::Vector3 &tau1, - const ModuleBase::Vector3 &tau2) { + const bool need_derivatives = request.force || request.stress; + auto evaluate_pair = [&](double r, + double R0_sum, + double C6_product, + double r_sqr, + int it1, + int ia1, + const ModuleBase::Vector3& tau1, + const ModuleBase::Vector3& tau2) { const double tmp_exp = exp(-para_.damping() * (r / R0_sum - 1)); - const double tmp_factor = C6_product / pow(r_sqr, 3) / r / (1 + tmp_exp) - * (-6 / r + tmp_exp / (1 + tmp_exp) * para_.damping() / R0_sum); - const ModuleBase::Vector3 dr = tau2 - tau1; - stress_ += tmp_factor / 2 - * ModuleBase::Matrix3(dr.x * dr.x, dr.x * dr.y, dr.x * dr.z, - dr.y * dr.x, dr.y * dr.y, dr.y * dr.z, - dr.z * dr.x, dr.z * dr.y, dr.z * dr.z); + const double tmp_damp_recip = 1.0 + tmp_exp; + + result.energy -= C6_product / pow(r_sqr, 3) / tmp_damp_recip / 2.0; + + if (!need_derivatives) + { + return; + } + + const double tmp_factor = C6_product / pow(r_sqr, 3) / r / tmp_damp_recip + * (-6.0 / r + + tmp_exp / tmp_damp_recip * para_.damping() / R0_sum); + + if (request.force) + { + result.force[ucell_.itia2iat(it1, ia1)] += tmp_factor * (tau1 - tau2); + } + + if (request.stress) + { + const ModuleBase::Vector3 dr = tau2 - tau1; + result.stress += tmp_factor / 2.0 + * ModuleBase::Matrix3(dr.x * dr.x, + dr.x * dr.y, + dr.x * dr.z, + dr.y * dr.x, + dr.y * dr.y, + dr.y * dr.z, + dr.z * dr.x, + dr.z * dr.y, + dr.z * dr.z); + } }; - index_loops(stress); - stress_ *= para_.scaling() / ucell_.omega; - ModuleBase::timer::end("Vdwd2", "stress"); + index_loops(evaluate_pair); + + result.energy *= para_.scaling(); + + if (request.force) + { + std::for_each(result.force.begin(), result.force.end(), [&](ModuleBase::Vector3& force) { + force *= para_.scaling() / ucell_.lat0; + }); + result.has_force = true; + } + + if (request.stress) + { + result.stress *= para_.scaling() / ucell_.omega; + result.has_stress = true; + } + + ModuleBase::timer::end("Vdwd2", "evaluate"); } } // namespace vdw diff --git a/source/source_hamilt/module_vdw/vdwd2.h b/source/source_hamilt/module_vdw/vdwd2.h index 4e1e836d57..12b119b03c 100644 --- a/source/source_hamilt/module_vdw/vdwd2.h +++ b/source/source_hamilt/module_vdw/vdwd2.h @@ -16,21 +16,20 @@ class Vdwd2 : public Vdw { public: - Vdwd2(const UnitCell &unit_in) : Vdw(unit_in) {} + Vdwd2(const UnitCell& unit_in) : Vdw(unit_in) {} ~Vdwd2() = default; - Vdwd2Parameters ¶meter() { return para_; } - const Vdwd2Parameters ¶meter() const { return para_; } + Vdwd2Parameters& parameter() { return para_; } + const Vdwd2Parameters& parameter() const { return para_; } private: Vdwd2Parameters para_; - void cal_energy() override; - void cal_force() override; - void cal_stress() override; + void evaluate_impl(const VdwRequest& request, VdwResult& result) override; - template void index_loops(F &&f) + template + void index_loops(F&& f) { int xidx = para_.period().x / 2; int yidx = para_.period().y / 2; @@ -44,7 +43,8 @@ class Vdwd2 : public Vdw = sqrt(para_.C6().at(ucell_.atoms[it1].ncpp.psd) * para_.C6().at(ucell_.atoms[it2].ncpp.psd)) / pow(ucell_.lat0, 6); const double R0_sum - = (para_.R0().at(ucell_.atoms[it1].ncpp.psd) + para_.R0().at(ucell_.atoms[it2].ncpp.psd)) / ucell_.lat0; + = (para_.R0().at(ucell_.atoms[it1].ncpp.psd) + para_.R0().at(ucell_.atoms[it2].ncpp.psd)) + / ucell_.lat0; if (!R0_sum) { ModuleBase::WARNING_QUIT("Input", "R0_sum can not be 0"); @@ -61,21 +61,23 @@ class Vdwd2 : public Vdw { for (ilat_loop.z = -zidx; ilat_loop.z <= zidx; ++ilat_loop.z) { - if ((!(ilat_loop.x || ilat_loop.y || ilat_loop.z)) && (it1 == it2) && (ia1 == ia2)) + if ((!(ilat_loop.x || ilat_loop.y || ilat_loop.z)) && (it1 == it2) + && (ia1 == ia2)) + { continue; + } const ModuleBase::Vector3 tau2 = ucell_.atoms[it2].tau[ia2] + ilat_loop * ucell_.latvec; const double r_sqr = (tau1 - tau2).norm2(); const double r = sqrt(r_sqr); - // calculations happen in f f(r, R0_sum, C6_product, r_sqr, it1, ia1, tau1, tau2); } } - } // end for ilat_loop - } // end for ia2 - } // end for ia1 - } // end for it2 - } // end for it1 + } + } + } + } + } } }; diff --git a/source/source_hamilt/module_vdw/vdwd3.cpp b/source/source_hamilt/module_vdw/vdwd3.cpp index 47d035311e..1db27de3e1 100644 --- a/source/source_hamilt/module_vdw/vdwd3.cpp +++ b/source/source_hamilt/module_vdw/vdwd3.cpp @@ -22,6 +22,8 @@ void Vdwd3::init() lat_[2] = ucell_.a3 * ucell_.lat0; std::vector at_kind = atom_kind(); + iz_.clear(); + xyz_.clear(); iz_.reserve(ucell_.nat); xyz_.reserve(ucell_.nat); for (size_t it = 0; it != ucell_.ntype; it++) { @@ -82,10 +84,10 @@ std::vector Vdwd3::atom_kind() return atom_kind; } -void Vdwd3::cal_energy() +void Vdwd3::evaluate_energy(double& energy) { - ModuleBase::TITLE("Vdwd3", "cal_energy"); - ModuleBase::timer::start("Vdwd3", "cal_energy"); + ModuleBase::TITLE("Vdwd3", "evaluate_energy"); + ModuleBase::timer::start("Vdwd3", "evaluate_energy"); init(); int ij = 0; @@ -257,51 +259,53 @@ void Vdwd3::cal_energy() { pbc_three_body(iz_, lat_, xyz_, rep_cn_, cc6ab, eabc); } - energy_ = (-para_.s6() * e6 - para_.s18() * e8 - eabc) * 2; - ModuleBase::timer::end("Vdwd3", "cal_energy"); + energy = (-para_.s6() * e6 - para_.s18() * e8 - eabc) * 2.0; + ModuleBase::timer::end("Vdwd3", "evaluate_energy"); } -void Vdwd3::cal_force() +void Vdwd3::evaluate_impl(const VdwRequest& request, VdwResult& result) { - ModuleBase::TITLE("Vdwd3", "cal_force"); - ModuleBase::timer::start("Vdwd3", "cal_force"); - init(); - - force_.clear(); - force_.resize(ucell_.nat); - - std::vector> g; - g.clear(); - g.resize(ucell_.nat); - ModuleBase::matrix smearing_sigma(3, 3); - - pbc_gdisp(g, smearing_sigma); - - for (size_t iat = 0; iat != ucell_.nat; iat++) { - force_[iat] = -2.0 * g[iat]; -} + if (!request.force && !request.stress) + { + evaluate_energy(result.energy); + return; + } - ModuleBase::timer::end("Vdwd3", "cal_force"); -} + ModuleBase::TITLE("Vdwd3", "evaluate"); + ModuleBase::timer::start("Vdwd3", "evaluate"); -void Vdwd3::cal_stress() -{ - ModuleBase::TITLE("Vdwd3", "cal_stress"); - ModuleBase::timer::start("Vdwd3", "cal_stress"); init(); - std::vector> g; - g.clear(); - g.resize(ucell_.nat); + std::vector> gradient(ucell_.nat); ModuleBase::matrix smearing_sigma(3, 3); + pbc_gdisp(gradient, smearing_sigma, result.energy); - pbc_gdisp(g, smearing_sigma); + if (request.force) + { + result.force.resize(ucell_.nat); + for (int iat = 0; iat < ucell_.nat; ++iat) + { + result.force[iat] = -2.0 * gradient[iat]; + } + result.has_force = true; + } + + if (request.stress) + { + result.stress = ModuleBase::Matrix3(2.0 * smearing_sigma(0, 0), + 2.0 * smearing_sigma(0, 1), + 2.0 * smearing_sigma(0, 2), + 2.0 * smearing_sigma(1, 0), + 2.0 * smearing_sigma(1, 1), + 2.0 * smearing_sigma(1, 2), + 2.0 * smearing_sigma(2, 0), + 2.0 * smearing_sigma(2, 1), + 2.0 * smearing_sigma(2, 2)) + / ucell_.omega; + result.has_stress = true; + } - stress_ = ModuleBase::Matrix3(2.0 * smearing_sigma(0, 0), 2.0 * smearing_sigma(0, 1), 2.0 * smearing_sigma(0, 2), - 2.0 * smearing_sigma(1, 0), 2.0 * smearing_sigma(1, 1), 2.0 * smearing_sigma(1, 2), - 2.0 * smearing_sigma(2, 0), 2.0 * smearing_sigma(2, 1), 2.0 * smearing_sigma(2, 2)) - / ucell_.omega; - ModuleBase::timer::end("Vdwd3", "cal_stress"); + ModuleBase::timer::end("Vdwd3", "evaluate"); } void Vdwd3::get_c6(int iat, int jat, double nci, double ncj, double &c6) @@ -735,8 +739,13 @@ void Vdwd3::get_dc6_dcnij(int mxci, int mxcj, double cni, double cnj, int izi, i } } -void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::matrix &smearing_sigma) +void Vdwd3::pbc_gdisp(std::vector>& g, + ModuleBase::matrix& smearing_sigma, + double& energy) { + double e6 = 0.0; + double e8 = 0.0; + double eabc = 0.0; std::vector c6save(ucell_.nat * (ucell_.nat + 1)), dc6_rest_sum(ucell_.nat * (ucell_.nat + 1) / 2), dc6i(ucell_.nat), cn(ucell_.nat); pbc_ncoord(cn); @@ -788,6 +797,9 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m t8 = std::pow(r / (para_.rs18() * r0), -para_.alp8()); damp8 = 1.0 / (1.0 + 6.0 * t8); + e6 += c6 * damp6 / r6 * 0.5; + e8 += 3.0 * c6 * r42 * damp8 / r8 * 0.5; + // d(r^(-6))/d(tau) drij[linii][taux + rep_vdw_[0]][tauy + rep_vdw_[1]][tauz + rep_vdw_[2]] += (-para_.s6() * (6.0 / (r7)*c6 * damp6) @@ -840,6 +852,9 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m t8 = std::pow(r / (para_.rs18() * r0), -para_.alp8()); damp8 = 1.0 / (1.0 + 6.0 * t8); + e6 += c6 * damp6 / r6; + e8 += 3.0 * c6 * r42 * damp8 / r8; + // d(r^(-6))/d(r_ij) drij[linij][taux + rep_vdw_[0]][tauy + rep_vdw_[1]][tauz + rep_vdw_[2]] += -para_.s6() * (6.0 / (r7)*c6 * damp6) - para_.s18() * (24.0 / (r9)*c6 * r42 * damp8); @@ -894,6 +909,9 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m t6 = r6 + std::pow(r0, 6); t8 = r8 + std::pow(r0, 8); + e6 += c6 / t6 * 0.5; + e8 += 3.0 * c6 * r42 / t8 * 0.5; + // d(1/r^(-6)+r0^6)/d(r) drij[linii][taux + rep_vdw_[0]][tauy + rep_vdw_[1]][tauz + rep_vdw_[2]] += -para_.s6() * c6 * 6.0 * r4 * r / (t6 * t6) * 0.5 @@ -939,6 +957,9 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m t6 = r6 + std::pow(r0, 6); t8 = r8 + std::pow(r0, 8); + e6 += c6 / t6; + e8 += 3.0 * c6 * r42 / t8; + drij[linij][taux + rep_vdw_[0]][tauy + rep_vdw_[1]][tauz + rep_vdw_[2]] += -para_.s6() * c6 * 6.0 * r4 * r / (t6 * t6) - para_.s18() * c6 * 24.0 * r42 * r7 / (t8 * t8); @@ -1027,6 +1048,7 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m ang = 0.375 * (rij2 + rjk2 - rik2) * (rij2 - rjk2 + rik2) * (-rij2 + rjk2 + rik2) / (geomean3 * geomean2) + 1.0 / geomean3; + eabc += ang * c9 * damp9; dc6_rest = ang * damp9; dfdmp = 2.0 * alp9 * std::pow(0.75 * r0av, alp9) * damp9 * damp9; @@ -1149,6 +1171,7 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m ang = 0.375 * (rij2 + rjk2 - rik2) * (rij2 - rjk2 + rik2) * (-rij2 + rjk2 + rik2) / (geomean3 * geomean2) + 1.0 / geomean3; + eabc += ang * c9 * damp9 / 2.0; dc6_rest = ang * damp9 / 2.0; dfdmp = 2.0 * alp9 * std::pow(0.75 * r0av, alp9) * damp9 * damp9; @@ -1269,6 +1292,7 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m ang = 0.375 * (rij2 + rjk2 - rik2) * (rij2 - rjk2 + rik2) * (-rij2 + rjk2 + rik2) / (geomean3 * geomean2) + 1.0 / geomean3; + eabc += ang * c9 * damp9 / 2.0; dc6_rest = ang * damp9 / 2.0; dfdmp = 2.0 * alp9 * std::pow(0.75 * r0av, alp9) * damp9 * damp9; @@ -1394,6 +1418,7 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m ang = 0.375 * (rij2 + rjk2 - rik2) * (rij2 - rjk2 + rik2) * (-rij2 + rjk2 + rik2) / (geomean3 * geomean2) + 1.0 / geomean3; + eabc += ang * c9 * damp9 / 6.0; dc6_rest = ang * damp9 / 6.0; dfdmp = 2.0 * alp9 * std::pow(0.75 * r0av, alp9) * damp9 * damp9; @@ -1533,6 +1558,8 @@ void Vdwd3::pbc_gdisp(std::vector> &g, ModuleBase::m } } } // end iat + + energy = (-para_.s6() * e6 - para_.s18() * e8 - eabc) * 2.0; } } // namespace vdw diff --git a/source/source_hamilt/module_vdw/vdwd3.h b/source/source_hamilt/module_vdw/vdwd3.h index f65f8f66b4..207b9f94c2 100644 --- a/source/source_hamilt/module_vdw/vdwd3.h +++ b/source/source_hamilt/module_vdw/vdwd3.h @@ -32,9 +32,9 @@ class Vdwd3 : public Vdw std::vector rep_vdw_; std::vector rep_cn_; - void cal_energy() override; - void cal_force() override; - void cal_stress() override; + void evaluate_impl(const VdwRequest& request, VdwResult& result) override; + + void evaluate_energy(double& energy); void init(); @@ -53,7 +53,9 @@ class Vdwd3 : public Vdw const std::vector &cc6ab, double &eabc); - void pbc_gdisp(std::vector> &g, ModuleBase::matrix &smearing_sigma); + void pbc_gdisp(std::vector>& g, + ModuleBase::matrix& smearing_sigma, + double& energy); void get_dc6_dcnij(int mxci, int mxcj, double cni, double cnj, int izi, int izj, int iat, int jat, double &c6check, double &dc6i, double &dc6j); diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp index 1b8c7d3b50..8f66b68b2c 100644 --- a/source/source_hamilt/module_vdw/vdwd4.cpp +++ b/source/source_hamilt/module_vdw/vdwd4.cpp @@ -219,90 +219,74 @@ void Vdwd4::compute(double& energy_ha, #endif } -void Vdwd4::cal_energy() +void Vdwd4::set_force_from_gradient(const std::vector& gradient_ha_bohr, + VdwResult& result) const { - ModuleBase::TITLE("Vdwd4", "cal_energy"); - ModuleBase::timer::start("Vdwd4", "cal_energy"); - - double energy_ha = 0.0; - compute(energy_ha, nullptr, nullptr); - - // DFT-D4 returns Hartree; ABACUS vdW energies are stored in Ry. - energy_ = 2.0 * energy_ha; - - ModuleBase::timer::end("Vdwd4", "cal_energy"); -} - -void Vdwd4::set_force_from_gradient(const std::vector& gradient_ha_bohr) -{ - force_.clear(); - force_.resize(ucell_.nat); + result.force.resize(ucell_.nat); for (int iat = 0; iat < ucell_.nat; ++iat) { // DFT-D4 returns dE/dR in Ha/Bohr; ABACUS forces are -dE/dR in Ry/Bohr. - force_[iat].x = -2.0 * gradient_ha_bohr[3 * iat + 0]; - force_[iat].y = -2.0 * gradient_ha_bohr[3 * iat + 1]; - force_[iat].z = -2.0 * gradient_ha_bohr[3 * iat + 2]; + result.force[iat].x = -2.0 * gradient_ha_bohr[3 * iat + 0]; + result.force[iat].y = -2.0 * gradient_ha_bohr[3 * iat + 1]; + result.force[iat].z = -2.0 * gradient_ha_bohr[3 * iat + 2]; } - has_force_cache_ = true; + result.has_force = true; } -void Vdwd4::set_stress_from_sigma(const std::array& sigma_ha) +void Vdwd4::set_stress_from_sigma(const std::array& sigma_ha, + VdwResult& result) const { // Tentative mapping consistent with the current D3 convention. // Confirm sign, transposition and volume normalization by finite-strain tests. - stress_ = ModuleBase::Matrix3(2.0 * sigma_ha[0], 2.0 * sigma_ha[1], 2.0 * sigma_ha[2], - 2.0 * sigma_ha[3], 2.0 * sigma_ha[4], 2.0 * sigma_ha[5], - 2.0 * sigma_ha[6], 2.0 * sigma_ha[7], 2.0 * sigma_ha[8]) - / ucell_.omega; - - has_stress_cache_ = true; + result.stress = ModuleBase::Matrix3(2.0 * sigma_ha[0], + 2.0 * sigma_ha[1], + 2.0 * sigma_ha[2], + 2.0 * sigma_ha[3], + 2.0 * sigma_ha[4], + 2.0 * sigma_ha[5], + 2.0 * sigma_ha[6], + 2.0 * sigma_ha[7], + 2.0 * sigma_ha[8]) + / ucell_.omega; + result.has_stress = true; } -void Vdwd4::cal_force() +void Vdwd4::evaluate_impl(const VdwRequest& request, VdwResult& result) { - ModuleBase::TITLE("Vdwd4", "cal_force"); - ModuleBase::timer::start("Vdwd4", "cal_force"); + ModuleBase::TITLE("Vdwd4", "evaluate"); + ModuleBase::timer::start("Vdwd4", "evaluate"); - if (!has_force_cache_ || !has_stress_cache_) + double energy_ha = 0.0; + if (request.force || request.stress) { - double energy_ha = 0.0; std::vector gradient(3 * ucell_.nat, 0.0); std::array sigma; sigma.fill(0.0); - // Request sigma together with the gradient. The DFT-D4 C API computes - // sigma internally for gradient calculations anyway, so keep it and - // avoid a second expensive D4 call when ABACUS subsequently requests stress. + // The DFT-D4 C API evaluates energy, gradient and sigma together. + // Keep all requested quantities from this single call. compute(energy_ha, &gradient, &sigma); - set_force_from_gradient(gradient); - set_stress_from_sigma(sigma); - } - - ModuleBase::timer::end("Vdwd4", "cal_force"); -} - -void Vdwd4::cal_stress() -{ - ModuleBase::TITLE("Vdwd4", "cal_stress"); - ModuleBase::timer::start("Vdwd4", "cal_stress"); - if (!has_stress_cache_) + if (request.force) + { + set_force_from_gradient(gradient, result); + } + if (request.stress) + { + set_stress_from_sigma(sigma, result); + } + } + else { - double energy_ha = 0.0; - std::vector gradient(3 * ucell_.nat, 0.0); - std::array sigma; - sigma.fill(0.0); - - // DFT-D4 may require a valid gradient buffer when sigma is requested. - compute(energy_ha, &gradient, &sigma); - set_force_from_gradient(gradient); - set_stress_from_sigma(sigma); + compute(energy_ha, nullptr, nullptr); } - ModuleBase::timer::end("Vdwd4", "cal_stress"); + // DFT-D4 returns Hartree; ABACUS vdW energies are stored in Ry. + result.energy = 2.0 * energy_ha; + + ModuleBase::timer::end("Vdwd4", "evaluate"); } } // namespace vdw diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h index 2c580388b6..f035db3d0b 100644 --- a/source/source_hamilt/module_vdw/vdwd4.h +++ b/source/source_hamilt/module_vdw/vdwd4.h @@ -25,15 +25,10 @@ class Vdwd4 : public Vdw double cutoff_cn_ = 0.0; // Bohr, coordination-number cutoff double total_charge_ = 0.0; // e, total system charge (sum zv*na - nelec) - bool has_force_cache_ = false; - bool has_stress_cache_ = false; + void evaluate_impl(const VdwRequest& request, VdwResult& result) override; - void set_force_from_gradient(const std::vector& gradient_ha_bohr); - void set_stress_from_sigma(const std::array& sigma_ha); - - void cal_energy() override; - void cal_force() override; - void cal_stress() override; + void set_force_from_gradient(const std::vector& gradient_ha_bohr, VdwResult& result) const; + void set_stress_from_sigma(const std::array& sigma_ha, VdwResult& result) const; void build_structure(std::vector& numbers, std::vector& positions, diff --git a/source/source_lcao/FORCE_STRESS.cpp b/source/source_lcao/FORCE_STRESS.cpp index 43e78b82a0..6634a42303 100644 --- a/source/source_lcao/FORCE_STRESS.cpp +++ b/source/source_lcao/FORCE_STRESS.cpp @@ -6,6 +6,7 @@ #include "source_io/module_parameter/parameter.h" // new #include "source_base/timer.h" +#include "source_base/tool_quit.h" #include "source_cell/module_neighbor/sltk_grid_driver.h" #include "source_estate/elecstate_lcao.h" #include "source_estate/module_pot/H_TDDFT_pw.h" // Taoni add 2025-02-20 @@ -66,6 +67,7 @@ Force_Stress_LCAO::~Force_Stress_LCAO() } template void Force_Stress_LCAO::getForceStress(UnitCell& ucell, + const vdw::VdwResult* vdw_result, const bool isforce, const bool isstress, const bool istestf, @@ -347,23 +349,31 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, // jiyy add 2019-05-18, update 2021-05-02 ModuleBase::matrix force_vdw; ModuleBase::matrix stress_vdw; - auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp); - if (vdw_solver != nullptr) + if (vdw_result != nullptr) { if (isforce) { + if (!vdw_result->has_force || vdw_result->force.size() != static_cast(nat)) + { + ModuleBase::WARNING_QUIT("Force_Stress_LCAO::getForceStress", + "The cached vdW force is unavailable or has an invalid size."); + } force_vdw.create(nat, 3); - const std::vector>& force_vdw_temp = vdw_solver->get_force(); - for (int iat = 0; iat < ucell.nat; ++iat) + for (int iat = 0; iat < nat; ++iat) { - force_vdw(iat, 0) = force_vdw_temp[iat].x; - force_vdw(iat, 1) = force_vdw_temp[iat].y; - force_vdw(iat, 2) = force_vdw_temp[iat].z; + force_vdw(iat, 0) = vdw_result->force[iat].x; + force_vdw(iat, 1) = vdw_result->force[iat].y; + force_vdw(iat, 2) = vdw_result->force[iat].z; } } if (isstress) { - stress_vdw = vdw_solver->get_stress().to_matrix(); + if (!vdw_result->has_stress) + { + ModuleBase::WARNING_QUIT("Force_Stress_LCAO::getForceStress", + "The cached vdW stress is unavailable."); + } + stress_vdw = vdw_result->stress.to_matrix(); } } @@ -553,7 +563,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, } #endif // VDW force of vdwd2 or vdwd3 - if (vdw_solver != nullptr) + if (vdw_result != nullptr) { fcs(iat, i) += force_vdw(iat, i); } @@ -675,7 +685,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, ModuleIO::print_force(GlobalV::ofs_running, ucell, "IMP_SOL FORCE", fsol, false); // this->print_force("IMP_SOL FORCE",fsol,1,ry); } - if (vdw_solver != nullptr) + if (vdw_result != nullptr) { ModuleIO::print_force(GlobalV::ofs_running, ucell, "VDW FORCE", force_vdw, false); // this->print_force("VDW FORCE",force_vdw,1,ry); @@ -746,7 +756,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, + sigmahar(i, j); // hartree stress // VDW stress from linpz and jiyy - if (vdw_solver != nullptr) + if (vdw_result != nullptr) { scs(i, j) += stress_vdw(i, j); } @@ -816,7 +826,7 @@ void Force_Stress_LCAO::getForceStress(UnitCell& ucell, ModuleIO::print_stress("EWALD STRESS", sigmaewa, screen, ry, GlobalV::ofs_running); ModuleIO::print_stress("cc STRESS", sigmacc, screen, ry, GlobalV::ofs_running); ModuleIO::print_stress("XC STRESS", sigmaxc, screen, ry, GlobalV::ofs_running); - if (vdw_solver != nullptr) + if (vdw_result != nullptr) { ModuleIO::print_stress("VDW STRESS", stress_vdw, screen, ry, GlobalV::ofs_running); } diff --git a/source/source_lcao/FORCE_STRESS.h b/source/source_lcao/FORCE_STRESS.h index 86d6f57b4c..a6e1caa913 100644 --- a/source/source_lcao/FORCE_STRESS.h +++ b/source/source_lcao/FORCE_STRESS.h @@ -18,6 +18,11 @@ #include "source_lcao/setup_dm.h" // mohan add 2025-11-03 #include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 +namespace vdw +{ +struct VdwResult; +} + template class Force_Stress_LCAO @@ -32,6 +37,7 @@ class Force_Stress_LCAO ~Force_Stress_LCAO(); void getForceStress(UnitCell& ucell, + const vdw::VdwResult* vdw_result, const bool isforce, const bool isstress, const bool istestf, diff --git a/source/source_pw/module_ofdft/of_stress_pw.cpp b/source/source_pw/module_ofdft/of_stress_pw.cpp index a481057507..6b988baf78 100644 --- a/source/source_pw/module_ofdft/of_stress_pw.cpp +++ b/source/source_pw/module_ofdft/of_stress_pw.cpp @@ -1,6 +1,7 @@ #include "of_stress_pw.h" #include "source_base/timer.h" +#include "source_base/tool_quit.h" #include "source_hamilt/module_vdw/vdw.h" #include "source_io/module_output/output_log.h" @@ -9,6 +10,7 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, ModuleBase::matrix& kinetic_stress, UnitCell& ucell, + const vdw::VdwResult* vdw_result, ModuleSymmetry::Symmetry* p_symm, const pseudopot_cell_vl& locpp, Structure_Factor* p_sf, @@ -80,8 +82,15 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, // nlcc stress_cc(sigmaxcc, this->rhopw, ucell, p_sf, true, locpp.numeric, pelec->charge); - // vdw term - stress_vdw(sigmavdw, ucell); + // vdW term prepared before SCF for this ionic configuration. + if (vdw_result != nullptr) + { + if (!vdw_result->has_stress) + { + ModuleBase::WARNING_QUIT("OF_Stress_PW::cal_stress", "The cached vdW stress is unavailable."); + } + sigmavdw = vdw_result->stress.to_matrix(); + } for (int ipol = 0; ipol < 3; ipol++) { @@ -119,13 +128,3 @@ void OF_Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, ModuleBase::timer::end("OF_Stress_PW", "cal_stress"); return; } - -void OF_Stress_PW::stress_vdw(ModuleBase::matrix& sigma, UnitCell& ucell) -{ - auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp); - if (vdw_solver != nullptr) - { - sigma = vdw_solver->get_stress().to_matrix(); - } - return; -} diff --git a/source/source_pw/module_ofdft/of_stress_pw.h b/source/source_pw/module_ofdft/of_stress_pw.h index d5d5d5feb1..17fc50f477 100644 --- a/source/source_pw/module_ofdft/of_stress_pw.h +++ b/source/source_pw/module_ofdft/of_stress_pw.h @@ -5,6 +5,11 @@ #include "source_pw/module_pwdft/vl_pw.h" #include "source_pw/module_pwdft/stress_func.h" +namespace vdw +{ +struct VdwResult; +} + class OF_Stress_PW : public Stress_Func { public: @@ -15,16 +20,13 @@ class OF_Stress_PW : public Stress_Func void cal_stress(ModuleBase::matrix& sigmatot, ModuleBase::matrix& kinetic_stress, UnitCell& ucell, + const vdw::VdwResult* vdw_result, ModuleSymmetry::Symmetry* p_symm, const pseudopot_cell_vl& locpp, Structure_Factor* p_sf, K_Vectors* p_kv); protected: - // call the vdw stress - void stress_vdw(ModuleBase::matrix& smearing_sigma, - UnitCell& ucell); // force and stress calculated in vdw together. - const elecstate::ElecState* pelec = nullptr; ModulePW::PW_Basis* rhopw = nullptr; }; diff --git a/source/source_pw/module_pwdft/forces.cpp b/source/source_pw/module_pwdft/forces.cpp index 386785cf9c..9330dac3f0 100644 --- a/source/source_pw/module_pwdft/forces.cpp +++ b/source/source_pw/module_pwdft/forces.cpp @@ -13,6 +13,7 @@ #include "source_base/mathzone.h" #include "source_base/timer.h" #include "source_base/tool_threading.h" +#include "source_base/tool_quit.h" #include "source_estate/module_pot/efield.h" #include "source_estate/module_pot/gatefield.h" #include "source_hamilt/module_ewald/H_Ewald_pw.h" @@ -28,6 +29,7 @@ template void Forces::cal_force(UnitCell& ucell, ModuleBase::matrix& force, + const vdw::VdwResult* vdw_result, const elecstate::ElecState& elec, const ModulePW::PW_Basis* const rho_basis, ModuleSymmetry::Symmetry* p_symm, @@ -87,18 +89,19 @@ void Forces::cal_force(UnitCell& ucell, // force due to core charge this->cal_force_scc(forcescc, rho_basis, elec.vnew, elec.vnew_exist, locpp->numeric, ucell); - ModuleBase::matrix stress_vdw_pw; //.create(3,3); ModuleBase::matrix force_vdw; force_vdw.create(nat, 3); - auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp); - if (vdw_solver != nullptr) + if (vdw_result != nullptr) { - const std::vector>& force_vdw_temp = vdw_solver->get_force(); + if (!vdw_result->has_force || vdw_result->force.size() != static_cast(this->nat)) + { + ModuleBase::WARNING_QUIT("Forces::cal_force", "The cached vdW force is unavailable or has an invalid size."); + } for (int iat = 0; iat < this->nat; ++iat) { - force_vdw(iat, 0) = force_vdw_temp[iat].x; - force_vdw(iat, 1) = force_vdw_temp[iat].y; - force_vdw(iat, 2) = force_vdw_temp[iat].z; + force_vdw(iat, 0) = vdw_result->force[iat].x; + force_vdw(iat, 1) = vdw_result->force[iat].y; + force_vdw(iat, 2) = vdw_result->force[iat].z; } if (PARAM.inp.test_force) { @@ -153,7 +156,7 @@ void Forces::cal_force(UnitCell& ucell, force(iat, ipol) = forcelc(iat, ipol) + forceion(iat, ipol) + forcenl(iat, ipol) + forcecc(iat, ipol) + forcescc(iat, ipol); - if (vdw_solver != nullptr) // linpz and jiyy added vdw force, modified by zhengdy + if (vdw_result != nullptr) // linpz and jiyy added vdw force, modified by zhengdy { force(iat, ipol) += force_vdw(iat, ipol); } diff --git a/source/source_pw/module_pwdft/forces.h b/source/source_pw/module_pwdft/forces.h index 3a229fc396..5a8531fdda 100644 --- a/source/source_pw/module_pwdft/forces.h +++ b/source/source_pw/module_pwdft/forces.h @@ -16,6 +16,11 @@ class pseudopot_cell_vnl; +namespace vdw +{ +struct VdwResult; +} + // forward declaration so that the dH module (out_mat_dh_vl) can reuse cal_force_loc namespace hamilt { template class Veff; } @@ -42,6 +47,7 @@ class Forces void cal_force(UnitCell& ucell, ModuleBase::matrix& force, + const vdw::VdwResult* vdw_result, const elecstate::ElecState& elec, const ModulePW::PW_Basis* const rho_basis, ModuleSymmetry::Symmetry* p_symm, diff --git a/source/source_pw/module_pwdft/stress_pw.cpp b/source/source_pw/module_pwdft/stress_pw.cpp index 4b26b7ec9f..dae145c5c4 100644 --- a/source/source_pw/module_pwdft/stress_pw.cpp +++ b/source/source_pw/module_pwdft/stress_pw.cpp @@ -1,6 +1,7 @@ #include "stress_pw.h" #include "source_base/timer.h" +#include "source_base/tool_quit.h" #include "source_base/global_variable.h" // use GlobalC #include "source_hamilt/module_vdw/vdw.h" #include "source_io/module_output/output_log.h" @@ -10,6 +11,7 @@ template void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, UnitCell& ucell, + const vdw::VdwResult* vdw_result, Plus_U &dftu, // mhan add 2025-11-07 const pseudopot_cell_vl& locpp, const pseudopot_cell_vnl& nlpp, @@ -116,8 +118,15 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, this->stress_us(sigmanl, rho_basis, nlpp, ucell); } - // vdw term - stress_vdw(sigmavdw, ucell); + // vdW term prepared before SCF for this ionic configuration. + if (vdw_result != nullptr) + { + if (!vdw_result->has_stress) + { + ModuleBase::WARNING_QUIT("Stress_PW::cal_stress", "The cached vdW stress is unavailable."); + } + sigmavdw = vdw_result->stress.to_matrix(); + } // DFT+U and DeltaSpin stress if (PARAM.inp.dft_plus_u || PARAM.inp.sc_mag_switch) @@ -182,16 +191,6 @@ void Stress_PW::cal_stress(ModuleBase::matrix& sigmatot, return; } -template -void Stress_PW::stress_vdw(ModuleBase::matrix& sigma, UnitCell& ucell) -{ - auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp); - if (vdw_solver != nullptr) - { - sigma = vdw_solver->get_stress().to_matrix(); - } - return; -} template class Stress_PW; #if ((defined __CUDA) || (defined __ROCM)) diff --git a/source/source_pw/module_pwdft/stress_pw.h b/source/source_pw/module_pwdft/stress_pw.h index 3defa09128..d77665d916 100644 --- a/source/source_pw/module_pwdft/stress_pw.h +++ b/source/source_pw/module_pwdft/stress_pw.h @@ -7,6 +7,11 @@ #include "source_lcao/module_dftu/dftu.h" // mohan add 2025-11-07 #include "source_lcao/module_ri/conv_coulomb_pot_k.h" +namespace vdw +{ +struct VdwResult; +} + template class Stress_PW : public Stress_Func { @@ -16,6 +21,7 @@ class Stress_PW : public Stress_Func // calculate the stress in PW basis void cal_stress(ModuleBase::matrix& smearing_sigmatot, UnitCell& ucell, + const vdw::VdwResult* vdw_result, Plus_U &dftu, // mhan add 2025-11-07 const pseudopot_cell_vl& locpp, const pseudopot_cell_vnl& nlpp, @@ -27,10 +33,6 @@ class Stress_PW : public Stress_Func const psi::Psi , Device>* d_psi_in = nullptr); protected: - // call the vdw stress - void stress_vdw(ModuleBase::matrix& smearing_sigma, - UnitCell& ucell); // force and stress calculated in vdw together. - // the stress from the non-local pseudopotentials in uspp // which is due to the dependence of the Q function on the atomic position void stress_us(ModuleBase::matrix& sigma,