实验二：分子与表面的对接优化及分子动力学模拟 - 副本

计算材料学》实验讲义实验二：分子与表面的对接、优化及分子动力学模拟一、前言分子力学优化分子力学方法又称力场方法，其基于非常简单的经典力学模型，忽略了电子运动，把体系能量看作是原子核坐标的函数，其贡献来自诸如键伸缩、单键键角的张合以及旋转等等。该方法从本质上说是能量最小值方法，即在原子间相互作用势的作用下，通过改变粒子分布的几何位型，以能量最小为判据，从而获得体系的最佳结构。因此，确定分子间的相互作用势是进行分子力学优化的关键，在分子力学中用力场来描述分子中各原子间的相互作用。所谓力场是指描述各种形式的相互作用对分子能量影响的函数，其有关参数、常数和表达式通常称为力场。一般力场的表达式为E—E+E+E+E+E+（2-6）stretchbendtorsionvdwelec式中，E.为键的伸缩能，E.为键的弯曲能，二者均采用谐振子模型；E.为键的stretchbendtorsion扭曲势，它采用傅立叶级数形式来描述；E、E为非键作用项，分别表示范德华相互作用vdwelec和静电相互作用。分子模拟所使用的力场，从最初的单原子体系扩展到多原子分子、聚合物、生化分子体系。力场也从简单的非键相互作用，扩展到复杂的形式。每个力场针对特殊目的有所侧重，各有优缺点和使用范围。在模拟中计算中选择合适的力场尤为重要，也是决定计算结果成败的关键。对于全原子模拟而言，人们越来越重视力场的发展，概括的讲，可以把力场的发展趋势归结为三点：第一是朝着通用的方式发展，几乎覆盖所有的原子类型；第二是重点强调和提高特定应用范围内的性质预测；第三是在适当的研究范围内追求结果的精确性，预测的性质包括分子结构、构型性质、振动频率、生成热等。目前常见的力场主要包括AMBER（针对蛋白质、核酸等生化分子）、OPLS（针对多肽核酸和有机溶剂的液体体系）、CHARMM（针对有机分子、溶液、聚合物、蛋白质等）、Tripos（有机小分子及生物大分子）、YETI（含金属的生物小分子）、Universal（主族元素化合物，有机分子，金属配合物）、Dreiding（主族元素小分子）、GROMOS（生物大分子）、UFF（—般分子）、CVFF（—般分子）、CFF（有机小分子）和COMPASS（有机和无机分子）。2.COMPASS力场COMPASS（Condensed-phaseOptimizedMolecularPotentialsforAtomisticSimulationStudies）是第一个能够在很大范围内，准确模拟和预测单个分子或者凝聚态物质的结构、构象、振动频率、热力学性质的从头算力场，之所以称之为从头算力场，是因为这个力场中多数参数都是通过量子力学计算得到的。COMPASS力场可研究的对象包括最基本的有机、无机小分子、聚合物等，还可以模拟计算由金属离子、金属、金属氧化物等形成的新型材料虽然不同的体系需要采用不同的模型，但COMPASS力场的参数是不变的，理论上讲。可以研究包括不同界面和材料的复杂体系。对于所有的无机和有机共价分子（包括聚合物），COMPASS力场中的总是能写作成键形式、非成键形式能量的总和。对于成键形式，势能表达式如下：E=E[k(b-b)2+k(b-b)3+k(b-b)4]+工[k(9-9)2+k(9-9)3+k(9-9)4]pot203040203040TOC\o"1-5"\h\z+E[k(1-cos©)+k(1-cos2©)+k(1-cos3©)]+》k咒2+Ek(b-b)(b'-b')+Ek(b-b)(9-9)12320000Xb,b'b,9(b-b)[kcos©+kcos2©+kcos3©]+乙(9-9)[kcos©+kcos2©+kcos3©]01230123+乙k(9-9)(9'-9')+乙k(9-9)(9'-9')cos©0000b，99,9,©对于非键形式，包括描述经典相互作用的库仑项和范德华相互作用的Lennard-Jones势，经典qqr0r0相互作用表示为：E=丫空j;范德华相互作用表示为：E=丫E[2（^）9-3（#”。elervdWLJrr与其它力场相比，COMPASS力场采用了更复杂的一套函数形式对相互作用势进行了更加准确的描述。自1998年力场产生以来，已经验证和预测了大量的科学实验数据和材料的结构性质。该力场的摘要编写者为中国学者孙准教授。3扮子动力学模拟分子动力学模拟方法的基本思想是把物质看成由原子和分子组成的粒子系统（many-bodysystems），从该体系的某一假定的位能模型出发，并假定体系粒子的运动遵循经典力学或量子力学描述的规律，若已知粒子的所有受力作用，则可以求解出运动方程而得到系统中全体粒子在相空间中的轨道，然后统计得到系统的热力学参数、结构和输运特性等。MD方法的基本原理是：在各态历经假设（ErgodicHypothesis）下，任何力学量的系综平均等效于该力学量对时间的平均，而力学量时间的平均可以从经典运动方程所决定的运动轨迹得到。具体地说，对一个由N个原子构成的简单系统其势能项可表示为U二U+UVDWint其中U为分子间和分子内的非键范德华作用，U为分子内部势能（包括键伸缩、键角弯VDWint曲、键扭转等）。依据经典力学，系统中任意原子i所受力可表示为势能的梯度F二—vu二—ii（dF二—vu二—ii—i+——j+——kU.Ox5y5z丿iii假定体系从某一势能模型出发，给每个原子赋予一定的初速度，借助牛顿运动方程aC）=巴imi可求得体系中各原子的加速度，进而将牛顿运动方程对时间积分，可求得各原子的速度和位d2-d--r=v=adt2idtiiv=v0+at）iii1r=ro+vot+a12iii2i式中，尸和V分别为原子的位置和速度，上标“0”为各物理量的初始值。重复以上步骤，可得到一个时间序列，从而获取各时间下系统中分子运动的位置、速度和加速度等信息；在各态历经假设的基础上，即可得任一物理量Q的系综平均，从而获取体系的有关性质:Q；：=QC)=limlj次C)dtst0因此，从计算数学的角度看，MD方法是微分方程的数值解问题。二、实验目的（1）理解分子力学方法及分子动力学方法的原理。（2）掌握分子与固体表面对接的方法。（3）掌握分子与固体表面体系优化及分子动力学模拟的基本步骤。三、实验内容

