动力学理论+++
2.3 多柔体系统动力学建模
2.3.1 柔性体上点的位置向量、速度和加速度 1．柔性体系统中的坐标系
图2.5 柔性体上节点P 的位置
柔性体系统中的坐标系如图2.5所示，包括惯性坐标系（r
e ）和动坐标系（b
e ）。前者不随时间而变化，后者是建立在柔性体上，用于描述柔性体的运动。动坐标系可以相对惯性坐标系进行有限的移动和转动。动坐标系在惯性坐标系中的坐标（移动、转动）称为参考坐标。
与刚体不同，柔性体是变形体，体内各点的相对位置时时刻刻都在变化，只靠动坐标系不能准确描述该柔性体在惯性坐标系中的位置，因此，引入弹性坐标来描述柔性体上各点相对动坐标系统的变形。这样柔性体上任一点的运动就是动坐标系的“刚性”运动与弹性变形的合成运动。由于柔体上各点之间有相对运动，所以动坐标系的选择不是采用连体坐标系，而需要采用随着柔性体形变而变化的坐标系，即“浮动坐标系”。
在研究多柔体系统时，合适的坐标系是非常重要的。在确定浮动坐标系时有两点准则：1、便于方程建立求解；2、柔性体刚体运动与变形运动的耦合尽量小。目前常见的浮动坐标系大致有如下5种，局部附着框架、中心惯性主轴框架、蒂斯拉德框架、巴克凯恩斯框架以及刚体模态框架。采用何种需因实际情况而定。 2．柔性体上任意点的位置向量、速度和加速度
在分析刚体平面运动的时候，把复杂的刚体平面运动分解为几种简单的运动。在对柔性体的运动，尤其是在小变形的情况下，也可以采用类似的方法。如某柔性体从位置L1运动到位置L2，其间运动可以分解为：刚性移动－>刚性转动－>变形运动。对于柔性体上任意一点P ，其位置向量为：
r
e b
e P
i
B i
C i
r P
i r P
i s 0
P i s P
i u
0()
p p r r A s u =++
（2-158）
r 为P 点在惯性坐标系中的向量；0r 为浮动坐标系原点在惯性坐标系中的向
量；根据式2.2- 14，A 为方向余弦矩阵；p s 为柔性体未变形时P 点在浮动坐标系中的向量;p u 为相对变形向量，p u 可以用不同的方法离散化，与讨论平面问题相同，对于点P ，该单元的变形采用模态坐标来描述，有：
p p f
u q =Φ
（2-159）
式（2-159）中，p Φ为点P 满足里兹基向量要求的假设变形模态矩阵，q 为变形的广义坐标。
柔性体上任一点的速度向量及加速度向量可以对式求对时间一阶导数和二阶导数得到：
0()P p p p f
A s u A q =+++Φ&&&&r r
（2-160）
0()2p p p p f p f
r r A s u A q A q =+++Φ+Φ&&&&&&&&&&
（2-161）
2．多柔体系统的能量
（1）动能和质量矩阵
考虑节点P 变形前后的位置、方向和模态，柔性体的广义坐标可以表示为：
[(1,,)][]T T i x
y z q i M ψ
θφ
===ξr ψ
q L （2-176）
速度表达式（2-160）在系统广义坐标式（2-176）的时间导数ξ&中表示为：
G ()
r e B r B
()
b e P
s P
P
u P '
图2.6 柔性
[()]P p p P v I
A s u
B A ξ
=-+Φ&% （2-177） 柔性体的动能为：
1122T T GBT GB
P P P P P P V P
T v vdV m v v I ρωω=
≈+∑?
（2-178）
其中，P m 和P I 分别为节点P 的节点质量和节点惯性张量；GB
P
P B ωψ=&，为点B 相对于全局坐标基的角速度在局部坐标基中的斜方阵表示。将式（2.3-20）
和关系式P P B ωψ=&代入式（2-178），得到动能的广义表达式：
1()2
T T M ξξξ=&&
（2-179）
上式中的质量矩阵()M ξ为33?维的方阵，表示为：
()tt
tr tm T
tr rr
rm T T
tm
rm
mm M M M M M M M M M M ξ????=?
?????
（2-180）
其中下标,,t r m 分别表示平动、旋转和模态自由度。质量矩阵的六个独立分量分别表示为：
12337889
456
[][[]][]
tt tr j j tm T T rr j j j ij i j T rm j j mm M I
M A q B
M A M B q q q B
M B q M τττττττττττ?
=?=-+??=??=-+-??=+?=??
%
（2-181）
其中9个惯性时不变矩阵列表如下：
式(2.3-24)中可以明显看出质量矩阵与模态坐标显性相关，而且由于引入转换矩阵A 和B ，质量矩阵也与系统的方向坐标显性相关。质量矩阵中的9个惯性时不变矩阵19~ττ可通过计算有限元模型的N 个节点信息在预处理过程中一次性得到，从而简化运动微分方程的求解。节点信息包括：每个节点的质量P m ，未变形时的位置矢量P s 以及模态矩阵P Φ。
（2）势能和刚度矩阵
势能一般分为重力势能和弹性势能两部分，可用下列二次项表示：
1
()2
T g W W K ξξξ
=+ （2-182）
在弹性势能中，K 是对应于模态坐标q 的结构部件的广义刚度矩阵，通常为
常量。重力势能g W
表示为：
[()]T g P B P P W
W
W r gdW r A s q gdW
ρρ=?=++Φ??
（2-183）
其中g 表示重力加速度矢量，重力g f 可对g W 求导得；
[][()][]W g T g P P W T P W dW W dW dW ρρρ????
????==+??
????
??????
???g
A f s Φq g ξξA Φg
（2-184）
（3）能量损失和阻尼矩阵
阻尼力依赖于广义模态速度并可以从下列二次项中推导得出：
12
T
q Dq Γ=&& （2-185）
上式称为Rayleigh 能量损耗函数。矩阵D 包含阻尼系数ij d
，它是常值对称
阵。当引入正交模态振型时，阻尼矩阵可用对角线为模态阻尼率i c
的对角阵来表示。对于每一个正交模态的阻尼率都可以取不同值，而且还能以该模态的临界
阻尼cr i c 的比值形式给出。
3．多柔体动力学方程
柔性体的运动方程从下列拉格朗日方程导出：
()[]T
d L L dt λ???Γ??-++-=????????
ΨQ 0ξξξξΨ=0
&& （2-186）
其中，Ψ为约束方程；
λ为对应于约束方程的拉氏乘子；
ξ为如式（2.3-19）定义的广义坐标； Q 为投影到ξ上的广义力；
L 为拉格朗日项，定义为L T W =-，T 和W 分别表示动能和势能； Γ表示能量损耗函数。
将求得的,,T W Γ代入式（2.3-29），得到最终的运动微分方程为：
1[][]2T T g M M M K f D Q ξξξξξξλξξ
??ψ+-++++=??&&&&&&& （2-187） 其中，ξ,ξ,ξ&&&为柔性体的广义坐标及其时间导数；
,M M
&为柔性体的质量矩阵及其对时间的导数；
??M
ξ
为质量矩阵对柔性体广义坐标的偏导数，它是一个(6)(6)(6)M M M +?+?+维张量，M 为模态数。
4.1.3 多刚体系统运动学
对于多体系统的运动学分析，传统的理论力学是以刚体位置、速度和加速度的微分关系以及矢量合成原理为基础进行分析的，而计算多体系统动力学中的运动学分析则是以系统中连接物体与物体的运动副为出发点，所进行的位置、速度和加速度分析都是基于与运动副对应的约束方程来进行的。
基于约束的多体系统运动学，首先寻求与系统中运动副等价的位置约束代数方程，再由位置约束方程的导数得到速度、加速度的约束代数方程，对这些约束方程进行数值求解，可得到广义位置坐标及相应的速度和加速度坐标，最后根据坐标变换就可以由系统广义坐标及相应导数得到系统中任何一点的位置、速度和加速度。
由于机械系统在二维空间运动时，广义坐标、约束方程、问题规模以及问题求解都相对简单，故本节先讨论二维多体系统运动学以解释多体系统运动学基本理论，在此基础上再给出三维多体系统的运动学方程。
1．约束方程（位置方程）
设一个平面机构由nb 个刚性构件组成。在机构所在平面上建立一个全局坐标系xoy ，机构在该坐标系中运动；再为机构上每个构件i 建立各自的连体坐标系i i i y o x '''，可由连体坐标系的运动确定构件的运动。选定构件i 连体坐标系原点i o '的全局坐标T i i i y x r ],[=和连体坐标系相对于全局坐标系的转角i φ组成构件i 的笛卡尔广义坐标矢量T i i i i y x q ],,[φ≡，如图2.2所示。由nb 个刚性构件组成的系统
的广义坐标数nb nc ?=3，则系统广义坐标矢量可表示为T
T nb T T q q q q ],...,,[21
=。
y
x o
x
i
i
φ
图2.2 平面笛卡尔广义坐标
一个实际的机械系统，系统中构件与支架或构件与构件之间存在运动副的联接，这些运动副可以用系统广义坐标表示为代数方程。设表示运动副的约束方程数为nh ，则用系统广义坐标矢量表示的运动学约束方程组为：
0)](),...,(),([)(21=ΦΦΦ=ΦT K
nh K K K q q q q (2-8)
这里给出的是定常完整约束情况。如果约束方程与时间相关，则自变量中显含时间项，这种约束被称为非定常约束；更一般的约束方程含有不可积速度项的不等式或关系式，这种约束称为非完整约束。一般的运动学约束是定常完整约束。
对于一个有nc 个广义坐标和nh 个约束方程的机械系统，若nh nc >，且这nh 个约束方程是独立、相容的，则系统自由度nh nc DOF -=。为使系统具有确定运动，可以有二种方法：
（1）为系统添加与系统自由度DOF 相等的附加驱动约束； （2）对系统施加力的作用。
在（1）情况下，系统实际自由度为零，被称为是在运动学上确定的，在此情况下求解系统运动过程中的位置、速度和加速度的分析是运动学分析，运动学分析本身不涉及作用力或反作用力问题。但是对于运动学上确定的系统，可以求解系统中约束反力，即已知运动求作用力，这是动力学逆问题。
在（2）情况下，系统有着大于零的自由度，但是在外力作用力，对于具有确定构型和特定初始条件的系统，其动力学响应是确定的，这种情况下求解系统运动过程中的位置、速度和加速度的分析，称为动力学分析。在这种情况下，特殊地，如果外力与时间无关，可以求解系统的静平衡位置，这就是静平衡分析问题。
考虑运动学分析，为使系统具有确定运动，也就是要使系统实际自由度为零，为系统施加等于自由度（nh nc -）的驱动约束：
0),(=Φt q D (2-9)
在一般情况下，驱动约束是系统广义坐标和时间的函数。驱动约束在其集合内部及其与运动学约束合集中必须是独立和相容的，在这种条件下，驱动系统运动学上是确定的，将作确定运动。
由式（2-8）表示的系统运动学约束和式（2-9）表示的驱动约束组合成系统所受的全部约束：
0),(),(),(=???
???ΦΦ=Φt q t q t q D K (2-10)
式（2-10）为nc 个广义坐标的nc 个非线性方程组，其构成了系统位置方程。求解式（2-10），就可得到系统在任意时刻的广义坐标位置)(t q 。
2．速度和加速度方程
对式(2-10)运用链式微分法则求导，得到速度方程：
0),(),(),,(=Φ+Φ=Φt q q t q t q q t q &&&
(2-11)
若令),(t q t Φ-=υ，则速度方程为
0),(),,(=-Φ=Φυq t q t q q q &&& (2-12) 如果
q
Φ是非奇异的，可以求解式(2-12)得到各离散时刻的广义坐标速度q &。
对式(2-11)运用链式微分法则求导，可得加速度方程
0),(),(2)),((),(),,,(=Φ+Φ+Φ+Φ=Φt q q t q q q t q q t q t q q q tt qt q q q &&&&&&&&&& (2-13)
若令tt qt q q q q q
Φ-Φ-Φ-=&&&2)(η，则加速度方程为 0
),,(),(),,,(=-Φ=Φt q q q t q t q q q q &&&&&&&&η
(2-14)
如果q Φ是非奇异的，可以求解式(2-14)得到各离散时刻的广义坐标加速度
q &&。在速度方程(2-12)和加速度方程(2-14)中出现的矩阵q Φ，称为雅可比矩阵，
雅可比矩阵是约束多体系统运动学和动力学分析中最重要的矩阵。如果Φ的维数为m ，q 维数为n ，那么q Φ维数为n m ?矩阵，其定义为j i j i q q ?Φ?=Φ),()(。这里q Φ为nc nc ?的方阵。
对式(2-12)中的υ和式(2-14)中的η进行计算时，会涉及到二阶导数，在实
际的数值求解中，并不是实时地调用求导算法来进行计算，而是先根据具体的约束类型，导出二阶导数以及雅可比矩阵的表示式，在计算中只需代入基本的数据即可。
3．坐标变换与任意点运动
在确定系统中构件上任意点的运动时，常要求将构件上点从连体坐标系变换到全局坐标系中，现讨论连体坐标系与全局坐标系的坐标变换及构件上任意点运动。
设矢量s ρ在全局坐标系xoy 和某连体坐标系i i i y o x '''中分别表示为：
?????='=''T
y x T y x s s s s s s ]
,[]
,[ (2-15) 若任意点P 在全局坐标系xoy 和连体坐标系i i i y o x '
''中坐标如图2-3所示，则
存在如下坐标变换关系：
P P P s A r s r r '+=+= (2-16)
其中，P
r 为点P 在全局坐标系中的坐标，r 为连体坐标系原点O '在全局坐
标系中的坐标，P s 为矢量s ρ在全局坐标系中坐标，P
s '为矢量s ρ在连体坐标系中的坐标，A 为旋转变换矩阵，其形式为：
??
?
?
??-==φφφφ
φcos sin sin cos )(A A (2-17) A 对时间的导数为：
B A d d A φφφφφφφφ&&&&≡??????---==sin cos cos sin (2-18) 根据式(2-16)，我们可以得到以连体坐标系表示的构件上的任一点的全局坐标。
x
y
r
P r P
s
P
o 'r
P
r P
s P
o
y
x
x
z
y
z
二维空间坐标变换 三维空间坐标变换 式(2.2-9)对时间求导数，可得任意点的速度变换公式： P P P s B r s A r r '+='+=φ&&&&& (2-19) 式(2.2-12)对时间求导数，可得任意点的加速度变换公式： P P P P P s A s B r s B s B r r '-'+='+'+=2φφφφ&&&&&&&&&&&&& (2-20) 对于一个平面机构来说，进行运动学分析时，先是选定最大集的广义坐标，再分别根据式(2-10)、(2-12)和(2-14）求解机构在各离散时刻的广义坐标位置、广义坐标速度和广义坐标加速度。对于任意一个由连体坐标系确定的构件上的点，可以根据式(2-16)、(2-19)和(2-20)求解其位置、速度和加速度。
《装备制造技术》 2007年第 10期 收稿日期 :2007-08-21 作者简介 :申海兰 , 24岁 , 女 , 河北人 , 在读研究生 , 研究方向为微机电系统。 分子动力学模拟方法概述 申海兰 , 赵靖松 (西安电子科技大学机电工程学院 , 陕西西安 710071 摘要 :介绍了分子动力学模拟的基本原理及常用的原子间相互作用势 , 如Lennard-Jones 势 ; 论述了几种常用的有限差分算法 , 如 Verlet 算法 ; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。关键词 :分子动力学模拟 ; 原子间相互作用势 ; 有限差分算法 ; 系综中图分类号 :O3 文献标识码 :A 文章编号 :1672-545X(200710-0029-02 从统计物理学中衍生出来的分子动力学模拟方法 (molec- ular dynamics simulation , M DS , 实践证明是一种描述纳米科技 研究对象的有效方法 , 得到越来越广泛的重视。所谓分子动力学模拟 , 是指对于原子核和电子所构成的多体系统 , 用计算机模拟原子核的运动过程 , 从而计算系统的结构和性质 , 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动 [1]。它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段 , 称之为“计算机实验” 手段 [2], 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
根据模拟对象的不同 , 将它分为平衡态分子动力学模拟 (EM DS (和非平衡态分子动力学模拟 (NEM DS 。其中 , EM DS 是分子动力学模拟的基础 ; NEM DS 适用于非线性响应系统的模拟 [3]。下面主要介绍 EM DS 。 1分子动力学方法的基本原理 计算中根据以下基本假设 [4]: (1 所有粒子的运动都遵循经典牛顿力学规律。 (2 粒子之间的相互作用满足叠加原理。 显然这两条忽略了量子效应和多体作用 , 与真实物理系统存在一定差别 , 仍然属于近似计算。 假设 N 为模拟系统的原子数 , 第 i 个原子的质量为 m i , 位置坐标向量为 r i , 速度为 v i =r ? i , 加速度为 a i =r ?? i , 受到的作用力为 F i , 原子 i 与原子 j 之间距离为 r ij =r i -r j , 原子 j 对原子 i 的作用力为 f ij , 原子 i 和原子 j 相互作用势能为 ! (r ij , 系统总的势能为 U (r 1, r 2, K r N = N i =1! j ≠ i ! " (r ij , 所有的物理量都是随时 间变化的 , 即 A=A (t , 控制方程如下 : m i r ?? i =F i =j ≠ i
分子动力学的模拟过程 分子动力学模拟作为一种应用广泛的模拟计算方法有其自身特定的模拟步骤,程序流程也相对固定。本节主要就分子动力学的模拟步骤和计算程序流程做一些简单介绍。 1. 分子动力学模拟步驟 分子动力学模拟是一种在微观尺度上进行的数值模拟方法。这种方法既可以得到一些使用传统方法,热力学分析法等无法获得的微观信息,又能够将实际实验研究中遇到的不利影响因素回避掉,从而达到实验研宄难以实现的控制条件。 分子动力学模拟的步骤为: (1)选取所要研究的系统并建立适当的模拟模型。 (2)设定模拟区域的边界条件,选取粒子间作用势模型。 (3)设定系统所有粒子的初始位置和初始速度。 (4)计算粒子间的相互作用力和势能,以及各个粒子的位置和速度。 (5)待体系达到平衡,统计获得体系的宏观特性。 分子动力学模拟的主要对象就是将实际物理模型抽象后的物理系统模型。因此,物理建模也是分子动力学模拟的一个重要的环节。而对于分子动力学模拟,主要还是势函数的选取,势函数是分子动力学模拟计算的核心。这是因为分子动力学模拟主要是计算分子间作用力,计算粒子的势能、位置及速度都离不开势函数的作用。系统中粒子初始位置的设定最好与实际模拟模型相符,这样可以使系统尽快达到平衡。另外,粒子的初始速度也最好与实际系统中分子的速度相当,这样可以减少计算机的模拟时间。 要想求解粒子的运动状态就必须把运动方程离散化,离散化的方法有经典Verlet算法、蛙跳算法(Leap-frog)、速度Veriet算法、Gear预估-校正法等。这些算法有其各自的优势,选取时可按照计算要求选择最合适的算法。 统计系统各物理量时,便又涉及到系统是选取了什么系综。只有知道了模拟系统采用的系综才能釆用相对应的统计方法更加准确,有效地进行统计计算,减少信息损失。 2. 分子动力学模拟程序流程 具体到分子动力学模拟程序的具体流程,主要包括: (1)设定和模拟相关的参数。 (2)模拟体系初始化。 (3)计算粒子间的作用力。 (4)求解运动方程。 (5)循环计算,待稳定后输出结果。 分子动力学模拟程序流程图如2.3所示。
分子动力学 分子动力学方法是一种计算机模拟实验方法，是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。它是对理论计算和实验的有力补充。 分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和Born-Oppenheimer绝热近似下，分子动力学的这一种假设是可行的[1]。所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。要进行分子动力学模拟就必须知道原子间的相互作用势。在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。然而采用经验势必然丢失了局域电子结构之间存在的强相关作用信息，即不能得到原子动力学过程中的电子性质[1]。 事实上，分子动力学就是模拟原子系统的趋衡过程。实际上，分子动力学方法就是确定某一描述与初始条件、边值关系的数值解。我们假定系统经过M步长之后达到稳定，而这
一稳定状态正是我们所求的。 １、分子动力学的算法分析 首先，我们假定我们研究的系统服从 Newton 方程所确定的描述，即： )(1 )(.. t F m t r = （1） 式中r(t)表征原子在t 时刻的位置矢量 F(t)表征原子在t 时刻所受到的力，它与所有原子的位置矢有关 m 表征原子的质量。 如果我们给定初始条件，即方程（１）的定解条件r(0)和v(0)，那么方程（１）的解就可以确定。６０年代中期发展了大量的分子动力学算法，如两步差分算法[2]、预测-校正算法 [3] 、中心差分算法[4]、蛙跳算法[5]等等。为了方便导出它们， 我们以Euler 一步法[6] 来讨论之。我们令)()(.. t r t v =（表征粒子 的速度），则有： ) ()()(1 )()(... . t v t r t F m t r t v === （２） 记??? ? ??????=? ? ? ???=)()(1)()()()(. t v t F m t f t r t v t w （３）
分子动力学方法模拟基本步骤 1.第一步 即模型的设定，也就是势函数的选取。势函数的研究和物理系统上对物质的描述研究息息相关。最早是硬球势，即小于临界值时无穷大，大于等于临界值时为零。常用的是LJ势函数，还有EAM势函数，不同的物质状态描述用不同的势函数。 模型势函数一旦确定，就可以根据物理学规律求得模拟中的守恒量。 2 第二步 给定初始条件。运动方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始条件。如：verlet算法需要两组坐标来启动计算，一组零时刻的坐标，一组是前进一个时间步的坐标或者一组零时刻的速度值。 一般意思上讲系统的初始条件不可能知道，实际上也不需要精确选择代求系统的初始条件，因为模拟实践足够长时，系统就会忘掉初始条件。当然，合理的初始条件可以加快系统趋于平衡的时间和步伐，获得好的精度。 常用的初始条件可以选择为：令初始位置在差分划分网格的格子上，初始速度则从玻尔兹曼分布随机抽样得到；令初始位置随机的偏离差分划分网格的格子上，初始速度为零；令初始位置随机的偏离差分划分网格的格子上，初始速度也是从玻尔兹曼分布随机抽样得到。 第三步 趋于平衡计算。在边界条件和初始条件给定后就可以解运动方程，进行分子动力学模拟。但这样计算出的系统是不会具有所要求的系统的能量，并且这个状态本身也还不是一个平衡态。 为使得系统平衡，模拟中设计一个趋衡过程，即在这个过程中，我们增加或者从系统中移出能量，直到持续给出确定的能量值。我们称这时的系统已经达到平衡。这段达到平衡的时间成为驰豫时间。 分子动力学中，时间步长的大小选择十分重要，决定了模拟所需要的时间。为了减小误差，步长要小，但小了系统模拟的驰豫时间就长了。因此根据经验选择适当的步长。如，对一个具有几百个氩气Ar分子的体系，lj势函数，发现取h为0.01量级，可以得到很好的相图。这里选择的h是没有量纲的，实际上这样选择的h对应的时间在10-14s的量级呢。如果模拟1000步，系统达到平衡，驰豫时间只有10-11s。 第四步 宏观物理量的计算。实际计算宏观的物理量往往是在模拟的最后揭短进行的。它是沿相空间轨迹求平均来计算得到的（时间平均代替系综平均）
第五章 分子动力学 第一节 Verlet 算法 牛顿方程 i i i m f dt r d 2 2 记 N r r r R ,,21 N N m f m f m f G ,,221 1 方程写为 2 d R G dt v v 三点公式 2 42 111122n n n n n n n R R R G R R v v v v v v v r 如果给出初始条件0R 和1R ，可求解方程，但常常给出的初 始条件是00,v R ， 那么 02 0012 G v R R (为什么？ 因为dv G dt r ，所以，0000 ()'(')t v t v dt G t v t G r r v ;， 所以，210000000 '(')R R dt v t G R v G r r r r r ;) 方法的优点： 保持时间反演不变性，即令 n n , 方程形式 不变 （尽管误差会破坏这一对称性）
如果问题与v 无关，计算精度相当高 方法的缺点： n v v 必须用到1n R v （为什么是缺点？） 另一方案 2 221112! ()2 n n n n n n n n R R v G v v G G v v v v v 缺点：失去时间反演不变性 第二节 多体问题的基本方法 （阅读材料） 全同粒子，概率分布为 N r r r W R W 21, 物理量平均值 1i i A A R W R dR dR dr Z Z W R dR v v v v v v v 分子动力学 1 lim dt t A A n 个粒子处于 n r r ,1的分布密度函数 N n n n r d r d R W n N N Z r r r 121!!1, !! n N N 来自N 个粒子中取n 个的组合数 例如：N n 是1
分子动力学模拟 分子动力学就是一门结合物理,数学与化学的综合技术。分子动力学就是一套分子模拟方法,该方法主要就是依靠牛顿力学来模拟分子体系的运动,以在由分子体系的不同状态构成的系统中抽取样本,从而计算体系的构型积分,并以构型积分的结果为基础进一步计算体系的热力学量与其她宏观性质。 这门技术的发展进程就是: 1980年:恒压条件下的动力学方法(Andersenの方法、Parrinello-Rahman法) 1983年:非平衡态动力学方法(Gillan and Dixon) 1984年:恒温条件下的动力学方法(能势‐フーバーの方法) 1985年:第一原理分子动力学法(→カー?パリネロ法) 1991年:巨正则系综的分子动力学方法(Cagin and Pettit)、 最新的巨正则系综,即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触,此时系统不仅同热源有能量交换,而且可以同粒子源有粒子的交换,最后达到平衡,这种系综称巨正则系综。 进行分子动力学模拟的第一步就是确定起始构型,一个能量较低的起始构型就是进行分子模拟的基础,一般分子的其实构型主要就是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度,这一速度就是根据玻尔兹曼分布随机生成,由于速度的分布符合玻尔兹曼统计,因此在这个阶段,体系的温度就是恒定的。另外,在随机生成各个原子的运动速度之后须进行调整,使得体系总体在各个方向上的动量之与为零,即保证体系没有平动位移。 由上一步确定的分子组建平衡相,在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子与分子中的原子开始根据初始速度运动,可以想象其间会发生吸引、排斥乃至碰撞,这时就根据牛顿力学与预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算,在这个过程中,体系总能量不变,但分子内部势能与动能不断相互转化,从而体系的温度也不断变化,在整个过程中,体系会遍历势能面上的各个点,计算的样本正就是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能,进而计算构型积分。 作用势的选择与动力学计算的关系极为密切,选择不同的作用势,体系的势能面会有不同的形状,动力学计算所得的分子运动与分子内部运动的轨迹也会不同,进而影响到抽样的结果与抽样结果的势能计算,在计算宏观体积与微观成分关系的时候主要采用刚球模型的二体势,计算系统能量,熵等关系时早期多采用Lennard-Jones、morse势等双体势模型,对于金属计算,主要采用morse势,但就是由于通过实验拟合的对势容易导致柯西关系,与实验不符,因此在后来的模拟中有人提出采用EAM等多体势模型,或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但就是对于二体势模型,多体势往往缺乏明确的表达式,参量很多,模拟收敛速度很慢,给应用带来很大困难,因此在一般应用中,通过第一性原理计算结果拟合势函数的L-J,morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想就是赋予分子体系初始运动状态之后,利用分子的自然运动在相空间中抽取样本进行统计计算,时间步长就就是抽样的间隔,因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞,体系数据溢出;太短的时间步长会降低模拟过程搜索相空间的能力,因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但就是通常情况下,体系各自由度中运动周期最短的就是各个化学键的振动,而这种运动对计算某些宏观性质并不产生影响,因此就产生了屏蔽分子内部振动或其她无关运动的约束动力学,约束动力学可以有效地增长分子动力学模拟时间步长,提高搜索相空间的能
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法） 1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。 由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。 进入生产相之后体系中的分子和分子中的原子开始根据初始速度运动，可以想象其间会发生吸引、排斥乃至碰撞，这时就根据牛顿力学和预先给定的粒子间相互作用势来对各个例子的运动轨迹进行计算，在这个过程中，体系总能量不变，但分子内部势能和动能不断相互转化，从而体系的温度也不断变化，在整个过程中，体系会遍历势能面上的各个点，计算的样本正是在这个过程中抽取的。 用抽样所得体系的各个状态计算当时体系的势能，进而计算构型积分。 作用势的选择与动力学计算的关系极为密切，选择不同的作用势，体系的势能面会有不同的形状，动力学计算所得的分子运动和分子内部运动的轨迹也会不同，进而影响到抽样的结果和抽样结果的势能计算，在计算宏观体积和微观成分关系的时候主要采用刚球模型的二体势，计算系统能量，熵等关系时早期多采用Lennard-Jones、morse势等双体势模型，对于金属计算，主要采用morse势，但是由于通过实验拟合的对势容易导致柯西关系，与实验不符，因此在后来的模拟中有人提出采用EAM等多体势模型，或者采用第一性原理计算结果通过一定的物理方法来拟合二体势函数。但是对于二体势模型，多体势往往缺乏明确的表达式，参量很多，模拟收敛速度很慢，给应用带来很大困难，因此在一般应用中，通过第一性原理计算结果拟合势函数的L-J，morse等势模型的应用仍非常广泛。 分子动力学计算的基本思想是赋予分子体系初始运动状态之后，利用分子的自然运动在相空间中抽取样本进行统计计算，时间步长就是抽样的间隔，因而时间步长的选取对动力学模拟非常重要。太长的时间步长会造成分子间的激烈碰撞，体系数据溢出；太短的时间步长会降低模拟过程搜索相空间的能力，因此一般选取的时间步长为体系各个自由度中最短运动周期的十分之一。但是通常情况下，体系各自由度中运动周期最短的是各个化学键的振动，而这种运动对计算某些宏观性质并不产生影响，因此就产生了屏蔽分子内部振动或其他无关运动的约束动力学，约束动力学可以有效地增长分子动力学模拟时间步长，提高搜索相空间的能
vasp做分子动力学的好处，由于vasp是近些年开发的比较成熟的软件，在做电子scf速度方面有较好的优势。 缺点：可选系综太少。 尽管如此，对于大多数有关分子动力学的任务还是可以胜任的。 主要使用的系综是NVT和NVE。 下面我将对主要参数进行介绍！ 一般做分子动力学的时候都需要较多原子，一般都超过100个。 当原子数多的时候，k点实际就需要较少了。有的时候用一个k点就行，不过这都需要严格的测试。通常超过200个原子的时候，用一个k点，即Gamma点就可以了。 INCAR： EDIFF 一般来说，用1E-4或者1E-5都可以，这个参数只是对第一个离子步的自洽影响大一些，对于长时间的分子动力学的模拟，精度小一点也无所谓，但不能太小。 IBRION=0 分子动力学模拟 IALGO=48 一般用48，对于原子数较多，这个优化方式较好。 NSW=1000 多少个时间步长。 POTIM=3 时间步长,单位fs,通常1到3. ISIF=2 计算外界的压力. NBLOCK= 1 多少个时间步长，写一次CONTCAR，CHG和CHGCAR，PCDAT. KBLOCK=50 NBLOCK*KBLOCK个步长写一次XDATCAR. ISMEAR=-1 费米迪拉克分布. SIGMA =0.05 单位:电子伏 NELMIN=8 一般用6到8,最小的电子scf数.太少的话,收敛的不好. LREAL=A APACO=10 径向分布函数距离,单位是埃. NPACO=200 径向分布函数插的点数. LCHARG=F 尽量不写电荷密度,否则CHG文件太大. TEBEG=300 初始温度. TEEND=300 终态温度。不设的话，等于TEBEG. SMASS -3 NVE ensemble；-1 用来做模拟退火；大于0 NVT 系综。 ///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// 1)收敛判据的选择 结构弛豫的判据一般有两种选择：能量和力。这两者是相关的，理想情况下，能量收敛到基态，力也应该是收敛到平衡态的。但是数值计算过程上的差异导致以二者为判据的收敛速度差异很大，力收敛速度绝大部分情况下都慢于能量收敛速度。这是因为力的计算是在能量的基础上进行的，能量对坐标的一阶导数得到力。计算量的增大和误差的传递导致力收敛慢。 到底是以能量为收敛判据，还是以力为收敛判据呢？关心能量的人，觉得以能量
分子动力学在材料科学中的应用 摘要：本文综述了几种常见条件下的分子动力学模拟方法以及分子动力学模拟的最新发展趋势.介绍用分子动力学模拟方法研究固休的休相结构,表面问题,界面问题以及薄膜形成过程等方面的研究成果。 关键词：分子动力学; 计算机模拟; 材料科学 1引言 分子动力学（Molecular Dyanmica,简称MD）用于计算以固体、液体、气体为模型的单个分子运动，它是探索各种现象本质和某些新规律的一种强有力的计算机模拟方法，具有沟通宏观特性与微观结构的作用，对于许多在理论分析和实验观察上难以理解的现象可以做出一定的解释［1］。MD方法不要求模型过分简化，可以基于分子（原子、离子）的排列和运动的模拟结果直接计算求和以实现宏观现象中的数值估算。可以直接模拟许多宏观现象，取得和实验相符合或可以比较的结果，还可以提供微观结构、运动以及它们和体系宏观性质之间关系的极其明确的图象［2］。MD以其不带近似、跟踪粒子轨迹、模拟结果准确［3］，而倍受研究者的关注，在物理、化学、材料、摩擦学等学科及纳米机械加工中得到广泛而成功的应用。本文主要评述MD方法在材料科学中的应用. 目前在材料微观结构的研究中，由于实验条件的限制，使得许多重要的微观结构的信息难以得到，如，对于由液态金属快速凝固的非晶转变过程，其微观结构的瞬时变化根本无法用实验仪器去测量。理论分析、实验测定及模拟计算已成为现代材料科学研究的3种主要方法[2]。20世纪90年代以来，由于计算机科学和技术的飞速发展，模拟计算的地位日渐突显。计算机模拟可以提供实验上尚无法获得或很难获得的信息。虽然计算机模拟不能完全取代实验，但可以用来指导实验的进行，从而促进理论和实践的发展，所以有必要对这一领域进行介绍。 2 分子动力学基本原理 分子动力学将连续介质看成由N个原子或分子组成的粒子系统，各粒子之间的作用力可以通过量子力学势能函数求导得出，忽略量子效应后，运用经典牛顿力学建立系统粒子运动数学模型，通过数值求解得到粒子在相空间的运动轨迹，然后由统计物理学原理得出该系统相应的宏观动态、静态特性。图1所示是MD
1、分子动力学简介： 分子动力学方法是一种计算机模拟的实验方法，是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。它是对理论计算和实验的有力补充。广泛应用于材料科学、生物物理和药物设计等。经典MD模拟，其系统规模在一般的计算机上也可达到数万个原子，模拟时间为纳秒量级。 分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和绝热近似（Born-Oppenheimer）下，分子动力学的这一种假设是可行的。所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。要进行分子动力学模拟就必须知道原子间的相互作用势。 在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。然而采用经验势必然丢失了局域电子结构之间存在的强相互作用的信息，即不能得到原子动力学过程中的电子性质。
2、分子模拟的三步法和大致分类 三步法： 第一步：建模。包括几何建模，物理建模，化学建模，力学建模。初始条件的设定，这里要从微观和宏观两个方面进行考虑。 第二步：过程。这里就是体现所谓分子动力学特点的地方。包括对运动方程的积分的有效算法。对实际的过程的模拟算法。关键是分清楚平衡和非平衡，静态和动态以及准静态情况。 第三步：分析。这里是做学问的关键。你需要从以上的计算的结果中提取年需要的特征，说明你的问题的实质和结果。因此关键是统计、平均、定义、计算。比如温度、体积、压力、应力等宏观量和微观过程量是怎么联系的。 大致分类： 2.1电子模拟（量化计算，DFT） 量子化学计算 一般处理几个到几十个原子 常见软件：GAUSSIAN，NWCHEM等 密度泛函（DFT） 可以算到上百个原子 常见软件：V ASP 2.2分子模拟（分子动力学，蒙特卡洛） 2.2.1分子级别的模拟 以分子的运动为主要模拟对象。采用经验性的分子间作用函数
《材料计算设计基础》 学号： 流水号： 姓名： 完成日期：
分子动力学模拟及其在材料中的研究进展 摘要：本文综述了分子动力学模拟技术的发展，介绍了分子动力学的分类、运动方程的求解、初始条件和边界条件的选取、平衡系综及其控制、感兴趣量的提取以及分子动力学模拟在材料中的研究进展。 关键词：分子动力学模拟平衡态系综金属材料感兴趣量径向分布函数 引言 科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。这种优点使分子动力学模拟在金属材料研究中显得非常有吸引力。 分子动力学MD (Molecular Dynamics)模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。MD模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。该模拟技术主要涉及粒子运动的动力学问题,与蒙特卡罗模拟方法(简称MC)相比,分子动力学是一种“确定性方法”, 它所计算的是时间平均,而MC进行的是系综平均。然而按照统计力学各态历经假设,时间平均等价于系综平均。因此,两种方法严格的比较计算能给出几乎相同的结果。 经典的分子动力学方法是Alder等于1957年提出并首先在“硬球”液体模型下应用，发现了由Kirkwood在1939年根据统计力学预言的“刚性球组成的集合系统会发生有液相到结晶相的转变”。后来人们称这种相变为Alder相变。Rahman
分子动力学模拟位错和界面的相互作用
学校代码10530 学号200910081121 分类号密级 硕士学位论文 分子动力学模拟位错和界面的相互作用 学位申请人 周银库 指导教师 陈尚达副教授 学院名称 材料与光电物理学院 学科专业 材料科学与工程 研究方向 金属薄膜的力学性能
二零一二年五月 Molecular dynamics simulations of interaction between dislocations and interfaces Candidate Yinku Zhou Supervisor Shangda Chen (Associate Professor) College Faculty of Materials, Optoelectronics and Physics Program Material Science and Engineering
Specialization Mechanical Properties of metal film Degree Engineering Master University Xiangtan University Date May, 2012
湘潭大学 学位论文原创性声明 本人郑重声明：所呈交的论文是本人在导师的指导下独立进行研究所取得的研究成果。除了文中特别加以标注引用的内容外，本论文不包含任何其他个人或集体已经发表或撰写的成果作品。对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。本人完全意识到本声明的法律后果由本人承担。 作者签名：日期：年月日 学位论文版权使用授权书 本学位论文作者完全了解学校有关保留、使用学位论文的规定，同意学校保留并向国家有关部门或机构送交论文的复印件和电子版，允许论文被查阅和借阅。本人授权湘潭大学可以将本学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存和汇编本学位论文。 涉密论文按学校规定处理。 作者签名：日期：年月日 导师签名：日期：年月日
分子动力学模拟基础知识 ? Molecular Dynamics Simulation o MD: Theoretical Background Newtonian Mechanics and Numerical Integration The Liouville Operator Formalism to Generating MD Integration Schemes o Case Study 1: An MD Code for the Lennard-Jones Fluid Introduction The Code, mdlj.c o Case Study 2: Static Properties of the Lennard-Jones Fluid (Case Study 4 in F&S) o Case Study 3: Dynamical Properties: The Self-Diffusion Coefficient ? Ensembles o Molecular Dynamics at Constant Temperature Velocity Scaling: Isokinetics and the Berendsen Thermostat Stochastic NVT Thermostats: Andersen, Langevin, and Dissipative Particle Dynamics The Nosé-Hoover Chain Molecular Dynamics at Constant Pressure: The Berendsen Barostat Molecular Dynamics Simulation We saw that the Metropolis Monte Carlo simulation technique generates a sequence of states with appropriate probabilities for computing ensemble averages (Eq. 1). Generating states probabilitistically is not the only way to explore phase space. The idea behind the Molecular Dynamics (MD) technique is that we can observe our dynamical system explore phase space by solving all particle equations of motion . We treat the particles as classical objects that, at least at this stage of the course, obey Newtonian mechanics. Not only does this in principle provide us with a properly weighted sequence of states over which we can compute ensemble averages, it additionally gives us time-resolved information, something that Metropolis Monte Carlo cannot provide. The ``ensemble averages'' computed in traditional MD simulations are in practice time averages : (99) The ergodic hypothesis partially requires that the measurement time, , i , in the system. The price we pay for this extra information is that we must at least access if not store particle velocities in addition to positions, and we must compute interparticle forces in addition to potential energy. We will introduce and explore MD in this section.
分子动力学(MD) 1 分子动力学（MD）基础 1.1 MD分类 1.2 MD简介 1.3 MD适用范围 2 分子动力学运动方程数值求解 2.1 基础知识 2.1.1 运动方程 2.1.2 空间描述 2.1.3 最小作用量原理 2.1.4 拉格朗日(Lagrange)方程 2.1.5 哈密顿(Hamilton)方程 2.2 粒子运动方程的数值解法 2.2.1 Verlet算法 2.2.2 欧拉(Euler)预测—矫正公式 2.2.3 Gear预测—矫正方法 3 分子动力学原胞与边界条件 3.1 分子动力学原胞 3.2 边界条件 3.2.1 自由表面边界 3.2.2 固定边界 3.2.3 柔性边界 3.2.4 周期性边界 4 势函数与分子力场 4.1 势函数 4.1.1 两体势 4.1.2 多体势 4.2 分子力场 4.2.1 分子力场函数的构成
4.2.2 常用力场函数和分类 5 分子动力学模拟的基本步骤 5.1 设定模拟所采用的模型 5.2 给定初始条件 5.3 趋于平衡计算 5.4 宏观物理量的计算 6 平衡态分子动力学模拟 6.1 系综 6.2 微正则系综的分子动力学模拟6.3 正则系综的分子动力学模拟
1 分子动力学（MD）基础 1.1MD分类 微正则系综（VNE） 正则系综（VNP） 平衡态MD 等温等压系综（NPT） 经典MD 等焓等压系综（NPH） 巨正则系综（ＶＴμ） 非平衡态ＭＤ 量子MD 1.2分子动力学(MD)简介 分子动力学是在原子、分子水平上求解多体问题的重要的计算机模拟方法。分子动力学方法为确定性模拟方法，广泛地用于研究经典的多粒子体系的研究中,是按该体系内部的内禀动力学规律来计算并确定位形的转变。 分子动力学方法是通过建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性, 从而得到系统的宏观性质。 在分子动力学中，粒子的运动行为是通过经典的Newton运动方程所描述。系统的所有粒子服从经典力学的运动规律，它的动力学方程就是从经典力学的运动方程——拉格朗日(lagrange)方程和哈密顿(Hamilton)方程导出。 1.3适用范围 原则上，分子动力学方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其它的微观粒子。 实际上，分子动力学模拟方法和随机模拟方法一样都面临着两个基本限制：
分子动力学模拟 The Standardization Office was revised on the afternoon of December 13, 2020
分子动力学模拟 分子动力学是一门结合物理，数学和化学的综合技术。分子动力学是一套分子模拟方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的热力学量和其他宏观性质。 这门技术的发展进程是： 1980年：恒压条件下的动力学方法（Andersenの方法、Parrinello-Rahman法）1983年：非平衡态动力学方法(Gillan and Dixon) 1984年：恒温条件下的动力学方法（能势‐フーバーの方法） 1985年：第一原理分子动力学法（→カー?パリネロ法） 1991年：巨正则系综的分子动力学方法(Cagin and Pettit). 最新的巨正则系综，即为组成系综的系统与一温度为T、化学势为μ的很大的热源、粒子源相接触，此时系统不仅同热源有能量交换，而且可以同粒子源有粒子的交换，最后达到平衡，这种系综称巨正则系综。 进行分子动力学模拟的第一步是确定起始构型，一个能量较低的起始构型是进行分子模拟的基础，一般分子的其实构型主要是来自实验数据或量子化学计算。在确定起始构型之后要赋予构成分子的各个原子速度，这一速度是根据玻尔兹曼分布随机生成，由于速度的分布符合玻尔兹曼统计，因此在这个阶段，体系的温度是恒定的。另外，在随机生成各个原子的运动速度之后须进行调整，使得体系总体在各个方向上的动量之和为零，即保证体系没有平动位移。由上一步确定的分子组建平衡相，在构建平衡相的时候会对构型、温度等参数加以监控。
【专业】计算物理【研究方向】分子动力学模拟【学术讲坛】 1、分子动力学简介：分子动力学方法是一种计算机模拟实验方法，是研究凝聚态系统的有力工具。该技术不仅可以得到原子的运动轨迹，还可以观察到原子运动过程中各种微观细节。它是对理论计算和实验的有力补充。广泛应用于材料科学、生物物理和药物设计等。经典MD模拟，其系统规模在一般的计算机上也可达到数万个原子，模拟时间为纳秒量级。2006年进行了三千二百亿个原子的模拟（IBM lueGene/L）。分子动力学总是假定原子的运动服从某种确定的描述，这种描叙可以牛顿方程、拉格朗日方程或哈密顿方程所确定的描述，也就是说原子的运动和确定的轨迹联系在一起。在忽略核子的量子效应和Born-Oppenheimer绝热近似下，分子动力学的这一种假设是可行的。所谓绝热近似也就是要求在分子动力学过程中的每一瞬间电子都处于原子结构的基态。要进行分子动力学模拟就必须知道原子间的相互作用势。在分子动力学模拟中，我们一般采用经验势来代替原子间的相互作用势，如Lennard-Jones势、Mores势、EAM原子嵌入势、F-S多体势。然而采用经验势必然丢失了局域电子结构之间存在的强相关作用信息，即
不能得到原子动力学过程中的电子性质。详细介绍请见附件。 2、分子模拟的三步法和大致分类三步法： 第一步：建模。包括几何建模，物理建模，化学建模，力学建模。初始条件的设定，这里要从微观和宏观两个方面进行考虑。 第二步：过程。这里就是体现所谓分子动力学特点的地方。包括对运动方程的积分的有效算法。对实际的过程的模拟算法。关键是分清楚平衡和非平衡，静态和动态以及准静态情况。 第三步：分析。这里是做学问的关键。你需要从以上的计算的结果中提取年需要的特征，说明你的问题的实质和结果。因此关键是统计、平均、定义、计算。比如温度、体积、压力、应力等宏观量和微观过程量是怎么联系的。有了这三步，你就可以做一个好的分子动力学专家了。推而广之，其实所谓的介观模拟，蒙特卡罗模拟、有限元模拟都是一个道理。大致分类：?量子化学计算?2.1电子模拟（量化计算，DFT）常见软件：VASP?可以算到上百个原子?密度泛函（DFT）?常见软件：GAUSSIAN，NWCHEM等?一般处理几个到几十个原子 2.2分子模拟（分子动力学，蒙特卡洛）
分子动力学模拟： 对于原子核和电子组成的多体体系，求解运动方程（哈密顿，牛顿，拉格朗日），用经典和量子化方法求解粒子的运动状态。 MC方法：系综（抽样）平均法分子动力学：时间平均 一优点：遇到的不利影响因素回避掉,从而达到实验研宄难以实现的控制条件。核心算法：粒子的运动状态就必须把运动方程离散化，离散化的方法有经典Verlet算法、蛙跳算法(Leap-frog)、速度Veriet算法、Gear预估-校正法等。 缺点：元胞体积和形状不变，不含有自由电子，对金属体系计算不理想。 注意：一般而言，MD模拟时间足够长，初始条件不会影响计算结果，但是会加大构型平衡的计算时间。 二步骤： 1.选取所要研究的系统并建立适当的模拟模型。 2.设定区域的边界条件，选取粒子间相互作用势模型；要注意观察PBC边界条 件的使用，以及计算格子和建模的晶格子之间的关系。体系是单胞沿不同方向重复叠合而组成。但模拟时只保留基本单元，由平移对称矩阵计算得到其他原子的空间坐标。最小近邻的截断半径。 3.设定系统所有粒子的初始位置和初始速度； 4.计算粒子间相互作用力和势能，以及各个粒子的位置和速度；最好与实际模 型相符，以减少达到平衡的时间。势场参数调整，最小近邻的截断半径。 对势：LJ势（惰性气体，过渡金属，柔韧材料），Born-lande势（离子晶体），Morse势，Johnson势（金属）发展到三体势，缺点是导致Cauchy关系，即不能描述晶体的弹性性质。
多体势：80年代以后，EAM势等（晶体对势+核嵌入电子云嵌入能），多用于金属。 5.待体系达到平衡后，构型积分获得体系的宏观性质。选取合适的系综，控制 温度和压力的变化。 控温方法：
分子动力学方法 一、引言 计算机模拟中的另一类确定性模拟方法，即统计物理中的所谓合于动力学方法（Molecular Dynamics Method）。这种方法是按该体系内部的内禀动力学规律来计算并确定位形的转变。它首先需要建立一组分子的运动方程，并通过直接对系统中的一个个分子运动方程进行数值求解，得到每个时刻各个分子的坐标与动量，即在相空间的运动轨迹，再利用统计计算方法得到多体系统的静态和动态特性，从而得到系统的宏观性质。在这样的处理过程中我们可以看出：MD方法中不存在任何随机因素。在MD方法处理过程中方程组的建立是通过对物理体系的微观数学描述给出的。在这个微观的物理体系中，每个分子都各自服从经典的牛顿力学。每个分子运动的内禀动力学是用理论力学上的哈密顿量或者拉格朗日量来描述，也可以直接用牛顿运动方程来描述。确定性方法是实现Boltzman的统计力学途径。这种方法可以处理与时间有关的过程，因而可以处理非平衡态问题。但是使用该方法的程序较复杂，讨算量大，占内存也多、本节将介绍分子动力学方法及其应用。 原则上，MD方法所适用的微观物理体系并无什么限制。这个方法适用的体系既可以是少体系统，也可以是多体系统；既可以是点粒子体系，也可以是具有内部结构的体系；处理的微观客体既可以是分子，也可以是其他的微观粒子。 实际上，MD模拟方法和随机模拟方法一样都面临着两个基本限制：一个是有限观测时间的限制；另一个是有限系统大小的限制。通常人们感兴趣的是体系在热力学极限下（即粒子数日趋于无穷时）的性质。但是计算机模拟允许的体系大小要比热力学极限小得多，因此可能会出现有限尺寸效应。为了减小有限尺寸效应，人们往往引入周期性、全反射、漫反射等边界条件。当然边界条件的引入显然会影响体系的某些性质。 对于MD方法，向然的系综是微正则系综，这时能量是运动常量。然而，当我们想要研究温度和（或）压力是运动常量的系统时，系统不再是封闭的。例如当温度为常量的系统可以认为系统是放置在一个热俗中。当然，在MD方法中我们只是在想像中将系统放入热浴中。实际上，在模拟计算中具体所采取的做法是对一些自由度加以约束。例如在恒温体系的情况下，体系的平均动能是一个不变量。这时我们可以设计一个算法，使平均动能被约束在一个给定值上。由于这个约束，我们并不是在真正处理一个正则系综，而实际上仅仅是复制了这个系综的位形部分。只要这一约束不破坏从一个状态到另一个状态的马尔科夫特性，这种做法就是正确的。不过其动力学性质可能会受到这一约束的影响。 自五十年代中期开始，MD方法得到了广泛的应用。它与蒙特卡洛方法一起已经成为计算机模拟的重要方法。应用MD方法取得了许多重要成果，例如气体或液体的状态方程、相变问题、吸附问题等，以及非平衡过程的研究。其应用已从化学反应、生物学的蛋白质到重离子碰撞等广泛的学科研究领域。 二、分子运动方程的数值求解 采用MD方法时，必须对一组分于运动微分方程做数值求解。从计算数学的角度来看，这个求解是一个初值问题。实际上计算数学为了求解这种问题己经发展了许多的算法，但并不是所有的这些算法都可以用来解决物理问题。下面我们先以一个一维谐振子为例，来看一下如何用计算机数值计算方法求解初值问题。一维谐振子的经典哈密顿量为 （2.1） 这里的哈密顿量（即能量）为守恒量。假定初始条件为x(p)、p(0),则它的哈密顿方程是对时间的一阶微分方程 （2.2） 现在我们要用数值积分方法计算在相空间中的运动轨迹(X(t)、p(t)) 。我们采用有限差分法，将微分方程变为有限差分方程，以便在计算机上做数值求解，并得到空间坐标和动量随时间的演化关系。首先，