Skip to content

Fix Pauli-to-Spinor Conversion in LCAO Non-Collinear Calculations #7513

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dyzheng wants to merge 3 commits into
base: develop
Choose a base branch
from
Open

Fix Pauli-to-Spinor Conversion in LCAO Non-Collinear Calculations #7513

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
eb94066
fix: correct Pauli-to-spinor Hamiltonian conversion for nspin=4
dyzheng Jun 23, 2026
4e5b62a
fix: correct Pauli-to-spinor conversion in DFT+U and DeltaSpin for ns…
dyzheng Jun 23, 2026
4738c6a
chore: remove .py and .md test files from PR
dyzheng Jun 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 13 additions & 5 deletions source/source_lcao/module_gint/gint_common.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,8 +168,15 @@ void merge_hr_part_to_hR(const std::vector<hamilt::HContainer<double>>& hr_gint_
std::vector<int> row_set = {0, 0, 1, 1};
std::vector<int> 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<int> clx_i = {1, 0, 0, -1};
std::vector<int> clx_j = {0, 1, -1, 0};
std::vector<int> 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();
Expand Down Expand Up @@ -203,17 +210,18 @@ void merge_hr_part_to_hR(const std::vector<hamilt::HContainer<double>>& hr_gint_
+ std::complex<double>(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);
for (int irow = 0; irow < upper_mat->get_row_size(); ++irow)
{
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));
}
}
}
Expand Down
20 changes: 18 additions & 2 deletions source/source_lcao/module_operator_lcao/dftu_force_stress.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,9 +148,25 @@ void DFTU<OperatorLCAO<TK, TR>>::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<double> (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<double> 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)
}
}

Expand Down
14 changes: 10 additions & 4 deletions source/source_lcao/module_operator_lcao/dftu_lcao.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -612,7 +612,14 @@ void hamilt::DFTU<hamilt::OperatorLCAO<std::complex<double>, std::complex<double
assert(vu.size() == vu_tmp.size());
#endif

