Critic2专题

Critic2 位移求力与力常数实操

从生成Si16位移结构、实际QE单点求力，到原生力常数、声子和振动热力学的完整终端流程。

更新于 2026-10-30约 7 分钟阅读按任务与环境选择适用配置

阅读本篇目录

没有 FORCE_CONSTANTS 时，先生成位移结构，用电子结构软件计算它们的力，再把真实输出交给Critic2。本例用固定晶格Si16、PBE/ONCV赝势和QE单点实际完成这一流程，不需要先准备一份外部力常数。

1. 获取输入并选运行方式#

unzip critic2-dft-forces-wave146-*.zip -d critic2-dft-forces
cd critic2-dft-forces
bash run.sh list

Critic2模块为 critic2/1.4（内部1.4.155，GCC/GFortran13.3、OpenMP、LAPACK/BLAS）；QE模块为 qe/7.6-8173m（InteloneAPI2024.2、MPI+OpenMP、MKL）。对照与绘图采用 phono3py/4.8.2 独立Python3.11.16环境内的Phonopy4.8.1与Symfc1.7.3。完整DFT阶段需要已分配至少4个CPU：单MPI进程×4OpenMP线程；Critic2后处理仅使用其中1核、1线程。普通计算按账户资源使用Slurm，包内提供 run.slurm 模板。

从DFT开始完整重算，在新解压目录运行：

bash run-all.sh

它依次执行：生成有限位移和随机位移结构→写入明确的QE电子计算参数→实际执行6个SCF→将新计算的力和新位移记录传递给Critic2→构造FC2→计算频率及热力学→对照和画图。其中2个SCF用于本篇的参考结构和单个位移结构，另4个用于随机配对与回归。

若只想学习有限位移后处理，在另一份新解压目录运行：

bash run.sh 2
head -n 6 work/fd-read/FORCE_CONSTANTS

此命令读取包内附带的真实QE输出，不重新执行DFT。已有结果目录不会被覆盖，重复尝试请使用新解压目录。

2. 生成实际位移结构#

对应的独立原生命令由 bash run.sh 1执行；完整输入：

crystal ../../data/unit.vasp
units bohr
vibrations calculator qe
vibrations memory 64 MB
vibrations create_displacements 2 2 2 distance .02 template si-*.scf.in dataset critic2_disp.txt

原胞为2个Si，2 2 2生成16原子超胞。UNITS BOHR使 DISTANCE .02明确表示 0.02bohr≈0.010584Å。利用当前结构的对称性，本例只需1个独立位移；si-000.scf.in为未位移参考结构，si-001.scf.in为位移结构。对其他结构，独立位移数量由对称性决定，不能预先固定为1。

critic2_disp.txt记录晶格、原子、超胞矩阵、位移长度和文件映射。求力与后处理期间一起保留它，不能只留下某个输出文件。

实际Si位移与能量和力响应的比较

3. 配置并运行真实QE计算#

生成器写出的QE文件是通用结构模板。完整脚本中的 prepare_qe_runtime.py保留实际几何和原子顺序，并明确写入本例PBE赝势、截止能、电子k网格与收敛阈值；不会直接执行通用模板里的XDM默认设置。

位移结构实际使用的完整QE输入如下：

&CONTROL
 calculation='scf', prefix='si', pseudo_dir='../../pseudo', outdir='./tmp', tprnfor=.true., tstress=.true., verbosity='high'
/
&SYSTEM
 ibrav=0, nat=16, ntyp=1, ecutwfc=40, ecutrho=160, nbnd=40, nosym=.true., noinv=.true.
/
&ELECTRONS
 conv_thr=1d-10, mixing_beta=0.5, electron_maxstep=150
/
ATOMIC_SPECIES
Si 28.0855 Si_ONCV_PBE-1.0.upf
ATOMIC_POSITIONS crystal
Si    0.00137821     0.00000000     0.00000000 
Si    0.50000000     0.00000000     0.00000000 
Si    0.00000000     0.50000000     0.00000000 
Si    0.50000000     0.50000000     0.00000000 
Si    0.00000000     0.00000000     0.50000000 
Si    0.50000000     0.00000000     0.50000000 
Si    0.00000000     0.50000000     0.50000000 
Si    0.50000000     0.50000000     0.50000000 
Si    0.12500000     0.12500000     0.12500000 
Si    0.62500000     0.12500000     0.12500000 
Si    0.12500000     0.62500000     0.12500000 
Si    0.62500000     0.62500000     0.12500000 
Si    0.12500000     0.12500000     0.62500000 
Si    0.62500000     0.12500000     0.62500000 
Si    0.12500000     0.62500000     0.62500000 
Si    0.62500000     0.62500000     0.62500000 

