分子动力学简介(简明)

1、 什么是分子模拟 分子模拟是在分子模型的基础上用计算机做实验，“计算机实验” 通过模拟微观粒子的运动来计算宏观性质温度压力黏度传递性质表面张力分子间的作用模型牛顿力学量子力学统计力学等 分子模拟具有理论和实验的双重性质 分子模拟不能完全取代实验理论实验模拟理论的正确性模拟参数的正确性模拟方法的选择理论的更新与自然界相比的准确程度尺度（米）时间（秒）10-910-710-510-310-1510-910-610-31电子模拟（量化计算，DFT）分子模拟（分子动力学，蒙特卡洛）颗粒方法流体力学原子结构薛定谔方程模拟电子云能量性质，化学键等信息 量子化学计算一般处理几个到几十个原子常见软件：GAUS

2、SIAN，NWCHEM等 密度泛函（DFT）可以算到上百个原子上百个原子常见软件：VASP 分子水平的模拟 以分子的运动为主要模拟对象 采用经验性的分子间作用函数模拟微粒之间的作用 一般情况下不考虑电子转移效应，因而不能准确模拟不能准确模拟化学成键作用 发展最早 1950s，Alder，劳伦斯利物默实验室，分子动力学模拟32个原子 1950s，Metropolis，洛斯阿洛莫斯实验室，蒙特卡洛模拟32个原子 分子级别的模拟应用的领域很广 广泛应用于化学，物理，生物，化工，材料，机械，治药等领域 简单易学 蒙特卡洛是一种优化方法通过蒙特卡洛算法来寻求能量最优点随机方法 通过系综平均来求取宏观性质

3、 模拟的是平衡状态平衡状态，不涉及时间效应（KMC除外） 优点是可以跨越时间因素，缺点是得不到有关时间信息的性质 可以模拟平衡状态，也可以模拟中间状态 可以获得有关时间的信息 受时间的限制，无法模拟缓慢过程分子体系（几百几亿）求解牛顿运动方程宏观性质( )dU rm rdr= -& & 同时考虑原子核的运动（牛顿力学）和电子的运动（量子力学） 能同时准确模拟物理作用和化学键作用 目前来说CPMD可以处理的体系还很小（几十个原子） 将分子基团（几个或者几十上百个原子）当成单个的微粒来处理 微粒之间的作用也是通过类似于分子动力学的未能函数来描述 可以模拟更长的时间跨度电子原子核原子

4、量子级别模拟分子级别模拟CG级别模拟分子动力学对势能函数的依赖性：所有从分子动力学计算出来得到的宏观性质最终都取决于势能模型( )dU rm rdr= -& &分子动力学的核心：牛顿运动方程势能（位能）模型：( )( )U rf r=rUr例：甲烷，某些惰性气体质点处理Ur方阱模型Ur阶梯模型键的振动键角扭矩分子内部各原子（基团）之间的范德华力、静电力一般要计算1-4（相隔超过两个键的原子或基团对）15432qqq分子之间的范德华力分子之间的静电力CCCHHHHHHHH 10根键长作用 18个键角作用 8个扭矩作用 27个范德华力作用 27个静电作用20()1exp() eu

5、lDa ll=- Morse类键长模型能量阱深参数键长平衡键长参数 胡克类键长模型20()()2ku lll=-键长平衡键长参数 胡克类键角模型20()()2ku lqq=-键角平衡键角参数0( ) 1cos()2NnnVunwwg=+-扭矩障碍参数扭动360度所经过的能量最低点的次数126126( )4U rrrsse骣=-桫Lennard-Jones模型Ur 混合规则ABAB1()2A BABsx ss=+A BABez e e=通常都取1 一般情况下只考虑点电荷之间的作用力 不考虑极化作用所带来的长程项的作用114ABNNi jijijqqUrpe=邋分子A的点电荷总数分子B的点电荷总数

6、 用有限的微观分子体系模拟实际宏观体系的必要手段3214342121计算周期边界条件下两个微粒之间的作用宏观体系微观体系主盒影子影子影子 Ewald方法 Particle Mesh Ewald 方法（PME) 计算量很大BAC( , , )( , )( , )( , )U A B CU A BU A CU B C+ 处理方法(2)(3)( )( , ,)totalijijikjkUUrUr r r=+邋邋 通过量子化学模拟回归得到点电荷范德华力键长、键角、扭矩力 实验数据回归键长键角范德华力初始化能量优化平衡数据产出避免局部分子重叠，并不是动力学模拟 根据所有分子的当前坐标计算个分子的受力（位

