From a40b37e957d18cd553832666349a3108a8d0aa13 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 16:52:36 +0800 Subject: [PATCH 01/11] change the output of step of electron evolve, starts from 1 --- examples/22_rt-tddft/01_H2_length_gauge/INPUT | 2 +- source/source_esolver/esolver_ks_lcao_tddft.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/22_rt-tddft/01_H2_length_gauge/INPUT b/examples/22_rt-tddft/01_H2_length_gauge/INPUT index 99d7386a82b..bd68e58862e 100644 --- a/examples/22_rt-tddft/01_H2_length_gauge/INPUT +++ b/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 diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp index 61e43797e13..033dbfa58d3 100644 --- a/source/source_esolver/esolver_ks_lcao_tddft.cpp +++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp @@ -249,7 +249,7 @@ template void ESolver_KS_LCAO_TDDFT::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; } From 8b80e5d5acf4f5acb4c9c68892b0f0b7ff60fcf9 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 17:32:02 +0800 Subject: [PATCH 02/11] update dipole outputs --- .../source_io/module_ctrl/ctrl_output_td.cpp | 2 +- source/source_io/module_dipole/dipole_io.h | 2 + .../source_io/module_dipole/write_dipole.cpp | 139 ++++++++++++------ 3 files changed, 93 insertions(+), 50 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_output_td.cpp b/source/source_io/module_ctrl/ctrl_output_td.cpp index 65d0c9ab965..e120285a8f9 100644 --- a/source/source_io/module_ctrl/ctrl_output_td.cpp +++ b/source/source_io/module_ctrl/ctrl_output_td.cpp @@ -37,7 +37,7 @@ void ctrl_output_td(const UnitCell& ucell, { 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()); + ModuleIO::write_dipole(ucell, rho_save[is], rhopw, is, istep, ss_dipole.str(), GlobalV::ofs_running); } } diff --git a/source/source_io/module_dipole/dipole_io.h b/source/source_io/module_dipole/dipole_io.h index b97810bf757..f0dc086f8a7 100644 --- a/source/source_io/module_dipole/dipole_io.h +++ b/source/source_io/module_dipole/dipole_io.h @@ -3,6 +3,7 @@ #include "source_basis/module_pw/pw_basis.h" +#include #include namespace ModuleIO @@ -13,6 +14,7 @@ void write_dipole(const UnitCell& ucell, 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); diff --git a/source/source_io/module_dipole/write_dipole.cpp b/source/source_io/module_dipole/write_dipole.cpp index 56448365bfe..28913513879 100644 --- a/source/source_io/module_dipole/write_dipole.cpp +++ b/source/source_io/module_dipole/write_dipole.cpp @@ -1,15 +1,30 @@ #include "source_base/parallel_reduce.h" -#include "source_estate/module_charge/charge.h" #include "source_io/module_dipole/dipole_io.h" -#include "source_lcao/module_rt/evolve_elec.h" -// fuxiang add 2017-03-15 +// Calculate and write dipole moment data for RT-TDDFT calculations +// +// Dipole moment is a measure of the separation of positive and negative charges +// in a system. In electronic structure calculations, we compute three components: +// +// 1. Electronic dipole moment (P_elec): +// Formula: P_elec = -int(r * rho(r) dr) +// where rho(r) is the electron density +// +// 2. Ionic dipole moment (P_ion): +// Formula: P_ion = sum_atoms (Z_v * r_atom) +// where Z_v is the valence charge of the atom +// +// 3. Total dipole moment (P_tot): +// Formula: P_tot = P_elec + P_ion +// +// The total dipole moment norm is |P_tot| = sqrt(P_tot_x^2 + P_tot_y^2 + P_tot_z^2) void ModuleIO::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) { ModuleBase::TITLE("ModuleIO", "write_dipole"); @@ -24,79 +39,93 @@ void ModuleIO::write_dipole(const UnitCell& ucell, ofs.open(fn.c_str(), std::ofstream::app); if (!ofs) { - ModuleBase::WARNING("ModuleIO", "Can't create Charge File!"); + ModuleBase::WARNING("ModuleIO", "Can't create Dipole File!"); } } + ofs_running << " Write dipole data to file: " << fn << std::endl; + + // Calculate modulus of reciprocal lattice vectors + // bmod[i] = |b_i| where b_i are reciprocal lattice vectors + // Used for coordinate transformation from fractional to Cartesian double bmod[3]; for (int i = 0; i < 3; i++) { bmod[i] = prepare(ucell, i); + if (bmod[i] < 1e-10) + { + ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "bmod[" + std::to_string(i) + "] is zero or too small!"); + } } -#ifndef __MPI - double dipole_elec_x = 0.0, dipole_elec_y = 0.0, dipole_elec_z = 0.0; - double lat_factor_x = ucell.lat0 * ModuleBase::BOHR_TO_A / rhopw->nx; - double lat_factor_y = ucell.lat0 * ModuleBase::BOHR_TO_A / rhopw->ny; - double lat_factor_z = ucell.lat0 * ModuleBase::BOHR_TO_A / rhopw->nz; - - for (int k = 0; k < rhopw->nz; k++) + // Validate grid dimensions to prevent division by zero + if (rhopw->nx == 0 || rhopw->ny == 0 || rhopw->nz == 0) { - for (int j = 0; j < rhopw->ny; j++) - { - for (int i = 0; i < rhopw->nx; i++) - { - 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; - } - } + ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "Grid dimensions nx, ny, or nz cannot be zero!"); } - dipole_elec_x *= ucell.omega / static_cast(rhopw->nxyz); - dipole_elec_y *= ucell.omega / static_cast(rhopw->nxyz); - dipole_elec_z *= ucell.omega / static_cast(rhopw->nxyz); - Parallel_Reduce::reduce_pool(dipole_elec_x); - Parallel_Reduce::reduce_pool(dipole_elec_y); - Parallel_Reduce::reduce_pool(dipole_elec_z); + if (rhopw->nxyz == 0) + { + ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "Total grid points nxyz cannot be zero!"); + } - ofs << istep+1 << " " << dipole_elec_x << " " << dipole_elec_y << dipole_elec_z; -#else + if (rhopw->nplane == 0) + { + ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "nplane cannot be zero!"); + } + // Calculate electronic dipole moment + // P_elec = -int(r * rho(r) dr) + // Discretized: P_elec[i] = -sum(r_grid[i] * rho(r_grid)) * (omega / nxyz) double dipole_elec[3] = {0.0, 0.0, 0.0}; + // Loop over local grid points (parallel decomposition) for (int ir = 0; ir < rhopw->nrxx; ++ir) { + // Convert 1D index to 3D indices int i = ir / (rhopw->ny * rhopw->nplane); int j = ir / rhopw->nplane - i * rhopw->ny; int k = ir % rhopw->nplane + rhopw->startz_current; + + // Convert to fractional coordinates: r_i = i / N_i double x = (double)i / rhopw->nx; double y = (double)j / rhopw->ny; double z = (double)k / rhopw->nz; + // Accumulate: P_elec -= rho * r (negative sign from electron charge) dipole_elec[0] -= rho_save[ir] * x; dipole_elec[1] -= rho_save[ir] * y; dipole_elec[2] -= rho_save[ir] * z; } + // Reduce across MPI processes to get global sum Parallel_Reduce::reduce_pool(dipole_elec[0]); Parallel_Reduce::reduce_pool(dipole_elec[1]); Parallel_Reduce::reduce_pool(dipole_elec[2]); + + // Convert to Cartesian coordinates and normalize + // Conversion factor: lat0 / bmod[i] transforms fractional to Cartesian + // Volume normalization: omega / nxyz accounts for grid spacing for (int i = 0; i < 3; ++i) { dipole_elec[i] *= ucell.lat0 / bmod[i] * ucell.omega / rhopw->nxyz; } - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Electronic dipole moment P_elec_x(t)", dipole_elec[0]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Electronic dipole moment P_elec_y(t)", dipole_elec[1]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Electronic dipole moment P_elec_z(t)", dipole_elec[2]); - - ofs << std::setprecision(precision) << istep+1 << " " << dipole_elec[0] << " " << dipole_elec[1] << " " - << dipole_elec[2] << std::endl; - + // Output electronic dipole moment + ofs_running << " Electronic dipole moment" << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running, "P_elec_x(t)", dipole_elec[0]); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_elec_y(t)", dipole_elec[1]); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_elec_z(t)", dipole_elec[2]); + + // Write to file: step index and three dipole components + ofs << std::setprecision(precision) << istep+1 + << " " << dipole_elec[0] + << " " << dipole_elec[1] + << " " << dipole_elec[2] << std::endl; + + // Calculate ionic dipole moment + // P_ion = sum_{atom_types} sum_{atoms} (Z_v * tau) + // where tau is the atomic position in fractional coordinates double dipole_ion[3] = {0.0}; double dipole_sum = 0.0; @@ -111,28 +140,36 @@ void ModuleIO::write_dipole(const UnitCell& ucell, } dipole_ion[i] += sum * ucell.atoms[it].ncpp.zv; } - dipole_ion[i] *= ucell.lat0 / bmod[i]; //* ModuleBase::FOUR_PI / ucell.omega; + // Convert to Cartesian coordinates + dipole_ion[i] *= ucell.lat0 / bmod[i]; } - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Ionic dipole moment P_ion_x(t)", dipole_ion[0]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Ionic dipole moment P_ion_y(t)", dipole_ion[1]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Ionic dipole moment P_ion_z(t)", dipole_ion[2]); + // Output ionic dipole moment + ofs_running << " Ionic dipole moment" << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running, "P_ion_x(t)", dipole_ion[0]); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_ion_y(t)", dipole_ion[1]); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_ion_z(t)", dipole_ion[2]); + // Calculate total dipole moment + // P_tot = P_elec + P_ion double dipole[3] = {0.0}; for (int i = 0; i < 3; ++i) { dipole[i] = dipole_ion[i] + dipole_elec[i]; } - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total dipole moment P_tot_x(t)", dipole[0]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total dipole moment P_tot_y(t)", dipole[1]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total dipole moment P_tot_z(t)", dipole[2]); + // Output total dipole moment + ofs_running << " Total dipole moment" << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running, "P_tot_x(t)", dipole[0]); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_tot_y(t)", dipole[1]); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_tot_z(t)", dipole[2]); + // Calculate and output total dipole moment norm + // |P_tot| = sqrt(P_tot_x^2 + P_tot_y^2 + P_tot_z^2) dipole_sum = sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]); - ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Total dipole moment norm |P_tot(t)|", dipole_sum); - -#endif + ofs_running << " Total dipole moment norm" << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running, "|P_tot(t)|", dipole_sum); if (GlobalV::MY_RANK == 0) { @@ -144,6 +181,10 @@ void ModuleIO::write_dipole(const UnitCell& ucell, return; } +// Calculate the modulus of a reciprocal lattice vector +// Input: cell - unit cell containing reciprocal lattice G +// dir - direction index (0=x, 1=y, 2=z) +// Output: bmod - |b_dir| where b_dir is the reciprocal lattice vector double ModuleIO::prepare(const UnitCell& cell, int& dir) { double bvec[3] = {0.0}; @@ -172,4 +213,4 @@ double ModuleIO::prepare(const UnitCell& cell, int& dir) } bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); return bmod; -} +} \ No newline at end of file From c6afecfea81f866210698745e667eb90e8e2f5c6 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 17:48:09 +0800 Subject: [PATCH 03/11] =?UTF-8?q?=E6=94=B9=E8=BF=9B=20dipole=5Fio=20?= =?UTF-8?q?=E6=A8=A1=E5=9D=97=EF=BC=9A=E4=BC=98=E5=8C=96=E6=80=A7=E8=83=BD?= =?UTF-8?q?=E5=92=8C=E4=BB=A3=E7=A0=81=E7=BB=93=E6=9E=84?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 1. 增强 write_dipole 函数： - 添加 ofs_running 参数，支持自定义输出流 - 移除未使用的头文件引用（charge.h, evolve_elec.h） - 统一电子偶极矩计算算法，移除 ifdef MPI 分支 - 预计算倒网格维度倒数，优化性能 - 添加除零保护检查 - 添加 OpenMP 并行化，使用 reduction 累加 - 提取公共打印函数 printDipoleMoment 2. 改进 prepare 函数： - 使用 switch 语句替代 if-else 链 - 移除中间变量，直接返回结果 3. 更新 dipole_io.h： - 添加 UnitCell 头文件引用 - 调整参数顺序，添加 ofs_running 在前 4. 添加详细注释： - 函数文档注释 - 物理公式说明 - 代码实现细节 5. 优化代码风格： - 使用小写常量名（small_value） - 优化变量作用域 - 改进错误处理 --- source/source_io/module_dipole/dipole_io.h | 1 + .../source_io/module_dipole/write_dipole.cpp | 189 ++++++++++-------- 2 files changed, 108 insertions(+), 82 deletions(-) diff --git a/source/source_io/module_dipole/dipole_io.h b/source/source_io/module_dipole/dipole_io.h index f0dc086f8a7..1a86f82c885 100644 --- a/source/source_io/module_dipole/dipole_io.h +++ b/source/source_io/module_dipole/dipole_io.h @@ -2,6 +2,7 @@ #define DIPOLE_IO_H #include "source_basis/module_pw/pw_basis.h" +#include "source_cell/unitcell.h" #include #include diff --git a/source/source_io/module_dipole/write_dipole.cpp b/source/source_io/module_dipole/write_dipole.cpp index 28913513879..f1ebb703c75 100644 --- a/source/source_io/module_dipole/write_dipole.cpp +++ b/source/source_io/module_dipole/write_dipole.cpp @@ -1,6 +1,21 @@ #include "source_base/parallel_reduce.h" #include "source_io/module_dipole/dipole_io.h" +namespace ModuleIO +{ + +// Helper function to print dipole moment components +// name: descriptive name for the dipole moment type +// px, py, pz: dipole moment components in x, y, z directions +void printDipoleMoment(std::ofstream& ofs_running, const std::string& name, + double px, double py, double pz) +{ + ofs_running << " " << name << std::endl; + ModuleBase::GlobalFunc::OUT(ofs_running, "P_x(t)", px); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_y(t)", py); + ModuleBase::GlobalFunc::OUT(ofs_running, "P_z(t)", pz); +} + // Calculate and write dipole moment data for RT-TDDFT calculations // // Dipole moment is a measure of the separation of positive and negative charges @@ -11,35 +26,45 @@ // where rho(r) is the electron density // // 2. Ionic dipole moment (P_ion): -// Formula: P_ion = sum_atoms (Z_v * r_atom) -// where Z_v is the valence charge of the atom +// Formula: P_ion = sum_{atom_types} sum_{atoms} (Z_v * tau) +// where Z_v is the valence charge and tau is atomic position // // 3. Total dipole moment (P_tot): // Formula: P_tot = P_elec + P_ion // // The total dipole moment norm is |P_tot| = sqrt(P_tot_x^2 + P_tot_y^2 + P_tot_z^2) -void ModuleIO::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) +// +// Parameters: +// - ucell: unit cell containing atomic structure and lattice information +// - rho_save: electron density on real-space grid +// - rhopw: plane wave basis information including grid dimensions +// - is: spin channel index +// - istep: current time step +// - fn: output file name +// - ofs_running: output stream for logging +// - precision: floating-point precision for output +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) { ModuleBase::TITLE("ModuleIO", "write_dipole"); time_t start, end; std::ofstream ofs; + // Open output file on master process only if (GlobalV::MY_RANK == 0) { start = time(NULL); - ofs.open(fn.c_str(), std::ofstream::app); if (!ofs) { - ModuleBase::WARNING("ModuleIO", "Can't create Dipole File!"); + ModuleBase::WARNING_QUIT("ModuleIO", "Cannot create dipole file: " + fn); } } @@ -48,38 +73,40 @@ void ModuleIO::write_dipole(const UnitCell& ucell, // Calculate modulus of reciprocal lattice vectors // bmod[i] = |b_i| where b_i are reciprocal lattice vectors // Used for coordinate transformation from fractional to Cartesian + double small_value = 1e-10; double bmod[3]; - for (int i = 0; i < 3; i++) + + for (int i = 0; i < 3; ++i) { bmod[i] = prepare(ucell, i); - if (bmod[i] < 1e-10) + if (bmod[i] < small_value) { - ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "bmod[" + std::to_string(i) + "] is zero or too small!"); + ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", + "bmod[" + std::to_string(i) + "] is too small or zero"); } } // Validate grid dimensions to prevent division by zero - if (rhopw->nx == 0 || rhopw->ny == 0 || rhopw->nz == 0) - { - ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "Grid dimensions nx, ny, or nz cannot be zero!"); - } - - if (rhopw->nxyz == 0) - { - ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "Total grid points nxyz cannot be zero!"); - } - - if (rhopw->nplane == 0) + if (rhopw->nx == 0 || rhopw->ny == 0 || rhopw->nz == 0 || + rhopw->nxyz == 0 || rhopw->nplane == 0) { - ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", "nplane cannot be zero!"); + ModuleBase::WARNING_QUIT("ModuleIO::write_dipole", + "Invalid grid parameters: nx, ny, nz, nxyz, or nplane is zero"); } // Calculate electronic dipole moment // P_elec = -int(r * rho(r) dr) // Discretized: P_elec[i] = -sum(r_grid[i] * rho(r_grid)) * (omega / nxyz) double dipole_elec[3] = {0.0, 0.0, 0.0}; - - // Loop over local grid points (parallel decomposition) + + // Precompute inverse grid dimensions for performance + double inv_nx = 1.0 / static_cast(rhopw->nx); + double inv_ny = 1.0 / static_cast(rhopw->ny); + double inv_nz = 1.0 / static_cast(rhopw->nz); + + // Loop over local grid points (parallel decomposition with OpenMP) + // Use reduction for thread-safe accumulation + #pragma omp parallel for reduction(-:dipole_elec[:3]) schedule(static) for (int ir = 0; ir < rhopw->nrxx; ++ir) { // Convert 1D index to 3D indices @@ -88,9 +115,10 @@ void ModuleIO::write_dipole(const UnitCell& ucell, int k = ir % rhopw->nplane + rhopw->startz_current; // Convert to fractional coordinates: r_i = i / N_i - double x = (double)i / rhopw->nx; - double y = (double)j / rhopw->ny; - double z = (double)k / rhopw->nz; + // Using multiplication instead of division for better performance + double x = static_cast(i) * inv_nx; + double y = static_cast(j) * inv_ny; + double z = static_cast(k) * inv_nz; // Accumulate: P_elec -= rho * r (negative sign from electron charge) dipole_elec[0] -= rho_save[ir] * x; @@ -106,29 +134,26 @@ void ModuleIO::write_dipole(const UnitCell& ucell, // Convert to Cartesian coordinates and normalize // Conversion factor: lat0 / bmod[i] transforms fractional to Cartesian // Volume normalization: omega / nxyz accounts for grid spacing + double vol_factor = ucell.omega / static_cast(rhopw->nxyz); for (int i = 0; i < 3; ++i) { - dipole_elec[i] *= ucell.lat0 / bmod[i] * ucell.omega / rhopw->nxyz; + dipole_elec[i] *= ucell.lat0 / bmod[i] * vol_factor; } // Output electronic dipole moment - ofs_running << " Electronic dipole moment" << std::endl; - ModuleBase::GlobalFunc::OUT(ofs_running, "P_elec_x(t)", dipole_elec[0]); - ModuleBase::GlobalFunc::OUT(ofs_running, "P_elec_y(t)", dipole_elec[1]); - ModuleBase::GlobalFunc::OUT(ofs_running, "P_elec_z(t)", dipole_elec[2]); + printDipoleMoment(ofs_running, "Electronic dipole moment", + dipole_elec[0], dipole_elec[1], dipole_elec[2]); // Write to file: step index and three dipole components - ofs << std::setprecision(precision) << istep+1 + ofs << std::setprecision(precision) << istep + 1 << " " << dipole_elec[0] - << " " << dipole_elec[1] - << " " << dipole_elec[2] << std::endl; + << " " << dipole_elec[1] + << " " << dipole_elec[2] << std::endl; // Calculate ionic dipole moment // P_ion = sum_{atom_types} sum_{atoms} (Z_v * tau) // where tau is the atomic position in fractional coordinates double dipole_ion[3] = {0.0}; - double dipole_sum = 0.0; - for (int i = 0; i < 3; ++i) { for (int it = 0; it < ucell.ntype; ++it) @@ -145,10 +170,8 @@ void ModuleIO::write_dipole(const UnitCell& ucell, } // Output ionic dipole moment - ofs_running << " Ionic dipole moment" << std::endl; - ModuleBase::GlobalFunc::OUT(ofs_running, "P_ion_x(t)", dipole_ion[0]); - ModuleBase::GlobalFunc::OUT(ofs_running, "P_ion_y(t)", dipole_ion[1]); - ModuleBase::GlobalFunc::OUT(ofs_running, "P_ion_z(t)", dipole_ion[2]); + printDipoleMoment(ofs_running, "Ionic dipole moment", + dipole_ion[0], dipole_ion[1], dipole_ion[2]); // Calculate total dipole moment // P_tot = P_elec + P_ion @@ -159,58 +182,60 @@ void ModuleIO::write_dipole(const UnitCell& ucell, } // Output total dipole moment - ofs_running << " Total dipole moment" << std::endl; - ModuleBase::GlobalFunc::OUT(ofs_running, "P_tot_x(t)", dipole[0]); - ModuleBase::GlobalFunc::OUT(ofs_running, "P_tot_y(t)", dipole[1]); - ModuleBase::GlobalFunc::OUT(ofs_running, "P_tot_z(t)", dipole[2]); + printDipoleMoment(ofs_running, "Total dipole moment", + dipole[0], dipole[1], dipole[2]); // Calculate and output total dipole moment norm // |P_tot| = sqrt(P_tot_x^2 + P_tot_y^2 + P_tot_z^2) - dipole_sum = sqrt(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2]); - + double dipole_sum = sqrt(dipole[0] * dipole[0] + + dipole[1] * dipole[1] + + dipole[2] * dipole[2]); ofs_running << " Total dipole moment norm" << std::endl; ModuleBase::GlobalFunc::OUT(ofs_running, "|P_tot(t)|", dipole_sum); + // Close file and report timing on master process if (GlobalV::MY_RANK == 0) { end = time(NULL); ModuleBase::GlobalFunc::OUT_TIME("write_dipole", start, end); ofs.close(); } - - return; } // Calculate the modulus of a reciprocal lattice vector -// Input: cell - unit cell containing reciprocal lattice G -// dir - direction index (0=x, 1=y, 2=z) -// Output: bmod - |b_dir| where b_dir is the reciprocal lattice vector -double ModuleIO::prepare(const UnitCell& cell, int& dir) +// Input: +// cell - unit cell containing reciprocal lattice G +// dir - direction index (0=x, 1=y, 2=z) +// Output: +// bmod - |b_dir| where b_dir is the reciprocal lattice vector +double prepare(const UnitCell& cell, int& dir) { double bvec[3] = {0.0}; - double bmod = 0.0; - if (dir == 0) - { - bvec[0] = cell.G.e11; - bvec[1] = cell.G.e12; - bvec[2] = cell.G.e13; - } - else if (dir == 1) + + // Select the appropriate reciprocal lattice vector components + switch (dir) { - bvec[0] = cell.G.e21; - bvec[1] = cell.G.e22; - bvec[2] = cell.G.e23; + case 0: // x-direction + bvec[0] = cell.G.e11; + bvec[1] = cell.G.e12; + bvec[2] = cell.G.e13; + break; + case 1: // y-direction + bvec[0] = cell.G.e21; + bvec[1] = cell.G.e22; + bvec[2] = cell.G.e23; + break; + case 2: // z-direction + bvec[0] = cell.G.e31; + bvec[1] = cell.G.e32; + bvec[2] = cell.G.e33; + break; + default: + ModuleBase::WARNING_QUIT("ModuleIO::prepare", "Invalid direction index"); } - else if (dir == 2) - { - bvec[0] = cell.G.e31; - bvec[1] = cell.G.e32; - bvec[2] = cell.G.e33; - } - else - { - ModuleBase::WARNING_QUIT("ModuleIO::prepare", "direction is wrong!"); - } - bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); - return bmod; -} \ No newline at end of file + + // Calculate and return the magnitude of the reciprocal lattice vector + return sqrt(bvec[0] * bvec[0] + bvec[1] * bvec[1] + bvec[2] * bvec[2]); +} + +} // namespace ModuleIO \ No newline at end of file From d25290f9616fda4839b3d16af16ca706e40b28b3 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:12:48 +0800 Subject: [PATCH 04/11] delete uselss #ifdef __LCAO --- source/source_io/module_ctrl/ctrl_output_td.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_output_td.cpp b/source/source_io/module_ctrl/ctrl_output_td.cpp index e120285a8f9..9588b67e417 100644 --- a/source/source_io/module_ctrl/ctrl_output_td.cpp +++ b/source/source_io/module_ctrl/ctrl_output_td.cpp @@ -29,7 +29,6 @@ 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) { @@ -82,7 +81,6 @@ void ctrl_output_td(const UnitCell& ucell, } } } -#endif // __LCAO } template void ctrl_output_td(const UnitCell&, From 057a4e0261ed9468dbd59ba85043e060bd9ee2dc Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:29:21 +0800 Subject: [PATCH 05/11] Refactor: migrate write_dipole from ctrl_output_td to ctrl_output_fp This refactoring extends the dipole output functionality to all DFT solvers by moving it from the TDDFT-specific module to the common output module. Changes: 1. source/source_io/module_ctrl/ctrl_output_fp.cpp - Add include for dipole_io.h - Add dipole output functionality after out_xc_r section - Now supports out_dipole parameter for all FP-based solvers 2. source/source_io/module_ctrl/ctrl_output_td.cpp - Remove duplicate dipole output code (migrated to ctrl_output_fp) - Remove unnecessary dipole_io.h include - Renumber comments (1) for current, (3) for restart 3. source/source_io/module_dipole/write_dipole.cpp - Rename printDipoleMoment to print_dipole_moment (snake_case convention) Benefits: - OFDFT, KSDFT (PW/LCAO), SDFT, and TDDFT all gain dipole output capability - Single implementation point reduces maintenance burden - Consistent behavior across all solver types - Better code organization following the base class aggregation pattern --- source/source_io/module_ctrl/ctrl_output_fp.cpp | 12 ++++++++++++ source/source_io/module_ctrl/ctrl_output_td.cpp | 14 +------------- source/source_io/module_dipole/write_dipole.cpp | 8 ++++---- 3 files changed, 17 insertions(+), 17 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 8e3991c8a14..13ba27ebc44 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/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, is, istep, ss_dipole.str(), GlobalV::ofs_running); + } + } + ModuleBase::timer::end("ModuleIO", "ctrl_output_fp"); } diff --git a/source/source_io/module_ctrl/ctrl_output_td.cpp b/source/source_io/module_ctrl/ctrl_output_td.cpp index 9588b67e417..6d7fbcc6658 100644 --- a/source/source_io/module_ctrl/ctrl_output_td.cpp +++ b/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" @@ -29,18 +28,7 @@ void ctrl_output_td(const UnitCell& ucell, { ModuleBase::TITLE("ModuleIO", "ctrl_output_td"); - // (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(), GlobalV::ofs_running); - } - } - - // (2) Write current information + // (1) Write current information const elecstate::ElecStateLCAO>* pelec_lcao = dynamic_cast>*>(pelec); diff --git a/source/source_io/module_dipole/write_dipole.cpp b/source/source_io/module_dipole/write_dipole.cpp index f1ebb703c75..f9343c27a91 100644 --- a/source/source_io/module_dipole/write_dipole.cpp +++ b/source/source_io/module_dipole/write_dipole.cpp @@ -7,7 +7,7 @@ namespace ModuleIO // Helper function to print dipole moment components // name: descriptive name for the dipole moment type // px, py, pz: dipole moment components in x, y, z directions -void printDipoleMoment(std::ofstream& ofs_running, const std::string& name, +void print_dipole_moment(std::ofstream& ofs_running, const std::string& name, double px, double py, double pz) { ofs_running << " " << name << std::endl; @@ -141,7 +141,7 @@ void write_dipole(const UnitCell& ucell, } // Output electronic dipole moment - printDipoleMoment(ofs_running, "Electronic dipole moment", + print_dipole_moment(ofs_running, "Electronic dipole moment", dipole_elec[0], dipole_elec[1], dipole_elec[2]); // Write to file: step index and three dipole components @@ -170,7 +170,7 @@ void write_dipole(const UnitCell& ucell, } // Output ionic dipole moment - printDipoleMoment(ofs_running, "Ionic dipole moment", + print_dipole_moment(ofs_running, "Ionic dipole moment", dipole_ion[0], dipole_ion[1], dipole_ion[2]); // Calculate total dipole moment @@ -182,7 +182,7 @@ void write_dipole(const UnitCell& ucell, } // Output total dipole moment - printDipoleMoment(ofs_running, "Total dipole moment", + print_dipole_moment(ofs_running, "Total dipole moment", dipole[0], dipole[1], dipole[2]); // Calculate and output total dipole moment norm From 473c08d06180358f221661c57161175633166b75 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:32:18 +0800 Subject: [PATCH 06/11] Cleanup: remove unused 'is' parameter from write_dipole The 'is' (spin channel index) parameter was not used in the write_dipole function: - Electron dipole moment uses rho_save (already selected by spin) - Ionic dipole moment uses ucell (spin-independent) - File naming is handled by the caller through 'fn' parameter Changes: 1. dipole_io.h - Remove 'is' parameter from function declaration 2. write_dipole.cpp - Remove 'is' parameter from function definition 3. ctrl_output_fp.cpp - Update function call to remove 'is' argument --- source/source_io/module_ctrl/ctrl_output_fp.cpp | 2 +- source/source_io/module_dipole/dipole_io.h | 1 - source/source_io/module_dipole/write_dipole.cpp | 2 -- 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 13ba27ebc44..710f4b87307 100644 --- a/source/source_io/module_ctrl/ctrl_output_fp.cpp +++ b/source/source_io/module_ctrl/ctrl_output_fp.cpp @@ -206,7 +206,7 @@ void ctrl_output_fp(UnitCell& ucell, { std::stringstream ss_dipole; ss_dipole << global_out_dir << "dipole_s" << is + 1 << ".txt"; - ModuleIO::write_dipole(ucell, chr.rho_save[is], pw_rhod, is, istep, ss_dipole.str(), GlobalV::ofs_running); + ModuleIO::write_dipole(ucell, chr.rho_save[is], pw_rhod, istep, ss_dipole.str(), GlobalV::ofs_running); } } diff --git a/source/source_io/module_dipole/dipole_io.h b/source/source_io/module_dipole/dipole_io.h index 1a86f82c885..389d7681f64 100644 --- a/source/source_io/module_dipole/dipole_io.h +++ b/source/source_io/module_dipole/dipole_io.h @@ -12,7 +12,6 @@ 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, diff --git a/source/source_io/module_dipole/write_dipole.cpp b/source/source_io/module_dipole/write_dipole.cpp index f9343c27a91..d09317ef518 100644 --- a/source/source_io/module_dipole/write_dipole.cpp +++ b/source/source_io/module_dipole/write_dipole.cpp @@ -38,7 +38,6 @@ void print_dipole_moment(std::ofstream& ofs_running, const std::string& name, // - ucell: unit cell containing atomic structure and lattice information // - rho_save: electron density on real-space grid // - rhopw: plane wave basis information including grid dimensions -// - is: spin channel index // - istep: current time step // - fn: output file name // - ofs_running: output stream for logging @@ -46,7 +45,6 @@ void print_dipole_moment(std::ofstream& ofs_running, const std::string& name, 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, From 225a6cf1bb7f7903cdaae1445457f36c3be4be26 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:45:16 +0800 Subject: [PATCH 07/11] update index.rst --- docs/advanced/output_files/index.rst | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/advanced/output_files/index.rst b/docs/advanced/output_files/index.rst index 75bb501517d..0754799ea78 100644 --- a/docs/advanced/output_files/index.rst +++ b/docs/advanced/output_files/index.rst @@ -7,4 +7,5 @@ Detailed Introduction of the Output Files output-specification running_scf.log - \ No newline at end of file + output_dipole.md + From 613277d5f5e9fc5ad1a63136c220cf5680d3937b Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:49:08 +0800 Subject: [PATCH 08/11] add back #ifdef __LCAO --- source/source_io/module_ctrl/ctrl_output_td.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/source/source_io/module_ctrl/ctrl_output_td.cpp b/source/source_io/module_ctrl/ctrl_output_td.cpp index 6d7fbcc6658..204c755da1b 100644 --- a/source/source_io/module_ctrl/ctrl_output_td.cpp +++ b/source/source_io/module_ctrl/ctrl_output_td.cpp @@ -28,6 +28,7 @@ void ctrl_output_td(const UnitCell& ucell, { ModuleBase::TITLE("ModuleIO", "ctrl_output_td"); +#ifdef __LCAO // (1) Write current information const elecstate::ElecStateLCAO>* pelec_lcao = dynamic_cast>*>(pelec); @@ -52,6 +53,7 @@ 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); } +#endif // __LCAO // (3) Output file for restart if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0 From 74108d7438cd60bb673f00fd816b619e725eaf03 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:49:37 +0800 Subject: [PATCH 09/11] add output_dipole.md --- docs/advanced/output_files/output_dipole.md | 67 +++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 docs/advanced/output_files/output_dipole.md diff --git a/docs/advanced/output_files/output_dipole.md b/docs/advanced/output_files/output_dipole.md new file mode 100644 index 00000000000..064cbfd4efb --- /dev/null +++ b/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. From 2afd4a8ec3ed9a2643b81fcb99bfd919188cdb80 Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 21:51:33 +0800 Subject: [PATCH 10/11] add INPUT --- examples/18_md/01_lcao_gamma_Si8/INPUT | 35 ++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 examples/18_md/01_lcao_gamma_Si8/INPUT diff --git a/examples/18_md/01_lcao_gamma_Si8/INPUT b/examples/18_md/01_lcao_gamma_Si8/INPUT new file mode 100644 index 00000000000..017c3ed8c8c --- /dev/null +++ b/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 From 64fdf6815307bc6500f0d4dd0c92fc05b9b30ade Mon Sep 17 00:00:00 2001 From: abacus_fixer Date: Sat, 16 May 2026 22:03:58 +0800 Subject: [PATCH 11/11] fix --- source/source_io/module_ctrl/ctrl_output_td.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/source/source_io/module_ctrl/ctrl_output_td.cpp b/source/source_io/module_ctrl/ctrl_output_td.cpp index 204c755da1b..d1d2a9aeb5c 100644 --- a/source/source_io/module_ctrl/ctrl_output_td.cpp +++ b/source/source_io/module_ctrl/ctrl_output_td.cpp @@ -53,8 +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); } -#endif // __LCAO - // (3) Output file for restart if (PARAM.inp.out_freq_td > 0) // default value of out_freq_td is 0 { @@ -71,6 +69,7 @@ void ctrl_output_td(const UnitCell& ucell, } } } +#endif // __LCAO } template void ctrl_output_td(const UnitCell&,