计算机在材料科学与工程中的应用-分子动力学.ppt

《计算机在材料科学与工程中的应用-分子动力学.ppt》由会员分享，可在线阅读，更多相关《计算机在材料科学与工程中的应用-分子动力学.ppt（69页珍藏版）》请在三一办公上搜索。

1、计算机在材料科学与工程中的应用计算材料学(分子动力学)郝 艳,分子动力学,1.介绍分子动力学的基本原理及步骤 2.分子动力学在Materials Studio（MS）中的实现 重点讲解Discover、Forcite的应用实例,原则上,第一原理方法在理论上已经能解决所有问题，但计算量太大，计算机资源有限，原子数目较多时，如高分子、蛋白质、原子簇以及研究表面问题、功能材料或材料的力学性能等，实际上难以完成计算。为此，发展了分子力学(Molecular Mechanics,MM)与分子动力学(Molecular Dynamics,MD)方法，MM与MD是经典力学方法，针对的最小结构单元不再是电子而

2、是原子。因原子的质量比电子大很多，量子效应不明显，可近似用经典力学方法处理。,20 世纪 30 年代,Andrews 最早提出分子力学的基本思想；40 年代以后得到发展,并用于有机小分子研究。90年代以来,随着计算机技术和计算算法的发展而得到迅猛发展和广泛应用。在材料科学,用于材料的优化设计、结构与力学性能、热加工性能预报、界面相互作用、纳米材料结构与性能研究等；在化学领域，用于表面催化与催化机理、溶剂效应、原子簇的结构与性质研究等；在生物科学和药物设计领域的应用也十分普及，如蛋白质的多级结构与性质，病毒、药物作用机理、特效药物的大通量筛选与快速开发等。,MM和MD的应用，通常称作分子模拟(m

3、olecular simulation,molecular modeling)或 分子设计(molecular design)。MM是确定分子结构的方法，是分子的静态势函数，利用分子势能随结构的变化而变化的性质，确定分子势能极小时的平衡结构(stationary point)；而实际过程通常是在一定温度和一定压力下发生的，为了更切实际地了解体系运动和演化的过程，必须考虑体系中原子的运动，并与温度T和时间t建立联系；MD含温度与时间,还可得到如材料的玻璃化转变温度、热容、晶体结晶过程、输送过程、膨胀过程、动态弛豫(relax)以及体系在外场作用下的变化过程等。,分子动力学模拟应用实例,晶体结构确

4、定及性能预测饱和脂肪酸大分子晶体结构确定及性质预测（西班牙巴赛罗那大学）,晶体中的分子以强烈的氢键形成二聚体而整齐排布,分子模拟方法与粉末衍射实验分析方法结合分析复杂分子的晶体结构行之有效,分子动力学模拟应用实例,晶体形貌研究吸附剂分子对Al2O3晶体各个面生长速率的影响（加拿大Alberta大学）,吸附能的排列次序与实验观察到的各个面的生长速率倒数成正比,两个具有最低吸附能的晶面在晶体生长过程中其主导作用，并最终决定晶体的宏观形貌,分子模拟方法可实现晶体生长形貌预测及控制,什么是分子动力学,经典分子动力学将原子视为经典粒子，通过求解各粒子的运动方程得到不同时刻粒子的空间位置、运动状态，从而统

5、计出材料的宏观行为特性。,MD是用计算机方法来表示统计力学，用来研究不能用解析方法来解决的复合体系的平衡和力学性质，是理论与实验的桥梁。,分子动力学的基本思想,经典力学定律 分子动力学中处理的体系的粒子遵从牛顿方程，即 式中 是粒子所受的力，为粒子的质量，是原子i的加速度，原子i所受的力 可以直接用势能函数对坐标 的一阶导数，即，其中U为势能函数（简称势函数或力场），因此对N个粒子体系的每个粒子有,求解这组方程要通过数值方法，即给出体系中每个粒子的初始坐标和速度，从而产生一系列的位置与速度，即为任意时刻的坐标和速度。,分子动力学方法工作框图,势函数或力场,轨迹:分子动力学整个过程中的坐标和速度