7、能函数） 根据受力更新分子的坐标 在此过程中收集用来计算宏观性质的有关信息读入模型参数，模拟控制参数去除某些可能存在的原子重叠去除某些严重扭曲的键长、键角、扭矩等方法最速下降法牛顿拉夫森方法其他一般优化几千到几万步 Verlet法简单易行，但是有精度损失Leap-frog法，Verlet法的变种 速度Verlet法和Verlet法相比，可以同时获得位置、速度、加速度，而且没有精度损失Beeman算法，速度计算精度更高，但是计算量大预测-校正法算法选择精度和运算量的折衷体系原子体系刚性分子体系非刚性分子，刚性键非刚性分子涉及到的分子运动平动平动、转动平动、转动、扭矩平动、转动、扭矩、振动建议步长

8、10 fs5 fs2 fs-1 fs 不同性质的作用力采取不同的步长 最常见的方法：r-RESPA 成功用于多种体系，简单流体，有机分子体系，蛋白质等 和普通速度Verlet法相比，对于复杂分子体系（高分子，生物分子等），r-RESPA可以大大加快运算速度，步长可以达到1-2 fs 为保证原子各自运动时分子的整体性而添加的分子内部的约束条件算法 最常见的是针对氢原子，因为氢原子很高的振动频率 Verlet算法：SHAKE 速度Verlet算法：RATTLE NVE，最简单的分子动力学E=K+U恒定体积，而且和外界没有能量或者物质交换，只有动能和势能之间的转换通过监视各能量项的变化，可以检查程序

