溶剂化自由能的计算

背景

将溶质从气相转移到溶剂的自由能称为溶剂化自由能，这个概念对于理解大量的化学和物理过程非常重要。例如，它可以决定溶解度，即当溶剂与溶质蒸汽接触时，溶剂溶解在溶剂中的量。了解一种化合物在两种或两种以上溶剂中的相对溶解度是设计溶剂萃取工艺（一种主要的液相分离技术）的关键。相对溶解度最常使用的参数是log P，它量化了化合物在水和正辛醇中的分配比例。log P 为正值表明溶质优先在正辛醇中溶剂化，这是非极性化合物的典型特征。log P 被广泛用于预测药物的吸附程度以及体内农药和环境毒素的活性。

可以使用热力学积分方法直接模拟出溶剂化自由能。在该方法中，溶质与溶剂的相互作用通过若干步从零逐渐增加到完全相互作用（反之亦然）。对于每个耦合强度，进行模拟并计算溶剂化自由能的导数。在计算所有导数之后，可以得到自由能作为积分。由于每一步都需要长时间的模拟计算（通常为100 ps），自由能计算非常昂贵。

实际上，将通过三次计算运行执行溶剂化自由能计算：

从真空中的溶质分子开始，电荷逐渐减少到零，同时保持所有其他相互作用不变。这个过程中的自由能变化称为理想贡献。

然后，通过几个步骤打开范德华相互作用，将无电荷分子与溶剂耦合。这个过程中的自由能变化称为范德华贡献。

最后，在溶质分子上重新引入电荷，溶质分子现在浸没在溶剂中。这一步的自由能是静电贡献。然后，总溶剂化自由能为理想、范德华和静电贡献之和。

介绍

本教程介绍如何执行溶剂化自由能计算。可以计算丙酸在正辛醇中的溶剂化自由能。丙酸是一种液体脂肪酸，天然存在于汗液、乳制品和细菌发酵产物中。它常以丙酸盐的形式作为面包中的防霉剂和香水中的成分。

本教程计算辛醇中的溶剂化自由能。结合水的实验溶剂化自由能-6.36 kcal/mol (Wolfenden, 1981)，可以确定丙酸的log P 值。

注意：由于本教程中的计算需要非常长时间的模拟运行，因此提供了长时间输出构型的计算结果以供比较。示例工程SolvationFreeEnergy.stp位于Materials Studio安装目录下的share\Examples\Projects\Forcite文件夹中。

对于非管理员的Windows用户，可以将SolvationFreeEnergy.stp工程和相关的SolvationFreeEnergy_Files文件夹复制到具有写入权限的位置。然后打开SolvationFreeEnergy.stp工程的新副本。

本教程包括如下部分：

开始

绘制丙酸和正辛醇的分子构型

含有丙酸分子的正辛醇无定形晶胞的构建

平衡结构构型

运行溶剂化自由能计算

注意：为了和本教程中的参数保持一致，可以使用Settings Organizer对话框将工程中所有参数都设置为BIOVIA的默认值。

1、开始

①　首先启动Materials Studio并创建一个新工程。

②　打开New Project对话框，输入SolvationFreeEnergy作为工程名，单击OK按钮。

③　新工程将以SolvationFreeEnergy为工程名显示于Project Explorer中。

2、绘制丙酸和正辛醇的分子构型

①　首先，构建丙酸的分子结构：

从菜单栏中选择File | New...，打开New Document对话框。选择3D Atomistic，单击OK按钮。

将文档重命名为Propionic acid。

②　使用Sketch工具条上的工具，绘制丙酸分子，确保单击调整氢原子（加氢）和结构整理按钮。在文档中单击鼠标右键，然后选择Display Style，打开Display Style对话框。

在Atom选项卡中，Display style选择CPK并关闭对话框。

③　为了为后续计算准备分子结构，应使用电荷组优化几何结构。一旦定义了电荷组，它们就可以在后续计算中使用，而无需重新计算。

④　从菜单栏中选择Modules | Forcite | Calculation，打开Forcite Calculation对话框。将Task更改为Geometry Optimization。

