理论课讲义04.mm methods dftffmd分子模拟与设计

1、分子模拟与设计Given by:周成冈 Ph. D., Professor of Chemistry中国地质大学(武汉)材料与化学学院Faculty of Materials Science and Chemistry 2013 Spring基于波函数的计算基于波函数的计算方法Ab initio HF后HF方法Semiempirical原理基于量子力学(求y)基于量子力学(求y)基于量子力学(求y)相关能不考虑考虑不考虑算法严格从头算，无任何经验参数严格从头算，无任何经验参数大量经验参数和近似速度非常慢最慢非常快精度较为精确最精确定性基本准确可靠性非常可靠最可靠误差随意性大计算性质基态基态、激

2、发态基态、激发态(ZINDO/S)处理体系100401000密度泛函理论从头算(HF, p-HF, SE)：基于复杂的多电子波函数，存在3N个变量(N为电子数，每个电子包含3个空间变量)，体系能量为波函数的本征值。密度泛函理论(Density Functional Theory)：电子密度取代波函数 (电子密度为3个空间变量的函数)，体系的基态能量为电子密度的泛函(Hohenberg-Kohn第一定理)。DFT通过Kohn-Sham框架实现：复杂的多体问题被简化成一个没有相互作用的电子在有效势场中运动的问题。有效势包括所有的相互作用效应，即Hartree势、交换势(Pauli原理决定的相互作用

3、势)和关联势(某电子对整个电荷分布的影响所产生的势)。DFT典型程序DFT的地位与作用近十几年来，DFT同分子动力学方法相结合，有大量的新发展(First-Principles MD)在材料设计、合成、模拟计算和评价诸多方面得到了广泛应用，是目前最流行的计算方法(量化计算界至少70%以上的人采用)已成为计算凝聚态物理、计算材料科学和计算量子化学的重要基础和核心技术逐渐进入工业技术领域DFT的优点DFT vs. HF: 降低了计算体系的维数波函数电子密度提供了第一性原理或从头算的计算框架，在计算精度和计算量之间达到了较好的权衡在这个框架下可以发展各式各样的结构和能带计算方法，例如材料电子结构和几

4、何结构、固体和液态金属中的相变等。这些方法已发展成为用量子力学方法计算力的精确的分子动力学方法。DFT的应用范围电子基态能量与原子位置之间的关系可以用来确定分子或晶体的结构，从而得到各种性质当原子不处在它的平衡位置时，给出作用在原子上的力因此，DFT可以解决物理和化学中的大多数问题，如力、光、电、磁、热等电离势、电子亲和能振动光谱化学反应生物分子的结构催化活性位置的特性Hohenberg-Kohn定理相互作用多体系统的粒子数密度是决定系统物理性能的基本变量，即系统所有物理性质都是由粒子数密度决定的，包括总能量和系统态密度。Hohenberg-Kohn能量泛函形式：能量是密度的唯一泛函，仅当 取

5、严格的基态密度时能量泛函取极小值，并代表电子系统的基态密度。H-K能量泛函与变分方程其中T0: 无相互作用系统的动能，T是相互作用系统的动能VH: Hartree势EXC:交换关联能，包含了Hartree近似以外所有的多体效应变分方程：限制性条件是总粒子数守恒：Kohn-Sham自洽方程组将相互作用多体系统基态问题在形式上严格转化为在有效势场中运动的独立粒子的基态问题。处理KS方程中的交换关联项EXC, VXC的处理需要用到局域密度近似(Local Density Approximation)或广义梯度近似(Generalized Gradient Approximation)KS方程求解流程

6、图rin(r)r(r)=rat(r)求解、Vxc、Veff求解Kohn-Sham方程得到i由i构造rout(r)比较rin与 rout(r)计算总能EtotNoYesrin与rout混合初始结构精度控制NoYes输出结果： Etot、 i、 r(r)Vxc、Veff、En(k)、N(E)交换关联HF(ab initio)工作流程假设初始分子轨道fi(0)代入Fock算子求解得到fi(1)判断fi(1)-fi(0)非常接近：得解否则：令fi(1)=fi(0)SCF自洽场过程LCAO-MO获得F: 原子位置调整结构基函数/基组STF, GTF单z：STO-3G双z：6-31G三z：6-311G扩展