以甲基乙酸甲酯与Fe表面对接过程为例，展开详述。优化分子结构。构建分子。首先构建二甲基乙酸甲酯分子，并命名为molecule。分配力场。利用分子力学方法对分子结构进行优化，首先需要给构建的分子分配力场，选择discover模块〒中的setup，出现如下对话框：选择Energy，确定Forcefield选择为compass；选择Typing，点击Calculate，给定分子每个原子力场参数，即分子内及分子间的相互作用势；点击Automation，确认相关选项均选择No；选择Non-Bond，Applysettingsto选择vdW&Coulomb，Summationmethod选择Atombased，Quality选择Medium，点击More，cutoffdistance选择9.5A；点击JobControl选择MyComputero(3)优化分子。选择discover模块T中的Minimizer，出现如下对话框：丽DiscoverIinuizationMethod:|SmartMinirniz已「二]Convergencelevel:|MediumMaximumiterations:厂丽DiscoverIinuizationMethod:|SmartMinirniz已「二]Convergencelevel:|MediumMaximumiterations:厂Optimizecell:Minimize]Files...Help5000Method选择SmartMinimizer,Convergencelevel选择Medium,点击Minimize,开始进行优化。优化完毕后，在Project中会出现文件名为moleculeDiscoMin的文件夹，文件夹内名为molecule.xsd的文件即为优化后的分子结构。如下图所示：构造表面并优化表面。（1）构造表面。首先导入Fe晶体结构，选择菜单栏Build一Surfaces一CleaveSurface按钮，出现如下对话框：将Cleaveplane（hk1）从（-100）改为（001）,并按下Tab键确认。将Thickness从1.0改为3.5，并按下Tab键确认(注意：表面的厚度要大于cut-offdistance的值)。点击cleave，在Project中出现Fe(001).xsd的文件，即为Fe晶体的(001)面。如下图所示。RT1EditConstraintsKelpRT1EditConstraintsKelp(2)分配力场。利用分子力学方法对Fe(001)面进行优化，首先需要给构建的分子分配力场，选择discover模块中的setup，具体参数设置与优化分子结构时的一致，此处不再重复。优化表面。通过鼠标拖拽选中表面两层以下的所有原子，选择菜单栏Modify—Constraints，在Editconstraints对话框中，选中FixCartesianPosition和选择FixFractionalPosition，并关闭对话框。如下图所示:Atom|Measurement|LatticepFixCartesianposition17X17Y|7Z層Fixfractionalposition选择discover模块丁中的Minimizer，具体参数设置同分子结构优化参数设置。最后点击Minimize，完成优化。(4)构造超晶胞。使优化好的Fe(001)处于活动窗口，选择Build一Symmetry一Supercell，将U、V均改为8(注意：构造的超晶胞的边长应该大于2倍的cutoffdistance)，如下图所

