GPAW激发态直接优化与稳定性分析

本篇用一个小水分子，练习激发态的轨道直接优化。先得到基态轨道，再把一个自旋通道的HOMO电子提升到LUMO，分别使用DO-MOM和DO-GMF优化电子态，检查电子Hessian的最低曲率方向，最后用真实原子位移核对解析力。

这里优化的是固定核位置下的电子轨道。没有进行水分子几何优化，也没有计算振动频率或实验吸收谱。

1 下载与运行顺序#

下载激发态直接优化算例，保存为gpaw-direct-excited.zip。登录节点解压后用sbatch job.slurm提交；已进入分配好的CPU资源时，依次执行：

python3 -m zipfile -e gpaw-direct-excited.zip .
cd gpaw-direct-excited
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

程序使用1个MPI进程、1个线程。计算顺序为：基态、DO-MOM激发态、DO-GMF激发态、H1沿x方向位移−0.003 Å、+0.003 Å的两组DO-MOM计算；两个中心激发态还分别进行Davidson电子稳定性分析。结果在results/，图片在plots/。

Hessian分析需要额外的能量和轨道梯度计算，通常比普通基态计算慢。同样输入再次运行会复用已完成结果；改变结构或方法时请在新目录计算。

2 水分子和基态初猜#

两个O—H键长均为0.9575 Å，H—O—H角为104.51°，四周加入4 Å真空，使用非周期边界。原子顺序为O、H1、H2，H1最初位于O的+x方向。

import numpy as np
from ase import Atoms
from gpaw import GPAW
from gpaw.directmin.etdm_lcao import LCAOETDM

angle = np.deg2rad(104.51)
bond = 0.9575
atoms = Atoms('OH2', positions=[
    [0, 0, 0], [bond, 0, 0],
    [bond * np.cos(angle), bond * np.sin(angle), 0]], pbc=False)
atoms.center(vacuum=4.0)

calc = GPAW(
    legacy_gpaw=True, mode='lcao', basis='dzp', h=0.22,
    xc='PBE', spinpol=True, nbands='nao', symmetry='off',
    eigensolver=LCAOETDM(representation='u-invar'),
    occupations={'name': 'fixed-uniform'},
    mixer={'backend': 'no-mixing'},
    convergence={'energy': 1e-6, 'density': 1e-6,
                 'eigenstates': 4e-9},
    maxiter=240, txt='ground.txt')
atoms.calc = calc
ground_energy = atoms.get_potential_energy()
calc.write('ground.gpw', mode='all')

nbands='nao'使LCAO表示保留与原子轨道数相同的轨道空间；no-mixing与直接轨道优化配合。这里明确使用自旋极化接口，以便只在一个自旋通道内提升电子。

3 DO-MOM：直接优化与最大重叠法配合#

从已收敛基态读取轨道，生成HOMO→LUMO的占据改变：

from gpaw.directmin.tools import excite
from gpaw.mom import prepare_mom_calculation

calc = GPAW('ground.gpw', legacy_gpaw=True, txt='DO-MOM.txt')
atoms = calc.get_atoms()
atoms.calc = calc
numbers = excite(calc, 0, 0, spin=(0, 0))
calc.set(eigensolver=LCAOETDM(
    searchdir_algo={'name': 'l-sr1p'},
    linesearch_algo={'name': 'max-step', 'max_step': 0.2},
    representation='u-invar', need_init_orbs=False))
prepare_mom_calculation(calc, atoms, numbers)
excited_energy = atoms.get_potential_energy()
forces = atoms.get_forces()
calc.write('DO-MOM.gpw', mode='all')

excite中的两个0分别表示从HOMO取走电子、放入LUMO；spin=(0,0)表示在同一个自旋通道内操作。它们不是两个原子的编号。

本例得到的是混合自旋开壳层行列式。上、下自旋电子数相同，只能说明净自旋投影为0，不能据此称它为已自旋纯化的单重态。

DO-MOM通过最大重叠法辅助保留所需电子态。换成电荷转移很强的体系时，仍需检查最终轨道和态的性质，不能只凭迭代结束判断找到了目标态。

4 DO-GMF：指定电子鞍点阶数#

同样从基态初猜出发，把求解器改成广义模跟随：

calc = GPAW('ground.gpw', legacy_gpaw=True, txt='DO-GMF.txt')
atoms = calc.get_atoms()
atoms.calc = calc
numbers = excite(calc, 0, 0, spin=(0, 0))
calc.set(
    eigensolver=LCAOETDM(
        searchdir_algo={'name': 'l-bfgs-p_gmf'},
        linesearch_algo={'name': 'max-step', 'max_step': 0.2},
        partial_diagonalizer={
            'name': 'Davidson', 'sp_order': 1, 'seed': 42,
            'logfile': 'GMF-inner-Davidson.txt'},
        representation='u-invar', update_ref_orbs_counter=1000,
        need_init_orbs=False),
    occupations={'name': 'mom', 'numbers': numbers,
                 'use_fixed_occupations': True})
