sd cg diis

claude/plan/sd_diis_refactor_plan.md

@@ -0,0 +1,162 @@
+# SD 和 DIIS 优化算法重构计划
+
+## 目标
+将 SD (Steepest Descent) 和 DIIS 优化算法重构为类式设计，与 LBFGS/RFO 保持一致的代码风格。
+
+---
+
+## 当前问题分析
+
+### 1. SD.py 问题
+- **函数式设计**：直接是 `SD()` 函数，不是类
+- **参数硬编码**：`max_step_size=0.2`, `maxiterations=128` 作为函数参数
+- **没有继承 JobABC**：无法使用 `_init_params()` 等工具方法
+- **单位转换不一致**：使用 `g_au` 转换，而 LBFGS 不使用
+- **DIIS 集成粗糙**：每 10 次迭代在循环内部 `from .DIIS import DIIS`
+- **没有轨迹输出**：缺少 `write_xyz` 功能
+
+### 2. DIIS.py 问题
+- **函数式设计**：`DIIS()` 是函数而非类
+- **职责不清**：DIIS 应该是加速器组件，不是独立优化器
+- **单位转换问题**：同样使用 `g_au`
+
+### 3. optimization.py 问题
+- **SD 调用不一致**：`SD(self.atoms, output=self.output)` 没有传 `paras` 参数
+
+---
+
+## 重构设计
+
+### 1. DIISAccelerator 类（加速器模式）
+
+将 DIIS 重新定位为**加速器组件**，可被任何优化器复用：
+
+```python
+@dataclass
+class DIISParams:
+    memory: int = 10          # 历史向量数量
+    min_vectors: int = 3      # 最少向量数才能外推
+
+class DIISAccelerator:
+    def __init__(self, params: DIISParams = None)
+    def store(self, x: np.ndarray, error: np.ndarray) -> None
+    def can_extrapolate(self) -> bool
+    def extrapolate(self) -> Optional[np.ndarray]
+    def reset(self) -> None
+```
+
+### 2. SD 类设计（继承 JobABC）
+
+```python
+@dataclass
+class SDParams:
+    max_step: float = 0.2           # 最大步长 (Angstrom)
+    max_iter: int = 256             # 最大迭代次数
+    diis_enabled: bool = True       # 是否启用 DIIS 加速
+    diis_interval: int = 10         # DIIS 触发间隔
+    diis_memory: int = 10           # DIIS 历史向量数
+    write_traj: bool = False        # 是否写轨迹
+    traj_every: int = 1             # 轨迹写入间隔
+    verbose: int = 1                # 日志详细程度
+
+class SD(JobABC):
+    def __init__(self, atoms, output, paras=None)
+    def run(self) -> Atoms
+    def _clip_step(self, step) -> np.ndarray
+    def _try_diis_acceleration(self, forces) -> Optional[np.ndarray]
+    def _build_iter_message(...) -> List[str]
+    def _log_iter(self, info_message) -> None
+    def _check_convergence(self) -> bool
+```
+
+---
+
+## 需要修改的文件
+
+| 文件 | 修改内容 |
+|------|---------|
+| `maple/function/dispatcher/optimization/algorithm/DIIS.py` | 重构为 DIISAccelerator 类 + 保留旧函数兼容 |
+| `maple/function/dispatcher/optimization/algorithm/SD.py` | 完全重写为类式设计 |
+| `maple/function/dispatcher/optimization/algorithm/__init__.py` | 更新导出 |
+| `maple/function/dispatcher/optimization/optimization.py` | 更新 SD 调用方式 |
+
+---
+
+## 详细实现步骤
+
+### Step 1: 重构 DIIS.py
+
+1. 创建 `DIISParams` dataclass
+2. 创建 `DIISAccelerator` 类
+   - `store()`: 存储位置和误差向量
+   - `extrapolate()`: 执行 DIIS 外推（构建 B 矩阵，求解线性系统）
+   - `reset()`: 重置历史
+3. 保留旧的 `OptimizationStorage` 和 `DIIS()` 函数（标记为 deprecated）
+
+### Step 2: 重构 SD.py
+
+1. 创建 `SDParams` dataclass（参考 LBFGSParams）
+2. 创建 `SD` 类继承 `JobABC`
+3. 实现核心方法：
+   - `__init__()`: 初始化参数，创建 DIISAccelerator
+   - `_clip_step()`: 限制步长
+   - `_try_diis_acceleration()`: 尝试 DIIS 加速
+   - `_build_iter_message()`: 构建迭代日志（与 LBFGS 格式一致）
+   - `_check_convergence()`: 检查收敛条件
+   - `run()`: 主优化循环
+4. 添加 `write_xyz()` 轨迹输出
+5. **移除 `g_au` 单位转换**，与 LBFGS 保持一致
+
+### Step 3: 更新 __init__.py
+
+```python
+from .DIIS import DIISAccelerator, DIISParams, OptimizationStorage, DIIS
+from .SD import SD, SDParams
+```
+
+### Step 4: 更新 optimization.py
+
+```python
+elif self.commandcontrol.get('method').lower() == 'sd':
+    from .algorithm import SD
+    opt = SD(self.atoms, output=self.output, paras=self.commandcontrol)
+    return opt.run()
+```
+
+---
+
+## 关键设计决策
+
+### 1. DIIS 作为加速器而非独立优化器
+- DIIS 本身不计算搜索方向，只做位置外推
+- 作为组件可被 SD、CG 等优化器复用
+
+### 2. 移除 `g_au` 单位转换
+- LBFGS/RFO 直接使用 ASE 返回的单位 (eV, eV/Å)
+- SD/DIIS 应该保持一致
+
+### 3. 保留向后兼容性
+- 保留旧的 `DIIS()` 函数和 `OptimizationStorage` 类
+- 使用 DeprecationWarning 提示迁移
+
+### 4. 统一的收敛检查
+- 四项标准：max_f, rms_f, max_dp, rms_dp
+- 与 LBFGS 使用相同的阈值属性
+
+---
+
+## 验证方案
+
+1. **运行测试**：使用简单分子（如 H2O）测试 SD 优化
+2. **对比收敛**：验证 DIIS 加速确实提高收敛速度
+3. **检查输出**：确保日志格式与 LBFGS 一致
+4. **回归测试**：确保现有功能不受影响
+
+---
+
+## 预期输出文件
+
+优化完成后将生成：
+- `{base}_traj.xyz` - 完整轨迹
+- `{base}_opt.xyz` - 最终优化结构
+- `{base}_opt_traj.xyz` - 详细轨迹（当 write_traj=True）