示：分子与表面对接及优化。(1)加入真空层。加入真空层后，表面变为一个3D周期单元，底部是表面而其上保持真空状态。它允许你在真空层中将分子放在表面上。选择菜单栏Build一Crystal一BuildVacuumSlab，出现如下对话框：将VacuumThickness改为30A,SlabPosition改为0.05A,点Build完毕。得到如下3D周期性结构。

(2)将分子对接到表面。首先打开优化好的甲基乙酸甲酯分子的3D文档，双击分子结构以选中分子，点击鼠标右键，选择copy；然后打开刚构造好的Fe(001)面3D结构，点击鼠标右键，选择paste，就将分子拷贝到了表面上。按下Shift+Alt键，同时点击鼠标右键，可拖动分子，将分子定位到表面上。如下图所示：'•!*B#.r栓席『'•!*B#.r栓席『[.一1\.4(4)优化该结构。选中所有Fe原子，选择菜单栏Modify—Constraints，在Editconstraints对话框中，选中FixCartesianPosition和选择FixFractionalPosition，并关闭对话框。选择discover模块"中的Minimizer,Method选择SmartMinimizer,Convergencelevel选择Medium,点击Minimize，开始进行优化。优化完毕后，得到如下结构。4扮子与固体表面相互作用的分子动力学模拟(1)打开上述优化好的3D结构，从工具栏中选择，然后选择Dynamics，或者从菜单上选择Modules-Discover-Dynamics，出现如下对话框：

将dynamicstime改为50ps,Save改为Full,Frameoutputevery改为500,点击Run.计算需要花费几分钟的时间，运行完毕后，出现名为Fe(001)DiscoDynamics的文件夹。5•计算分子与固体表面的相互作用能分子与固体表面的相互作用能可通过以下公式计算：E二E-(E+E)interacitontotalsurfacemolecule其中，Ettl为分子与表面的总能量，Ef是除去分子后表面的能量。E是除去表面后totalsurfacemolecule分子的能量。打开动力学模拟后的3D文件，选中体系所有原子，并且在菜单上选择Modify-Constraints，取消原子固定。首先计算整个体系的能量，选择discover中的setup，在Energy标签下，按下Calculate按钮。计算结束后，在输出文档中查找TotalPotentialEnergy(-47519.078738kcal/mol)。(3)计算固