7、：*, +局域密度近似(LDA)假设能量泛函为：将能量密度泛函在某个定域 中展开：只取第一项，就是局域密度近似(LDA)，即交换关联能为电子密度的泛函，可写成定域积分的形式：其中 为单粒子多体交换关联能。LDA是指将非均匀电子体系近似为均匀电子体系，体系密度是恒定的。 表示均匀电子体系的密度，该密度等于原来非均匀电子体系r处的电子密度。如果考虑自旋(即LSDA)，则交换关联能可写为：这样一个简单的近似往往能给出很好的结果。这直接导致了后来DFT的广泛应用。尽管L(S)DA获得了巨大的成功，但是也有许多不足之处，如系统地高估结合能和HOMO能量、低估分子的键长(或键能)和晶体的晶格参数、处理富电

8、子体系不准确等。在L(S)DA基础上的成功改进有广义梯度近似(GGA)。广义梯度近似(GGA)GGA近似下，交换关联能是电子密度及其梯度的泛函。对于非常高的电子密度，交换能起主导作用，GGA的非局域性更适合处理密度的非均匀性，可以大大改进交换关联能计算结果。对较轻的元素的结果一般与实验符合得较好(优于LDA)GGA比LDA计算慢很多，但有许多相对优势：计算解离能和内聚能更接近实验值(LDA高估)对分子键长和固体格常数计算更准确(LDA低估)CB、MB、HB和VB等键能计算更准确(优于LDA)可以很好的描述分子体系的热力学性质(优于LDA)常见的泛函形式：LDAPWC：近期发展PWC泛函是在对V

9、WN泛函的某些错误校正后的结果，是DMol3模块的默认泛函。以DMol3模块为例，常用的两种LDA方法是VWN和PWCVWN：最常用的LSD (Local spin density)相关势函数，用来拟合电子气的精确数值结果。可精确预测共价体系结构、频率和相关能量，但往往高估键能、不能处理弱健体系，如氢键等。GGA的泛函形式GGA方法实际上是一种经验性的描述函数。对于不同的研究体系，计算所使用的泛函的精确度实际上是不一样的。如果要求计算的结果准确可靠，则需要对相关的函数进行查阅，或者从相关文献中查找所使用的泛函。以Perdew为首的一派认为“泛函必须以一定的物理规律为基础”，如标度关系、渐进行为

10、等。基于这种理念构造的PBE泛函是现在用的最广泛的GGA泛函之一。以Becke为首的一派认为“一切都是合法的”，所以人们可以选择任何可能的泛函形式，通常泛函的参数由拟合大量的计算数据得到。构造GGA交换相关泛函的方法分为两个流派：DMol3中的泛函PW91Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992).BOPTsuneda, T.; Suzumura, T.; Hirao, K. J. Chem. Phys., 110, 10664 (1999).BPBecke, A. D. J. Chem. Phys., 88, 2547 (19

11、88).Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992).VWN-BPVosko, S. J.; Wilk, L.; Nusair, M. Can. J. Phys., 58, 1200-1211 (1980). Becke, A. D. J. Chem. Phys., 88, 2547 (1988).Perdew, J. P.; Wang, Y. Phys. Rev. B, 45, 13244 (1992).PBEPerdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett., 77

12、, 3865 (1996).RPBEHammer, B.; Hansen, L. B.; Norskov, J. K. Phys. Rev. B, 59, 7413 (1999).BLYPBecke, A. D. J. Chem. Phys., 88, 2547 (1988). Lee, C.; Yang, W.Parr, R. G. Phys. Rev. B, 37, 786 (1988).HCTHBoese, A. D.; Handy, N. C. J. Chem. Phys., 114, 5497 (2001).常用的GGA泛函HCTH：对无机物和氢键体系的热力学计算结果较为可靠PBE(

13、Perdew-Burke-Enzerhof)：交换项与Becke相类似，相关项与Perdew-Wang函数相接近。具有较强的物理背景，主要用于固体计算，可靠的数值计算性能。 DFT中最常用；RPBE：热力学计算更可靠P91, BP, BLYP, BOP：一般是Becke交换函数(B88)与Perdew-Wang相关函数(BP)或者Lee-Yang-Parr相关函数(BLYP)组合使用。常用的为PW91和B3LYPVWN-BP：COSMO或者COSMO-RS计算推荐使用。缺点：计算反应能垒的数据和实验相比仍任会有所低估。关于DFT与泛函的补充知识请参阅文件密度泛函的发展与应用(中科大杨金龙?)，

