diff --git a/source/source_lcao/module_gint/gint_common.cpp b/source/source_lcao/module_gint/gint_common.cpp index 987fcbe5a0..1de8969fd0 100644 --- a/source/source_lcao/module_gint/gint_common.cpp +++ b/source/source_lcao/module_gint/gint_common.cpp @@ -168,8 +168,15 @@ void merge_hr_part_to_hR(const std::vector>& hr_gint_ std::vector row_set = {0, 0, 1, 1}; std::vector col_set = {0, 1, 0, 1}; //construct complex matrix + // Pauli-to-spinor conversion: H = V_0*I + B_x*sigma_x + B_y*sigma_y + B_z*sigma_z + // sigma_y = [[0,-i],[i,0]], so H_{up,down} = B_x - i*B_y, H_{down,up} = B_x + i*B_y + // coefficient = clx_i + i*clx_j for each Pauli channel: + // is=0 (up,up): V_0 + B_z => coeff on B_z = +1 => clx_i=1, clx_j=0 + // is=1 (up,down): B_x - i*B_y => coeff on B_y = -i => clx_i=0, clx_j=-1 + // is=2 (down,up): B_x + i*B_y => coeff on B_y = +i => clx_i=0, clx_j=+1 + // is=3 (down,down): -(V_0 - B_z) => coeff on V_0 = -1 => clx_i=-1, clx_j=0 std::vector clx_i = {1, 0, 0, -1}; - std::vector clx_j = {0, 1, -1, 0}; + std::vector clx_j = {0, -1, 1, 0}; for (int is = 0; is < 4; is++){ if(!PARAM.globalv.domag && (is==1 || is==2)) continue; hR_tmp->set_zero(); @@ -203,9 +210,10 @@ void merge_hr_part_to_hR(const std::vector>& hr_gint_ + std::complex(clx_i[is], clx_j[is]) * mat_nspin2->get_value(irow, icol); } } - //fill the lower triangle matrix - //When is=0 or 3, the real part does not need conjugation; - //when is=1 or 2, the small matrix is not Hermitian, so conjugation is not needed + //fill the lower triangle matrix at -R by conjugate transpose of upper at R + // This ensures H(-R) = H(R)^dagger, required for Hermiticity of H(k). + // For real matrices (is=0,3), conj has no effect. + // For complex matrices (is=1,2), conj is essential. if (iat1 < iat2) { auto lower_mat = lower_ap->find_matrix(-R_index); @@ -213,7 +221,7 @@ void merge_hr_part_to_hR(const std::vector>& hr_gint_ { for (int icol = 0; icol < upper_mat->get_col_size(); ++icol) { - lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol); + lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol)); } } } diff --git a/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp b/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp index 38c96025fa..50b157c03e 100644 --- a/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp @@ -148,9 +148,25 @@ void DFTU>::cal_force_stress(const bool cal_force, this->cal_v_of_u(occ, tlp1, u_value, &VU[0], eu_tmp); if(this->nspin == 4) { - for (int i = 0; i < VU.size(); i++) + // For nspin=4, VU from cal_v_of_u is in Pauli basis: + // block 0 = V_0, block 1 = V_x, block 2 = V_y, block 3 = V_z + // The force formula is F = -Tr(VU * dDM/dR). + // Since DM is stored as BaseMatrix (real, spinor basis), + // we need VU in spinor basis as well. The conversion is: + // V_{up,up} = 0.5 * (V_0 + V_z) + // V_{down,down} = 0.5 * (V_0 - V_z) + // V_{up,down} = 0.5 * (V_x - i*V_y) → Re = 0.5*V_x, Im = -0.5*V_y + // V_{down,up} = 0.5 * (V_x + i*V_y) → Re = 0.5*V_x, Im = +0.5*V_y + // For real DM, only the real part of VU contributes to force. + // We convert VU in-place: block 0→V_uu, 1→V_ud(Re), 2→V_du(Re), 3→V_dd + const int m_size2 = tlp1 * tlp1; + std::vector VU_pauli = VU; // save Pauli-basis VU + for (int m = 0; m < m_size2; m++) { - VU[i] /= 2.0; + VU[m] = 0.5 * (VU_pauli[m] + VU_pauli[m + 3 * m_size2]); // V_uu = 0.5*(V_0+V_z) + VU[m + m_size2] = 0.5 * VU_pauli[m + m_size2]; // Re(V_ud) = 0.5*V_x + VU[m + 2 * m_size2] = 0.5 * VU_pauli[m + m_size2]; // Re(V_du) = 0.5*V_x + VU[m + 3 * m_size2] = 0.5 * (VU_pauli[m] - VU_pauli[m + 3 * m_size2]); // V_dd = 0.5*(V_0-V_z) } } diff --git a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp index 2b8e4f3506..92fac0ffab 100644 --- a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp @@ -612,7 +612,14 @@ void hamilt::DFTU, std::complex transfer from double to std::complex + // Pauli-to-spinor conversion for DFT+U potential: + // V = V_0*I + V_x*sigma_x + V_y*sigma_y + V_z*sigma_z + // sigma_y = [[0,-i],[i,0]], so: + // V_{up,up} = 0.5*(V_0 + V_z) + // V_{down,down} = 0.5*(V_0 - V_z) + // V_{up,down} = 0.5*(V_x - i*V_y) <- note: minus sign from sigma_y + // V_{down,up} = 0.5*(V_x + i*V_y) <- note: plus sign from sigma_y + // This is consistent with the convention in gint_common.cpp merge_hr_part_to_hR(). const int m_size = int(sqrt(vu.size()) / 2); const int m_size2 = m_size * m_size; vu.resize(vu_tmp.size()); @@ -627,9 +634,8 @@ void hamilt::DFTU, std::complex type, but here we use double type for test - vu[index[1]] = 0.5 * (vu_tmp[index[1]] + std::complex(0.0, 1.0) * vu_tmp[index[2]]); - vu[index[2]] = 0.5 * (vu_tmp[index[1]] - std::complex(0.0, 1.0) * vu_tmp[index[2]]); + vu[index[1]] = 0.5 * (vu_tmp[index[1]] - std::complex(0.0, 1.0) * vu_tmp[index[2]]); + vu[index[2]] = 0.5 * (vu_tmp[index[1]] + std::complex(0.0, 1.0) * vu_tmp[index[2]]); } } } diff --git a/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp b/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp index 955c6b4d26..0ca04bb76f 100644 --- a/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp @@ -274,43 +274,95 @@ void DeltaSpin>::cal_force_IJR(const int& iat1, } double tmp[3] = {0.0}; // calculate the local matrix - for (int is = 1; is < nspin; is++) + // For nspin=4, the constraint force is F = lambda · dM/dR where: + // Mx = DM_ud + DM_du, My = -i*(DM_ud - DM_du), Mz = DM_uu - DM_dd + // For real DM (dDM_ud = dDM_du), My contribution vanishes, so: + // F = lambda_x * (dDM_ud + dDM_du) + lambda_z * (dDM_uu - dDM_dd) + // We convert lambda from Pauli basis to spinor basis: + // lambda_uu = lambda_z, lambda_dd = -lambda_z + // lambda_ud = lambda_x, lambda_du = lambda_x + if (nspin == 4) { - const double lambda_tmp = nspin==2?lambda[2]:lambda[is-1]; - const double* dm_pointer = dmR_pointer->get_pointer(); - for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + // lambda in spinor basis: (lambda_uu, lambda_ud, lambda_du, lambda_dd) + const double lambda_spinor[4] = {lambda[2], lambda[0], lambda[0], -lambda[2]}; + for (int is = 0; is < 4; is++) { - const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; - for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + const double lambda_tmp = lambda_spinor[is]; + if (std::abs(lambda_tmp) < 1e-15) continue; + const double* dm_pointer = dmR_pointer->get_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) { - const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; #ifdef __DEBUG - assert(nlm1.size() == nlm2.size()); + assert(nlm1.size() == nlm2.size()); #endif - const int length = nlm1.size() / 4; - const int lmax = sqrt(length); - int index = 0; - for(int l = 0; lget_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) { - for (int m = 0; m < 2*l+1; m++) + const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; +#ifdef __DEBUG + assert(nlm1.size() == nlm2.size()); +#endif + const int length = nlm1.size() / 4; + const int lmax = sqrt(length); + int index = 0; + for(int l = 0; l * - // force2 = - VU * * } - force1[0] += tmp[0]; - force1[1] += tmp[1]; - force1[2] += tmp[2]; - force2[0] -= tmp[0]; - force2[1] -= tmp[1]; - force2[2] -= tmp[2]; + for (int m = 0; m < 2*l+1; m++) + { + index = l*l + m; + tmp[0] = lambda_tmp * nlm1[index + length] * nlm2[index] * dm_pointer[step_trace[is]]; + tmp[1] = lambda_tmp * nlm1[index + length * 2] * nlm2[index] * dm_pointer[step_trace[is]]; + tmp[2] = lambda_tmp * nlm1[index + length * 3] * nlm2[index] * dm_pointer[step_trace[is]]; + force1[0] += tmp[0]; + force1[1] += tmp[1]; + force1[2] += tmp[2]; + force2[0] -= tmp[0]; + force2[1] -= tmp[1]; + force2[2] -= tmp[2]; + } } + dm_pointer += npol; } - dm_pointer += npol; + dm_pointer += (npol - 1) * col_indexes.size(); } - dm_pointer += (npol - 1) * col_indexes.size(); } } } @@ -345,48 +397,102 @@ void DeltaSpin>::cal_stress_IJR(const int& iat1, step_trace[3] = col_indexes.size() + 1; } // calculate the local matrix - for (int is = 1; is < nspin; is++) + // For nspin=4, convert lambda from Pauli basis to spinor basis + if (nspin == 4) { - const double lambda_tmp = nspin==2?lambda[2]:lambda[is-1]; - const double* dm_pointer = dmR_pointer->get_pointer(); - for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + const double lambda_spinor[4] = {lambda[2], lambda[0], lambda[0], -lambda[2]}; + for (int is = 0; is < 4; is++) { - const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; - for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + const double lambda_tmp = lambda_spinor[is]; + if (std::abs(lambda_tmp) < 1e-15) continue; + const double* dm_pointer = dmR_pointer->get_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) { - const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; #ifdef __DEBUG - assert(nlm1.size() == nlm2.size()); + assert(nlm1.size() == nlm2.size()); #endif - const int length = nlm1.size() / 4; - const int lmax = sqrt(length); - double tmp = lambda_tmp * dm_pointer[step_trace[is]]; - int index = 0; - for(int l = 0; lget_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) { - for (int m = 0; m < 2*l+1; m++) + const std::vector& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; +#ifdef __DEBUG + assert(nlm1.size() == nlm2.size()); +#endif + const int length = nlm1.size() / 4; + const int lmax = sqrt(length); + double tmp = lambda_tmp * dm_pointer[step_trace[is]]; + int index = 0; + for(int l = 0; l