分子动力学-ZHJ

势函数
原子间相互作用势根据来源可分为经典势和第一 性原理势 经典势又可根据使用范围分为原子间相互作用势 和分子间相互作用势 根据势函数的形式由可分为对势和多体势 势函数的研究和开发是分子动力学发展的最重要 的任务之一
v为速度矢量， mi为粒子质量
注：这些方程的求解一般来说需要通过数值方法进行；这些数值解产生一 系列的位置与速度对{xn,vn}；n 表示一系列的离散的时间，t＝nΔt, Δt表示时 间间隔(时间步长)；求解此方程组，须给出体系每个粒子的初始坐标和速度
分子动力学方法工作框图
分子动力学适用范围
内涵：分子动力学方法只考虑多体系统中原子核 的运动，而电子的运动不予考虑，量子效应忽略。 适用：经典近似在很宽的材料体系都较精确 不适用：涉及电荷重新分布的化学反应、键的形 成与断裂、解离、极化以及金属离子的化学键都 不适用，也不适用于低温 这几个不适用的领域可用量子力学方法替代，也 有结合了经典分子动力学和量子力学的新方法， 叫做第一性原理分子动力学，主要有CPMD和 BOMD
t
2
) [ri (t t ) ri (t )]/ t
Verlet 算法优缺点
Verlet 算法简单，存储要求适度 缺点 位置r(t+δt)要通过小项与非常大的两项2r(t) 与r(t-δt)的差得到，这容易造成精度损失 方程中没有显式速度项，在下一步的位置 没得到之前，难以得到速度项 不是自启动算法，新位置必须由t 时刻与前 一时刻t－δt 的位置得到
该过程中，动能、势能相互转化，当动能、势能、总能量只在平均 值附近涨落时，体系就达到平衡了（对微正则系综）
时间步长
时间步长Δt的选取非常重要
不合适的时间步长可能会导致模拟的失败 或结果错误，或者造成模拟的效率太低
时间步长的选取参考原子或分子特征运动 频率
势函数
势函数是描述原子（分子）间相互作用的 函数 原子间相互作用控制着原子间的相互作用 行为，从根本上决定材料的所有性质，这 种作用具体由势函数来描述 在分子动力学模拟中，势函数的选取对模 拟的结果起着决定性的作用 在选取势函数时，不能任意选取一个函数 就进行模拟，要认真考察势函数是否合适
边界条件
三维周期性边界条件 模拟大块固体或液体 时，必需选取三维周 期性边界条件 在计算粒子受力时， 采样区的边长应当至 少大于两倍的rc ，使 粒子i 不能同时与j 粒 子和它的镜像粒子相 作用
周期性边界条件
边界条件
二维周期性边界条件 对于薄膜情况，可使 用二维周期边界条件， 认为薄膜在x-y平面 内无限扩展(存在周 期性边界条件)，而 在z 方向受到限制(不 赋予周期边界条件)
ri (t t ) ri (t ) t vi (t t / 2)
Leap-frog 算法优缺点
Leap-frog 算法相比Verlet 算法有两个优点： 它包括显速度项，并且计算量稍小 它也有明显的缺陷：位置与速度不是同步 的，这意味着在位置一定时，不可能同时 计算动能对总能量的贡献
分子动力学模拟
identation 2010年8月
-Ref. 张跃等 <计算材料学基础> [北京航空航 天大学出版社]
第一节 概论
概述 分子动力学思想 分子动力学方法工作框图 分子动力学适用范围
概述
什么是分子动力学 一种确定性方法 按照该体系内部的内禀动力学规律来确定 位形的转变，跟踪系统中每个粒子的个体 运动 根据统计物理规律，给出微观量（分子的 坐标、速度）与宏观可观测量（温度、压 力、比热容、弹性模量等）的关系 研究材料性能
Leap-frog 算法
vi (t t / 2) v i (t t / 2) Fi (t ) t mi
Fi (t ) v i (t t / 2) v i (t t / 2) t mi
ri (t t ) ri (t ) vi (t t / 2) t
Verlet 算法
在分子动力学中，积分运动方程运用最广泛的方法是 Verlet 算法，这种算法运用t 时刻的位置和速度及t-δ t时 刻的位置，计算出t+δ t时刻的位置r(t+δ t)
ri (t t ) ri (t ) tv i (t ) 1 t 2ai (t )+ 2 ri (t t ) ri (t ) tv i (t ) 1 t 2ai (t ) 2
薄膜周期性边界条件示意图
边界条件
非周期性边界条件 有时，仅对系统的一部分感兴趣(比如表面)，这 时可把系统分为两个部分：表面部分应用自由边 界条件，而另一部分可以应用周期性边界条件 有时，又需要用到固定边界条件，比如在有些单 向加载的模拟中 有时还要采用以上介绍的几种边界条件的结合应 用，即混合边界条件 在具体的应用中要根据模拟的对象和目的来选定 合适的边界条件
初始体系的设置
体系对初始条件和其他计算条件具有敏感性 初始化要求给每个粒子指定初始坐标和速度 实验结构已知：初始坐标和速度可从晶体结 构中得到，须进一步最小化来弛豫应力
指定的初始矢量未必对应于所使用的势函数的最小值
实验结构未知：需要根据已知结构来搭建一 个结构，这时需要进行最小化过程
初始体系的设置
力的计算方法
截断半径法
预先选一个截断半径rc ，只计算 以截断半径为球体内的粒子间的 作用力；而与粒子之间的距离超 过截断半径时，不考虑它们的作 用 当系统中粒子数量很大时，还可 以利用邻域列表法判断粒子的分 布情况，进一步节省机时
计算粒子1 受力时，只计算截断半 径以内的2、3、4 三个粒子对粒子1的 作用力之和，其他粒子的作用忽略不 计 计算过程中，每隔一定步数需要更 新此表
算法的选取
分子动力学常见算法：有Verlet算法、 Leap-frog 算法、Velocity Verlet 算法、 Gear 算法、Tucterman 和Berne 多时间步 长算法
选取算法时，其评判标准中最重要的一点 是能量守恒
第三节 运动方程的数值求解
Verlet 算法 Verlet 算法优缺点 Leap-frog 算法 Leap-frog 算法优缺点
分子动力学思想
经典力学定律
Fi mi ai (t )
Fi (t)为粒子所受的力，mi 为粒子的质 量，ai (t)为原子i 的加速
Fi mi ai (t )
原子i 所受的力Fi(t)可以直接用势能函数对坐标 ri 的一阶导数表示，从而
v i U mi Fi t ri r (t ) v (t )
初值问题
为了进行分子动力学模拟，建立系统的初 始构形是必要的 初始构形可以通过实验数据，或理论模型， 或两者的结合来获得 给每个原子赋初速度也是必要的，它可以 从模拟的温度下的Maxwell-Boltzmann 分布 来任意选取
2 mi 1 1 mi vix p (vix ) ( ) 2 exp[ ] 2 k BT 2 k BT
势函数的选取原则
势函数的选取原则
通过文献调研，查明它的出处、应用 范围 自己验证势函数的好坏，通过模拟材 料一些已知的性质来验证，然后再进 行使用
力的计算方法
背景：分子动力学模拟中所花费的计算时 间，90%以上是用来计算作用在原子上的力， 所用时间大致正比于原子数目的平方 减少计算量的方法 短程力：截断半径法，只计算力程以内的 作用力，大大减少了计算量 长程力：库仑力，需找到近似处理办法减 少计算量，常用Ewald 求和法
概述
分子动力学的原理&步骤
原理：多体问题的严格求解，需要建立并求解体 系的薛定谔方程，根据波恩-奥本海默近似，原子 核的运动可以用经典动力学方法处理 薛定谔方程 牛顿运动方程（简化计算） 步骤
• 建立一个由N 个粒子（分子）组成的模型体系 • 解N 个粒子（分子）组成的模型体系的牛顿运动方程直至 平衡 • 平衡后，进行材料性能的计算，对模拟结果进行分析
初值问题
Maxwell-Boltzmann 方程给出了质量为mi 的原子i 在温度T 下沿x 方向速度为vix 的概 率 Maxwell-Boltzmann 分布是一种Gaussian 分布，它可以用随机数发生器得到 大多数随机数发生器产生的随机数均匀分 布在0 ～ 1之间，通过变换得到Gaussian 分布
初始速度矢量根据伪随机数进行设置，使体系的 总动能与目标温度对应
利用能量均分定理，考虑体系总的自由度数NF 给某个速度分量设置麦克斯韦分布
1 3N Ek mi vi2 N F k BT / 2 2 i 1
弛豫：体系的初始坐标和初始速度设置以后，在 进行模拟体系的性质以前，首先必须使体系进行 趋于平衡的过程
势函数
表示原子（分子）间相互作用的函数，也 称力场 原子间相互作用控制着原子间相互作用行 为，从根本上决定材料的所有性质，具体 由势函数来描述 在原子间相互作用势的现实模型是作用在 粒子上的力将随着此粒子的位置或和它有 相互作用的粒子的位置的变化而改变 采用的势函数的准确与否，将直接影响到 模拟结果的精确度
第四节 边界条件和初值
边界条件 初值问题
边界条件
模拟能力限制，不能模拟大量分子，只能模拟有限空 间中的有限个分子：有限空间边界 边界条件的选取 对于分子动力学模拟来说，合适的边界条件的选 取需要考虑两个方面的问题。第一，为减小计算量， 模拟的单元应尽可能小，同时模拟原胞还应足够大， 以排除任何可能的动力学扰动而造成对结果的影响， 此外模拟的原胞必须足够大还可以满足统计学处理的 可靠性要求；第二，还要从物理角度考虑体积变化、 应变相容性及环境的应力平衡等实际耦合问题。
第二节 主要技术
分子动力学运行流程图 初始体系的设置 时间步长 势函数 势函数的选取原则 力的计算方法 算法的选取
分子动力学运行流程图
进行分子动力学运算的几 个必备步骤：
首先建立计算模型 设定计算模型的初始坐标和 初始速度 选定合适的时间步长 选取合适的原子间相互作用 势函数，便于进行力的计算 选择合适的算法、边界条件 和外界条件 计算 对计算数据进行统计处理
Fi (t ) ri (t t ) 2ri (t ) ri (t t ) t mi
2
差分方程误差为 (δt)4数量级，速度 表达式未出现在 Verlet中，但可用 位置求得，其误差 为(δt)3数量级
v i (t ) [ri (t t ) ri (t t )]/2 t Or v i (t
邻域列表法
力的计算方法
Байду номын сангаасwald 加和法
Ewald 加和法是Ewald 于1921 年在研究离 子晶体的能量时发展的 在这个方法中，一个粒子与模拟的盒子中 所有其他原子以及周期性元胞中的镜像粒 子作用
力的计算方法
对一个位于晶格点n[=(nxL,nyL,nzL), nx, ny, nz 是整 数]的电荷所受势能为
初值问题
初始速度经常被校正以满足总动量为零 分别计算沿三方向的动量总和，然后用每 一方向的总动量除以总质量，得到一速度 值 用每个原子的速度减去此速度值，即保持 总动量为零
第五节 势函数
势函数 对势 适用于金属和合金的多体势 共价晶体的作用势 有机分子中的作用势(力场) 分子间作用势
N N qi q j 1 U 2 n i 1 j 1 4 0 rij n N N qi q j 1 U ' 2 n 0 i 1 j 1 4 0 rij n
问题：这个式子的求和收敛得非常慢，事实上它是条件 收敛的，Ewald 在求和时采用了一个技巧，即将上式转 变为两个收敛很快的级数，每个点电荷周围用两个大小 相等、符号相反的电荷分布包围（上图）