CP2K 分子红外与拉曼光谱
H2O与CO2的真实结构、完整输入、有限位移、模式与IR/Raman活动数据。
阅读本篇目录
计算分子红外和拉曼,需要先把结构优化到驻点,再计算位移前后的力、偶极矩和极化率。只设置 INTENSITIES TRUE 不会自动替你打开这些输出。本篇用H₂O与CO₂分开说明非线性和线性分子。
下载完整输入、原始输出、频率表和绘图脚本,解压后进入 run.sh 所在目录。本文在鸢算科技P100上用CP2K2026.2双精度 cp2k.psmp 实算;构建采用GNU13.3、MPICH5.0.1、CUDA12.9、OpenBLAS0.3.33、ScaLAPACK2.2.3。分析环境为ASE3.29.0。
1 提交与查看结果#
sbatch submit.slurm water
模板申请四张P100、四个MPI进程,每进程一个CPU核和一个OpenMP线程。NPROC_REP 1让不同进程计算不同位移,脚本按MPI全局编号分配GPU。在已分配的计算资源里也可以执行 bash run.sh water。套餐用户用 yskit 生成对应分区的提交脚本,再保留这条运行命令;不要在登录节点直接运行。
脚本依次执行“优化结构 → 提取最终XYZ → 有限位移力、偶极与极化率 → 振动模式 → 数据提取 → 画图”。inputs为可编辑输入,reference为本文结果,work为本次新结果。每个条件请新解压一份,脚本不会覆盖已有work。
2 先优化,再算振动#
采用PBE、DZVP-MOLOPT-SR-GTH、GTH-PBE赝势、600/50 Ry网格和12 Å非周期MT盒。振动对数值噪声敏感,电子收敛EPS_SCF为1E−8;优化最大力1E−5 Hartree/Bohr,约0.000514 eV/Å,比一般结构优化更紧。
完整输入:water-opt#
&GLOBAL
PROJECT water
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../../DATA/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../../DATA/GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-12
&END QS
&MGRID
CUTOFF 600
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-8
MAX_SCF 100
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&END DFT
&SUBSYS
&CELL
ABC 12 12 12
PERIODIC NONE
&END CELL
&COORD
O 5.9999998739 6.0000000010 5.9992568887
H 5.9999999865 6.7657755697 6.5949200315
H 5.9999998417 5.2342244900 6.5949201085
&END COORD
&KIND O
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&KIND H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&END SUBSYS
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL
&MOTION
&GEO_OPT
OPTIMIZER BFGS
MAX_ITER 100
MAX_FORCE 1.0E-5
RMS_FORCE 5.0E-6
MAX_DR 1.0E-4
RMS_DR 5.0E-5
&END GEO_OPT
&END MOTION
脚本取 water-pos-1.xyz 的最后一帧写成optimized.xyz。振动输入使用该文件,而不是重新读取最初的未优化坐标。
完整输入:water-ir-raman#
&GLOBAL
PROJECT water
RUN_TYPE VIBRATIONAL_ANALYSIS
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../DATA/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../DATA/GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-12
&END QS
&MGRID
CUTOFF 600
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-8
MAX_SCF 100
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&END DFT
&PROPERTIES
&LINRES
EPS 1.0E-8
MAX_ITER 100
PRECONDITIONER FULL_SINGLE_INVERSE
&POLAR ON
DO_RAMAN .TRUE.
PERIODIC_DIPOLE_OPERATOR .FALSE.
&END POLAR
&END LINRES
&END PROPERTIES
&SUBSYS
&CELL
ABC 12 12 12
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME optimized.xyz
COORD_FILE_FORMAT XYZ
&END TOPOLOGY
&KIND O
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&KIND H
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q1
&END KIND
&END SUBSYS
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL
&VIBRATIONAL_ANALYSIS
NPROC_REP 1
DX 0.005
FULLY_PERIODIC .FALSE.
INTENSITIES .TRUE.
THERMOCHEMISTRY .TRUE.
TC_TEMPERATURE 298.15
TC_PRESSURE 101325
&PRINT
&MOLDEN_VIB ON
&END MOLDEN_VIB
&END PRINT
&END VIBRATIONAL_ANALYSIS
DX 0.005单位Bohr,约0.002646 Å。三个原子的中心有限位移共18个正负方向,分配给四个MPI进程;它不是18个独立材料体系。FULLY_PERIODIC FALSE配合分子边界处理刚体自由度。IR需要 DFT/PRINT/MOMENTS,Raman需要 PROPERTIES/LINRES/POLAR,非周期分子采用局域偶极算符。
3 查看原生频率和活性#
grep -E 'VIB\|Frequency|VIB\|IR int|VIB\|Raman' work/cases/water-ir-raman/output.log
cat work/analysis/water-ir-raman.csv

