Fix: Incorrect negative signs in length gague RT-TDDFT

source/module_elecstate/potentials/H_TDDFT_pw.cpp

@@ -1,6 +1,7 @@
 #include "H_TDDFT_pw.h"
-#include "module_base/math_integral.h"
+
 #include "module_base/constants.h"
+#include "module_base/math_integral.h"
 #include "module_base/timer.h"
 #include "module_hamilt_lcao/module_tddft/evolve_elec.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
@@ -27,7 +28,7 @@ std::vector<int> H_TDDFT_pw::ttype;
 int H_TDDFT_pw::tstart;
 int H_TDDFT_pw::tend;
 double H_TDDFT_pw::dt;
-//cut dt for integral
+// cut dt for integral
 double H_TDDFT_pw::dt_int;
 int H_TDDFT_pw::istep_int;
 // space domain parameters
@@ -36,7 +37,7 @@ int H_TDDFT_pw::istep_int;
 double H_TDDFT_pw::lcut1;
 double H_TDDFT_pw::lcut2;
 
-//velocity gauge
+// velocity gauge
 ModuleBase::Vector3<double> H_TDDFT_pw::At;
 
 // time domain parameters
@@ -48,7 +49,7 @@ std::vector<double> H_TDDFT_pw::gauss_phase;
 std::vector<double> H_TDDFT_pw::gauss_sigma; // time(a.u.)
 std::vector<double> H_TDDFT_pw::gauss_t0;
 std::vector<double> H_TDDFT_pw::gauss_amp; // Ry/bohr
-std::vector<int> H_TDDFT_pw::gauss_ncut; //cut for integral
+std::vector<int> H_TDDFT_pw::gauss_ncut;   // cut for integral
 
 // trapezoid
 int H_TDDFT_pw::trape_count;
@@ -58,7 +59,7 @@ std::vector<double> H_TDDFT_pw::trape_t1;
 std::vector<double> H_TDDFT_pw::trape_t2;
 std::vector<double> H_TDDFT_pw::trape_t3;
 std::vector<double> H_TDDFT_pw::trape_amp; // Ry/bohr
-std::vector<int> H_TDDFT_pw::trape_ncut; //cut for integral
+std::vector<int> H_TDDFT_pw::trape_ncut;   // cut for integral
 
 // Trigonometric
 int H_TDDFT_pw::trigo_count;
@@ -67,20 +68,22 @@ std::vector<double> H_TDDFT_pw::trigo_omega2; // time(a.u.)^-1
 std::vector<double> H_TDDFT_pw::trigo_phase1;
 std::vector<double> H_TDDFT_pw::trigo_phase2;
 std::vector<double> H_TDDFT_pw::trigo_amp; // Ry/bohr
-std::vector<int> H_TDDFT_pw::trigo_ncut; //cut for integral
+std::vector<int> H_TDDFT_pw::trigo_ncut;   // cut for integral
 
 // Heaviside
 int H_TDDFT_pw::heavi_count;
 std::vector<double> H_TDDFT_pw::heavi_t0;
 std::vector<double> H_TDDFT_pw::heavi_amp; // Ry/bohr
 
