分子动力学
本教程展示了如何使用MLatom进行分子动力学模拟
首先我们演示如何使用输入文件/命令行进行分子动力学模拟，此外还有更加通用的使用 Python API 的教程。您可以下载分子的初始XYZ坐标文件 h2_md_kreg_init.xyz，初始XYZ速度文件 h2_md_kreg_init.vxyz 和模型文件 h2_kreg_energies.unf。
输入文件 h2_md_kreg.inp 展示了如何使用Nosé-Hoover恒温器模拟NVT系综下氢分子的分子动力学：
# h2_md_kreg.inp
MD # 1. requests molecular dynamics
initConditions=user-defined # 2. use user-defined initial conditions
initXYZ=h2_md_kreg_init.xyz # 3. file with initial geometry; Unit: Angstrom
initVXYZ=h2_md_kreg_init.vxyz # 4. file with initial velocity; Unit: Angstrom/fs
dt=0.3 # 5. time step; Unit: fs
trun=30 # 6. total time; Unit: fs
thermostat=Nose-Hoover # 7. use Nose-Hoover thermostat
ensemble=NVT # 8. NVT ensemble
temperature=300 # 9. Run MD at 300 Kelvin
MLmodelType=KREG # 10. KREG model is used
MLprog=MLatomF # 11. use KREG implemented in the Fortran part of MLatom
MLmodelIn=h2_kreg_energies.unf # 12. file with the trained model
初始坐标以XYZ格式(Å)定义:
2
H 0.0 0.0 1.0
H 0.0 0.0 0.0
初始速度也以简化的XYZ格式定义(Å/fs):
2
0.0 0.0 -0.05
0.0 0.0 0.05
运行分子动力学模拟。
mlatom h2_md_kreg.inp &> h2_md_kreg.out
模拟结束后，分子的XYZ坐标轨迹将被保存在 traj.xyz 文件中，XYZ速度将被保存在 traj.vxyz 文件中动能、势能和总能（单位Hartree）将被分别保存在 traj.ekin, traj.epot 和 traj.etot 中,能量梯度（单位Hartree/Å）将被保存在 traj.grad 中。
备注
在MD中，通过输入文件或命令行，平移动量和角动量默认被移除。用户不能像Python API那样保持角动量或改变自由度。为确保正确的自由度，MLatom会自动检测分子的线性。
初始条件
除了使用上述方式自行定义初始条件，用户还可以通过输入文件生成随机速度。 initTemperature 用于控制分子的动能( \(E_\text{kin}=N_\text{DOF}k_\text{B}T_\text{init}/2\) )。
MD # 1. requests molecular dynamics
initConditions=random # 2. generate random velocities
initXYZ=h2_md_kreg_init.xyz # 3. file with initial geometry; Unit: Angstrom
initTemperature=300 # 4. initial temperature
dt=0.3 # 5. time step; Unit: fs
trun=30 # 6. total time; Unit: fs
thermostat=Nose-Hoover # 7. use Nose-Hoover thermostat
ensemble=NVT # 8. NVT ensemble
temperature=300 # 9. Run MD at 300 Kelvin
MLmodelType=KREG # 10. KREG model is used
MLprog=MLatomF # 11. use KREG implemented in the Fortran part of MLatom
MLmodelIn=h2_kreg_energies.unf # 12. file with the trained model
用户也可以使用简正模采样来生成初始条件。我们提供了两个选项：Wigner分布(initConditions=wigner)和量子谐波Boltzmann分布(initConditions=quantum-harmonic-boltzmann)。
MD # 1. requests molecular dynamics
method=AIQM1 # 2. use AIQM1
initConditions=wigner # 3. use Wigner distribution
#initConditions=harmonic-quantum-boltzmann # 3. use harmonic quantum Boltzmann distribution
normalModeFile=ethanol_freq.json # 4. molecule with frequencies and normal modes
initTemperature=300 # 5. initial temperature
trun=30 # 6. total time; Unit: fs
dt=0.3 # 7. time step; Unit: fs
normalModeFile 指定保存分子对象的json文件，该文件应该包含简正模采样所需的信息，即频率和简正模。在做 频率计算 时可自动生成该文件（检查 freq[你的分子数].json 文件）。
选择系综和恒温器
MLatom目前支持NVE和NVT系综，可以分别使用 ensemble=NVE 和 ensemble=NVT 关键字定义。在NVT系综下，用户可以选择Andersen恒温器( thermostat=Andersen )和Nosé-Hoover恒温器( thermostat=Nose-Hoover )。
选择方法或模型
MLatom为用户提供了多种方法或模型用于MD。
上述例子是由用户训练的KREG模型，并使用以下关键字定义:
...
MLmodelType=KREG # 10. KREG model is used
MLprog=MLatomF # 11. use KREG implemented in the Fortran part of MLatom
MLmodelIn=h2_kreg_energies.unf # 12. file with the trained model
关于如何训练此类模型，请参阅具体的手册和教程，例如 MACE。
以上几行可以替换为一些可被自动识别的方法，如ANI-1ccx或AIQM1，举个简单的例子:
...
ANI-1ccx # or AIQM1
定义QM方法的一般方法是使用 method 和对应的QM程序(MLatom使用这些程序的接口)。例如，可以选择从头算方法(如HF或MP2)或DFT方法(如B3LYP/6-31G*):
...
method=B3LYP/6-31G* # 2. request running DFT optimization with B3LYP
qmprog=PySCF # 3. request using PySCF for B3LYP calculations; qmprog=Gaussian can be also used
备注
原则上，我们建议使用AIQM1，因为该方法通常是一定成本下最准确的方法。但目前仅适用于含有CHNO元素的化合物。另外，为了使其更好地得到运用，建议为其QM部分安装MNDO程序。XACS云平台使用的是Sparrow，它目前只提供数值梯度，限制了它在分子动力学中的使用。因此，改用ANI-1ccx是XACS云平台的一个很好的替代方案(它对基态中性闭壳层化合物有额外的限制)。
使用PyAPI
在PyAPI中运行分子动力学包括4个步骤： 生成初始条件, 初始化模型, 初始化恒温器 和 运行分子动力学。此处我们使用预训练的KREG模型来进行 \(\text{H}_2\) 的MD模拟。
初始条件
MLatom提供了一个函数 ml.generate_initial_conditions 来生成MD的初始几何构型和速度，该函数返回一个 ml.data.molecular_database 类，类中的 ml.data.molecule 对象包含这些属性。下面的示例使用用户提供的初始条件。它们可以是先前计算的初始条件，也可以是其他软件生成的初始条件。
import mlatom as ml
# Use user-defined initial conditions
mol = ml.data.molecule.from_xyz_file('h2_md_kreg_init.xyz')
init_cond_db = ml.generate_initial_conditions(molecule=mol,
generation_method='user-defined',
file_with_initial_xyz_coordinates='h2_md_kreg_init.xyz',
file_with_initial_xyz_velocities='h2_md_kreg_init.vxyz')
init_mol = init_cond_db[0]
用户还可以通过 generation_method='random' 生成随机速度。 initial_temperature (单位：开尔文)可以用来控制分子的动能。动能由 \(N_\text{DOF}k_\text{B}T_\text{init}/2\) 计算，其中 \(N_\text{DOF}\) 是分子的自由度， \(k_\text{B}\) 是玻尔兹曼常数，而 \(T_\text{init}\) 是用户设置的初始温度。用户也可以使用 initial_kinetic_energy 参数直接设置初始动能（单位是Hartree）。
import mlatom as ml
# Generate random velocities given geometry
mol = ml.data.molecule.from_xyz_file('h2_md_kreg_init.xyz')
init_cond_db = ml.generate_initial_conditions(molecule=mol,
generation_method='random',
initial_temperature=300, # Set initial temperature in K
#initial_kinetic_energy=4.75E-4, # Or set initial kinetic energy in Hartree
number_of_initial_conditions=1)
备注
在MLatom中，如果定义 generation_method='random' ，平移动量和角动量(速度)默认被删除(对非线性分子 \(N_\text{DOF}=3N_\text{atoms}-6\) ）。如果想要保留角动量，可以使用 eliminate_angular_momentum=False ，( \(N_\text{DOF}=3N_\text{atoms}-3\) )。默认情况下，对于双原子分子， \(N_\text{DOF}=3N_\text{atoms}-5\) ，而对于其他线性分子( \(N_\text{atoms}>2\) )， \(N_\text{DOF}=3N_\text{atoms}-6\) (但MLatom会检查分子的线性并输出警示信息)。您可以使用 degrees_of_freedom 参数来覆盖MLatom中自由度的默认设置。
除了随机速度之外，MLatom还提供了 准经典分子动力学 所需的两种简正模采样方法。一种是基于Newton-X中实现的非锐化Wigner分布（参见该软件 论文 中的公式20）。这些代码在Python中被修改和重写。另一种是基于 此文 描述的谐波量子玻尔兹曼分布。下面我们将展示如何使用简正模采样，带有分子频率和简正模的文件可在此处下载： ethanol_freq.json
import mlatom as ml
eqmol = ml.data.molecule()
eqmol.load('ethanol_freq.json',format='json')
# Sample in Wigner distribution
init_cond_db = ml.generate_initial_conditions(molecule=eqmol,
generation_method='wigner',
initial_temperature=300, # Set initial temperature in K
number_of_initial_conditions=1)
# Sample in harmonic quantum Boltzmann distribution
init_cond_db = ml.generate_initial_conditions(molecule=eqmol,
generation_method='harmonic-quantum-boltzmann',
initial_temperature=300, # Set initial temperature in K
number_of_initial_conditions=1)
备注
简正模采样需要频率和简正模。确保分子对象具有这些属性(通常只需运行 频率计算 )。对于过渡态，我们建议使用谐波量子玻尔兹曼分布(确保频率计算后只有一个虚频率)。反应坐标是单独处理的。
方法或模型
用户应选择一种模型或方法来计算MD的力。这里提供了一个预训练的ML模型 h2_kreg_energies.unf ，也可以选择其他模型或方法，参见以下示例中被注释掉的内容：
# Initializing model
model = ml.models.kreg(model_file='h2_kreg_energies.unf',ml_program='MLatomF')
# or another ML model, e.g., MACE:
#mymodel = ml.models.mace(model_file='mace.pt')
# or choose one of the predifined (automatically recognized) methods
#model = ml.models.methods(method='ANI-1ccx')
# or QM method, e.g., B3LYP with Gaussian
#model = ml.models.methods(method='B3LYP/6-31G*', program='Gaussian')
恒温器
若不使用NVE系综(微正则系综，系统具有恒定粒子N、恒定体积V和恒定能量E)，则应指定恒温器。这里，我们在NVT（正则系综，系统具有恒定的粒子N、恒定的体积V且温度在T周围波动）中运行MD，需要先初始化Nosé-Hoover恒温器。
# Initializing thermostat
nose_hoover = ml.md.Nose_Hoover_thermostat(temperature=300, molecule=init_mol)
备注
除了线性分子外，Nosé-Hoover恒温器的自由度基本上与初始条件下的自由度相同。任意线性分子默认的自由度为 \(3N_\text{atoms}-5\) 。用户仍然可以使用 degrees_of_freedom 参数来覆盖MLatom中的默认设置。
Nosé-Hoover恒温器的高级设置，如Nosé-Hoover链长度( nose_hoover_chain_length )，多个时间步长( multiple_time_step )，Yoshida-Suzuki步数(number_of_yoshida_suzuki_steps)和Nosé-Hoover链频率( nose_hoover_chain_frequency )也可以更改。更多信息请参考使用手册。
Andersen恒温器也可以在MLatom中使用。
andersen = ml.md.Andersen_thermostat(temperature=300)
备注
Andersen恒温器中的随机碰撞适用于所有 \(3N_\text{atoms}\) 自由度。
您可以使用 gamma 来修改Andersen恒温器中的碰撞率( \(\text{fs}^{-1}\) )，其默认值为0.2。
动力学
最后一步是运行分子动力学。用户应提供最长模拟时间(以fs为单位)和时间步长(以fs为单位)。在这里，最长模拟时间为30 fs，时间步长为0.3 fs。
# Run molecular dynamics
dyn = ml.md(model=model,
molecule_with_initial_conditions=init_mol,
thermostat=nose_hoover,
ensemble='NVT',
time_step=0.3,
maximum_propagation_time=30.0,
dump_trajectory_interval=1 # e.g., you can dump every time step (by default nothing is dumped!)
filename='mytraj.h5' # file where you want to save the traj
format='h5md' # format of dumped traj (supported formats: plain_text, h5md, json)
)
轨迹被保存为 ml.md 的属性，能够以纯文本文件和h5md格式的文件输出。
# Dump trajectory
traj = dyn.molecular_trajectory
traj.dump(filename='traj', format='plain_text')
traj.dump(filename='traj.h5', format='h5md')
汇总
综合上述4个步骤，整个工作流程即从设置初始条件到获得动力学轨迹文件的代码如下所示:
import mlatom as ml
# Use user-defined initial conditions
mol = ml.data.molecule.from_xyz_file('h2_md_kreg_init.xyz')
init_cond_db = ml.generate_initial_conditions(molecule=mol,
generation_method='user-defined',
file_with_initial_xyz_coordinates='h2_md_kreg_init.xyz',
file_with_initial_xyz_velocities='h2_md_kreg_init.vxyz')
init_mol = init_cond_db[0]
# Initializing model
model = ml.models.kreg(model_file='h2_kreg_energies.unf',ml_program='MLatomF')
# Initializing thermostat
nose_hoover = ml.md.Nose_Hoover_thermostat(temperature=300, molecule=init_mol)
# Run molecular dynamics
dyn = ml.md(model=model,
molecule_with_initial_conditions=init_mol,
thermostat=nose_hoover,
ensemble='NVT',
time_step=0.3,
maximum_propagation_time=30.0)
# Dump trajectory
traj = dyn.molecular_trajectory
traj.dump(filename='traj', format='plain_text')
traj.dump(filename='traj.h5', format='h5md')
print(f"Number of steps in the trajectory: {len(traj.steps)}")
XYZ坐标保存在 traj.xyz 中，XYZ速度保存在 traj.vxyz 中。动能保存在 traj.ekin 中，势能保存在 traj.epot 中，总能量保存在 traj.etot 中，能量梯度保存在 traj.grad 中。轨迹共有101步(包括初始条件)，最后一步的几何构型为:
2
H 0.0000000000000 0.0000000000000 0.8764323118216
H 0.0000000000000 0.0000000000000 0.1235676881784
停止函数
用户可以在MD中定义一个停止函数，以便在满足某些条件时停止(仅可通过Python API使用)。例如，下面展示了如何在键长小于0.9 Å时停止MD。
import mlatom as ml
# Use user-defined initial conditions
mol = ml.data.molecule.from_xyz_file('h2_md_kreg_init.xyz')
init_cond_db = ml.generate_initial_conditions(molecule=mol,
generation_method='user-defined',
file_with_initial_xyz_coordinates='h2_md_kreg_init.xyz',
file_with_initial_xyz_velocities='h2_md_kreg_init.vxyz')
init_mol = init_cond_db[0]
# Initializing model
model = ml.models.kreg(model_file='h2_kreg_energies.unf',ml_program='MLatomF')
# Initializing thermostat
nose_hoover = ml.md.Nose_Hoover_thermostat(temperature=300, molecule=init_mol)
# User-defined stop function
def stop_function(molecule):
dist_matrix = molecule.get_internuclear_distance_matrix()
# Get bond length
bond_length = dist_matrix[0,1]
if bond_length < 0.9:
return True
else:
return False
# Run molecular dynamics
dyn = ml.md(model=model,
molecule_with_initial_conditions=init_mol,
thermostat=nose_hoover,
ensemble='NVT',
time_step=0.3,
maximum_propagation_time=30.0,
stop_function=stop_function,
stop_function_kwargs={})
# Dump trajectory
traj = dyn.molecular_trajectory
traj.dump(filename='traj_stop_function', format='plain_text')
traj.dump(filename='traj_stop_function.h5', format='h5md')
print(f"Number of steps in the trajectory: {len(traj.steps)}")
轨迹有4步，最后一步的几何构型是:
2
H 0.0000000000000 0.0000000000000 0.9417784217628
H 0.0000000000000 0.0000000000000 0.0582215782372
备注
每一步中，分子都会被传递给停止函数来检查是否需要停止MD。用户定义的停止函数中的第一个参数必须接受一个 molecule 对象。
问题或建议
如果您有进一步的问题、批评和建议，我们很乐意通过 电子邮件 、 Slack (首选)或微信(请发送电子邮件请求将您添加到XACS用户支持组)以英文或中文接收您的意见。