From 3379cde62fff51c63a75364d1d79a6ea35916fb9 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Fri, 9 Oct 2026 15:23:39 +0800 Subject: [PATCH] Fix: Output dipole moments in Cartesian coordinates for arbitrary lattices --- docs/advanced/output_files/output_dipole.md | 6 +- source/source_io/module_dipole/dipole_io.h | 2 - .../source_io/module_dipole/write_dipole.cpp | 98 ++++--------------- 3 files changed, 24 insertions(+), 82 deletions(-) diff --git a/docs/advanced/output_files/output_dipole.md b/docs/advanced/output_files/output_dipole.md index 064cbfd4efb..dfb57930d76 100644 --- a/docs/advanced/output_files/output_dipole.md +++ b/docs/advanced/output_files/output_dipole.md @@ -46,7 +46,9 @@ 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 +- `dipole_x, dipole_y, dipole_z`: The Cartesian x, y, z components of the electronic dipole moment for the corresponding spin channel, in units of elementary charge times Bohr ($e a_0$). + +These components refer to the Cartesian axes, including for nonorthogonal or rotated cells; they are not components along the lattice vectors. The electronic dipole includes the negative sign of the electron charge. Example output: @@ -64,4 +66,4 @@ During the calculation, the dipole moment is also printed in the `running_*.log` - 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. +The log reports Cartesian components in the same units. The total dipole moment is the sum of the electronic and ionic dipoles, and its norm is calculated from these Cartesian components. The `dipole_s${spin}.txt` files contain only the electronic contribution. diff --git a/source/source_io/module_dipole/dipole_io.h b/source/source_io/module_dipole/dipole_io.h index 3911b5f0375..c38f5f93b7c 100644 --- a/source/source_io/module_dipole/dipole_io.h +++ b/source/source_io/module_dipole/dipole_io.h @@ -17,8 +17,6 @@ void write_dipole(const UnitCell& ucell, std::ofstream& ofs_running, const int& precision = 11); -double prepare(const UnitCell& cell, int& dir); - } // namespace ModuleIO #endif diff --git a/source/source_io/module_dipole/write_dipole.cpp b/source/source_io/module_dipole/write_dipole.cpp index c8a31671c0b..74020c0492e 100644 --- a/source/source_io/module_dipole/write_dipole.cpp +++ b/source/source_io/module_dipole/write_dipole.cpp @@ -67,31 +67,14 @@ void write_dipole(const UnitCell& ucell, 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) - { - bmod[i] = prepare(ucell, i); - if (bmod[i] < small_value) - { - 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 || 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}; + // Accumulate the electronic first moment in fractional coordinates. + double dipole_elec_frac[3] = {0.0, 0.0, 0.0}; // Precompute inverse grid dimensions for performance double inv_nx = 1.0 / static_cast(rhopw->nx); @@ -100,7 +83,7 @@ void write_dipole(const UnitCell& ucell, // 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) +#pragma omp parallel for reduction(- : dipole_elec_frac[ : 3]) schedule(static) for (int ir = 0; ir < rhopw->nrxx; ++ir) { // Convert 1D index to 3D indices @@ -114,25 +97,22 @@ void write_dipole(const UnitCell& ucell, 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[ir] * x; - dipole_elec[1] -= rho[ir] * y; - dipole_elec[2] -= rho[ir] * z; + // The negative sign accounts for the electron charge. + dipole_elec_frac[0] -= rho[ir] * x; + dipole_elec_frac[1] -= rho[ir] * y; + dipole_elec_frac[2] -= rho[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] * vol_factor; - } + Parallel_Reduce::reduce_pool(dipole_elec_frac[0]); + Parallel_Reduce::reduce_pool(dipole_elec_frac[1]); + Parallel_Reduce::reduce_pool(dipole_elec_frac[2]); + + // Lattice vectors are rows: r_cart = r_frac * latvec * lat0. + // Apply the full matrix to retain shear, rotation, and axis signs. + const double vol_factor = ucell.omega / static_cast(rhopw->nxyz); + const ModuleBase::Vector3 elec_frac(dipole_elec_frac[0], dipole_elec_frac[1], dipole_elec_frac[2]); + const ModuleBase::Vector3 dipole_elec = elec_frac * ucell.latvec * (ucell.lat0 * vol_factor); // Output electronic dipole moment print_dipole_moment(ofs_running, "Electronic dipole moment", dipole_elec[0], dipole_elec[1], dipole_elec[2]); @@ -142,9 +122,8 @@ void write_dipole(const UnitCell& ucell, << 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}; + // Accumulate Z_v * taud, then transform to Cartesian coordinates. + ModuleBase::Vector3 ion_frac(0.0, 0.0, 0.0); for (int i = 0; i < 3; ++i) { for (int it = 0; it < ucell.ntype; ++it) @@ -154,11 +133,10 @@ void write_dipole(const UnitCell& ucell, { sum += ucell.atoms[it].taud[ia][i]; } - dipole_ion[i] += sum * ucell.atoms[it].ncpp.zv; + ion_frac[i] += sum * ucell.atoms[it].ncpp.zv; } - // Convert to Cartesian coordinates - dipole_ion[i] *= ucell.lat0 / bmod[i]; } + const ModuleBase::Vector3 dipole_ion = ion_frac * ucell.latvec * ucell.lat0; // Output ionic dipole moment print_dipole_moment(ofs_running, "Ionic dipole moment", dipole_ion[0], dipole_ion[1], dipole_ion[2]); @@ -189,40 +167,4 @@ void write_dipole(const UnitCell& ucell, } } -// 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}; - - // Select the appropriate reciprocal lattice vector components - switch (dir) - { - 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"); - } - - // 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