理论课讲义-04.p分子模拟与设计

分子模拟与设计Given

by:周成冈

Ph.

D.,Professor

of

Chemistry中国地质大学( )材料与化学学院Faculty

of

Materials

Science

and

Chemistry2013

Spring基于波函数的计算基于波函数的计算方法Ab

initio

HF后HF方法Semiempirical原理基于量子力学(求)基于量子力学(求)基于量子力学(求)相关能不考虑考虑不考虑算法严格从头算，无任何经验参数严格从头算，无任何经验参数大量经验参数和近似速度非常慢最慢非常快精度较为精确最精确定性基本准确可靠性非常可靠最可靠误差随意性大计算性质基态基态、激发态基态、激发态(ZINDO/S)处理体系～100～40>1000密度泛函理论从头算(HF,

p-HF,

SE)：基于复杂的多电子波函数，存在

3N个变量(N为电子数，每个电子包含3个空间变量)，体系能量为波函数的本征值。密度泛函理论(DensityFunctionalTheory)：电子密度取代波函数(电子密度为3个空间变量的函数)，体系的基态能量为电子密度的泛函(Hohenberg-Kohn第一定理)。DFT通过Kohn-Sham框架实现：复杂的多体问题被简化成一个没有相互作用的电子在有效势场中运动的问题。有效势包括所有的相互作用效应，即Hartree势、交换势

(Pauli原理决定的相互作用势)和关联势(某电子对整个电荷分布的影响所产生的势)。DFT典型程序DFT的地位与作用近十几年来，DFT同分子动力学方法相结合，有大量的新发展( -Principles

MD)在材料设计、、模拟计算和评价诸多方面得到了广泛应用，是目前最流行的计算方法(量化计算界至少70%以上的人采用)已成为计算凝聚态物理、计算材料科学和计算量子化学的重要基础和 技术逐渐进入工业技术领域DFT的优点DFT

vs.

HF:

降低了计算体系的维数[波函数→电子密度]提供了第一性原理或从头算的计算框架，在计算精度和计算量之间达到了较好的权衡在这个框架下可以发展各式各样的结构和能带计算方法，例如材料电子结构和几何结构、固体和液态金属中的相变等。这些方法已发展成为用量子力学方法计算力的精确的分子动力学方法。DFT的应用范围电子基态能量与原子位置之间的关系可以用来确定分子或晶体的结构，从而得到各种性质当原子不处在它的平衡位置时，给出作用在原子上的力因此，DFT可以解决物理和化学中的大多数问题，如力、光、电、磁、热等电离势、电子亲和能

振动光谱化学反应生物分子的结构催化活性位置的特性…Hohenberg-Kohn定理相互作用多体系统的粒子数密度是决定系统物理性能的基本变量，即系统所有物理性质都是由粒子数密度决定的，包括总能量和系统态密度。Hohenberg-Kohn能量泛函形式：

VextE

,V

T

Vee

T

Vee

drV

r

r

F

drV

r

r

取严格的基态密度时能量泛函取极小值，并代表电子系统的基态密度。能量是密度的唯一泛函，仅当

H-K能量泛函与变分方程E

,V

T

T0

VH

Exc

d其中T0

iii,

j

,V2m

i

H2iT0:无相互作用系统的动能，T是相互作用系统的动能VH:Hartree势EXC:交换关联能，包含了Hartree近似以外所有的多体效应Exc

T

T0

Vee

VH

变分方程：

E

,V

dr

r

0限制性条件是总粒子数守恒：N

dr

r

Kohn-Sham自洽方程组

2212effeffi2mdrVxc

r

e

r

r

r

2

V

r

r

r

i

i

iV

r

V

r

Vc

r

Vxc

r

,Vc

r

r

r

Exc,

r

i将相互作用多体系统基态问题在形式上严格转化为在有效势场中运动的独立粒子的基态问题。处理KS方程中的交换关联项[EXC,VXC]的处理需要用到局域密度近似(LocalDensity

Approximation)或广义梯度近似(Generalized

Gradient

Approximation)KS方程求解流程图求解φ、Vxc、Veff求解Kohn-Sham方程得到ψi由ψi构造out(r)比较in与out(r)计算总能EtotNoYesin与out混合in(r)

