【翻译】IRC 的“最后手段”——沿虚频方向移动结构！
上篇文章介绍了更换算法或使用 GRRM 等常规解决方案，但即便如此，也可能仍无法顺利完成 IRC 计算。
这时，可以尝试一招“最后手段”：手动沿虚频方向移动结构，然后进行结构优化。
只要你能理解 log 文件中的内容，沿虚振动方向移动结构其实并不难；但对于不熟悉的人来说，可能还是有些挑战。
因此，本文将介绍一种通过手动改变结构以解决 IRC 问题的方法。
进行 TS 的振动计算
首先，需要准备好进行过渡态（TS）振动计算的文件。通常，在执行 opt=ts 时会一并加入 freq 关键字。
这次我们以以下论文中的 TS3 为例来讲解：
“Enzyme-catalysed [6+4] cycloadditions in the biosynthesis of natural products”
Zhang, B. et al. Nature 2019, 568, 122–138.
如何读取 log 文件中的振动模式
在执行完频率（freq）计算的 log 文件中，可以找到如下片段：
Frequencies -- ...用关键词“Frequencies —”搜索即可快速定位。
每一个振动模式下方，都会列出每个原子在 x、y、z 三个方向上的位移量。与 GaussView 中的显示内容对应查看会更直观。
图示参考如下：
频率段截图
位移量截图只需将这些位移值乘以一个系数（例如 -1.0 到 1.0）并加到原始坐标中，就可以得到沿虚频方向偏移后的新结构。
接下来我们说明如何保存这些偏移后的结构。
方法一：使用 GaussView
如果您拥有 GaussView，可依次选择：
Results > Vibrations
→ 勾选 Manual Displacement
→ 移动滑动条来调整位移幅度（可正向或反向）
当结构偏移到合适位置时，点击旁边的 “Save Structure” 即可保存该结构。
方法二：使用程序脚本
若您没有 GaussView，仅有 Gaussian，也可以通过编写简单脚本来完成坐标偏移。
以下为一个用 python 编写的示例程序（原博主未提供，我让AI写了一个）：
import re
import numpy as np
import os
def parse_gaussian_output(filename):
with open(filename, 'r') as f:
lines = f.readlines()
# 提取原始坐标
coords = []
atoms = []
for i, line in enumerate(lines):
if 'Standard orientation:' in line:
coords = []
atoms = []
j = i + 5
while '-----' not in lines[j]:
parts = lines[j].split()
atoms.append(parts[1])
coords.append([float(parts[3]), float(parts[4]), float(parts[5])])
j += 1
# 提取虚频位移向量
displacements = []
for i, line in enumerate(lines):
if 'Frequencies --' in line:
freq = float(line.split()[2])
if freq < 0:
disp = []
for j in range(len(atoms)):
disp_line = lines[i + 5 + j]
parts = disp_line.split()
disp.append([float(parts[2]), float(parts[3]), float(parts[4])])
displacements.append(disp)
break # 只处理第一个虚频
return atoms, np.array(coords), np.array(displacements[0])
def write_gjf(atoms, coords, filename, charge=0, multiplicity=1, method='B3LYP', basis='6-31G(d)', job_type='opt freq'):
with open(filename, 'w') as f:
f.write(f"%chk={os.path.splitext(filename)[0]}.chk\n")
f.write(f"# {job_type} {method}/{basis}\n\n")
f.write(f"Generated by script\n\n")
f.write(f"{charge} {multiplicity}\n")
for atom, coord in zip(atoms, coords):
f.write(f"{atom} {coord[0]:.6f} {coord[1]:.6f} {coord[2]:.6f}\n")
f.write('\n')
def generate_structures(atoms, coords, displacement, step=0.1, output_dir='structures'):
if not os.path.exists(output_dir):
os.makedirs(output_dir)
factors = np.arange(-1.0, 1.0 + step, step)
for factor in factors:
new_coords = coords + factor * displacement
base_filename = os.path.join(output_dir, f'structure_{factor:.1f}')
gjf_filename = f"{base_filename}.gjf"
write_gjf(atoms, new_coords, gjf_filename)
# 使用示例
if __name__ == '__main__':
gaussian_output_file = 'your_gaussian_output.out' # 替换为你的 Gaussian 输出文件名
atoms, coords, displacement = parse_gaussian_output(gaussian_output_file)
generate_structures(atoms, coords, displacement)程序会将虚振动的位移向量按 0.1 为步长，从 -1.0 倍至 1.0 倍依次加到坐标上。输出到structures文件夹下，格式为.gjf。
使用方法如下：
- 安装python环境
pip install numpypython test.py
与实际 IRC 计算结果对比
笔者亲测发现，使用虚频方向偏移结构后再进行结构优化，与通过 IRC 路径计算再优化所得结构在大多数情况下没有显著差异。
其中的技巧在于：加到坐标中的位移量控制在约 0.1 倍左右。
若使用 1.0 倍或更大的位移，很容易导致结构严重扭曲，反而增加优化难度与耗时。
本质上，我们只是轻轻把结构从过渡态“山峰”上推一小步，接下来它就会自然地“滑”向谷底（即产物或反应物）。
总结
在论文的实验部分或“Computational Details”中，经常可以见到如下语句：
The intrinsic reaction coordinate (IRC) calculations were performed at the same level of theory to assure the transition states could connect the reactants and expected product.
本征反应坐标（IRC）计算在同一理论水平上进行，以确保过渡态能够连接反应物和预期产物。
—— J Phys Org Chem. 2019, e4035.
Intrinsic reaction coordinate (IRC) calculations were used to characterize transition state structures.
本征反应坐标（IRC）计算用于确定过渡态结构的特征。
—— Tetrahedron Lett. 2012, 53, 6919–6922.
由此可见，在完成过渡态结构优化后，进行 IRC 计算是不可或缺的环节。
但如果尝试各种方法仍无法成功，也可以考虑手动沿虚频方向移动结构后再优化作为备选策略——前提是必须在文中明确说明 IRC 计算未成功这一事实。
关于如何在论文中写出“IRC 计算失败的正当理由”，可参考上篇文章提供的例子。