claude/plan/sdcg_fusion_plan.md

@@ -0,0 +1,197 @@
+# SDCG Fusion Optimizer Plan
+
+## Overview
+Fuse Steepest Descent (SD) and Conjugate Gradient (CG) into a single phased optimizer `SDCG`, with optional GDIIS acceleration. SD provides robust initial descent; CG (PRP+ variant) provides faster convergence near the minimum.
+
+## File Changes
+
+### 1. NEW: `maple/function/dispatcher/optimization/algorithm/SDCG.py`
+The main implementation file, replacing SD.py as the active optimizer.
+
+**SDCGParams dataclass:**
+```python
+@dataclass
+class SDCGParams:
+    # General
+    max_step: float = 0.2
+    max_iter: int = 256
+    write_traj: bool = False
+    traj_every: int = 1
+    verbose: int = 1
+
+    # Phase control
+    sd_enabled: bool = True         # Enable SD phase
+    cg_enabled: bool = True         # Enable CG phase
+    sd_max_iter: int = 50           # Max SD iterations before forced CG switch
+    cg_switch_fmax: float = 0.0     # Switch to CG when max_f < threshold (0.0 = auto)
+
+    # CG parameters
+    cg_restart_threshold: float = 0.2  # Powell restart: |f_k·f_{k-1}| >= threshold * ||f_k||^2
+    cg_beta_method: str = "prp+"       # CG variant (prp+ only for now)
+
+    # GDIIS parameters
+    diis_enabled: bool = True
+    diis_store_every: int = 5
+    diis_min_snapshots: int = 3
+    diis_memory: int = 6
+```
+
+**SDCG class (inherits JobABC):**
+- Phase state machine: `_phase = "sd"` or `"cg"`
+- SD phase: identical to current SD.py logic (step = max_step * forces, clip)
+- CG phase: PRP+ conjugate gradient with Powell restart
+- GDIIS: works in both phases, reset on phase transition
+- Barzilai-Borwein step scale estimation (carried over from SD)
+
+**PRP+ CG Algorithm:**
+```
+d_0 = f_0  (first CG direction = steepest descent)
+For k >= 1:
+    beta = max(dot(f_k, f_k - f_{k-1}) / dot(f_{k-1}, f_{k-1}), 0)  # PRP+
+
+    # Powell restart condition
+    if |dot(f_k, f_{k-1})| >= 0.2 * dot(f_k, f_k):
+        beta = 0  # restart to SD direction
+
+    d_k = f_k + beta * d_{k-1}
+
+    # Step size: use max_step scaling (same as SD), clip per-atom
+    step = max_step * d_k / max(|d_k|)  # normalize so max displacement = max_step
+    x_{k+1} = x_k + step
+```
+
+**SD → CG Transition Logic:**
+1. Auto mode (`cg_switch_fmax = 0.0`): threshold = `0.5 * initial_max_f` (computed at iter 0)
+2. Explicit mode: user sets `cg_switch_fmax` to desired value
+3. Iteration-based: always switch after `sd_max_iter` SD steps
+4. Transition triggers on whichever condition is met first
+5. On transition: reset GDIIS, log phase change, start CG from current forces
+
+**User Control Modes:**
+- `#opt(method=sd)` → `sd_enabled=True, cg_enabled=False` (SD only)
+- `#opt(method=cg)` → `sd_enabled=False, cg_enabled=True` (CG only)
+- `#opt(method=sdcg)` → `sd_enabled=True, cg_enabled=True` (default fusion)
+- `diis_enabled=false` → disable GDIIS in any mode
+
+### 2. MODIFY: `maple/function/dispatcher/optimization/algorithm/SD.py`
+Convert to backward-compatibility stub:
+```python
+from .SDCG import SDCG, SDCGParams
+
+# Backward compatibility aliases
+SD = SDCG
+SDParams = SDCGParams
+```
+Keep the legacy function-based `DIIS()` reference in DIIS.py unchanged.
+
+### 3. MODIFY: `maple/function/dispatcher/optimization/algorithm/__init__.py`
+```python
+from .LBFGS import LBFGS, LBFGSParams
+from .RFO import RFO, RFOParams
+from .DIIS import DIIS, DIISAccelerator, DIISParams, OptimizationStorage
+from .SDCG import SDCG, SDCGParams
+from .SD import SD, SDParams  # backward compat
+```
+
+### 4. MODIFY: `maple/function/dispatcher/optimization/optimization.py`
+Add routing for `sdcg` and `cg`, update `sd` routing:
+```python
+elif method in ('sd', 'sdcg', 'cg'):
+    from .algorithm import SDCG
+    opt = SDCG(self.atoms, output=self.output, paras=self.commandcontrol)
+    return opt.run()
+```
+
+### 5. MODIFY: `maple/function/read/command_control.py`
+Add `"sdcg"` to IMPLEMENTATION_MAP:
+```python
+"opt": {"lbfgs", "rfo", "sd", "cg", "sdcg", ""},
+```
+
+### 6. NEW: `example/opt/sdcg/inp1.inp`
+Test input for SDCG fusion mode (same molecule as SD tests):
+```
+#model=uma
+#opt(method=sdcg)
+#device=gpu0
+```
+
+### 7. NEW: `example/opt/cg/inp1.inp`
+Test input for CG-only mode:
+```
+#model=uma
+#opt(method=cg)
+#device=gpu0
+```
+
+## Implementation Details
+
+### Phase Transition in `run()`:
+```python
+# In the main loop:
+if self._phase == "sd" and self.params.cg_enabled:
+    should_switch = False
+    if sd_iter_count >= self.params.sd_max_iter:
+        should_switch = True
+    if max_f < cg_switch_threshold:
+        should_switch = True
+    if should_switch:
+        self._switch_to_cg(iteration)
+```
+
+### CG Step Computation:
+```python
+def _cg_step(self, forces: np.ndarray) -> np.ndarray:
+    if self._cg_prev_forces is None:
+        # First CG step = SD direction
+        direction = forces.copy()
+    else:
+        # PRP+ beta
+        df = forces - self._cg_prev_forces
+        beta = np.dot(forces.ravel(), df.ravel()) / max(np.dot(self._cg_prev_forces.ravel(), self._cg_prev_forces.ravel()), 1e-20)
+        beta = max(beta, 0.0)  # PRP+ clamp
+
+        # Powell restart
+        if abs(np.dot(forces.ravel(), self._cg_prev_forces.ravel())) >= self.params.cg_restart_threshold * np.dot(forces.ravel(), forces.ravel()):
+            beta = 0.0
+
+        direction = forces + beta * self._cg_prev_direction
+
+    # Normalize and scale
+    max_disp = np.abs(direction).max()
+    if max_disp > 0:
+        step = self.params.max_step * direction / max_disp
+    else:
+        step = direction
+    step = self._clip_step(step)
+
+    # Store for next iteration
+    self._cg_prev_forces = forces.copy()
+    self._cg_prev_direction = direction.copy()
+
+    return step
+```
+
+### GDIIS Integration:
+- Same `_try_diis_acceleration()` method as current SD.py
+- Called in both SD and CG phases at `diis_store_every` intervals
+- On phase transition: `self.diis.reset()` + `self._sd_step_counter = 0`
+- Step validation identical: reject if energy rises AND force increases
+
+## Files Touched (Summary)
+| File | Action |
+|------|--------|
+| `algorithm/SDCG.py` | CREATE - main implementation |
+| `algorithm/SD.py` | MODIFY - backward-compat stub |
+| `algorithm/__init__.py` | MODIFY - add SDCG imports |
+| `optimization/optimization.py` | MODIFY - add routing |
+| `read/command_control.py` | MODIFY - add "sdcg" |
+| `example/opt/sdcg/inp1.inp` | CREATE - test input |
+| `example/opt/cg/inp1.inp` | CREATE - test input |
+
+## Verification
+1. Run SDCG default mode: `example/opt/sdcg/inp1.inp` - should show SD phase then CG phase transition
+2. Run SD-only mode: `example/opt/sd/inp1.inp` with `method=sd` - should behave identically to current SD
+3. Run CG-only mode: `example/opt/cg/inp1.inp` - should skip SD phase entirely
+4. Run with DIIS disabled: add `diis_enabled=false` - should show no GDIIS attempts
+5. Check .out files for phase transition logging, convergence, and trajectory output