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 + 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. 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 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; } diff --git a/source/source_io/module_ctrl/ctrl_output_fp.cpp b/source/source_io/module_ctrl/ctrl_output_fp.cpp index 8e3991c8a14..710f4b87307 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, 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 65d0c9ab965..d1d2a9aeb5c 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" @@ -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>* pelec_lcao = dynamic_cast>*>(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 { diff --git a/source/source_io/module_dipole/dipole_io.h b/source/source_io/module_dipole/dipole_io.h index b97810bf757..389d7681f64 100644 --- a/source/source_io/module_dipole/dipole_io.h +++ b/source/source_io/module_dipole/dipole_io.h @@ -2,7 +2,9 @@ #define DIPOLE_IO_H #include "source_basis/module_pw/pw_basis.h" +#include "source_cell/unitcell.h" +#include #include namespace ModuleIO @@ -10,9 +12,9 @@ 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); diff --git a/source/source_io/module_dipole/write_dipole.cpp b/source/source_io/module_dipole/write_dipole.cpp index 56448365bfe..d09317ef518 100644 --- a/source/source_io/module_dipole/write_dipole.cpp +++ b/source/source_io/module_dipole/write_dipole.cpp @@ -1,105 +1,157 @@ #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 -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, - const int& precision) + +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 print_dipole_moment(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 +// 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_{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) +// +// 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 +// - 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& 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 Charge File!"); + ModuleBase::WARNING_QUIT("ModuleIO", "Cannot create dipole file: " + fn); } } + 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 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); - } - -#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++) - { - for (int j = 0; j < rhopw->ny; j++) + if (bmod[i] < small_value) { - 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", + "bmod[" + std::to_string(i) + "] is too small or 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); - - ofs << istep+1 << " " << dipole_elec_x << " " << dipole_elec_y << dipole_elec_z; -#else + // Validate grid dimensions to prevent division by zero + if (rhopw->nx == 0 || rhopw->ny == 0 || rhopw->nz == 0 || + rhopw->nxyz == 0 || rhopw->nplane == 0) + { + 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}; - + + // 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 int i = ir / (rhopw->ny * rhopw->nplane); int j = ir / rhopw->nplane - i * rhopw->ny; int k = ir % rhopw->nplane + rhopw->startz_current; - double x = (double)i / rhopw->nx; - double y = (double)j / rhopw->ny; - double z = (double)k / rhopw->nz; - + + // Convert to fractional coordinates: r_i = i / N_i + // 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; 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 + 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; } - 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]); + // Output electronic dipole moment + print_dipole_moment(ofs_running, "Electronic dipole moment", + dipole_elec[0], dipole_elec[1], dipole_elec[2]); - ofs << std::setprecision(precision) << istep+1 << " " << dipole_elec[0] << " " << dipole_elec[1] << " " - << dipole_elec[2] << std::endl; + // 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; - for (int i = 0; i < 3; ++i) { for (int it = 0; it < ucell.ntype; ++it) @@ -111,65 +163,77 @@ 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 + print_dipole_moment(ofs_running, "Ionic dipole moment", + dipole_ion[0], dipole_ion[1], 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]); - - 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); + // Output total dipole moment + print_dipole_moment(ofs_running, "Total dipole moment", + dipole[0], dipole[1], dipole[2]); -#endif + // Calculate and output total dipole moment norm + // |P_tot| = sqrt(P_tot_x^2 + P_tot_y^2 + P_tot_z^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; } -double ModuleIO::prepare(const UnitCell& cell, int& dir) +// 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 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) - { - bvec[0] = cell.G.e21; - bvec[1] = cell.G.e22; - bvec[2] = cell.G.e23; - } - else if (dir == 2) - { - bvec[0] = cell.G.e31; - bvec[1] = cell.G.e32; - bvec[2] = cell.G.e33; - } - else + + // Select the appropriate reciprocal lattice vector components + switch (dir) { - ModuleBase::WARNING_QUIT("ModuleIO::prepare", "direction is wrong!"); + 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"); } - bmod = sqrt(pow(bvec[0], 2) + pow(bvec[1], 2) + pow(bvec[2], 2)); - return bmod; + + // 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