(r)=Σat(r)

初始结构精度控制NoYes输出结果：Etot、ψi、(r)Vxc、Veff、En(k)、N(E)交换关联HF(ab

initio)工作流程i假设初始分子轨道

(0)代入Fock算子i求解得到

(1)判断i

-i(1)

(0)非常接近：得解(1)否则：令i

=i(0)SCF自洽场过程LCAO-MO获得F:

原子位置调整结构基函数/基组STF,

GTF单：STO-3G双：6-31G三：6-311G扩展：**,++局域密度近似(LDA)假设能量泛函为：G

drgr

,其中gr

是r处能量密度泛函将能量密度泛函在某个定域

r

中展开：gr

g0

r

1

只取第一项，就是局域密度近似(LDA)，即交换关联能为

电子密度的泛函，可写成定域积分的形式：

xcLDAxcEdr

r

r

其中xc

为单粒子多体交换关联能。LDA是指将非均匀电子体系近似为均匀电子体系，体系密度是恒定的。

r

表示均匀电子体系的密度，该密度等于原来非均匀电子体系r处的电子密度。如果考虑自旋(即LSDA)，则交换关联能可写为：xcLSDAxcE

,

dr

r

,

r

r

这样一个简单的近似往往能给出很好的结果。这直接导致了后来DFT的广泛应用。尽管L(S)DA获得了巨大的成功，但是也有许多不足之处，如系统地高估结合能和HOMO能量、低估分子的键长(或键能)和晶体的晶格参数、处理富电子体系确等。在L(S)DA基础上的成功改进有广义梯度近似(GGA)。广义梯度近似(GGA)GGA近似下，交换关联能是电子密度及其梯度的泛函。对于非常高的电子密度，交换能起主导作用，GGA的非局域性更适合处理密度的非均匀性，可以大大改进交换关联能计算结果。对分子键长和固体格常数计算更准确(LDA低估)计算解离能和内聚能更接近实验值(LDA高估)对较轻的元素的结果一般与实验符合得较好(优于LDA)CB、MB、HB和VB等键能计算更准确(优于LDA)可以很好的描述分子体系的热力学性质(优于LDA)GGA比LDA计算慢很多，但有许多相对优势：常见的泛函形式：LDAPWC：近期发展PWC泛函是在对VWN泛函的某些错误校正后的结果，是DMol3模块的默认泛函。以DMol3模块为例，常用的两种LDA方法是VWN和PWCVWN：最常用的LSD

(Local

spindensity)相关势函数，用来拟合电子气的精确数值结果。可精确

共价体系结构、频率和相关能量，但往往高估键能、不能处理弱健体系，如氢键等。GGA的泛函形式GGA方法实际上是一种经验性的描述函数。对于不同的研究体系，计算所使用的泛函的精确度实际上是不一样的。如果要求计算的结果准确可靠，则需要对相关的函数进行查阅，或者从相关文献中查找所使用的泛函。构造GGA交换相关泛函的方法分为两个流派：以Becke为首的一派认为“一切都是合法的”，所以人们可以选择任何可能的泛函形式，通常泛函的参数由拟合大量的计算数据得到。以Perdew为首的一派认为“泛函必须以一定的物理规律为基础”，如标度关系、渐进行为等。基于这种理念构造的

PBE泛函是现在用的最广泛的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(1988).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,

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.; ,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(Perdew-Burke-Enzerhof)：交换项与Becke相类似，相关项与Perdew-Wang函数相接近。具有较强的物理背景，主要用于固体计算，可靠的数值计算性能。[DFT中最常用；RPBE：热力学计算更可靠]P91,BP,

BLYP,

BOP：一般是Becke交换函数(B88)与Perdew-Wang相关函数(BP)或者Lee--Parr相关函数

(BLYP)组合使用。[常用的为PW91和B3LYP]VWN-BP：COSMO或者COSMO-RS计算

