计算材料科学的模型和算法和多尺度关联

• 例: COMPASS-98力场(condensed-phase optimized molecular potentials for atomistic simulation studies)的表达式如下
键伸缩：Eb = ∑ [k2 (b − b0 ) 2 + k3 (b − b0 )3 + k4 (b − b0 ) 4 ]
r
θ E
r
分子势函数曲面——势能面 示意图
θ
r
过渡状态
E
产物 反应物
• 1、分子力场 、
• 分子片段力场的函数表达式中包含自变量和力场参数 • 其中自变量为分子的结构参数，独立参数为键长、键 角和二面角，如图
b b θ θ θ’ φ
b’
有的还使用一个非独立参数： 面外弯曲角
χ
• 而 力场参数 一般通过与实验数据 和 从头算数据进 行最小二乘法拟合来确定
b ,θ b ,ϕ
' + ∑ (b ' − b0 )[k1 cos ϕ + k 2 cos 2ϕ + k 3 cos 3ϕ ] b ',ϕ
+ ∑ (θ − θ 0 )[k1 cos ϕ + k 2 cos 2ϕ + k 3 cos 3ϕ ]
θ ,ϕ
键面外弯曲：Eχ = ∑ k χ χ 2
χ
Biblioteka Baidu
显然是一个 + ∑ kϕθθ ' cosϕ (θ − θ 0 )(θ ' − θ 0' ) 非谐性力场
b
• 每个k是一独立的力场参数,下标“0”代表参考(~平衡)结构参数：
键弯曲：Eθ = ∑ [k2 (θ − θ 0 ) 2 + k3 (θ − θ 0 )3 + k4 (θ − θ 0 ) 4 ]
θ
二面角：Eϕ = ∑{k1[1 − cos(ϕ − ϕ 0,1 )] + k2 [1 − cos(2ϕ − ϕ 0, 2 ) + k3[1 − cos(3ϕ − ϕ 0,3 )]}
范德华势：EvdW
0 6 r 0 9 rij ij = ∑ ε ij 2 − 3 r rij i, j ij
即Lennard-Jones 9-6函数
其它力场范德华势较多采用L-J 12-6函数：
EvdW r 0 12 r 0 6 ij ij = ∑ ( D0 ) ij − 2 r rij i, j ij
(i=1,2,3,…3N-6)
• 1) 用各种方法构造出一个分子的任意结构，得到初始的结构参 数{xi0, i=1,2,3,…3N-6} • 2) 进行坐标变换，即根据需要将原子的直角坐标转变成内坐标 (键长、键角、二面角)或反之
3) 建立分子体系的势能表达式E(x) 4) 计算E(x)随各坐标的一阶、二阶导数 5) 计算接近数学条件的坐标增量 6) 得到新的结构参数{xi1, i=1,2,3,…3N-6} 重复 4)、5)、6)，直至最后两次得到的体系势能之差 或总体势能梯度的均方根值达到预定精度范围为止 • MM计算中，有时还要考虑体系所处的外压条件，此 时，可通过压力因子的定义来调节原子坐标 • • • • •
• 非键合势函数中，静电相互作用表示分子中各原子静 电荷的库仑相互作用对势能的贡献 • 不同的力场，静电相互作用表达式基本相同 • 范德华势也大都采用Lennard-Jones函数，但函数中的 指数有所不同。如COMPASS-98的非键合势函数为： 静电相互作用： Eelec = ∑
i> j
qi q j rij
• 势函数形式很多，目前已被广泛使用的力场有如CFF、 MM2、MM3、MM4、MMFF、AMBER、CHARMM、 DREIDING、UFF和COMPASS等 • 形式虽多, 但一般总表达为分子内与分子间势能之和: V总=V键合+V非键合 • 分子内势能(键合)包括键伸缩、键角弯曲和二面角扭转势能 • 分子间势能(非键合)包括范德华势和静电势, 有的还包括H键： V键合=V键伸缩+V键角弯曲+V二面角扭转 V非键合=V范德华+V静电+V氢键 • 键合势函数中，一些力场还包含交叉项，使精度更高 • 交叉项的含义：如键长变化时，键角弯曲势能随键长 的不同而不同，等
基本思想
• 事先构造出简单体系(如链段、官能团等各种不同结构的小 片段)的势能函数, 简称 势函数 或 力场(force field) • 将势函数建成数据库，在形成较大分子的势函数时， 从数据库中检索到结构相同的片段，组合成大体系的 势函数 • 利用分子势能随原子位置的变化有极小值的性质，确 定大分子的结构即为分子力学(MM) • 利用势函数，建立并求解与温度和时间有关的牛顿运 动方程，得到一定条件下体系的结构随时间的演化关 系即为分子动力学(MD) • 理论方法的核心 核心是构造势函数 核心 • 势函数：势能与原子位置的关系。且往往是不知道的 • 需要通过其他方法，如量子化学方法及实验数据获得
• 2、分子力学方法 、
• MM是确定分子结构的方法 • 利用分子势能随结构的变化而变化的性质，确定分子 势能极小时的平衡结构(stationary point) • 物理模型: 视原子为质点，视化学键为弹簧，而弹力常 数完全由数据库中的分子力场来确定 • 因此是直接用势函数研究问题，不考虑原子的动能 • 不考虑动能所对应的结构，相当于体系处于 T=0 K 时 的结果
• 当然，在建立分子的势能函数时，还有一些更细致的 问题要考虑，如： – 势能展开项的截断 – 周期结构的处理 – 多组分混合物体系 – 含有离子的体系 或 – 金属中的离子 等 • 不同的方法或程序中，分别都有更详细的讨论。参考：
[德] D. 罗伯. 计算材料学. 北京：化学工业出版社, 2002, 9 俞庆森,朱龙观. 分子设计导论.北京：高等教育出版社,2000 杨小震. 分子模拟与高分子材料. 北京：科学出版社，2002 熊家炯主编. 材料设计. 天津：天津大学出版社，2000 Sun H，Ren P, Fried J R. The COMPASS Force Field: Parameterization and Validation for polyphosphazenes. Computational and Theoretical Polymer Science, 1998, 8(1/2): 229 • Sun H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications - Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem., 1998, 102: 7338 • • • • •
• 3、分子动力学方法 、
• 分子力场是分子的静态势函数 • 而实际过程通常是在一定温度和一定压力下发生的 • 为了更切实际地了解体系运动和演化的过程，必须考 虑体系中原子的运动，并与温度T和时间t建立联系 • 我们知道，温度是原子分子热运动剧烈程度的量度 • 根据统计热力学，对于n个原子的体系，体系的温度T 与各原子的运动速率vi的关系为：
2
(i=1,2,3,…3N-6)
进行结构优化 具体步骤： • 1) 用各种方法构造出一个分子的任意结构，得到初始 的结构参数{xi0, i=1,2,3,…3N-6} • 2) 进行坐标变换，即根据需要将原子的直角坐标转变 成内坐标{ri、θi、φi}或反之
∂2 E ( x) ∂E ( x ) = 0, >0 2 ∂xi ∂xi
• 计算过程一般为：在一定的统计系综下 • 1) 由原子位置和连接方式，从数据库调用力场参数并 形成体系势函数 • 2) 由给定温度计算体系动能以及总能量 • 3) 计算各原子的势能梯度, 得到原子在力场中所受的力 即dp/dt=mdv/dt=ma=F
dp ∂H ( x, p ) =− , dt ∂x
H ( x, p ) = K ( p ) + E ( x )
• 与量子力学不同，经典力学对Hamilton量不进行算符 化处理，也不建立和求解本征方程，而是建立并求解 经典运动方程： dp ∂H ( x, p )
dt = − ∂x dx ∂H ( x, p ) = dt ∂p
1 n k BT = mi vi2 ∑ 3n i =1
• 又因体系中各原子的速率为vi时，动量pi=mivi，对应总 n 动能K(p)为： 1 2 2 2
K ( p) = ∑
i =1
2 mi
( p i , x + p i , y + pi , z )
• 势能由力场确定为E(x)，因此体系的Hamilton量H为：
• 由力场首先构造并得到分子的势函数E(x) • 利用在E(x)的极小点处，E(x)随各原子独立的空间坐标 {xi, i=1,2,3,…3N-6} (与内坐标{ri、θi、φi}等价) 的一阶 微分=0以及全部二阶微分>0的数学条件：
∂E ( x) =0 ∂xi ∂ E ( x) >0 2 ∂xi
ϕ
' 交叉项：E x = ∑ kbb (b − b0 )(b' − b0 ) + ∑ kθθ (θ − θ 0 )(θ ' − θ 0' ) b ,b '
θ ,θ
+ ∑ k bθ (b − b0 )(θ − θ 0 ) + ∑ (b − b0 )[k1 cosϕ + k 2 cos 2ϕ + k 3 cos 3ϕ ]
dx ∂H ( x, p ) = dt ∂p
计算材料科学中的模型、 计算材料科学中的模型、算法和 多尺度关联 系列讲座 与 讨论
西北工业大学科研处 西北工业大学理学院 西北工业大学材料学院 西北工业大学高性能计算中心 联合举办
III、纳米、亚纳米尺度
一、第一原理方法及其在材料科学中的应用 二、密度泛函方法及其应用
三、从原子分子到纳米尺度
——分子力学、分子动力学方法及其应用
ϕ ,θ ,θ '
• 力场参数 最小二乘法确定 力场参数{k}最小二乘法确定
基本思想 • 如：R-COOH基团 • 1) 由ab initio (构型优化方法)计算出 平衡结构，得到 参考结构参数{bi0, θi0, φi0}
C
O O H
• 2) 用伪随机数方法将{bi0, θi0, φi0}人为改变成若干(n)组 非平衡结构参数 {bi, θi, φi}l (l=1,2, …, n) • 3) 用各{bi, θi, φi}l 结构参数分别进行ab initio计算，得 到{bi, θi, φi}l 结构参数下对应的能量El (l=1,2, …, n) • 4) 将El (l=1,2, …, n)和{bi, θi, φi}l代入势能表达式 • 5) 用最小二乘法拟合，确定力场参数{k}
从原子分子到纳米尺度
——分子力学、分子动力学方法及其应用 Molecular Mechanics / Molecular Dynamics
• 一、 MM、MD理论基础 • 二、 MM、MD计算程序 • 三、 MM、MD方法的应用
一、 MM、MD基础理论 、 基础理论
• 原则上, 第一原理方法在理论上已经能解决所有问题 • 但计算量太大，计算机资源有限，原子数目较多时，如高分子、 蛋白质、原子簇以及研究表面问题、功能材料或材料的力学性 能等，实际上难以完成计算 • 为此，发展了分子力学(Molecular Mechanics, MM)与分子动力 学(Molecular Dynamics, MD)方法 • 它们的应用，又称分子模拟(molecular simulation, molecular modeling) 或 分子设计(molecular design) • MM与MD是经典力学方法，针对的最小结构单元不再是电子 电子而 电子 是原子 • 因原子的质量比电子大很多，量子效应不明显，可近似用经典 力学方法处理 • 20 世纪 30 年代, Andrews 最早提出分子力学(MM)的基本思想； 40 年代以后得到发展, 并用于有机小分子研究。90年代以来得 到迅猛发展和广泛应用