Skip to content

Commit d5c2916

Browse files
committed
fix: correct Pauli-to-spinor Hamiltonian conversion for nspin=4
Fix two bugs in LCAO non-collinear Hamiltonian construction: 1. Wrong sign in off-diagonal elements: H_{up,down} = B_x + i*B_y (wrong) should be B_x - i*B_y (correct), and vice versa for H_{down,up}. Fixed by correcting clx_j coefficients in merge_hr_part_to_hR(). 2. Missing complex conjugate in lower triangle fill: H(-R) used transpose instead of conjugate transpose, breaking Hermiticity for complex matrices. Fixed by using std::conj() when filling lower triangle. These errors caused the non-collinear Hamiltonian to be the complex conjugate of the correct result, leading to incorrect spin textures in nspin=4 calculations. The PW code path was not affected. Add test case and verification script to validate: - H(R=0) Hermiticity: max|H - H^dagger| < 1e-10 - Off-diagonal phase: Im(H_{up,down}) < 0 for m||+y direction See tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md for details.
1 parent 022948c commit d5c2916

6 files changed

Lines changed: 542 additions & 5 deletions

File tree

source/source_lcao/module_gint/gint_common.cpp

Lines changed: 13 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -168,8 +168,15 @@ void merge_hr_part_to_hR(const std::vector<hamilt::HContainer<double>>& hr_gint_
168168
std::vector<int> row_set = {0, 0, 1, 1};
169169
std::vector<int> col_set = {0, 1, 0, 1};
170170
//construct complex matrix
171+
// Pauli-to-spinor conversion: H = V_0*I + B_x*sigma_x + B_y*sigma_y + B_z*sigma_z
172+
// sigma_y = [[0,-i],[i,0]], so H_{up,down} = B_x - i*B_y, H_{down,up} = B_x + i*B_y
173+
// coefficient = clx_i + i*clx_j for each Pauli channel:
174+
// is=0 (up,up): V_0 + B_z => coeff on B_z = +1 => clx_i=1, clx_j=0
175+
// is=1 (up,down): B_x - i*B_y => coeff on B_y = -i => clx_i=0, clx_j=-1
176+
// is=2 (down,up): B_x + i*B_y => coeff on B_y = +i => clx_i=0, clx_j=+1
177+
// is=3 (down,down): -(V_0 - B_z) => coeff on V_0 = -1 => clx_i=-1, clx_j=0
171178
std::vector<int> clx_i = {1, 0, 0, -1};
172-
std::vector<int> clx_j = {0, 1, -1, 0};
179+
std::vector<int> clx_j = {0, -1, 1, 0};
173180
for (int is = 0; is < 4; is++){
174181
if(!PARAM.globalv.domag && (is==1 || is==2)) continue;
175182
hR_tmp->set_zero();
@@ -203,17 +210,18 @@ void merge_hr_part_to_hR(const std::vector<hamilt::HContainer<double>>& hr_gint_
203210
+ std::complex<double>(clx_i[is], clx_j[is]) * mat_nspin2->get_value(irow, icol);
204211
}
205212
}
206-
//fill the lower triangle matrix
207-
//When is=0 or 3, the real part does not need conjugation;
208-
//when is=1 or 2, the small matrix is not Hermitian, so conjugation is not needed
213+
//fill the lower triangle matrix at -R by conjugate transpose of upper at R
214+
// This ensures H(-R) = H(R)^dagger, required for Hermiticity of H(k).
215+
// For real matrices (is=0,3), conj has no effect.
216+
// For complex matrices (is=1,2), conj is essential.
209217
if (iat1 < iat2)
210218
{
211219
auto lower_mat = lower_ap->find_matrix(-R_index);
212220
for (int irow = 0; irow < upper_mat->get_row_size(); ++irow)
213221
{
214222
for (int icol = 0; icol < upper_mat->get_col_size(); ++icol)
215223
{
216-
lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol);
224+
lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol));
217225
}
218226
}
219227
}
Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
INPUT_PARAMETERS
2+
suffix verify_soc
3+
nbands 40
4+
5+
calculation scf
6+
ecutwfc 10
7+
scf_thr 1.0e-4
8+
scf_nmax 200
9+
out_chg 0
10+
11+
smearing_method gaussian
12+
smearing_sigma 0.01
13+
14+
mixing_type pulay
15+
mixing_beta 0.2
16+
mixing_restart 1e-3
17+
mixing_dmr 1
18+
mixing_gg0 1.1
19+
20+
ks_solver scalapack_gvx
21+
basis_type lcao
22+
gamma_only 0
23+
noncolin 1
24+
lspinorb 1
25+
26+
#Parameter DFT+U
27+
dft_plus_u 1
28+
orbital_corr 2
29+
hubbard_u 5.0
30+
onsite_radius 5.0
31+
pseudo_dir ../../PP_ORB
32+
orbital_dir ../../PP_ORB
33+
34+
# Output H(R) for verification
35+
out_mat_hs2 1
Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
K_POINTS
2+
0
3+
Gamma
4+
2 1 1 0 0 0
Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
ATOMIC_SPECIES
2+
Fe 1.000 Fe.upf
3+
4+
NUMERICAL_ORBITAL
5+
Fe_gga_6au_100Ry_4s2p2d1f.orb
6+
7+
LATTICE_CONSTANT
8+
8.190
9+
10+
LATTICE_VECTORS
11+
1.00 0.50 0.50
12+
0.50 1.00 0.50
13+
0.50 0.50 1.00
14+
ATOMIC_POSITIONS
15+
Direct
16+
17+
Fe
18+
0.0
19+
2
20+
0.00 0.00 0.00 mag 1.0 1.0 1.0
21+
0.51 0.51 0.51 mag 1.0 1.0 1.0
Lines changed: 223 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,223 @@
1+
# LCAO 非共线哈密顿量复共轭错误修复
2+
3+
## 0. 问题描述
4+
5+
在 LCAO 基组下,nspin=4(非共线自旋)计算中,从 Pauli 基的实数有效势
6+
`(V_0, B_x, B_y, B_z)` 构建旋量基复数哈密顿量的过程中存在两处错误:
7+
8+
1. **非对角元相位符号错误**`H_{↑,↓}``H_{↓,↑}` 的虚部符号相反
9+
2. **下三角填充缺少复共轭**`H(-R)` 使用转置而非共轭转置,破坏 Hermiticity
10+
11+
**影响**:非共线哈密顿量与正确结果互为复共轭(仅对非对角元),导致自洽后的
12+
波函数自旋纹理与磁化方向不一致。PW 代码路径不受影响。
13+
14+
**与 DM Fourier 变换符号修复的区别**:此 bug 位于哈密顿量构建环节
15+
(grid integration → H(R)),与 DM 的 k→R Fourier 变换无关。DM 的 Fourier
16+
符号修复(`density_matrix.cpp``dmr_complex.cpp` 中的 `sinp → -sinp`
17+
是另一个独立的问题。
18+
19+
---
20+
21+
## 1. 正确公式
22+
23+
### 1.1 Pauli 基 → 旋量基转换
24+
25+
有效势在 Pauli 基中为 4 个实数矩阵:
26+
```
27+
v(0) = V_0 (标量势,正比于电荷密度 n)
28+
v(1) = B_x (磁场 x 分量)
29+
v(2) = B_y (磁场 y 分量)
30+
v(3) = B_z (磁场 z 分量)
31+
```
32+
33+
旋量基哈密顿量为 `H = V_0 * I + B_x * σ_x + B_y * σ_y + B_z * σ_z`,其中:
34+
35+
```
36+
σ_x = [[0, 1], [1, 0]]
37+
σ_y = [[0, -i], [i, 0]] ← 关键:σ_y 的非对角元含 ±i
38+
σ_z = [[1, 0], [0, -1]]
39+
```
40+
41+
展开得 4 个旋量矩阵元:
42+
43+
```
44+
H_{↑,↑} = V_0 + B_z (纯实数)
45+
H_{↑,↓} = B_x - i*B_y (虚部为负)
46+
H_{↓,↑} = B_x + i*B_y (虚部为正)
47+
H_{↓,↓} = V_0 - B_z (纯实数)
48+
```
49+
50+
### 1.2 实空间 Hermiticity
51+
52+
哈密顿量矩阵元 `H_{μν}(R) = <φ_{μ0}|H|φ_{νR}>` 需满足:
53+
54+
```
55+
H_{μν}(R) = H*_{νμ}(-R)
56+
```
57+
58+
`H(-R) = H†(R)`。在代码中,对 atom pair `(iat1, iat2)`,若已知
59+
`(iat1, iat2, R)` 处的矩阵 `H_upper`,则 `(iat2, iat1, -R)` 处的矩阵为:
60+
61+
```
62+
H_lower(-R) = H_upper†(R) = conj(H_upper^T(R))
63+
```
64+
65+
---
66+
67+
## 2. 代码中的错误
68+
69+
### 2.1 错误位置
70+
71+
文件:`source/source_lcao/module_gint/gint_common.cpp`
72+
函数:`merge_hr_part_to_hR()`
73+
74+
### 2.2 错误 1:非对角元相位
75+
76+
代码用系数 `(clx_i + i*clx_j)` 组合 Pauli 分量:
77+
78+
```cpp
79+
// 错误代码:
80+
std::vector<int> clx_i = {1, 0, 0, -1};
81+
std::vector<int> clx_j = {0, 1, -1, 0};
82+
// ↑ ↑
83+
// is=1: +i is=2: -i
84+
// 应为: -i 应为: +i
85+
```
86+
87+
这导致:
88+
```
89+
H_{↑,↓} = B_x + i*B_y (错误:虚部为正)
90+
H_{↓,↑} = B_x - i*B_y (错误:虚部为负)
91+
```
92+
93+
与正确公式恰好互为复共轭。
94+
95+
### 2.3 错误 2:下三角填充
96+
97+
```cpp
98+
// 错误代码:只转置,不取共轭
99+
lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol);
100+
```
101+
102+
对实数矩阵 (is=0,3),转置 = 共轭转置,无影响。
103+
对复数矩阵 (is=1,2),转置 ≠ 共轭转置,导致 `H(-R) ≠ H†(R)`
104+
105+
---
106+
107+
## 3. 修复
108+
109+
```cpp
110+
// 修复后:
111+
std::vector<int> clx_j = {0, -1, 1, 0};
112+
// ↑ ↑
113+
// is=1: -i is=2: +i ← 正确
114+
115+
// 下三角填充:使用共轭转置
116+
lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol));
117+
```
118+
119+
---
120+
121+
## 4. 测试方案
122+
123+
### 4.1 测试体系
124+
125+
使用 bcc Fe 二聚体 + DFT+U,磁化方向沿 (1,1,1),开启 SOC:
126+
127+
```
128+
STRU:
129+
Fe 0.0 0.00 0.00 0.00 mag 1.0 1.0 1.0
130+
Fe 0.51 0.51 0.51 0.51 mag 1.0 1.0 1.0
131+
```
132+
133+
该体系具有稳定的磁矩 (~3.47 μB/atom),自洽后可产生显著的
134+
非对角哈密顿量矩阵元。
135+
136+
### 4.2 验证条件
137+
138+
**条件 1:H(R=0) 的 Hermiticity**
139+
140+
```
141+
max|H(R=0) - H†(R=0)| < 1e-10
142+
```
143+
144+
修复前:~1e-2(因下三角缺少共轭)
145+
修复后:~1e-19(机器精度)
146+
147+
**条件 2:非对角元相位符号**
148+
149+
对磁化方向 m = (mx, my, mz),XC 势的非对角元为:
150+
```
151+
H_{↑,↓} ∝ mx - i*my
152+
```
153+
154+
当 my > 0 时,Im(H_{↑,↓}) < 0。
155+
156+
测试中 m 沿 (1,1,1),故 Im(H_{↑,↓}) 应为负值。
157+
158+
修复前:Im(H_{↑,↓}) > 0(错误)
159+
修复后:Im(H_{↑,↓}) < 0(正确)
160+
161+
**条件 3:自旋纹理一致性(定性)**
162+
163+
对角化 H(k=0),计算占据态的自旋期望值 `<σ>`。
164+
修复后,`<σ>` 应与输入磁化方向 (1,1,1) 一致。
165+
修复前,`<σ>` 指向共轭后的方向。
166+
167+
### 4.3 测试文件
168+
169+
| 文件 | 说明 |
170+
|------|------|
171+
| `tests/03_NAO_multik/verify_hamiltonian_convention/` | 测试目录 |
172+
| `check_hamiltonian_convention.py` | Python 验证脚本 |
173+
| `tests/03_NAO_multik/verify_dm_symmetry_soc/` | 复用已有的 SOC 测试算例 |
174+
175+
### 4.4 运行方式
176+
177+
```bash
178+
# 运行 SCF 计算(需要 out_mat_hs2=1 输出 H(R))
179+
cd tests/03_NAO_multik/verify_dm_symmetry_soc
180+
mpirun -np 4 abacus > log.txt 2>&1
181+
182+
# 运行验证脚本
183+
python3 ../verify_hamiltonian_convention/check_hamiltonian_convention.py OUT.verify_soc/
184+
```
185+
186+
预期输出:
187+
```
188+
TEST 1: H(R=0) Hermiticity: H = H^dagger
189+
max|H - H^dagger| = 1.73e-19
190+
[PASS]
191+
192+
TEST 2: Off-diagonal phase for m along y
193+
Mean Im(H_{up,down}) = -1.59e-03
194+
[PASS] Im(H_{up,down}) < 0, consistent with correct convention
195+
```
196+
197+
---
198+
199+
## 5. PW vs LCAO 对比
200+
201+
| 代码路径 | 文件 | 状态 |
202+
|---------|------|------|
203+
| PW | `source/source_pw/module_pwdft/kernels/veff_op.cpp` | **正确** |
204+
| LCAO | `source/source_lcao/module_gint/gint_common.cpp` | **已修复** |
205+
206+
PW 代码直接实现了正确的公式:
207+
```cpp
208+
// veff_op.cpp: 正确
209+
sup = out*(in[0]+in[3]) + out1*(in[1] - i*in[2]); // V_0+V_z, V_x-i*V_y
210+
sdown = out1*(in[0]-in[3]) + out*(in[1] + i*in[2]); // V_0-V_z, V_x+i*V_y
211+
```
212+
213+
LCAO 代码通过 `merge_hr_part_to_hR()` 间接转换,修复前符号错误。
214+
215+
---
216+
217+
## 6. 影响范围
218+
219+
- 所有 nspin=4 的 LCAO 计算(包括 SOC 和非共线无 SOC)
220+
- CPU 和 GPU 路径均受影响(GPU 代码调用同一 `merge_hr_part_to_hR` 函数)
221+
- regular 和 meta-GGA 均受影响
222+
223+
nspin=1,2(共线自旋)不受影响,因为此时非对角元为零。

0 commit comments

Comments
 (0)