diff --git a/source/source_hamilt/module_vdw/test/vdw_test.cpp b/source/source_hamilt/module_vdw/test/vdw_test.cpp index faf62ba0139..ee021d66f88 100644 --- a/source/source_hamilt/module_vdw/test/vdw_test.cpp +++ b/source/source_hamilt/module_vdw/test/vdw_test.cpp @@ -623,7 +623,9 @@ class vdwd4Test: public testing::Test {0.3, 0.25, 0.25} }}}}; construct_ucell(structure,ucell); + ucell.atoms[0].ncpp.zv = 4.0; + input.nelec = 8.0; input.vdw_method = "d4"; input.vdw_d4_xc = "pbe"; input.vdw_d4_model = "d4"; @@ -646,6 +648,15 @@ TEST_F(vdwd4Test, D4GetEnergy) EXPECT_NEAR(ene, -0.04998837990336073, 1E-10); } +TEST_F(vdwd4Test, D4GetEnergyForChargedSystem) +{ + input.nelec = 7.0; + + auto vdw_solver = vdw::make_vdw(ucell, input); + const double ene = vdw_solver->get_energy(); + EXPECT_NEAR(ene, -0.04359451765256733, 1E-10); +} + TEST_F(vdwd4Test, D4GetForce) { auto vdw_solver = vdw::make_vdw(ucell, input); diff --git a/source/source_hamilt/module_vdw/vdwd4.cpp b/source/source_hamilt/module_vdw/vdwd4.cpp index 5eca0c63d32..1b8c7d3b50d 100644 --- a/source/source_hamilt/module_vdw/vdwd4.cpp +++ b/source/source_hamilt/module_vdw/vdwd4.cpp @@ -86,6 +86,13 @@ Vdwd4::Vdwd4(const UnitCell& unit_in, const std::string& xc_name, const Input_pa cutoff_disp2_ = cutoff_to_bohr(input.vdw_cutoff_radius, input.vdw_radius_unit); cutoff_disp3_ = std::min(40.0, cutoff_disp2_); cutoff_cn_ = length_to_bohr(input.vdw_cn_thr, input.vdw_cn_thr_unit); + + double valence_charge = 0.0; + for (int it = 0; it < ucell_.ntype; ++it) + { + valence_charge += ucell_.atoms[it].ncpp.zv * ucell_.atoms[it].na; + } + total_charge_ = valence_charge - input.nelec; } void Vdwd4::build_structure(std::vector& numbers, @@ -163,7 +170,7 @@ void Vdwd4::compute(double& energy_ha, ucell_.nat, numbers.data(), positions.data(), - nullptr, + &total_charge_, lattice.data(), periodic.data()); check_dftd4_error(error, "dftd4_new_structure"); diff --git a/source/source_hamilt/module_vdw/vdwd4.h b/source/source_hamilt/module_vdw/vdwd4.h index 0cb46c2ca8c..2c580388b67 100644 --- a/source/source_hamilt/module_vdw/vdwd4.h +++ b/source/source_hamilt/module_vdw/vdwd4.h @@ -23,6 +23,7 @@ class Vdwd4 : public Vdw double cutoff_disp2_ = 0.0; // Bohr, two-body dispersion cutoff double cutoff_disp3_ = 0.0; // Bohr, three-body ATM cutoff 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;