Materials Studio官方教程：Forcite——溶剂化自由能的计算【1】

4、平衡结构构型

在开始溶剂化自由能计算之前，先平衡初始结构构型。Amorphous Cell计算执行了几何优化，但结构可能仍然包含需要动力学弛豫的应力集中位置和密度不均匀区域。选择恒温方法进行动力学计算是一种消除体系多余热量的方法。NHL恒温方法是平衡和生产构型运行的良好恒温方法。

使得PA_nOct.xsd为当前文档。打开Forcite Calculation对话框，从Task中选择Dynamics，单击More...按钮，打开Forcite Dynamics对话框。

选择NVT作为Ensemble。将Total simulation time设置为10 ps。在Thermostat选项卡中，选择Thermostat为NHL，将Q ratio设置为1.0。关闭对话框。

注意：默认Q ratio为0.01，可确保抑制Nosé恒温方法可能发生的温度振荡，从而实现更快的平衡。NHL恒温方法是Nosé恒温方法的延伸，具有自己的阻尼机制。这意味着可以在NHL中使用更高的Q ratio。当Q ratio为1时，恒温方法的影响很小，可以在平衡结构和生成构型的计算中使用相同的值。

在Energy选项卡中，对于Electrostatic和van der Waals求和方法均选择Group based。

单击Run按钮运行计算任务，关闭对话框。

将在Project Explorer中创建一个新文件夹PA_nOct Forcite Dynamics。在计算开始后很短的时间，就会出现图表和状态文本文件，记录计算过程。平衡构型的计算需要几分钟的时间才能完成，取决于计算机的性能。

提示：可以通过并行运行动力学计算来减少计算所需的时间，这需要使用多个CPU核和许可证。

在实际计算中，需要运行更长时间构型平衡计算，比如100 ps，然后用大约相同时间的恒压方法进行构型平衡，以调整晶胞的尺寸。

注意：如果不想等待计算结束就进行之后的步骤，可以从实例文件夹中导入已经计算好的结构，在Examples\Projects\Forcite\SolvationFreeEnergy.stp工程中。该结构已经过100 ps NVT动力学平衡和100 ps NPT动力学平衡。由于使用的结构的差异，结果可能会有所不同。

继续下个步骤之前，请保存所有打开的文档并清除工作区文件。

从菜单栏中选择File | Save Project，然后选择Window | Close All。

5、运行溶剂化自由能计算

对PA_nOct Forcite Dynamics文件夹中平衡运行的最终结构，或示例文件夹中的等效结构，继续接下来的计算。

在Forcite dynamics文件夹中打开3D Atomistic文件（非轨迹文件）PA_nOct.xsd。

此3D原子文档包含动力学计算的最终结构。

选择Modules | Forcite | Calculation，打开Forcite Calculation对话框。选择Task为Solvation Free Energy，单击More...按钮，打开Forcite Solvation Free Energy对话框。

可以把总溶剂化自由能计算为Ideal、van der Waals和Electrostatic三个贡献项的总和。

每个贡献项都需要单独计算溶剂化自由能。可以对所有贡献项使用相同的输入结构，并按任意顺序运行计算。也可以在一次计算中运行所有三个贡献项，这将在所有运行中使用相同的设置顺序执行。

本教程将依次地运行各个贡献项计算，从理想贡献开始，始终使用默认的热力学积分算法。要减少本教程中计算所需的时间，请将平衡的时间步数减少到1000，生产构型的时间步数减少到5000。在真实的模拟中，将运行更长的时间。示例文件夹中的模拟使用了50000步平衡和100000步构型生产。本教程将耦合参数从0更改为1，这意味着电荷将从全COMPASS电荷减少到零。也可以反向运行耦合参数，结果平均值保持不变。对于理想贡献，可以减少耦合参数步数。

将Contribution设置为Ideal。Equilibration steps设置为1000，Production steps设置为5000。耦合参数的Steps设置为5。

可以将Dynamics和Geometry Optimization设置以及Energy选项卡设置保留为在前面步骤中设置的值。在运行计算之前，通过创建一个集合，指定结构中的哪个分子是溶质分子。

在结构中，双击溶质分子中的任何原子以选择分子。单击Add以创建名为SoluteAtoms的集合。

单击Run按钮开始计算。

将在Dynamics文件夹中创建一个新的文件夹PA_nOct Forcite FreeEnergy。较短时间后，将会出现一个状态文本文档，记录计算的进度。理想贡献的计算只需要几分钟就可以完成，因为它仅涉及一个空盒子中的单个分子。计算完成后，将返回报告文档、图表文档和数据表文档。

打开报告文档PA_nOct.txt并向下滚动到包含溶剂化自由能结果的部分。

丙酸溶剂化自由能的理想贡献约为25.5 kcal/mol。由于这是去除分子上电荷的自由能损耗，因此接近于优化后结构的静电能的相反数，如上所述。可以通过查找丙酸几何优化运行报告文档中的静电能（约-26.2 kcal/mol）来验证这一点。

注意：本教程中使用的自由能值为更长时间计算后的结果，可在示例文件夹中读取。由于本教程中的模拟时间要短得多，因此结果可能会有所不同。

溶剂化自由能通过模拟中计算的自由能导数的积分来确定。图表显示了该积分与耦合参数之间的函数关系。

打开图表文档PA_nOct FreeEnergy.xcd。确认耦合参数为1处的值与报告的溶剂化自由能匹配。

图表中的数据也包含在数据表文件的列A(lambda)中。实际自由能导数值包含在标有<d E/d lambda>的列中。