6、称为轨迹。,经典分子动力学的适用范围,分子动力学方法只考虑多体系统中原子核的运动，而电子的运动不予考虑，量子效应忽略。,1,300K时，KBT=2.5J/mol临界频率=6.251012 S-1,因此，经典分子动力学不适用：相对的高频率的运动；涉及电荷重新分布的化学反应、键的形成与断裂、解离、极化以及金属离子的化学键等等。,分子动力学模拟实施步骤,分子动力学运行流程图,时间步长,参考原子或分子特征运动频率来选取,在应用软件的实际操作中，需要设置的参数还很多，包括温度、压力、力场、算法等等。,第一性原理分子动力学 用分子动力学方法讨论材料的结构、相变及力学性质，已经被广泛的研究，其势函数的选取有

7、很多种，诸如Lennard-Jones势分子动力学、Morse势分子动力学方法等。针对不同的材料，构建介观条件下的对势，取决于对材料介观结构的深刻理解，这给势函数的构建带来一定的困难，从而给经典分子动力学的模拟带来困难。随着密度泛函理论的发展，第一性原理方法可解决这些困难。1985年，Car和Parrinello 提出第一性原理分子动力学方法，第一性原理分子动力学方法的势场直接来源于电子结构的计算，而不是经验势。,重点模块涉及模块,Discover：是Materials Studio的分子力学计算引擎，适用于很大范围的分子和材料。Forcite：先进的分子力学和分子动力学模拟程序。,Mater

8、ials Visualizer：MS的核心模块，提供建模、分 析和可视化的工具。Amorphous Cell：构建复杂无定形模型并预测关键性 质，一般与Discover、Forcite连用。,分子动力学在Materials Studio（MS）中的实现,本课程所使用的软件Materials Studio 7.0,Discover是Materials Studio的分子力学计算引擎为Amorphous Cell、Forcite等模块提供了基础计算方法。,计算最低能量构象 给出不同系综下体系结构的动力学轨迹 得到各类结构参数 热力学性质 力学性质 动力学量 振动强度,Discover Setup菜

