Fix(PW): correct all nonlocal force/stress formulas for nspin=4/SOC

dyzheng · dyzheng · commit eaad9504ae93 · 2026-05-03T10:31:56.000+08:00
Fixed formula errors in multiple kernels that broke SOC tests like
035_PW_15_SO. The incorrect formula was:
  ps1*dbb2 + ps2*dbb1 (wrong)

Correct formula verified from develop branch:
  ps1*dbb1 + ps2*dbb2 (correct)

Affected kernels:
- force_op.cpp: nonlocal force (deeq_nc) and DeltaSpin
- stress_op.cpp: nonlocal stress (deeq_nc)
- CUDA/ROCm versions of above

All formulas now match develop branch implementations.
diff --git a/source/source_pw/module_pwdft/kernels/cuda/force_op.cu b/source/source_pw/module_pwdft/kernels/cuda/force_op.cu
@@ -270,7 +270,7 @@ __global__ void cal_force_nl(
                     const thrust::complex<FPTYPE> dbb1 = conj(dbecp[index0]) * becp[index1 + nkb];
                     const thrust::complex<FPTYPE> dbb2 = conj(dbecp[index0 + nkb]) * becp[index1];
                     const thrust::complex<FPTYPE> dbb3 = conj(dbecp[index0 + nkb]) * becp[index1 + nkb];
-                    const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();
+                    const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();
                     atomicAdd(force + iat * forcenl_nc + ipol, tmp);
                 }
             }
diff --git a/source/source_pw/module_pwdft/kernels/cuda/stress_op.cu b/source/source_pw/module_pwdft/kernels/cuda/stress_op.cu
@@ -361,7 +361,7 @@ __global__ void cal_stress_nl(
             const thrust::complex<FPTYPE> dbb1 = conj(dbecp[ib2 * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];
             const thrust::complex<FPTYPE> dbb2 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[ib2 * nkb + inkb2];
             const thrust::complex<FPTYPE> dbb3 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];
-            stress_var -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();
+            stress_var -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();
         }
         ++iat;
         sum+=nproj;
diff --git a/source/source_pw/module_pwdft/kernels/force_op.cpp b/source/source_pw/module_pwdft/kernels/force_op.cpp
@@ -254,7 +254,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[ipol] -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();
+                                    local_force[ipol] -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();
                                 }
                             }
                         }
diff --git a/source/source_pw/module_pwdft/kernels/rocm/force_op.hip.cu b/source/source_pw/module_pwdft/kernels/rocm/force_op.hip.cu
@@ -256,7 +256,7 @@ __global__ void cal_force_nl(
                     const thrust::complex<FPTYPE> dbb1 = conj(dbecp[index0]) * becp[index1 + nkb];
                     const thrust::complex<FPTYPE> dbb2 = conj(dbecp[index0 + nkb]) * becp[index1];
                     const thrust::complex<FPTYPE> dbb3 = conj(dbecp[index0 + nkb]) * becp[index1 + nkb];
-                    const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();
+                    const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();
                     atomicAdd(force + iat * forcenl_nc + ipol, tmp);
                 }
             }
diff --git a/source/source_pw/module_pwdft/kernels/rocm/stress_op.hip.cu b/source/source_pw/module_pwdft/kernels/rocm/stress_op.hip.cu
@@ -322,7 +322,7 @@ __global__ void cal_stress_nl(
             const thrust::complex<FPTYPE> dbb1 = conj(dbecp[ib2 * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];
             const thrust::complex<FPTYPE> dbb2 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[ib2 * nkb + inkb2];
             const thrust::complex<FPTYPE> dbb3 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];
-            stress_var -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();
+            stress_var -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();
         }
         ++iat;
         sum+=nproj;
diff --git a/source/source_pw/module_pwdft/kernels/stress_op.cpp b/source/source_pw/module_pwdft/kernels/stress_op.cpp
@@ -232,7 +232,7 @@ struct cal_stress_nl_op<FPTYPE, base_device::DEVICE_CPU>
                                 const std::complex<FPTYPE> dbb1 = conj(dbecp[ib2 * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];
                                 const std::complex<FPTYPE> dbb2 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[ib2 * nkb + inkb2];
                                 const std::complex<FPTYPE> dbb3 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];
-                                local_stress -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();
+                                local_stress -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();
                             }
                         } // end ip
                     }     // ia