14、里面介绍了：密度泛函的“雅各梯”密度泛函的选用新密度泛函的发展方法以及大量的参考文献课间故事2003 Nobel Prize for Chemistry诺贝尔化学奖授予两位医生？阿格雷(P. Agre, 1949-)，Johns-Hopkins大学医学院，医生、分子生物学家。麦金农(R. MacKinnon, 1956-), Rockefeller大学Howard-Hughesy医学研究中心教授。Molecular channels let us enter the chemistry of the cell.瑞典皇家科学院2003年10月日宣布，将该年度诺贝尔化学奖授予美国科学家彼得阿格雷和

15、罗德里克麦金农，分别表彰他们发现细胞膜水通道，以及对离子通道结构和机理研究作出的开创性贡献。力场方法(Force Field)HF, p-HF, SE：波函数 vs. 能量DFT：电子密度 vs. 能量FF：经典力学模型(如谐振子) vs. 能量基本的BO近似：原子核的运动与电子的运动可以看成是独立的分子是一组靠各种作用力(过近：互相排斥；过远：互相吸引。莫尔斯势？)维系在一起的原子集合，真实的分子结构是这些作用力达到平衡状态的表现。什么是FF？基本思想(D. Andrews, 1930)：分子内部的化学键都有“自然”的键长值和键角值。分子要调整它的几何构象，以使其键长值和键角值尽可能接近自然

16、值，同时也使非键作用(范德华力)处于最小的状态，给出原子核位置的最佳排布。函数(组)的参数由实验或第一原理计算得到的训练组拟合得到，可以较为精确地推广至预测类似的分子结构和能量等性质(可移植性: Transferability)一个或一组可以描述原子间相互作用简单解析函数(组) 简单代数函数取代了薛定谔方程！一个典型的FF应该包括：FF代数方程：势函数FF参数化及其可靠性训练组(Training Sets)：实验上的热力学和光谱学数据，或者精确的第一原理计算的数据量要充足！参数化(Parameterization)：势函数参数的的置定过程，通常称之为训练，最后获得的是一系列固定参数。例如：各类

17、键长、键角的“本征值”一般取自晶体学、电子衍射或其他的谱学数据，键伸缩和角变力常数主要由振动光谱数据确定，等等。不同FF的函数形式不同，参数不同。FF计算结果的准确性和可靠性主要取决于势函数和参数。典型的FF (I)传统力场AMBER：由Kollman课题组开发的力场，是目前使用比较广泛的一种力场，适合处理生物大分子。CHARMM：由Karplus课题组开发，对小分子体系到溶剂化的大分子体系都有很好的拟合。CVFF：可以用于计算无机体系MMX：包括MM2和MM3，是目前应用最为广泛的一种力场，主要针对有机小分子ReaxFF：反应力场，可以预测键的生成和断裂，具有“预测”反应的特性第二代力场势能

18、函数形式比传统力场要更加复杂，涉及的力场参数更多，计算量也更大，相应地更加准确。CFF：力场家族，包括了CFF91、PCFF、CFF95等很多力场，可以进行从有机小分子、生物大分子到分子筛等诸多体系的计算。COMPASS：MSI公司，擅长进行高分子体系计算。MMFF94：是目前最准确的力场之一通用力场ESFF：MSI公司，可以进行有机、无机分子的计算UFF：以计算周期表上所有元素的参数Dreiding：适用于有机小分子、大分子、主族元素代表性程序：LAMMPS, GROMACS, DISCOVER等FF的应用范围和优劣主要应用于有机领域(现在逐渐包括无机)尤其是在大分子体系例如生物过程、药物设

19、计等领域应用非常广泛优点：最快的计算方法，可处理的体系最大(1M)，在适当范围内给出的计算精度与第一原理的精度相当缺点：对能量的处理比较粗糙，可拓展性、可移植性需要仔细鉴定典型的FF应用NMR or X-ray structure refinementProtein structure predictionProtein folding kinetics and mechanicsConformational dynamicsGlobal optimizationDNA/RNA simulationsMembrane proteins/lipid layers simulationsCompu

