《物理问题的计算机模拟方法(1)—分子动力学》由会员分享，可在线阅读，更多相关《物理问题的计算机模拟方法(1)—分子动力学（22页珍藏版）》请在金锄头文库上搜索。
1、硕士研究生课程硕士研究生课程物理问题的计算机模拟方法物理问题的计算机模拟方法讲义讲义适用专业适用专业: : 凝聚态物理、材料物理与化学、理论物理、光学工程学学 时：时：3040 学时参考教材：参考教材：1德D.W.Heermann 著，秦克诚译，理论物理中的计算机模拟方法，北京大学出版社，1996。2荷 Frenkel 第三步：计算粒子 i 的新位置ir),(),(ziyixiiiiizyxzyx r第四步：计算粒子在新旧两个位置系统的能量之差),.,.,(),.,.,(2121NiNiUUUrrrrrrrr第五步：由的大小判断粒子 i 是否从 ri运动到:Uir若，则 ri；0Uir若，则再
2、选一个随机数 R(0 rc 的相互作用忽略不计，代价是忽略了背景。5积分格式 jiijij jiiN iirruU)()()()()(FrF于是，运动方程可写为(2-1) jiiji iirdtdmdtd)()( , Fvvr（1）递推公式在第一章中求解运动方程时，我们是直接求解的关于粒子位置坐标的二阶微分方程，得到的递推公式需要知道最初的两点位置才能启动计算，但在实际计算中，我们常常是给出最初的位置和速度，于是，我们可通过选取一定的差分格式，有运动方程(16)得到关于位置和速度的递推公式。或 已知 )1()()1()(n in in in i vvrr )()()()( htvtvhtrtr
3、 00 )0()0(vvrr递推公式的具体形式取决于差分各式的选取。（2）时间步长选择时间步长的原则：在保证计算精度的前提下，尽量节省计算时间。实例：Ar 原子系统，LJ 势时间步长取为 (=100ps=10fs)sh14210)(10约化时间步长（3）约束条件保证能量、动量或角动量守恒。（4）减小计算误差的技巧数值计算不可避免有误差，与误差有关的因素主要有：差分格式时间步长切断半径最小影像约定等等6. 计算热力学量dttVttAttAttNt)(),(, )()(lim00vr（1） 若系统的 NVE 恒定，则（微正则系综平均microcannonical ensemble average）
4、NVEAA（2） 若系统的 NVT 恒定，则（正则系综平均cannonical ensemble average）NVTAA在实际计算中，往往还需要要求系统的总动量 P=0 或恒定，因整个系统不受外界的作用。（3） 温度与能量均分定理)3(21ckNNkTEN总粒子数，Nc约束条件个数, 一般情况下，N NcP=0 Nc = 3（4） 位势切断误差与修正g(r) 对关联函数（pair correlation function） 粒子数密度VN当原点位置上一个粒子时，在 r 附近 dr 内发现有一个粒子的几rr dg)(率。对气体或液体， （各向同性）)()(rr平均每个粒子的能量cccrrrV
5、VdrrrgruNUdrrrgrudrrrgrudrrrgruddrrrgrudrgruNU02202022)()(2)()(2)()(2)()(2 )()(21)()(21r切断误差修正为：crcdrrrgruU2)()(26计算机模拟的组织（1） 初始化给定粒子的初始位置、初始速度等；（2） 趋衡从初始转台开始，按运动方程要求的规律，从非平衡态趋向平衡态；（3） 投产计算物理量的统计平均值。 2.2 微正则系综分子动力学方法微正则系综分子动力学方法微正则系综：NVE 恒定，且 P = 0 （分子动力学方法的自然选择）1Verlet 算法运动方程的显示中心差分格式mttttttiiii/)(
6、)()(2)(2Frrrttttttiii2/)()()(rrv由可得到如下递推公式（2-2） 2/ )(/211211hrrvmhFrrrn in in in in in in i此即 Verlet 算法算法（A2）, 其特点是：（1） 要启动此算法，必须已知两点，所以又称为二步法；10, iirr（2）和不能同步算出，第 n+1 步算出的第 n 步的速度。n irn iv2Verlet 算法的速度形式22 22)(21)()( )()( 21)()()(tmttttdttdtdttdttti iiii iiFvrrrrr2)()()()()()()()()()( tmtttttt tmtt
7、ttttmtttt FFvvFvvFvv（向后差分）（向前差分）由此可得如下递推公式：（2-3） 2/ )(2/1121mFFhvvmFhhvrrn in in in in in in in i此即 Verlet 算法的速度形式算法的速度形式（A3）, 其特点是：（1） 启动此算法，必须已知两点；00, iivr（2）和同步计算。 （的计算需要知道）1n ir1n iv1n iF1n ir此算法比 A2 算法更稳定。3趋恒阶段的能量调整由于系统初始状态的能量离我们要模拟系统的能量（或温度）有一定差异，这就需要在系统趋于平衡的阶段进行调整。最常用的方法之一就是进行速度标度：11n in ivv标
8、度因子：（2-4）2/12/)33(irefmvkTNTref 为设定的系统的温度值。能量设定值： refrefkTNE)33(21能量参考值： 2 21imvE2EEref2)(21irefvmE要注意的问题：（1） 一般不需要每一步都进行标度，可每隔若干步（比如 50 步）调整一次；（2） 当能量达到所设定的值后，停止标度，在以下的模拟时间内能量保持不变。4实例氩原子系统， N=256LJ 势: 612 4)(ijijijrrru作用力： 814221)(48)()(ijijji iij ijxrrxxxrurF 取时间单位：sm122/12 10348 取长度单位： = 0.3405nm
9、取质量单位：m=6.6338210-26 kg(氩原子质量) 8142814221482148)(ijijjiijijji ijxrrxxmmrrxxrF 使用上述单位，可得到约化单位下的作用力表达式： 8*14*11)()(ijijjiijxrrxxrF标度因子：2/12*2/12*2/12* 222/1216/)1( 16/)1( 48/)13( /)33(iirefirefirefvTNvkTNvmmkTNmvkTN约化温度： kTTkTrefref /*Kk8.119取时间步长 h = 210-14s = 20fs, 相当于约化时间 0.064约化温度：T* = 2.53 303K,
10、T* = 0.722 86.5K约化数密度：* = N/L3, N = 256* = 0.636 L = 7.83 , * = 0.83134 L = 6.75当 N = 64 时，* = 0.83134 L = 4.25, 所以，取 rc = 2.5 是错的，因为要求 rc L/2 = 2.13。 （习题 3.5）计算源程序见 p129, 程序 PL1 2.3 正则系综分子动力学方法正则系综分子动力学方法正则系综：NVT 恒定，且 P = 0 如何保持温度 T 恒定？1速度标度方法11n in ivv2/12/)33(irefmvkTN算法算法 A4:(1) 规定初始位置和初始速度；0 ir
11、0 iv(2) 计算； 2/ )(2/1121mFFhvvmFhhvrrn in in in in in in in i(3) 计算；2/12/)33(irefmvkTN(4) 对速度进行标度：，返回（2） 。11n in ivv说明：（1） 在前面的微正则系综的模拟中也对速度进行了标度，但其目的是使总能量达到所需要得值，而且不必每步进行速度标度，可隔若干步标度一次，一旦总能量达到所设定的值，即可停止标度；而在正则系综的模拟中，每步都必须进行标度，以始终保持温度恒定。（2） 通过一个广义位势引进能量涨落，连同细致耦合的一种特殊选择，会导致速度标度机制。 （见书中的证明，公式 3.59 有误）（
12、3） 由于控制温度的反馈环内的时间延迟，温度将在一定程度上涨落。2随机方法（Anderson 热浴）体系与一强加了温度的热浴相耦合，与热浴的耦合由偶尔作用于随机选择的粒子上的脉冲力表示。在开始模拟前首先要选择系统与热浴的耦合强度，这一耦合强度由随机碰撞的频率 决定。模拟步骤：（1）规定初始位置和初始速度；0 ir0 iv（2）计算； 2/ )(2/1121mFFhvvmFhhvrrn in in in in in in in i（3）粒子 i 在时间间隔 h 内发生碰撞的几率为 h，产生随机数，如果 h, 则粒子 i 经历一次随即碰撞，否则，返回（2） ；（4）粒子 i 从温度为 T 的麦克斯
13、韦玻尔兹曼分布中获得一个新的速度，返回（2） 。 2.4 等温等压系综分子动力学方法等温等压系综分子动力学方法1等压等焓系综N P H 恒定，P = 0 (P:压强，P：总动量)(PE为外压强，平衡时，PE =P)VPEHEVPEHE若系统与外界绝热，则 ，即（等焓过程）VPEE0H要保持压强不变，体积则必须变化，要保持压强不变，体积则必须变化，为此，将体积 V 作为一个新的动力学变量，构造一个新的拉氏函数：VPVMUmUTUTVVE iiVV22 21)(21),(&rrrrL引入标度坐标：或 LrVriii/3/1iiLr（，体积变化缓慢）iiii iiLLLmpvr&)(0L&， （共轭动量） iiiiimLrmr& L2VMV&VPVMLUmLVVE ii222 21)(21),(&L系统的哈密顿量为：（2-5）VPVMLUmLVVHE ii222 21)(21),(&运动方程为：（对粒子） （2-6）iiHmLdtd )(2&（对活塞） （2-7）VHVMdtd )(&将（2-5）代入（2-6）得i iiiiiiLFr rUUmLLmL 22&LL LmFii i&2VVdtdV dVdV dtdV dtdLL&3/23/13/131 VV VV