9、是否基本正确21/2(3)NiicBimKNNk T=-p p体系总动能和体系温度的关系原子总数体系总的约束条件数自由度模拟开始时，体系中各微粒的初始速度也是可以根据这个关系计算微粒 i 的动量微粒 i 的质量BPVN k TW=+压力与分子间的作用势能的关系理想气体分子间的作用造成的剩余压力1111( )33iNNiiiiiiWU= -籽=邋r rrrr frrr f维利系数W的计算扩散系数二元径向分布函数2( )()NNijij iVg rNd=-邋rrrr1( )(0)6iiDtt=-rrrrxtx通过斜率求取D实际位置，不是PBC后的位置 速度调节 简单scale v=v*(T0/T(

10、t)1/2 Berendsen法 增加调节参数 重新随机设置速度每隔一定步数，按照控制温度随机抽取一组高斯速度分布热量交换Nose-Hoover法将交换热源也当成体系的一部分进行积分产生严格的NVT系综控制体积活塞类似于温度控制通过调节体系的体积来调节压力压力的控制比温度更难，因为压力和能量的微分量（力）成正比，每一步的波动更大 气体的模拟 小分子体系，不需要复杂的势能模型 几百到几千个分子，分子分布稀疏，大部分是短程作用 一般用一台微机就可以处理，计算时间几分钟几小时 简单的液体，不涉及太多的界面性质 小分子体系，势能模型不是很复杂 几百个分子，可能涉及到静电作用，可能需要长程校正 用微机也

11、可以处理，计算时间一般几小时几天 必要性 体系越来越大 模拟时间越来越长 解决办法 制造更快的处理器 并行计算机 例子：50000原子的生物体系，1ns模拟 单个处理器：12天 16个并行处理器：1天或者 Message Passing Interface 90年代初制定和完善的一套并行语法 支持Fortran, C, C+ 简单易学 并行效率需要1小时需要1/2小时完美的并行效率 处理器的速度远远超过数据传输的速度，大量的时间花在处理器之间的信息传递上了 CPU的速度几乎是几何级数增长 内存的速度是代数级数增长 加快数据传输，尽量减少花在数据传输上的时间 数据传输硬件上的进步 算法上做文章

12、每一个处理器负责处理一部分原子 每一步计算每一个处理器都要接受其它处理器负责处理的原子的相关信息 信息传输量大，使用的处理器越多并行效率效率越低 一般适合处理5-10万左右微粒的体系N=20n=15n=610n=1620n=11151201481291961613155271741131018实际情况 按照体系的实际物理位置按区域划分每个处理器的处理范围 每一步计算每一个处理器只需要和相邻的处理器交换信息 数据传输量小，并行效率高，适合处理大型体系（超过10万微粒） 算法比较复杂（边界的处理）N=201201481291961613155271741131018 主要针对与生物和化学软材料体系

13、 优点 程序设计水平高，计算效率高，号称可以有效并行到上千个处理器 兼容多种输入和输出文件格式，有很好的分析辅助软件VMD 有很好的维护服务 不需安装 免费 缺点 万一需要自己安装的话比较麻烦/Research/namd/ 主要针对生物体系，也适当兼容一般化学分子 优点 有很好的内置势能模型 自定义新模型和新分子很方便 有很完善的维护网站 缺点 计算效率不高（收敛到16个处理器），运算速度慢 $400 主要针对生物体系，也包含部分化学体系 优点 势能模型更新很快 自定义新模型比较方便 维护服务很好 缺点 运

14、算速度慢，计算效率低 $600/ 一般性分子动力学软件，对生物体系略有偏重 优点 支持多种模型 免费 缺点 仍在开发中，某些方面还不完善/tinker/ 一般性分子模拟软件 优点 兼容当前大多数的势能模型 编程水平高，计算效率高（比NAMD差，强于其他所有类似软件） 可以模拟软材料和固体物理系统 免费 缺点 维护差/sjplimp/lammps.html 一般性分子模拟软件 优点 界面友好 计算效率高（有两个版本供选择，适合于不同大小的体系） 维护服务很好 缺点 兼

15、容性不好 100英镑http:/www.cse.clrc.ac.uk/msi/software/DL_POLY/ 主要针对生物体系，也适当照顾一般化学体系 优点 算法好，计算效率高 界面友好 维护服务好 免费软件 缺点 兼容性不好/纳米颗粒一般是人工合成物质纳米管，纳米球，纳米棒等等纳米颗粒在生物医学方面的应用(PNAS, v99, p6645, 2002)药物输送，基因输送 纳米颗粒的生物安全性 (Science, v300,p243, 2003)C60在鱼类脑部的累计 (Environmental Health Perspectives, v112

16、, p1058, 2004)C60导致老鼠的肺部病症 (Toxicol. Sci. v77, p117, p126, 2004)用分子动力学研究C60和基因分子的作用MD引擎 NAMD静电力-PME恒温 (300 K)衡压 (1 bar)时间步长 = 2 fs离子-中和体系模拟时间 4-20 nsDNA模型 AMBER 99 双链, 12-20 碱基对 单链, 12-20 碱基C60, sp2 碳原子用LJ 模型模拟溶剂 4600到6000水分子 TIP3P模型水，离子水，离子水，离子水，离子三个方向施加周期边界条件1 ns 的动画起始结构 C60和DNA能相互吸引，形成稳定的杂化结构 二者的

17、最初作用一般发生在1-2纳秒之内C60和B-DNA之间的距离和时间的函数C60和A-DNA之间的距离与时间的函数 C60和DNA分子之间典型的结合能大约在15-40 kcal/mol, 取决于二者的相对作用位置和结合点 C60和DNA分子之间的结合能远大于热波动kT (0.6 kcal/mol) 同时也远大于在同等模拟环境下C60和C60之间的结合能 (7.5 kcal/mol)C60吸附于DNA侧面C60吸附于DNA尾部 C60和B-DNA形成稳定的结合体，一般的结合能在20-40kcal/mol 比较稳定的结合部位在DNA的尾部或者侧面（minor grooves）C60进入DNA尾部 C

18、60和A型DNA的侧面和尾部结合，而且，C60可以进入A型DNA的尾部如果没有C60吸附，尾部的碱基对形成氢键的原子对之间的距离 穿入DNA尾部的C60破坏了DNA尾部的碱基对的氢键C60进入DNA尾部后，尾部的碱基对形成氢键的原子对之间的距离变形的角度 吸附在A型DNA侧面的DNA使得DNA碱基对之间的相对结构变形。本应平行的碱基对变形形成了一个大约20度的角 C60和单链DNA之间的结合导致单链DNA的严重变形 C60被单链DNA的某些部分包围，疏水的碱基平面和C60紧密接触 C60和单链DNA的结合能大约在20到30 kcal/mol C60和DNA之间结合的驱动力是疏水作用DNA亲水部位的水密度分布C60和DNA结合部位的水密度分布C60进入损伤部