1、学习目标学习目标v 了解计算化学的主要方法和基本概念了解计算化学的主要方法和基本概念v 能够读懂相关领域的文献，了解常见分子模拟软件及能够读懂相关领域的文献，了解常见分子模拟软件及使用使用v 掌握基本的计算化学研究思路和方法。掌握基本的计算化学研究思路和方法。v 了解分子模拟能干什么，我们可以作那些工作，作到了解分子模拟能干什么，我们可以作那些工作，作到什么精度，为在自己的研究领域开展分子模拟方面的什么精度，为在自己的研究领域开展分子模拟方面的研究打一个基础。研究打一个基础。 分子模拟具有理论和实验的双重性质 分子模拟不能完全取代实验理论实验模拟理论的正确性模拟参数的正确性模拟方法的选择理论的
2、更新系综系综 系综（ensemble）代表一大群相类似的体系的集合。对一类相同性质的体系，其微观状态微观状态（比如每个粒子的位置和速度）仍然可以大不相同。（实际上，对于一个宏观体系，所有可能的微观状态数是天文数字。） 统计物理的一个基本假设（各态历经各态历经假设）是：对于一个处于平衡的体系，物理量的时间平均，等于对对应系综里所有体系进行平均的结果。 体系的平衡态的物理性质可以对不同的微观状态求和来得到。系综的概念是由约西亚威拉德吉布斯（J. Willard Gibbs）在1878年提出的。 常用系综常用系综v微正则系综微正则系综 (microcanonical ensemble)v正则系综正则
3、系综 (canonical ensemble) v巨正则系综巨正则系综 (grand canonical ensemble)v等温等压系综等温等压系综 (isothermal-isobaric ensemble)常用系综常用系综微正则系综微正则系综 (microcanonical ensemble)：系综里的每个体系具有相同的：系综里的每个体系具有相同的能量（通常每个体系的粒子数和体积也是相同的）。能量（通常每个体系的粒子数和体积也是相同的）。正则系综正则系综 (canonical ensemble)：系综里的每个体系都可以和其他体系交：系综里的每个体系都可以和其他体系交换能量（每个体系的粒子
4、数和体积仍然是固定且相同的），但是系综里换能量（每个体系的粒子数和体积仍然是固定且相同的），但是系综里所有体系的能量总和是固定的。系综内各体系有相同的温度。所有体系的能量总和是固定的。系综内各体系有相同的温度。巨正则系综巨正则系综 (grand canonical ensemble)：正则系综的推广，每个体系都：正则系综的推广，每个体系都可以和其他体系交换能量和粒子，但系综内各体系的能量总和以及粒子可以和其他体系交换能量和粒子，但系综内各体系的能量总和以及粒子数总和都是固定的。（系综内各体系的体积相同。）系综内各个体系有数总和都是固定的。（系综内各体系的体积相同。）系综内各个体系有相同的温度和
5、化学势。相同的温度和化学势。等温等压系综等温等压系综 (isothermal-isobaric ensemble)：正则系综的推广，体系间：正则系综的推广，体系间可交换能量和体积，但能量总和以及体积总和都是固定的。（系综内各可交换能量和体积，但能量总和以及体积总和都是固定的。（系综内各体系有相同的粒子数。）正如它的名字，系综内各个体系有相同的温度体系有相同的粒子数。）正如它的名字，系综内各个体系有相同的温度和压强。和压强。 分子动力学方法工作框图分子动力学方法工作框图给定t 时刻的坐标和速度以及其他动力学信息，那么就可计算出 t t 时刻的坐标和速度。程序构成方式 输入指定运算条件的参数（初始
6、温度，粒子数，密度，时间步长） ； 体系初始化（选定初始坐标和初始速度） ； 计算作用在所有粒子上的力； 解牛顿运动方程（第 和第 步构成了模拟的核心，重复这两步，直到体系的演化到指 定的时间） ； 计算并输出物理量的平均值，完成模拟。积分算法优劣的判据积分算法优劣的判据 分子动力学中一个好的积分算法的判据主要包括： 计算速度快； 需要较小的计算机内存； 允许使用较长的时间步长； 表现出较好的能量守恒。分子动力学的适用范围分子动力学方法只考虑多体系统中原子核的运动，而电子的运动不予考虑，量子效应忽略。经典近似在很宽的材料体系都较精确；但对于涉及电荷重新分布的化学反应、键的形成与断裂、解离、极化
7、以及金属离子的化学键都不适用，此时需要使用量子力学方法。经典分子动力学方法（MD）也不适用于低温，因为量子物理给出的离散能级之间的能隙比体系的热能大，体系被限制在一个或几个低能态中。当温度升高或与运动相关的频率降低（有较长的时间标度）时，离散的能级描述变得不重要，在这样的条件下，更高的能级可以用热激发描述。 电子运动具有更高的特征频率，必须用量子力学以及量子经典理论联合处理。这些技术近年来取得了很大进步。在这些方法中，体系中化学反应部分用量子理论处理，而其他部分用经典模型处理。QM/MM分子力学、分子动力学方法及其应用分子力学、分子动力学方法及其应用Molecular Mechanics /
8、Molecular Dynamics 一、 MM、MD理论基础 二、 MM、MD计算程序 三、 MM、MD方法的应用一、一、 MM、MD基础理论基础理论 原则上, 第一原理方法在理论上已经能解决所有问题 但计算量太大，计算机资源有限，原子数目较多时，如高分子、蛋白质、原子簇以及研究表面问题、功能材料或材料的力学性能等，实际上难以完成计算 为此，发展了分子力学(Molecular Mechanics, MM)与分子动力学(Molecular Dynamics, MD)方法 它们的应用，又称分子模拟(molecular simulation, molecular modeling) 或 分子设计(
9、molecular design) MM与MD是经典力学方法，针对的最小结构单元不再是电子电子而是原子 因原子的质量比电子大很多，量子效应不明显，可近似用经典力学方法处理 20 世纪 30 年代, Andrews 最早提出分子力学(MM)的基本思想；40 年代以后得到发展, 并用于有机小分子研究。90年代以来得到迅猛发展和广泛应用基本思想基本思想 事先构造出简单体系(如链段、官能团等各种不同结构的小片段)的势能函数, 简称 势函数 或 力场(force field) 将势函数建成数据库，在形成较大分子的势函数时，从数据库中检索到结构相同的片段，组合成大体系的势函数 利用分子势能随原子位置的变化
10、有极小值的性质，确定大分子的结构即为分子力学(MM) 利用势函数，建立并求解与温度和时间有关的牛顿运动方程，得到一定条件下体系的结构随时间的演化关系即为分子动力学(MD) 理论方法的核心核心是构造势函数 势函数：势能与原子位置的关系。且往往是不知道的 需要通过其他方法，如量子化学方法及实验数据获得rrrE过渡状态 E反应物产物分子势函数曲面势能面示意图 1、分子力场、分子力场 分子片段力场的函数表达式中包含自变量和力场参数 其中自变量为分子的结构参数，独立参数为键长、键角和二面角，如图 而 力场参数 一般通过与实验数据 和 从头算数据进行最小二乘法拟合来确定bbb有的还使用一个非独立参数：面外
11、弯曲角 势函数形式很多，目前已被广泛使用的力场有如CFF、MM2、MM3、MM4、MMFF、AMBER、CHARMM、DREIDING、UFF和COMPASS等 形式虽多, 但一般总表达为分子内与分子间势能之和:V总总=V键合键合+V非键合非键合 分子内势能(键合)包括键伸缩、键角弯曲和二面角扭转势能 分子间势能(非键合)包括范德华势和静电势, 有的还包括H键：V键合=V键伸缩+V键角弯曲+V二面角扭转V非键合=V范德华+V静电+V氢键 键合势函数中，一些力场还包含交叉项，使精度更高 交叉项的含义：如键长变化时，键角弯曲势能随键长的不同而不同，等 例例: COMPASS-98力场(conden
12、sed-phase optimized molecular potentials for atomistic simulation studies)的表达式如下 每个k是一独立的力场参数,下标“0”代表参考(平衡)结构参数：bbbbkbbkbbkE)()()(404303202键伸缩：)()()(404303202kkkE键弯曲：)3cos(1 )2cos(1 )cos(1 3 , 032, 021 , 01kkkE二面角：2kE键面外弯曲：,0000)()(bbbbxkbbbbkE交叉项：,3210003cos2coscos)()(bbbkkkbbbbk, 32103cos2coscos)(
13、bkkkbb,32103cos2coscos)(kkk,00)(cosk 力场参数力场参数k最小二乘法确定最小二乘法确定 基本思想 如：R-COOH基团 1) 由ab initio (构型优化方法)计算出 平衡结构，得到 参考结构参数bi0, i0, i0 2) 用伪随机数方法将bi0, i0, i0人为改变成若干(n)组非平衡结构参数 bi, i, il (l=1,2, , n) 3) 用各bi, i, il 结构参数分别进行ab initio计算，得到bi, i, il 结构参数下对应的能量El (l=1,2, , n) 4) 将El (l=1,2, , n)和bi, i, il代入势能表
14、达式 5) 用最小二乘法拟合，确定力场参数kCOOH 非键合势函数中，静电相互作用表示分子中各原子静电荷的库仑相互作用对势能的贡献 不同的力场，静电相互作用表达式基本相同 范德华势也大都采用Lennard-Jones函数，但函数中的指数有所不同。如COMPASS-98的非键合势函数为：静电相互作用：jiijjirqqEelec范德华势：jiijijijijijvdWrrrrE,609032即Lennard-Jones 9-6函数其它力场范德华势较多采用L-J 12-6函数：jiijijijijijvdWrrrrDE,6012002)( 当然，在建立分子的势能函数时，还有一些更细致的问题要考虑，
15、如： 势能展开项的截断 周期结构的处理 多组分混合物体系 含有离子的体系 或 金属中的离子 等 不同的方法或程序中，分别都有更详细的讨论。参考：德 D. 罗伯. 计算材料学. 北京：化学工业出版社, 2002, 9 俞庆森,朱龙观. 分子设计导论.北京：高等教育出版社,2000杨小震. 分子模拟与高分子材料. 北京：科学出版社，2002熊家炯主编. 材料设计. 天津：天津大学出版社，2000Sun H，Ren P, Fried J R. The COMPASS Force Field: Parameterization and Validation for polyphosphazenes.
16、Computational and Theoretical Polymer Science, 1998, 8(1/2): 229Sun H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications - Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem., 1998, 102: 7338 2、分子力学方法、分子力学方法 MM是确定分子结构的方法 利用分子势能随结构的变化而变化的性质，确定分子势能极小时的平衡结构(
17、stationary point) 物理模型: 视原子为质点，视化学键为弹簧，而弹力常数完全由数据库中的分子力场来确定 因此是直接用势函数研究问题，不考虑原子的动能 不考虑动能所对应的结构，相当于体系处于 T=0 K 时的结果 由力场首先构造并得到分子的势函数E(x) 利用在E(x)的极小点处，E(x)随各原子独立的空间坐标xi, i=1,2,3,3N-6 (与内坐标ri、i、i等价) 的一阶微分=0以及全部二阶微分0的数学条件：0)(0)(22iixxExxE(i=1,2,3,3N-6) 进行结构优化 具体步骤： 1) 用各种方法构造出一个分子的任意结构，得到初始的结构参数xi0, i=1,
18、2,3,3N-6 2) 进行坐标变换，即根据需要将原子的直角坐标转变成内坐标ri、i、i或反之 1) 用各种方法构造出一个分子的任意结构，得到初始的结构参数xi0, i=1,2,3,3N-6 2) 进行坐标变换，即根据需要将原子的直角坐标转变成内坐标(键长、键角、二面角)或反之 3) 建立分子体系的势能表达式E(x) 4) 计算E(x)随各坐标的一阶、二阶导数 5) 计算接近数学条件的坐标增量 6) 得到新的结构参数xi1, i=1,2,3,3N-6 重复 4)、5)、6)，直至最后两次得到的体系势能之差或总体势能梯度的均方根值达到预定精度范围为止 MM计算中，有时还要考虑体系所处的外压条件，
19、此时，可通过压力因子的定义来调节原子坐标0)( , 0)(22iixxExxE(i=1,2,3,3N-6) 3、分子动力学方法、分子动力学方法 分子力场是分子的静态势函数 而实际过程通常是在一定温度和一定压力下发生的 为了更切实际地了解体系运动和演化的过程，必须考虑体系中原子的运动，并与温度T和时间t建立联系 我们知道，温度是原子分子热运动剧烈程度的量度 根据统计热力学，对于n个原子的体系，体系的温度T与各原子的运动速率vi的关系为：niiivmnTk12B31 又因体系中各原子的速率为vi时，动量pi=mivi，对应总动能K(p)为：niziyixiipppmpK12,2,2,)(21)(
20、势能由力场确定为E(x)，因此体系的Hamilton量H为：)()(),(xEpKpxH 与量子力学不同，经典力学对Hamilton量不进行算符化处理，也不建立和求解本征方程，而是建立并求解经典运动方程：ppxHtxxpxHtp),(dd ),(dd 计算过程一般为：在一定的统计系综下 1) 由原子位置和连接方式，从数据库调用力场参数并形成体系势函数 2) 由给定温度计算体系动能以及总能量 3) 计算各原子的势能梯度, 得到原子在力场中所受的力 即dp/dt=mdv/dt=ma=F 1) 由原子位置和连接方式，调用力场参数并形成体系势函数 2) 由给定温度计算体系动能以及总能量 3) 计算各原
21、子的势能梯度, 得到原子在力场中所受的力 即dp/dt=mdv/dt=ma=F 4) 对每个原子，在一定时间间隔内，用牛顿方程求解其运动行为：ppxHtxxpxHtp),(dd ,),(dd 22ddddtxmtvmamFiiiiiii 5) 显示体系能量和结构 6) 取下一时间间隔，返回步骤 1) 不断循环反复，可设定循环次数或强行终止计算 其中, 不仅时间间隔可以根据需要取不同大小(一般1 fs=10-15 s),温度可以任意设定，而且还可以在循环过程中逐渐改变温度，即研究体系的退火(annealing)行为 有时需要进一步考虑外场的作用，如压力、电场、磁场、重力场等 从原则上讲，这些问题
22、都不难解决，而且还在进一步发展中，不再赘述 此外，即使是分子力学与分子动力学方法，也受计算量的限制，所处理的体系不可能太大，好的计算机可达到数十万百万个原子的规模 (2050 nm以下/ MD 106108个原子罗伯p110) 对于更大(50 nm、微米)尺度的问题，人们也在发展介观尺度(mesoscale)的方法，并且已取得一些成果 分子动力学方法也有一些本质的缺陷 如：势函数精度(客观性) 势函数形式在每次计算中都不变, 故不能模拟如分子在高温下, 结构发生断裂的热裂解过程 (国内引入开关函数已开始做) 此时, 可用QM处理化学变化或e转移, 余用MM，称QM/MM二、二、 MM、MD计算
23、程序计算程序 因程序结构相对简单，国内外有不少, 且还在不断产生 如：Chem3D, ChemOffice, Cerius2, Material Studio, Alchemy 2000, Sybyl, Biosym, HyperChem x.x, Spartan x.x, Chemgraf, Bilder, Script, COGS, Gaussian-03; MP(Molecular Properties, 杨小震)等 输入：输入：体系模型建模(Builder)，选择力场、系综 MM: 无特殊输入 MD: 循环次数(万几10万)、温度、时间间隔(fs), 外场,温度变化速率，其它性质 等
24、输出：输出： MM: 势能-结构曲线(数值)，动态结构图, 其它性质 MD: (同上)，及 其它MD性质 注：注：第一原理及MM一般只能得到势能极小点结构，而MD可越过一些小的势垒，甚至可达到最小点三、三、 MM、MD方法的应用方法的应用 领域：高分子、生命科学、药物设计、催化、半导体 其它功能材料、结构材料等 分子力学是用计算机在原子水平上模拟给定分子模型的结构与性质，进而得到分子的各种物理性质与化学性质，如结构参数、振动频率、构象能量、相互作用能量、偶极矩、密度、摩尔体积、汽化焓等 分子动力学方法能实时将分子的动态行为显示到计算机屏幕上, 便于直观了解体系在一定条件下的演变过程 MD含温度
25、与时间, 因此还可得到如材料的玻璃化转变温度、热容、晶体结晶过程、输送过程、膨胀过程、动态弛豫(relax)以及体系在外场作用下的变化过程等 分子模拟较早应用于高分子问题的研究 应用范围主要包括: 物理性质、结构、构象与弹性、晶体结构、力学性能、玻璃态与玻璃化转变、光谱性质、非线性光学性质、电性质、共混与分子间相互作用等 应用于生物科学和药物设计也十分普及，如蛋白质的多级结构与性质，病毒、药物作用机理、特效药物的大通量筛选与快速开发等 在化学领域，用于表面催化与催化机理、溶剂效应、原子簇的结构与性质研究等 在材料科学, 用于材料的优化设计、结构与力学性能、热加工性能预报、界面相互作用、纳米材料
26、结构与性能研究等罗伯，计算材料学，p129: 已有很多研究者采用分子动力学在原子尺度上进行了材料性质及结构的模拟 在材料科学领域，相变及晶格缺陷结构的研究一直是普遍关注的热点 在这方面，人们关注的焦点集中在同相和异相界面、位错、裂纹、界面偏析、失配(misfit)位错、次级晶界位错以及聚合物中的结构筹等问题 而且，原子尺度上的材料合成和设计已经取得实质性进展 也正是在这一点上，应力被分别地引入到金属、聚合物和陶瓷材料的具体研究当中 已发表的工作包括： 经典例子和基本原理、外界环境条件、失配位错的结构、界面结构及其迁移率、晶格位错的结构和能量、偏析效应、微裂纹结构、薄膜和表面、异相界面、材料设计、原子论和有限元讨论及大作业讨论题目：分子力学研究进展分子动力学研究进展常用计算化学软件使用介绍气体小分子在聚合物中扩散系数的计算举例说明分子模拟在化学领域的应用作业要求：（1）以论文综述的形式查阅并撰写有关计算机模拟（分子模拟）在化学领域中的应用进展。（2）要有参考文献的出处，避免抄袭拷贝现有文献综述，一定要自己归纳总结，提炼各综述的精华。（3）WORD电子版形式提交到： (4)文件名：题目-姓名-学号（5）气体小分子在聚合物中扩散系数的计算（要有计算模拟的主要流程图，分阶段描述，各阶段的目的）
