Update write_dipole.cpp

docs/advanced/output_files/index.rst

@@ -7,4 +7,5 @@ Detailed Introduction of the Output Files
 
     output-specification
     running_scf.log
-
+    output_dipole.md
+

docs/advanced/output_files/output_dipole.md

@@ -0,0 +1,67 @@
+# Outputting Dipole Moment
+
+ABACUS can output the dipole moment by adding the keyword `out_dipole` in the INPUT file:
+
+```
+out_dipole             1
+```
+
+## Supported Calculations
+
+This feature is available for all types of DFT calculations that use charge density:
+
+- **KSDFT** (Kohn-Sham DFT)
+  - Plane wave (PW) basis
+  - Linear combination of atomic orbitals (LCAO) basis
+- **SDFT** (Stochastic DFT)
+- **OFDFT** (Orbital-Free DFT)
+- **TDDFT** (Time-Dependent DFT)
+
+## Input Parameters
+
+| Parameter | Type | Description | Default |
+|-----------|------|-------------|---------|
+| `out_dipole` | Integer | Whether to output dipole moment | 0 |
+
+- `out_dipole = 0`: Disable dipole output
+- `out_dipole = 1`: Enable dipole output
+
+## Output Files
+
+When `out_dipole` is set to 1, ABACUS will generate files named `dipole_s${spin}.txt` for each spin channel in the output directory.
+
+For spin-polarized calculation (nspin=2):
+- `dipole_s1.txt`
+- `dipole_s2.txt`
+
+For non-spin-polarized calculation (nspin=1):
+- `dipole_s1.txt`
+
+## Output Format
+
+The dipole output file contains one line for each ionic/electronic step:
+
+```
+step_index dipole_x dipole_y dipole_z
+```
+
+- `step_index`: The current step number (starts from 1)
+- `dipole_x, dipole_y, dipole_z`: The x, y, z components of the dipole moment
+
+Example output:
+
+```
+1 -0.00123456 0.00234567 -0.00345678
+2 -0.00123457 0.00234568 -0.00345679
+...
+```
+
+## Additional Information
+
+During the calculation, the dipole moment is also printed in the `running_*.log` file, including:
+- Electronic dipole moment
+- Ionic dipole moment
+- Total dipole moment
+- Total dipole moment norm
+
+The dipole moment calculation includes both electronic and ionic contributions. The total dipole moment is the sum of electronic and ionic dipoles.

examples/18_md/01_lcao_gamma_Si8/INPUT

@@ -0,0 +1,35 @@
+INPUT_PARAMETERS
+#Parameters (1.General)
+suffix                 Si_nve
+calculation            md
+nbands                 20
+symmetry               0
+pseudo_dir             ../../../tests/PP_ORB
+orbital_dir            ../../../tests/PP_ORB
+
+#Parameters (2.Iteration)
+ecutwfc                30
+scf_thr                1e-5
+scf_nmax               100
+
+#Parameters (3.Basis)
+basis_type             lcao
+ks_solver              genelpa
+gamma_only             1
+
+#Parameters (4.Smearing)
+smearing_method        gaussian
+smearing_sigma         0.001
+
+#Parameters (5.Mixing)
+mixing_type            broyden
+mixing_beta            0.3
+chg_extrap             second-order
+
+#Parameters (6.MD)
+md_type                nve
+md_nstep               10
+md_dt                  1
+md_tfirst              300
+
+out_dipole             1

examples/22_rt-tddft/01_H2_length_gauge/INPUT

@@ -20,7 +20,7 @@ smearing_method      gauss
 
 #Parameters (5.MD Parameters)
 md_type              nve
-md_nstep             1000
+md_nstep             5
 md_dt                0.005
 md_tfirst            0
 

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -249,7 +249,7 @@ template <typename TR, typename Device>
 void ESolver_KS_LCAO_TDDFT<TR, Device>::print_step()
 {
     std::cout << " -------------------------------------------" << std::endl;
-    std::cout << " STEP OF ELECTRON EVOLVE : " << unsigned(totstep) << std::endl;
+    std::cout << " STEP OF ELECTRON EVOLVE : " << unsigned(totstep)+1 << std::endl;
     std::cout << " -------------------------------------------" << std::endl;
 }
 

