diff --git a/source/source_io/module_ctrl/ctrl_output_pw.cpp b/source/source_io/module_ctrl/ctrl_output_pw.cpp index 58c88461641..ac525789c3e 100644 --- a/source/source_io/module_ctrl/ctrl_output_pw.cpp +++ b/source/source_io/module_ctrl/ctrl_output_pw.cpp @@ -227,7 +227,8 @@ void ModuleIO::ctrl_scf_pw(const int istep, { // float type has not been implemented auto* onsite_p = projectors::OnsiteProjector::get_instance(); onsite_p->cal_occupations(reinterpret_cast, Device>*>(stp.template get_psi_t()), - pelec->wg); + pelec->wg, + inp.nspin); } ModuleBase::timer::end("ModuleIO", "ctrl_scf_pw"); diff --git a/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp b/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp index 307dd312103..0aa0afa3de4 100644 --- a/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp +++ b/source/source_lcao/module_deltaspin/lambda_loop_helper.cpp @@ -353,6 +353,14 @@ void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) const int nspin = sc.get_nspin(); const auto& Mi = sc.get_Mi(); const auto& atomLabel = sc.get_atomLabels(); + std::vector atom_labels_iat(nat); + for (const auto& it : sc.get_atomCounts()) + { + for (int ia = 0; ia < it.second; ++ia) + { + atom_labels_iat[sc.get_iat(it.first, ia)] = atomLabel[it.first] + std::to_string(ia + 1); + } + } std::vector mag_x(nat, 0.0); std::vector mag_y(nat, 0.0); std::vector mag_z(nat, 0.0); @@ -369,7 +377,7 @@ void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) { mag_z[iat] = Mi[iat].z; } - table << atomLabel << mag_z; + table << atom_labels_iat << mag_z; ofs_running << table.str() << std::endl; } else if (nspin == 4) @@ -387,7 +395,7 @@ void print_Mi(const SpinConstrain& sc, std::ostream& ofs_running) mag_y[iat] = Mi[iat].y; mag_z[iat] = Mi[iat].z; } - table << atomLabel << mag_x << mag_y << mag_z; + table << atom_labels_iat << mag_x << mag_y << mag_z; ofs_running << table.str() << std::endl; } } @@ -406,6 +414,14 @@ void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) const int nspin = sc.get_nspin(); const auto& lambda = sc.get_sc_lambda(); const auto& atomLabel = sc.get_atomLabels(); + std::vector atom_labels_iat(nat); + for (const auto& it : sc.get_atomCounts()) + { + for (int ia = 0; ia < it.second; ++ia) + { + atom_labels_iat[sc.get_iat(it.first, ia)] = atomLabel[it.first] + std::to_string(ia + 1); + } + } std::vector mag_force_x(nat, 0.0); std::vector mag_force_y(nat, 0.0); std::vector mag_force_z(nat, 0.0); @@ -422,7 +438,7 @@ void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) { mag_force_z[iat] = lambda[iat].z * ModuleBase::Ry_to_eV; } - table << atomLabel << mag_force_z; + table << atom_labels_iat << mag_force_z; ofs_running << table.str() << std::endl; } else if (nspin == 4) @@ -440,7 +456,7 @@ void print_Mag_Force(const SpinConstrain& sc, std::ostream& ofs_running) mag_force_y[iat] = lambda[iat].y * ModuleBase::Ry_to_eV; mag_force_z[iat] = lambda[iat].z * ModuleBase::Ry_to_eV; } - table << atomLabel << mag_force_x << mag_force_y << mag_force_z; + table << atom_labels_iat << mag_force_x << mag_force_y << mag_force_z; ofs_running << table.str() << std::endl; } } diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h index 1af2c116273..343cef7bbe3 100644 --- a/source/source_lcao/module_deltaspin/spin_constrain.h +++ b/source/source_lcao/module_deltaspin/spin_constrain.h @@ -361,7 +361,7 @@ class SpinConstrain /// check atomCounts void check_atomCounts() const { state_.check_atomCounts(); } /// get iat - int get_iat(int itype, int atom_index) { return state_.get_iat(itype, atom_index); } + int get_iat(int itype, int atom_index) const { return state_.get_iat(itype, atom_index); } /// set nspin void set_nspin(int nspin) { state_.set_nspin(nspin); } /// get nspin diff --git a/source/source_lcao/module_operator_lcao/operator_lcao.cpp b/source/source_lcao/module_operator_lcao/operator_lcao.cpp index f0560e69492..c1829b25d10 100644 --- a/source/source_lcao/module_operator_lcao/operator_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/operator_lcao.cpp @@ -78,6 +78,14 @@ void OperatorLCAO::init(const int ik_in) { this->hR->set_zero(); } } + // propagate current_spin to next operator so all nodes in the chain + // share the same spin state, set by HamiltLCAO::updateHk via set_current_spin. + // This is done before processing this node so that children receive the + // correct spin regardless of any local toggling that may happen inside + // this operator's contributeHR(). + if (this->next_op != nullptr) { + dynamic_cast*>(this->next_op)->current_spin = this->current_spin; + } switch (this->cal_type) { case calculation_type::lcao_overlap: { // cal_type=lcao_overlap refer to overlap matrix operators, which are diff --git a/source/source_pw/module_pwdft/onsite_proj.h b/source/source_pw/module_pwdft/onsite_proj.h index a4181e4b4bb..0ede3a6e8c0 100644 --- a/source/source_pw/module_pwdft/onsite_proj.h +++ b/source/source_pw/module_pwdft/onsite_proj.h @@ -75,7 +75,9 @@ namespace projectors const ModuleBase::matrix& ekb); /// @brief calculate and print the occupations of all lm orbitals - void cal_occupations(const psi::Psi, Device>* psi, const ModuleBase::matrix& wg_in); + void cal_occupations(const psi::Psi, Device>* psi, + const ModuleBase::matrix& wg_in, + const int nspin_in); int get_size_becp() const { return size_becp; } std::complex* get_becp() const { return becp; } diff --git a/source/source_pw/module_pwdft/onsite_proj_overlap.cpp b/source/source_pw/module_pwdft/onsite_proj_overlap.cpp index 989a6092148..aa2e280a0e9 100644 --- a/source/source_pw/module_pwdft/onsite_proj_overlap.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_overlap.cpp @@ -143,7 +143,8 @@ void projectors::OnsiteProjector::overlap_proj_psi( template void projectors::OnsiteProjector::cal_occupations( const psi::Psi, Device>* psi_in, - const ModuleBase::matrix& wg_in) + const ModuleBase::matrix& wg_in, + const int nspin_in) { ModuleBase::timer::start("OnsiteProj", "cal_occupation"); this->tabulate_atomic(0); @@ -151,6 +152,7 @@ void projectors::OnsiteProjector::cal_occupations( // loop over k-points to calculate Mi of \sum_{k,i,l,m} const int nbands = psi_in->get_nbands(); + const int npol = psi_in->get_npol(); for(int ik = 0; ik < psi_in->get_nk(); ik++) { psi_in->fix_k(ik); @@ -159,16 +161,19 @@ void projectors::OnsiteProjector::cal_occupations( this->tabulate_atomic(ik); } // std::cout << __FILE__ << ":" << __LINE__ << " nbands = " << nbands << std::endl; - this->overlap_proj_psi( - nbands * psi_in->get_npol(), - psi_in->get_pointer()); + this->overlap_proj_psi(nbands * npol, psi_in->get_pointer()); const std::complex* becp_p = this->get_h_becp(); // becp(nbands*npol , nkb) // mag = wg * \sum_{nh}becp * becp int nkb = this->tot_nproj; - //nkb = 18; - //std::cout << "at " << __FILE__ << ": " << __LINE__ << " output nbands: " << nbands << std::endl; - //std::cout << "at " << __FILE__ << ": " << __LINE__ << " output nkb: " << nkb << std::endl; + // nspin=2 (npol=1): the spin-up and spin-down channels are separate + // k-points. Store spin-up occupancy in the up-up Pauli block (occ[0]) + // and spin-down occupancy in the down-down block (occ[3]) so that + // print_orb_chg() yields: + // Charge = occ[0] + occ[3], Mag(z) = occ[0] - occ[3] + // nspin=1 (npol=1): no spin polarization, split the occupancy evenly + // between occ[0] and occ[3] so that the printed magnetization is zero. + // nspin=4 (npol=2): both spinor components are interleaved per band. for(int ib = 0;ib::cal_occupations( for(int ih = 0; ih < nh; ih++) { const int occ_index = (begin_ih + ih) * 4; - const int index = ib*2*nkb + begin_ih + ih; - occs[occ_index] += weight * conj(becp_p[index]) * becp_p[index]; - occs[occ_index + 1] += weight * conj(becp_p[index]) * becp_p[index + nkb]; - occs[occ_index + 2] += weight * conj(becp_p[index + nkb]) * becp_p[index]; - occs[occ_index + 3] += weight * conj(becp_p[index + nkb]) * becp_p[index + nkb]; + if (npol == 1) + { + const int index = ib * nkb + begin_ih + ih; + const double occ = weight * (conj(becp_p[index]) * becp_p[index]).real(); + if (nspin_in == 2 && this->isk_ && this->isk_[ik] == 1) + { + occs[occ_index + 3] += occ; + } + else if (nspin_in == 1) + { + occs[occ_index] += 0.5 * occ; + occs[occ_index + 3] += 0.5 * occ; + } + else + { + occs[occ_index] += occ; + } + } + else + { + const int index = ib * 2 * nkb + begin_ih + ih; + occs[occ_index] += weight * conj(becp_p[index]) * becp_p[index]; + occs[occ_index + 1] += weight * conj(becp_p[index]) * becp_p[index + nkb]; + occs[occ_index + 2] += weight * conj(becp_p[index + nkb]) * becp_p[index]; + occs[occ_index + 3] += weight * conj(becp_p[index + nkb]) * becp_p[index + nkb]; + } } begin_ih += nh; } @@ -213,7 +239,8 @@ void projectors::OnsiteProjector::overlap_proj_ template void projectors::OnsiteProjector::cal_occupations( const psi::Psi, base_device::DEVICE_CPU>*, - const ModuleBase::matrix&); + const ModuleBase::matrix&, + const int); #if ((defined __CUDA) || (defined __ROCM)) template @@ -227,5 +254,6 @@ void projectors::OnsiteProjector::overlap_proj_ template void projectors::OnsiteProjector::cal_occupations( const psi::Psi, base_device::DEVICE_GPU>*, - const ModuleBase::matrix&); + const ModuleBase::matrix&, + const int); #endif diff --git a/source/source_pw/module_pwdft/onsite_proj_print.cpp b/source/source_pw/module_pwdft/onsite_proj_print.cpp index adb6e1246c0..6e97d133a81 100644 --- a/source/source_pw/module_pwdft/onsite_proj_print.cpp +++ b/source/source_pw/module_pwdft/onsite_proj_print.cpp @@ -32,11 +32,13 @@ void print_orb_chg( /*align=*/{/*value*/FmtTable::Align::RIGHT, /*title*/FmtTable::Align::LEFT}); // parameters for mag output int occ_index = 0; + std::vector atom_labels_iat(ucell->nat); for(int iat=0; iatnat; iat++) { const int it = ucell->iat2it[iat]; std::string atom_label = atom_labels[it]; int ia = ucell->iat2ia[iat]; + atom_labels_iat[iat] = atom_label + std::to_string(ia+1); GlobalV::ofs_running << FmtCore::format("%-20s", atom_label+std::to_string(ia+1)) << std::endl; std::vector sum(4, 0.0); int current_l = 1; @@ -72,7 +74,7 @@ void print_orb_chg( GlobalV::ofs_running << std::endl; // Print magnetism table - print_mag_table(atom_labels, mag_x, mag_y, mag_z); + print_mag_table(atom_labels_iat, mag_x, mag_y, mag_z); } void print_mag_table( diff --git a/tests/03_NAO_multik/scf_deltaspin2/result.ref b/tests/03_NAO_multik/scf_deltaspin2/result.ref index fc858235e90..ea1d55e48c8 100644 --- a/tests/03_NAO_multik/scf_deltaspin2/result.ref +++ b/tests/03_NAO_multik/scf_deltaspin2/result.ref @@ -1,4 +1,5 @@ -etotref -6762.435776188675 -etotperatomref -3381.217888094338 -totalforceref 63.230574 -totalstressref 2916.957427 +etotref -6789.4013179079911424 +etotperatomref -3394.7006589540 +totalforceref 67.888296 +totalstressref 6225.500296 +totaltimeref 2.43 diff --git a/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref b/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref index 9dd83d8e706..dd48e508258 100644 --- a/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref +++ b/tests/17_DS_DFTU/12_PW_DS_S2_Z/result.ref @@ -1,3 +1,3 @@ -etotref -6366.569118260046 -etotperatomref -3183.2845591300 -totaltimeref 1.97 +etotref -6369.198268154857 +etotperatomref -3184.5991340774 +totaltimeref 4.25 diff --git a/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref b/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref index 654cb15b3ee..72fb5c43b87 100644 --- a/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref +++ b/tests/17_DS_DFTU/18_PW_DFTU_DS_S2_Z/result.ref @@ -1,3 +1,3 @@ -etotref -6355.9855588350255857 -etotperatomref -3177.9927794175 -totaltimeref 2.89 +etotref -6360.5555597255606699 +etotperatomref -3180.2777798628 +totaltimeref 5.86 diff --git a/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref b/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref index e40be3f01bb..c3b83965fd3 100644 --- a/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref +++ b/tests/17_DS_DFTU/19_PW_DFTU_DS_S4_XY/result.ref @@ -1,3 +1,3 @@ -etotref -6360.5554588729937677 -etotperatomref -3180.277729436497 -totaltimeref 1.0 +etotref -6360.5555339529537378 +etotperatomref -3180.2777669765 +totaltimeref 9.56 diff --git a/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref b/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref index b06c7ee03a6..ae5a5a853dd 100644 --- a/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref +++ b/tests/17_DS_DFTU/21_PW_DFTU_DS_S4_Z/result.ref @@ -1,3 +1,3 @@ -etotref -6360.5554655414534864 -etotperatomref -3180.2777327707267 -totaltimeref 1.0 +etotref -6360.5554079455796455 +etotperatomref -3180.2777039728 +totaltimeref 6.65 diff --git a/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref b/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref index 8d127ea7298..eb51c776c6e 100644 --- a/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref +++ b/tests/17_DS_DFTU/41_PW_DS_S4_Thr10_XY/result.ref @@ -1,3 +1,3 @@ -etotref -6369.198254647004 -etotperatomref -3184.599127323502 -totaltimeref 1.0 +etotref -6369.198255055176 +etotperatomref -3184.5991275276 +totaltimeref 9.00