从分子动力学计算振动光谱
本教程将展示如何从 分子动力学 轨迹计算振动光谱（红外光谱和功率谱）。
我们还制作了该教程的视频，你可以在 bilibili 上找到它。
理论
红外光谱可以通过对偶极矩的自相关函数做傅里叶变换得到：
\(P(\nu)\propto\frac{\nu\text{tanh}(h\nu/2k_{\text{B}}T)}{n(\nu)}\int\langle\mathbf{\mu}(\tau)\mathbf{\mu}(t+\tau) \rangle _{\tau}e^{-i2\pi\nu t} \text{d}t\),
其中， \(\langle \mathbf{\mu}(\tau)\mathbf{\mu}(t+\tau) \rangle _{\tau}\) 是偶极矩的自相关函数， \(\nu\) 是频率， \(h\) 是普朗克常数， \(k_{\text{B}}\) 是玻尔兹曼常数， \(T\) 是温度， \(n(\nu)\) 是折射率，在气相中，该值为1。在做快速傅里叶变换之前，我们会先用Hann窗口函数处理自相关函数，用户可以调节自相关深度（autocorrelation depth）来控制谱峰的宽度。用户也可以在自相关函数后用零值填充（zero padding）来提高谱图的精度。
类似地，功率谱也可以从速度的自相关函数得到：
\(P(\nu)\propto\int \langle\mathbf{v}(\tau)\mathbf{v}(t+\tau)\rangle_{\tau}e^{-i2\pi\nu t} \text{d}t\),
其中， \(\langle\mathbf{v}(\tau)\mathbf{v}(t+\tau)\rangle_{\tau}\) 是速度的自相关函数。
从分子动力学轨迹计算振动光谱的理论部分可以参考以下两篇文章： DOI: 10.1039/C3CP44302G 和 DOI: 10.1039/C6CP05849C 。
MLatom中的计算细节可以参考 这篇文章（DOI: 10.1039/D3CP03515H） 。使用此功能时请引用该文章
备注
MLatom中的一些QM方法（例如使用Gaussian和MNDO程序的方法）会在用户计算单点能时自动提供偶极矩。除了AIQM1，其它所有的预训练机器学习模型都不能提供偶极矩。AIQM1提供的是ODM2*级别的偶极矩，且只有当使用MNDO（而非Sparrow）时有效。所以在计算红外光谱前请确认使用的方法是否能提供偶极矩。所有的方法都能用来计算功率谱，因为其只需要速度即可得到。
使用输入文件/命令行
接下来我们将展示如何通过输入文件或命令行计算振动光谱。输入文件中展示了如何从h5md格式的分子动力学轨迹计算乙醇的红外光谱。在这里我们提供了乙醇的轨迹，它是用AIQM1在300K的Nose-Hoover温标下得到的。轨迹的步长为0.5飞秒。分子动力学轨迹可以在这里下载： ethanol_traj.h5 （h5md格式）， ethanol_traj.xyz （xyz坐标）， ethanol_traj.vxyz (xyz速度)和 ethanol_traj.dp (偶极矩)。
MD2vibr # 1. requests vibrational spectra from MD trajectory
trajH5MDin=ethanol_traj.h5 # 2. file with MD trajectory in h5md format
dt=0.5 # 3. time step of trajectory; unit: fs
start_time=0 # 4. use trajectory starting from 0 fs
end_time=10000 # 5. use trajectory ending at 10000 fs
autocorrelationDepth=1024 # 6. autocorrelation depth; unit: fs
zeropadding=1024 # 7. zero padding length; unit: fs
output=ir # 8. output infrared spectrum
你会在文件夹里找到两个文件： ir.png 和 ir.npy 。前者是红外光谱的图片，后者是谱图的numpy数组。计算得到的红外光谱如下所示：
关键词 output 可以用来控制计算哪种振动光谱。 output=ir 计算红外光谱， output=ps 计算功率谱。下面是计算功率谱的输入文件。
MD2vibr # 1. requests vibrational spectra from MD trajectory
trajH5MDin=ethanol_traj.h5 # 2. file with MD trajectory in h5md format
dt=0.5 # 3. time step of trajectory; unit: fs
start_time=0 # 4. use trajectory starting from 0 fs
end_time=10000 # 5. use trajectory ending at 10000 fs
autocorrelationDepth=1024 # 6. autocorrelation depth; unit: fs
zeropadding=1024 # 7. zero padding length; unit: fs
output=ps # 8. output power spectrum
计算得到的功率谱如下所示：
除了h5md格式的轨迹，还可以直接提供偶极矩（用来计算红外光谱）和速度（用来计算功率谱）的文本文件，需要使用 trajdpin 或 trajVXYZin 关键词。下面是两种情况的例子。
MD2vibr # 1. requests vibrational spectra from MD trajectory
trajdpin=ethanol_traj.dp # 2. file with dipole moments in plain text
dt=0.5 # 3. time step of trajectory; unit: fs
start_time=0 # 4. use trajectory starting from 0 fs
end_time=10000 # 5. use trajectory ending at 10000 fs
autocorrelationDepth=1024 # 6. autocorrelation depth; unit: fs
zeropadding=1024 # 7. zero padding length; unit: fs
output=ir # 8. output infrared spectrum
MD2vibr # 1. requests vibrational spectra from MD trajectory
trajVXYZin=ethanol_traj.vxyz # 2. file with velocities in plain text
dt=0.5 # 3. time step of trajectory; unit: fs
start_time=0 # 4. use trajectory starting from 0 fs
end_time=10000 # 5. use trajectory ending at 10000 fs
autocorrelationDepth=1024 # 6. autocorrelation depth; unit: fs
zeropadding=1024 # 7. zero padding length; unit: fs
output=ps # 8. output power spectrum
备注
原则上我们推荐使用AIQM1方法。通常它在一定的计算代价下是最准确的。AIQM1的偶极矩是ODM2*级别的。但是偶极矩只有在使用MNDO程序时才能获得。我们推荐安装MNDO以解锁AIQM1的全部功能。XACS云计算平台使用的是Sparrow，目前不提供偶极矩，而且它只提供数值梯度，这限制了它在MD中的应用。AIQM1的另一个限制是它只能用于CHNO元素。ANI-1ccx是XACS平台上一个好的替代方案，但是它同样不能计算偶极矩，而且它只能用于计算中性闭壳层的基态体系。
其它关键词
关键词 start_time 和 end_time 可以控制选择轨迹的某一部分计算谱图。默认情况下，MLatom会着用整条轨迹。 autocorrelationDepth 控制自相关深度，它会影响峰的宽度； zeropadding 控制零填充的长度，它会影响谱图的分辨率。这两个设置的默认值是1024飞秒。MLatom会用谱图的最大值对谱图进行归一化。如果想要获得原始谱图，可以使用 normalize_intensity=0 。
使用API
接下来将展示如何用MLatom的API计算振动光谱。
import mlatom as ml
# Loading MD trajectory in h5md format
traj = ml.data.molecular_trajectory()
traj.load(filename='ethanol_traj.h5',format='h5md')
moldb = ml.data.molecular_database()
moldb.molecules = [step.molecule for step in traj.steps]
# Loading MD trajectory in plain text
# moldb = ml.data.molecular_database.from_xyz_file('ethanol_traj.xyz')
# moldb.add_xyz_vectorial_properties_from_file('ethanol_traj.vxyz','xyz_velocities')
# with open('ethanol_traj.dp','r') as f:
# lines = f.readlines()
# for ii in range(len(moldb)):
# moldb[ii].dipole_moment = [float(each) for each in lines[ii].strip().split()]
# Calculate vibrational spectra
vibr = ml.vibrational_spectrum(molecular_database=moldb,dt=0.5)
# Infrared spectrum
freqs_ir,ints_ir = vibr.plot_infrared_spectrum(filename='ethanol_ir.png', # file to save the figure
autocorrelation_depth=1024, # autocorrelation depth; unit: fs
zero_padding=1024, # zero padding; unit: fs
lb=0, # start time; unit: fs
ub=10000, # end time; unit: fs
normalize=True, # normalize the spectrum
return_spectrum=True, # get the spectrum array from the function
format='modulus') # use the modulus of the dipole moments
# Power spectrum
freqs_ps,ints_ps = vibr.plot_power_spectrum(filename='ethanol_ps.png', # file to save the figure
autocorrelation_depth=1024, # autocorrelation depth; unit: fs
zero_padding=1024, # zero padding; unit: fs
lb=0, # start time; unit: fs
ub=10000, # end time; unit: fs
normalize=True, # normalize the spectrum
return_spectrum=True) # use the modulus of the dipole moments
# Save the spectra
np.save('ethanol_ir_saved.npy',np.array((freqs_ir,ints_ir)))
np.save('ethanol_ps_saved.npy',np.array((freqs_ps,ints_ps)))
关键词与输入文件有些许不同。 lb 和 ub 可以控制轨迹的开始时间和结束时间。 autocorrelation_depth 和 zero_padding 分别控制自相关深度和零填充长度。 normalize=True 会让谱图归一化。使用 return_spectrum=True 可以让函数返回图谱的数组。在红外光谱中，还可以使用 format='modulus' 或 format='vector' 来选择使用偶极矩的模或XYZ向量。
问题或建议
如果您有进一步的问题、批评和建议，我们很乐意通过 电子邮件 、 Slack (首选)或微信(请发送电子邮件请求将您添加到XACS用户支持组)以英文或中文接收您的意见。