Refactor: Remove cal_ux from UnitCell.

source/Makefile.Objects

@@ -226,6 +226,7 @@ OBJS_ELECSTAT=elecstate.o\
     H_Hartree_pw.o\
     H_TDDFT_pw.o\
     pot_xc.o\
+    cal_ux.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
       elecstate_lcao_cal_tau.o\

source/module_cell/test/CMakeLists.txt

@@ -102,7 +102,7 @@ add_test(NAME cell_parallel_kpoints_test
 AddTest(
   TARGET cell_unitcell_test
   LIBS parameter  ${math_libs} base device cell_info  symmetry
-  SOURCES unitcell_test.cpp ../../module_io/output.cpp 
+  SOURCES unitcell_test.cpp ../../module_io/output.cpp ../../module_elecstate/cal_ux.cpp
 
 )
 

source/module_cell/test/support/mock_unitcell.cpp

@@ -8,10 +8,6 @@
    to avoid using UnitCell functions because there is GLobalC, which will bring
    endless compile troubles like undefined behavior"
 */
-void UnitCell::cal_ux() {}
-bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
-    return true;
-}
 void UnitCell::set_iat2iwt(const int& npol_in) {}
 UnitCell::UnitCell() {
     Coordinate = "Direct";

source/module_cell/test/unitcell_test.cpp

@@ -4,6 +4,7 @@
 #include "module_parameter/parameter.h"
 #undef private
 #include "module_cell/unitcell.h"
+#include "module_elecstate/cal_ux.h"
 
 #include "memory"
 #include "module_base/global_variable.h"
@@ -97,7 +98,7 @@ Magnetism::~Magnetism()
  *   - UpdateVel
  *     - update_vel(const ModuleBase::Vector3<double>* vel_in)
  *   - CalUx
- *     - cal_ux(): calculate magnetic moments of cell
+ *     - cal_ux(UnitCell& ucell): calculate magnetic moments of cell
  *   - ReadOrbFile
  *     - read_orb_file(): read header part of orbital file
  *   - ReadOrbFileWarning
@@ -610,7 +611,7 @@ TEST_F(UcellTest, JudgeParallel)
 {
     ModuleBase::Vector3<double> b(1.0, 1.0, 1.0);
     double a[3] = {1.0, 1.0, 1.0};
-    EXPECT_TRUE(ucell->judge_parallel(a, b));
+    EXPECT_TRUE(elecstate::judge_parallel(a, b));
 }
 
 TEST_F(UcellTest, Index)
@@ -1034,7 +1035,7 @@ TEST_F(UcellTest, CalUx1)
     ucell->atoms[0].m_loc_[0].set(0, -1, 0);
     ucell->atoms[1].m_loc_[0].set(1, 1, 1);
     ucell->atoms[1].m_loc_[1].set(0, 0, 0);
-    ucell->cal_ux();
+    elecstate::cal_ux(*ucell);
     EXPECT_FALSE(ucell->magnet.lsign_);
     EXPECT_DOUBLE_EQ(ucell->magnet.ux_[0], 0);
     EXPECT_DOUBLE_EQ(ucell->magnet.ux_[1], -1);
@@ -1050,7 +1051,7 @@ TEST_F(UcellTest, CalUx2)
     ucell->atoms[1].m_loc_[0].set(1, 1, 1);
     ucell->atoms[1].m_loc_[1].set(0, 0, 0);
     //(0,0,0) is also parallel to (1,1,1)
-    ucell->cal_ux();
+    elecstate::cal_ux(*ucell);
     EXPECT_TRUE(ucell->magnet.lsign_);
     EXPECT_NEAR(ucell->magnet.ux_[0], 0.57735, 1e-5);
     EXPECT_NEAR(ucell->magnet.ux_[1], 0.57735, 1e-5);

source/module_cell/unitcell.cpp

@@ -473,70 +473,6 @@ void UnitCell::bcast_atoms_tau() {
 #endif
 }
 
-void UnitCell::cal_ux() {
-    double amag, uxmod;
-    int starting_it = 0;
-    int starting_ia = 0;
-    bool is_paraller;
-    // do not sign feature in teh general case
-    magnet.lsign_ = false;
-    ModuleBase::GlobalFunc::ZEROS(magnet.ux_, 3);
-
-    for (int it = 0; it < ntype; it++) {
-        for (int ia = 0; ia < atoms[it].na; ia++) {
-            amag = pow(atoms[it].m_loc_[ia].x, 2)
-                   + pow(atoms[it].m_loc_[ia].y, 2)
-                   + pow(atoms[it].m_loc_[ia].z, 2);
-            if (amag > 1e-6) {
-                magnet.ux_[0] = atoms[it].m_loc_[ia].x;
-                magnet.ux_[1] = atoms[it].m_loc_[ia].y;
-                magnet.ux_[2] = atoms[it].m_loc_[ia].z;
-                starting_it = it;
-                starting_ia = ia;
-                magnet.lsign_ = true;
-                break;
-            }
-        }
-        if (magnet.lsign_) {
-            break;
-}
-    }
-    // whether the initial magnetizations is parallel
-    for (int it = starting_it; it < ntype; it++) {
-        for (int ia = 0; ia < atoms[it].na; ia++) {
-            if (it > starting_it || ia > starting_ia) {
-                magnet.lsign_
-                    = magnet.lsign_
-                      && judge_parallel(magnet.ux_, atoms[it].m_loc_[ia]);
-            }
-        }
-    }
-    if (magnet.lsign_) {
-        uxmod = pow(magnet.ux_[0], 2) + pow(magnet.ux_[1], 2)
-                + pow(magnet.ux_[2], 2);
-        if (uxmod < 1e-6) {
-            ModuleBase::WARNING_QUIT("cal_ux", "wrong uxmod");
-        }
-        for (int i = 0; i < 3; i++) {
-            magnet.ux_[i] *= 1 / sqrt(uxmod);
-        }
-        //       std::cout<<"    Fixed quantization axis for GGA: "
-        //<<std::setw(10)<<ux[0]<<"  "<<std::setw(10)<<ux[1]<<"
-        //"<<std::setw(10)<<ux[2]<<std::endl;
-    }
-    return;
-}
-
-bool UnitCell::judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
-    bool jp = false;
-    double cross;
-    cross = pow((a[1] * b.z - a[2] * b.y), 2)
-            + pow((a[2] * b.x - a[0] * b.z), 2)
-            + pow((a[0] * b.y - a[1] * b.x), 2);
-    jp = (fabs(cross) < 1e-6);
-    return jp;
-}
-
 //==============================================================
 // Calculate various lattice related quantities for given latvec
 //==============================================================

