Chem3D中的分子力场优化

Posted by ZH DENG on September 13, 2016

分子力场简介

在研究中,我们时常需要进行分子模拟;而分子模拟的基础是准确计算原子之间的相互作用:包括组成同一分子的原子之间的成键相互作用,和不同分子间的范德华相互作用;有的分子间还有氢键相互作用。描述原子间的这些相互作用有两种方式,一个是通过量子化学计算,另外一种方式就是采用分子力场计算。

虽然量子化学计算分子结构和原子、分子间相互作用比较准确,但是需要消耗大量的硬件资源和时间;而分子力场计算则很快,因为分子力场并不计算电子相互作用——它是对分子结构进行简化得到的一种模型。在这个模型中,分子力场使用牛顿经典力学来处理原子之间的相互作用:把组成分子的原子看成是由弹簧连接起来的球,然后用简单的数学函数来描述球与球之间的相互作用。如将氢分子看做有弹簧连接的两个球,我们可以用简单的胡克定律描述两个氢原子间的能量:E=k(b-b0)2。其中,b表示两氢原子间距离,b0表示平衡时原子间距,k为键能系数,b0和k称为力场参数。更复杂一点可以用四次Taylor展开式表达:E=k1(b-b0)2+K2(b-b0)3+K3(b-b0)4,更多的参数可以获得对成键分子的更精确的描述。

但是,分子力学对分子结构和原子间相互作用描述的准确性这依赖于所用的参数。而这些参数通常拟合自实验数据,或者量子化学结果。它属于经验描述,显然质量要低一些,但是由于计算速度快,适合于描述上千个乃至百万个原子的模拟,而在这些情况下,我们无法采用量子力学计算,因此,只能采用分子力场进行模拟。

简单的势能面扫描

在Chem3D中,最为常用的分子力场就是MM2,它适合对一些小分子进行模拟。要学习MM2力场的使用,我们首先要理解势能面与构象之间的关系。我们知道,一个分子的能量很大程度上取决于构象,而我们可以将能量与构象联系起来,绘制出能量随构象参数变化的图形——势能面。有了势能面,我们就可以对分子模型直观的分析。

势能面使得我们能够直观地研究贩子能量随构象的变化。Chem3D提供了基于分子力场的势能面扫描。相应的操作在 Calculations → Dihedral Driver 内。

Dihedral Driver最多能够绘制关于两个二面角参数的势能面,对于每个二面角参数,它采用MM2力场进行能量优化,最后得到势能面图像。

例如:1)二维势能曲线 绘制乙烷的势能面

选择乙烷的碳碳键,点击 Calculations → Dihedral Driver → Single Angle Plot 。然后计算开始。等到计算完成后,我们会得到一条二维曲线,它描述了乙烷分子能量随H-C-C-H二面角变化的情况。

dgs

dgs

用鼠标选择在曲线上的某一点,就可以查看对应的二面角角度及能量,Model Window中也会显示相应的构象。用鼠标在曲线上连续拖动,就可以在Model Window观察构象的连续改变过程。

dgs

2)三维势能面 绘制(R)-3-methylbutan-2-amine的势能面

选择相邻的碳碳键和碳氮键,点击 Calculations→Dihedral Driver→Double Angle Plot。

dgs

dgs

计算完成后,我们会得到三维势能面,只是能量高低换做颜色变化来描述了。默认状况下为黑白(若要调整颜色,在显示框内右键选择即可),从右侧的Legend bar 我们能知道,白色代表能量高,黑色代表能量低。查看构象及对应的能量,需要用鼠标选择图像中的一个点。连续拖动就可以在Model Window观察构象的连续改变过程。a为能量最高的构象,b为能量最低的构象。

dgs

dgs

dgs

但是,对于二面角的简单旋转并未考虑其他可能受到旋转操作影响的构象因素,这使得我们绘制出的势能面可能不准确。这时,就需要进行二次计算,来保证计算结果的准确性:在Dihedral Driver Chart内任意部分右键,点击Recompute with Minimization即可。当得到满意的图表以后,在Dihedral Driver Chart内任意部分右键,可选择Copy data,Copy Picture,Save Picture等多种导出或保存方式。

MM2分子动力学简介

MM2 分子动力学采用牛顿力学来模拟分子中的原子的热运动,随着温度变化来调整原子的动能。

dgs

在工具栏中,MM2 Dynamics是MM2下面有绿色箭头的按钮,在它的右边从左到右依次是:

Step Interval :指定两次模拟之间的间隔,以fs为单位。通常1-2fs就能给出较好的结果。

Frame Interval :指定帧和统计数据收集的时间间隔,也以fs为单位。通常10-20fs能给出了相当平滑的帧序列。

Heating/Cooling Rate:指定模型的动能加减速率。

Target Temperature :动力学计算的目标温度。

这里有一个极端的例子,氯化钠的分解:

实例操作

D-A反应的选择性

我们都知道,D-A反应会得到内型和外型产物;通常情况下后者的热稳定性更高。接下来,我们将试着利用Chem3D来计算两种过渡态的能量并加以比较。

1.在ChemDraw panel中分别绘制丁烯二酸酐和环戊二烯,然后再用单键工具连接需要成键的两对碳原子。注意,当用单键连接以后,程序会自动应用碳四价规则,这需要我们手动编辑这四个碳原子的结构以恢复到超价状态;

dgs

2.首先画出外型过渡态。我们关心的是新形成的六元环过渡态,所以需要手动设置这六根键的参数。由于过渡态有对称面,将这六根键分为3组。首先是C(1)-C(2)键和C(6)-C(7)键,这两根键在过渡态中的键级处于1-2之间,我们将它们的键长都调整为1.4埃左右(碳碳单键键长约为1.3埃,碳碳双键键长约为1.5埃)。

dgs

3.接下来是C(1)-C(5)键和C(2)-C(3)键,这两根键在过渡态中的键级也处于1-2之间,也调整键长为约1.4埃。

dgs

4.最后是新形成的C(3)-C(7)键和C(5)-C(6)键,它们的键级在0-1之间,键长调整为约2.0埃。

dgs

5.点击MM2 Dynamics,然后再使用MM2 Minimize,得出计算结果和优化后的构型。

6.那么如何得到内型过渡态呢?重画一遍?No。选择原来丁烯二酸酐中的一个桥头碳原子,点击菜单栏Structure→Invert,翻转其构型,另一个碳原子如法炮制,就得到内型过渡态。接下来的步骤如前所述,得到计算结果。

dgs

dgs

7.计算结果是内型过渡态能量为535.3kcal/mol,而外型过渡态能量为493.2kcal/mol,外型过渡态的热力学稳定性更高。由于使用的是牛顿力学体系,计算结果的可靠性并不高。但是可以预见的是,内型过渡态空间上比外型过渡态更为拥挤,能量也应该更高。在本例中,计算结果与我们的理论判断相符合。