// TR == std::complex<double> transfer from double to std::complex<double>
// 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());
Expand All @@ -627,9 +634,8 @@ void hamilt::DFTU<hamilt::OperatorLCAO<std::complex<double>, std::complex<double
index[3] = m2 * m_size + m1 + m_size2 * 3;
vu[index[0]] = 0.5 * (vu_tmp[index[0]] + vu_tmp[index[3]]);
vu[index[3]] = 0.5 * (vu_tmp[index[0]] - vu_tmp[index[3]]);
// vu should be std::complex<double> type, but here we use double type for test
vu[index[1]] = 0.5 * (vu_tmp[index[1]] + std::complex<double>(0.0, 1.0) * vu_tmp[index[2]]);
vu[index[2]] = 0.5 * (vu_tmp[index[1]] - std::complex<double>(0.0, 1.0) * vu_tmp[index[2]]);
vu[index[1]] = 0.5 * (vu_tmp[index[1]] - std::complex<double>(0.0, 1.0) * vu_tmp[index[2]]);
vu[index[2]] = 0.5 * (vu_tmp[index[1]] + std::complex<double>(0.0, 1.0) * vu_tmp[index[2]]);
}
}
}
Expand Down
224 changes: 165 additions & 59 deletions source/source_lcao/module_operator_lcao/dspin_force_stress.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -274,43 +274,95 @@ void DeltaSpin<OperatorLCAO<TK, TR>>::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<double>& 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<double>& nlm2 = nlm2_all.find(col_indexes[iw2l])->second;
const std::vector<double>& nlm1 = nlm1_all.find(row_indexes[iw1l])->second;
for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol)
{
const std::vector<double>& 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; l<lmax; l++)
const int length = nlm1.size() / 4;
const int lmax = sqrt(length);
int index = 0;
for(int l = 0; l<lmax; l++)
{
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 - 1) * col_indexes.size();
}
}
}
else
{
// nspin=1 or nspin=2: original logic
for (int is = 1; is < nspin; is++)
{
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 std::vector<double>& 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<double>& 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<lmax; l++)
{
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 = - VU * <d phi_{I,R1}/d R1|chi_m> * <chi_m'|phi_{J,R2}>
// force2 = - VU * <phi_{I,R1}|d chi_m/d R0> * <chi_m'|phi_{J,R2>}
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();
}
}
}
Expand Down Expand Up @@ -345,48 +397,102 @@ void DeltaSpin<OperatorLCAO<TK, TR>>::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<double>& 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<double>& nlm2 = nlm2_all.find(col_indexes[iw2l])->second;
const std::vector<double>& nlm1 = nlm1_all.find(row_indexes[iw1l])->second;
for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol)
{
const std::vector<double>& 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; l<lmax; l++)
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<lmax; l++)
{
for (int m = 0; m < 2*l+1; m++)
{
index = l*l + m;
stress[0]
+= tmp * (nlm1[index + length] * dis1.x * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.x);
stress[1]
+= tmp * (nlm1[index + length] * dis1.y * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.y);
stress[2]
+= tmp * (nlm1[index + length] * dis1.z * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.z);
stress[3] += tmp
* (nlm1[index + length * 2] * dis1.y * nlm2[index]
+ nlm1[index] * nlm2[index + length * 2] * dis2.y);
stress[4] += tmp
* (nlm1[index + length * 2] * dis1.z * nlm2[index]
+ nlm1[index] * nlm2[index + length * 2] * dis2.z);
stress[5] += tmp
* (nlm1[index + length * 3] * dis1.z * nlm2[index]
+ nlm1[index] * nlm2[index + length * 3] * dis2.z);
}
}
dm_pointer += npol;
}
dm_pointer += (npol - 1) * col_indexes.size();
}
}
}
else
{
// nspin=1 or nspin=2: original logic
for (int is = 1; is < nspin; is++)
{
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 std::vector<double>& 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<double>& 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<lmax; l++)
{
index = l*l + m;
stress[0]
+= tmp * (nlm1[index + length] * dis1.x * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.x);
stress[1]
+= tmp * (nlm1[index + length] * dis1.y * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.y);
stress[2]
+= tmp * (nlm1[index + length] * dis1.z * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.z);
stress[3] += tmp
* (nlm1[index + length * 2] * dis1.y * nlm2[index]
+ nlm1[index] * nlm2[index + length * 2] * dis2.y);
stress[4] += tmp
* (nlm1[index + length * 2] * dis1.z * nlm2[index]
+ nlm1[index] * nlm2[index + length * 2] * dis2.z);
stress[5] += tmp
* (nlm1[index + length * 3] * dis1.z * nlm2[index]
+ nlm1[index] * nlm2[index + length * 3] * dis2.z);
for (int m = 0; m < 2*l+1; m++)
{
index = l*l + m;
stress[0]
+= tmp * (nlm1[index + length] * dis1.x * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.x);
stress[1]
+= tmp * (nlm1[index + length] * dis1.y * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.y);
stress[2]
+= tmp * (nlm1[index + length] * dis1.z * nlm2[index] + nlm1[index] * nlm2[index + length] * dis2.z);
stress[3] += tmp
* (nlm1[index + length * 2] * dis1.y * nlm2[index]
+ nlm1[index] * nlm2[index + length * 2] * dis2.y);
stress[4] += tmp
* (nlm1[index + length * 2] * dis1.z * nlm2[index]
+ nlm1[index] * nlm2[index + length * 2] * dis2.z);
stress[5] += tmp
* (nlm1[index + length * 3] * dis1.z * nlm2[index]
+ nlm1[index] * nlm2[index + length * 3] * dis2.z);
}
}
dm_pointer += npol;
}
dm_pointer += npol;
dm_pointer += (npol - 1) * col_indexes.size();
}
dm_pointer += (npol - 1) * col_indexes.size();
}
}
}
Expand Down
Loading
Loading