1、三、分子动力学方法 赖文生清华大学材料系 技科楼 2315室* 13.1 常规等容方法 (NVE系综 ) 分子动力学模拟是一种用来计算一个经典多体体系的平衡和传递性质的方法。经典这个词意味着组成粒子的核心运动遵守经典力学定律。真实实验：准备试样，将试样同一测试仪器连结，在一定的时间间隔内测量一下感兴趣的性质。如果测定服从统计噪音，那么测定的次数越多，结果就越准确。分子动力学模拟：准备试样：选择一个由 N个粒子组成的模型体系，解体系的牛顿运动方程直到体系的性质不再随时间改变（体系达到平衡）。平衡以后，就进行实际的测定。模拟和实验容易犯的错误很相似：如试样未准备好，测定时间太短，体系在实验过程中经
2、历一不可逆的变化，或者并未测定相要测定的量。Date 2温度n 在分子动力学模拟中，测定一个可观测的量，必须能将这一量表达为体系中粒子的位置与动量的函数。例如，对温度的方便定义是采用在所有自由度上的等配分能，以使它们以二次型进入哈密尔顿函数。n 实际过程中，温度是由测量体系的总动能除以体系的自由度数Nf(对总动能固定的体系 Nf=3N-3， N是粒子数 )得到。由于体系总动能涨落 ,体系的瞬时温度也随之改变n 温度的相对涨落的数量级在 左右。因此要达到温度的一个准确的估计 ,应对许多涨落进行平均才行。 Date 3程序对分子动力学模拟的最好介绍就是考察一个程序。程序按以下方式构成。1）读入指定
3、运算条件的参数（如初始温度、粒子数、密度、时间等）。2）体系初始化（即选定初始位置和速度）。3）计算作用于所有粒子上的力4）解牛顿运动方程。这一步和上一步构成了模拟的核心。重复这两步直至我们计算体系的演化到指定的时间长度。5）中心循环完成之后，计算并打印测定的平均值，模拟结束。Date 4虚拟算法Program md 简单的 MD程序call init 初始化 t=0do while (t .lt. tmax) MD循环call force(f,en) 计算力call integrate(f,en) 积分运动方程 t = t + deltcall sample 抽样平均 enddostopen
4、dDate 5初始化3.1.1 初始化要启动模拟，应对体系中所有的粒子赋予初始位置和速度。所选粒子位置应与要模拟的结构相容。通常，将粒子初始放在立方晶格上来完成。然后对每一个粒子的分速度赋予一个区间 (-0.5,0.5)上随机分布的值。接着，改变所有粒子的速度以确保总动量为零，并标定速度以使平均动能为指定值。热平衡时有式中 v是给定粒子速度的 分量。瞬时温度 T(t):显然 , 通过因子 fs = T/T(t)1/2来调整所有粒子的速度以使得瞬时温度与指定温度相一致。Date 6分子动力学的初始化subroutine init MD程序的初始化sumv=0sumv2=0 do i=1, npa
5、rtx(i)=lattice_pos(i) 将粒子放在晶格上v(i) = (ranf( ) -0.5) 给定随机速度 sumv=sumv + v(i) 质心速度sumv2=sumv2 + v(i)*v(i) 动能 enddosumv = sumv/npart 质心速度sumv2 = sumv2/npart 均方速度fs = sqrt(3*temp/sumv2) 速度的标定因子do i = 1, npartv(i) = (v(i)-sumv)*fs 设定所需的动能和质心速度零 xm(i) = x(i) v(i)*dt 先前时间步的位置enddoreturnend Date 7力的计算3.1.2
6、力的计算我们以 Lennard-Jones体系为例，来计算粒子间的互作用力。第 i个原子受的力为：如果粒子对 i， j足够近以至可发生作用的话，必须计算它们之间的力及其对势能的贡献。假设希望计算 x方向上的力Date 8力的计算subroutine force(f,en) 确定力和能量en = 0do i=1, npartf(i)=0 设置力为零enddodo i=1, npartdo j = i+1, npartxr = x(i)-x(j)xr = xr box*nint(xr/box) 周期性边界条件 r2 = xr*xr if ( r2 .lt. rc2) then 检验截断距离 r2i = 1/r2r6i = r2i*3ff = 48*r2i*r6i(r6i-0.5) Lennard-Jones势 Date 9力的计算（续）f(i) = f(i) + ff*xr 更新力 f(j) = f(j) ff*xren = en + 4*r6i*(r6i-1)-ecut 更新能量 , ecut = 4(rc-12-rc-6)endifenddoenddoreturnend Date 10