CELL_PARAMETERS bohr
    0.000000000000    10.261212901569    10.261212901569 
   10.261212901569     0.000000000000    10.261212901569 
   10.261212901569    10.261212901569     0.000000000000 

K_POINTS automatic
2 2 2 0 0 0

TPRNFOR输出原子受力；CONV_THR=1d-10的单位是Ry，较紧的电子收敛用于小位移求力。NOSYM与NOINV保持位移计算的显式电子k采样；截止能40/160Ry和2³电子k网格是本例固定设置，正式研究应自行检查收敛。这里没有结构优化，不能把该晶格当作已经针对这套计算参数优化完成。

完整重算后检查：

grep 'convergence has been achieved' qe-work/si-001/out.out
grep -A18 'Forces acting on atoms' qe-work/si-001/out.out

附带的源结果在 reference/qe-work/，读取示例：

grep 'convergence has been achieved' reference/qe-work/si-001/out.out
grep -A18 'Forces acting on atoms' reference/qe-work/si-001/out.out

QE详细输出中还有non-local等受力分量；应使用 Forces acting on atoms后紧跟的16行总力，不能把后面的某一分量或多组分量一起当作总力。

4. 用原生READ_FORCES构造FC2#

bash run.sh 2的完整输入如下，完整重算则由脚本自动传递新的dataset与新力：

crystal ../../data/unit.vasp
units bohr
vibrations calculator qe
vibrations memory 64 MB
vibrations read_forces dataset critic2_disp.txt
vibrations write_fc2 FORCE_CONSTANTS full
vibrations print summary
vibrations calcq file points.dat
q 0 0 0
q .5 .5 0
q .5 .5 .5
q .13 .27 .31
ENDCALCQ
vibrations thermo mesh 6 6 6 shift .5 .5 .5 temperatures 0 300 1000 noxdebye file thermo.dat

READ_FORCES从dataset推断 si-*.scf.out，核对输出结构与预期位移，然后构造力常数。它会按x、y、z分别扣除每份总力的原子平均值，即整体受力漂移；本例finite reader 不会自动减去si-000的参考力。参考结构用来检查未位移残余力，两个步骤不要混淆。

WRITE_FC2导出的单位跟随 CALCULATOR QE，本例为 Ry/bohr²。传入通常采用eV/Å²的Phonopy数据流程时，需要乘 13.605693/0.52917721²换算；不能因文件名同叫FORCE_CONSTANTS就按VASP的单位读取。

head -n 6 work/fd-read/FORCE_CONSTANTS
cat work/fd-read/points.dat
cat work/fd-read/thermo.dat

实际Γ点三个光学模均为 522.346462cm⁻¹；Γ点声学模约10⁻⁵cm⁻¹。6³移位q网格300K得到Fvib=6.0731717kJ/mol原胞、S=41.5678402J/K/mol原胞、Cv=39.5566105J/K/mol原胞。Fvib只是振动贡献，原胞含2Si，不能直接把它当作每1Si化学式单元的值。

5. 读取能量与作图#

源数据的位移能量增量为 0.775797meV，由实际力常数得到的 0.5 u^T Phi u为 0.775343meV。两者来自不同的实际输出关系，不强制写成相等。

只处理附带数据并画图：

bash analyze.sh
ls figures-reference

同一DFT设置下有限位移与配对回归的离散q点频率

这张图是X、L两个离散q点的6模频谱，横轴为模式编号，并非连续声子路径。连续路径、DOS与振动热力学的操作见力常数与声子色散和声子态密度与振动热力学。本例用于学习流程，未证明截止能、k/q网格、超胞和位移幅度的材料物性收敛。QE参数含义可参看官方PW输入说明。

完整重算的数值与图位于 analysis/CHECKS.json 和 figures/；只分析附带真实数据时，对应为 analysis-reference/CHECKS.json 与 figures-reference/。

需要根据你的体系选择资源？联系讨论计算需求。实际可用软件与权限，以当前账户和运行环境为准。