9、单Energy,Commpass:是第一个把有机分子体系与无机分子体系统一的分子立场。适合于共价分子体系，包括大多数常见有机物、无机物和聚合物、金属、金属氧化物和金属卤化物。Pcff：polymer consistent force field，基于CFF91力场发展而来，适用于聚合物及有机物。可用于聚碳酸酯类、多糖类等聚合物、无机和有机材料，包括约20种金属(Li,K,Cr,Mo,W,Fe,Na,Ni,Pd,Pt,Cu,Ag,Au,Al,Sn,Pb）、糖类、脂类和核酸，以及惰性气体(He,Ne,Kr,Xe)。Cvff：Consistent valence force field，一致性价力场

10、，最初以生化分子为主，经不断强化，可适用于计算多肽、蛋白质与大量的有机分子。,Discover SetupNon-bond,设置相互作用力：范德华力（vdW）库仑力（Coulomb）范德华和库仑力（vdW&Coulomb）,计算方法（Summation method）：基于原子的总量（Atom Based）适合于孤立体系，对于周期性体系计算量较小，但是准确性较差 基于电子群的总量（Group based）适合于周期性和非周期性体系，计算的准确性好一些，计算量最小 基于指定数量层（Cell Mutipole）不采用截断方法（No Cutoff）适合于周期性体系，计算最为准确，但计算量最大,精度（

11、Quality）：Corse（粗略）Medium（中等）Fine（精细）Ultra-fine（超精细）,Job Control菜单,Discover Minimization,优化方法（Method）：综合法(Smart Minimizer)：综合以下三种方法。最速下降法(Steepest Descent)：计算简单，速度快，但在极小值附近收敛性不够好，造成移动方向正交，适用于优化的最初阶段。共轭梯度法(Conjugate gradient)：收敛性较好，但对分子起始结构要求较高，因此常与最速下降法联合使用，先用最速下降法优化，再用共轭梯度法优化至收敛。牛顿法(Newton)：计算量较大，当微

12、商小时收敛快，也常与最速下降法联合使用。,Convergence level：收敛精度水平Maximum interation：最大迭代数Optimizer cell：选中的话表示优化晶胞参数和原子位置,Discover MolecularDynamics,Ensemble：系综Temperature：目标温度Pressure：给系统所施加的压力Number of steps：整个动力学所运行的总步数Time step：每一动力学步骤所花费的时间（单步长时间）Dynamics time：Number of stepsTime step（总模拟时间）Trajectory SaveCoordina

13、tes表示保存坐标Final Structure表示只保存最终结构Full表示保存所有。Frame output every：若输入5000，则表示每5000步输出一帧（即晶体结构）。,字母含义如下：N=固定粒子数V=固定体积E=固定能量T=固定温度P=固定压强H=固定焓,系综（Ensemble）,速率法（Velocity Scale）：系统温度和粒子的速度直接相关，可以通过调整粒子的速度使系统温度维持在目标值。Nose and Andersen：控制热力学温度并生成正确的统计系综的概率，热力学温度不变，允许模拟系统交换能量的“热浴”。Berendsen:包括之间的热能交换系统和热浴，允许指定

14、衰减常数的值。注：Nose法适用于自相关研究，速率法和Andersen法不适应于不连续轨迹和速度的自相关研究。,Setup菜单,计算任务（Task）：单点能计算 几何优化 动力学模拟 淬火模拟 退火模拟 剪切模拟 限制性剪切模拟 内聚能密度计算 力学特性计算精度控制（Qualigy）,几何优化的参数设置,ABNR法：改良后的Newton-Raphion法，常用于生物分子体系。,淬火模拟的参数设置,退火模拟的参数设置,内聚能密度计算的参数设置Setup菜单-Cohesive Energy Densit/More 计算分子内相互作用 输出Study Table文件 Study table中包括输入

15、的结构文件,力学特性计算的参数设置Setup菜单-Mechanical Properties/More 优化结构 明确应变模式中产生的应变数目 推荐使用偶数值（2-100）指出结构最大的形变 值在之间较为合理 应变模式（Strain Pattern）应变张量矩阵，由结构对称性决定,COMPASS,力场,适合于共价分子体系，包括,大多数常见有机物、无机物和聚合物、,金属、金属氧化物和金属卤化物,。,GroupA,共价模型，有机物、聚合物和气体分子,GroupB,离子模型，金属、金属氧化物、金属卤化物、沸石,（,0K,）,COMPASS26,增加二苯醚、苯腈、叠氮化物、醇类新参数，改进甲烷的参数,

16、COMPASS27,加入二硫键基团的参数,COMPASS,加入硫酸根、磺酸根基团的参数,COMPASS,力场,适合于共价分子体系，包括,大多数常见有机物、无机物和聚合物、,金属、金属氧化物和金属卤化物,。,GroupA,共价模型，有机物、聚合物和气体分子,GroupB,离子模型，金属、金属氧化物、金属卤化物、沸石,（,0K,）,COMPASS26,增加二苯醚、苯腈、叠氮化物、醇类新参数，改进甲烷的参数,COMPASS27,加入二硫键基团的参数,COMPASS,加入硫酸根、磺酸根基团的参数,COMPASS,力场,适合于共价分子体系，包括,大多数常见有机物、无机物和聚合物、,金属、金属氧化物和金属

17、卤化物,。,GroupA,共价模型，有机物、聚合物和气体分子,GroupB,离子模型，金属、金属氧化物、金属卤化物、沸石,（,0K,）,COMPASS26,增加二苯醚、苯腈、叠氮化物、醇类新参数，改进甲烷的参数,COMPASS27,加入二硫键基团的参数,COMPASS,加入硫酸根、磺酸根基团的参数,Universal:元素周期表的完整覆盖。适用于计算有机、主族无机分子和金属络合物的几何结构、构象能量差异。对于有机金属体系或其他力场不包含相关参数的体系，推荐使用该力场。Dreiding:基于杂化规则的力常数和几何结构参数。适用于计算有机、生物和主族无机分子的几何结构、构象能、分子间结合能和晶体堆

18、积。COMPASS26 增加二苯醚、苯腈、叠氮化物、醇类新参数，改进甲烷的参数COMPASS27 加入二硫键基团的参数pcff30：pcff力场会有所更新，最新的版本始终命名为pcff。于此同时，较早的一个版本也会保留下来，用于验证和已有工作的继续，例如pcff30。cvff_nocross_nomorse:当体系能量较高时，Morse函数会允许成键原子分离至不合理的距离。当体系结构远离平衡时，交叉项可能是不稳定的。,气体在聚合体中扩散的测量,目的：介绍如何使用力场方法来计算气体在材料中的扩散系数。模块：Materials Visualizer,Discover,COMPASS,Amorpho

19、us Cell,背景 气体在有机溶剂，聚合体或沸石中的扩散率可以通过分子动力学模拟来计算，同时也可以计算气体在材料中的均方位移。这可以让你计算气体的自扩散系数，并进而可以研究全扩散系数。当进行分子动力学计算的时候，可以分析温度，压力，密度，渗透尺度和结构对扩散的影响。简介 在本教程中，我们将通过构建一个包括氧和二甲基硅氧烷（PDMS）的无定形晶胞中计算氧气在该聚合物。当构建了晶胞以后，将进行分子动力学模拟并计算氧分子的均方位移。虽然本教程中的时间尺度限制了计算，但还是可以用来熟悉相关的方法。本教程基于Charati 和Stern(1998)年发表的一篇研究气体在硅聚合物中扩散的文章。,选择此文

20、件夹存放数据,生成了名称为diffusivity的Project,写入diffusivity,这样就产生了新的 Materials Studio project，开始了Materials Studio 运行,.新建一个Materials Studio工作任务,2.建立初始结构 第一步是构建并优化氧分子和PDMS 聚合物来构建无定形原胞。从菜单栏中选择Build/Build Polymers/Homopolymer 来显示Homopolymer 对话框。,把库Library改成硅氧烷siloxanes，把重复单元Repeat unit改成二甲基硅化物dimeth_siloxane。,在Homop

21、olymer 对话框中选取Advanced。选上Random，点击Build。关闭Homopolymer 对话框。一个名为Polydimeth_siloxane.xsd 的新的3D 自动文档会打开。,在Project Explorer 中，右键点击project root 并选择新的3D Atomistic Document。右键点击3D Atomistic.xsd 并选择重命名。把名字改成Oxygen 并点击回车。,现在可以勾画出氧分子。激活oxygen.xsd。点击Sketch Atom 按钮，从下拉菜单中选择oxygen。在3D Viewer上左键单击，然后松开左键，移动鼠标以形成一根键

22、。鼠标移到一定距离，键不能再伸长。双击左键，完成构建。把鼠标移到键上面，它会变成浅蓝色，这时左键点击一下变为双键，O2分子完成构建。注意，在这些操作中，鼠标状态为。不能点。完成O2分子的构建后，点，避免产生新的原子。,你需要对氧分子命名一下，不然，MS Modeling 就会用默认的名字。在Properties Explorer 中，把Filter 改成Molecule。双击Name，输入oxygen，点击OK。注意核对 ChemicalFormula中是否显示O2。,一个经验力场计算（能量最小化或分子动力学）中花费最大的部分是非键参数的确定（库仑相互作用和范德华力）。涉及力场的计算会用各种方

23、法来计算非键参数，随所研究系统的尺度和类型而变化。不过对范德华力默认的方法是原子级模拟，对库仑相互作用则是Ewald加和模拟。对某些聚合物，可以用一组原子而不是单个原子来逼近非键参数。这种方法叫作charge groups。本教程中你会从头到尾用到这个方法。这种方法可以在不损害精度的情况下加速计算。现在聚合体将自动用charge groups 来计算，如果要显示的话，点击Display Style 对话框。激活Polydimeth_siloxane.xsd 文档。右键点击3D 原子文档，选取DisplayStyle。在Display Style 对话框中，把Color by 选项改成Charg

24、e Group。,在Charges 对话框中指明氧分子是用charge group 的。激活oxygen.xsd。从菜单栏中选取Modify/Charges 来显示Charges 对话框，选择Charge Groups条目，点击Calculate。,在优化两个分子的几何结构之前，必须要让Discover 知道用charge goups 来进行非键计算，而不是用默认选项。在Job Control中选My Computer。,现在可以开始优化两个几何结构了。点击工具条上的Discover 按钮，然后从下拉列表中选择Minimizer。激活oxygen.xsd。点击Discover Minimiza

25、tion 对话框中的Minimize 按钮。,现在任务浏览器显示出来了，并且在Project Explorer 中创建了一个新目录oxygen Disco Min。当计算完成时，最小化的结构会被存放到这个新目录下。,激活Polydimeth_siloxane.xsd，点击Minimize 按钮。计算结束后最小化的结果被返回到Polydimeth_siloxane Disco Min/Polydimeth_siloxane.xsd 中。关闭Discover Minimization 对话框。,现在有了两个优化的几何结构。,在File中点击Save Project。从菜单栏中选择Windows|C

26、lose All。在Project Explorer 中打开最小化的结构oxygen Disco Min/oxygen and Polydimeth_siloxane Disco Min/Polydimeth_siloxane.xsd。,3 建一个无定形的晶胞当你建好两个结构后，就可以用Amorphous Cell 模块来把它们往一个晶胞中成倍地复制。在工具栏上选择Amorphous Cell 按钮，然后从下拉列表中选择Construction。将会显示Amorphous Cell 对话框。,第一步是指明组成晶胞的分子。激活oxygen.xsd，点击Add 按钮。对Polydimeth_sil

27、oxane.xsd 重复同样操作。,氧分子和PDMS 各十个被添加到晶胞中去。不过，我们想建的是包含个氧分子和八聚PDMS的晶胞。,在Constituent molecules 部分，点击Number cell for oxygen，把它改为4。对Polydimeth_siloxane.xsd 作同样操作，不过把数值改为8。把Number of configurations 从10 改为1，把Target density of the final configurations 从1 改为0.95。不选上the Refine configurations following construct

28、复选框。,在Amorphous Cell Construction 对话框中选择Setup 条目。在Job Control 部分，不选上Automatic 并在文本区域输入cell，点击Construct。当Amorphous Cell 构建了一个结构后，默认是把这个结构与组成分子列表中的第一个分子取相同的名字。本例中，你要把它改成cell。,在Project Explorer 中出现了一个新的名为 AC Constr 的文件夹。当计算结束时，会产生一个包含不规则晶胞的轨迹文档cell.xtd。关闭Amorphous Cell Construction 对话框。双击cell.xtd。这个文档中

29、包含了一个有八聚PDMS 和4 个氧分子的周期性晶胞。,3.晶胞的弛豫 当一个无规则晶胞生成时，分子可能不是等价地分布在晶胞中，这样就造成了真空区。为了矫正这个，要进行能量最小化来优化晶胞。最小化过后，要进行分子动力学模拟来平衡晶胞。在能量最小化之前，清空工作区。选择File|Save Project，接着再从菜单栏中选取Windows|Close All。双击 Project Explorer 中的cell.xtd。当一个包含周期性结构的3D 原子文档被打开时，那些非键的设定会重新变成默认值。文档cell.xtd 中也有周期性结构，因此在打开之后要把非键的设定从默认值改回来。从菜单栏中选择M

30、odules|Discover|Setup 来显示Discover Setup 对话框，从中选取Non-Bond条目。把Apply settings to 改成 vdW&Coulomb。把Summation method 改成 Group Based。,在构建无规则晶胞时，都要用能量最小化和分子动力学来进行结构弛豫。,关闭Discover Setup 对话框。现在你已经准备好对整个晶胞进行能量最小化了。由于本教程中时间有限，只能进行2000步的优化计算。在实际计算中，因该把整个优化运行完全。点击工具条上的Discover 按钮，然后从下拉列表中选择Minimizer。在Discover Min

31、imization对话框中，把Maximum iterations从5000改为2000。点击Minimize。关闭Discover Minimization 对话框。,任务结束后，最终的结构保存在文件夹cell Disco Min 中。现在要用分子动力学模拟继续进行弛豫。从菜单栏中选取Modules|Discover|Dynamics。将会显示Discover Molecular Dynamics 对话框。,把Ensemble 改为NVT。把温度改为300。把Number of steps 从5000 改为2000，把Trajectory Save 选项改为 Final Structure，

32、点击Run。,字母含义如下：N=固定粒子数V=固定体积E=固定能量T=固定温度P=固定压强H=固定焓,注：要平衡一个准备进行扩散计算的晶胞，NPT 系综是最好的选择。本教程中采用最快的NVT 系综。在一个实际的模拟中，你很可能需要至少50ps 来平衡晶胞。这与系统的大小有关。系统越大，平衡所需时间越长。对NVT 系综来说，当即时更新的图表文档中的能量固定不变时，系统就平衡了。在平衡过程中你也要根据速度来调节温度。,分子动力学的运行和分析 当系统平衡以后，要计算要分子在晶胞中的均方位移，你需要很多帧来分析氧原子往哪里移动。因此现在要再运行另外一个分子动力学模拟并生成一个可以用Discover A

33、nalysis 工具来分析的轨迹文档。之前，运行了一个NVT 系综，现在用NVE 系综。因为就方法而言，NVE 动力学不会被系统的热力学过程干扰。激活cell.xtd。,在Discover Molecular Dynamics 对话框中，把Ensemble 改为NVE。运行的步数也要增加。把Number of steps 改为5000。把Trajectory Save 选项改为 Full。把Frame output every改为250。把Trajectory Save 选项选成Full 意味着轨迹文件不仅输出坐标，还包含其它信息，如温度，能量，速度和晶格参数。有些动力学分析函数只需要坐标作为

34、输入，但均方位移需要全部的输出信息。关于分析函数需要什么样的轨迹输出可以参阅Discover Analysis dialog 帮助主题。按下Run 按钮。关闭Discover Molecular Dynamics 对话框。,计算完成后，就可以开始分析输出文件了。激活cell.xtd。点击Animation 工具条上的Play 按钮。轨迹从1 到20 帧循环，可以让你观察分子动力学模拟过程。当动画结束后，按Stop 按钮。,为了计算氧分子的均方位移，你要把它们同聚合物分子区分开来。这可以通过把它们定义成一组来达到。要选取所有的氧原子，按住ALT 键，双击其中一个。不过，如果一个氧原子在聚合体内部

35、，你要把它同其它氧原子区分开来。最简单的方法是用它们的力场类型来标记它们，只选中那些对应一种特定力场的氧原子。先使氧分子与聚合物清楚地区分。,显示O2的扩散更清楚,选中氧分子中的一个氧原子，右键点选context menu 中的Label。在Label 对话框中，选择ForcefieldType 特性并点击Apply。氧原子被标记为o1o。,从菜单栏中选取Edit|Atom Selection，会显示Atom Selection 对话框。,按OK，关闭对话框。在3D trajectory document上双击左键，去除对O2的选中。,点击Forcite Analysis 对话框中的可用选项C

36、hoose sets箭头，选择oxygen。点击Analyze。关闭Forcite Analysis 对话框。,5.输出数据并计算扩散系数 本教程的最后一部分包括一种电子表格或画图软件的使用。你可以用它来检验均方位移的计算是否正确，在此基础上再来计算扩散系数。复制并粘贴图表文档到你的电子表格中。右键点击plot，并从context menu 中选取Copy。打开新的电子表格，右键点击它并选择Paste。,在你的电子表格中有14列数。第一列是时间，它每隔一列重复出现一次。另外的列里包含所有均方位移的x-，y-和z-分量。在本次计算中你只要前面两列。,提示：在实际计算中，你要检查计算结果是否可靠。

37、你可以画出log(MSD)对log(time)的曲线。如果你的计算收敛了，那么你将得到一条直线。不然，你就要重新计算了。,要计算氧在PDMS 中的扩散率，你要画出MSD 对时间的曲线，拟合后计算斜率。画出MSD 对时间的曲线。线性拟合成直线y=ax+b，记下斜率a。扩散系数由下面式子给出：,其中N 是系统中扩散原子的数目。上式中的微分近似用MSD 对时间微分的比率来代替，也就是曲线的斜率a。由于MSD 的值已经对扩散原子数N 作了平均，所以公式可以简化为：,实验中氧扩散到PDMS 中的扩散系数约为5.7 x 10-5 cm2 s-1，而报道的计算结果在2 x 10-5 cm2 s-1与 8 x

38、 10-5 cm2 s-1 之间(Charati and Stern,1998;Hofmann et al.,2000)。你的计算结果很可能和这个有很大差别，因为你取的聚合体太短，晶格尺度太小，运算时间太短以至于Einstein扩散没有时间发生。参考文献Charati,S.G.;Stern,S.A.Diffusion of Gases in Silicone Polymers:Molecular Dynamics Simulations,Macromolecules,31,5529-5538(1998).Hofmann,D.;Fritz,L.;Ulbrich,J.;Schepers,C;Boehning,M.Detailed-atomistic molecular modeling of small molecule diffusion andsolution processes in polymeric membrane materials,Macromol.TheorySimul.,9,293-327(2000).,Thank you！,