使用。缺点：计算反应能垒的数据和实验相比仍任会有所低估。[关于DFT与泛函的补充知识]请参阅文件《密度泛函的发展与应用》(中杨金龙?)，里面介绍了：密度泛函的“雅各梯”密度泛函的选用新密度泛函的发展方法以及大量的参考文献课间故事—2003

Nobel

Prize

for

Chemistry化学奖授予两位医生？？阿格雷(P.Agre,1949-)，Johns-Hopkins

大学医学院，医生、分子生物学家。麦金农(R.MacKinnon,1956-),Rockefeller

大学Howard-Hughesy医学教授。Molecular

channels

let

us

enterthe

chemistry

of

the

cell.瑞典宣布，将该年度诺贝尔化学奖授予科学家彼得·阿格雷和罗德里克·麦金农，分别表彰他们发现细胞膜水通道，以及对离子通道结构和机理研究作出的开创性贡献。20世纪50年代中期，科学家发现，细胞膜中存在着某种通道，只允许水分子出入，人们称之为水通道

(Aquaporins)。道蛋的特2000白的殊结1998通道离子道上过滤

液的情在水滤嘴，而一如接上溶液有个蛋白质结化学真的是一门中心学科！深入到原子层次的研究工作只能通过模拟来进行。这里的机理解析工作就是由力场方法进行的。力场方法(Force

Field)HF,

p-HF,

SE：波函数vs.

能量DFT：电子密度vs.

能量FF：经典力学模型(如谐振子)vs.

能量基本的BO近似：原子核的运动与电子的运动可以看成是独立的分子是一组靠各种作用力(过近：互相排斥；过远：互相吸引。

势？)维系在一起的原子集合，真实的分子结构是这些作用力达到平衡状态的表现。FF？基本思想(D.Andrews,1930)：分子的化学键都有“自然”的键长值和键角值。分子要调整它的几何构象，以使其键长值和键角值尽可能接近自然值，同时也使非键作用(范德华力)处于最小的状态，给出原子核位置的最佳排布。一个典型的FF应该包括：一个或一组可以描述原子间相互作用简单解析函数(组)

[简单代数函数取代了薛定谔方程！！]函数(组)的参数由实验或第一原理计算得到的训练组拟合得到，可以较为精确地推广至类似的分子结构和能量等性质(可移植性:Transferability)FF代数方程：势函数FF参数化及其可靠性训练组(Training

Sets)：实验上的热力学和光谱学数据，或者精确的第一原理计算的数据[量要充足！]例如：各类键长、键角的“本征值”一般取自晶体学、电子衍射或其他的谱学数据，键伸缩和角变力常数主要由振动光谱数据确定，等等。参数化(Parameterization)：势函数参数的的置定过程，通常称之为训练，最后获得的是一系列固定参数。不同FF的函数形式不同，参数不同。FF计算结果的准确性和可靠性主要取决于势函数和参数。典型的FF

(I)传统力场AMBER：由Kollman课题组开发的力场，是目前使用比较广泛的一种力场，适合处理生物大分子。CHARMM：由Karplus课题组开发，对小分子体系到溶剂化的大分子体系都有很好的拟合。CVFF：可以用于计算无机体系MMX：包括MM2和MM3，是目前应用最为广泛的一种力场，主要针对有机小分子–ReaxFF：反应力场，可以键的生成和断裂，具有“”反应的特性第二代力场–势能函数形式比传统力场要更加复杂，涉及的力场参数 ，计算量也更大，相应地更加准确。PCFF、CFF95等很生物大分子到分子–CFF：力场多力场，可筛等诸多体COMPASSMMFF94：高分子体系计算。之一通用力场ESFF：MSUFF：以计无机分子的计算参数–Dreiding：适用于有机小分子、大分子、主族元素代表性程序：LAMMPS,

GROMACS,

