diff --git a/source/module_elecstate/cal_ux.cpp b/source/module_elecstate/cal_ux.cpp index 441aacceae..f8e4b1c82a 100644 --- a/source/module_elecstate/cal_ux.cpp +++ b/source/module_elecstate/cal_ux.cpp @@ -4,62 +4,91 @@ namespace elecstate { void cal_ux(UnitCell& ucell) { + if (PARAM.inp.nspin != 4) { return; } - double amag, uxmod; + const double absolute_mag_thr = 1.0e-6; + + double amag = 0.0; + double uxmod = 0.0; + int starting_it = 0; int starting_ia = 0; - bool is_paraller; + bool is_paraller = false; + // 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++) { + for (int it = 0; it < ucell.ntype; it++) + { + for (int ia = 0; ia < ucell.atoms[it].na; ia++) + { + // m_loc_: local magnetization vector for each atom 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; - } - } + + // find the first atom (it,ia) whose magnetism is not zero + // compute ux + if (amag > absolute_mag_thr) + { + 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 any atom has magnetism, then break the for iteration + 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: " - //< starting_it || ia > starting_ia) + { + ucell.magnet.lsign_ + = ucell.magnet.lsign_ + && judge_parallel(ucell.magnet.ux_, ucell.atoms[it].m_loc_[ia]); + } + } + } + + // if all of the atoms have the same parallel magnetism direction, + // then set the direction to a unit vector + 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 < absolute_mag_thr) + { + ModuleBase::WARNING_QUIT("elecstate::cal_ux", "wrong uxmod"); + } + + // reset the magnetism for each direction + for (int i = 0; i < 3; i++) + { + ucell.magnet.ux_[i] *= 1 / sqrt(uxmod); + } + } + return; } bool judge_parallel(double a[3], ModuleBase::Vector3 b) { @@ -71,4 +100,4 @@ bool judge_parallel(double a[3], ModuleBase::Vector3 b) { jp = (fabs(cross) < 1e-6); return jp; } -} \ No newline at end of file +} diff --git a/source/module_elecstate/magnetism.cpp b/source/module_elecstate/magnetism.cpp index b41f5261ad..4cf77c3ff0 100644 --- a/source/module_elecstate/magnetism.cpp +++ b/source/module_elecstate/magnetism.cpp @@ -22,6 +22,8 @@ void Magnetism::compute_magnetization(const double& omega, const double* const * rho, double* nelec_spin) { + assert(nxyz>0); + if (PARAM.inp.nspin==2) { this->tot_magnetization = 0.00; @@ -57,36 +59,53 @@ void Magnetism::compute_magnetization(const double& omega, // noncolliear : else if(PARAM.inp.nspin==4) { - for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] = 0.00; -} + for(int i=0;i<3;i++) + { + this->tot_magnetization_nc[i] = 0.00; + } + this->abs_magnetization = 0.00; for (int ir=0; irtot_magnetization_nc[i] += rho[i+1][ir]; -} + for(int i=0;i<3;i++) + { + this->tot_magnetization_nc[i] += rho[i+1][ir]; + } this->abs_magnetization += std::abs(diff); } #ifdef __MPI Parallel_Reduce::reduce_pool(this->tot_magnetization_nc, 3); Parallel_Reduce::reduce_pool(this->abs_magnetization); #endif - for(int i=0;i<3;i++) {this->tot_magnetization_nc[i] *= omega/ nxyz; -} + for(int i=0;i<3;i++) + { + this->tot_magnetization_nc[i] *= omega/ nxyz; + } + this->abs_magnetization *= omega/ nxyz; - GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"<<'\t'<tot_magnetization_nc[0]<<'\t'<tot_magnetization_nc[1]<<'\t'<tot_magnetization_nc[2]<<'\n'; + + GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)" + <<'\t'<tot_magnetization_nc[0] + <<'\t'<tot_magnetization_nc[1] + <<'\t'<tot_magnetization_nc[2]<<'\n'; + ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization); } return; } + bool Magnetism::judge_parallel(double a[3], ModuleBase::Vector3 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); + double cross=0.0; + 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; } diff --git a/source/module_elecstate/magnetism.h b/source/module_elecstate/magnetism.h index 408ba9d039..8a841c7006 100644 --- a/source/module_elecstate/magnetism.h +++ b/source/module_elecstate/magnetism.h @@ -13,7 +13,7 @@ class Magnetism Magnetism(); ~Magnetism(); - // notice : becast is done in unitcell + // notice : bcast (MPI operation) is done in unitcell double *start_magnetization=nullptr; // tot_magnetization : majority spin - minority spin (nelup - neldw). @@ -28,10 +28,9 @@ class Magnetism double* nelec_spin = nullptr); - ModuleBase::Vector3 *m_loc_=nullptr;//magnetization for each element along c-axis - double *angle1_=nullptr; //angle between c-axis and real spin std::vector - double *angle2_=nullptr; //angle between a-axis and real spin std::vector projection in ab-plane - //void cal_ux(const int ntype); + ModuleBase::Vector3 *m_loc_=nullptr; //magnetization for each element along c-axis + double *angle1_=nullptr; //angle between c-axis and real spin std::vector + double *angle2_=nullptr; //angle between a-axis and real spin std::vector projection in ab-plane double ux_[3]={0.0}; bool lsign_=false; diff --git a/source/module_esolver/esolver_fp.cpp b/source/module_esolver/esolver_fp.cpp index c92c2ebd34..1f31453ebd 100644 --- a/source/module_esolver/esolver_fp.cpp +++ b/source/module_esolver/esolver_fp.cpp @@ -258,6 +258,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) { ModuleBase::TITLE("ESolver_FP", "before_scf"); + // if the cell has changed if (ucell.cell_parameter_updated) { // only G-vector and K-vector are changed due to the change of lattice @@ -266,6 +267,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pw_rho->collect_local_pw(); this->pw_rho->collect_uniqgg(); + // if double grid used in USPP, update related quantities in dense grid if (PARAM.globalv.double_grid) { this->pw_rhod->initgrids(ucell.lat0, ucell.latvec, pw_rhod->nx, pw_rhod->ny, pw_rhod->nz); @@ -273,35 +275,40 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pw_rhod->collect_uniqgg(); } + // reset local pseudopotentials this->locpp.init_vloc(ucell, this->pw_rhod); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL"); this->pelec->omega = ucell.omega; + // perform symmetry analysis if (ModuleSymmetry::Symmetry::symm_flag == 1) { ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY"); } + // reset k-points kv.set_after_vc(PARAM.inp.nspin, ucell.G, ucell.latvec); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS"); } + //---------------------------------------------------------- // charge extrapolation + //---------------------------------------------------------- if (ucell.ionic_position_updated) { this->CE.update_all_dis(ucell); - this->CE.extrapolate_charge(&(this->Pgrid), + this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, - &(this->sf), + &this->sf, GlobalV::ofs_running, GlobalV::ofs_warning); } //---------------------------------------------------------- - // about vdw, jiyy add vdwd3 and linpz add vdwd2 + //! calculate D2 or D3 vdW //---------------------------------------------------------- auto vdw_solver = vdw::make_vdw(ucell, PARAM.inp, &(GlobalV::ofs_running)); if (vdw_solver != nullptr) @@ -309,19 +316,22 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) this->pelec->f_en.evdw = vdw_solver->get_energy(); } - // calculate ewald energy + //---------------------------------------------------------- + //! calculate ewald energy + //---------------------------------------------------------- if (!PARAM.inp.test_skip_ewald) { this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac); } //---------------------------------------------------------- - //! cal_ux should be called before init_scf because - //! the direction of ux is used in noncoline_rho + //! set direction of magnetism, used in non-collinear case //---------------------------------------------------------- elecstate::cal_ux(ucell); + //---------------------------------------------------------- //! output the initial charge density + //---------------------------------------------------------- if (PARAM.inp.out_chg[0] == 2) { for (int is = 0; is < PARAM.inp.nspin; is++) @@ -339,7 +349,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep) } } + //---------------------------------------------------------- //! output total local potential of the initial charge density + //---------------------------------------------------------- if (PARAM.inp.out_pot == 3) { for (int is = 0; is < PARAM.inp.nspin; is++) diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp index b123a354fe..1e6cc8d827 100644 --- a/source/module_esolver/esolver_ks_lcao.cpp +++ b/source/module_esolver/esolver_ks_lcao.cpp @@ -160,7 +160,9 @@ void ESolver_KS_LCAO::before_all_runners(UnitCell& ucell, const Input_pa #ifdef __EXX // 5) initialize exx // PLEASE simplify the Exx_Global interface - if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "relax" || PARAM.inp.calculation == "cell-relax" + if (PARAM.inp.calculation == "scf" + || PARAM.inp.calculation == "relax" + || PARAM.inp.calculation == "cell-relax" || PARAM.inp.calculation == "md") { if (GlobalC::exx_info.info_global.cal_exx) diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp index 5426ce2294..dc7b522c43 100644 --- a/source/module_esolver/lcao_others.cpp +++ b/source/module_esolver/lcao_others.cpp @@ -21,15 +21,16 @@ #include "module_base/formatter.h" #include "module_elecstate/elecstate_lcao.h" #include "module_elecstate/module_dm/cal_dm_psi.h" -#include "module_hamilt_general/module_ewald/H_Ewald_pw.h" -#include "module_hamilt_general/module_vdw/vdw.h" + #include "module_hamilt_lcao/hamilt_lcaodft/LCAO_domain.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/op_exx_lcao.h" #include "module_hamilt_lcao/hamilt_lcaodft/operator_lcao/operator_lcao.h" #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h" + #include "module_io/read_wfc_nao.h" #include "module_io/write_elecstat_pot.h" #include "module_io/write_wfc_nao.h" + #ifdef __EXX #include "module_io/restart_exx_csr.h" #endif @@ -263,6 +264,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) this->psi, this->pelec); } + //========================================================= // cal_ux should be called before init_scf because // the direction of ux is used in noncoline_rho @@ -271,14 +273,15 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) // pelec should be initialized before these calculations this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm); + // self consistent calculations for electronic ground state if (cal_type == "get_pchg") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting partial charge"); - IState_Charge ISC(this->psi, &(this->pv)); + IState_Charge chr_i(this->psi, &(this->pv)); if (PARAM.globalv.gamma_only_local) { - ISC.begin(this->GG, + chr_i.begin(this->GG, this->chr.rho, this->pelec->wg, this->pelec->eferm.get_all_ef(), @@ -306,7 +309,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - ISC.begin(this->GK, + chr_i.begin(this->GK, this->chr.rho, this->chr.rhog, this->pelec->wg, @@ -342,10 +345,10 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) else if (cal_type == "get_wf") { std::cout << FmtCore::format("\n * * * * * *\n << Start %s.\n", "getting wave function"); - IState_Envelope IEP(this->pelec); + IState_Envelope wavefunc(this->pelec); if (PARAM.globalv.gamma_only_local) { - IEP.begin(ucell, + wavefunc.begin(ucell, this->psi, this->pw_rhod, this->pw_wfc, @@ -367,7 +370,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - IEP.begin(ucell, + wavefunc.begin(ucell, this->psi, this->pw_rhod, this->pw_wfc, @@ -391,7 +394,7 @@ void ESolver_KS_LCAO::others(UnitCell& ucell, const int istep) } else { - ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); + ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::others", "CALCULATION type not supported"); } ModuleBase::timer::tick("ESolver_KS_LCAO", "others"); diff --git a/source/module_esolver/pw_others.cpp b/source/module_esolver/pw_others.cpp index c36ddae6e7..5461bed0a0 100644 --- a/source/module_esolver/pw_others.cpp +++ b/source/module_esolver/pw_others.cpp @@ -58,7 +58,8 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) const std::string cal_type = PARAM.inp.calculation; - if (cal_type == "test_memory") { + if (cal_type == "test_memory") + { Cal_Test::test_memory(ucell.nat, ucell.ntype, ucell.GGT, @@ -66,7 +67,9 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) this->pw_wfc, this->p_chgmix->get_mixing_mode(), this->p_chgmix->get_mixing_ndim()); - } else if (cal_type == "gen_bessel") { + } + else if (cal_type == "gen_bessel") + { Numerical_Descriptor nc; nc.output_descriptor(ucell, this->psi[0], @@ -75,7 +78,9 @@ void ESolver_KS_PW::others(UnitCell& ucell, const int istep) PARAM.inp.bessel_descriptor_tolerence, this->kv.get_nks()); ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "GENERATE DESCRIPTOR FOR DEEPKS"); - } else { + } + else + { ModuleBase::WARNING_QUIT("ESolver_KS_PW::others", "CALCULATION type not supported"); }