寻找分子在表面上的 低能构型

背景

小分子的势能面可能非常复杂，有许多局部能量极小值和一个全局能量最小值。有几种方法可用于确定全局最小值，包括蒙特卡罗算法和不同形式的分子动力学。淬火分子动力学通过附加的几何优化步骤执行标准分子动力学计算，在该步骤中，对轨迹文件中的每一帧执行几何优化。分子动力学可以对许多不同的低能构型进行有效采样。

介绍

本教程使用Forcite Plus对由小有机分子邻氯苯酚和金属氧化物二氧化钛表面组成的体系执行分子动力学。本教程由金红石相的TiO2晶体，切割二氧化钛表面，再将分子放置在表面上，对其进行优化，并运行淬火分子动力学。可以通过数据表寻找能量最低的构象，最后计算结合能。

本教程包括如下部分：

开始

氯苯酚结构的绘制与优化

构建并优化TiO2表面结构

表面上的分子的平衡位置计算

使用淬火动力学对构型进行采样

计算结合能

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

1、开始

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

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

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

2、氯苯酚结构的绘制与优化

绘制邻氯苯酚分子结构。

创建一个新的3D原子文档3D Atomistic Document。绘制一个邻氯苯酚分子，随后使用Adjust Hydrogen和Clean工具。

在Project Explorer，将3D原子文档重命名为chlorophenol.xsd。

所绘制的分子的结构并不精确，因此需要使用经典模拟工具进行几何优化。经典模拟工具需要一个力场，用不同类型原子之间的距离、角度等函数的总和来描述原子之间的力。

这里将使用Forcite经典模拟引擎和COMPASS力场。

单击Modules工具条上的Forcite按钮

在下拉列表中选择Calculation，或在菜单栏中选择Modules | Forcite | Calculation。

将打开Forcite Calculation对话框。

在该对话框中，可以选择执行何种计算任务，计算精度Quality根据所选择的计算任务控制不同的参数。在Energy选项卡中，可以选择力场、电荷计算方法，并自定义对非键项的处理方式。

将Task修改为Geometry Optimization。在Energy选项卡中，将Forcefield修改为COMPASSIII。

COMPASS是一个已经计算了电荷的、经过充分验证的力场。因此，当选择COMPASS时，Charges选项将自动更改为由力场分配Forcefield assigned。非键设置默认为基于原子截断Atom based。这适用于小分子，但当稍后引入TiO2表面时，会将其更改为Ewald。

单击Run按钮。

这将启动计算任务并将其显示在Job Explorer中。计算任务将在Job Explorer中创建一个名为chlorophenol Forcite GeomOpt的新文件夹。计算任务完成后，将显示提示，计算结果文件夹中包含6个文件。

chlorophenol.xsd：包含初始结构优化后的几何结构

chlorophenol - Calculation：包含计算任务设置的xml状态文件。单击此文件将打开Forcite Calculation对话框，其中包含计算中指定的设置。

chlorophenol Convergence.xcd：能量和梯度法线的变化图。

chlorophenol Energies.xcd：焓变图。急剧下降表明初始绘制的较粗糙的几何构型已成功优化。

Status.txt：包含实时更新状态。

chlorophenol.txt：初始设置的文本版本，以及初始和最终结构的能量项分解。

3、构建并优化TiO2表面结构

Materials Studio的结构库中包含大量各种结构，其中包含许多常见的金属氧化物。存在三种晶型的二氧化钛。在这里，使用金红石相。

在Project Explorer中，右键单击工程根目录，从弹出的快捷菜单中选择Import...，打开Import Document对话框。导航至Structures\metal-oxides并导入TiO2_rutile.xsd文件。

在切割表面之前，优化晶体结构的几何构型。

在Forcite Calculation对话框的Setup选项卡中，单击More...按钮，打开Forcite Geometry Optimization对话框。选择Optimize cell复选框，关闭对话框。

