add two method of orb_mag

Author: dyzheng

source/source_esolver/esolver_ks_lcao.cpp

@@ -21,6 +21,7 @@
 #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
 #include "source_lcao/LCAO_set.h" // mohan add 20251111
 #include "source_lcao/module_orbital_mag/orbital_mag.h"
+#include "source_lcao/module_orbital_mag/orbital_mag_dm.h"
 
 namespace ModuleESolver
 {
@@ -582,6 +583,35 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(UnitCell& ucell, const int istep, const
             GlobalV::ofs_running << " M_y = " << M_orb.y << std::endl;
             GlobalV::ofs_running << " M_z = " << M_orb.z << std::endl;
         }
+
+        //! 5) Calculate orbital magnetic moment using density matrix method if requested
+        if (PARAM.inp.out_orbital_mag_dm)
+        {
+            hamilt::OrbitalMagDM<TK> orbital_mag_dm(ucell,
+                                                    this->kv,
+                                                    *this->dmat.dm,
+                                                    hamilt_lcao->getHR(),
+                                                    hamilt_lcao->getSR(),
+                                                    &this->pv,
+                                                    this->orb_);
+            ModuleBase::Vector3<double> M_itin = orbital_mag_dm.calculate_orbital_moment();
+            ModuleBase::Vector3<double> M_local = orbital_mag_dm.calculate_local_moment();
+            ModuleBase::Vector3<double> M_total = M_itin + M_local;
+
+            GlobalV::ofs_running << "\n ORBITAL MAGNETIC MOMENT (DM method, Bohr magneton):" << std::endl;
+            GlobalV::ofs_running << " --- Itinerant contribution ---" << std::endl;
+            GlobalV::ofs_running << " M_x = " << M_itin.x << std::endl;
+            GlobalV::ofs_running << " M_y = " << M_itin.y << std::endl;
+            GlobalV::ofs_running << " M_z = " << M_itin.z << std::endl;
+            GlobalV::ofs_running << " --- Local (on-site) contribution ---" << std::endl;
+            GlobalV::ofs_running << " M_x = " << M_local.x << std::endl;
+            GlobalV::ofs_running << " M_y = " << M_local.y << std::endl;
+            GlobalV::ofs_running << " M_z = " << M_local.z << std::endl;
+            GlobalV::ofs_running << " --- Total orbital moment ---" << std::endl;
+            GlobalV::ofs_running << " M_x = " << M_total.x << std::endl;
+            GlobalV::ofs_running << " M_y = " << M_total.y << std::endl;
+            GlobalV::ofs_running << " M_z = " << M_total.z << std::endl;
+        }
     }
 
     ModuleBase::timer::tick("ESolver_KS_LCAO", "after_scf");

source/source_io/cal_r_overlap_R.cpp

@@ -516,6 +516,55 @@ ModuleBase::Vector3<double> cal_r_overlap_R::get_psi_r_psi(const ModuleBase::Vec
     return temp_prp;
 }
 
+ModuleBase::Vector3<std::complex<double>> cal_r_overlap_R::get_psi_L_psi(
+    const ModuleBase::Vector3<double>& R1,
+    const int& T1,
+    const int& L1,
+    const int& m1,
+    const int& N1,
+    const ModuleBase::Vector3<double>& R2,
+    const int& T2,
+    const int& L2,
+    const int& m2,
+    const int& N2)
+{
+    ModuleBase::Vector3<double> origin_point(0.0, 0.0, 0.0);
+    double factor = sqrt(ModuleBase::FOUR_PI / 3.0);
+    const ModuleBase::Vector3<double>& distance = R2 - R1;
+
+    // grad_r_beta[beta] = ∂/∂R · ⟨φ_μ|(r-R_B)_beta|φ_ν⟩  (a Vector3<double> for each beta)
+    // Following the same sign convention as get_psi_r_psi:
+    //   mA2=0 → z component (positive factor)
+    //   mA2=1 → x component (negative factor)
+    //   mA2=2 → y component (negative factor)
+
+    // grad of ⟨φ|(r-R_B)_x|φ⟩ w.r.t. R
+    ModuleBase::Vector3<double> grad_rx
+        = -1.0 * factor
+          * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, 1, m2);
+
+    // grad of ⟨φ|(r-R_B)_y|φ⟩ w.r.t. R
+    ModuleBase::Vector3<double> grad_ry
+        = -1.0 * factor
+          * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, 2, m2);
+
+    // grad of ⟨φ|(r-R_B)_z|φ⟩ w.r.t. R
+    ModuleBase::Vector3<double> grad_rz
+        = factor
+          * center2_orb21_r[T1][T2][L1][N1][L2].at(N2).cal_grad_overlap(origin_point, distance, m1, 0, m2);
+
+    // L_x = i * [∂/∂R_y ⟨r_z⟩ - ∂/∂R_z ⟨r_y⟩]
+    // L_y = i * [∂/∂R_z ⟨r_x⟩ - ∂/∂R_x ⟨r_z⟩]
+    // L_z = i * [∂/∂R_x ⟨r_y⟩ - ∂/∂R_y ⟨r_x⟩]
+    const std::complex<double> imag_i(0.0, 1.0);
+
+    std::complex<double> Lx = imag_i * (grad_rz.y - grad_ry.z);
+    std::complex<double> Ly = imag_i * (grad_rx.z - grad_rz.x);
+    std::complex<double> Lz = imag_i * (grad_ry.x - grad_rx.y);
+
+    return ModuleBase::Vector3<std::complex<double>>(Lx, Ly, Lz);
+}
+
 void cal_r_overlap_R::get_psi_r_beta(const UnitCell& ucell,
                                      std::vector<std::vector<double>>& nlm,
                                      const ModuleBase::Vector3<double>& R1,

source/source_io/cal_r_overlap_R.h

@@ -5,6 +5,7 @@
 #include "source_base/sph_bessel_recursive.h"
 #include "source_base/vector3.h"
 #include "source_base/ylm.h"
+#include <complex>
 #include "source_basis/module_ao/ORB_atomic_lm.h"
 #include "source_basis/module_ao/ORB_gaunt_table.h"
 #include "source_basis/module_ao/ORB_read.h"
@@ -46,6 +47,18 @@ class cal_r_overlap_R
       const int& m2,
       const int& N2
     );
