MD and SCF results are not same with same STRU

[dyzheng]

Describe the bug

in #6991 , I changed test case md_out_sync and find this error.

Expected behavior

No response

To Reproduce

No response

Environment

No response

Additional Context

No response

Task list for Issue attackers (only for developers)

• Verify the issue is not a duplicate.
• Describe the bug.
• Steps to reproduce.
• Expected behavior.
• Error message.
• Environment details.
• Additional context.
• Assign a priority level (low, medium, high, urgent).
• Assign the issue to a team member.
• Label the issue with relevant tags.
• Identify possible related issues.
• Create a unit test or automated test to reproduce the bug (if applicable).
• Fix the bug.
• Test the fix.
• Update documentation (if necessary).
• Close the issue and inform the reporter (if applicable).

[dyzheng]

---

问题描述

现象：相同结构和输入参数，SCF 和 MD 计算的第一个电子步：

• E_band 相同 ✓
• E_hartree 相同 ✓
• E_xc 不同 ✗

---

能量公式（代码排查结果）

总能量：

// fp_energy.cpp:22-27
etot = eband + deband + (etxc - etxcc) + hartree_energy + descf + ...

各项定义：

能量项	公式	使用的 rho	计算时机
E_band	Σ wg × ε	无依赖	psiToRho 之前
deband	-∫ ρ × V_eff dr	rho_final	Symmetry_rho 之后
E_xc	∫ ρ ε_xc(ρ) dr	rho_final	update_pot
E_hartree	Σ_G |ρ(G)|² / G²	rho_final	update_pot

关键发现：

• deband 和 E_xc 都使用 Symmetry_rho 之后的电荷密度
• deband = -∫ ρ × V_eff dr，其中 V_eff = V_H + V_xc + ...
• 所以 deband 也包含与 XC 相关的贡献

---

电荷密度状态追踪

iter_init():
    save_rho_before_sum_band() → rho_save = rho_in

hamilt2density_single():
    psiToRho(ψ) → rho = rho_out
        ↓
    Symmetry_rho::begin():  ← 【关键差异点】
        SCF (symm_flag=1): rho = rho_symm  ← 修改
        MD (symm_flag=0): rho = rho_out    ← 不修改
        ↓
    rho_final 确定
        ↓
    deband = -∫ rho_final × V_eff dr  ← 使用 rho_final

iter_finish():
    update_pot():
        E_hartree = f(rho_final)  ← 使用 rho_final
        E_xc = f(rho_final)       ← 使用 rho_final

---

根因

symmetry 参数默认值不同：

计算类型	symmetry	symm_flag	Symmetry_rho
SCF	1	1	执行，rho 被对称化
MD	0	0	不执行，rho 保持原样

代码证据：

// read_input_item_system.cpp:145-162
if (para.input.symmetry == "default") {
    para.input.symmetry = "1";  // SCF 默认
}
if (para.input.calculation == "md") {
    para.input.symmetry = "0";  // MD 强制
}

---

完整影响链

symmetry 参数差异
    ↓
Symmetry_rho 执行条件不同
    ↓
rho_final 不同：
    SCF: rho_out → rho_symm (对称化)
    MD: rho_out (保持原样)
    ↓
deband 不同（使用 rho_final）
E_xc 不同（使用 rho_final）

---

为什么 E_band 相同

• E_band = Σ wg × ε（波函数本征值加权求和）
• 计算时机：calEBand() 在 psiToRho() 之前
• 不依赖电荷密度

---

为什么 E_hartree 相同但 E_xc 不同

E_hartree 公式：
$$E_H = \frac{1}{2} \sum_G \frac{4\pi |\rho(G)|^2}{G^2}$$

• G 空间积分
• 对称化保持 |ρ(G)|² 总量
• 对局域分布不敏感

E_xc 公式：
$$E_{xc} = \int \rho(r) \epsilon_{xc}(\rho(r)) dr$$

• r 空间局域积分
• 对称化修改 ρ(r) 的局域分布
• 对局域分布敏感

---

解决方案

# INPUT 文件
symmetry = 0  # SCF 也关掉对称性，与 MD 一致

---

关键代码位置

文件	功能
read_input_item_system.cpp:145-162	symmetry 默认值设置
symmetry_rho.cpp:19-33	rho_symmetry 实现
esolver_ks_lcao.cpp:774-782	Symmetry_rho 调用和 deband 计算
elecstate_energy.cpp:149-151	deband 公式
pot_xc.cpp:38-42	etxc 计算
hsolver_lcao.cpp:78-82	calEBand 时机

---