作为一个三维周期性结构，可以选择优化晶胞和原子位置，或仅优化原子位置。在本例中，同时优化晶胞和原子位置。

在Properties Explorer中，将Filter改为Lattice 3D。

二氧化钛的金红石相属于P42/MNM空间群。Forcite可以在带有对称性约束的条件下优化结构。

单击Run按钮，关闭Forcite Calculation对话框。

将创建一个新文件夹，TiO2_rutile Forcite GeomOpt。当计算结束时，新文件夹包含与氯苯酚文件夹中类似的一系列文件。另外，由于这是一个周期性结构，Forcite也给出了晶胞参数和密度的变化。

注意：TiO2_rutile.txt输出文件会报告某些能源贡献缺少参数。这些不会影响本教程的结果，但可以通过手动分配力场类型，然后删除结构中的所有化学键来避免这些提示消息。聚合物与金属氧化物表面相互作用（Polymer interactions with a metal oxide surface）教程的第2步使用这种方法。

现在已经优化了结构，切割晶体并暴露活性110表面。

提示：如果不知道活性表面，可以使用Morphology模块进行搜索。

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

打开优化过的TiO2_rutile.xsd文件。选择Build | Surfaces | Cleave Surface，打开Cleave Surface对话框。将Cleave plane更改为1 1 0，按下TAB键。

在TiO2_rutile.xsd中，将显示一个蓝色盒子，以蓝色实线显示切割平面，点线表示切割深度。

使用微调按钮将Fractional Top增加到0.4，Fractional Thickness增加到3.0。单击Cleave按钮，并关闭对话框。

将打开一个新的3D原子文档TiO2_rutile (1 1 0).xsd。表面由白色矩形表示，包含钛原子和氧原子。

现在要优化表面。由于已经优化了块体晶体结构，因此可以固定底层的几何构型，因为它们表示块体结构。

单击3D Viewer Selection Mode工具。重新定位结构的方向，使得表面在屏幕顶端。

选择底部的Ti和O原子混合层以及其正上方和下方的O原子层部分。从菜单栏中选择Modify | Constraints。

打开Edit Constraints对话框。

勾选Fix Cartesian position复选框，关闭对话框。

单击3D Viewer以取消选中所选部分。右键单击3D Viewer，从弹出的快捷菜单中选择Display Style，打开Display Style对话框。在Atom选项卡的Coloring部分，将Color by选项更改为Constraint。

结构应与下图类似，被固定的原子显示为红色，未被固定的原子显示为灰色。

在Display Style对话框，将Color by选项重新改回Element，关闭对话框。

现在将优化表面结构。

打开Forcite Calculation对话框，单击Run按钮。

将再次显示包含结果文件的新文件夹。本部分的最后一步是确定表面的大小。表面必须足够大，以容纳有机分子，同时尽量减少与镜像分子之间的相互作用。因此，首先使用距离监视工具计算氯酚分子的长度尺寸。

打开优化过的chlorophenol.xsd文件。单击Measure/Change按钮

的下拉箭头，选择Distance。选择Cl原子及与之对位的H原子。

将显示一个值约为5.6 Å的距离测量。假设每个分子周围有2 Å区域，则分子间相互作用可以达到最小化，这表明单位晶胞大小为5.6 + 2 × 2或约为9.6 Å。现在，需要确定必须生成多少个当前表面的副本才能实现此目的。

使得优化好的金红石相表面为当前文档。在Properties Explorer中，将Filter更改为Lattice 2D。

U和V的晶胞参数分别约为3.0和6.3 Å。因此，3 × 2的超晶胞对于容纳有机分子是足够的。

从菜单栏中选择Build | Symmetry | Supercell。将U更改为3，V更改为2。单击Create Supercell按钮，关闭对话框。

已经建立了表面结构。在继续下一步之前，保存工程。

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

4、表面上的分子的平衡位置计算

下一步是在表面上放置分子，并优化结构。在工程根目录文件夹中创建一个新文件，并将其用于计算。

