Fix: Output dipole moments in Cartesian coordinates for arbitrary lattices in RT-TDDFT

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.

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

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<double>(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<double>(j) * inv_ny;
         double z = static_cast<double>(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<double>(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<double>(rhopw->nxyz);
+    const ModuleBase::Vector3<double> elec_frac(dipole_elec_frac[0], dipole_elec_frac[1], dipole_elec_frac[2]);
+    const ModuleBase::Vector3<double> 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<double> 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<double> 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