fix: correct Pauli-to-spinor conversion in DFT+U and DeltaSpin for nspin=4

Fix three critical bugs in non-collinear (nspin=4) LCAO calculations:

1. DFT+U transfer_vu (dftu_lcao.cpp): Fix sign error in Pauli-to-spinor
   conversion. The off-diagonal elements had wrong imaginary part sign:
   - Before: V_{up,down} = 0.5*(V_x + i*V_y)  (wrong)
   - After:  V_{up,down} = 0.5*(V_x - i*V_y)  (correct, from sigma_y)

2. DFT+U force/stress (dftu_force_stress.hpp): Convert VU from Pauli basis
   to spinor basis before force calculation. The old code incorrectly mixed
   Pauli-basis VU with spinor-basis DM.

3. DeltaSpin force/stress (dspin_force_stress.hpp): Convert lambda from
   Pauli basis to spinor basis. The constraint force F = lambda·dM/dR
   requires proper Pauli-to-spinor conversion:
   - lambda_spinor = (lambda_z, lambda_x, lambda_x, -lambda_z)
   for (uu, ud, du, dd) components.

These fixes ensure consistent Pauli-to-spinor conversion across all modules:
- H construction (gint_common.cpp): already fixed
- DFT+U Hamiltonian: fixed in this commit
- DFT+U force/stress: fixed in this commit
- DeltaSpin force/stress: fixed in this commit

Verified by scf_u_spin4 test (nspin=4 + DFT+U): SCF converges correctly.

Author: dyzheng

source/source_lcao/module_operator_lcao/dftu_force_stress.hpp

@@ -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)
             }
         }
 

source/source_lcao/module_operator_lcao/dftu_lcao.cpp

@@ -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());
@@ -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]]);
         }
     }
 }

source/source_lcao/module_operator_lcao/dspin_force_stress.hpp

@@ -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();
         }
     }
 }
@@ -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();
         }
     }
 }