一种准确计算气体分子绝对自由能的方法

摘要

本发明公开了一种准确计算气体分子绝对自由能的方法，包括以下步骤：S1、降低气体分子中柔性二面角势能面的能垒；S2、气体分子绝对自由能的计算。本发明实现了分子中柔性二面角自动识别和势能面扫描，实现了拟合二面角势能面，并自动修改二面角参数，以降低势能面至特定值，实现了新的WHAM程序，可以从biased和unbiased的势能reweight得到实际的自由能，设计了一套可行的计算流程，实现气体分子在逐步去除各类能量时，势能的平缓变化；将上述所有流程衔接在一起，实现了气体分子绝对自由能的自动化计算。

说明书

技术领域

本发明属于分子动力学模拟技术领域，具体涉及一种准确计算气体分子绝对自由能的方法。

背景技术

近年来，分子模拟技术的应用越来越广泛，对于自由能计算的需求也越来越大。当前，由于晶体结构的稳定较好，采样空间相对较小，已有一些方法可以比较准确地预测固体分子的绝对自由能，可应用于药物分子的稳定实验晶型的预测和筛选。

而对于液体和气体分子的绝对自由能，还没有比较好的准确预测方法。这使得涉及到固液(如熔点、溶解度)、固气(如升华焓)等的性质还没有非常有效的预测方法。

气体分子和液体分子绝对自由能的计算，当前还没有非常有效的办法：

对于气体分子，由于其势能面极其复杂，模拟时很容易陷入局部能量最低点，从而影响自由能计算结果的收敛性。

对于液体分子，情况更加复杂，不仅需要克服液体分子本身的势能面的影响，还要克服分子快速扩散带来的采样问题。因此，要准确计算液体分子的绝对自由能就更难实现。

发明内容

针对上述技术问题，本发明的目的在于提供一种准确计算气体分子绝对自由能的方法，能准确计算固体和气体的绝对自由能，通过溶剂化自由能的计算，就可以间接计算出液体分子的绝对自由能。

为实现上述目的，本发明提供如下技术方案：

一种准确计算气体分子绝对自由能的方法，包括以下步骤：

S1、降低气体分子中柔性二面角势能面的能垒；

S1.1准备好要计算分子的坐标文件，并通过相关工具抓取力场参数；

S1.2调用yoda程序，加载分子坐标和力场，搜索分子中的所有柔性二面角；

S1.3对找到的所有柔性二面角进行势能面扫描；

S1.4对势能面中存在能垒大于5kcal/mol的二面角势能面进行拟合，修改对应二面角参数，降低能垒；

S1.5重复步骤S1.3和S1.4，直至所有柔性角的势能面中不存在大于5kcal/mol的能垒；

S1.6生成新的力场参数文件；

S2、气体分子绝对自由能的计算；

S2.1使用新的力场参数，建一个20nm×20nm×20nm的模拟盒子，在一系列温度下，进行NVT模拟，使用cut-off方式计算静电作用力，步长1fs，模拟10ns；

S2.2选取两个参考温度下的平衡结构，作为自由能模拟的初始结构；

S2.3逐步限制一个中心重原子的位置，力常数从0变到10000kJ/mol/nm

S2.4逐步限制其它重原子的位置，力常数从0变到10000kJ/mol/nm

S2.5逐步限制非重原子的位置，力常数从0变到10000kJ/mol/nm

S2.6逐步增大所有原子的位置限制，力常数从10000变到200000kJ/mol/nm

S2.7逐步去除所有原子的分子内(静电、范德华、键、角、二面角)相互作用；

S2.8上述模拟正常结束后，全部使用原始力场参数对轨迹进行rerun，得到unbiased的势能；

S2.9使用修改后的WHAM程序，加载各个λ或T窗口biased和unbiased能量，得到最终的reweight后的自由能；

S2.10使用reweight后的Δf(T)和PSCP自由能，以及相关理论推导值，计算最终的绝对自由能随温度的变化。

其中，步骤S1.1所述的相关工具，为antechamber或CHARMM-GUI。