-void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
+void H_TDDFT_pw::cal_fixed_v(double* vl_pseudo)
 {
     ModuleBase::TITLE("H_TDDFT_pw", "cal_fixed_v");
 
-    //skip if velocity_gague
-    if(stype==1) {return;
-}
+    // skip if velocity_gague
+    if (stype == 1)
+    {
+        return;
+    }
 
     // time evolve
     H_TDDFT_pw::istep++;
@@ -104,7 +107,7 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
     for (auto direc: module_tddft::Evolve_elec::td_vext_dire_case)
     {
         std::vector<double> vext_space(this->rho_basis_->nrxx, 0.0);
-        double vext_time = cal_v_time(ttype[count],true);
+        double vext_time = cal_v_time(ttype[count], true);
 
         if (module_tddft::Evolve_elec::out_efield && GlobalV::MY_RANK == 0)
         {
@@ -118,15 +121,17 @@ void H_TDDFT_pw::cal_fixed_v(double *vl_pseudo)
 
         cal_v_space(vext_space, direc);
         for (size_t ir = 0; ir < this->rho_basis_->nrxx; ++ir)
+        {
             vl_pseudo[ir] += vext_space[ir] * vext_time;
+        }
         count++;
     }
 
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_fixed_v");
     return;
 }
 
-void H_TDDFT_pw::cal_v_space(std::vector<double> &vext_space, int direc)
+void H_TDDFT_pw::cal_v_space(std::vector<double>& vext_space, int direc)
 {
     ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space");
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space");
@@ -145,7 +150,7 @@ void H_TDDFT_pw::cal_v_space(std::vector<double> &vext_space, int direc)
     return;
 }
 
-void H_TDDFT_pw::cal_v_space_length(std::vector<double> &vext_space, int direc)
+void H_TDDFT_pw::cal_v_space_length(std::vector<double>& vext_space, int direc)
 {
     ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length");
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length");
@@ -190,21 +195,22 @@ double H_TDDFT_pw::cal_v_space_length_potential(double i)
     double vext_space = 0.0;
     if (i < lcut1)
     {
-        vext_space = ((i - lcut1) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0;
+        vext_space = -((i - lcut1) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut1) * this->ucell_->lat0;
     }
     else if (i >= lcut1 && i < lcut2)
     {
-        vext_space = -i * this->ucell_->lat0;
+        vext_space = i * this->ucell_->lat0;
     }
     else if (i >= lcut2)
     {
-        vext_space = ((i - lcut2) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0;
+        vext_space = -((i - lcut2) * (lcut2 - lcut1) / (lcut1 + 1.0 - lcut2) - lcut2) * this->ucell_->lat0;
     }
     return vext_space;
 }
+
 int H_TDDFT_pw::check_ncut(int t_type)
 {
-    int ncut=0;
+    int ncut = 0;
     switch (t_type)
     {
     case 0:
@@ -233,12 +239,13 @@ int H_TDDFT_pw::check_ncut(int t_type)
     }
     return ncut;
 }
+
 void H_TDDFT_pw::update_At()
 {
     std::cout << "calculate electric potential" << std::endl;
     // time evolve
     H_TDDFT_pw::istep++;
-    
+
     // judgement to skip vext
     if (!module_tddft::Evolve_elec::td_vext || istep > tend || istep < tstart)
     {
@@ -249,41 +256,46 @@ void H_TDDFT_pw::update_At()
     trape_count = 0;
     trigo_count = 0;
     heavi_count = 0;
-    //parameters for integral
+    // parameters for integral
     int ncut = 1;
     bool last = false;
     double out = 0.0;
 
     for (auto direc: module_tddft::Evolve_elec::td_vext_dire_case)
     {
         last = false;
-        //cut the integral space and initial relvant parameters
+        // cut the integral space and initialize relevant parameters
         ncut = check_ncut(ttype[count]);
-        istep_int = istep*ncut;
-        dt_int = dt/double(ncut);
+        istep_int = istep * ncut;
+        dt_int = dt / double(ncut);
 
-        //store vext_time for each time point, include the first and last point
-        double* vext_time = new double[ncut+1]();
-        for(int i=0; i<=ncut; i++)
+        // store vext_time for each time point, include the first and last point
+        std::vector<double> vext_time(ncut + 1, 0.0); // Use std::vector to manage memory
+        for (int i = 0; i <= ncut; i++)
         {
-            //if this is the last point, type_count++
-            if(i==ncut)last=true;
-            vext_time[i] = cal_v_time(ttype[count],last);
+            // if this is the last point, type_count++
+            if (i == ncut)
+            {
+                last = true;
+            }
+            vext_time[i] = cal_v_time(ttype[count], last);
             istep_int++;
         }
-        ModuleBase::Integral::Simpson_Integral(ncut+1, vext_time, dt_int, out);
-        //update At value for its direction
+        // Call the Simpson's rule integration using std::vector data
+        ModuleBase::Integral::Simpson_Integral(ncut + 1, vext_time.data(), dt_int, out);
+
+        // update At value for its direction
         switch (stype)
         {
         case 1:
-            At[direc-1] -= out;
+            At[direc - 1] -= out;
             break;
         default:
             std::cout << "space_domain_type of electric field is wrong" << std::endl;
             break;
         }
 
-        //output Efield
+        // output Efield
         if (module_tddft::Evolve_elec::out_efield && GlobalV::MY_RANK == 0)
         {
             std::stringstream as;
@@ -293,8 +305,7 @@ void H_TDDFT_pw::update_At()
                 << vext_time[0] * ModuleBase::Ry_to_eV / ModuleBase::BOHR_TO_A << std::endl;
             ofs.close();
         }
-        //total count++
-        delete[] vext_time;
+        // total count++
         count++;
     }
     return;
@@ -345,8 +356,10 @@ double H_TDDFT_pw::cal_v_time_Gauss(const bool last)
 
     double gauss_t = (istep_int - t0 * ncut) * dt_int;
     vext_time = cos(omega * gauss_t + phase) * exp(-gauss_t * gauss_t * 0.5 / (sigma * sigma)) * amp;
-    if(last) {gauss_count++;
-}
+    if (last)
+    {
+        gauss_count++;
+    }
 
     return vext_time;
 }
@@ -376,8 +389,10 @@ double H_TDDFT_pw::cal_v_time_trapezoid(const bool last)
     }
 
     vext_time = vext_time * amp * cos(omega * istep_int * dt_int + phase);
-    if(last) {trape_count++;
-}
+    if (last)
+    {
+        trape_count++;
+    }
 
     return vext_time;
 }
@@ -394,8 +409,10 @@ double H_TDDFT_pw::cal_v_time_trigonometric(const bool last)
     const double timenow = istep_int * dt_int;
 
     vext_time = amp * cos(omega1 * timenow + phase1) * sin(omega2 * timenow + phase2) * sin(omega2 * timenow + phase2);
-    if(last) {trigo_count++;
-}
+    if (last)
+    {
+        trigo_count++;
+    }
 
     return vext_time;
 }
@@ -405,11 +422,14 @@ double H_TDDFT_pw::cal_v_time_heaviside()
     double t0 = *(heavi_t0.begin() + heavi_count);
     double amp = *(heavi_amp.begin() + heavi_count);
     double vext_time = 0.0;
-    if (istep < t0) {
+    if (istep < t0)
+    {
         vext_time = amp;
-    } else if (istep >= t0) {
+    }
+    else if (istep >= t0)
+    {
         vext_time = 0.0;
-}
+    }
     heavi_count++;
 
     return vext_time;

source/module_io/write_dipole.cpp

@@ -37,21 +37,26 @@ void ModuleIO::write_dipole(const double* rho_save,
 
 #ifndef __MPI
     double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0;
+    double lat_factor_x = GlobalC::ucell.lat0 * 0.529177 / rhopw->nx;
+    double lat_factor_y = GlobalC::ucell.lat0 * 0.529177 / rhopw->ny;
+    double lat_factor_z = GlobalC::ucell.lat0 * 0.529177 / rhopw->nz;
+
     for (int k = 0; k < rhopw->nz; k++)
     {
         for (int j = 0; j < rhopw->ny; j++)
         {
             for (int i = 0; i < rhopw->nx; i++)
             {
-                dipole_elec_x += rho_save[i * rhopw->ny * rhopw->nz + j * rhopw->nz + k] * i
-                                 * GlobalC::ucell.lat0 * 0.529177 / rhopw->nx;
-                dipole_elec_y += rho_save[i * rhopw->ny * rhopw->nz + j * rhopw->nz + k] * j
-                                 * GlobalC::ucell.lat0 * 0.529177 / rhopw->ny;
-                dipole_elec_z += rho_save[i * rhopw->ny * rhopw->nz + j * rhopw->nz + k] * k
-                                 * GlobalC::ucell.lat0 * 0.529177 / rhopw->nz;
+                int index = i * rhopw->ny * rhopw->nz + j * rhopw->nz + k;
+                double rho_val = rho_save[index];
+
+                dipole_elec_x += rho_val * i * lat_factor_x;
+                dipole_elec_y += rho_val * j * lat_factor_y;
+                dipole_elec_z += rho_val * k * lat_factor_z;
             }
         }
     }
+
     dipole_elec_x *= GlobalC::ucell.omega / static_cast<double>(rhopw->nxyz);
     dipole_elec_y *= GlobalC::ucell.omega / static_cast<double>(rhopw->nxyz);
     dipole_elec_z *= GlobalC::ucell.omega / static_cast<double>(rhopw->nxyz);
@@ -103,7 +108,6 @@ void ModuleIO::write_dipole(const double* rho_save,
             double sum = 0;
             for (int ia = 0; ia < GlobalC::ucell.atoms[it].na; ++ia)
             {
-
                 sum += GlobalC::ucell.atoms[it].taud[ia][i];
             }
             dipole_ion[i] += sum * GlobalC::ucell.atoms[it].ncpp.zv;
@@ -138,7 +142,7 @@ void ModuleIO::write_dipole(const double* rho_save,
     return;
 }
 
-double ModuleIO::prepare(const UnitCell &cell, int &dir)
+double ModuleIO::prepare(const UnitCell& cell, int& dir)
 {
     double bvec[3] = {0.0};
     double bmod = 0.0;

tests/integrate/601_NO_TDDFT_H2_len_current/refcurrent_total.dat

@@ -1,4 +1,4 @@
 0 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
 1 0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
-2 4.2068003950947370e-04 -3.6201227751435700e-05 -5.7838914079840259e-05
-3 -1.7465510685220912e-04 -9.7934211131431987e-05 -1.5648656468267278e-04
+2 -4.2482488063038011e-04 3.6208093141272442e-05 5.7880862345902912e-05
+3 1.7506474069499989e-04 9.7936761289276890e-05 1.5650847511644919e-04

tests/integrate/601_NO_TDDFT_H2_len_current/result.ref

@@ -1,6 +1,6 @@
-etotref -30.91255703869593
-etotperatomref -15.4562785193
-totalforceref 0.459640
-totalstressref 0.945730
+etotref -30.91255706614114
+etotperatomref -15.4562785331
+totalforceref 0.459938
+totalstressref 0.946457
 CompareCurrent_pass 0
 totaltimeref 1.72

tests/integrate/601_NO_TDDFT_H2_len_gauss/result.ref

@@ -1,5 +1,5 @@
-etotref -31.36940219279552
-etotperatomref -15.6847010964
-totalforceref 0.574354
-totalstressref 5.532849
+etotref -30.45477750879012
+etotperatomref -15.2273887544
+totalforceref 0.574814
+totalstressref 4.287440
 totaltimeref 1.57

tests/integrate/601_NO_TDDFT_H2_len_gauss_dire/result.ref

@@ -1,5 +1,5 @@
-etotref -31.36928614218948
-etotperatomref -15.6846430711
-totalforceref 0.572958
-totalstressref 5.533528
+etotref -30.45466274119030
+etotperatomref -15.2273313706
+totalforceref 0.573436
+totalstressref 4.285714
 totaltimeref 1.55

tests/integrate/601_NO_TDDFT_H2_len_heavi/result.ref

@@ -1,5 +1,5 @@
-etotref -38.85688376698273
-etotperatomref -19.4284418835
-totalforceref 0.564906
-totalstressref 77.567005
+etotref -22.97718184915547
+etotperatomref -11.4885909246
+totalforceref 0.564908
+totalstressref 76.073961
 totaltimeref 1.56

tests/integrate/601_NO_TDDFT_H2_len_hhg/result.ref

@@ -1,5 +1,5 @@
-etotref -82.62629897647498
-etotperatomref -41.3131494882
-totalforceref 0.748884
-totalstressref 501.373588
+etotref 20.38731422853930
+etotperatomref 10.1936571143
+totalforceref 0.738256
+totalstressref 489.377111
 totaltimeref 1.56

tests/integrate/601_NO_TDDFT_H2_len_trape/result.ref

@@ -1,5 +1,5 @@
-etotref -30.92452072277197
-etotperatomref -15.4622603614
-totalforceref 0.574664
-totalstressref 1.256903
+etotref -30.89965242483882
+etotperatomref -15.4498262124
+totalforceref 0.574662
+totalstressref 1.034216
 totaltimeref 1.44