20、ter aided drug design分子动力学(Molecular Dynamics)Kinetics vs. Dynamics 速度 vs. 演化单点能(Single-Point Energy)势能面上的任意单点几何优化(Geometry Optimization)走向势能面局部最小点过渡态(Transition State)势能面的鞍点分子动力学？MD一种计算经典多体体系的平衡和性质传递的方法，组成粒子根据已知的微观粒子间的相互作用，遵循牛顿经典力学(Classic MD)或量子力学规律(ab initio MD)进行运动，从而模拟出分子体系动态演化行为。通过求解出多体体系的运动方程

21、，得到全体粒子在相空间中的轨道，然后统计得到系统的性质，通过体系的微观性质来求算其宏观性质，属于微观尺度的模拟技术。MD包含温度(T)、压力(P)和时间(t)，可以模拟体系在实际条件约束下随时间推进的微观动态过程，进而通过统计力学方法得到体系的平衡参数、力学性质、热力学性质等。典型的分子动力学程序有：LAMMPS (C); VASP (AI)MD的运行流程(以NVE为例)给定每个分子的初始位置ri(0)和速度vi(0)计算每个分子的受力Fi和加速度ai求出每个分子运动一个时间步长后的新位置和速度统计系统的热力学性质及其它物理量统计性质不变？打印结果，结束结构移动所有分子到新位置并赋予当前时刻的

22、速度YESNO分子动力学的运行设置Structure?NVE or NVT?Temperature and ThermostatTime step?Total simulation time/Number of steps?结构：周期边界条件为什么要采用周期边界条件(Periodic Boundary Conditions)？通常选取小数体系(几十个到数千个分子)作为研究对象，但表面和内部的分子受力差别较大，将会产生表面效应。受限于实际计算模拟能力(计算机硬件限制)，不可能模拟一个真实体系(10nm60k)，因此，必须利用有限空间中有限个分子来代表真实体系。PBC：为减小表面效应，同时构造一个

23、准无穷大体积，以更能精确地代表宏观系统。PBC基本原理令基本元胞(最小影像)完全等同重复无穷多次， 当有一个粒子穿过基本MD元胞壁时， 就让这个粒子以相同的速度穿过此表面对面的壁重新进入该MD元胞内。利用PBC，可以从模拟有限微观系统来预测宏观系统的演化行为和性质。相空间(Phase Space)经典系综理论为描述系统微观状态而引入的一个基本概念(Gibbs, 1902)，表示某系统所有可能状态的空间，是一个六维假想空间，其中动量和空间各占三维。相空间中的每一点，代表系统一种可能的运动状态，系统的不同运动状态，可以用相空间中不同的点来代表。系统的状态随时间变化，相空间中代表运动状态的点将在相空

24、间中移动，形成“轨迹”。理论上，时间足够时可以遍历相空间的所有态。统计系综(Ensemble)(统计)系综指处在一定宏观条件下，大量性质和结构完全相同的系统的集合，是采用统计方法描述系统的统计规律性时引入的一个最基本的概念。相空间中的不同点代表了系统的不同运动状态，也可以把它们看作是性质和结构相同而初始条件不同的系统所处的状态，这些系统的集合就是系综。统计系综的各个系统反映了实际系统在不同时刻的面貌。常用系综：微正则系综、正则系综、巨正则系综、等温等压系综。系综与分子动力学微正则系综(Microcanonical Ensemble, NVE)：孤立、保守的系统的统计系综，系统与外界不交换能量，

25、粒子数守恒，体积也不发生变化。例如：孤立体系正则系综(Canonical Ensemble, NVT)：系统的粒子数(N)、体积(V)和温度(T)都保持不变且总动量守恒的系综。例如：恒温热浴等温等压系综(IsothermalIsobaric Ensemble, NPT)：系统的粒子数(N)、压强(P)和温度(T)恒定，系综内的每个体系可和其他体系进行能量(E)和体积(V)交换，但系综内各体系的能量总和以及体积总和固定。例如：相变过程温度 (Temperature)温度与初速度：初始速度矢量根据伪随机数进行设置，使体系的总动能与目标温度对应。利用能量均分定理，根据Maxwell-Boltzmann分布来任意设置初始速度：弛豫(Relax)：初始坐标与速度设置后，在模拟体系性质之前，必须先使得体系进行趋于平衡的过程。在弛豫过程中，动能与势能相互转化，当动能、势能与总能量只在平均值附近涨落时，体系就达到平衡(NVE)质量mi的原子i在温度T下沿x方向速度为vix的概率温控子(Thermostat, NVT)NVT：恒温热浴，模拟过程中必须保证体系温度恒定(T恒定，总能可