打开数据表文件PA_nOct.std。通过单击首行标题A选择第一列。按住CTRL键并单击行标题D以选择第四列。单击Quick Plot按钮。

该图显示了耦合参数为1处从约50 kcal/mol的值到0的直线。自由能对应于线下的面积，在这种情况下是初始值的一半。

在继续进行其他贡献项的计算之前，请保存所有打开的文档并清除工作区。

从菜单栏中选择File | Save Project，然后选择Window | Close All。

继续讨论静电对溶剂化自由能的贡献。由于静电的自由能导数近似为线性，因此可以使用与理想贡献相同数量的耦合参数步数。

打开与上一步相同的输入结构文档PA_nOct.xsd。

在Forcite Solvation Free Energy对话框中，选择Contribution为Electrostatic。

单击Run按钮开始计算。

将在Dynamics文件夹中创建一个新文件夹，即PA_nOct Forcite FreeEnergy (2)。较短时间后，将会出现一个状态文本文档，记录计算的进度。静电贡献的计算需要更长的时间，因为它涉及系统中的所有分子。

注意：Job Explorer中报告的进度是单个动力学运行的进度，而不是整个溶剂化自由能运行的进度。每次平衡或生产运行后，该进度将重置为0。

提示：溶剂化自由能计算中的每个贡献项的计算是独立的。如果有适当的许可证和可用的CPU核数，这些许可证和CPU可以同时运行。

如前所述，完成后返回报告、图表和数据表文件。报告文件包括静电贡献项计算结果。

打开报告文档PA_nOct.txt并向下滚动到包含溶剂化自由能数据的部分。

丙酸在辛醇中的静电贡献约为-29.4 kcal/mol，接近理想贡献，接近优化后分子的静电能。理想贡献和静电贡献之和约为-3.9 kcal/mol，表明丙酸和正辛醇之间的静电相互作用储存了约3.9 kcal/mol的能量。

注意：本教程中使用的自由能值为更长时间计算后的结果，可在示例文件夹中读取。由于本教程中的模拟时间要短得多，因此结果可能会有所不同。

在继续进行最终贡献的计算之前，请保存所有打开的文档并清除工作区。

从菜单栏中选择File | Save Project，然后选择Window | Close All。

范德华贡献项是最难确定的，因为自由能导数是非线性的。这意味着需要有足够数量的耦合参数步数来确定曲率。可以使用与静电贡献项相同数量的平衡和构型生产步骤，但增加耦合参数的步数。

打开与上一步相同的输入结构文档PA_nOct.xsd。

在Forcite Solvation Free Energy对话框中，选择Contribution为van der Waals。对于耦合参数Steps，设置为10并关闭对话框。

单击Run按钮，启动计算并关闭对话框。

将在Project Explorer中创建一个新文件夹PA_nOct Forcite FreeEnergy (3)。由于使用的步骤是静电贡献项的两倍，因此此运行所需的时间大约是静电贡献的两倍。

完成后，打开生成的图表文档。

打开图表文档PA_nOct FreeEnergy.xcd。

注意，自由能与耦合参数之间的函数是非单调的。溶质分子最初浸入溶剂后，自由能变为正值，表明在辛醇溶剂中形成空腔的能量。一旦形成空腔，与溶剂的吸引色散作用占主导地位，自由能变为负值。最终值应为负值，约为-2.8 kcal/mol。

打开文本文档PA_nOct.txt，确认范德华项对溶剂化自由能的贡献约为-3 kcal/mol。

现已得到了溶剂化自由能的所有贡献。

总计为：

25.5 - 29.4 - 2.8 = -6.7 kcal/mol

提示：通过将Contribution选择为ALL，可以在单次计算中运行所有贡献项。示例文件夹中给出的计算就是这样进行的。这样的计算返回单个图表文档，其中包含所有的图像，和单个数据表，其中每个贡献都有单独的表单。

最后，要计算丙酸的log P，可以使用水作为溶剂代替辛醇重复上述计算。利用丙酸的实验水合自由能−6.36 kcal/mol (Wolfenden, 1981)。由于该值（负值）的绝对值略小于辛醇的计算值，因此可以预测丙酸可溶于辛醇，并且预计log P为正值。

log P通过下式获得：

这里：

0.434对应于10log e。

R是气体常数（1.987×10⁻³ kcal/mol/K）。

T是模拟中的温度（298 K）。

Ax是溶剂x中的溶剂化自由能。

代入上式，得出log P=0.25，与实验值0.246 (Vitas-M Lab)非常一致。

从菜单栏中选择File | Save Project。

本教程到此结束。

参考文献：

Wolfenden, R.; Andersson, L.; Cullis, P. M.; Southgate, C. C. B., "Affinities of amino acid side chains for solvent water", Biochemistry, 20(4), 849–855 (1981).

Log P data obtained from Vitas-M Lab (ID: STL168039), https://vitasmlab.biz/finded-

stk=STL168039.

Materials Studio是久负盛名计算模拟软件，问世20余年来，经过不断地迭代优化，使其功能异常强大，极易上手，初学者只需通过简单的参数设置和点击鼠标就能完成DFT计算。其计算可靠性久经考验，备受Nature、Science等顶级期刊认可。

华算科技和Materials Studio官方代理深圳浦华系统联合推出Materials Studio建模、计算、分析课程。课程专为零基础学员设计，沿着理论讲解、模型搭建、性质计算、结果分析层层递进讲解，带你快速入门DFT计算。课程极度注重实操，全程线上直播，提供无限回放，课程群在线答疑。（详情点击下方图片跳转）

识别下方二维码报名，或者联系手机13005427160。