准经典分子动力学
准经典分子动力学也称为准经典轨迹（QCT），与经典动力学相比，其考虑了零点能（ZPE），在化学反应的研究中被广泛使用（例如，参见Houk等人的文章 PNAS 2012, 109, 12860 ， J. Am. Chem. Soc. 2023, 145, 14446 ）。在本教程中，我们展示了如何执行此类模拟的 示例 ，示例内容为重现上述的PNAS论文。另请参阅有关 分子动力学 的相关教程。
以下文章提供了实现细节（在使用此功能时，请与其他必需的引用一起引用）：
Yi-Fan Hou, Lina Zhang, Quanhao Zhang, Fuchun Ge, Pavlo O. Dral. Physics-informed active learning for accelerating quantum chemical simulations. arXiv 2024, DOI: 10.48550/arXiv.2404.11811.
背景
进行准经典分子动力学模拟的关键是确保初始条件（几何结构和速度）包含所需的零点能（ZPE）。MLatom支持两种简正模采样方法来生成这些条件：
Wigner采样（参见 对应论文 中的公式20）仅从最低振动能级进行采样。此方法在 面跳跃动力学 和光谱生成中也被广泛使用。
sampling from the harmonic quantum Boltzmann distribution (implemented following this work, also called “thermal sampling” in the manual to its early implementation in the VENUS program). In contrast to Wigner sampling, it samples from higher vibrational levels. This sampling is commonly used for performing QCT investigation of transition states: for all the real normal modes that are perpendicular to the reaction coordinate, the quantum number \(n_{i}\) of normal mode \(i\) is first sampled from a harmonic quantum Boltzmann distribution:
\(p(n_{i})=e^{- n_{i} h v_{i} / k_{\text{B}} T} (1 - e^{-h v_{i} / k_{\text{B}} T})\)
其中 \(h\) 是普朗克常数， \(v_{i}\) 是振动频率， \(k_{B}\) 是玻尔兹曼常数， \(T\) 是温度。坐标 \(Q_{i}\) 和动量 \(P_{i}\) 由以下公式生成：
\(Q_{i}=A_{i} \text{cos} (2 \pi R_{i})\),
\(P_{i}=-\omega_{i}A_{i} \text{cos} (2 \pi R_{i})\),
此处， \(\omega_{i}=2 \pi v_{i}\) ， \(A_{i}= \sqrt{2E_{i}} / \omega_{i}\) 和 \(R_{i}\) 是 \([0,1]\) 内的均匀随机数。反应坐标（对应虚频的简正模）固定，质量加权动量 \(P\) 的计算公式为 \(\pm \sqrt{-2 k_{B} T \text{ln} (1-R)}\) ，其中 \(R\) 也是 \([0,1]\) 内的一个均匀随机数。
示例
下面我们将向您展示如何在MLatom中运行准经典分子动力学模拟，并重现 PNAS 2012, 109, 12860 和 arXiv.2404.11811 中1,3-丁二烯和乙烯的Diels-Alder反应的结果。
一般步骤如下：
优化构型，此例中我们进行过渡态优化
进行 频率计算 ，生成以JSON格式转储的分子对象文件（此处我们提供了该文件
da_eqmol.json，是过渡态的）使用分子文件中的简正模（此例中是
da_eqmol.json）来生成初始条件使用这些条件进行分子动力学模拟
分析结果。此处，我们将计算穿越过渡区的平均时间和C-C键形成的平均时间间隔。
请参阅以下完整版本，其中包含Python API，以及一个简要版本的输入文件（仅显示如何进行简正模采样和传播轨迹）。
输入文件示例
生成初始条件并运行单个轨迹可以使用以下输入文件示例：
MD # 1. requests molecular dynamics
method=UB3LYP/6-31G* qmprog=gaussian # 2. use DFT method (if you have Gaussian). You can also use, e.g., AIQM1 or ANI-1ccx, or B3LYP and qmprog=PySCF...
initConditions=harmonic-quantum-boltzmann # 3. use harmonic quantum Boltzmann distribution
#initConditions=wigner # 3. also possible - use Wigner sampling
normalModeFile=da_eqmol.json # 4. molecule with frequencies and normal modes
initTemperature=300 # 5. initial temperature
trun=150 # 6. total time; Unit: fs
dt=0.5 # 7. time step; Unit: fs
您将获得多个与此轨迹相关的输出文件（例如，可使用VMD查看的xyz坐标文件，以及h5md轨迹文件等）。
上面的示例将花费很长的时间（一个轨迹30分钟），但如果使用机器学习，速度会快得多。这里提供了一个针对该反应预训练的模型 da_energy_ani.npz ，使用以下输入文件进行相同的模拟只需不到一分钟：
MD
MLmodelType=ANI
MLmodelIn=da_energy_ani.npz
initConditions=harmonic-quantum-boltzmann
normalModeFile=da_eqmol.json
initTemperature=300
trun=150
dt=0.5
Python API的示例
一个使用Python API的更完整的示例：
import mlatom as ml
import os
import numpy as np
# Load molecule with frequencies and normal modes
eqmol = ml.data.molecule()
eqmol.load('da_eqmol.json',format='json')
# Generate initial conditions
if not os.path.exists('init_cond_db.json'):
init_cond_db = ml.generate_initial_conditions(molecule=eqmol,
generation_method='harmonic-quantum-boltzmann',
initial_temperature=298,
number_of_initial_conditions=128)
init_cond_db.dump('init_cond_db.json',format='json')
else:
init_cond_db = ml.data.molecular_database.load('init_cond_db.json',format='json')
# Reverse the velocities to get the backward trajectories
re_init_cond_db = init_cond_db.copy()
for mol in re_init_cond_db:
for atom in mol:
atom.xyz_velocities = -atom.xyz_velocities
# Define method
method = ml.models.methods(method='UB3LYP/6-31G*',program='Gaussian',nthreads=8)
# Run dynamics
for imol in range(128):
if not os.path.exists(f'forward_traj{imol+1}.xyz'):
dyn = ml.md(model=method,
molecule_with_initial_conditions = init_cond_db[imol],
ensemble='NVE',
time_step=0.5,
maximum_propagation_time=150.0)
traj = dyn.molecular_trajectory
traj.dump(filename=f'forward_traj{imol+1}',format='plain_text')
if not os.path.exists(f'backward_traj{imol+1}.xyz'):
dyn = ml.md(model=method,
molecule_with_initial_conditions = re_init_cond_db[imol],
ensemble='NVE',
time_step=0.5,
maximum_propagation_time=150.0)
traj = dyn.molecular_trajectory
traj.dump(filename=f'backward_traj{imol+1}',format='plain_text')
# Analyze trajectories
def bond_length(molecule):
d1 = molecule.internuclear_distance(0,10)
d2 = molecule.internuclear_distance(6,11)
return d1, d2
def check_traj(forward_traj,backward_traj):
forward_label = 'ts'
backward_label = 'ts'
forward_time = -1
backward_time = -1
forward_distances_list = []
backward_distances_list = []
for imol in range(len(forward_traj)):
d1,d2 = bond_length(forward_traj[imol])
min_dist = min(d1,d2)
max_dist = max(d1,d2)
distances = [d1,d2]
forward_distances_list.append(distances)
if min_dist > 5.0:
forward_label = 'reactant'
break
if max_dist < 1.6:
forward_label = 'product'
break
if (min_dist > 2.52 or max_dist < 2.02) and forward_time < 0:
forward_time = imol*0.5
for imol in range(len(backward_traj)):
d1,d2 = bond_length(backward_traj[imol])
min_dist = min(d1,d2)
max_dist = max(d1,d2)
distances = [d1,d2]
backward_distances_list.append(distances)
if min_dist > 5.0:
backward_label = 'reactant'
break
if max_dist < 1.6:
backward_label = 'product'
break
if (min_dist > 2.52 or max_dist < 2.02) and backward_time < 0:
backward_time = imol*0.5
# Check time gap
if (forward_label=='product' and backward_label=='reactant') or (forward_label=='reactant' and backward_label=='product'):
if forward_label == 'product':
distances_list = forward_distances_list
else:
distances_list = backward_distances_list
bond1_formed = False
bond2_formed = False
for istep in range(len(distances_list)):
distances = distances_list[istep]
# print(distances)
if min(distances) < 1.6 and not bond1_formed:
gap = istep * 0.5
bond1_formed = True
if max(distances) < 1.6 and not bond2_formed:
gap = istep * 0.5 - gap
bond2_formed = True
break
else:
gap = None
return forward_label,backward_label,forward_time+backward_time,gap
reactive = 0
nonreactive = 0
traverse_time_list = []
gap_list = []
for imol in range(128):
forward_traj = ml.data.molecular_database.from_xyz_file(f'forward_traj{imol+1}.xyz')
backward_traj = ml.data.molecular_database.from_xyz_file(f'backward_traj{imol+1}.xyz')
forward_label, backward_label, traverse_time, gap = check_traj(forward_traj,backward_traj)
if (forward_label == 'product' and backward_label == 'reactant') or (forward_label == 'reactant' and backward_label == 'product'):
reactive += 1
# print(f'{forward_time} + {backward_time} = {forward_time+backward_time}')
traverse_time_list.append(traverse_time)
gap_list.append(gap)
else:
nonreactive += 1
print(f"Number of reactive trajectories: {reactive}")
print(f"Number of non-reactive trajectories: {nonreactive}")
print(f"Average time: {np.mean(traverse_time_list)} fs")
print(f"Median time: {np.median(traverse_time_list)} fs")
print(f"Standard deviation: {np.std(traverse_time_list)} fs")
print(f"Average time gap: {np.mean(gap_list)} fs")
print(f"Median time gap: {np.median(gap_list)} fs")
print(f"Standard deviation: {np.std(gap_list)} fs")
然而，每个DFT轨迹大约需要30分钟，128个轨迹将需要超过2天的时间。使用机器学习模型可以大大加速模拟过程，因此我们提供了这个反应的预训练模型 da_energy_ani.npz 。只需将DFT方法替换为这个ANI模型即可：
method = ml.models.ani(model_file='da_energy_ani.npz',device='cpu') # If you want to use gpu, set device='cuda'
此处提供了我们的初始条件（ init_cond_db.json ），以便您复现我们的结果。模拟大约需要15分钟（如果您使用MLatom的 simulation.run_in_parallel 函数并行运行，速度会快得多），您将能够获得以下输出：
Number of reactive trajectories: 121
Number of non-reactive trajectories: 7
Average time: 53.35537190082645 fs
Median time: 53.0 fs
Standard deviation: 11.622406049624619 fs
Average time gap: 3.479338842975207 fs
Median time gap: 3.0 fs
Standard deviation: 2.885126215935735 fs