|
1 | | -# LCAO 力/应力中 DMR 配对约定推导(DM/DMK/DMR 链路) |
2 | | - |
3 | | -- 日期:2026-08-09 |
4 | | -- 相关代码:`source/source_estate/module_dm/`(DMK/DMR)、 |
5 | | - `source/source_hamilt/module_hcontainer/func_folding.cpp`(folding_HR)、 |
6 | | - `source/source_lcao/module_operator_lcao/`(overlap/ekinetic/nonlocal 力与应力) |
7 | | -- 配套文档:`docs/dm_dmk_dmr_investigation.md`(排查)、 |
8 | | - `docs/dm_dmk_dmr_action_plan.md`(整改与测试方案) |
9 | | - |
10 | | -本文档锁定 DMK→DMR 的傅里叶符号约定,以及力/应力计算中实空间密度矩阵 |
11 | | -与两中心积分(含 nonlocal 的 R2−R1 配对)的收缩对应关系,作为代码注释 |
12 | | -(`density_matrix.cpp` 中 "See Sec. 3 ...")所指的推导依据。 |
13 | | - |
14 | | -## 1. 指标与 R 格矢约定 |
15 | | - |
16 | | -- 波函数系数 `C_{μn}(k)`:μ = 轨道(行),n = 能带(列),存储为 |
17 | | - `psi::Psi` 的 `wfc(ib, iw)`。 |
18 | | -- 密度矩阵元素 `D_{μν}`:μ = 行(bra)轨道,ν = 列(ket)轨道。 |
19 | | -- 格矢 R:`HContainer::AtomPair(iat1, iat2, R)` 的矩阵块为 |
20 | | - `O_{μν}(R)`,μ∈iat1(0 胞)、ν∈iat2(R 胞);格矢满足 |
21 | | - `dtau = τ(iat2) + R·L − τ(iat1)`(`UnitCell::cal_dtau`)。 |
22 | | -- k、R 均为**直接坐标**,相位 `k·R` 前乘 `2π`。 |
23 | | - |
24 | | -## 2. k 空间与实空间傅里叶对 |
25 | | - |
26 | | -- 正变换(算符,`folding_HR`,`func_folding.cpp`): |
| 1 | +# DMR pairing conventions in LCAO force/stress (DM/DMK/DMR chain) |
| 2 | + |
| 3 | +- Date: 2026-08-09 |
| 4 | +- Related code: `source/source_estate/module_dm/` (DMK/DMR), |
| 5 | + `source/source_hamilt/module_hcontainer/func_folding.cpp` (folding_HR), |
| 6 | + `source/source_lcao/module_operator_lcao/` (overlap/ekinetic/nonlocal |
| 7 | + force and stress) |
| 8 | +- This document fixes the DMK->DMR Fourier sign convention and the |
| 9 | + contraction mapping between the real-space density matrix and the two-center |
| 10 | + integrals (including the nonlocal R2-R1 pairing) used by the force/stress |
| 11 | + code. It is the derivation referenced by the code comments (the |
| 12 | + "See Sec. 3 ..." comments in `density_matrix.cpp`). |
| 13 | + |
| 14 | +## 1. Index and R-lattice conventions |
| 15 | + |
| 16 | +- Wave-function coefficients `C_{μn}(k)`: μ = orbital (row), n = band |
| 17 | + (column), stored as `psi::Psi`'s `wfc(ib, iw)`. |
| 18 | +- Density-matrix elements `D_{μν}`: μ = row (bra) orbital, ν = column (ket) |
| 19 | + orbital. |
| 20 | +- Lattice vector R: the matrix block of `HContainer::AtomPair(iat1, iat2, R)` |
| 21 | + is `O_{μν}(R)` with μ ∈ iat1 (home cell) and ν ∈ iat2 (cell R); the vector |
| 22 | + satisfies `dtau = τ(iat2) + R·L − τ(iat1)` (`UnitCell::cal_dtau`). |
| 23 | +- k and R are both in **direct coordinates**; the phase `k·R` carries a |
| 24 | + factor of `2π`. |
| 25 | + |
| 26 | +## 2. k-space and real-space Fourier pair |
| 27 | + |
| 28 | +- Forward transform (operators, `folding_HR`, `func_folding.cpp`): |
27 | 29 |
|
28 | 30 | ``` |
29 | 31 | O(k) = Σ_R e^{+ikR} O(R) |
30 | 32 | ``` |
31 | 33 |
|
32 | | -- 逆变换(密度矩阵,`DensityMatrix_Tools::cal_DMR*`,`density_matrix.cpp`): |
| 34 | +- Inverse transform (density matrix, `DensityMatrix_Tools::cal_DMR*`, |
| 35 | + `density_matrix.cpp`): |
33 | 36 |
|
34 | 37 | ``` |
35 | 38 | D(R) = Σ_k e^{-ikR} DMK(k) |
36 | 39 | ``` |
37 | 40 |
|
38 | | -- 归一化:k 点权重 `w_k` 已内嵌在 `DMK`(`wg = w_k·occ`, |
39 | | - `occupy.cpp`),因此 `cal_DMR` **不**再乘 `1/Nk`。 |
| 41 | +- Normalization: the k-point weights `w_k` are embedded in `DMK` |
| 42 | + (`wg = w_k·occ`, `occupy.cpp`), so `cal_DMR` must **not** multiply by |
| 43 | + `1/Nk` again. |
40 | 44 | - `DMK(μ,ν;k) = Σ_n f_nk C*_{μn}(k) C_{νn}(k) = (C f C†)^T = D_std^T` |
41 | | - (`cal_dm_psi.cpp`,`zgemm('N','T')` + 先取共轭)。 |
| 45 | + (`cal_dm_psi.cpp`, `zgemm('N','T')` on pre-conjugated coefficients). |
42 | 46 |
|
43 | | -## 3. 正逆变换互为逆(Sec. 3 约定) |
| 47 | +## 3. Forward and inverse transforms are mutual inverses (Sec. 3 convention) |
44 | 48 |
|
45 | | -对同一套 R 定义,`e^{+ikR}`(folding_HR)与 `e^{-ikR}`(cal_DMR)构成 |
46 | | -自洽的正逆傅里叶对。若把 `density_matrix.cpp` 中的 `-sinp` 改回 `+sinp`, |
47 | | -则 `D(R)` 实际存的是正变换 `Σ_k e^{+ikR} DMK(k)`,与 `folding_HR` 不再 |
48 | | -互为逆变换,非对角/非对称消费方(如 nspin=4 nonlocal 力)即出错。 |
49 | | -该符号由单元测试 `T1_fourier_round_trip` 锁定(见 |
50 | | -`docs/dm_dmk_dmr_action_plan.md` §T1)。 |
| 49 | +For the same set of R vectors, `e^{+ikR}` (folding_HR) and `e^{-ikR}` |
| 50 | +(cal_DMR) form a self-consistent forward/inverse Fourier pair. If the |
| 51 | +`-sinp` in `density_matrix.cpp` were changed back to `+sinp`, then `D(R)` |
| 52 | +would actually store the forward transform `Σ_k e^{+ikR} DMK(k)`, which is no |
| 53 | +longer the inverse of `folding_HR`; non-diagonal / non-symmetric consumers |
| 54 | +(such as the nspin=4 nonlocal force) then become wrong. This sign is locked |
| 55 | +by the unit test `T1_fourier_round_trip` (see below). |
51 | 56 |
|
52 | | -## 4. 力/应力的一般形式 |
| 57 | +## 4. General form of the force/stress |
53 | 58 |
|
54 | | -Feynman–Hellmann 型受力(overlap/ekinetic/nonlocal 等): |
| 59 | +The Feynman–Hellmann-type forces (overlap/ekinetic/nonlocal etc.) read: |
55 | 60 |
|
56 | 61 | ``` |
57 | 62 | F = Σ_k Tr[ D(k) · ∂O(k)/∂τ ] |
58 | 63 | = Σ_k Σ_{μν} D_{μν}(k) ∂O_{νμ}(k)/∂τ |
59 | 64 | ``` |
60 | 65 |
|
61 | | -代入 `∂O(k) = Σ_R e^{+ikR} ∂O(R)`: |
| 66 | +Substituting `∂O(k) = Σ_R e^{+ikR} ∂O(R)`: |
62 | 67 |
|
63 | 68 | ``` |
64 | 69 | F = Σ_R Σ_{μν} ∂O_{νμ}(R) D_{μν}(−R) (*) |
65 | | - 其中 D(−R) := Σ_k e^{+ikR} DMK(k) |
| 70 | + with D(−R) := Σ_k e^{+ikR} DMK(k) |
66 | 71 | ``` |
67 | 72 |
|
68 | | -代码存储的是 `D(R) = Σ_k e^{-ikR} DMK(k)`,两者通过 §5 的闭合迹论证等价。 |
| 73 | +The code stores `D(R) = Σ_k e^{-ikR} DMK(k)`; the two are equivalent through |
| 74 | +the closed-trace argument of Sec. 5. |
69 | 75 |
|
70 | | -## 5. 闭合迹共轭保护(Sec. 5 配对) |
| 76 | +## 5. Closed-trace conjugate protection (Sec. 5 pairing) |
71 | 77 |
|
72 | | -代码在两中心积分路径中,对每个 `(iat1,iat2,R)` 同指标收缩: |
| 78 | +In the two-center integral path the code contracts same-index pairs for every |
| 79 | +`(iat1,iat2,R)`: |
73 | 80 |
|
74 | 81 | ``` |
75 | 82 | W_code = Σ_{μνR} D(μ,ν;R) ∂O(μ,ν;R) |
76 | 83 | ``` |
77 | 84 |
|
78 | | -把 `D(R) = Σ_k e^{-ikR} DMK(k)` 代入: |
| 85 | +Substituting `D(R) = Σ_k e^{-ikR} DMK(k)`: |
79 | 86 |
|
80 | 87 | ``` |
81 | 88 | W_code = Σ_k Σ_{μν} DMK(μ,ν;k) [Σ_R e^{-ikR} ∂O(μ,ν;R)] |
82 | 89 | = Σ_k Tr( DMK(k) · ∂O'(−k) ) |
83 | 90 | ``` |
84 | 91 |
|
85 | | -其中 `∂O'(−k) := Σ_R e^{-ikR} ∂O(R)` 是 folding 在 −k 处的取值。 |
86 | | -若 `∂O` 在 HContainer 中**全方向成对**(对每个 `(iat1,iat2,R)` 存在 |
87 | | -`(iat2,iat1,−R)`,且 `∂O(iat2,iat1,−R) = ∂O(iat1,iat2,R)†`),则 |
88 | | -`∂O'(k) = Σ_R e^{+ikR} ∂O(R)` 逐 k 厄米,`∂O'(−k) = ∂O'(k)†`。于是 |
| 92 | +where `∂O'(−k) := Σ_R e^{-ikR} ∂O(R)` is the folded operator evaluated at |
| 93 | +`−k`. If `∂O` is stored in the HContainer with **full-direction pairing** |
| 94 | +(for every `(iat1,iat2,R)` there is `(iat2,iat1,−R)` with |
| 95 | +`∂O(iat2,iat1,−R) = ∂O(iat1,iat2,R)†`), then `∂O'(k) = Σ_R e^{+ikR} ∂O(R)` |
| 96 | +is Hermitian for each k and `∂O'(−k) = ∂O'(k)†`. Hence |
89 | 97 |
|
90 | 98 | ``` |
91 | 99 | W_code = Σ_k w_k Tr( DMK(k) · ∂O'(k)† ) |
92 | 100 | = Σ_k w_k conj( Tr( DMK(k) ∂O'(k) ) ) |
93 | | - = Σ_k w_k Tr( DMK(k) ∂O'(k) ) [厄米矩阵对的迹为实数] |
| 101 | + = Σ_k w_k Tr( DMK(k) ∂O'(k) ) [trace of Hermitian pair is real] |
94 | 102 | = Σ_k w_k Re Tr( DMK(k) ∂O'(k) ) |
95 | 103 | ``` |
96 | 104 |
|
97 | | -即**总和受共轭保护**,而非"逐项相等"。保护成立的三条件: |
| 105 | +That is, the **total is protected by the conjugate** rather than by |
| 106 | +"term-by-term equality". The protection holds under three conditions: |
98 | 107 |
|
99 | | -1. 收缩是**全轨道闭合 Frobenius 迹**(不按元素取用); |
100 | | -2. 配对算符**逐 k 厄米**(HContainer 全方向成对存储); |
101 | | -3. 只取**实部**。 |
| 108 | +1. The contraction is a **closed Frobenius trace over all orbitals** |
| 109 | + (no element-wise consumption); |
| 110 | +2. The paired operator is **Hermitian per k** (HContainer full-direction |
| 111 | + paired storage); |
| 112 | +3. Only the **real part** is kept. |
102 | 113 |
|
103 | | -破坏任一条件的消费方即暴露于 O(1) 配对误差。 |
| 114 | +Any consumer that breaks one of these conditions is exposed to an O(1) |
| 115 | +pairing error. |
104 | 116 |
|
105 | | -### 5.1 overlap / ekinetic(两中心积分路径) |
| 117 | +### 5.1 overlap / ekinetic (two-center integral path) |
106 | 118 |
|
107 | | -`operator_fs_utils.hpp` 的 `cal_force_stress_2center` 对 |
108 | | -`(iat1,iat2,R)` 直接取 `D(μ,ν;R) ∂O(μ,ν;R)` 同指标收缩,且对 |
109 | | -`(iat2,iat1,−R)` 同样各计一次(全 R 集迭代),故 |
110 | | -`finalize_force_stress(force_factor = 1.0)`(无需因子 2)。 |
| 119 | +`cal_force_stress_2center` in `operator_fs_utils.hpp` contracts |
| 120 | +`D(μ,ν;R) ∂O(μ,ν;R)` with same indices for every `(iat1,iat2,R)`, and |
| 121 | +`(iat2,iat1,−R)` is visited once as well (full R-set iteration), hence |
| 122 | +`finalize_force_stress(force_factor = 1.0)` (no factor 2 needed). |
111 | 123 |
|
112 | | -### 5.2 nonlocal(β 投影,R2−R1 配对) |
| 124 | +### 5.2 nonlocal (β projection, R2−R1 pairing) |
113 | 125 |
|
114 | | -`nonlocal_fs.cpp::cal_force_stress` 以中心原子 `iat0` 为基准枚举两个 |
115 | | -邻原子:`iat1`(胞 `R1`)、`iat2`(胞 `R2`),取用的 DMR 块是 |
| 126 | +`nonlocal_fs.cpp::cal_force_stress` enumerates two neighbor atoms relative to |
| 127 | +the central atom `iat0`: `iat1` (cell `R1`) and `iat2` (cell `R2`), and |
| 128 | +consumes the DMR block |
116 | 129 |
|
117 | 130 | ``` |
118 | 131 | dmR->find_matrix(iat1, iat2, R2 − R1) |
119 | 132 | ``` |
120 | 133 |
|
121 | | -即配对结构为 `R_vector = R2 − R1`(`nonlocal_fs.cpp` 中 |
122 | | -`R_index2 − R_index1`),不是简单同 R 闭合迹。该路径的受力公式为 |
| 134 | +i.e. the pairing structure is `R_vector = R2 − R1` (`R_index2 − R_index1` in |
| 135 | +`nonlocal_fs.cpp`), which is not a simple same-R closed trace. The force on |
| 136 | +this path reads |
123 | 137 |
|
124 | 138 | ``` |
125 | 139 | F(iat0) = Σ_{iat1,iat2} <∂β_{iat1,R1}/∂τ_iat0 | D_{iat1,iat2}(R2−R1) | β_{iat2,R2}> |
126 | 140 | ``` |
127 | 141 |
|
128 | | -逐项按 `(iat1,iat2,R2−R1)` 取用复数 DMR 块(nspin=4 时由 |
129 | | -`cal_DMR_full` 提供,不丢虚部)。该路径**不满足** §5 的同 R 闭合迹 |
130 | | -保护条件,其正确性依赖: |
131 | | -1. `cal_DMR_full` 采用 `e^{-ikR}`(与 `folding_HR` 的 `e^{+ikR}` 互逆); |
132 | | -2. DMR 块在 `(iat1,iat2,R2−R1)` 上的取值与 β 投影的指标顺序一致; |
133 | | -3. 数值上由 SOC 有限差分力测试锁定(见 `docs/dm_dmk_dmr_action_plan.md` §T4)。 |
134 | | - |
135 | | -若未来引入半集存储优化(只存 `(iat1,iat2,R)` 不存 `(iat2,iat1,−R)`), |
136 | | -§5 的保护论证与 `force_factor=1.0` 约定同时失效,须恢复 factor=2 并 |
137 | | -重跑 T2/T3/T8 守卫。 |
| 142 | +Each term consumes the complex DMR block element-wise by |
| 143 | +`(iat1,iat2,R2−R1)` (for nspin=4 the block is provided by `cal_DMR_full` |
| 144 | +without dropping the imaginary part). This path **does not** satisfy the |
| 145 | +same-R closed-trace protection of Sec. 5; its correctness relies on: |
| 146 | + |
| 147 | +1. `cal_DMR_full` using `e^{-ikR}` (mutually inverse with `folding_HR`'s |
| 148 | + `e^{+ikR}`); |
| 149 | +2. the DMR block taken at `(iat1,iat2,R2−R1)` matching the index order of |
| 150 | + the β projection; |
| 151 | +3. numerically locked by the SOC finite-difference force integration test |
| 152 | + (`tests/integrate/240_NO_KP_15_SO_FD`). |
| 153 | + |
| 154 | +If a half-set storage optimization is introduced in the future (storing only |
| 155 | +`(iat1,iat2,R)` without `(iat2,iat1,−R)`), the Sec. 5 protection argument |
| 156 | +and the `force_factor = 1.0` convention both break; the factor must be |
| 157 | +restored to 2 and the T2/T3/T8 guards re-run. |
0 commit comments