从计算机在化学中的应用一第七章分子模拟基础第七章分子模拟基础7.1引言分子模拟(MolecularModeling)这个概念近年来才开始使用的，这个研究领域在最近十儿年里随计算机的发展而迅速发展，并且成为化学领域中与实验仪器一样成为日常的研究工具：通过这个研究工具可以从更深的层次理解化学反应及过程，确定或预测分子结构以及化学体系的稳定性，估算不同状态之间的能量，在原子水平上解释反应途径和机理，了解分子结构与性质之间的关系，预测设计有价值的分子等。在这个研究领域中也常使用理论化学（TheoreticalChemistry）和计算化学（ComputationalChemistry)的概念。一般来说，理论化学等价于量子力学：计算化学则除了量子力学外还包括分子力学、最小化、模拟、构象分析和其他以计算机为基础的用以理解和预测分子体系性质的方法。分子模拟泛指用于模拟分子或分子体系性质的方法，定位于表述和处理基手三维结构的分子结构和性质。近年来发展了许多分子模拟的计算机软件，其一般都包括以下几个方面：1）建立和显示分子；2）优化分子结构：3）研究分子的反应性（如从轨道能量、组合系数、节点性质等根据前线轨道理论研究分子不同取代基的相对反应性、反应选择性以及亲核、亲电试剂的反应位置等)；4）计算产生和显示分子轨道、电荷分布、静电势图等；5）估计化学反应途径和机理；6）研究分子的动态性质：7）研究分子的结构-性质(活性)关系(QSPR/QSAR)。下面我们首先介绍几个分子模拟中常用的概念。7.1.1坐标系对于分子模拟程序来说，无疑确定原子和分子在体系中的位置是非常重要的。一般采用两种方法来表示原子/分子的位置，一种是笛卡儿坐标，另一种是内坐标。笛卡儿坐标比较好理解，即对每一个原子都用一套XYZ来表示其在体系中的位置，但笛卡儿坐标常常并不是十分方便的，因为它不能直接显示我们所研究的分子结构的常用数据如键长、键角等，因此在研究中也常使用内坐标。内坐标是采用一个原子对体系内其它原子相对位置进行定位的方法，常被写成乙矩阵的形式。每一个原子对应于Z矩阵的一行，我们以完全交叉构象的乙烷为例，在Z矩阵第一行，定义原子1为C原子；原子2也为C原子，与原子1的距离为1.54A；原子3为H原子，与原子1之间的键长为1A，3-1-2的键角为109.5°；原子4为H原子，与原子1之间的键长为1A，4-1-2的键角为109.5，4-1-2-3的二面角为120°：后面的原子与第4个原子相同，用与前面原子之间的键长、键角和二面角定位量子力学计算程序一般使用内坐标，分子力学程序一般使用笛卡儿坐标，现在的程序一般都可以使用两种中的任意一种。7 - 1
7 1 7.1 (Molecular Modeling) (Theoretical Chemistry) (Computational Chemistry) 1) 2) 3) ( ) 4) 5) 6) 7) - ( ) (QSPR/QSAR) 7.1.1 / XYZ Z Z Z 1 C 2 C 1 1.54Å 3 H 1 1Å 3-1-2 109.5° 4 H 1 1Å 4-1-2 109.5° 4-1-2-3 120° 4
双计算机在化学中的应用月一第七章分子模拟基础1c2C1.5483H1.0109.5214H1.0109.52120315H1.0231109.5-1206H1.02109.51803LH1.0109.512068H1.0-12062109.5167.1.2势能面分子模拟中默认采用Born-Oppenheimer近似，因此分子基电子态的能量可以看作只是核坐标的函数，核的移动必然引起能量的变化。当然不同的移动方式引起的能量变化大小不同，如乙烷C-C键长从平衡态移动0.1A需要3kcal/mol的能量：而两个Ar之间非共价移动1A只需要0.1kcal/mol。体系能量的变化可以看成是在一个多维面上的运动，这个面称之为势能面(energysurface)。对势能面上的点我们最感兴趣的是不动点(stationarypoint)，即对内坐标或笛卡儿坐标的一阶微商为零的点，在不动点，所有原子所受到的力为零。极小点是一种不动点，它代表稳定的结构：另一种不动点是鞍点，代表过渡态的结构。通过数学方法可以在势能面上找出这些不动点。100右图为两个变量的势能面图。A和C点为极小点，这意味着在结0构A和C为稳定结构，如果C比ableA稳定，C为该势能面的最小点N在势能面上进行的构型优化可以0找到象A和C这样的极小点。B是从A到最小能量途径上的极大点，称作鞍点(SaddlePoint)，在B点所代表的结构中所有原子所受的力仍然为零，B点代表从A到C40802060100变化的过渡态(TransitionState)。连variable 1(degrees)接反应物和产物之间的最低能量途径称为反应坐标。如正工烷沿中心两个C原子间的C-C键旋转的反应途径有一系列的不动点，A、C、E和G为极小点，B、D和F为极大点。除少数简单体系外，大多数势能面是非常复杂的，如肌红蛋白(约150residues)7-2
7 2 1 C 2 C 1.54 3 H 1.0 1 109.5 2 4 H 1.0 1 109.5 2 120 3 5 H 1.0 1 109.5 2 -120 3 6 H 1.0 2 109.5 1 180 3 7 H 1.0 2 109.5 1 120 6 8 H 1.0 2 109.5 1 -120 6 8 6 2 7 5 1 3 4 7.1.2 Born-Oppenheimer C-C 0.1Å 3kcal/mol Ar 1Å 0.1kcal/mol (energy surface) (stationary point) A C A C C A C A C B A C (Saddle Point) B B A C (Transition State) C C-C A C E G B D F ( 150 residues) variable 1(degrees) variable 2(degrees) A B C
从计算机在化学中的应用一第七章分子模拟基础的构象有4x1034种之多。因此，单一的稳定结构并不能完全反映一个分子的性质，必须考虑热力学Boltzman分布：除小分子外，要找到极小点有时是不可能的。RSotu.eoxeuCH:o0-180-120-60120180Angel(°)7.1.3分子图形我们在第三章已经介绍了分子图形的一些绘制方法和软件，分子图形的优点是可以方便地得到一些分子结构的信息如键长、键角等。分子显示在计算机显示器屏幕上可以采用stick，tube,ball&stick和space-filling等类型的图形(如下图的aspirin分子图形)。大分子体系中要显示所有的原子是困难的事，因此除使用上面的类型外，还可使用cartoon、ribbon等类型的显示方法。下图为松弛肽的cartoon和ribbon图7-3
— 7 3 4 1034 Boltzman -180 -120 -60 0 60 120 180 0 3 6 9 • • • • • • • G F E D C B A H H CH3 H H CH3 H H CH3 H3C H H H H CH3 H3C H H H H H3C H H CH3 H H CH3 H CH3 H H H CH3 H CH3 H H CH3 H H CH3 Energy(kcal·mol-1 ) Angel(°) 7.1.3 stick, tube, ball&stick space-filling ( aspirin ) cartoon ribbon cartoon ribbon
从计算机在化学中的应用一第七章分子模拟基础7.1.4表面分子模拟中常需要处理两个或更多的分子之间的非共价相互作用问题，这些问题涉及到范德华表面(vanderWaalssurface)和分子表面(molecularsurface)。范德华表面是原子范德华球的简单空间重叠，可基于CPK或spacefilling模型。如果使用一个小的探针分子接近大分子的范德华表面，由于探针分子球的有限体积使一些缝隙成为死空间（探针球在大分子表面滚动时无法进入的缝隙）。分子表面是由探针球在分子范德华表面上滚动时的内向部分。它包括两类表面元素，接触表面(contactsurface)对应于指探针真正接触到范德华表面的部分，凹表面(re-entrantsurface)对应于那些由于缝隙太小而探针无法进入的部分。分子表面常用水作为探针，球半径为1.4A。可接近表面(accessiblesurface)的概念也被广泛使用，它的含义是探针分子在范德华表面上滚动时其球心所形成的面。分子表面的计算方法是由Connolly等人提出的，显示表面的方法很多，下面是一些显示类型。7 - 4
7 4 7.1.4 (van der Waals surface) (molecular surface) CPK space filling ( ) (contact surface) (re-entrant surface) 1.4Å (accessible surface) Connolly
双计算机在化学中的应用一第七章分子模拟基础7.2量子力学模型7.2.1Schrodinger方程化学反应必然要涉及到电子，而对于电子性质的计算，最常用的就是量子力学方法。讨论量子力学问题的起点自然从Schrodinger方程开始，单粒子定态Schrodinger方程：HY=EY(7.1)其中 E 为粒子能量，Halmiton 算符H=v2 +P2ma2a2拉普拉斯算符：V2=4ax?y+?求解Schrodinger方程为偏微分本征方程的问题。为本征函数，E为本征值Schrodinger方程只能对少数几个体系可以有精确解，如势箱中的粒子、谐振子、环中粒子、球中粒子和氢原子等。解这些问题的共同之处是对方程的可能解引入边界条件，因此对于在势箱中的粒子，波函数在边界处应为0，环中粒子的波函数必须有2元的周期性。另外一个在求解Schrodinger方程时需要的条件是波函数必须归一化，满足此条件的波函数被称作是归一化的。通常Schrodinger方程的解要满足正交化(orthogonal)，同时满足正交和归一条件的波函数称为正交归一的(orthonaomal)7.2.2氢原子和类氢离子在只包含一个电子的原子中，原子单位下Hamilton算符具有如下形式H=--v?_Z(7.2)2r因为原子是球对称的，可将Schrodinger方程转为极坐标系求解，r为到核(原点)的距离，为与z轴的夹角，Φ为在xy平面上的投影与x轴夹角。解可以写成只与r有关的径向函数R(r)部分和与e有关的角度函数Y(O，)部分的积：(7.3)Ynlm=Rn(r)Yim(0,Φ)波函数一般称作轨道(orbital)，其由3个量子数n,l,m决定，量子数可采用的值如下：主量子数n：0,1,2..：角量子数l:0,1...,(n-1)：磁量子数m：-1(1-1)..0.(1-1)l径向函数：(n-1-1)!(7.4)p'L(p)Ru(r)exp2n[(n + 1)]]3nnap=2Zrlnao，ao为Bohr半径。7- 5
7 5 7.2 7.2.1 Schrödinger Schrödinger Schrödinger H E (7.1) E Halmiton V m H 2 2 2 2 2 2 2 2 2 2 x y z Schrödinger E Schrödinger 0 2 Schrödinger Schrödinger (orthogonal) (orthonaomal) 7.2.2 Hamilton r Z H 2 2 1 (7.2) Schrödinger r ( ) z xy x r R(r) Y( , ) nlm=Rnl(r)Ylm( , ) (7.3) (orbital) 3 n, l, m n 0, 1, 2 l 0, 1, , (n-1) m -l,-(l-1), 0 (l-1),l ( ) 2 exp 2 [( )!] 2 ( 1)! ( ) 2 1 1 2 1 3 3 0 l n l nl L n n l n l na Z R r (7.4) =2Zr/na0 a0 Bohr