3_分子动力学方法_微正则系综

2020/1/301三、分子动力学方法赖文生清华大学材料系技科楼2315室2020/1/302§3.1常规等容方法(NVE系综)分子动力学模拟是一种用来计算一个经典多体体系的平衡和传递性质的方法。经典这个词意味着组成粒子的核心运动遵守经典力学定律。真实实验：准备试样，将试样同一测试仪器连结，在一定的时间间隔内测量一下感兴趣的性质。如果测定服从统计噪音，那么测定的次数越多，结果就越准确。分子动力学模拟：准备试样：选择一个由N个粒子组成的模型体系，解体系的牛顿运动方程直到体系的性质不再随时间改变（体系达到平衡）。平衡以后，就进行实际的测定。模拟和实验容易犯的错误很相似：如试样未准备好，测定时间太短，体系在实验过程中经历一不可逆的变化，或者并未测定相要测定的量。2020/1/303温度在分子动力学模拟中，测定一个可观测的量，必须能将这一量表达为体系中粒子的位置与动量的函数。例如，对温度的方便定义是采用在所有自由度上的等配分能，以使它们以二次型进入哈密尔顿函数。实际过程中，温度是由测量体系的总动能除以体系的自由度数Nf(对总动能固定的体系Nf=3N-3，N是粒子数)得到。由于体系总动能涨落,体系的瞬时温度也随之改变温度的相对涨落的数量级在左右。因此要达到温度的一个准确的估计,应对许多涨落进行平均才行。TkmvB21212NifBiNktmvtT12)()(fN/12020/1/304程序对分子动力学模拟的最好介绍就是考察一个程序。程序按以下方式构成。1）读入指定运算条件的参数（如初始温度、粒子数、密度、时间等）。2）体系初始化（即选定初始位置和速度）。3）计算作用于所有粒子上的力4）解牛顿运动方程。这一步和上一步构成了模拟的核心。重复这两步直至我们计算体系的演化到指定的时间长度。5）中心循环完成之后，计算并打印测定的平均值，模拟结束。2020/1/305虚拟算法Programmd简单的MD程序callinit初始化t=0dowhile(t.lt.tmax)MD循环callforce(f,en)计算力callintegrate(f,en)积分运动方程t=t+deltcallsample抽样平均enddostopend2020/1/306初始化3.1.1初始化要启动模拟，应对体系中所有的粒子赋予初始位置和速度。所选粒子位置应与要模拟的结构相容。通常，将粒子初始放在立方晶格上来完成。然后对每一个粒子的分速度赋予一个区间(-0.5,0.5)上随机分布的值。接着，改变所有粒子的速度以确保总动量为零，并标定速度以使平均动能为指定值。热平衡时有式中v是给定粒子速度的分量。瞬时温度T(t):显然,通过因子fs=[T/T(t)]1/2来调整所有粒子的速度以使得瞬时温度与指定温度相一致。mTkvB/2NifiBNtmvtTk12,)()(2020/1/307分子动力学的初始化subroutineinitMD程序的初始化sumv=0sumv2=0doi=1,npartx(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)速度的标定因子doi=1,npartv(i)=(v(i)-sumv)*fs设定所需的动能和质心速度零xm(i)=x(i)–v(i)*dt先前时间步的位置enddoreturnend2020/1/308力的计算3.1.2力的计算我们以Lennard-Jones体系为例，来计算粒子间的互作用力。第i个原子受的力为：如果粒子对i，j足够近以至可发生作用的话，必须计算它们之间的力及其对势能的贡献。假设希望计算x方向上的力NiNijijijNiNijijtotrrruE16121)11(4)()5.01(48)(,)()()(6122rrrxrfrrurxrrurfxx得到NiNijijijijijitottotiidrrdurrrEEf1)(2020/1/309力的计算subroutineforce(f,en)确定力和能量en=0doi=1,npartf(i)=0设置力为零enddodoi=1,npartdoj=i+1,npartxr=x(i)-x(j)xr=xr–box*nint(xr/box)周期性边界条件r2=xr*xrif(r2.lt.rc2)then检验截断距离r2i=1/r2r6i=r2i**3ff=48*r2i*r6i(r6i-0.5)Lennard-Jones势2020/1/3010力的计算（续）f(i)=f(i)+ff*xr更新力f(j)=f(j)–ff*xren=en+4*r6i*(r6i-1)-ecut更新能量,ecut=4(rc-12-rc-6)endifenddoenddoreturnend2020/1/3011运动方程积分3.1.3运动方程积分采用Verlet算法，这一算法不仅是最简单的，而且通常是最好的。其推导可通过对在t时刻某粒子的坐标泰勒展开获得类似地，将这两个方程相加，可得或(3.1)新位置的估计值有一t4级别的误差，式中t是分子动力学模拟中的时间步长。注意Verlet算法并不使用速度来计算新的位置。)(!32)()()()(432tOrttmtfttvtrttr)(!32)()()()(432tOrttmtfttvtrttr)()()(2)()(42tOtmtftrttrttr2)()()(2)(tmtfttrtrttr2020/1/3012运动方程积分(续)将上面两个泰勒展开式相减,可以得到速度(可从轨迹知识得到)。或(3.2)速度的这一表达式的精度只在t2级别左右。由上可以看出，如果原子的前两个位置是已知的，则根据式(3.1)和式(3.2)就可进行相关的计算。然而，在典型的分子动力学模拟中，只有初始位置和初始速度是给定的。所以，在利用递归的Verlet算法开始计算之前，首先要找到原子第二个位置的计算方法。一般地，应用普通运动学方程可得)()(2)()(3tOttvttrttr)(2)()()(2tOtttrttrtv)(2)()()(2tftttvtrttriiii2020/1/3013积分动动方程subroutineintegrate(f,en)积分动动方程sumv=0sumv2=0doi=1,npartMD循环xx=2*x(i)-xm(i)+dt*dt*f(i)/2Verlet算法(3.1)vi=(xx-xm(i))/(2*dt)速度(3.2)sumv=sumv+vi质心速度sumv2=sumv2+vi*vi总动量xm(i)=x(i)更新先前步的位置x(i)=xx更新当前步的位置enddotemp=sumv2/(3*npart)瞬时温度etot=en+sumv2/2总能量returnend2020/1/3014§3.2等压方法（Isobaricmethod）上节我们讨论的分子动力学是一种研究在体积V中N个粒子的经典体系的自然的时间演绎的方法。在这些模拟中，总能量E是运动的常数。若我们假设时间平均等效于系综平均，那么在一般的MD模拟中得到的时间平均等价于在微正则系综（恒NVE）的系综。在其他系综（如NVT或NpT）中进行模拟更方便，更接近实际情况。乍一看，在不是微正则系综的其他系综中进行MD模拟是不可能。但已提出了两种相当不同的解决这一问题的方法。一种是基于这样一种思想，即其他系综的动力学模拟可以通过混合牛顿力学的MD和一定的MonteCarlo移动来实现。另一种方法完全源自动力学，它基于对体系拉格朗日运动方程的重新推导。2020/1/3015经典粒子服从牛顿力学、分析力学？*牛顿力学可以解决经典粒子的运动问题：),...,2,1(,NiFzmFymFxmiiiziiyiixii牛顿力学用矢量形式建立运动基本方程（矢量力学）。在处理约束系统、变形体的动力学等问题时有些困难。分析力学（标量力学）属于牛顿力学的范畴，但是又高于牛顿力学。ripHqqHpiiii,...,2,1*分析力学也可以解决经典粒子的运动问题（哈密顿方程）：qi是粒子i的广义坐标，pi是与之共轭的广义动量。2020/1/3016•牛顿力学用力、加速度等矢量表示质点运动规律，着眼于体系内每个质点受到的力及其产生的效果。•分析力学通过动能、势能等标量表示运动规律，并能将着眼点由质点转为体系整体（能量是标量，是体系的整体性质）。牛顿力学、分析力学？＝动能－势能),,,...,,,(),,,...,,,(),,,...,,,;,,,...,,,(111111111111NNNNNNNNNNNNzyxzyxVzyxzyxTzyxzyxzyxzyxLL),...,2,1(,NizLzLdtdyLyLdtdxLxLdtdiiiiii对于保守力学体系，力可以从势函数得到。如果定义一个Lagrange函数如下：),...,2,1(,NiFzmFymFxmiiiziiyiixii2020/1/3017Lagrange运动方程的优点是：不紧紧局限于笛卡儿坐标系，可以推广到任意坐标系。为此，引入广义坐标、广义速度Lagrange运动方程),...,2,1(,NizLzLdtdyLyLdtdxLxLdtdiiiiii广义坐标：足以能够确定体系内质点位置的任意一组独立参量。用符号qi(i=1,2,…,f）表示。有时候也称为独立坐标。广义坐标的数目与体系的自由度相等。所谓完整约束是指约束方程中不存在速度的约束:自由度的数目取决于体系受到的独立的完整约束的个数（f＝3N－m），因为每个独立的完整约束都会降低1个自由度。mitrrrFNi,2,10);,,(21Nitqqqrrfii,2,1);,,(212020/1/3018广义速度是广义坐标对时间的微商：Lagrange运动方程dtdqqii),...,2,1(,fiqLqLdtdii),...,2,1(fiqLpii，完整保守体系的力学运动状态可以用Largrange变量（时间、广义坐标、广义速度）完全描述。Hamilton提出，同样可以用时间、广义坐标及其共轭广义动量描述。Lagrange运动方程可以用广义坐标、广义速度表示：共轭广义动量的定义：f,...,2,1ipHqqHpiiii2020/1/3019Hamilton（正则）运动方程VTqqTLqqTLqqLLqppqHfiiififiiiiifiiiii1111)(),(＝Hamilton函数的定义如下：),...,2,1(,,ifipqHqpHiii对于稳定完整系统，其Hamilton函数就等于体系的机械能。所谓“正则”（canonical）是指“变量的微商已经解出在方程的另一边”。对称性诠释？