使得优化好的表面结构为当前文档，从菜单栏中选择File | Save As...，打开Save As对话框。导航至工程根目录文件夹，将文件名更改为ChloroTiO2_110，单击Save按钮。

将在工程根目录文件夹中保存表面的一个副本。现在需要将有机分子粘贴到表面上。

选择优化后的chlorophenol.xsd文件。选择3D Viewer Selection Mode工具，选择距离测量，按下DELETE键。按下CTRL + C键，将当前文件重新改回ChloroTiO2_110.xsd文件。按下CTRL + V键。

将氯苯酚粘贴在包含表面的文件中。现在将重新调整分子取向，使其使垂直于表面。

按住SHIFT键和ALT键，单击鼠标右键，然后定位氯苯酚，使其悬停在表面上。

提示：如果使用带滚轮或鼠标中键的鼠标，可以通过按住SHIFT键和鼠标中键或滚轮来平移片段。

初始起始结构与下图类似。

注意：注意确保氯苯酚分子不会穿透表面。在尝试优化体系之前，请使用CPK显示样式验证这一点。

当研究分子与表面的不同相互作用时，固定表面原子。

使用3D Viewer Selection Mode工具选择TiO2表面中的所有原子。打开Edit Constraints对话框，清除并重新选择Fix Cartesian position复选框。关闭对话框。

现在优化结构。

单击Forcite Calculation对话框上的Run按钮。

这将创建一个包含优化结构的新结果文件夹。分子在表面上的取向大致是平的。

这是结构的一个能量极小值构型，但这只是对最小能量吸附位置的猜测。要计算全局能量最小值，可以尝试大约20种不同的起始构型，以确保已包含尽可能多的相空间。

或者，可以使用淬火分子动力学方法对许多不同的构型进行采样。继续计算之前，请保存项目并关闭所有文档。

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

5、使用淬火动力学对构型进行采样

将运行分子动力学模拟。

重新打开优化过的ChloroTiO2_110.xsd文件。

在Forcite Calculation对话框中，将Task更改为Quench，并单击More...按钮。

将打开Forcite Quench Dynamics对话框。

将可以访问淬火的基本选项。但是，可以通过其他动力学选项指定温度等。

单击Dynamics options后面的More...按钮。

打开Forcite Dynamics对话框。

注意：Frame output every XXX steps选项取决于Forcite Quench Dynamics对话框上的Quench every选项。

设置Temperature为350 K。

这个温度平衡了动能过大的体系（分子从表面解吸）和动能不足的体系（分子在表面移动）。可能需要在不同的温度下执行多个计算以确定正确的温度。对于本教程，350 K是一个合适的温度。

使用默认的5000步数以350 K运行动力学，进行5 ps的模拟时间。目前，每5000步执行一次淬火或几何优化。对于此模拟，它将只执行一次淬火。因此，将淬火次数值更改为每250步一次。

关闭Forcite Dynamics对话框。在Forcite Quench Dynamics对话框中，将Quench every更改为250步，然后关闭该对话框。在Forcite Calculation对话框中，单击Run按钮并关闭该对话框。

将在Project Explorer中创建另一个新文件夹ChloroTiO2_110 Forcite Quench，其中包含实时更新的结果文件以及以下文件：

ChloroTiO2_110.xsd：输入结构

ChloroTiO2_110 Quench Energy.xcd：包含淬火结构能量的图表

ChloroTiO2_110 Quench.std：包含淬火动力学结果的研究表，包括能量的分项

ChloroTiO2_110 Quench.xtd：包含淬火结构的轨迹文件

ChloroTiO2_110.xtd：包含原始非淬火结构的轨迹文件

ChloroTiO2_110.txt：包含初始设置和结果总结的文本文档

首先，研究初始的非淬火轨迹。

在Project Explorer中，双击以打开ChloroTiO2_110.xtd。

