Evaluate vdW corrections once per ionic step

Author: Growl1234

source/source_esolver/esolver_double_xc.cpp

@@ -1,7 +1,6 @@
 #include "esolver_double_xc.h"
 #include "source_hamilt/module_xc/xc_functional.h"
 #include "source_hamilt/module_ewald/H_Ewald_pw.h"
-#include "source_hamilt/module_vdw/vdw.h"
 #ifdef __MLALGO
 #include "source_lcao/module_deepks/LCAO_deepks.h"
 #include "source_lcao/module_deepks/LCAO_deepks_interface.h"
@@ -117,13 +116,9 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     ESolver_KS_LCAO<TK,TR>::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
@@ -382,6 +377,7 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
     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,

source/source_esolver/esolver_fp.cpp

@@ -187,11 +187,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

source/source_esolver/esolver_fp.h

@@ -14,6 +14,12 @@
 #include "source_pw/module_pwdft/structure_factor.h" // structure factor
 
 #include <fstream>
+#include <memory>
+
+namespace vdw
+{
+struct VdwResult;
+}
 
 
 //! The First-Principles (FP) Energy Solver Class
@@ -45,6 +51,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<const vdw::VdwResult> vdw_result_;
+
     //! These pointers will be deleted in the free_pointers() function every ion step.
     elecstate::ElecState* pelec = nullptr; ///< Electronic states
 

source/source_esolver/esolver_ks_lcao.cpp

@@ -242,7 +242,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
 
     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,

source/source_esolver/esolver_ks_pw.cpp

@@ -311,7 +311,7 @@ void ESolver_KS_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
     this->stp.update_psi_d();
 
     // Calculate forces
-    ff.cal_force(ucell, force, *this->pelec, this->pw_rhod, &ucell.symm,
+    ff.cal_force(ucell, force, this->get_vdw_result(), *this->pelec, this->pw_rhod, &ucell.symm,
                  &this->sf, this->solvent, &this->dftu, &this->locpp, &this->ppcell, 
                  &this->kv, this->pw_wfc, this->stp.template get_psi_d<T, Device>());
 }
@@ -324,7 +324,7 @@ void ESolver_KS_PW<T, Device>::cal_stress(UnitCell& ucell, ModuleBase::matrix& s
     // mohan add 2025-10-12
     this->stp.update_psi_d();
 
-    ss.cal_stress(stress, ucell, this->dftu, this->locpp, this->ppcell, this->pw_rhod,
+    ss.cal_stress(stress, ucell, this->get_vdw_result(), this->dftu, this->locpp, this->ppcell, this->pw_rhod,
                   &ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->stp.template get_psi_d<T, Device>());
 
     // external stress

source/source_esolver/esolver_of.cpp

@@ -546,7 +546,8 @@ void ESolver_OF::cal_force(UnitCell& ucell, 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);
 }
 
 /**
@@ -562,6 +563,6 @@ void ESolver_OF::cal_stress(UnitCell& ucell, 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

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<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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,15 +486,20 @@ 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);
 }
 
 TEST_F(vdwd3Test, D3bjGetForce)
 {
     input.vdw_method = "d3_bj";
     auto vdw_solver = vdw::make_vdw(ucell, input);
-    std::vector<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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,15 +609,20 @@ 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);
 }
 
 TEST_F(vdwd3abcTest, D3bjGetForce)
 {
     input.vdw_method = "d3_bj";
     auto vdw_solver = vdw::make_vdw(ucell, input);
-    std::vector<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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,15 +742,20 @@ 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);
 }
 
 TEST_F(vdwd4Test, D4SGetForce)
 {
     input.vdw_d4_model = "d4s";
     auto vdw_solver = vdw::make_vdw(ucell, input);
-    std::vector<ModuleBase::Vector3<double>> 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<ModuleBase::Vector3<double>>& 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);