diff --git a/source/Makefile.Objects b/source/Makefile.Objects index 9a5243f1e3..0d3fddbf71 100644 --- a/source/Makefile.Objects +++ b/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\ diff --git a/source/module_cell/test/CMakeLists.txt b/source/module_cell/test/CMakeLists.txt index 9bd09c154f..6558197058 100644 --- a/source/module_cell/test/CMakeLists.txt +++ b/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 ) diff --git a/source/module_cell/test/support/mock_unitcell.cpp b/source/module_cell/test/support/mock_unitcell.cpp index b91430987a..973e347f31 100644 --- a/source/module_cell/test/support/mock_unitcell.cpp +++ b/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 b) { - return true; -} void UnitCell::set_iat2iwt(const int& npol_in) {} UnitCell::UnitCell() { Coordinate = "Direct"; diff --git a/source/module_cell/test/unitcell_test.cpp b/source/module_cell/test/unitcell_test.cpp index a9fb6503eb..4aefb70870 100644 --- a/source/module_cell/test/unitcell_test.cpp +++ b/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* 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 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); diff --git a/source/module_cell/unitcell.cpp b/source/module_cell/unitcell.cpp index 653e632094..71f460b39a 100755 --- a/source/module_cell/unitcell.cpp +++ b/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: " - //< 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 //============================================================== diff --git a/source/module_cell/unitcell.h b/source/module_cell/unitcell.h index 552bea559c..6fcf87a259 100644 --- a/source/module_cell/unitcell.h +++ b/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 b); std::vector> atom_mulliken; //[nat][nspin] int n_mag_at; diff --git a/source/module_elecstate/CMakeLists.txt b/source/module_elecstate/CMakeLists.txt index 78efe4ead0..8d6f2e5f69 100644 --- a/source/module_elecstate/CMakeLists.txt +++ b/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) diff --git a/source/module_elecstate/cal_ux.cpp b/source/module_elecstate/cal_ux.cpp new file mode 100644 index 0000000000..40e9bdde93 --- /dev/null +++ b/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: " + //< 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; +} +} \ No newline at end of file diff --git a/source/module_elecstate/cal_ux.h b/source/module_elecstate/cal_ux.h new file mode 100644 index 0000000000..2d00770fb2 --- /dev/null +++ b/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 b); + +} + +#endif \ No newline at end of file diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index f18957c01c..f0f9345506 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/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 #ifdef __EXX @@ -593,7 +594,7 @@ void ESolver_KS_LCAO::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::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(); diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp index 92be005770..3c60fb2526 100644 --- a/source/module_esolver/esolver_ks_lcao_tddft.cpp +++ b/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(); diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp index d32c4790f9..f081f5596d 100644 --- a/source/module_esolver/esolver_ks_pw.cpp +++ b/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::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::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(); diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp index 8b9019e40c..a5a68ba14a 100644 --- a/source/module_esolver/esolver_of.cpp +++ b/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 diff --git a/source/module_esolver/esolver_of_tool.cpp b/source/module_esolver/esolver_of_tool.cpp index 9d655d8e61..85c6b640fe 100644 --- a/source/module_esolver/esolver_of_tool.cpp +++ b/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(); diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp index b643ffe4c6..9f5b51481c 100644 --- a/source/module_esolver/lcao_before_scf.cpp +++ b/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::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 diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index b4c05b2af7..75f4164e2f 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/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::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 diff --git a/source/module_hamilt_general/module_xc/test/test_xc3.cpp b/source/module_hamilt_general/module_xc/test/test_xc3.cpp index 3b859d2969..13255831e1 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc3.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc3.cpp @@ -48,7 +48,7 @@ class XCTest_GRADCORR : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - ucell.cal_ux(); + elecstate::cal_ux(ucell); chr.rho = new double*[4]; chr.rho[0] = new double[5]; diff --git a/source/module_hamilt_general/module_xc/test/test_xc5.cpp b/source/module_hamilt_general/module_xc/test/test_xc5.cpp index 2eed4b3508..67e81b4208 100644 --- a/source/module_hamilt_general/module_xc/test/test_xc5.cpp +++ b/source/module_hamilt_general/module_xc/test/test_xc5.cpp @@ -44,7 +44,7 @@ class XCTest_VXC : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - ucell.cal_ux(); + elecstate::cal_ux(ucell); ucell.omega = 1; chr.rhopw = &(rhopw); @@ -142,7 +142,7 @@ class XCTest_VXC_Libxc : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - ucell.cal_ux(); + elecstate::cal_ux(ucell); ucell.omega = 1; chr.rhopw = &(rhopw); @@ -240,7 +240,7 @@ class XCTest_VXC_meta : public XCTest ucell.tpiba = 1; ucell.magnet.lsign_ = true; - ucell.cal_ux(); + elecstate::cal_ux(ucell); ucell.omega = 1; chr.rhopw = &(rhopw); diff --git a/source/module_hamilt_general/module_xc/test/xc3_mock.h b/source/module_hamilt_general/module_xc/test/xc3_mock.h index 8ecccd932d..f0848e8bf7 100644 --- a/source/module_hamilt_general/module_xc/test/xc3_mock.h +++ b/source/module_hamilt_general/module_xc/test/xc3_mock.h @@ -190,15 +190,18 @@ Charge::~Charge(){}; Magnetism::Magnetism(){}; Magnetism::~Magnetism(){}; -void UnitCell::cal_ux() +namespace elecstate { - magnet.lsign_ = false; + void cal_ux(UnitCell& ucell) + { + ucell.magnet.lsign_ = false; - magnet.ux_[0] = 0; - magnet.ux_[1] = 1; - magnet.ux_[2] = 2; + ucell.magnet.ux_[0] = 0; + ucell.magnet.ux_[1] = 1; + ucell.magnet.ux_[2] = 2; - magnet.lsign_ = true; + ucell.magnet.lsign_ = true; + }; } #ifdef __LCAO diff --git a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp index f2f36620c4..417c812059 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/forces_cc.cpp @@ -15,6 +15,7 @@ #include "module_hamilt_general/module_ewald/H_Ewald_pw.h" #include "module_hamilt_general/module_surchem/surchem.h" #include "module_hamilt_general/module_vdw/vdw.h" +#include "module_elecstate/cal_ux.h" #ifdef _OPENMP #include @@ -70,9 +71,10 @@ void Forces::cal_force_cc(ModuleBase::matrix& forcecc, } else { - if (PARAM.inp.nspin == 4) { - ucell_in.cal_ux(); -} + if (PARAM.inp.nspin == 4) + { + elecstate::cal_ux(ucell_in); + } const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &ucell_in); // etxc = std::get<0>(etxc_vtxc_v); diff --git a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp index 26d5e66c37..898657febc 100644 --- a/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp +++ b/source/module_hamilt_pw/hamilt_pwdft/stress_func_cc.cpp @@ -4,6 +4,7 @@ #include "module_base/math_integral.h" #include "module_base/timer.h" #include "module_hamilt_pw/hamilt_pwdft/global.h" +#include "module_elecstate/cal_ux.h" #ifdef USE_LIBXC #include "module_hamilt_general/module_xc/xc_functional_libxc.h" @@ -63,8 +64,10 @@ void Stress_Func::stress_cc(ModuleBase::matrix& sigma, } else { - if(PARAM.inp.nspin==4) { GlobalC::ucell.cal_ux(); -} + if(PARAM.inp.nspin==4) + { + elecstate::cal_ux(GlobalC::ucell); + } const auto etxc_vtxc_v = XC_Functional::v_xc(rho_basis->nrxx, chr, &GlobalC::ucell); // etxc = std::get<0>(etxc_vtxc_v); // may delete? // vtxc = std::get<1>(etxc_vtxc_v); // may delete?