3、含有丙酸分子的正辛醇无定形晶胞的构建

①　接下来，使用Amorphous Cell模块创建溶剂化自由能计算的输入结构。该结构包含一个丙酸分子和几个溶剂分子。

②　从菜单栏中选择Modules | Amorphous Cell | Calculation，打开Amorphous Cell Calculation对话框。

③　在Composition列表中，单击Molecule列，从下拉列表中选择优化后的结构，即Propionic acid Forcite GeomOpt文件夹中的Propionic acid.xsd文件。确保Loading是1。

④　单击下一行，选择优化后的结构n-Octanol.xsd。将正辛醇Loading设置为100。

⑤　正辛醇密度的实验值是0.824 g/cm3。以该密度构建体系，盒子长度约为3 nm。在Amorphous Cell计算中使用与在Forcite中相同的能量设置。

⑥　在Setup选项卡中，将Density定义为0.824。单击Options...按钮，打开Amorphous Cell Options对话框。勾选Include non backbone torsions复选框，关闭对话框。

⑦　在Energy选项卡中，从Forcefield下拉列表中选择COMPASSIII。选择Group based作为Electrostatic 的求和方法。

⑧　在Project Explorer中选择工程根目录，单击Run按钮，开始计算任务。关闭对话框。

⑨　注意：由于输入结构上存在有效的电荷组，Amorphous Cell不会重新将其进行计算；不需要取消自动计算电荷组。范德华相互作用的基于电荷组截断Group based的求和方法在Amorphous Cell中是不可用的，可以用基于原子截断Atom based的求和方法来代替。

⑩　将在Project Explorer中创建一个新文件夹Propionic acid AC Construct。计算完成后，该文件夹的轨迹文档中将包含溶剂化分子结构。Amorphous Cell总是输出一个轨迹文档，从而方便存储多帧构型。但是，在下例中，处理结构文档更方便，因此先复制结构。

⑪　选中工程根目录，创建一个新的3D原子结构文档3D Atomistic Document，将其重命名为PA_nOct。

⑫　右键单击Amorphous Cell结构Propionic acid.xtd，选择Copy。右键单击新文档PA_nOct.xsd并选择Paste。

⑬　在继续之后的步骤之前，保存并关闭所有打开的文件。

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

⑤　在Energy选项卡上，选择COMPASSIII作为Forcefield。确保Charges设置为Forcefield assigned。将Electrostatic和van der Waals求和方法更改为Group based。

⑥　单击Forcefield后的More...按钮，打开Forcite Preparation Options对话框。单击Charge groups后的More...按钮，打开Forcite Charge Groups对话框。

对于Method，选择Divide-and-conquer。关闭对话框。

⑦　单击Run按钮，开始计算任务。

将在Project Explorer中创建一个新文件夹Propionic acid Forcite GeomOpt。当计算结束后，在新文件夹中存储了优化过的结构。

输出结构分配了COMPASS电荷，并分成中性的电荷组。可以通过用Charge Group给分子着色，并用Charge标记原子来验证这一点。

对正辛醇使用相同的过程：

⑨　选中工程根目录，创建一个新的3D原子结构文档3D Atomistic Document，将其重命名为n-Octanol。

⑩　绘制一个正辛醇分子，保留其显示模式为默认设置。

⑪　如前所述，使用电荷组执行几何优化。这将使用与上面相同的设置，因此可以立即开始计算。

⑫　确保n-Octanol.xsd是当前文档，单击Forcite Calculation对话框上的Run按钮以启动计算。关闭对话框。

⑬　这将在Project Explorer中创建一个新文件夹n-Octanol Forcite GeomOpt。计算完成后，该文件夹中包含分配了电荷组的优化后结构。

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

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。

⑩　在结构中，双击溶质分子中的任何原子以选择分子。单击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>的列中。

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选项卡设置保留为在前面步骤中设置的值。在运行计算之前，通过创建一个集合，指定结构中的哪个分子是溶质分子。

打开数据表文件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−3 kcal/mol/K）。

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

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

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

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