Fix DFT-D4 calculations for charged systems

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);

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<int>& 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");

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;