步骤S2.1所述的一系列温度为100～400K范围。

步骤S2.2选取两个参考温度为200K和300K。

与现有技术相比，本发明的有益效果是：

(1)实现了分子中柔性二面角自动识别和势能面扫描。

(2)实现了拟合二面角势能面，并自动修改二面角参数，以降低势能面至特定值。

(3)实现了新的WHAM程序，可以从biased和unbiased的势能reweight得到实际的自由能。

(4)设计了一套可行的计算流程，实现气体分子在逐步去除各类能量时，势能的平缓变化。

(5)将上述所有流程衔接在一起，实现了气体分子绝对自由能的自动化计算。

附图说明

图1为实施例丁烷分子结构示意图；

图2为实施例的丁烷分子初始C_C1_C2_C3二面角势能面图；

图3为实施例扫描得到的势能面；

图4是本发明气体分子绝对自由能的计算流程；

图5使用原力场得到的绝对自由能曲线；

图6为本实施例使用新力场经WHAM程序reweight得到的绝对自由能曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本实施例以丁烷(butane)分子绝对自由能的计算步骤，如图1所示的丁烷结构式：

(1)准备好丁烷分子的坐标(pdb)文件，使用antechamber抓取对应的GAFF2力场参数；

(2)调用yoda程序，加载分子坐标和力场，搜索分子中的所有柔性二面角；

(3)对找到的所有柔性二面角进行势能面扫描；

丁烷分子只有一个柔性二面角(C_C1_C2_C3)如图2所示。

(4)对势能面中存在能垒大于5kcal/mol的二面角势能面进行拟合，修改对应二面角参数，降低能垒；

拟合函数：

f(x)＝k1*(1+cosx)+k2*(1-cos2x)+k3*(1+cos3x)+k4*(1+cos4x)+k6*(1+cos6x)+b

(5)重复步骤(3)和(4)，直至所有柔性角的势能面中不存在大于5kcal/mol的能垒；修改C_C1_C2_C3二面角参数后，扫描得到的势能面，如图3所示；

(6)生成新的力场参数文件；

气体分子绝对自由能的计算流程如图4所示。

(7)使用新的力场参数，建一个20nm×20nm×20nm的模拟盒子，在一系列温度(100～400K)下，进行NVT模拟，使用cut-off方式计算静电作用力，步长1fs，模拟10ns；

通过这一系列温度下模拟得到的势能数据H(T)，可以积分得到无量纲的自由能随温度的变化Δf(T)：

(8)选取两个参考温度(200K和300K，用于收敛性验证)下的平衡结构，作为自由能模拟的初始结构；

(9)逐步限制一个中心重原子的位置，力常数从0变到10000kJ/mol/nm

(10)逐步限制其它重原子的位置，力常数从0变到10000kJ/mol/nm

(11)逐步限制非重原子的位置，力常数从0变到10000kJ/mol/nm

(12)逐步增大所有原子的位置限制，力常数从10000变到200000kJ/mol/nm

(13)逐步去除所有原子的分子内(静电、范德华、键、角、二面角)相互作用(ΔA

(14)上述模拟正常结束后，全部使用原始力场参数对轨迹进行rerun，得到unbiased的势能；

通过上述模拟可以算出总的PSCP自由能变化(ΔG

(15)使用修改后的WHAM程序，加载各个λ或T窗口biased和unbiased能量，得到最终的reweight后的自由能；

用于reweight的WHAM程序的原理：

标准WHAM的计算公式：

(16)使用reweight后的Δf(T)和PSCP自由能，以及相关理论推导值，计算最终的绝对自由能随温度的变化。

最终绝对自由能的计算公式如下：

结果如图5和图6所示。

从计算结果可以看出：新的力场的结果结合reweight后的自由能，收敛性明显提升。这种提升在复杂一些的分子上表现更明显。

尽管已经示出和描述了本发明的实施例，对于本领域的普通技术人员而言，可以理解在不脱离本发明的原理和精神的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。