DISCOVER等FF的应用范围和优劣主要应用于有机领域(现在逐渐包括无机)尤其是在大分子体系例如生物过程、药物设计等领域应用非常广泛优点：最快的计算方法，可处理的体系最大(>1M)，在适当范围内给出的计算精度与第一原理的精度相当缺点：对能量的处理比较粗糙，可拓展性、可移植性需要仔细鉴定[典型的FF应用]NMR

or

X-raystructure

refinementProtein

structure

predictionProteinfolding

kineticsandmechanicsConformationaldynamicsGlobal

optimizationDNA/RNA

simulationsMembrane

proteins/lipid

layerssimulationsComputer

aideddrug

design分子动力学(Molecular

Dynamics)Kinetics

vs.

Dynamics[速度vs.

演化]单点能(Single-PointEnergy)势能面上的任意单点几何优化(Geometry

Optimization)势能面局部最小点过渡态(Transition

State)势能面的鞍点分子动力学？MD一种计算经典多体体系的平衡和性质传递的方法，组成粒子根据已知的微观粒子间的相互作用，遵循牛顿经典力学

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

or

NVT?Temperature

andThermostatTime

step?Total

simulationti mber

of

steps?结构：周期边界条件为什么要采用周期边界条件(PeriodicBoundaryConditions)？通常选取小数体系(几十个到数千个分子)作为研究对象，但表面和的分子受力差别较大，将会产生表面效应。受限于实际计算模拟能力(计算机硬件限制)，不可能模拟一个真实体系(10nm>60k)，因此，必须利用有限空间中有限个分子来代表真实体系。PBC

：为减小表面效应，同时构造一个准无穷大体积，以更能精确地代表宏观系统。PBC基本原理令基本元胞(最小影像)完全等同重复无穷多次，当有一个粒子穿过基本MD元胞壁时，就让这个粒子以相同的速度穿过此表面对面的壁重新进入该MD元胞内。利用PBC，可以从模拟有限微观系统来宏观系统的演化行为和性质。相空间(Phase

Space)经典系综理论为描述系统微观状态而引入的一个基本概念

(Gibbs,1902)，表示某系统所有可能状态的空间，是一个六维假想空间，其中动量和空间各占三维。相空间中的每一点，代表系种可能的运动状态，系统的不同运动状态，可以用相空间中不同的点来代表。系统的状态随时间变化，相空间中代表运动状态的点将在相空间中移动，形成“轨迹”。理论上，时间足够时可以遍历相空间的所有态。统计系综(Ensemble)(统计)系综指处在一定宏观条件下，大量性质和结构完全相同的系统的集合，是采用统计方法描述系统的统计规律性时引入的一个最基本的概念。相空间中的不同点代表了系统的不同运动状态，也可以把它们看作是性质和结构相同而初始条件不同的系统所处的状态，这些系统的集合就是系综。统计系综的各个系统反映了实际系统在不同时刻的面貌。常用系综：微正则系综、正则系综、巨正则系综、等温等压系综。系综与分子动力学微正则系综(MicrocanonicalEnsemble,NVE)：孤立、保守的系统的统计系综，系统与外界不交换能量，粒子数守恒，体积也不发生变化。[例如：孤系]正则系综(Canonical

Ensemble,

NVT)：系统的粒子数(N)、体积(V)和温度(T)都保持不变且总动量守恒的系综。[例如：恒温热浴]othermal–Isobaric

Ensemble,

NPT)：系压强(P)和温度(T)恒定，系综内的每进行能量(E)和体积(V)交换，但系综内各体及体积总和固定。[例如：相变过程]温度(Temperature)温度与初速度：初始速度矢量根据伪随机数进行设置，使体系的总动能与目标温度对应。利用能量均分定理，根据

Maxwell-Boltzmann分布来任意设置初始速度：m

v21/

2mP(vix

)

i

2

kBT

exp

i ix

2kBT

质量mi的原子i在温度T下沿x方向速度为vix的概率弛豫(Relax)：初始坐标与速度设置后，在模拟体系性质之前，必须先使得体系进行趋于平衡的过程。在弛豫过程中，动能与势能相互转化，当动能、势能与总能量只在平均值附近涨落时，体系就达到平衡(NVE)温控子(Thermostat,