4 换成CO₂:保留线性几何,释放振动自由度#
新解压一份,执行 sbatch submit.slurm co2。CO₂优化时只允许三个原子的z坐标调整,保持x=y=6 Å的严格直线;在本例正交固定盒中,XY约束与笛卡尔方向一致。优化结束后,振动输入完全移除该约束,让所有原子三个方向都参与有限位移。若把约束一直带进振动,可能会漏掉弯曲模式。
完整输入:co2-opt#
&GLOBAL
PROJECT water
RUN_TYPE GEO_OPT
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../../DATA/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../../DATA/GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-12
&END QS
&MGRID
CUTOFF 600
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-8
MAX_SCF 100
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&END DFT
&SUBSYS
&CELL
ABC 12 12 12
PERIODIC NONE
&END CELL
&COORD
C 6 6 6
O 6 6 4.82
O 6 6 7.18
&END COORD
&KIND O
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&KIND C
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q4
&END KIND
&END SUBSYS
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL
&MOTION
&GEO_OPT
OPTIMIZER BFGS
MAX_ITER 100
MAX_FORCE 1.0E-5
RMS_FORCE 5.0E-6
MAX_DR 1.0E-4
RMS_DR 5.0E-5
&END GEO_OPT
&CONSTRAINT
&FIXED_ATOMS
LIST 1 2 3
COMPONENTS_TO_FIX XY
&END FIXED_ATOMS
&END CONSTRAINT
&END MOTION
完整输入:co2-ir-raman#
&GLOBAL
PROJECT water
RUN_TYPE VIBRATIONAL_ANALYSIS
PRINT_LEVEL MEDIUM
&END GLOBAL
&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME ../DATA/BASIS_MOLOPT
POTENTIAL_FILE_NAME ../DATA/GTH_POTENTIALS
&QS
EPS_DEFAULT 1.0E-12
&END QS
&MGRID
CUTOFF 600
REL_CUTOFF 50
&END MGRID
&SCF
EPS_SCF 1.0E-8
MAX_SCF 100
SCF_GUESS ATOMIC
&OT
MINIMIZER DIIS
PRECONDITIONER FULL_SINGLE_INVERSE
&END OT
&END SCF
&XC
&XC_FUNCTIONAL PBE
&END XC_FUNCTIONAL
&END XC
&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON
&PRINT
&MOMENTS ON
PERIODIC .FALSE.
REFERENCE COM
&END MOMENTS
&END PRINT
&END DFT
&PROPERTIES
&LINRES
EPS 1.0E-8
MAX_ITER 100
PRECONDITIONER FULL_SINGLE_INVERSE
&POLAR ON
DO_RAMAN .TRUE.
PERIODIC_DIPOLE_OPERATOR .FALSE.
&END POLAR
&END LINRES
&END PROPERTIES
&SUBSYS
&CELL
ABC 12 12 12
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME optimized.xyz
COORD_FILE_FORMAT XYZ
&END TOPOLOGY
&KIND O
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q6
&END KIND
&KIND C
BASIS_SET DZVP-MOLOPT-SR-GTH
POTENTIAL GTH-PBE-q4
&END KIND
&END SUBSYS
&PRINT
&FORCES ON
&END FORCES
&END PRINT
&END FORCE_EVAL
&VIBRATIONAL_ANALYSIS
NPROC_REP 1
DX 0.005
FULLY_PERIODIC .FALSE.
INTENSITIES .TRUE.
THERMOCHEMISTRY .TRUE.
TC_TEMPERATURE 298.15
TC_PRESSURE 101325
&PRINT
&MOLDEN_VIB ON
&END MOLDEN_VIB
&END PRINT
&END VIBRATIONAL_ANALYSIS
非线性H₂O应有3N−6=3个内部模式,线性CO₂应有3N−5=4个。CO₂两条弯曲模式几乎简并,不能因为频率重合就只保留一条。本文CO₂参考采用八个MPI/八张P100,公开脚本采用四个MPI/四张P100;改变的是位移的并行分配,仍计算相同18个位移。

5 本例得到什么#
| 分子/模式 | 频率/cm⁻¹ | IR强度/km·mol⁻¹ | Raman活性/Å⁴·amu⁻¹ |
|---|---|---|---|
| H₂O 模式1 | 1612.1373 | 69.6924 | 6.04528 |
| H₂O 模式2 | 3718.1052 | 19.546 | 48.6739 |
| H₂O 模式3 | 3822.7933 | 95.6427 | 14.0597 |
| CO₂ 模式1 | 670.1025 | 25.1708 | 1.95968e-10 |
| CO₂ 模式2 | 670.1033 | 25.1707 | 1.41208e-10 |
| CO₂ 模式3 | 1294.7412 | 1.60565e-05 | 15.0776 |
| CO₂ 模式4 | 2354.8074 | 515.215 | 3.39661e-07 |
CO₂对称伸缩在IR中很弱,反对称伸缩在非共振Raman中很弱;水的三个内部模式均为正频率。结果来自有限基组和网格下的谐近似,不乘经验缩放系数。Raman列是活性,不是加入激光频率、温度因子和仪器响应后的实验散射强度;两种图的纵轴分别归一化,不能互相比较峰高绝对值。
提取与重画#
module load ase/3.29.0
python analyze.py work
python plot_results.py work
ls work/analysis
ls work/figures
CSV给出模式编号、频率、IR强度和Raman活性;JSON保留质量、最终结构和热修正。原生 water-VIBRATIONS-1.mol 是Molden振动文件,可以用支持该格式的分子查看器检查模式位移。项目名前缀water是本包统一输出名,不代表CO₂输入成了水。
只重画本文结果,把上述命令中的work换成reference。脚本用实际模式的线位置与活性叠加高斯曲线,展示宽度σ=10 cm⁻¹;这是绘图展宽,不是计算得到的声子寿命。