Original file line number	Diff line number	Diff line change
`@@ -270,7 +270,7 @@ __global__ void cal_force_nl(`
`270`	`270`	`const thrust::complex<FPTYPE> dbb1 = conj(dbecp[index0]) * becp[index1 + nkb];`
`271`	`271`	`const thrust::complex<FPTYPE> dbb2 = conj(dbecp[index0 + nkb]) * becp[index1];`
`272`	`272`	`const thrust::complex<FPTYPE> dbb3 = conj(dbecp[index0 + nkb]) * becp[index1 + nkb];`
`273`		`- const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();`
	`273`	`+ const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();`
`274`	`274`	`atomicAdd(force + iat * forcenl_nc + ipol, tmp);`
`275`	`275`	`}`
`276`	`276`	`}`
Original file line number	Diff line number	Diff line change
`@@ -361,7 +361,7 @@ __global__ void cal_stress_nl(`
`361`	`361`	`const thrust::complex<FPTYPE> dbb1 = conj(dbecp[ib2 * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];`
`362`	`362`	`const thrust::complex<FPTYPE> dbb2 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[ib2 * nkb + inkb2];`
`363`	`363`	`const thrust::complex<FPTYPE> dbb3 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];`
`364`		`- stress_var -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();`
	`364`	`+ stress_var -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();`
`365`	`365`	`}`
`366`	`366`	`++iat;`
`367`	`367`	`sum+=nproj;`
Original file line number	Diff line number	Diff line change
`@@ -254,7 +254,7 @@ struct cal_force_nl_op<FPTYPE, base_device::DEVICE_CPU>`
`254`	`254`	`const std::complex<FPTYPE> dbb2 = conj(dbecp[index0 + nkb]) * becp[index1];`
`255`	`255`	`const std::complex<FPTYPE> dbb3 = conj(dbecp[index0 + nkb]) * becp[index1 + nkb];`
`256`	`256`
`257`		`- local_force[ipol] -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();`
	`257`	`+ local_force[ipol] -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();`
`258`	`258`	`}`
`259`	`259`	`}`
`260`	`260`	`}`
Original file line number	Diff line number	Diff line change
`@@ -256,7 +256,7 @@ __global__ void cal_force_nl(`
`256`	`256`	`const thrust::complex<FPTYPE> dbb1 = conj(dbecp[index0]) * becp[index1 + nkb];`
`257`	`257`	`const thrust::complex<FPTYPE> dbb2 = conj(dbecp[index0 + nkb]) * becp[index1];`
`258`	`258`	`const thrust::complex<FPTYPE> dbb3 = conj(dbecp[index0 + nkb]) * becp[index1 + nkb];`
`259`		`- const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();`
	`259`	`+ const FPTYPE tmp = - fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();`
`260`	`260`	`atomicAdd(force + iat * forcenl_nc + ipol, tmp);`
`261`	`261`	`}`
`262`	`262`	`}`
Original file line number	Diff line number	Diff line change
`@@ -322,7 +322,7 @@ __global__ void cal_stress_nl(`
`322`	`322`	`const thrust::complex<FPTYPE> dbb1 = conj(dbecp[ib2 * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];`
`323`	`323`	`const thrust::complex<FPTYPE> dbb2 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[ib2 * nkb + inkb2];`
`324`	`324`	`const thrust::complex<FPTYPE> dbb3 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];`
`325`		`- stress_var -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();`
	`325`	`+ stress_var -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();`
`326`	`326`	`}`
`327`	`327`	`++iat;`
`328`	`328`	`sum+=nproj;`
Original file line number	Diff line number	Diff line change
`@@ -232,7 +232,7 @@ struct cal_stress_nl_op<FPTYPE, base_device::DEVICE_CPU>`
`232`	`232`	`const std::complex<FPTYPE> dbb1 = conj(dbecp[ib2 * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];`
`233`	`233`	`const std::complex<FPTYPE> dbb2 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[ib2 * nkb + inkb2];`
`234`	`234`	`const std::complex<FPTYPE> dbb3 = conj(dbecp[(ib2+1) * nkb + inkb1]) * becp[(ib2+1) * nkb + inkb2];`
`235`		`- local_stress -= fac * (ps0 * dbb0 + ps1 * dbb2 + ps2 * dbb1 + ps3 * dbb3).real();`
	`235`	`+ local_stress -= fac * (ps0 * dbb0 + ps1 * dbb1 + ps2 * dbb2 + ps3 * dbb3).real();`
`236`	`236`	`}`
`237`	`237`	`} // end ip`
`238`	`238`	`} // ia`