第六章-分子动力学模拟课件

分子动力学模拟分子动力学模拟分子动力学模拟一、一、系综理论系综理论本章主要内容本章主要内容二、二、分子动力学方法分子动力学方法分子动力学模拟分子动力学模拟三、三、模拟细节模拟细节四、四、参量的计算参量的计算五、五、液态水的液态水的MDMD模拟模拟六、六、误差分析误差分析七、七、分子动力学模拟方法的应用分子动力学模拟方法的应用一、系综理论本章主要内容二、分子动力学方法分子动力学模拟三、分子动力学模拟（分子动力学模拟（分子动力学模拟（分子动力学模拟（molecular dynamics simulationmolecular dynamics simulation，简称，简称，简称，简称MDMD）方法首先是由方法首先是由方法首先是由方法首先是由AlderAlder和和和和WainwrightWainwright提出的，现已逐渐成为预测系统特性、提出的，现已逐渐成为预测系统特性、提出的，现已逐渐成为预测系统特性、提出的，现已逐渐成为预测系统特性、验证理论和改进模型的计算工具。验证理论和改进模型的计算工具。验证理论和改进模型的计算工具。验证理论和改进模型的计算工具。MDMD方法的方法的方法的方法的基本思想基本思想基本思想基本思想是把物质看成由原子和分子组成的粒子系统，是把物质看成由原子和分子组成的粒子系统，是把物质看成由原子和分子组成的粒子系统，是把物质看成由原子和分子组成的粒子系统，从该体系的某一假定的位能模型出发，并假定体系粒子的运动遵循经典从该体系的某一假定的位能模型出发，并假定体系粒子的运动遵循经典从该体系的某一假定的位能模型出发，并假定体系粒子的运动遵循经典从该体系的某一假定的位能模型出发，并假定体系粒子的运动遵循经典力学或量子力学描述的规律，若已知粒子的所有受力作用，则可以求解力学或量子力学描述的规律，若已知粒子的所有受力作用，则可以求解力学或量子力学描述的规律，若已知粒子的所有受力作用，则可以求解力学或量子力学描述的规律，若已知粒子的所有受力作用，则可以求解出运动方程而得到系统中全体粒子在相空间中的轨道，然后统计得到系出运动方程而得到系统中全体粒子在相空间中的轨道，然后统计得到系出运动方程而得到系统中全体粒子在相空间中的轨道，然后统计得到系出运动方程而得到系统中全体粒子在相空间中的轨道，然后统计得到系统的热力学参数、结构和输运特性等。也就是由体系的微观性质来求算统的热力学参数、结构和输运特性等。也就是由体系的微观性质来求算统的热力学参数、结构和输运特性等。也就是由体系的微观性质来求算统的热力学参数、结构和输运特性等。也就是由体系的微观性质来求算其宏观性质。其宏观性质。其宏观性质。其宏观性质。属于微观尺度的模拟技术。属于微观尺度的模拟技术。属于微观尺度的模拟技术。属于微观尺度的模拟技术。一、系综理论一、系综理论分子动力学模拟分子动力学模拟 分子动力学模拟（molecular dynau2.1 Newtonian mechanics2.1 Newtonian mechanics In the MD method,every new distribution is derived from the In the MD method,every new distribution is derived from the previous one by using the interactions between the particles.previous one by using the interactions between the particles.The interactions depends on the position of the particles.The interactions depends on the position of the particles.In that potential the particle feels a forceIn that potential the particle feels a forceAccording to Newtons second lawAccording to Newtons second law二、分子动力学方法二、分子动力学方法分子动力学模拟分子动力学模拟2.1 Newtonian mechanics u2.2 Potential energy functions 2.2 Potential energy functions 对一个由对一个由对一个由对一个由N N个原子构成的简单系统，其势能项由下式给出个原子构成的简单系统，其势能项由下式给出个原子构成的简单系统，其势能项由下式给出个原子构成的简单系统，其势能项由下式给出 式中右端第一项是外场（如电场、式中右端第一项是外场（如电场、式中右端第一项是外场（如电场、式中右端第一项是外场（如电场、磁场、声场等）对系统的作用；第二项磁场、声场等）对系统的作用；第二项磁场、声场等）对系统的作用；第二项磁场、声场等）对系统的作用；第二项是两体势即系统中每两个粒子间的相互是两体势即系统中每两个粒子间的相互是两体势即系统中每两个粒子间的相互是两体势即系统中每两个粒子间的相互作用；第三项是三体势，表示系统中每作用；第三项是三体势，表示系统中每作用；第三项是三体势，表示系统中每作用；第三项是三体势，表示系统中每三个粒子间的相互作用三个粒子间的相互作用三个粒子间的相互作用三个粒子间的相互作用 有效两体势有效两体势有效两体势有效两体势它包含多体效应，可很好地反映系统粒子间的相互作用。它包含多体效应，可很好地反映系统粒子间的相互作用。它包含多体效应，可很好地反映系统粒子间的相互作用。它包含多体效应，可很好地反映系统粒子间的相互作用。分子动力学模拟分子动力学模拟2.2 Potential energy functions1 1、Lennard-JoansLennard-Joans势势势势下面仅对简单系统的相互作用模型给予简介下面仅对简单系统的相互作用模型给予简介下面仅对简单系统的相互作用模型给予简介下面仅对简单系统的相互作用模型给予简介分子动力学模拟分子动力学模拟1、Lennard-Joans势下面仅对简单系统的相互作用模2.2.硬球势（硬球势（硬球势（硬球势（Hard-SphereHard-Sphere）3.3.软球势（软球势（软球势（软球势（Soft-SphereSoft-Sphere）通常，通常，通常，通常，v v 是为整数的参数。是为整数的参数。是为整数的参数。是为整数的参数。4.4.方阱势（方阱势（方阱势（方阱势（Square-WellSquare-Well）分子动力学模拟分子动力学模拟2.硬球势（Hard-Sphere）3.软球势（u2.3 Calculations of force,velocity,position2.3 Calculations of force,velocity,position The initial distribution of the Molecular dynamics simulation The initial distribution of the Molecular dynamics simulation is generated in a is generated in a random distributionrandom distribution.In simulation:In simulation:Each particle is also assigned an initial velocity Each particle is also assigned an initial velocity v vi i 分子动力学模拟分子动力学模拟2.3 Calculations of force,velThe force causes an accelerationThe force causes an accelerationWhich in turn modifies the initial velocity Which in turn modifies the initial velocity v vi i as asAnd modifies the initial position And modifies the initial position r ri i as asThe MD simulation can describe systems that evolve in time.The The MD simulation can describe systems that evolve in time.The new positions are derived from the Newtonian law of motions and new positions are derived from the Newtonian law of motions and therefore deterministic.therefore deterministic.分子动力学模拟分子动力学模拟The force causes an acceleratiu2.4 Equations of motion2.4 Equations of motion 为了在计算机上解运动方程，必须为微分方程建立一个为了在计算机上解运动方程，必须为微分方程建立一个为了在计算机上解运动方程，必须为微分方程建立一个为了在计算机上解运动方程，必须为微分方程建立一个有限差分格式有限差分格式有限差分格式有限差分格式，从差分方程中再导出，从差分方程中再导出，从差分方程中再导出，从差分方程中再导出位置和速度位置和速度位置和速度位置和速度的递推关系的递推关系的递推关系的递推关系式。这些算法是一步一步执行的，先算式。这些算法是一步一步执行的，先算式。这些算法是一步一步执行的，先算式。这些算法是一步一步执行的，先算t t 时刻的位置和速度，时刻的位置和速度，时刻的位置和速度，时刻的位置和速度，然后在此基础上计算然后在此基础上计算然后在此基础上计算然后在此基础上计算t+1t+1时刻的位置和速度。时刻的位置和速度。时刻的位置和速度。时刻的位置和速度。微分方程最为直接的离散化格式来自泰勒展开微分方程最为直接的离散化格式来自泰勒展开微分方程最为直接的离散化格式来自泰勒展开微分方程最为直接的离散化格式来自泰勒展开:1、有限差分方法、有限差分方法-预测校正法预测校正法分子动力学模拟分子动力学模拟2.4 Equations of motion 2、有限差分方法、有限差分方法-Verlet算法算法、VerletVerlet算法的一般形式算法的一般形式算法的一般形式算法的一般形式为了用数值方法求解微分方程，为了用数值方法求解微分方程，为了用数值方法求解微分方程，为了用数值方法求解微分方程，采用有限差分方法离散化得到差分格式采用有限差分方法离散化得到差分格式采用有限差分方法离散化得到差分格式采用有限差分方法离散化得到差分格式两式相加得两式相加得两式相加得两式相加得 分子动力学模拟分子动力学模拟2、有限差分方法-Verlet算法、Verlet算法的一般 上式提供了一个方法，从粒子在前两步（上式提供了一个方法，从粒子在前两步（上式提供了一个方法，从粒子在前两步（上式提供了一个方法，从粒子在前两步（t t和和和和t-ht-h）时刻）时刻）时刻）时刻的位置以及的位置以及的位置以及的位置以及t t时刻的作用力来得到粒子在时刻的作用力来得到粒子在时刻的作用力来得到粒子在时刻的作用力来得到粒子在t+ht+h时刻的位置。时刻的位置。时刻的位置。时刻的位置。速度对轨道计算没有关系，但对速度对轨道计算没有关系，但对速度对轨道计算没有关系，但对速度对轨道计算没有关系，但对估算动能及速度相关函数（用来研估算动能及速度相关函数（用来研估算动能及速度相关函数（用来研估算动能及速度相关函数（用来研究粒子的输运特性）非常有用。速究粒子的输运特性）非常有用。速究粒子的输运特性）非常有用。速究粒子的输运特性）非常有用。速度形式为度形式为度形式为度形式为分子动力学模拟分子动力学模拟 上式提供了一个方法，从粒子在前两步（t和t-、Verlet Verlet 蛙跳格式（蛙跳格式（蛙跳格式（蛙跳格式（leap-frogleap-frog）The velocity varies so one has to choose a reasonable The velocity varies so one has to choose a reasonable average value to be used.The velocity at the middle of average value to be used.The velocity at the middle of the step ought to be a good compromise,the step ought to be a good compromise,分子动力学模拟分子动力学模拟、Verlet 蛙跳格式（leap-frog）The ve、VerletVerlet速度相加形式速度相加形式速度相加形式速度相加形式分子动力学模拟分子动力学模拟、Verlet速度相加形式分子动力学模拟分子动力学模拟分子动力学模拟分子动力学模拟u3.1 3.1 初始化位型初始化位型 模拟开始时首先要初始化位型分布，既首先要给定微模拟开始时首先要初始化位型分布，既首先要给定微模拟开始时首先要初始化位型分布，既首先要给定微模拟开始时首先要初始化位型分布，既首先要给定微元中元中元中元中分子的初始位置和初始速度分子的初始位置和初始速度分子的初始位置和初始速度分子的初始位置和初始速度。、简单立方晶格、简单立方晶格、体心立方晶格、体心立方晶格、面心立方晶格、面心立方晶格abc三、模拟细节三、模拟细节分子动力学模拟分子动力学模拟3.1 初始化位型 模拟开始时首先abcabcabcabcabcabcabcabcabc分子动力学模拟分子动力学模拟abcabcabcabcabcabcabcabcabc分子动初始速度初始速度初始速度初始速度 模拟时，各粒子的初始速度按模拟时，各粒子的初始速度按模拟时，各粒子的初始速度按模拟时，各粒子的初始速度按麦克斯韦速度分布麦克斯韦速度分布麦克斯韦速度分布麦克斯韦速度分布取样。取样。取样。取样。Maxwells distribution law of velocityMaxwells distribution law of velocity 指平衡状态下理想气体分子指平衡状态下理想气体分子指平衡状态下理想气体分子指平衡状态下理想气体分子速度分布的统计规律。在平衡状速度分布的统计规律。在平衡状速度分布的统计规律。在平衡状速度分布的统计规律。在平衡状态下，分布在任一速率区间内的态下，分布在任一速率区间内的态下，分布在任一速率区间内的态下，分布在任一速率区间内的分子数与总分子数的比率为分子数与总分子数的比率为分子数与总分子数的比率为分子数与总分子数的比率为分子动力学模拟分子动力学模拟初始速度 模拟时，各粒子的初始u3.2 Periodic boundary conditions 3.2 Periodic boundary conditions 通常选取小数体系（几十个到数千个分子）作为研究通常选取小数体系（几十个到数千个分子）作为研究通常选取小数体系（几十个到数千个分子）作为研究通常选取小数体系（几十个到数千个分子）作为研究对象，但由于位于表面和内部的分子受力差别较大，将会对象，但由于位于表面和内部的分子受力差别较大，将会对象，但由于位于表面和内部的分子受力差别较大，将会对象，但由于位于表面和内部的分子受力差别较大，将会产生表面效应。产生表面效应。产生表面效应。产生表面效应。为了消除此效应，同时建造为了消除此效应，同时建造为了消除此效应，同时建造为了消除此效应，同时建造出一个准无穷大体积，使其可以出一个准无穷大体积，使其可以出一个准无穷大体积，使其可以出一个准无穷大体积，使其可以更精确地代表宏观系统，必须引更精确地代表宏观系统，必须引更精确地代表宏观系统，必须引更精确地代表宏观系统，必须引入周期性边界条件。入周期性边界条件。入周期性边界条件。入周期性边界条件。把小数体系看作一个中心元胞，此中心元胞被与中心元把小数体系看作一个中心元胞，此中心元胞被与中心元把小数体系看作一个中心元胞，此中心元胞被与中心元把小数体系看作一个中心元胞，此中心元胞被与中心元胞相同的其它元胞所包围，当某分子离开中心元胞时，该分胞相同的其它元胞所包围，当某分子离开中心元胞时，该分胞相同的其它元胞所包围，当某分子离开中心元胞时，该分胞相同的其它元胞所包围，当某分子离开中心元胞时，该分子将在整个点格上移动以致它从中心元胞对面重新进入这个子将在整个点格上移动以致它从中心元胞对面重新进入这个子将在整个点格上移动以致它从中心元胞对面重新进入这个子将在整个点格上移动以致它从中心元胞对面重新进入这个元胞。元胞。元胞。元胞。周期性边界条件周期性边界条件周期性边界条件周期性边界条件：分子动力学模拟分子动力学模拟3.2 Periodic boundary conditioD D基本元胞基本元胞A AB BF FE EHHG GC C分子动力学模拟分子动力学模拟D基本元胞ABFEHGC分子动力学模拟u3.3 Calculation of3.3 Calculation of interactions interactions i,j i,j 两个离子的相互作用势两个离子的相互作用势两个离子的相互作用势两个离子的相互作用势体系总势能可表示为：体系总势能可表示为：体系总势能可表示为：体系总势能可表示为：在计算势能时必须考虑相互作用的力程，对粒子数为在计算势能时必须考虑相互作用的力程，对粒子数为在计算势能时必须考虑相互作用的力程，对粒子数为在计算势能时必须考虑相互作用的力程，对粒子数为N N的模拟系统，原则上任何两个粒子间都有相互作用，在的模拟系统，原则上任何两个粒子间都有相互作用，在的模拟系统，原则上任何两个粒子间都有相互作用，在的模拟系统，原则上任何两个粒子间都有相互作用，在计算体系势能时须进行计算体系势能时须进行计算体系势能时须进行计算体系势能时须进行N N（N-1N-1）/2/2次运算。次运算。次运算。次运算。为了提高计算效率，实际模拟为了提高计算效率，实际模拟为了提高计算效率，实际模拟为了提高计算效率，实际模拟过程中通常在一个方便的力程上过程中通常在一个方便的力程上过程中通常在一个方便的力程上过程中通常在一个方便的力程上将位势截断，用以减少计算势能将位势截断，用以减少计算势能将位势截断，用以减少计算势能将位势截断，用以减少计算势能所消耗的时间。所消耗的时间。所消耗的时间。所消耗的时间。分子动力学模拟分子动力学模拟3.3 Calculation of interaction 对势能的最大贡献来自于粒子的近邻区域，位势截断对势能的最大贡献来自于粒子的近邻区域，位势截断对势能的最大贡献来自于粒子的近邻区域，位势截断对势能的最大贡献来自于粒子的近邻区域，位势截断常用的方法是常用的方法是常用的方法是常用的方法是球形截断球形截断球形截断球形截断，截断半径一般取，截断半径一般取，截断半径一般取，截断半径一般取2.52.5或或或或3.6 3.6，对，对，对，对截断距离之外分子间相互作用能按平均密度近似的方法进截断距离之外分子间相互作用能按平均密度近似的方法进截断距离之外分子间相互作用能按平均密度近似的方法进截断距离之外分子间相互作用能按平均密度近似的方法进行校正。行校正。行校正。行校正。分子动力学模拟分子动力学模拟 对势能的最大贡献来自于粒子的近邻区域，位势截断Cell methodCell methodNeighbor list method Neighbor list method 1、近邻表、近邻表-Verlet近邻表近邻表2、近邻表、近邻表-网格近邻表网格近邻表分子动力学模拟分子动力学模拟Cell methodNeighbor list metu3.4 3.4 标度与趋衡标度与趋衡 正则系综正则系综正则系综正则系综的的的的MDMD模拟，模拟过程中必须保证体系温度模拟，模拟过程中必须保证体系温度模拟，模拟过程中必须保证体系温度模拟，模拟过程中必须保证体系温度恒定，所以在每一时间步要对分子的速度进行标度使其满恒定，所以在每一时间步要对分子的速度进行标度使其满恒定，所以在每一时间步要对分子的速度进行标度使其满恒定，所以在每一时间步要对分子的速度进行标度使其满足正则系综的条件，标度因子为足正则系综的条件，标度因子为足正则系综的条件，标度因子为足正则系综的条件，标度因子为 ，其中，其中，其中，其中T T是指定的体系温度，是指定的体系温度，是指定的体系温度，是指定的体系温度，T*T*是每一时间步的瞬时动能温度。是每一时间步的瞬时动能温度。是每一时间步的瞬时动能温度。是每一时间步的瞬时动能温度。标度标度标度标度：分子动力学模拟分子动力学模拟3.4 标度与趋衡 正则系综的M 一种按上述方法建立的系统不会具有所要的能量，而且很可能一种按上述方法建立的系统不会具有所要的能量，而且很可能一种按上述方法建立的系统不会具有所要的能量，而且很可能一种按上述方法建立的系统不会具有所要的能量，而且很可能这个状态不对应一个平衡态，而所期望的体系参量是在平衡态通过这个状态不对应一个平衡态，而所期望的体系参量是在平衡态通过这个状态不对应一个平衡态，而所期望的体系参量是在平衡态通过这个状态不对应一个平衡态，而所期望的体系参量是在平衡态通过统计平均得到的，所以为使系统达到平衡需要一个趋衡阶段。统计平均得到的，所以为使系统达到平衡需要一个趋衡阶段。统计平均得到的，所以为使系统达到平衡需要一个趋衡阶段。统计平均得到的，所以为使系统达到平衡需要一个趋衡阶段。趋衡趋衡趋衡趋衡：对于正则系综，体系趋衡是通过增减体系能量来实现的。如果对于正则系综，体系趋衡是通过增减体系能量来实现的。如果对于正则系综，体系趋衡是通过增减体系能量来实现的。如果对于正则系综，体系趋衡是通过增减体系能量来实现的。如果在进行足够多的运算步数以后，系统能持续给出确定的平均动能和在进行足够多的运算步数以后，系统能持续给出确定的平均动能和在进行足够多的运算步数以后，系统能持续给出确定的平均动能和在进行足够多的运算步数以后，系统能持续给出确定的平均动能和平均势能的数值，那么就认为平衡已经建立。平均势能的数值，那么就认为平衡已经建立。平均势能的数值，那么就认为平衡已经建立。平均势能的数值，那么就认为平衡已经建立。分子动力学模拟分子动力学模拟 一种按上述方法建立的系统不会具有所要的能量，而u3.5 The MD units 3.5 The MD units 这里我们以能量这里我们以能量这里我们以能量这里我们以能量 （焦耳，（焦耳，（焦耳，（焦耳，J J）和粒子间距离）和粒子间距离）和粒子间距离）和粒子间距离 （米，（米，（米，（米，mm）为基本单位来对其它量进行约化，使其成为无量纲量。仅举为基本单位来对其它量进行约化，使其成为无量纲量。仅举为基本单位来对其它量进行约化，使其成为无量纲量。仅举为基本单位来对其它量进行约化，使其成为无量纲量。仅举例如下：例如下：例如下：例如下：粒子数密度：粒子数密度：粒子数密度：粒子数密度：温度：温度：温度：温度：能量：能量：能量：能量：压强：压强：压强：压强：时间：时间：时间：时间：力：力：力：力：分子动力学模拟分子动力学模拟3.5 The MD units 这里我们 A simulation run produces the raw data in form of a A simulation run produces the raw data in form of a very large very large disk filedisk file,The the complete state of the system,The the complete state of the system at each step.at each step.The disk processed after the simulation is finished.It The disk processed after the simulation is finished.It contains at least contains at least all the positions and velocities of all all the positions and velocities of all particlesparticles.This information is sufficient to calculate all the.This information is sufficient to calculate all the properties of the system.However,it is more economical properties of the system.However,it is more economical to calculate properties during the simulation and store to calculate properties during the simulation and store them in the than reading the calculating them afterwards.them in the than reading the calculating them afterwards.The The potentials and forcespotentials and forces certainly belong to this class certainly belong to this class because they are needed at each time step anyway in order because they are needed at each time step anyway in order to calculate the new position and velocities of the particles.to calculate the new position and velocities of the particles.四、参量的计算四、参量的计算分子动力学模拟分子动力学模拟 A simulation run produce1 1内能内能内能内能 模拟体系的内能由两部分组成模拟体系的内能由两部分组成模拟体系的内能由两部分组成模拟体系的内能由两部分组成2 2比热比热比热比热 涨落是体系中热力学量在不同时刻出现的波动。系统涨落是体系中热力学量在不同时刻出现的波动。系统涨落是体系中热力学量在不同时刻出现的波动。系统涨落是体系中热力学量在不同时刻出现的波动。系统比热可以通过能量的涨落导出，表示如下比热可以通过能量的涨落导出，表示如下比热可以通过能量的涨落导出，表示如下比热可以通过能量的涨落导出，表示如下分子动力学模拟分子动力学模拟1内能 模拟体系的内能由两部分组成2比热 涨3.3.压强压强压强压强 体系的压强可通过体系的压强可通过体系的压强可通过体系的压强可通过维里定理维里定理维里定理维里定理导出：导出：导出：导出：分子动力学模拟分子动力学模拟3.压强 体系的压强可通过维里定理导出：分子动力学模拟4 4扩散系数扩散系数扩散系数扩散系数 在模拟过程中，体系的迁移系数可以通过在模拟过程中，体系的迁移系数可以通过在模拟过程中，体系的迁移系数可以通过在模拟过程中，体系的迁移系数可以通过平衡相关函数平衡相关函数平衡相关函数平衡相关函数导导导导出。由出。由出。由出。由Green-KuboGreen-Kubo方法，计算扩散系数的公式如下：方法，计算扩散系数的公式如下：方法，计算扩散系数的公式如下：方法，计算扩散系数的公式如下：由于受计算机硬件的限制，在模拟过程中不可能对时间积由于受计算机硬件的限制，在模拟过程中不可能对时间积由于受计算机硬件的限制，在模拟过程中不可能对时间积由于受计算机硬件的限制，在模拟过程中不可能对时间积分至无穷。为了选择合适的积分上限，在模拟过程中应首先计分至无穷。为了选择合适的积分上限，在模拟过程中应首先计分至无穷。为了选择合适的积分上限，在模拟过程中应首先计分至无穷。为了选择合适的积分上限，在模拟过程中应首先计算各水分子的速度自相关函数，从水体系归一化的速度自相关算各水分子的速度自相关函数，从水体系归一化的速度自相关算各水分子的速度自相关函数，从水体系归一化的速度自相关算各水分子的速度自相关函数，从水体系归一化的速度自相关函数可以看出在什么长度的时间间隔时，速度自相关函数逐渐函数可以看出在什么长度的时间间隔时，速度自相关函数逐渐函数可以看出在什么长度的时间间隔时，速度自相关函数逐渐函数可以看出在什么长度的时间间隔时，速度自相关函数逐渐趋零，也就是说在什么时间范围内速度基本上不再相关，那么趋零，也就是说在什么时间范围内速度基本上不再相关，那么趋零，也就是说在什么时间范围内速度基本上不再相关，那么趋零，也就是说在什么时间范围内速度基本上不再相关，那么求扩散系数所选取的积分区间就可以确定。求扩散系数所选取的积分区间就可以确定。求扩散系数所选取的积分区间就可以确定。求扩散系数所选取的积分区间就可以确定。为速度自相关函数（为速度自相关函数（为速度自相关函数（为速度自相关函数（VACFVACF）分子动力学模拟分子动力学模拟4扩散系数 在模拟过程中，体系的迁移系数可以通过平衡5 5、粘滞系数、粘滞系数、粘滞系数、粘滞系数 利用平衡分子动力学模拟方法计算平衡系统的粘滞利用平衡分子动力学模拟方法计算平衡系统的粘滞利用平衡分子动力学模拟方法计算平衡系统的粘滞利用平衡分子动力学模拟方法计算平衡系统的粘滞系数的基础是系数的基础是系数的基础是系数的基础是Green-KoboGreen-Kobo理论。根据该理论，粘滞系数理论。根据该理论，粘滞系数理论。根据该理论，粘滞系数理论。根据该理论，粘滞系数由下式给定由下式给定由下式给定由下式给定 分子动力学模拟分子动力学模拟5、粘滞系数 利用平衡分子动力学模拟方法计算6 6、径向分布函数、径向分布函数、径向分布函数、径向分布函数-g g(r r)径向分布函数径向分布函数径向分布函数径向分布函数是与时间无关的关联性的量度，精确地说就是与时间无关的关联性的量度，精确地说就是与时间无关的关联性的量度，精确地说就是与时间无关的关联性的量度，精确地说就是与中心粒子间距为是与中心粒子间距为是与中心粒子间距为是与中心粒子间距为rr+dr rr+dr 的体积元内分子的平均数密度。的体积元内分子的平均数密度。的体积元内分子的平均数密度。的体积元内分子的平均数密度。在实际模拟过程中我们以元胞边在实际模拟过程中我们以元胞边在实际模拟过程中我们以元胞边在实际模拟过程中我们以元胞边长的长的长的长的1/21/2为半径作圆球，将其分为半径作圆球，将其分为半径作圆球，将其分为半径作圆球，将其分割为厚度相等的割为厚度相等的割为厚度相等的割为厚度相等的100100个球壳，计个球壳，计个球壳，计个球壳，计算分子间距算分子间距算分子间距算分子间距R R满足满足满足满足 （ijij，且，且，且，且i,j i,j 依次取值）的分子依次取值）的分子依次取值）的分子依次取值）的分子对个数对个数对个数对个数N(i)N(i)，i=1,2,.100.i=1,2,.100.分子动力学模拟分子动力学模拟6、径向分布函数-g(r)分子动力学模拟方法的基本步骤分子动力学模拟方法的基本步骤分子动力学模拟方法的基本步骤分子动力学模拟方法的基本步骤1.1.确定研究对象确定研究对象确定研究对象确定研究对象，如确定所研究的体系是属于微正则系综、，如确定所研究的体系是属于微正则系综、，如确定所研究的体系是属于微正则系综、，如确定所研究的体系是属于微正则系综、正则系综、巨正则系综呢还是属于等温等压系综等等；正则系综、巨正则系综呢还是属于等温等压系综等等；正则系综、巨正则系综呢还是属于等温等压系综等等；正则系综、巨正则系综呢还是属于等温等压系综等等；2.2.建立适当的位能模型建立适当的位能模型建立适当的位能模型建立适当的位能模型，如假定所研究的体系位势形式满足，如假定所研究的体系位势形式满足，如假定所研究的体系位势形式满足，如假定所研究的体系位势形式满足硬球势、软球势或勒纳德硬球势、软球势或勒纳德硬球势、软球势或勒纳德硬球势、软球势或勒纳德琼斯势等等；琼斯势等等；琼斯势等等；琼斯势等等；3.3.初始化位形和速度初始化位形和速度初始化位形和速度初始化位形和速度；4.4.建立体系粒子的运动方程并进行求解建立体系粒子的运动方程并进行求解建立体系粒子的运动方程并进行求解建立体系粒子的运动方程并进行求解，求解时需对位势进，求解时需对位势进，求解时需对位势进，求解时需对位势进行截断，并应用周期性边界条件和最小影像数约定，每一时行截断，并应用周期性边界条件和最小影像数约定，每一时行截断，并应用周期性边界条件和最小影像数约定，每一时行截断，并应用周期性边界条件和最小影像数约定，每一时间步对速度进行标度以使体系达到趋衡；间步对速度进行标度以使体系达到趋衡；间步对速度进行标度以使体系达到趋衡；间步对速度进行标度以使体系达到趋衡；5.5.通过统计平均来求解一些我们感兴趣的通过统计平均来求解一些我们感兴趣的通过统计平均来求解一些我们感兴趣的通过统计平均来求解一些我们感兴趣的热力学参量热力学参量热力学参量热力学参量等。等。等。等。五、液态水的五、液态水的MDMD模拟模拟分子动力学模拟分子动力学模拟分子动力学模拟方法的基本步骤1.确定研究对象，如确定所研1 1模拟参数模拟参数模拟参数模拟参数 本文将水系统视为正则系综进行本文将水系统视为正则系综进行本文将水系统视为正则系综进行本文将水系统视为正则系综进行MDMD模拟，模拟所用分子模拟，模拟所用分子模拟，模拟所用分子模拟，模拟所用分子数为数为数为数为256256，水分子初始位型按面心立方晶格分布，晶格所占空，水分子初始位型按面心立方晶格分布，晶格所占空，水分子初始位型按面心立方晶格分布，晶格所占空，水分子初始位型按面心立方晶格分布，晶格所占空间体积由水的密度导出。各分子起始速度按间体积由水的密度导出。各分子起始速度按间体积由水的密度导出。各分子起始速度按间体积由水的密度导出。各分子起始速度按MaxwellMaxwell分布取分布取分布取分布取样。模拟温度分别取样。模拟温度分别取样。模拟温度分别取样。模拟温度分别取=273K=273K、293K293K、313K313K、333K333K、353K353K、373K373K。模拟过程中为保证体系温度恒定，在每一时间步对分。模拟过程中为保证体系温度恒定，在每一时间步对分。模拟过程中为保证体系温度恒定，在每一时间步对分。模拟过程中为保证体系温度恒定，在每一时间步对分子的速度进行标度；计算体系位能时采用球形截断法，截断子的速度进行标度；计算体系位能时采用球形截断法，截断子的速度进行标度；计算体系位能时采用球形截断法，截断子的速度进行标度；计算体系位能时采用球形截断法，截断半径为半径为半径为半径为2.52.5，对截断距离之外分子间相互作用能按平均密度，对截断距离之外分子间相互作用能按平均密度，对截断距离之外分子间相互作用能按平均密度，对截断距离之外分子间相互作用能按平均密度近似的方法进行校正；运动方程的求解采用近似的方法进行校正；运动方程的求解采用近似的方法进行校正；运动方程的求解采用近似的方法进行校正；运动方程的求解采用Velocity Verlet Velocity Verlet AlgorithmAlgorithm算法；时间步长为算法；时间步长为算法；时间步长为算法；时间步长为210210-15-15s s，每次模拟的总步数为，每次模拟的总步数为，每次模拟的总步数为，每次模拟的总步数为10105 5步步步步 ，其中前，其中前，其中前，其中前5000050000步使体系趋衡，后步使体系趋衡，后步使体系趋衡，后步使体系趋衡，后5000050000步用来统计各种步用来统计各种步用来统计各种步用来统计各种平衡性质和结构性质；模拟过程中各参量均采用约化单位。平衡性质和结构性质；模拟过程中各参量均采用约化单位。平衡性质和结构性质；模拟过程中各参量均采用约化单位。平衡性质和结构性质；模拟过程中各参量均采用约化单位。分子动力学模拟分子动力学模拟1模拟参数 本文将水系统视为正则系综进行MD模拟，计算体系的内能有着重要的意义，不仅是可以通过能量求出计算体系的内能有着重要的意义，不仅是可以通过能量求出计算体系的内能有着重要的意义，不仅是可以通过能量求出计算体系的内能有着重要的意义，不仅是可以通过能量求出比热、表面张力等热力学量，更为重要的是内能（或势能）的比热、表面张力等热力学量，更为重要的是内能（或势能）的比热、表面张力等热力学量，更为重要的是内能（或势能）的比热、表面张力等热力学量，更为重要的是内能（或势能）的演变是我们在模拟过程中判断体系平衡的依据。体系是否达到演变是我们在模拟过程中判断体系平衡的依据。体系是否达到演变是我们在模拟过程中判断体系平衡的依据。体系是否达到演变是我们在模拟过程中判断体系平衡的依据。体系是否达到平衡直接关系到热力学量的真实性、稳定性和准确性。平衡直接关系到热力学量的真实性、稳定性和准确性。平衡直接关系到热力学量的真实性、稳定性和准确性。平衡直接关系到热力学量的真实性、稳定性和准确性。2 2水的内能和比热水的内能和比热水的内能和比热水的内能和比热分子动力学模拟分子动力学模拟 计算体系的内能有着重要的意义，不仅是可以通过能量求出表表表表1 1 1 1 水的势能模拟结果与文献值和实验值的比较水的势能模拟结果与文献值和实验值的比较水的势能模拟结果与文献值和实验值的比较水的势能模拟结果与文献值和实验值的比较方法方法方法方法实验值实验值实验值实验值3434本文本文本文本文MDMDMCMC3535MDMD2424势能（势能（势能（势能（kJ/molkJ/mol）-44.77*-44.77*-41.99-41.99-36.09-36.09-42.20-42.20模型模型模型模型 L-JL-JTIPTIPSPCSPC系综系综系综系综 NVTNVTNVTNVTEVNEVN表表表表2 2 2 2 水的比热的模拟结果和实验结果的比较水的比热的模拟结果和实验结果的比较水的比热的模拟结果和实验结果的比较水的比热的模拟结果和实验结果的比较温度（温度（温度（温度（KK）273273293293313313333333353353373373模拟值模拟值模拟值模拟值（J/gKJ/gK）4.01774.01773.97853.97853.95923.95923.94133.94133.93293.93293.92543.9254实验值实验值实验值实验值（J/gKJ/gK）4.20744.20744.18014.18014.17884.17884.18554.18554.19854.1985 分子动力学模拟分子动力学模拟表1 水的势能模拟结果与文献值和实验值的比较方法实验值33 3水的扩散系数和粘滞系数水的扩散系数和粘滞系数水的扩散系数和粘滞系数水的扩散系数和粘滞系数 表表表表3 3 3 3 水的扩散系数模拟结果与文献值和实验值的比较水的扩散系数模拟结果与文献值和实验值的比较水的扩散系数模拟结果与文献值和实验值的比较水的扩散系数模拟结果与文献值和实验值的比较方法方法方法方法实验值实验值实验值实验值3434本文本文本文本文MDMDMDMD3535MDMD3636扩散系数扩散系数扩散系数扩散系数(10(10-5-5cmcm2 2/s)/s)2.32.32.9622.9622.952.952.82.8模型模型模型模型 L-JL-JSPCSPCMCYMCY系综系综系综系综 NVTNVTNVTNVTNVTNVT分子动力学模拟分子动力学模拟3水的扩散系数和粘滞系数 表3 水的扩散系数模拟结果与文4 4水的径向分布函数水的径向分布函数水的径向分布函数水的径向分布函数 径向分布函数是对物质结构的反映，在研究物质特别是径向分布函数是对物质结构的反映，在研究物质特别是径向分布函数是对物质结构的反映，在研究物质特别是径向分布函数是对物质结构的反映，在研究物质特别是液态时有着重要的作用。液态时有着重要的作用。液态时有着重要的作用。液态时有着重要的作用。分子动力学模拟分子动力学模拟4水的径向分布函数 径向分布函数是对物质结构1 1、位能模型、位能模型、位能模型、位能模型 位能模型是对体系中分子或离子之间相互作用势的反映，位能模型是对体系中分子或离子之间相互作用势的反映，位能模型是对体系中分子或离子之间相互作用势的反映，位能模型是对体系中分子或离子之间相互作用势的反映，位能模型选择的正确与否决定了模拟的成败。但是许多情况位能模型选择的正确与否决定了模拟的成败。但是许多情况位能模型选择的正确与否决定了模拟的成败。但是许多情况位能模型选择的正确与否决定了模拟的成败。但是许多情况下由于研究对象的复杂性，使对问题的认识还不够深入，导下由于研究对象的复杂性，使对问题的认识还不够深入，导下由于研究对象的复杂性，使对问题的认识还不够深入，导下由于研究对象的复杂性，使对问题的认识还不够深入，导致所选择的位能模型和实际情况存在一定差距；同时受计算致所选择的位能模型和实际情况存在一定差距；同时受计算致所选择的位能模型和实际情况存在一定差距；同时受计算致所选择的位能模型和实际情况存在一定差距；同时受计算机硬件的限制，机硬件的限制，机硬件的限制，机硬件的限制，也不可能选择更接近实际的复杂模型，如也不可能选择更接近实际的复杂模型，如也不可能选择更接近实际的复杂模型，如也不可能选择更接近实际的复杂模型，如多体相互作用模型，从而导致误差出现。多体相互作用模型，从而导致误差出现。多体相互作用模型，从而导致误差出现。多体相互作用模型，从而导致误差出现。模型的名称ST2TIPSCITIPS2SPCBFPPC模型的类型经验公式经验公式经验公经验公式式从头计从头计算算经验公经验公式式经验公经验公式式半经验公半经验公式式偶极子及偶极子及交换作用交换作用表表表表1 1 1 1 关于纯水结构的理论模型关于纯水结构的理论模型关于纯水结构的理论模型关于纯水结构的理论模型六、六、误差分析误差分析分子动力学模拟分子动力学模拟1、位能模型 位能模型是对体系中分子或离子之间相互作用2 2、模拟分子数、模拟分子数、模拟分子数、模拟分子数 分子动力学模拟的主要任务就是要从微观来求算体系分子动力学模拟的主要任务就是要从微观来求算体系分子动力学模拟的主要任务就是要从微观来求算体系分子动力学模拟的主要任务就是要从微观来求算体系的宏观性质，所以模拟过程中参与模拟的分子数越多，模的宏观性质，所以模拟过程中参与模拟的分子数越多，模的宏观性质，所以模拟过程中参与模拟的分子数越多，模的宏观性质，所以模拟过程中参与模拟的分子数越多，模拟结果就越真实，越接近体系的宏观属性。例如拟结果就越真实，越接近体系的宏观属性。例如拟结果就越真实，越接近体系的宏观属性。例如拟结果就越真实，越接近体系的宏观属性。例如TanakaTanaka等等等等人采用人采用人采用人采用1.36101.36104 4个水分子对水体系进行了模拟，其扩散系数个水分子对水体系进行了模拟，其扩散系数个水分子对水体系进行了模拟，其扩散系数个水分子对水体系进行了模拟，其扩散系数更加接近实验值。但也应考虑到：模拟过程中由于分子受更加接近实验值。但也应考虑到：模拟过程中由于分子受更加接近实验值。但也应考虑到：模拟过程中由于分子受更加接近实验值。但也应考虑到：模拟过程中由于分子受力和势能的计算需要涉及微元内所有的分子对，所以计算力和势能的计算需要涉及微元内所有的分子对，所以计算力和势能的计算需要涉及微元内所有的分子对，所以计算力和势能的计算需要涉及微元内所有的分子对，所以计算量很大，非常消耗机时，一般占总计算量的量很大，非常消耗机时，一般占总计算量的量很大，非常消耗机时，一般占总计算量的量很大，非常消耗机时，一般占总计算量的80%80%左右。模左右。模左右。模左右。模拟体系分子数扩大一倍，总模拟时间并不仅增加一倍，而拟体系分子数扩大一倍，总模拟时间并不仅增加一倍，而拟体系分子数扩大一倍，总模拟时间并不仅增加一倍，而拟体系分子数扩大一倍，总模拟时间并不仅增加一倍，而是呈级数形式上升。是呈级数形式上升。是呈级数形式上升。是呈级数形式上升。要想解决大数粒子体系的模拟问题，还须引入一些模拟要想解决大数粒子体系的模拟问题，还须引入一些模拟要想解决大数粒子体系的模拟问题，还须引入一些模拟要想解决大数粒子体系的模拟问题，还须引入一些模拟技巧，如技巧，如技巧，如技巧，如Verlet Neighbour ListVerlet Neighbour List和和和和Link ListLink List等方法。等方法。等方法。等方法。分子动力学模拟分子动力学模拟2、模拟分子数 分子动力学模拟的主要任务就是要粒子数粒子数粒子数粒子数本文模拟本文模拟本文模拟本文模拟LCLCVerlet Verlet tabletableLLCLLC2562560.68sec0.68sec 13721372 0.45 sec0.45 sec0.072 sec0.072 sec0.037 sec0.037 sec20482048 0.57 sec0.57 sec0.114 sec0.114 sec0.058 sec0.058 sec40004000 0.094 sec0.094 sec69126912 1.59 sec1.59 sec0.529 sec0.529 sec0.164 sec0.164 sec表表表表2 2 2 2 不同方法每时间步所消耗不同方法每时间步所消耗不同方法每时间步所消耗不同方法每时间步所消耗CPUCPUCPUCPU时间的比较时间的比较时间的比较时间的比较分子动力学模拟分子动力学模拟粒子数本文模拟LCVerlet tableLLC2560.63 3、求解运动方程的算法、求解运动方程的算法、求解运动方程的算法、求解运动方程的算法 求解分子运动方程的方法有两种，一种是求解分子运动方程的方法有两种，一种是求解分子运动方程的方法有两种，一种是求解分子运动方程的方法有两种，一种是Verlet AlgorithmVerlet Algorithm，它是对二阶微分运动方程直接求解的一种方法；另一种则是它是对二阶微分运动方程直接求解的一种方法；另一种则是它是对二阶微分运动方程直接求解的一种方法；另一种则是它是对二阶微分运动方程直接求解的一种方法；另一种则是GearGear预测预测预测预测校正算法（校正算法（校正算法（校正算法（Gear Predictor-Corrector AlgorithmGear Predictor-Corrector Algorithm），），），），它是通过泰勒级数展开的方法实现的。它是通过泰勒级数展开的方法实现的。它是通过泰勒级数展开的方法实现的。它是通过泰勒级数展开的方法实现的。在时间步长较小的情况下，通过高阶在时间步长较小的情况下，通过高阶在时间步长较小的情况下，通过高阶在时间步长较小的情况下，通过高阶GearGear算法得到的结果算法得到的结果算法得