gmf_energy = atoms.get_potential_energy()
gmf_forces = atoms.get_forces()
calc.write('DO-GMF.gpw', mode='all')

附件中的DO-GMF分支重新读取自己的基态检查点后执行这一步。sp_order=1表示本例目标是电子轨道旋转能量面上的一阶鞍点。这个阶数不能对所有分子、所有激发一律设为1。

use_fixed_occupations=True使本分支固定给定占据；DO-GMF通过跟随指定的低曲率模寻找目标鞍点。本例采用LCAO实现，不能直接把这里的GMF参数照搬到PW模式。

5 实际能量和方法对照#

两种方法在本模型的中心结构总能差约4.5×10⁻⁹ eV。这个接近程度说明它们在本例找到了相同的能量解，不等于8.272 eV就是可直接对应实验的吸收峰。

6 怎样检查电子稳定性#

对中心激发态调用原生广义Davidson方法：

from gpaw.directmin.derivatives import Davidson

davidson = Davidson(calc.wfs.eigensolver,
                    'stability.txt', eps=1e-3, seed=42)
order = davidson.estimate_sp_order(
    calc, method='full-hess', target_more=2)
print(order, davidson.lambda_e, davidson.error_e)

这里用Hessian作用的数值计算和部分特征分解检查最低几个电子曲率方向，没有构造并完整对角化整个Hessian矩阵。接口名full-hess用于区别初猜轨道能量的对角近似。

本例两种方法均识别出一个负曲率方向，最低三个曲率的实际数值如下，保留原生轨道旋转单位：

这些是电子轨道变化的曲率，不是原子位移Hessian，也不是cm⁻¹振动频率。附件还保存每个特征对的残差，不能只看曲率符号而忽略求解是否收敛。

对于新的激发态，应先估计所需鞍点阶数，并结合电子态性质选择。不要把本例的阶数、最低曲率或样本种子当成通用参数。

7 用真实位移检查解析力#

保持O、H2和晶胞不变，分别把H1的x坐标移−0.003 Å和+0.003 Å，从中心DO-MOM激发态检查点继续优化电子态：

calc = GPAW('DO-MOM.gpw', legacy_gpaw=True, txt='displaced.txt')
atoms = calc.get_atoms()
atoms.positions[1, 0] += 0.003  # 另一组改为-0.003
atoms.calc = calc
energy_displaced = atoms.get_potential_energy()

本例直接沿用检查点保存的求解器与占据，不额外重置它们。真实计算得到：

差分使用F_x = −[E(+δ) − E(−δ)]/(2δ)。本例解析力很大，说明当前激发态几何远未达到受力很小的结构；不能因为解析力与差分吻合，就说结构优化已经完成。

8 密度与轨道数据#

前两图为水分子平面上的实际全电子总密度，第三图为激发态上、下自旋密度之差。净自旋投影为0时，局部自旋密度仍可以不为0。密度色标下限以下统一显示为最低颜色，原始数组保留。

每组fields.npz还保存完整LCAO系数、PAW重叠矩阵、占据和C* S Cᵀ度量。不同运行的轨道相位或等价轨道旋转可能不同，不能仅凭系数逐项不相同就判断计算错误；应结合密度、能量、力及PAW度量下的正交和电子数检查。

9 导出能量和电子曲率#

source env.sh
"$YS_PY" - <<'PYCODE'
from pathlib import Path
import json, csv, numpy as np
r = Path('results')
v = json.loads((r / 'VALIDATION.json').read_text())
with open('direct-excited.csv', 'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(['condition', 'energy_eV', 'H1_Fx_eV_A'])
    for row in v['results']:
        w.writerow([row['name'], row['energy_eV'], row['forces_eV_A'][1][0]])
for name in ['DO-MOM', 'DO-GMF']:
    d = np.load(r / name / 'stability.npz')
    table = np.column_stack([d['Hessian_eigenvalues_native'],
                             d['Hessian_residual_max_components']])
    np.savetxt(name + '-curvature.csv', table, delimiter=',',
               header='electronic_curvature_native,residual_max', comments='')
    print(name, table)
print(v['summary'])
PYCODE

只重画附件参考数据：

bash setup.sh
bash plot.sh --reference
ls reference-plots

自己的数据和图片分别在results/、plots/，附件参考数据与参考图片在reference/、reference-plots/。下载包不包含PAW数据集、GPW检查点或程序二进制。

10 环境与后续计算#

本例使用GPAW 26.7.0 legacy LCAO、Python 3.11.16，CPU构建采用GCC 13.3、HPC-X Open MPI 5.0.10rc2、MKL 2023.2和Libxc 7.0.0。setup.sh加载gpaw/26.7.0-8173m。

更换体系时，先确定目标电子态和自旋处理，检查基组、真空和网格，再选择初猜与电子优化算法。需要自旋纯化、复杂电荷转移态或激发态几何优化时，应进一步核对相应方法条件，不能只复制本例的一个占据改变。

接口参考：GPAW激发态直接优化说明、GPAW DO-GMF与稳定性分析说明。