计算物理：chapter3.2

1、分子动力学模拟分子动力学模拟(II)Molecular Dynamics Simulation基本思路基本思路根据粒子的运动规律根据粒子的运动规律( (如牛顿方程）对系统中的每个如牛顿方程）对系统中的每个粒子的运动方程进行数值求解，得到各个时刻每个粒子的运动方程进行数值求解，得到各个时刻每个粒子的坐标和动量粒子的坐标和动量( (即相空间中的运动轨迹即相空间中的运动轨迹) )，再利，再利用统计计算方法得到体系的静态和动态特征。用统计计算方法得到体系的静态和动态特征。基本步骤基本步骤1.初始化初始化2.求力求力3.积分运动方程积分运动方程4.计算物理量计算物理量Program MDCall ini

2、t ！初始化！初始化t=0Do while (t.lt.tmax) ！MD循环循环 call force(f,en) ！计算力！计算力 call integrate(f,en) ！积分运动方程！积分运动方程 t=t+deltt call sample ！抽样平均！抽样平均EnddoStopend积分运动方程积分运动方程vverletvLeapfrogvVelocity verletvGear Predictor-Corrector已知初始速度，坐标和力，通过数值求解运动方程已知初始速度，坐标和力，通过数值求解运动方程来求各时刻每个粒子的坐标和速度。来求各时刻每个粒子的坐标和速度。 已知已知t-

3、t和和t时刻的位置时刻的位置 由由t时刻位置，可得时刻位置，可得t时刻的受力时刻的受力 由由t时刻的受力及时刻的受力及t-t和和t时刻的位置，可得时刻的位置，可得t+t时刻的位置时刻的位置 由由t+t及及t-t时刻的位置，可得时刻的位置，可得t时刻的速度时刻的速度 已知已知t时刻位置、速度时刻位置、速度 由由t时刻的位置，可得时刻的位置，可得t时候受力时候受力 由由t时刻的受力和位置、速度，可得时刻的受力和位置、速度，可得t+t的位置和受力的位置和受力 由由t+t的受力和的受力和t时刻的速度和受力，可得时刻的速度和受力，可得t+t时刻的速度时刻的速度 已知已知t-t和和t时刻的位置时刻的位置

4、由由t-t和和t时刻的位置，可得时刻的位置，可得t-t/2时刻的速度和时刻的速度和t时刻的受力时刻的受力 由由t时刻的受力和时刻的受力和t- t/2时刻的速度，可得时刻的速度，可得t+t/2时刻的速度时刻的速度 由由t+t/2时刻的速度和时刻的速度和t时刻的位置，可得时刻的位置，可得t+t时刻的位置时刻的位置Subroutine Integrate(f,en) 积分运动方程积分运动方程(微正则系综微正则系综)Sumv=0Sumv2=0Do I=1,npart MD 循环循环xx=2*x(i)-xm(i)+dt*f(i) verlet 算法算法vi=(xx-xm(i)/(2*dt) 速度速度 s

5、umv=sumv+v(i) 质心速度质心速度 sum2=sumv2+v(i)*2 总动能总动能xm(i)=x(i) 更新上一步位置更新上一步位置x(i)=xx 更新当前步位置更新当前步位置EnddoTemp=sumv2/(3*npart) 瞬时温度瞬时温度Etot=en+sumv2/(2*npart) 总能量总能量ReturnendProgram MDCall init ！初始化！初始化t=0Do while (t.lt.tmax) ！MD循环循环 call force(f,en) ！计算力！计算力 call integrate(f,en) ！积分运动方程！积分运动方程 t=t+deltt c

6、all sample ！抽样平均！抽样平均EnddoStopend平衡态分子动力学模平衡态分子动力学模拟拟计算物理观测量计算物理观测量热力学热力学 ( (系综系综) )平均：某一时刻系宗中所有系统的平均平均：某一时刻系宗中所有系统的平均遍历性假设：遍历性假设：对无限长的轨迹，热力学平均动力学平均对无限长的轨迹，热力学平均动力学平均 或或 时间平均系综平均时间平均系综平均动力学平均动力学平均( (时间时间) )平均平均: :某一系统在所有时间的平均某一系统在所有时间的平均遍历性假设遍历性假设在在MD模拟中，物理量的时间平均表示为：模拟中，物理量的时间平均表示为： Time average:NNN

7、NNNNNNNensdrdprpHdrdprpHrpAA),(exp),(exp),(等温过程模拟等温过程模拟Andersen热浴热浴Velocity-scaling基本思路基本思路:粒子和虚粒子相互碰撞粒子和虚粒子相互碰撞连续两次随机碰撞的时间间隔满足泊松分布：连续两次随机碰撞的时间间隔满足泊松分布： Andersen热浴热浴tvvevtPvt/1;),(粒子碰撞后的速度由粒子碰撞后的速度由Maxwell分布抽样获得，分布抽样获得， Maxwell分布的对应温度为目标温度分布的对应温度为目标温度Andersen, JCP,72,2384(1980)T实现步骤：实现步骤：1）从初始位置、动量）

8、从初始位置、动量rN(0),pN(0)开始，并积开始，并积 分运动方程。分运动方程。2）根据一个粒子在一个时间步长）根据一个粒子在一个时间步长 t内经历一次内经历一次随机碰撞的几率为随机碰撞的几率为vdeltat，选择被碰撞的粒子。，选择被碰撞的粒子。3)被碰撞粒子的速度根据对应于温度被碰撞粒子的速度根据对应于温度T的的Maxwell分布随机抽取。分布随机抽取。Program MD_andersenCall init(temp) ！初始化！初始化Call force(f,en)T=0Do while (t.lt.tmax) ！MD循环循环 call integrate(1,f,en) ! 速度

9、速度verlet第一步第一步 call force(f,en) ！计算力！计算力 call integrate(2,f,en) ！速度！速度verlet第二步第二步 t=t+deltt call sample ！抽样平均！抽样平均EnddoendSubroutine integrate(switch, f,en,temp)If(switch.eq.1)then 速度速度Verlet法第一步法第一步Do I=1,npartX(i)=x(i)+dt*v(i)+dt*dt*f(i)/2 更新当前位置更新当前位置V(i)=v(i)+dt*f(i)/2 第一次更新速度第一次更新速度EnddoElse i

10、f(switch.eq.2)then 速度速度Verlet法第二步法第二步Do I=1,npartV(i)=v(i)+dt*f(i)/2. 第二次更新速度第二次更新速度EnddoDo I=1,npart Anderson热浴热浴If(rand().lt.nu*dt)v(i)=gauss(temp) 与热浴碰撞，给出粒子速度与热浴碰撞，给出粒子速度Enddo EndifReturnendVelocity-Scaling:32kBENK T粒子的动能：粒子的动能：23iiiBmvTNkiiikvmE221oldnewTT2oldnewvv2/1tanarg/tinsettTT等动等动能能Beren

11、dsen thermostat (weak coupling)：Berendsen et al. JCP, 81:3684(1984)(1)(0TTdttdT)(0TTtT1102TTt)()(tvtv等压过程模拟等压过程模拟Volume-scaling:)(1)(tPPdttdPbathP)(1bathPPPtinewirr3/1 jiijijkFrEVP2132其中：其中：压强由压强由Virial公式给出公式给出Berendsen et al. JCP, 81:3684(1984)分子动力学模拟应用举例分子动力学模拟应用举例v二维简单体系二维简单体系v生物大分子生物大分子二维简单体系二维简

12、单体系(Argon)6124)(rrrV相互作用：相互作用：单位：单位：长度：长度： ； 时间时间:(m 2/ )1/2 (对对Ar, 1.8x10-12s)速度速度:( /m)1/2 温度：温度： /kB (对对Ar, 120K)TVNT=0T=0.2T=0.4T=0.6T=0.8T=1.0Melting Melting 过程过程T=1.2T=1.4T=1.6T=2.0T=3.0位置涨落位置涨落作业：作业：三维简单体系三维简单体系(6x6x6)，计算，计算melting过程。过程。1. 优化结构优化结构2. Melting温度的密度依赖性温度的密度依赖性( (要求：写出计算程序要求：写出计算

13、程序(C(C或或FORTRAN)FORTRAN)，给出计算结果，给出计算结果) )DNA sequenceProtein sequence3D structureMolecular recognition例：蛋白质系统的分子动力学模拟例：蛋白质系统的分子动力学模拟初始结构初始结构PDBPDB数据库：数据库： X-rayX-ray NMR NMR软件建模软件建模(hyperchem)(hyperchem)力场参数力场参数kb is the spring constant of the bond.r0 is the bond length at equilibr

14、ium.Unique kb and r0 assigned for each bond pair, i.e. C-C, O-Hk is the spring constant of the bend.0 is the bond length at equilibrium.Unique parameters for angle bending are assigned to each bonded triplet of atoms based on their types (e.g. C-C-C, C-O-C, C-C-H, etc.)A controls the amplitude of th

15、e curven controls its periodicityf shifts the entire curve along the rotation angle axis (). The parameters are determined from curve fitting.Unique parameters for torsional rotation are assigned to each bonded quartet of atoms based on their types (e.g. C-C-C-C, C-O-C-N, H-C-C-H, etc.)A determines

16、the degree the attractivenessB determines the degree of repulsionq is the chargeC 12.01 0.616 mass pol.N 14.01 0.530C -C 310.0 1.525 k lengthC -O 570.0 1.229C -C -O 80.0 120.00 k angleX -C -O -X 4 11.20 180.0 2.00 divf Pk phi N 溶剂模型溶剂模型Alpah-helix折叠分子动力学模拟折叠分子动力学模拟Rc L/2The fraction of native contac

17、tsSimulation time (nsec)030020010040000.81.0常用的分子动力学模拟软件常用的分子动力学模拟软件AMBERCHARMMGROMACSTINKERNAMDAmber Amber是程序套件集合的名称，内含大约是程序套件集合的名称，内含大约60个程序个程序 现已有现已有Amber11 网址网址 / 主要程序有：主要程序有： sander leap antechamber ptraj和和carnal pmemd nmode mm_pbsa Amber 提供可视化软件包提供可视化软件包, 如如MOIL-View，

18、MD display, Interactive Essential DynamicsAmber 给其他程序提供参数支持给其他程序提供参数支持 AmberFFC: Using the Amber force fields in Accelrys programs Using the Amber force field in NAMD Using the Amber force field in CHARMM Using the Amber force field in GROMACS CHARMm CHARMm：Chemistry at HARvard Macromolecular Mechanics 网址：网址：http:/