From eb940669894c6fe8acc355b730c6ef45b71c42d2 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 23 Jun 2026 15:50:48 +0800 Subject: [PATCH 1/3] 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. --- .../source_lcao/module_gint/gint_common.cpp | 18 +- .../verify_hamiltonian_convention/INPUT | 35 +++ .../verify_hamiltonian_convention/KPT | 4 + .../verify_hamiltonian_convention/STRU | 21 ++ .../TEST_DESIGN.md | 223 ++++++++++++++++ .../check_hamiltonian_convention.py | 246 ++++++++++++++++++ 6 files changed, 542 insertions(+), 5 deletions(-) create mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/INPUT create mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/KPT create mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/STRU create mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md create mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py diff --git a/source/source_lcao/module_gint/gint_common.cpp b/source/source_lcao/module_gint/gint_common.cpp index 987fcbe5a0e..1de8969fd0d 100644 --- a/source/source_lcao/module_gint/gint_common.cpp +++ b/source/source_lcao/module_gint/gint_common.cpp @@ -168,8 +168,15 @@ void merge_hr_part_to_hR(const std::vector>& hr_gint_ std::vector row_set = {0, 0, 1, 1}; std::vector 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 clx_i = {1, 0, 0, -1}; - std::vector clx_j = {0, 1, -1, 0}; + std::vector 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(); @@ -203,9 +210,10 @@ void merge_hr_part_to_hR(const std::vector>& hr_gint_ + std::complex(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); @@ -213,7 +221,7 @@ void merge_hr_part_to_hR(const std::vector>& hr_gint_ { 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)); } } } diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/INPUT b/tests/03_NAO_multik/verify_hamiltonian_convention/INPUT new file mode 100644 index 00000000000..2a2a4f6b6c7 --- /dev/null +++ b/tests/03_NAO_multik/verify_hamiltonian_convention/INPUT @@ -0,0 +1,35 @@ +INPUT_PARAMETERS +suffix verify_soc +nbands 40 + +calculation scf +ecutwfc 10 +scf_thr 1.0e-4 +scf_nmax 200 +out_chg 0 + +smearing_method gaussian +smearing_sigma 0.01 + +mixing_type pulay +mixing_beta 0.2 +mixing_restart 1e-3 +mixing_dmr 1 +mixing_gg0 1.1 + +ks_solver scalapack_gvx +basis_type lcao +gamma_only 0 +noncolin 1 +lspinorb 1 + +#Parameter DFT+U +dft_plus_u 1 +orbital_corr 2 +hubbard_u 5.0 +onsite_radius 5.0 +pseudo_dir ../../PP_ORB +orbital_dir ../../PP_ORB + +# Output H(R) for verification +out_mat_hs2 1 diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/KPT b/tests/03_NAO_multik/verify_hamiltonian_convention/KPT new file mode 100644 index 00000000000..e769af76382 --- /dev/null +++ b/tests/03_NAO_multik/verify_hamiltonian_convention/KPT @@ -0,0 +1,4 @@ +K_POINTS +0 +Gamma +2 1 1 0 0 0 diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/STRU b/tests/03_NAO_multik/verify_hamiltonian_convention/STRU new file mode 100644 index 00000000000..10cb3175d61 --- /dev/null +++ b/tests/03_NAO_multik/verify_hamiltonian_convention/STRU @@ -0,0 +1,21 @@ +ATOMIC_SPECIES +Fe 1.000 Fe.upf + +NUMERICAL_ORBITAL +Fe_gga_6au_100Ry_4s2p2d1f.orb + +LATTICE_CONSTANT +8.190 + +LATTICE_VECTORS + 1.00 0.50 0.50 + 0.50 1.00 0.50 + 0.50 0.50 1.00 +ATOMIC_POSITIONS +Direct + +Fe +0.0 +2 +0.00 0.00 0.00 mag 1.0 1.0 1.0 +0.51 0.51 0.51 mag 1.0 1.0 1.0 diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md b/tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md new file mode 100644 index 00000000000..1fe49c367ba --- /dev/null +++ b/tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md @@ -0,0 +1,223 @@ +# LCAO 非共线哈密顿量复共轭错误修复 + +## 0. 问题描述 + +在 LCAO 基组下，nspin=4（非共线自旋）计算中，从 Pauli 基的实数有效势 +`(V_0, B_x, B_y, B_z)` 构建旋量基复数哈密顿量的过程中存在两处错误： + +1. **非对角元相位符号错误**：`H_{↑,↓}` 和 `H_{↓,↑}` 的虚部符号相反 +2. **下三角填充缺少复共轭**：`H(-R)` 使用转置而非共轭转置，破坏 Hermiticity + +**影响**：非共线哈密顿量与正确结果互为复共轭（仅对非对角元），导致自洽后的 +波函数自旋纹理与磁化方向不一致。PW 代码路径不受影响。 + +**与 DM Fourier 变换符号修复的区别**：此 bug 位于哈密顿量构建环节 +（grid integration → H(R)），与 DM 的 k→R Fourier 变换无关。DM 的 Fourier +符号修复（`density_matrix.cpp` 和 `dmr_complex.cpp` 中的 `sinp → -sinp`） +是另一个独立的问题。 + +--- + +## 1. 正确公式 + +### 1.1 Pauli 基 → 旋量基转换 + +有效势在 Pauli 基中为 4 个实数矩阵： +``` +v(0) = V_0 (标量势，正比于电荷密度 n) +v(1) = B_x (磁场 x 分量) +v(2) = B_y (磁场 y 分量) +v(3) = B_z (磁场 z 分量) +``` + +旋量基哈密顿量为 `H = V_0 * I + B_x * σ_x + B_y * σ_y + B_z * σ_z`，其中： + +``` +σ_x = [[0, 1], [1, 0]] +σ_y = [[0, -i], [i, 0]] ← 关键：σ_y 的非对角元含 ±i +σ_z = [[1, 0], [0, -1]] +``` + +展开得 4 个旋量矩阵元： + +``` +H_{↑,↑} = V_0 + B_z （纯实数） +H_{↑,↓} = B_x - i*B_y （虚部为负） +H_{↓,↑} = B_x + i*B_y （虚部为正） +H_{↓,↓} = V_0 - B_z （纯实数） +``` + +### 1.2 实空间 Hermiticity + +哈密顿量矩阵元 `H_{μν}(R) = <φ_{μ0}|H|φ_{νR}>` 需满足： + +``` +H_{μν}(R) = H*_{νμ}(-R) +``` + +即 `H(-R) = H†(R)`。在代码中，对 atom pair `(iat1, iat2)`，若已知 +`(iat1, iat2, R)` 处的矩阵 `H_upper`，则 `(iat2, iat1, -R)` 处的矩阵为： + +``` +H_lower(-R) = H_upper†(R) = conj(H_upper^T(R)) +``` + +--- + +## 2. 代码中的错误 + +### 2.1 错误位置 + +文件：`source/source_lcao/module_gint/gint_common.cpp` +函数：`merge_hr_part_to_hR()` + +### 2.2 错误 1：非对角元相位 + +代码用系数 `(clx_i + i*clx_j)` 组合 Pauli 分量： + +```cpp +// 错误代码： +std::vector clx_i = {1, 0, 0, -1}; +std::vector clx_j = {0, 1, -1, 0}; +// ↑ ↑ +// is=1: +i is=2: -i +// 应为: -i 应为: +i +``` + +这导致： +``` +H_{↑,↓} = B_x + i*B_y （错误：虚部为正） +H_{↓,↑} = B_x - i*B_y （错误：虚部为负） +``` + +与正确公式恰好互为复共轭。 + +### 2.3 错误 2：下三角填充 + +```cpp +// 错误代码：只转置，不取共轭 +lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol); +``` + +对实数矩阵 (is=0,3)，转置 = 共轭转置，无影响。 +对复数矩阵 (is=1,2)，转置 ≠ 共轭转置，导致 `H(-R) ≠ H†(R)`。 + +--- + +## 3. 修复 + +```cpp +// 修复后： +std::vector clx_j = {0, -1, 1, 0}; +// ↑ ↑ +// is=1: -i is=2: +i ← 正确 + +// 下三角填充：使用共轭转置 +lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol)); +``` + +--- + +## 4. 测试方案 + +### 4.1 测试体系 + +使用 bcc Fe 二聚体 + DFT+U，磁化方向沿 (1,1,1)，开启 SOC： + +``` +STRU: + Fe 0.0 0.00 0.00 0.00 mag 1.0 1.0 1.0 + Fe 0.51 0.51 0.51 0.51 mag 1.0 1.0 1.0 +``` + +该体系具有稳定的磁矩 (~3.47 μB/atom)，自洽后可产生显著的 +非对角哈密顿量矩阵元。 + +### 4.2 验证条件 + +**条件 1：H(R=0) 的 Hermiticity** + +``` +max|H(R=0) - H†(R=0)| < 1e-10 +``` + +修复前：~1e-2（因下三角缺少共轭） +修复后：~1e-19（机器精度） + +**条件 2：非对角元相位符号** + +对磁化方向 m = (mx, my, mz)，XC 势的非对角元为： +``` +H_{↑,↓} ∝ mx - i*my +``` + +当 my > 0 时，Im(H_{↑,↓}) < 0。 + +测试中 m 沿 (1,1,1)，故 Im(H_{↑,↓}) 应为负值。 + +修复前：Im(H_{↑,↓}) > 0（错误） +修复后：Im(H_{↑,↓}) < 0（正确） + +**条件 3：自旋纹理一致性（定性）** + +对角化 H(k=0)，计算占据态的自旋期望值 `<σ>`。 +修复后，`<σ>` 应与输入磁化方向 (1,1,1) 一致。 +修复前，`<σ>` 指向共轭后的方向。 + +### 4.3 测试文件 + +| 文件 | 说明 | +|------|------| +| `tests/03_NAO_multik/verify_hamiltonian_convention/` | 测试目录 | +| `check_hamiltonian_convention.py` | Python 验证脚本 | +| `tests/03_NAO_multik/verify_dm_symmetry_soc/` | 复用已有的 SOC 测试算例 | + +### 4.4 运行方式 + +```bash +# 运行 SCF 计算（需要 out_mat_hs2=1 输出 H(R)） +cd tests/03_NAO_multik/verify_dm_symmetry_soc +mpirun -np 4 abacus > log.txt 2>&1 + +# 运行验证脚本 +python3 ../verify_hamiltonian_convention/check_hamiltonian_convention.py OUT.verify_soc/ +``` + +预期输出： +``` +TEST 1: H(R=0) Hermiticity: H = H^dagger + max|H - H^dagger| = 1.73e-19 + [PASS] + +TEST 2: Off-diagonal phase for m along y + Mean Im(H_{up,down}) = -1.59e-03 + [PASS] Im(H_{up,down}) < 0, consistent with correct convention +``` + +--- + +## 5. PW vs LCAO 对比 + +| 代码路径 | 文件 | 状态 | +|---------|------|------| +| PW | `source/source_pw/module_pwdft/kernels/veff_op.cpp` | **正确** | +| LCAO | `source/source_lcao/module_gint/gint_common.cpp` | **已修复** | + +PW 代码直接实现了正确的公式： +```cpp +// veff_op.cpp: 正确 +sup = out*(in[0]+in[3]) + out1*(in[1] - i*in[2]); // V_0+V_z, V_x-i*V_y +sdown = out1*(in[0]-in[3]) + out*(in[1] + i*in[2]); // V_0-V_z, V_x+i*V_y +``` + +LCAO 代码通过 `merge_hr_part_to_hR()` 间接转换，修复前符号错误。 + +--- + +## 6. 影响范围 + +- 所有 nspin=4 的 LCAO 计算（包括 SOC 和非共线无 SOC） +- CPU 和 GPU 路径均受影响（GPU 代码调用同一 `merge_hr_part_to_hR` 函数） +- regular 和 meta-GGA 均受影响 + +nspin=1,2（共线自旋）不受影响，因为此时非对角元为零。 diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py b/tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py new file mode 100644 index 00000000000..04e7ed7c06e --- /dev/null +++ b/tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py @@ -0,0 +1,246 @@ +#!/usr/bin/env python3 +""" +Simple Hamiltonian convention check. +Reads the first R=0 block of H(R) and checks: +1. Hermiticity of R=0 block +2. Off-diagonal phase for m along y +""" +import sys +import re +import numpy as np + +def read_complex_value(s): + """Parse a complex value like '(-1.13e-02,0.00e+00)' or a real number.""" + s = s.strip().strip('()') + if ',' in s: + parts = s.split(',') + return complex(float(parts[0]), float(parts[1])) + else: + return complex(float(s), 0.0) + +def parse_hr_r0(filepath): + """Parse the R=(0,0,0) block from hrs1_nao.csr.""" + with open(filepath) as f: + lines = f.readlines() + + # Parse header + idx = 0 + nbasis = None + nR = None + for i, line in enumerate(lines): + stripped = line.strip() + if stripped.startswith('#') or stripped == '' or stripped.startswith('---'): + continue + if 'number of spin' in stripped or 'spin index' in stripped: + continue + if 'number of localized basis' in stripped: + nbasis = int(stripped.split()[0]) + continue + if 'number of Bravais' in stripped: + nR = int(stripped.split()[0]) + break + + # Find R=(0,0,0) block + found_r0 = False + for i, line in enumerate(lines): + stripped = line.strip() + # Look for "0 0 0 nnz" pattern + parts = stripped.split() + if len(parts) == 4: + try: + rx, ry, rz, nnz = int(parts[0]), int(parts[1]), int(parts[2]), int(parts[3]) + if rx == 0 and ry == 0 and rz == 0 and nnz > 0: + found_r0 = True + break + except ValueError: + continue + + if not found_r0: + print("R=(0,0,0) block not found!") + return None, None, None + + # Now parse values, column indices, row pointers after this line + idx = i + 1 + # Skip "# CSR values" + while idx < len(lines) and '# CSR values' not in lines[idx]: + idx += 1 + idx += 1 + + # Read nnz complex values + values = [] + while len(values) < nnz and idx < len(lines): + line = lines[idx].strip() + idx += 1 + if line.startswith('#'): + break + if not line: + continue + # Parse all complex values on this line + tokens = re.findall(r'\([^)]+\)', line) + if not tokens: + # Try space-separated values + tokens = line.split() + for t in tokens: + try: + values.append(read_complex_value(t)) + except ValueError: + break + if len(values) >= nnz: + break + + # Skip "# CSR column indices" + while idx < len(lines) and '# CSR column indices' not in lines[idx]: + idx += 1 + idx += 1 + + # Read nnz column indices + col_indices = [] + while len(col_indices) < nnz and idx < len(lines): + line = lines[idx].strip() + idx += 1 + if line.startswith('#'): + break + if not line: + continue + for t in line.split(): + try: + col_indices.append(int(t)) + except ValueError: + break + if len(col_indices) >= nnz: + break + + # Skip "# CSR row pointers" + while idx < len(lines) and '# CSR row pointers' not in lines[idx]: + idx += 1 + idx += 1 + + # Read nbasis+1 row pointers + row_ptrs = [] + while len(row_ptrs) < nbasis + 1 and idx < len(lines): + line = lines[idx].strip() + idx += 1 + if line.startswith('#'): + break + if not line: + continue + for t in line.split(): + try: + row_ptrs.append(int(t)) + except ValueError: + break + if len(row_ptrs) >= nbasis + 1: + break + + # Build dense matrix + H = np.zeros((nbasis, nbasis), dtype=complex) + for irow in range(nbasis): + if irow >= len(row_ptrs): + break + start = row_ptrs[irow] + end = row_ptrs[irow + 1] if irow + 1 < len(row_ptrs) else nnz + for k in range(start, min(end, nnz, len(values), len(col_indices))): + icol = col_indices[k] + if 0 <= icol < nbasis: + H[irow, icol] = values[k] + + return nbasis, H, nnz + +def main(): + if len(sys.argv) < 2: + print("Usage: python3 check_hamiltonian_convention.py ") + sys.exit(1) + + out_dir = sys.argv[1] + hr_file = None + dmr_file = None + + import os + for f in os.listdir(out_dir): + if f.startswith('hrs') and f.endswith('.csr'): + hr_file = os.path.join(out_dir, f) + if f.startswith('dmrs') and f.endswith('.csr'): + dmr_file = os.path.join(out_dir, f) + + if not hr_file: + print("ERROR: hrs*_nao.csr not found") + sys.exit(1) + + print(f"Parsing H(R=0) from: {hr_file}") + nbasis, H0, nnz = parse_hr_r0(hr_file) + if H0 is None: + print("Failed to parse H(R=0)") + sys.exit(1) + + print(f" nbasis = {nbasis}, nnz = {nnz}") + N = nbasis // 2 # number of spatial orbitals per spin + + # Test 1: Hermiticity of R=0 block + print("\n" + "=" * 70) + print("TEST 1: H(R=0) Hermiticity: H = H^dagger") + print("=" * 70) + herm_err = np.max(np.abs(H0 - H0.conj().T)) + print(f" max|H - H^dagger| = {herm_err:.2e}") + if herm_err < 1e-10: + print(" [PASS]") + else: + print(" [FAIL]") + + # Test 2: Off-diagonal spinor block phase + print("\n" + "=" * 70) + print("TEST 2: Off-diagonal phase for m along y") + print(" For m||+y: H_{up,down} should have Im < 0") + print(" (i.e., H(2i, 2j+1) should be -i*|B_y|)") + print("=" * 70) + + # Extract the up-down block: rows 0,2,4,..., cols 1,3,5,... + H_updown = np.zeros((N, N), dtype=complex) + for i in range(N): + for j in range(N): + H_updown[i, j] = H0[2*i, 2*j+1] + + # Check the sign of imaginary part + nonzero_mask = np.abs(H_updown) > 1e-10 + if np.any(nonzero_mask): + imag_parts = H_updown[nonzero_mask].imag + mean_imag = np.mean(imag_parts) + n_neg = np.sum(imag_parts < 0) + n_pos = np.sum(imag_parts > 0) + print(f" Non-zero H_{{up,down}} elements: {np.sum(nonzero_mask)}") + print(f" Mean Im(H_{{up,down}}) = {mean_imag:.6e}") + print(f" Negative: {n_neg}, Positive: {n_pos}") + if mean_imag < -1e-10: + print(" [PASS] Im(H_{{up,down}}) < 0, consistent with correct convention") + elif mean_imag > 1e-10: + print(" [FAIL] Im(H_{{up,down}}) > 0, indicates complex conjugation error") + else: + print(" [INFO] Im(H_{{up,down}}) ~ 0, magnetization may not be along y") + else: + print(" No significant off-diagonal elements found") + + # Test 3: Diagonal spinor block structure + print("\n" + "=" * 70) + print("TEST 3: Diagonal blocks H_{{up,up}} and H_{{down,down}}") + print(" For m||y with no SOC: H_{{up,up}} = H_{{down,down}} (real)") + print("=" * 70) + H_upup = np.zeros((N, N), dtype=complex) + H_downdown = np.zeros((N, N), dtype=complex) + for i in range(N): + for j in range(N): + H_upup[i, j] = H0[2*i, 2*j] + H_downdown[i, j] = H0[2*i+1, 2*j+1] + + diag_diff = np.max(np.abs(H_upup - H_downdown)) + upup_imag = np.max(np.abs(H_upup.imag)) + dndn_imag = np.max(np.abs(H_downdown.imag)) + print(f" max|H_{{up,up}} - H_{{down,down}}| = {diag_diff:.2e}") + print(f" max|Im(H_{{up,up}})| = {upup_imag:.2e}") + print(f" max|Im(H_{{down,down}})| = {dndn_imag:.2e}") + + if diag_diff < 0.1: + print(" [PASS] Diagonal blocks are nearly equal (expected for m perp to z)") + else: + print(" [INFO] Diagonal blocks differ (may indicate m has z-component or SOC)") + +if __name__ == '__main__': + main() From 4e5b62a0037d0ae054fdf8d50a8be8445e7d9fc9 Mon Sep 17 00:00:00 2001 From: dyzheng Date: Tue, 23 Jun 2026 21:17:32 +0800 Subject: [PATCH 2/3] fix: correct Pauli-to-spinor conversion in DFT+U and DeltaSpin for nspin=4 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- .../dftu_force_stress.hpp | 20 +- .../module_operator_lcao/dftu_lcao.cpp | 14 +- .../dspin_force_stress.hpp | 224 +++++++++++++----- 3 files changed, 193 insertions(+), 65 deletions(-) diff --git a/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp b/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp index 38c96025fa5..50b157c03e2 100644 --- a/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/dftu_force_stress.hpp @@ -148,9 +148,25 @@ void DFTU>::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 (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 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) } } diff --git a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp index 2b8e4f3506f..92fac0ffabd 100644 --- a/source/source_lcao/module_operator_lcao/dftu_lcao.cpp +++ b/source/source_lcao/module_operator_lcao/dftu_lcao.cpp @@ -612,7 +612,14 @@ void hamilt::DFTU, std::complex transfer from double to std::complex + // 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, std::complex type, but here we use double type for test - vu[index[1]] = 0.5 * (vu_tmp[index[1]] + std::complex(0.0, 1.0) * vu_tmp[index[2]]); - vu[index[2]] = 0.5 * (vu_tmp[index[1]] - std::complex(0.0, 1.0) * vu_tmp[index[2]]); + vu[index[1]] = 0.5 * (vu_tmp[index[1]] - std::complex(0.0, 1.0) * vu_tmp[index[2]]); + vu[index[2]] = 0.5 * (vu_tmp[index[1]] + std::complex(0.0, 1.0) * vu_tmp[index[2]]); } } } diff --git a/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp b/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp index 955c6b4d267..0ca04bb76fa 100644 --- a/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp +++ b/source/source_lcao/module_operator_lcao/dspin_force_stress.hpp @@ -274,43 +274,95 @@ void DeltaSpin>::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& 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& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const std::vector& 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; lget_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const std::vector& 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& 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 * - // force2 = - VU * * } - 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>::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& 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& nlm2 = nlm2_all.find(col_indexes[iw2l])->second; + const std::vector& nlm1 = nlm1_all.find(row_indexes[iw1l])->second; + for (int iw2l = 0; iw2l < col_indexes.size(); iw2l += npol) + { + const std::vector& 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; lget_pointer(); + for (int iw1l = 0; iw1l < row_indexes.size(); iw1l += npol) + { + const std::vector& 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& 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 Date: Tue, 23 Jun 2026 22:54:06 +0800 Subject: [PATCH 3/3] chore: remove .py and .md test files from PR --- .../TEST_DESIGN.md | 223 ---------------- .../check_hamiltonian_convention.py | 246 ------------------ 2 files changed, 469 deletions(-) delete mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md delete mode 100644 tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md b/tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md deleted file mode 100644 index 1fe49c367ba..00000000000 --- a/tests/03_NAO_multik/verify_hamiltonian_convention/TEST_DESIGN.md +++ /dev/null @@ -1,223 +0,0 @@ -# LCAO 非共线哈密顿量复共轭错误修复 - -## 0. 问题描述 - -在 LCAO 基组下，nspin=4（非共线自旋）计算中，从 Pauli 基的实数有效势 -`(V_0, B_x, B_y, B_z)` 构建旋量基复数哈密顿量的过程中存在两处错误： - -1. **非对角元相位符号错误**：`H_{↑,↓}` 和 `H_{↓,↑}` 的虚部符号相反 -2. **下三角填充缺少复共轭**：`H(-R)` 使用转置而非共轭转置，破坏 Hermiticity - -**影响**：非共线哈密顿量与正确结果互为复共轭（仅对非对角元），导致自洽后的 -波函数自旋纹理与磁化方向不一致。PW 代码路径不受影响。 - -**与 DM Fourier 变换符号修复的区别**：此 bug 位于哈密顿量构建环节 -（grid integration → H(R)），与 DM 的 k→R Fourier 变换无关。DM 的 Fourier -符号修复（`density_matrix.cpp` 和 `dmr_complex.cpp` 中的 `sinp → -sinp`） -是另一个独立的问题。 - ---- - -## 1. 正确公式 - -### 1.1 Pauli 基 → 旋量基转换 - -有效势在 Pauli 基中为 4 个实数矩阵： -``` -v(0) = V_0 (标量势，正比于电荷密度 n) -v(1) = B_x (磁场 x 分量) -v(2) = B_y (磁场 y 分量) -v(3) = B_z (磁场 z 分量) -``` - -旋量基哈密顿量为 `H = V_0 * I + B_x * σ_x + B_y * σ_y + B_z * σ_z`，其中： - -``` -σ_x = [[0, 1], [1, 0]] -σ_y = [[0, -i], [i, 0]] ← 关键：σ_y 的非对角元含 ±i -σ_z = [[1, 0], [0, -1]] -``` - -展开得 4 个旋量矩阵元： - -``` -H_{↑,↑} = V_0 + B_z （纯实数） -H_{↑,↓} = B_x - i*B_y （虚部为负） -H_{↓,↑} = B_x + i*B_y （虚部为正） -H_{↓,↓} = V_0 - B_z （纯实数） -``` - -### 1.2 实空间 Hermiticity - -哈密顿量矩阵元 `H_{μν}(R) = <φ_{μ0}|H|φ_{νR}>` 需满足： - -``` -H_{μν}(R) = H*_{νμ}(-R) -``` - -即 `H(-R) = H†(R)`。在代码中，对 atom pair `(iat1, iat2)`，若已知 -`(iat1, iat2, R)` 处的矩阵 `H_upper`，则 `(iat2, iat1, -R)` 处的矩阵为： - -``` -H_lower(-R) = H_upper†(R) = conj(H_upper^T(R)) -``` - ---- - -## 2. 代码中的错误 - -### 2.1 错误位置 - -文件：`source/source_lcao/module_gint/gint_common.cpp` -函数：`merge_hr_part_to_hR()` - -### 2.2 错误 1：非对角元相位 - -代码用系数 `(clx_i + i*clx_j)` 组合 Pauli 分量： - -```cpp -// 错误代码： -std::vector clx_i = {1, 0, 0, -1}; -std::vector clx_j = {0, 1, -1, 0}; -// ↑ ↑ -// is=1: +i is=2: -i -// 应为: -i 应为: +i -``` - -这导致： -``` -H_{↑,↓} = B_x + i*B_y （错误：虚部为正） -H_{↓,↑} = B_x - i*B_y （错误：虚部为负） -``` - -与正确公式恰好互为复共轭。 - -### 2.3 错误 2：下三角填充 - -```cpp -// 错误代码：只转置，不取共轭 -lower_mat->get_value(icol, irow) = upper_mat->get_value(irow, icol); -``` - -对实数矩阵 (is=0,3)，转置 = 共轭转置，无影响。 -对复数矩阵 (is=1,2)，转置 ≠ 共轭转置，导致 `H(-R) ≠ H†(R)`。 - ---- - -## 3. 修复 - -```cpp -// 修复后： -std::vector clx_j = {0, -1, 1, 0}; -// ↑ ↑ -// is=1: -i is=2: +i ← 正确 - -// 下三角填充：使用共轭转置 -lower_mat->get_value(icol, irow) = std::conj(upper_mat->get_value(irow, icol)); -``` - ---- - -## 4. 测试方案 - -### 4.1 测试体系 - -使用 bcc Fe 二聚体 + DFT+U，磁化方向沿 (1,1,1)，开启 SOC： - -``` -STRU: - Fe 0.0 0.00 0.00 0.00 mag 1.0 1.0 1.0 - Fe 0.51 0.51 0.51 0.51 mag 1.0 1.0 1.0 -``` - -该体系具有稳定的磁矩 (~3.47 μB/atom)，自洽后可产生显著的 -非对角哈密顿量矩阵元。 - -### 4.2 验证条件 - -**条件 1：H(R=0) 的 Hermiticity** - -``` -max|H(R=0) - H†(R=0)| < 1e-10 -``` - -修复前：~1e-2（因下三角缺少共轭） -修复后：~1e-19（机器精度） - -**条件 2：非对角元相位符号** - -对磁化方向 m = (mx, my, mz)，XC 势的非对角元为： -``` -H_{↑,↓} ∝ mx - i*my -``` - -当 my > 0 时，Im(H_{↑,↓}) < 0。 - -测试中 m 沿 (1,1,1)，故 Im(H_{↑,↓}) 应为负值。 - -修复前：Im(H_{↑,↓}) > 0（错误） -修复后：Im(H_{↑,↓}) < 0（正确） - -**条件 3：自旋纹理一致性（定性）** - -对角化 H(k=0)，计算占据态的自旋期望值 `<σ>`。 -修复后，`<σ>` 应与输入磁化方向 (1,1,1) 一致。 -修复前，`<σ>` 指向共轭后的方向。 - -### 4.3 测试文件 - -| 文件 | 说明 | -|------|------| -| `tests/03_NAO_multik/verify_hamiltonian_convention/` | 测试目录 | -| `check_hamiltonian_convention.py` | Python 验证脚本 | -| `tests/03_NAO_multik/verify_dm_symmetry_soc/` | 复用已有的 SOC 测试算例 | - -### 4.4 运行方式 - -```bash -# 运行 SCF 计算（需要 out_mat_hs2=1 输出 H(R)） -cd tests/03_NAO_multik/verify_dm_symmetry_soc -mpirun -np 4 abacus > log.txt 2>&1 - -# 运行验证脚本 -python3 ../verify_hamiltonian_convention/check_hamiltonian_convention.py OUT.verify_soc/ -``` - -预期输出： -``` -TEST 1: H(R=0) Hermiticity: H = H^dagger - max|H - H^dagger| = 1.73e-19 - [PASS] - -TEST 2: Off-diagonal phase for m along y - Mean Im(H_{up,down}) = -1.59e-03 - [PASS] Im(H_{up,down}) < 0, consistent with correct convention -``` - ---- - -## 5. PW vs LCAO 对比 - -| 代码路径 | 文件 | 状态 | -|---------|------|------| -| PW | `source/source_pw/module_pwdft/kernels/veff_op.cpp` | **正确** | -| LCAO | `source/source_lcao/module_gint/gint_common.cpp` | **已修复** | - -PW 代码直接实现了正确的公式： -```cpp -// veff_op.cpp: 正确 -sup = out*(in[0]+in[3]) + out1*(in[1] - i*in[2]); // V_0+V_z, V_x-i*V_y -sdown = out1*(in[0]-in[3]) + out*(in[1] + i*in[2]); // V_0-V_z, V_x+i*V_y -``` - -LCAO 代码通过 `merge_hr_part_to_hR()` 间接转换，修复前符号错误。 - ---- - -## 6. 影响范围 - -- 所有 nspin=4 的 LCAO 计算（包括 SOC 和非共线无 SOC） -- CPU 和 GPU 路径均受影响（GPU 代码调用同一 `merge_hr_part_to_hR` 函数） -- regular 和 meta-GGA 均受影响 - -nspin=1,2（共线自旋）不受影响，因为此时非对角元为零。 diff --git a/tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py b/tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py deleted file mode 100644 index 04e7ed7c06e..00000000000 --- a/tests/03_NAO_multik/verify_hamiltonian_convention/check_hamiltonian_convention.py +++ /dev/null @@ -1,246 +0,0 @@ -#!/usr/bin/env python3 -""" -Simple Hamiltonian convention check. -Reads the first R=0 block of H(R) and checks: -1. Hermiticity of R=0 block -2. Off-diagonal phase for m along y -""" -import sys -import re -import numpy as np - -def read_complex_value(s): - """Parse a complex value like '(-1.13e-02,0.00e+00)' or a real number.""" - s = s.strip().strip('()') - if ',' in s: - parts = s.split(',') - return complex(float(parts[0]), float(parts[1])) - else: - return complex(float(s), 0.0) - -def parse_hr_r0(filepath): - """Parse the R=(0,0,0) block from hrs1_nao.csr.""" - with open(filepath) as f: - lines = f.readlines() - - # Parse header - idx = 0 - nbasis = None - nR = None - for i, line in enumerate(lines): - stripped = line.strip() - if stripped.startswith('#') or stripped == '' or stripped.startswith('---'): - continue - if 'number of spin' in stripped or 'spin index' in stripped: - continue - if 'number of localized basis' in stripped: - nbasis = int(stripped.split()[0]) - continue - if 'number of Bravais' in stripped: - nR = int(stripped.split()[0]) - break - - # Find R=(0,0,0) block - found_r0 = False - for i, line in enumerate(lines): - stripped = line.strip() - # Look for "0 0 0 nnz" pattern - parts = stripped.split() - if len(parts) == 4: - try: - rx, ry, rz, nnz = int(parts[0]), int(parts[1]), int(parts[2]), int(parts[3]) - if rx == 0 and ry == 0 and rz == 0 and nnz > 0: - found_r0 = True - break - except ValueError: - continue - - if not found_r0: - print("R=(0,0,0) block not found!") - return None, None, None - - # Now parse values, column indices, row pointers after this line - idx = i + 1 - # Skip "# CSR values" - while idx < len(lines) and '# CSR values' not in lines[idx]: - idx += 1 - idx += 1 - - # Read nnz complex values - values = [] - while len(values) < nnz and idx < len(lines): - line = lines[idx].strip() - idx += 1 - if line.startswith('#'): - break - if not line: - continue - # Parse all complex values on this line - tokens = re.findall(r'\([^)]+\)', line) - if not tokens: - # Try space-separated values - tokens = line.split() - for t in tokens: - try: - values.append(read_complex_value(t)) - except ValueError: - break - if len(values) >= nnz: - break - - # Skip "# CSR column indices" - while idx < len(lines) and '# CSR column indices' not in lines[idx]: - idx += 1 - idx += 1 - - # Read nnz column indices - col_indices = [] - while len(col_indices) < nnz and idx < len(lines): - line = lines[idx].strip() - idx += 1 - if line.startswith('#'): - break - if not line: - continue - for t in line.split(): - try: - col_indices.append(int(t)) - except ValueError: - break - if len(col_indices) >= nnz: - break - - # Skip "# CSR row pointers" - while idx < len(lines) and '# CSR row pointers' not in lines[idx]: - idx += 1 - idx += 1 - - # Read nbasis+1 row pointers - row_ptrs = [] - while len(row_ptrs) < nbasis + 1 and idx < len(lines): - line = lines[idx].strip() - idx += 1 - if line.startswith('#'): - break - if not line: - continue - for t in line.split(): - try: - row_ptrs.append(int(t)) - except ValueError: - break - if len(row_ptrs) >= nbasis + 1: - break - - # Build dense matrix - H = np.zeros((nbasis, nbasis), dtype=complex) - for irow in range(nbasis): - if irow >= len(row_ptrs): - break - start = row_ptrs[irow] - end = row_ptrs[irow + 1] if irow + 1 < len(row_ptrs) else nnz - for k in range(start, min(end, nnz, len(values), len(col_indices))): - icol = col_indices[k] - if 0 <= icol < nbasis: - H[irow, icol] = values[k] - - return nbasis, H, nnz - -def main(): - if len(sys.argv) < 2: - print("Usage: python3 check_hamiltonian_convention.py ") - sys.exit(1) - - out_dir = sys.argv[1] - hr_file = None - dmr_file = None - - import os - for f in os.listdir(out_dir): - if f.startswith('hrs') and f.endswith('.csr'): - hr_file = os.path.join(out_dir, f) - if f.startswith('dmrs') and f.endswith('.csr'): - dmr_file = os.path.join(out_dir, f) - - if not hr_file: - print("ERROR: hrs*_nao.csr not found") - sys.exit(1) - - print(f"Parsing H(R=0) from: {hr_file}") - nbasis, H0, nnz = parse_hr_r0(hr_file) - if H0 is None: - print("Failed to parse H(R=0)") - sys.exit(1) - - print(f" nbasis = {nbasis}, nnz = {nnz}") - N = nbasis // 2 # number of spatial orbitals per spin - - # Test 1: Hermiticity of R=0 block - print("\n" + "=" * 70) - print("TEST 1: H(R=0) Hermiticity: H = H^dagger") - print("=" * 70) - herm_err = np.max(np.abs(H0 - H0.conj().T)) - print(f" max|H - H^dagger| = {herm_err:.2e}") - if herm_err < 1e-10: - print(" [PASS]") - else: - print(" [FAIL]") - - # Test 2: Off-diagonal spinor block phase - print("\n" + "=" * 70) - print("TEST 2: Off-diagonal phase for m along y") - print(" For m||+y: H_{up,down} should have Im < 0") - print(" (i.e., H(2i, 2j+1) should be -i*|B_y|)") - print("=" * 70) - - # Extract the up-down block: rows 0,2,4,..., cols 1,3,5,... - H_updown = np.zeros((N, N), dtype=complex) - for i in range(N): - for j in range(N): - H_updown[i, j] = H0[2*i, 2*j+1] - - # Check the sign of imaginary part - nonzero_mask = np.abs(H_updown) > 1e-10 - if np.any(nonzero_mask): - imag_parts = H_updown[nonzero_mask].imag - mean_imag = np.mean(imag_parts) - n_neg = np.sum(imag_parts < 0) - n_pos = np.sum(imag_parts > 0) - print(f" Non-zero H_{{up,down}} elements: {np.sum(nonzero_mask)}") - print(f" Mean Im(H_{{up,down}}) = {mean_imag:.6e}") - print(f" Negative: {n_neg}, Positive: {n_pos}") - if mean_imag < -1e-10: - print(" [PASS] Im(H_{{up,down}}) < 0, consistent with correct convention") - elif mean_imag > 1e-10: - print(" [FAIL] Im(H_{{up,down}}) > 0, indicates complex conjugation error") - else: - print(" [INFO] Im(H_{{up,down}}) ~ 0, magnetization may not be along y") - else: - print(" No significant off-diagonal elements found") - - # Test 3: Diagonal spinor block structure - print("\n" + "=" * 70) - print("TEST 3: Diagonal blocks H_{{up,up}} and H_{{down,down}}") - print(" For m||y with no SOC: H_{{up,up}} = H_{{down,down}} (real)") - print("=" * 70) - H_upup = np.zeros((N, N), dtype=complex) - H_downdown = np.zeros((N, N), dtype=complex) - for i in range(N): - for j in range(N): - H_upup[i, j] = H0[2*i, 2*j] - H_downdown[i, j] = H0[2*i+1, 2*j+1] - - diag_diff = np.max(np.abs(H_upup - H_downdown)) - upup_imag = np.max(np.abs(H_upup.imag)) - dndn_imag = np.max(np.abs(H_downdown.imag)) - print(f" max|H_{{up,up}} - H_{{down,down}}| = {diag_diff:.2e}") - print(f" max|Im(H_{{up,up}})| = {upup_imag:.2e}") - print(f" max|Im(H_{{down,down}})| = {dndn_imag:.2e}") - - if diag_diff < 0.1: - print(" [PASS] Diagonal blocks are nearly equal (expected for m perp to z)") - else: - print(" [INFO] Diagonal blocks differ (may indicate m has z-component or SOC)") - -if __name__ == '__main__': - main()