1、Molecular Dynamics分子分子动力学动力学 基础基础Chapter 34学时分子动力学？分子动力学？分子动力学是一门结合物理，数学和化学分子动力学是一门结合物理，数学和化学,研究处于原子、分子状态的固体或液体的微研究处于原子、分子状态的固体或液体的微观动力学特性的科学，一种计算机模拟实验观动力学特性的科学，一种计算机模拟实验方法，是研究凝聚态系统的有力工具。方法，是研究凝聚态系统的有力工具。分子动力学是一套分子模拟方法，该方法主分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽以
2、在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。学量和其他宏观性质。Molecular Dynamics 是从统计物理学发展起来的一种描述纳米科技这类处于原子、是从统计物理学发展起来的一种描述纳米科技这类处于原子、分子状态的固体或液体的动力学特性的微观现象的有效方法。分子状态的固体或液体的动力学特性的微观现象的有效方法。分子动力学总是假定原子的运动服从某种确定的描述，这种分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方
3、程、拉格朗日方程或哈密顿方程所确定的描描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和核子的量子效应和Born-Oppenheimer绝热近似下，分子动绝热近似下，分子动力学的这一种假设是可行的。力学的这一种假设是可行的。是对理论计算和实验的有力补充，广泛应用于材料科学、生是对理论计算和实验的有力补充，广泛应用于材料科学、生物物理和药物设计等。经典物物理和药物设计等。经典MD模拟，其系统规模在一般的计模拟，其系统规模在一般的计算机上也可达到数万个原子，模拟时间为纳秒量级。
4、算机上也可达到数万个原子，模拟时间为纳秒量级。分子动力学模拟？分子动力学模拟？Molecular Dynamics3.1分子动力学的基本原理及特点分子动力学模拟的分子动力学模拟的基本原理基本原理：建立一个粒子系统来模拟所研究的微观现象，建立一个粒子系统来模拟所研究的微观现象，系统中各粒子之间的相互作用根据量子力学来系统中各粒子之间的相互作用根据量子力学来确定。确定。对符合经典牛顿力学规律的大量粒子系统，通对符合经典牛顿力学规律的大量粒子系统，通过粒子动力学方程组的数值求解，决定各粒子过粒子动力学方程组的数值求解，决定各粒子在相空间的运动规律和轨迹。在相空间的运动规律和轨迹。然后按照统计物理原理
5、得出该系统相应的宏观然后按照统计物理原理得出该系统相应的宏观物理特性。物理特性。在分子动力学模拟中，通常假设分子的动力学行为在分子动力学模拟中，通常假设分子的动力学行为遵循经典运动方程，这在典型分子的遵循经典运动方程，这在典型分子的De Broglie 波长波长远小于下例公式表示的分子平均距离的平移运动中远小于下例公式表示的分子平均距离的平移运动中确实是一种满意的近似。确实是一种满意的近似。3121)/(3NVTmkhB（3-1）h:普朗克常数；普朗克常数；m:粒子质量；粒子质量；V:体积；体积；N:颗粒数。颗粒数。一、研究内容一、研究内容：l 分子动力学是在原子、分子水平上求解多体问分子动力
6、学是在原子、分子水平上求解多体问题的重要的计算机模拟方法，可以预测纳米尺题的重要的计算机模拟方法，可以预测纳米尺度上的材料动力学特性。度上的材料动力学特性。l 通过求解所有粒子的运动方程，分子动力学方通过求解所有粒子的运动方程，分子动力学方法可以用于模拟与原子运动路径相关的基本过法可以用于模拟与原子运动路径相关的基本过程。程。l 在分子动力学中，粒子的运动行为是通过经典在分子动力学中，粒子的运动行为是通过经典的的NewtonNewton运动方程所描述。运动方程所描述。分子动力学理论特点：分子动力学理论特点：能有效进行微观动态过程分析；能有效进行微观动态过程分析；计算的空间尺度可达到计算的空间尺
7、度可达到1010-9-9 m m；计算的时间尺度可达到计算的时间尺度可达到1010-15-15 s s；可以通过分子动力学理论可以通过分子动力学理论建立起宏观与微观之建立起宏观与微观之间的桥梁间的桥梁。仍存在的主要问题：v 真实所含微观粒子数量巨大，一般高达真实所含微观粒子数量巨大，一般高达10102323数量级；数量级；v 对具有枝链和环状结构等结构形式复杂的柔性分子的模拟计算还对具有枝链和环状结构等结构形式复杂的柔性分子的模拟计算还很困难。很困难。v 符合热力学极限的宏观系统由几千万亿亿个分子或原子组成（一符合热力学极限的宏观系统由几千万亿亿个分子或原子组成（一般高达般高达10102323
8、数量级），所含微观粒子数量巨大。实际计算中，分子数量级），所含微观粒子数量巨大。实际计算中，分子动力学方法受到有限观测时间和有限系统尺寸的限制。由于受计动力学方法受到有限观测时间和有限系统尺寸的限制。由于受计算机运算速度和内存空间的限制，计算机模拟允许的微观体系尺算机运算速度和内存空间的限制，计算机模拟允许的微观体系尺寸要比热力学极限小很多。寸要比热力学极限小很多。一粒黄豆那么大的水滴也含有约一粒黄豆那么大的水滴也含有约10102020个水分子，这数目对于目前计个水分子，这数目对于目前计算机来说实在是太大了，用当前最快速的计算机模拟它，恐怕也算机来说实在是太大了，用当前最快速的计算机模拟它，恐
9、怕也要运行几个世纪。计算机运算能力对分子数目的限制，出现的问要运行几个世纪。计算机运算能力对分子数目的限制，出现的问题是小样本系统的表面效应会掩盖其体效应，小系统的模拟不能题是小样本系统的表面效应会掩盖其体效应，小系统的模拟不能完全反映真实系统的性质和行为。完全反映真实系统的性质和行为。克服这一困难的办法：克服这一困难的办法：是对所选定的模拟单元施加周期性边界条件。是对所选定的模拟单元施加周期性边界条件。图中标有阴影的单元是选定体积为图中标有阴影的单元是选定体积为V V（边长为（边长为L)L)的立方体、内含的立方体、内含N N个粒个粒子的模拟单元。子的模拟单元。何谓周期性边界？何谓周期性边界？
10、就是在它周围想象存在着就是在它周围想象存在着无穷多个与模拟单元完全相同的单元，它们无穷多个与模拟单元完全相同的单元，它们像晶体元胞一样充满整个空间。每个单元内像晶体元胞一样充满整个空间。每个单元内部有数量相等、分布也相同的粒子，且相应部有数量相等、分布也相同的粒子，且相应的粒子具有相同的速度。上述周期性的安排，的粒子具有相同的速度。上述周期性的安排，使得一个粒子从单元的一面离开，必从相对使得一个粒子从单元的一面离开，必从相对的一面以同样的速度进入该单元，从而维持的一面以同样的速度进入该单元，从而维持各单元内粒子数不变。周期性边界条件的引各单元内粒子数不变。周期性边界条件的引入，使得每个小单元的
11、情况完全等价，但我入，使得每个小单元的情况完全等价，但我们面对的仍是一个无限系统，计算中只需处们面对的仍是一个无限系统，计算中只需处理和存储一个小单元的数据就够了，而不必理和存储一个小单元的数据就够了，而不必计算周期性边界条件带来的这个无限系统。计算周期性边界条件带来的这个无限系统。因此，因此，周期性边界条件的引入，使模拟运算周期性边界条件的引入，使模拟运算工作摆脱了巨大分子数的困境，并成功消除工作摆脱了巨大分子数的困境，并成功消除了为减小粒子数而带来的的有限尺寸效应了为减小粒子数而带来的的有限尺寸效应。当然边界条件的引入又会影响体系的某些性当然边界条件的引入又会影响体系的某些性质，但系统的有
12、限尺寸也有其好处，它使我质，但系统的有限尺寸也有其好处，它使我们有可能计算二阶的热力学性质。们有可能计算二阶的热力学性质。周期性边界条件周期性边界条件2.分子动力学模拟的基本步骤 在计算机上对分子系统的在计算机上对分子系统的MD模拟的实际步骤可以分为：模拟的实际步骤可以分为：首先是设定模拟所采用的模型；首先是设定模拟所采用的模型；包括几何建模，物理建模，化学建模，力学建模。包括几何建模，物理建模，化学建模，力学建模。给定初始条件；给定初始条件；初始条件的设定，这里要从微观和宏观两个方面进行考虑。初始条件的设定，这里要从微观和宏观两个方面进行考虑。计算计算 是体现所谓分子动力学特点的地方，包括对
13、运动方程的积分的有效算是体现所谓分子动力学特点的地方，包括对运动方程的积分的有效算法。对实际的过程的模拟算法，关键是分清楚平衡和非平衡，静态和法。对实际的过程的模拟算法，关键是分清楚平衡和非平衡，静态和动态以及准静态情况。动态以及准静态情况。分析分析 需要从以上的计算的结果中提取年需要的特征，说明你的问题的实质需要从以上的计算的结果中提取年需要的特征，说明你的问题的实质和结果，因此关键是统计、平均、定义、计算，比如温度、体积、压和结果，因此关键是统计、平均、定义、计算，比如温度、体积、压力、应力等宏观量和微观过程量是怎么联系的。力、应力等宏观量和微观过程量是怎么联系的。模拟模型的设定 假定两个
14、分子间的相互作用势为硬球势，其势函数表示为常用Lennard-Jones势函数：3.2 平衡态分子动力学模拟理论平衡态分子动力学模拟理论 核心核心是要计算系统中所有粒子的运动规律和是要计算系统中所有粒子的运动规律和轨迹。轨迹。计算中需要假设：计算中需要假设：l 所有粒子的运动都遵循经典牛顿力学规律；所有粒子的运动都遵循经典牛顿力学规律；l 粒子之间的相互作用满足叠加原理；粒子之间的相互作用满足叠加原理；3.2 平衡态分子动力学模拟 经典经典MD模拟方法的应用当中，存在着对两种系统模拟方法的应用当中，存在着对两种系统状态的状态的MD模拟。模拟。一种是对平衡态的一种是对平衡态的MD模拟，另一类是对
15、非平衡态模拟，另一类是对非平衡态的的MD模拟。模拟。对平衡态系综对平衡态系综MD模拟又可以分为如下类型：微正模拟又可以分为如下类型：微正则系综的则系综的MD（NVE）模拟，正则系综的）模拟，正则系综的MD（NVT）模拟，等温等压系综）模拟，等温等压系综MD（NPT）模拟和）模拟和等焓等压系综等焓等压系综MD（NPH）模拟等。）模拟等。在统计物理中的正则系综模拟是针对一个粒在统计物理中的正则系综模拟是针对一个粒子数子数N、体积、体积V、温度、温度T和总动量（和总动量（P=0）为守）为守恒量的系综（恒量的系综（NVT）。这种情况就如同一个）。这种情况就如同一个系统置于热浴之中，此时系统的能量可能有
16、系统置于热浴之中，此时系统的能量可能有涨落，但系统温度则已经保持恒定。在正则涨落，但系统温度则已经保持恒定。在正则系综的系综的MD模拟中施加的约束与微正则系综中模拟中施加的约束与微正则系综中的不一样。正则系综的不一样。正则系综MD方法是在运动方程组方法是在运动方程组上加上动能恒定（即温度恒定）的约束，而上加上动能恒定（即温度恒定）的约束，而不是像微正则系综的不是像微正则系综的MD模拟中对运动方程加模拟中对运动方程加上能量恒定的约束。上能量恒定的约束。1、基本运动方程2.典型势函数3.2.2 典型势函数(1).偶对势函数在分子动力学模拟计算中，势函数的选择极为重要，它决定了计算工作量、计算模型与
17、真实系统的近似程度。根据前述2个假设，通常选择偶对势函数uij(rij)即jiijijirurv)()(（3-7）式（3-6）变为ijijrrr（3-8）iijiFFijijijjiijrruFF)(6124)(ijijijijrrru(2).Lennard-Jones(LJ)势（3-9）（3-10）在处理流体润滑问题时，根据润滑流体分子相互作用的特性，常应用Lennard-Jones势函数：（3-11）(3).Stillinger-Weber作用势该作用势已被该作用势已被Kluge、Volz等人证实能很好地描述硅材料。硅等人证实能很好地描述硅材料。硅材料是现代电子技术基础。因此，该作用势对微
18、器件中的一材料是现代电子技术基础。因此，该作用势对微器件中的一些热问题分析将具有特别重要的意义。些热问题分析将具有特别重要的意义。)()()(2rhrBArUqprbcaabccabbcacabhhhrrrU),(323/1)cos()(cabbcrabrcabrhrahhStilliger-Weber作用势包括二体及三体相：作用势包括二体及三体相：（3-12）（3-13）（3-14）4.硬球势硬球势202)()(ijijijijttvr22/12)(ijijijijijijacabbt（3-15）（3-16）硬球的作用势在硬球直径 r=处呈尖锐的不连续性，超出这一范围，则无相互作用，而对于
19、r 范围，则作用势无穷大。相应的势函数公式为 0ijijru上述几种势函数与实际情况虽然有一定差别，但它们简单实用，在流体问题的模拟计算中已得到一定的应用。3.2.3 求解方法 根据给定的边界条件和初始值，求解方程组根据给定的边界条件和初始值，求解方程组（3-5）至（）至（3-10）就可求得各个粒子在）就可求得各个粒子在6n维维相空间的运动轨迹，然后再由统计理论得到该相空间的运动轨迹，然后再由统计理论得到该物理系统的客观特性。求解物理系统的客观特性。求解n个粒子的方程组个粒子的方程组计算工作量大，为了简化计算，一般采用每步计算工作量大，为了简化计算，一般采用每步只计算一次力的单步算法，这样可以
20、使计算省只计算一次力的单步算法，这样可以使计算省时、稳定和高精度。设粒子系统具有下列形式时、稳定和高精度。设粒子系统具有下列形式的初值的初值)(irfr)(00tiirr（3-17）在分子动力学模拟中，常采用在分子动力学模拟中，常采用Gear、Verlet、截断半径截断半径法法等几类算法。等几类算法。分子动力学模拟计算 求解方法 以一个一维谐振子为例，其经典哈密顿量为 一步法：有了初始条件x(0)、p(0)，就可以一步一步地使用前一时刻的坐标、动量值确定下一时刻的坐标、动量值。两步法 泰勒展开式的一般形式：二步法：所得的坐标和动量的递推公式分别为：截断半径法是为了克服计算繁杂，耗费机时而引入的
21、一种处理方法。其特点是预先选定一个截断半径rc只计算以截断半径为球体内的粒子间的作用力，而与粒子之间的距离超过截断半径时，则不考虑它们的作用。同时，系统中粒子数量很大时，还可以利用邻域列表法判断粒子的分步情况，进一步节省机时。Verlet算法的邻域列表法如图所示r1rc1234567rc:截断半径；截断半径；r1:邻域半径邻域半径如粒子1受的力只计算粒子2，3，4的作用力之和；而与粒子之间的距离超过截断半径时，则不考虑他们的作用（如5，6对1）。计算过程中每隔一定步数需要更新邻域表：即模拟计算的开始阶段，固定更新邻域表的步数为10到20步；此后计算采用自动调整法更新后rc内粒子与1间距的变化；
22、当粒子5，6与1的距离 rc或粒子7与1的距离 r1时，需要更新邻域表。3.2.4 边界条件与初值 由于计算机的运算能力有限，模拟计算系由于计算机的运算能力有限，模拟计算系统的粒子数不可能很大，这就会导致模拟统的粒子数不可能很大，这就会导致模拟系统粒子数少于真实系统的所谓系统粒子数少于真实系统的所谓“尺寸效尺寸效应应”问题。问题。为了选取尽可能多的粒子数而又不至于使为了选取尽可能多的粒子数而又不至于使计算工作量过分庞大，在统计物理中，对计算工作量过分庞大，在统计物理中，对于平衡态分子动力学模拟计算引用三维周于平衡态分子动力学模拟计算引用三维周期边界条件。期边界条件。图3-2.周期边界条件（阴影
23、区为采样区）模拟系统粒子数与真实模拟系统粒子数与真实系统存在系统存在“尺寸效应尺寸效应”问题。问题。平衡态分子动力学模拟平衡态分子动力学模拟计算采用三维周期边界计算采用三维周期边界条件，选取粒子数尽可条件，选取粒子数尽可能多，计算量不过分庞能多，计算量不过分庞大。大。采样区的周围是其镜像，采样区的周围是其镜像，根据采样区的边界条件根据采样区的边界条件计算采样区内的粒子运计算采样区内的粒子运动。动。3.3 非平衡态分子动力学模拟 基本方法是通过对模拟系统施加扰动，利用系统基本方法是通过对模拟系统施加扰动，利用系统的自然波动出现的响应，求解时间相关函数的粘的自然波动出现的响应，求解时间相关函数的粘
24、度等输运特性。这种模拟属于非线性相关系统的度等输运特性。这种模拟属于非线性相关系统的模拟，广泛应用于液体剪切流动模拟中。模拟，广泛应用于液体剪切流动模拟中。粘度是润滑薄膜摩擦行为中最主要的特性。计算粘度是润滑薄膜摩擦行为中最主要的特性。计算剪切粘度的计算方法有剪切粘度的计算方法有均匀非平衡态模拟均匀非平衡态模拟和和非均非均匀平衡态模拟匀平衡态模拟。二者的主要区别在于引起系统扰。二者的主要区别在于引起系统扰动的方法不同，从而周期边界条件不同。动的方法不同，从而周期边界条件不同。1.均匀非平衡态模拟均匀非平衡态模拟均匀非平衡态模拟是在方程中加入扰动项，求解响应之后均匀非平衡态模拟是在方程中加入扰动
25、项，求解响应之后确定粘度，方程的一般解如下：确定粘度，方程的一般解如下：)(tpmpq)(tqFp（3-26）（3-27）Evans等人近年来作了重大发展，采用三维周期边界条件，等人近年来作了重大发展，采用三维周期边界条件，提出了适用于大扰动，非牛顿流动特点的提出了适用于大扰动，非牛顿流动特点的Sllod算法，在宏算法，在宏观剪切流动模拟计算中取得了良好效果。观剪切流动模拟计算中取得了良好效果。Fig.Lee-Eduards 边界条件r图中 r为剪切应变率，L采样尺寸，t为时间，影线部分为采样区。rrvrmpriiiiiiiiapvpFpSllod算法求解均匀非平衡态系统（对等温剪切流动系统）
26、模拟的应用实例方程：（3-28）（3-29）式中a为等温系数，其计算式表示为：nijiijixiippppFp12.非均匀非平衡态分子动力学模拟非均匀非平衡态分子动力学模拟非均匀非平衡态分子动力学模拟是通过固体壁面扰动诱发的，采用下图所示的二维周期边界条件。流体的剪切流动仅在壁面运动产生，上下2个壁面保持恒温，周期边界条件仅在x、y方向成立。因壁面只影响到邻近壁面流体分子的运动，故为非均匀非平衡态系统。Fig.二维边界条件5.量子分子动力学方法 如激光加工、材料制造中的激光控制及激光冷却的进展，使得有必要分析和考虑带电粒子(电子、核)与光电场之间的相互作用中的非平衡现象，这些相互作用具有量子及
27、热能水平上的特征；假设原子系统简单地由离子及电子组成，而离子是一个含一定数量电离能的封闭壳(如下图)，离子的动能可通过分于动力学方法预测出，而电子的波函数则由时间依赖型Schrodinge方程求得，这种简化方法通常称为Borm-Oppenheimer近似。利用这一量子分子动力学方法，可以定性地研究光辐照与动能改变及原于分裂之间的关系，从而较好地理解光与物质之间相互作用的基本机制。思考题 分子动力学模拟适用于对什么体系进行描分子动力学模拟适用于对什么体系进行描述？述？分子动力学的基本原理是什么？分子动力分子动力学的基本原理是什么？分子动力学理论有哪些特点？存在的主要问题是什学理论有哪些特点？存在的主要问题是什么？么？
