Fix(PW): correct DFT+U force/stress spinor formula for npol=2

The vu array stores coefficients in order [H_↑↑, H_↓↑, H_↑↓, H_↓↓],
different from deeq_nc's [H_↑↑, H_↑↓, H_↓↑, H_↓↓]. The force/stress
formula must match: ps[2]*dbb1 + ps[1]*dbb2 for off-diagonal terms.

Fixed in CPU, CUDA, and ROCm kernels for both force and stress.

Author: dyzheng

source/source_pw/module_pwdft/kernels/cuda/force_op.cu

@@ -375,7 +375,7 @@ __global__ void cal_force_onsite(int wg_nc,
                 const thrust::complex<FPTYPE> dbb1 = conj(dbecp[inkb0]) * becp[inkb2 + nkb];
                 const thrust::complex<FPTYPE> dbb2 = conj(dbecp[inkb0 + nkb]) * becp[inkb2];
                 const thrust::complex<FPTYPE> dbb3 = conj(dbecp[inkb0 + nkb]) * becp[inkb2 + nkb];
-                const FPTYPE tmp = -fac * (ps[0] * dbb0 + ps[1] * dbb1 + ps[2] * dbb2 + ps[3] * dbb3).real();
+                const FPTYPE tmp = -fac * (ps[0] * dbb0 + ps[2] * dbb1 + ps[1] * dbb2 + ps[3] * dbb3).real();
                 atomicAdd(force + iat * forcenl_nc + ipol, tmp);
             }
         }

source/source_pw/module_pwdft/kernels/cuda/stress_op.cu

@@ -977,7 +977,7 @@ __global__ void cal_stress_onsite(
             const thrust::complex<FPTYPE> dbb1 = conj(dbecp[inkb1]) * becp[inkb2 + nkb];
             const thrust::complex<FPTYPE> dbb2 = conj(dbecp[inkb1 + nkb]) * becp[inkb2];
             const thrust::complex<FPTYPE> dbb3 = conj(dbecp[inkb1 + nkb]) * becp[inkb2 + nkb];
-            stress_var -= fac * (ps[0] * dbb0 + ps[1] * dbb1 + ps[2] * dbb2 + ps[3] * dbb3).real();
+            stress_var -= fac * (ps[0] * dbb0 + ps[2] * dbb1 + ps[1] * dbb2 + ps[3] * dbb3).real();
         }
         ++iat;
         sum+=nprojs;

source/source_pw/module_pwdft/kernels/force_op.cpp

@@ -352,7 +352,7 @@ struct cal_force_nl_op<FPTYPE, base_device::DEVICE_CPU>
                                     const std::complex<FPTYPE> dbb2 = conj(dbecp[index0 + nkb]) * becp[index1];
                                     const std::complex<FPTYPE> dbb3 = conj(dbecp[index0 + nkb]) * becp[index1 + nkb];
 
-                                    local_force[iforce] -= fac * (ps[0] * dbb0 + ps[1] * dbb1 + ps[2] * dbb2 + ps[3] * dbb3).real();
+                                    local_force[iforce] -= fac * (ps[0] * dbb0 + ps[2] * dbb1 + ps[1] * dbb2 + ps[3] * dbb3).real();
                                 }
                             }
                             else if(npol == 1)

source/source_pw/module_pwdft/kernels/rocm/force_op.hip.cu

@@ -361,7 +361,7 @@ __global__ void cal_force_onsite(int wg_nc,
                 const thrust::complex<FPTYPE> dbb1 = conj(dbecp[inkb0]) * becp[inkb2 + nkb];
                 const thrust::complex<FPTYPE> dbb2 = conj(dbecp[inkb0 + nkb]) * becp[inkb2];
                 const thrust::complex<FPTYPE> dbb3 = conj(dbecp[inkb0 + nkb]) * becp[inkb2 + nkb];
-                const FPTYPE tmp = -fac * (ps[0] * dbb0 + ps[1] * dbb1 + ps[2] * dbb2 + ps[3] * dbb3).real();
+                const FPTYPE tmp = -fac * (ps[0] * dbb0 + ps[2] * dbb1 + ps[1] * dbb2 + ps[3] * dbb3).real();
                 atomicAdd(force + iat * forcenl_nc + ipol, tmp);
             }
         }

source/source_pw/module_pwdft/kernels/rocm/stress_op.hip.cu

@@ -966,7 +966,7 @@ __global__ void cal_stress_onsite(
             const thrust::complex<FPTYPE> dbb1 = conj(dbecp[inkb1]) * becp[inkb2 + nkb];
             const thrust::complex<FPTYPE> dbb2 = conj(dbecp[inkb1 + nkb]) * becp[inkb2];
             const thrust::complex<FPTYPE> dbb3 = conj(dbecp[inkb1 + nkb]) * becp[inkb2 + nkb];
-            stress_var -= fac * (ps[0] * dbb0 + ps[1] * dbb1 + ps[2] * dbb2 + ps[3] * dbb3).real();
+            stress_var -= fac * (ps[0] * dbb0 + ps[2] * dbb1 + ps[1] * dbb2 + ps[3] * dbb3).real();
         }
         ++iat;
         sum+=nprojs;

source/source_pw/module_pwdft/kernels/stress_op.cpp

@@ -318,7 +318,7 @@ struct cal_stress_nl_op<FPTYPE, base_device::DEVICE_CPU>
                                 const std::complex<FPTYPE> dbb1 = conj(dbecp[inkb1]) * becp[nkb + inkb2];
                                 const std::complex<FPTYPE> dbb2 = conj(dbecp[nkb + inkb1]) * becp[inkb2];
                                 const std::complex<FPTYPE> dbb3 = conj(dbecp[nkb + inkb1]) * becp[nkb + inkb2];
-                                local_stress -= fac * (ps[0] * dbb0 + ps[1] * dbb1 + ps[2] * dbb2 + ps[3] * dbb3).real();
+                                local_stress -= fac * (ps[0] * dbb0 + ps[2] * dbb1 + ps[1] * dbb2 + ps[3] * dbb3).real();
                             }
                         } // end ip
                         break;