+    ModuleBase::Vector3<std::complex<double>> get_psi_L_psi(
+      const ModuleBase::Vector3<double>& R1,
+      const int& T1,
+      const int& L1,
+      const int& m1,
+      const int& N1,
+      const ModuleBase::Vector3<double>& R2,
+      const int& T2,
+      const int& L2,
+      const int& m2,
+      const int& N2
+    );
     void get_psi_r_beta(
       const UnitCell& ucell,
       std::vector<std::vector<double>>& nlm,

source/source_io/module_parameter/input_parameter.h

@@ -402,6 +402,7 @@ struct Input_para
     bool out_mat_r = false;               ///< jingan add 2019-8-14, output r(R) matrix.
     int out_wfc_lcao = 0;                 ///< output the wave functions in local basis.
     bool out_orbital_mag = false;         ///< output orbital magnetic moment (LCAO multi-k only)
+    bool out_orbital_mag_dm = false;      ///< output orbital magnetic moment using density matrix method
     bool out_dipole = false;              ///< output the dipole or not
     bool out_efield = false;              ///< output the efield or not
     int out_current = 0;                  ///< output the current or not

source/source_io/module_parameter/read_input_item_output.cpp

@@ -163,6 +163,18 @@ void ReadInput::item_output()
         };
         this->add_item(item);
     }
+    {
+        Input_Item item("out_orbital_mag_dm");
+        item.annotation = "output orbital magnetic moment using density matrix method (LCAO multi-k only)";
+        read_sync_bool(input.out_orbital_mag_dm);
+        item.check_value = [](const Input_Item& item, const Parameter& para) {
+            if (para.input.basis_type != "lcao" && para.input.out_orbital_mag_dm)
+            {
+                ModuleBase::WARNING_QUIT("ReadInput", "out_orbital_mag_dm is only for lcao");
+            }
+        };
+        this->add_item(item);
+    }
     {
         Input_Item item("out_proj_band");
         item.annotation = "output projected band structure";

source/source_io/test/CMakeLists.txt

@@ -279,6 +279,25 @@ AddTest(
 )
 endif()
 
+if(ENABLE_LCAO)
+AddTest(
+  TARGET MODULE_IO_cal_r_overlap_R_L_test
+  LIBS parameter base ${math_libs} device numerical_atomic_orbitals container orb
+  SOURCES
+  cal_r_overlap_R_L_test.cpp
+  ../cal_r_overlap_R.cpp
+  ../single_R_io.cpp
+  ../../source_cell/atom_spec.cpp
+  ../../source_cell/parallel_kpoints.cpp
+  ../../source_cell/test/support/mock_unitcell.cpp
+  ../../source_lcao/center2_orb.cpp
+  ../../source_lcao/center2_orb-orb11.cpp
+  ../../source_lcao/center2_orb-orb21.cpp
+  ../../source_lcao/center2_orb-orb22.cpp
+  ../orb_io.cpp
+)
+endif()
+
 AddTest(
   TARGET MODULE_IO_write_elf_logic_test
   SOURCES write_elf_logic_test.cpp