Fix: Add the correct RT-TDDFT electric field force acting on ions

Author: AsTonyshment

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -18,6 +18,9 @@ double H_TDDFT_pw::amp;
 double H_TDDFT_pw::bmod;
 double H_TDDFT_pw::bvec[3];
 
+// Used for calculating electric field force on ions
+vector<double> H_TDDFT_pw::global_vext_time = {0.0, 0.0, 0.0};
+
 int H_TDDFT_pw::stype; // 0 : length gauge  1: velocity gauge
 
 std::vector<int> H_TDDFT_pw::ttype;
@@ -121,11 +124,22 @@ void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
     trigo_count = 0;
     heavi_count = 0;
 
+    global_vext_time = {0.0, 0.0, 0.0};
+
+    if (PARAM.inp.td_vext_dire.size() != 1)
+    {
+        ModuleBase::WARNING_QUIT("H_TDDFT_pw::cal_fixed_v",
+                                 "For the sake of program stability, the feature of applying multiple electric fields "
+                                 "simultaneously has been temporarily disabled. Thank you for your understanding!");
+    }
+
     for (auto direc: PARAM.inp.td_vext_dire)
     {
         std::vector<double> vext_space(this->rho_basis_->nrxx, 0.0);
         double vext_time = cal_v_time(ttype[count], true);
 
+        global_vext_time[direc - 1] += vext_time;
+
         if (PARAM.inp.out_efield && GlobalV::MY_RANK == 0)
         {
             std::stringstream as;
@@ -479,7 +493,7 @@ void H_TDDFT_pw::prepare(const ModuleBase::Matrix3& G, int& dir)
     bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2));
 }
 
-void H_TDDFT_pw ::compute_force(const UnitCell& cell, ModuleBase::matrix& fe)
+void H_TDDFT_pw::compute_force(const UnitCell& cell, ModuleBase::matrix& fe)
 {
     int iat = 0;
     for (int it = 0; it < cell.ntype; ++it)
@@ -488,7 +502,9 @@ void H_TDDFT_pw ::compute_force(const UnitCell& cell, ModuleBase::matrix& fe)
         {
             for (int jj = 0; jj < 3; ++jj)
             {
-                fe(iat, jj) = ModuleBase::e2 * amp * cell.atoms[it].ncpp.zv * bvec[jj] / bmod;
+                // No need to multiply ModuleBase::e2, since the unit of force is Ry/Bohr
+                fe(iat, jj)
+                    = (std::abs(bmod) > 1e-10 ? global_vext_time[jj] * cell.atoms[it].ncpp.zv * bvec[jj] / bmod : 0);
             }
             ++iat;
         }

source/module_elecstate/potentials/H_TDDFT_pw.h

@@ -124,6 +124,7 @@ class H_TDDFT_pw : public PotBase
     static bool is_initialized; // static flag variable, used to ensure initialization only once
 
     static double amp;
+    static vector<double> global_vext_time;
 
     static double bmod;
     static double bvec[3];

source/module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.cpp

@@ -8,6 +8,7 @@
 #include "module_base/timer.h"
 #include "module_cell/module_neighbor/sltk_grid_driver.h"
 #include "module_elecstate/elecstate_lcao.h"
+#include "module_elecstate/potentials/H_TDDFT_pw.h"       // Taoni add 2025-02-20
 #include "module_elecstate/potentials/efield.h"           // liuyu add 2022-05-18
 #include "module_elecstate/potentials/gatefield.h"        // liuyu add 2022-09-13
 #include "module_hamilt_general/module_surchem/surchem.h" //sunml add 2022-08-10
@@ -271,10 +272,10 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
 
     //! atomic forces from E-field of rt-TDDFT
     ModuleBase::matrix fefield_tddft;
-    if (PARAM.inp.esolver_type == "TDDFT" && isforce)
+    if (PARAM.inp.esolver_type == "tddft" && isforce)
     {
         fefield_tddft.create(nat, 3);
-        elecstate::Efield::compute_force(ucell, fefield_tddft);
+        elecstate::H_TDDFT_pw::compute_force(ucell, fefield_tddft);
     }
 
     //! atomic forces from gate field
@@ -447,7 +448,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                     fcs(iat, i) += fefield(iat, i);
                 }
                 // E-field force of tddft
-                if (PARAM.inp.esolver_type == "TDDFT")
+                if (PARAM.inp.esolver_type == "tddft")
                 {
                     fcs(iat, i) += fefield_tddft(iat, i);
                 }
@@ -566,7 +567,7 @@ void Force_Stress_LCAO<T>::getForceStress(UnitCell& ucell,
                 ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD     FORCE", fefield, false);
                 // this->print_force("EFIELD     FORCE",fefield,1,ry);
             }
-            if (PARAM.inp.esolver_type == "TDDFT")
+            if (PARAM.inp.esolver_type == "tddft")
             {
                 ModuleIO::print_force(GlobalV::ofs_running, ucell, "EFIELD_TDDFT     FORCE", fefield_tddft, false);
                 // this->print_force("EFIELD_TDDFT     FORCE",fefield_tddft,1,ry);