source/module_cell/unitcell.h

@@ -19,8 +19,6 @@ class UnitCell {
 
     bool set_atom_flag; // added on 2009-3-8 by mohan
     Magnetism magnet;   // magnetism Yu Liu 2021-07-03
-    void cal_ux();
-    bool judge_parallel(double a[3], ModuleBase::Vector3<double> b);
     std::vector<std::vector<double>> atom_mulliken; //[nat][nspin]
     int n_mag_at;
 

source/module_elecstate/CMakeLists.txt

@@ -31,6 +31,7 @@ list(APPEND objects
     fp_energy.cpp
     magnetism.cpp
     occupy.cpp
+    cal_ux.cpp
 )
 
 if(ENABLE_LCAO)

source/module_elecstate/cal_ux.cpp

@@ -0,0 +1,68 @@
+#include "cal_ux.h"
+
+namespace elecstate {
+
+void cal_ux(UnitCell& ucell) {
+    double amag, uxmod;
+    int starting_it = 0;
+    int starting_ia = 0;
+    bool is_paraller;
+    // do not sign feature in teh general case
+    ucell.magnet.lsign_ = false;
+    ModuleBase::GlobalFunc::ZEROS(ucell.magnet.ux_, 3);
+
+    for (int it = 0; it < ucell.ntype; it++) {
+        for (int ia = 0; ia < ucell.atoms[it].na; ia++) {
+            amag = pow(ucell.atoms[it].m_loc_[ia].x, 2)
+                   + pow(ucell.atoms[it].m_loc_[ia].y, 2)
+                   + pow(ucell.atoms[it].m_loc_[ia].z, 2);
+            if (amag > 1e-6) {
+                ucell.magnet.ux_[0] = ucell.atoms[it].m_loc_[ia].x;
+                ucell.magnet.ux_[1] = ucell.atoms[it].m_loc_[ia].y;
+                ucell.magnet.ux_[2] = ucell.atoms[it].m_loc_[ia].z;
+                starting_it = it;
+                starting_ia = ia;
+                ucell.magnet.lsign_ = true;
+                break;
+            }
+        }
+        if (ucell.magnet.lsign_) {
+            break;
+        }
+    }
+    // whether the initial magnetizations is parallel
+    for (int it = starting_it; it < ucell.ntype; it++) {
+        for (int ia = 0; ia < ucell.atoms[it].na; ia++) {
+            if (it > starting_it || ia > starting_ia) {
+                ucell.magnet.lsign_
+                    = ucell.magnet.lsign_
+                      && judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]);
+            }
+        }
+    }
+    if (ucell.magnet.lsign_) {
+        uxmod = pow(ucell.magnet.ux_[0], 2) + pow(ucell.magnet.ux_[1], 2)
+                + pow(ucell.magnet.ux_[2], 2);
+        if (uxmod < 1e-6) {
+            ModuleBase::WARNING_QUIT("cal_ux", "wrong uxmod");
+        }
+        for (int i = 0; i < 3; i++) {
+            ucell.magnet.ux_[i] *= 1 / sqrt(uxmod);
+        }
+        //       std::cout<<"    Fixed quantization axis for GGA: "
+        //<<std::setw(10)<<ux[0]<<"  "<<std::setw(10)<<ux[1]<<"
+        //"<<std::setw(10)<<ux[2]<<std::endl;
+    }
+    return;
+}
+
+bool judge_parallel(double a[3], ModuleBase::Vector3<double> b) {
+    bool jp = false;
+    double cross;
+    cross = pow((a[1] * b.z - a[2] * b.y), 2)
+            + pow((a[2] * b.x - a[0] * b.z), 2)
+            + pow((a[0] * b.y - a[1] * b.x), 2);
+    jp = (fabs(cross) < 1e-6);
+    return jp;
+}
+}