.xtd文件是由动力学计算结果的帧构成轨迹文件。可以使用可以使用动画工具栏查看轨迹。

选择View | Toolbars | Animation，显示Animation工具条，单击Play按钮

播放轨迹文件，可以看到氯苯酚分子在表面上的运动。氯苯酚分子可能会离开表面。这是2D周期性边界条件导致的镜像。可以更改可视化以删除此镜像。

单击Stop按钮

单击Animation Mode按钮的下拉箭头，从下拉菜单中选择Options，打开Animation Options对话框。勾选Recalculate atom visibility everyframe复选框，关闭对话框。单击Play按钮。

氯苯酚分子将在表面附近移动。

单击Stop按钮结束轨迹播放。

接下来将查看淬火后的轨迹。

将ChloroTiO2_110 Quench.xtd打开为当前文档，单击Play按钮。

将产生若干低能量构型，但是其中哪个是能量最低的构型？可以通过查看Forcite文本输出文件或使用数据表确定。

停止播放轨迹，并使得ChloroTiO2_110 Quench.std数据表文件为当前文档。

数据表的第一列为结构，其余列为轨迹和能量分项。总能量列于D列中，标记为哈密顿量。

选择列D, Hamiltonian。单击Quick Plot按钮

将绘制总能量与行号之间的关系图。

单击Chart Viewer Selection Mode工具。在数据表可见的情况下，单击图表中的一个点。

在图表中选择一个点将高亮显示数据表中的相应行。可以用这种方法来识别低能构型。或者，可以使用数据表上的排序工具。

使得数据表为当前文档，再次选择Hamiltonian列。单击Sort Ascending按钮

这将对数据表中的行进行排序，顶端的是能量最低的行。可以通过双击A列中的结构进行查看。

双击以打开A列中的能量最低的帧。

还可以使用3D原子收集文档将所有结构重叠显示，以查看表面上是否有优先成键位置。

使数据表文件为当前文档。选择结构列A，将鼠标移动到此列中的任何单元格上，单击鼠标右键，然后从快捷菜单中选择Extract To Collection。在警告对话框中单击OK。

默认视图将显示所有分子的位置。

打开Display Style对话框，在Lattice选项卡上将Style更改为In-Cell。

这将2D表面上的所有分子重叠显示，因此可以看到是否存在优先方向或吸附位置。存在有一种有效成键模式。

6、计算结合能

可以使用以下公式计算分子和表面之间的结合能：

Binding Energy = Etotal - (Echlorophenol + ETiO2Surface)

如果将这些能量中的每一个取为最小值（换句话说，使用几何优化构型的能量），则可以在低温极限下获得一致的结合能结果。

Etotal是氯苯酚-TiO2体系的能量，从淬火计算中得到。

由于已经把数据表按能量增加的顺序排列，可以从第一行读取哈密顿量的值，大约是-23.2 kcal/mol。

同时从初始几何优化结果中，可获得氯苯酚的能量。打开chlorophenol Forcite GeomOpt/chlorophenol.txt，找到最后报告的总能量Total Energy值，大约是5.85 kcal/mol。

注意：由于表面中的原子受到约束，它们之间的力的能量在整个模拟过程中保持不变。由于这些恒定力不影响吸附质的整体运动，计算中不包括这些恒定力，也并未对其进行报告。因此，对于固定的TiO2表面，报告的能量为零。

由于TiO2表面的能量为零，可以得出结论，TiO2上氯苯酚的零温度结合能为-29.05 kcal/mol。要在有限温度下获得结合能，需要孤立氯苯酚分子和TiO2表面上的氯苯酚进行分子动力学计算。然后使用每次计算得到的平均能量，包括动能。

使用之前的淬火分子动力学方法，已经找到了该结构和表面的优先结合位点，并计算了结合能。现在可以进一步功能化这个有机分子，看看是否有不同的相互作用，或者可以改变表面。如果需要，也可以尝试运行更长时间来采样更多构型。