source/source_io/module_ctrl/ctrl_output_fp.cpp

@@ -1,6 +1,7 @@
 #include "ctrl_output_fp.h" // use ctrl_output_fp()
 
 #include "../module_output/cube_io.h"                                  // use write_vdata_palgrid
+#include "../module_dipole/dipole_io.h"                                // use write_dipole
 #include "source_estate/module_charge/symmetry_rho.h" // use Symmetry_rho
 #include "source_hamilt/module_xc/xc_functional.h"    // use XC_Functional
 #include "source_io/module_chgpot/write_elecstat_pot.h"             // use write_elecstat_pot
@@ -198,6 +199,17 @@ void ctrl_output_fp(UnitCell& ucell,
     }
 #endif
 
+    // 8) write dipole moment
+    if (PARAM.inp.out_dipole == 1)
+    {
+        for (int is = 0; is < nspin; ++is)
+        {
+            std::stringstream ss_dipole;
+            ss_dipole << global_out_dir << "dipole_s" << is + 1 << ".txt";
+            ModuleIO::write_dipole(ucell, chr.rho_save[is], pw_rhod, istep, ss_dipole.str(), GlobalV::ofs_running);
+        }
+    }
+
     ModuleBase::timer::end("ModuleIO", "ctrl_output_fp");
 }
 

source/source_io/module_ctrl/ctrl_output_td.cpp

@@ -1,7 +1,6 @@
 #include "ctrl_output_td.h"
 
 #include "source_base/parallel_global.h"
-#include "source_io/module_dipole/dipole_io.h"
 #include "source_io/module_parameter/parameter.h"
 #include "source_io/module_current/td_current_io.h"
 
@@ -30,18 +29,7 @@ void ctrl_output_td(const UnitCell& ucell,
     ModuleBase::TITLE("ModuleIO", "ctrl_output_td");
 
 #ifdef __LCAO
-    // (1) Write dipole information
-    for (int is = 0; is < PARAM.inp.nspin; ++is)
-    {
-        if (PARAM.inp.out_dipole == 1)
-        {
-            std::stringstream ss_dipole;
-            ss_dipole << PARAM.globalv.global_out_dir << "dipole_s" << is + 1 << ".txt";
-            ModuleIO::write_dipole(ucell, rho_save[is], rhopw, is, istep, ss_dipole.str());
-        }
-    }
-
-    // (2) Write current information
+    // (1) Write current information
     const elecstate::ElecStateLCAO<std::complex<double>>* pelec_lcao
         = dynamic_cast<const elecstate::ElecStateLCAO<std::complex<double>>*>(pelec);
 
@@ -65,7 +53,6 @@ void ctrl_output_td(const UnitCell& ucell,
     {
         ModuleIO::write_current(ucell, grid, istep, psi, pelec, kv, pv, orb, td_p->r_calculator, p_hamilt->getSR(), p_hamilt->getHR(), exx_nao);
     }
-
     // (3) Output file for restart
     if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0
     {

source/source_io/module_dipole/dipole_io.h

@@ -2,17 +2,19 @@
 #define DIPOLE_IO_H
 
 #include "source_basis/module_pw/pw_basis.h"
+#include "source_cell/unitcell.h"
 
+#include <fstream>
 #include <string>
 
 namespace ModuleIO
 {
 void write_dipole(const UnitCell& ucell,
                   const double* rho_save,
                   const ModulePW::PW_Basis* rhopw,
-                  const int& is,
                   const int& istep,
                   const std::string& fn,
+                  std::ofstream& ofs_running,
                   const int& precision = 11);
 
 double prepare(const UnitCell& cell, int& dir);