source/module_elecstate/cal_ux.h

@@ -0,0 +1,14 @@
+#ifndef CAL_UX_H
+#define CAL_UX_H
+
+#include "module_cell/unitcell.h"
+
+namespace elecstate {
+
+    void cal_ux(UnitCell& ucell);
+
+    bool judge_parallel(double a[3], ModuleBase::Vector3<double> b);
+
+}
+
+#endif

source/module_esolver/esolver_ks_lcao.cpp

@@ -40,6 +40,7 @@
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
+#include "module_elecstate/cal_ux.h"
 
 #include <memory>
 #ifdef __EXX
@@ -593,7 +594,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
 
             if (PARAM.inp.nspin == 4)
             {
-                ucell.cal_ux();
+                elecstate::cal_ux(ucell);
             }
 
             //! update the potentials by using new electron charge density
@@ -829,7 +830,7 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(UnitCell& ucell, const int istep, const
     {
         if (PARAM.inp.nspin == 4)
         {
-            ucell.cal_ux();
+            elecstate::cal_ux(ucell);
         }
         this->pelec->pot->update_from_charge(this->pelec->charge, &ucell);
         this->pelec->f_en.descf = this->pelec->cal_delta_escf();

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -29,6 +29,7 @@
 #include "module_hsolver/hsolver_lcao.h"
 #include "module_parameter/parameter.h"
 #include "module_psi/psi.h"
+#include "module_elecstate/cal_ux.h"
 
 //-----force& stress-------------------
 #include "module_hamilt_lcao/hamilt_lcaodft/FORCE_STRESS.h"
@@ -239,7 +240,7 @@ void ESolver_KS_LCAO_TDDFT::update_pot(UnitCell& ucell, const int istep, const i
     {
         if (PARAM.inp.nspin == 4)
         {
-            ucell.cal_ux();
+            elecstate::cal_ux(ucell);
         }
         this->pelec->pot->update_from_charge(this->pelec->charge, &ucell);
         this->pelec->f_en.descf = this->pelec->cal_delta_escf();

source/module_esolver/esolver_ks_pw.cpp

@@ -8,6 +8,7 @@
 #include "module_hamilt_general/module_ewald/H_Ewald_pw.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
+#include "module_elecstate/cal_ux.h"
 //-----force-------------------
 #include "module_hamilt_pw/hamilt_pwdft/forces.h"
 //-----stress------------------
@@ -298,7 +299,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     //! the direction of ux is used in noncoline_rho
     if (PARAM.inp.nspin == 4)
     {
-        ucell.cal_ux();
+        elecstate::cal_ux(ucell);
     }
 
     //! calculate the total local pseudopotential in real space
@@ -471,7 +472,7 @@ void ESolver_KS_PW<T, Device>::update_pot(UnitCell& ucell, const int istep, cons
     {
         if (PARAM.inp.nspin == 4)
         {
-            ucell.cal_ux();
+            elecstate::cal_ux(ucell);
         }
         this->pelec->pot->update_from_charge(this->pelec->charge, &ucell);
         this->pelec->f_en.descf = this->pelec->cal_delta_escf();

source/module_esolver/esolver_of.cpp

@@ -10,6 +10,7 @@
 #include "module_hamilt_general/module_ewald/H_Ewald_pw.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
+#include "module_elecstate/cal_ux.h"
 //-----force-------------------
 #include "module_hamilt_pw/hamilt_pwdft/forces.h"
 //-----stress------------------
@@ -315,7 +316,7 @@ void ESolver_OF::update_potential(UnitCell& ucell)
     // (1) get dL/dphi
     if (PARAM.inp.nspin == 4)
     {
-        ucell.cal_ux();
+        elecstate::cal_ux(ucell);
     }
 
     this->pelec->pot->update_from_charge(pelec->charge, &ucell); // Hartree + XC + external

source/module_esolver/esolver_of_tool.cpp

@@ -5,6 +5,7 @@
 #include "module_elecstate/potentials/gatefield.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_parameter/parameter.h"
+#include "module_elecstate/cal_ux.h"
 
 namespace ModuleESolver
 {
@@ -134,7 +135,7 @@ void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& uce
 
     if (PARAM.inp.nspin == 4) 
     {
-        ucell.cal_ux();
+        elecstate::cal_ux(ucell);
     }
     this->pelec->pot->update_from_charge(this->ptemp_rho_, &ucell);
     ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v();
@@ -172,9 +173,10 @@ void ESolver_OF::cal_dEdtheta(double** ptemp_phi, Charge* temp_rho, UnitCell& uc
 {
     double* dphi_dtheta = new double[this->pw_rho->nrxx];
 
-    if (PARAM.inp.nspin == 4) {
-        ucell.cal_ux();
-}
+    if (PARAM.inp.nspin == 4) 
+    {
+        elecstate::cal_ux(ucell);
+    }
     this->pelec->pot->update_from_charge(temp_rho, &ucell);
     ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v();
 

source/module_esolver/lcao_before_scf.cpp

@@ -3,6 +3,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_elecstate/cal_ux.h"
 //
 #include "module_base/timer.h"
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
@@ -246,7 +247,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     //=========================================================
     if (PARAM.inp.nspin == 4)
     {
-        ucell.cal_ux();
+        elecstate::cal_ux(ucell);
     }
 
     // Peize Lin add 2016-12-03

source/module_esolver/lcao_others.cpp

@@ -3,6 +3,7 @@
 #include "module_hamilt_lcao/hamilt_lcaodft/hamilt_lcao.h"
 #include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
+#include "module_elecstate/cal_ux.h"
 //
 #include "module_base/timer.h"
 #include "module_cell/module_neighbor/sltk_atom_arrange.h"
@@ -247,7 +248,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     //=========================================================
     if (PARAM.inp.nspin == 4)
     {
-        ucell.cal_ux();
+        elecstate::cal_ux(ucell);
     }
 
     // pelec should be initialized before these calculations