一、量子化学的简介在20世纪初,物理学家发现经典力学不能准确地解释小粒子之间的相互作用,例如分子和原子中的原子核、电子,但是这些粒子的相互作用却可以用量子力学描述。1926年薛定谔提出了关于描述微观粒子运动的“波动方程”,即Schrodinger方程,标志着现代量子时代的到来。时隔一年,1927年海特勒和伦敦用量子力学基本原理讨论了关于氢分子的结构问题,说明了两个氢原子能够结合成一个稳定的氢分子的原因,并且利用相当近似的计算方法,算出其结合能。由此,使人们认识到可以用量子力学原理讨论分子结构问题,从而逐渐形成了量子化学这一分支学科。量子化学是量子力学应用于化学形成的学科。从量子力学的观点看,化学变化是由原子核和电子组成的体系的各种稳定态或亚稳态之间的相互转换。量子化学就是通过求解体系的“波动方程”得到电子及核的运动状态,用以阐明各种谱图,反应规律,分子稳定性和反应活性等化学现象的一门学科。二、QTAIM 简介由于电子密度可以解释化学实验中观察到的现象,所以分子中原子量子理论被大家认识。20世纪60年代,Bader主要探究电子密度分布,并首次提出化学中最重要的性质即电子密度,他的想法为之后的“分子中原子”(AIM,Atoms in Molecules)理论提供了基础。通过研究的不断发展与完善,在其精密的量子力学理论基础上该理论被发展为“分子中原子”量子理论(QTAIM,Quantum Theory of Atoms in Molecules)。Bader早期对分子中电子密度分布的探讨与现代密度泛函理论(DFT)所得到的公式是一致的,现代DFT理论的计算能够进一步提高精确度,质量空前的电子密度图能够被快速描绘出。由于DFT的出现、计算机功率的惊人增长、以及算法的进步,导致了应用分子中原子量子理论的研究数量爆炸性增长,这些研究在固态物理学、材料科学、表面科学、X射线分析等一系列的科学领域中得到了非常广泛的应用。电子密度ρ(r)分布是QTAIM理论的核心,一个分子的电子密度分布的拓扑性质是由其电子密度梯度矢量场即一阶导数▽ρ(r)和电子密度的二阶导数Laplacian值▽2ρ(rb)决定的。若空间中某一点电子密度分布的一阶导数▽ρ(r)=0,该点被称为临界点,在该点建立Hessian(3*3)矩阵,可获得三个本征值(λ1,λ2,λ3),其中Laplacian值▽2ρ(rb)被定义为▽2ρ(rb)=λ1+λ2+λ3。如果 ▽2ρ>(rb)<0,表示在该区域 Laplacian值为负值,势能占优势,负电荷在该处集中,并且值越小,化学键的共价性越强;如果▽2ρ(rb)>0,表示该区域Laplacian值为正值,动能占优势,负电荷在该处分散,值越大,化学键的离子性越强;若▽2ρ(rb)接近于0时,一般形成的是较弱的化学键。而且,电子密度ρ(r)的海森矩阵的对角化给出有序特征值λ1<λ2<λ3的集合和相应的特征向量e1,e2,e3。临界点的性质由相应的海森矩阵的特征值(λ1、λ2和λ3)的有序集合揭示,临界点标记为(R,ω),其中R代表矩阵的秩,符号ω是三个特征值正负号数目的代数和。在临界点处ω只能等于3。临界点它有四种类型,核临界点(3,-3),一般来说它指的是对应核位置的局部最大值;当(R,ω)对应的是(3,-1)和(3,1),它们分别是键临界点(BCP)和环临界点(RCP);(3,+3)称为笼临界点(CCP)。特征向量e3表示在BCP处键路径的方向。电子积累的最佳方向和最不优先的方向分别是e2和e1。从数学的角度对真正空间下分子的本质存在形式进行捕捉得到相应的分子图,图中的相应关系表示的是临界点和相关键合路径的结合。当一个分子是独立稳定时,它会满足Poincare-Hopf关系式:n-b+r-c=1(其中n,b,r,c 分别代表NCPs,BCPs,RCPs,CCPs的数目)。结合Poincare-Hopf关系式从拓扑学角度我们课题组对分子的维度进行了重新定义:把包含CCP即{n,b r,c}的分子定义为3-DQT,含有RCP但不含CCP即{n,b,r}的分子定义为2-DQT,含有NCP和BCP即{n,b}的分子定义为1-DQT,只含有NCP即{n}孤立的单原子结构,定义为0-DQT。这种维度的划分更显细致,它和传统的欧几里得维度方法的区别在于,此种方法有效地区分了 2-D和1-D的维度。本论文采用QTAIM以及应力张量理论,使用四种不同的计算方法对酪氨酸-甘氨酸二肽构象异构体的分子拓扑结构,键性质和键临界点进行分析,探讨这四种计算方法的不同点。接着讨论结构效应和键性质导致的BSSE的差异。介绍一个新的概念螺旋度H,在这四种不同的计算方法中,使用MP2计算方法时发现其对扭转角度有不同于其他三种计算方法的独特响应。其次为了进一步发展理论的适用范围,把QTAIM以及应力张量理论作为一种新的指导思想,研究采用键路径集合组来分析六种不同的对位取代联苯;最后在应力张量理论中拓展新的基于三维矢量的化学键解释来分析整个键径,并且对QTAIM和应力张量理论下的键路径框架集进行对比,探讨二者之间的差异性。三、基于QTAIM和应力张量理论研究二肽构象在不同理论计算方法下的性质肽构象异构体由于其自身的灵活性会产生大量不同的异构体,这样就会导致即使同时使用较高理论水平和较高计算基组研究所有可能的几何异构体也是不现实的。随着体系规模的增大,高级方法不仅不利于扩展,而且理论水平越高,需要较大的基组才能获得可靠的结果。在量子化学中,由于采用的是有限基组近似,故在计算一个超分子比如二聚体AB的相互作用能时,会由于描述片段A的基组扩散到B上,反之同理。这样,描述片段A,B二聚体分子的基组就不处于同一个水平了--基函数的数目发生了一定的变化。这样会导致计算二聚体能量时所用的基函数数目增多,人为地降低了二聚体的能量,从而导致相互作用能的计算偏大,这种现象称为基组重叠误差。之前有学者研究过酪氨酸-甘氨酸二肽化合物最稳定的构象“bookl”,随后实验者发现BSSE效应在三种不同理论计算水平下的数值分别是 28.63(MP2),6.27(M06-2X)和 4.28(B3LYP-D3)kJ/mol。对于B3LYP理论水平来说因为其不考虑伦敦色散力影响,基本不考虑BSSE效应。由于BSSE导致分子不同部分之间产生吸引力,这可能使分子产生更多折叠和紧凑的构象。由于分子内的BSSE难以计算,所以本文的研究重点不是为了减少BSSE值,而是更好的理解分子内的BSSE影响。虽然不少研究也在探讨修正分子内的BSSE,但是还没有一种独特的方法可以真正做到。目前关于了解这种BSSE效应的方法通常是关注其能量最小值的位置以及从原子几何学中获得关于键合的信息。除QTAIM方法外,通常用几何结构’紧凑’程度来作为判断方式。在本节中将用QTAI M和应力张量分析法提供出较丰富的拓扑描述。通过跟踪键临界点(B CPs)和笼临界点(CCPs)随着二面角 Cβ(Tyr)-Cα(Tyr)-Ccarb(Tyr)-N(Gly)变化的数量可以得知分子之间的键合程度以及分子的紧凑程度。另外利用总局部能H(rb)以及BCP应力张量λ3σ来判断BCP的稳定程度。最后测量跟踪由于扭转键的扭曲而导致的键余长,这种由于键的扭转而具有的长度我们给它命名为键路径的螺旋度长度,符号表示为H;借此用来比较四种不同理论水平之间的不同特性。本工作使用Gaussian 09程序扫描四种理论水平下的不同二面角ψtyr对应的势能面,使用相同的计算基组6-31+G(d),扫描范围为80°到110°,以0.5°为步长。随后使用AIMAl]程序来辅助计算分子和键临界点的性质,所有的分子被证实没有非核吸引子临界点。由于BSSE的已知效应是分子内不同部分的吸引力导致分子产生过度结合的结构。选取在能量最小值θ=96°时不同理论计算水平下生成的3-DQT量子拓扑的分子图做研究,相对能的图展示在补充材料S1中。不同理论水平下的键合程度用BCP(b)的数量对扭转角度作图,见图2-2,由图可知MP2,M06-2X计算产生的分子中的BCP的数值的最高值包含36个BCP,而由B3LYP计算得分子中最多包含34个BCP,对于这四种理论计算水平来说B3LYP分子中则有着最少的BCP值,最小值为b=33。B3LYP-D3理论水平的BCP数量介于34-35之间。这种结果与采用非QTAIM法得到的结论一致,MP2理论计算水平与B3LYP相比往往会造成二肽构象结构更多的折叠。文中用CCP的数量随二面角角度的变化来佐证分子结构紧凑程度。在B3LYP理论水平下会产生开放的2-DQT分子结构,由B3LYP预测到的开放结构可能是由于分子内缺乏色散力作用或者是由于MP2和M06-2X计算过程中会存在大的BSSE影响。对于B3LYP-D3下产生的较紧凑结构,由于不考虑其BSSE效应的影响,造成此结果的主要原因是色散力。由能量最小值θ=96°附近的四个理论水平的局部总能量密度H(rb)表示的扭转C2-C3 BCP的强度增加顺序为MP2BPL。即使H*和H的值相等或者相似时,由于键路径p和q遍历不同的空间区域,所以真实的键路径p,q也是不同的。对键路径长度的表示方法进行相应的拓展,定义了无量纲表示法即 Hf*,Hf,Hf0min,Hf0min*,其方程如下 Hf=(H-BPL)/BPL,Hf0min=(H-Hθmin)/Hθmin。在探究不同表示方法的差异时,以扭转角θ作为绘图的横坐标,不同特征向量上扫描得到的键路径长度(H*,H),(Hf*Hf),(Hf0min,Hf0min*)作为纵坐标进行分析,详情见图3-1,3-2和3-3。首先比较不同对位取代基的联苯的相对能,随后将键路径长度随角度的变化与之做对比,见图3-1。其中BPL随角度变化的也绘制在其中,用黑色线条来表示。在图中用彩色图标来表示不同取代环境下联苯中的强扭转键C4-C7 BCP的螺旋度长度。从图中可以得出相对能随角度变化的趋势与特征向量e2所对应的键路径长度H相似,同时观察到在二面角θ=90.0°的位置时对位取代联苯的相对能和螺旋度长度均存在着较大的间隔,在θ=40.0°时,对位取代联苯的相对能和螺旋度长度之间的差距小于上述间距。当θ=40.0°时H的值处于最小值,此时H的变化趋势与相对能的变化趋势相似,但是此位置下对应的BPL值则与相对能的变化趋势不同,BPL值之间的差距比相对能之间的差距更大。对比二面角C4-C7 BCP的键路径长度H和H*见图3-1和补充材料4(a)可知二者不具有区分度,这是由于其没有明显的键曲率所导致的。非线性键合路径被定义为(BPL-GBL)/GBL,它更易发生在闭壳BCP中而不是开壳SCP。当键路径r伴随有非零的键曲率时,会导致键合路径长度H>H*这是因为在特征向量e2与e1的方向上相比是电子密度累积最多的方向,即最优选方向。H15---H20 BCP则与C4-C7 BCP恰好相反,在图3-2中显示出了较明显的键曲率,其特征向量键径长度H和H*显示出了不一样的变化,从图中能够得出对于不同取代联苯下的H值都大于BPL,由方程qi=ri+εie2可知,当椭圆率ellipticit ε不等于零时可推出如上结论。具有较短的]H值的相对能值也较低,这是因为对于一个含有较小椭圆率的BCP而言,它在随着扭转角变化时产生更小的阻碍作用。接下来对连接两个苯环之间的弱相互作用H15---H20 BCP进行重点的讨论。研究发现当二面角θ≈22.0°时该键会自动断裂,键路径长度]H和H*在作图时采用彩色图形表示,BPL用黑色图形表示,数据显示在0.0°≤θ≤20.0°的范围内显示出了与BPL变化完全不同的趋势,H值在这样的角度范围下逐渐升高,对于BPL的变化则不太明显。其中特殊的位置θ ≥ 17.0°时螺旋度长度H和H*的值均出现了分离的变化;相反,对于BPL而言则无论θ的值为多少都不会出现值得明显分离;比较键路径长度H和H*值可以得出在H15---H20 BCP破裂之前H>H*,当θ ≥ 0.0°时H---H BCP开始在最优选选e2向最近的RCP转变,RCP用红色的球体表示,见图Scheme 3-1。接着检查第一种分数形式的特征向量对应的键路径长度Hf和H*f,见图3-2。对于 C4-C7 BCP 和 H15---H20 BCP 二者的Hf和 Hf*数值具备可区分性。Hf值的数量级要远远大于键曲率所对应的数量级由图3-2可知,造成这个现象的原因是由于二面角的存在破坏了特征向量e1和e2相对于e3的相对取向。最后分析了第二种无量纲形式的特征向量对应的键路径长度Hfθmin 和 H fθmin*,从图3-3中可知对于弱相互作用的H15---H20 BCP而言其数值不可区分。这与]Hf和Hf*数值可区分性存在着对比。另外从此还得得到一个结论对于强扭转键C4-C7 BCP而言其]Hf和Hfθmin的变化有着很大的差别。通过上述研究提供了一种新的探究化学键的思想,在将来可以用键路径框架集合组]B={p,q,r}对化学键中三维方向的变化进行详细阐述。由于p,q,r各自有其对应的键路径长度H*,H和BPL,当C4-C7 B和H---H BCP处于旋转角的环境时由H*,H量化得到的键路径长度明显要比BPL在这范围内的变化要敏感的多,当分子受扭转角影响时实际的键长变化不是很明显,但是由于受扭转角的影响p和q路径相对于r的取向会发生改变。这点足以体现出当分子处于扭转环境中时将特征向量e1,e2和e3结合考虑来定性和定量的描述化学键性质。由于在实际的化学反应过程中分子中的化学键并不只是受拉伸和压缩的影响,而是包含一定的扭转,对于前者来说对化学键的影响则是只影响BPL的大小。因此键路径框架集合组的方法对于化学键进行深刻而详细的分析是有意义的。五、应力张量特征向量法拓展出新一代的QTAIM理论电子应力张量的特征向量对于化学反应性的预测是有效的。它特别适合用来理解化学反应过程中分子中的电子结构变化。量子应力张量σ(r),是直接与维里定理的Ehrenfest force相关的。它可以为低频正常模式下的分子往往会伴随着结构重排的现象而提供了物理解释。在这项工作中我们由Bader定义的应力张量去探究QTAIM分区中的应力张量的相关性质。在三维空间中流动的微小流体立方体Π(x,y,z,t),这是一个秩为3的张量场。张量场中有九个分量,其中处于对角线分量中的Πxx,Πyy,和Πzz所对应的是正应力,这些分量中的负值表示的是多维数据集的压缩。相反正值表示拉力或张力,其中更负/正值对应于正方体所受到的更大的压缩/张力。对角线化的应力张量σ(r),用Πxx,Πy y,和Πzz来量化,我们得到结论在QTAIM分区中应力张量的特征值λ1σ。λ2σ,λ3σ是与特征向量e1σ,e2σ,e3σ相对应的。用QTAIM和应力张量理论来解释特征值时是不同的。QTAIM方法定义最优选的方向是电荷密度最易积累的方向。与之不同的应力张量法则认为最优选的方向是最可被压缩的方向,即受到最小拉力的方向。对于应力张量而言,特征向量的特征值是有序的λ1σ<λ2σ<λ3σ,对于特征值λ3σ而言代表着纯拉伸作用,特征值λ1σ则代表着最具有最大的压缩作用。与之相比QTAIM理论中特征向量的特征值呈现λ1<λ2<λ3的排序,对于特征值λ2而言代表着比λ1更可变的作用。这便得到一个结论QTAIM方法中的特征值λ2和应力张量方法中代表压缩作用的λ1σ。具有可比性。通常而言情况下应力张量方法中的特征向量e1σ和e2σ和QTAIM方法中的特征向量e1和e2是不统一的,尤其是对于对称键,例如链接两个苯环的 C-C BCP。探讨在QTAIM理论中如何使用应力张量,要对应力张量的性质即特征向量e1σ和e2σ进行直接的探究,就不能仅仅只是从特征值λ1σ。和λ2σ的角度去考虑。首先从应力张量的椭圆率εσ=|λ2σ|/|λ1σ|-1和εσH=|λ1σ|/|λ2σ|-1入手进行分析,下角标’H’代表的是Hessian分子和分母排序。从相应的理论公式可以推出应力张量的椭圆率εσH≥ 0,ε≤0且ε≥0,见4-1和4-2。QTAIM理论中椭圆率ε提供了在BCP位置处垂直于键合路径的两个方向上的ρ(rb)的相对积累,λ1和λ2分别是相应的特征向量e1和e2的负特征值。之前实验假设电荷密度运动方向中QTAIM法中最不优选方向e1和最优选方向e2与应力张量中e1σ。和亡。存在着一一对映关系,但是之后的探究我们获悉QTAIM的特征向量和特征值与应力张量映射之间缺乏一对一的映射。出现这种混乱的原因是对于非对称键中远离几何中心的BCP而言QTAIM中的e2方向与应力张量中的e2σ方向会有偶然的重叠。伴随着QTAIM分区内的应力张量法在不对称BCP中的应用,发现QTAIM法中会导致▽·ρ(rb)=0,而应力张量中-▽·σ(rb)≠0。因此,在这项工作中要做的重点是把QTAIM和应力张量理论应用在整个键径上进行分析,用以避免这种仅仅只在键临界点BCP分析而得出的错误的结论。当出现▽·ρrb=0和-▽·σ(rb)=0不匹配的情况时将会导致甶应力张量的特征值决定的应力张量的性质对较小的变化变得很敏感。比如使用λ3σ<0的标志可以用来作为判断不稳定或者接近相变的有用工具。本研究的计算细节如下,用Gaussian 09B01对乙烯分子的势能面进行扫描。选用计算方法以及对应基组B3LYP/cc-pVQZ来计算每一步的单点能,扭转角的角度由二面角H3-C2-C1-H6来定义,角度范围是-180.0°≤θ≤180.0°,设定其势能面扫描步长为1.0°。使用AIMAll程序分析生成的波函数,获取相应的计算信息。最后使用Topviz程序可视化出相应的键路径p-,q-。由图4-2(b)可以发现三种不同形式的椭圆率在扭转角θ≈±150.0°的附近位置上会有一个明显的变化,但是这种变化并没有体现在相对能的图中。对于C1-H3 BCP的椭圆率变化中εσH和ε的变化很相似,就会产生一个一般性的思想,认为标量下的QTAIM和应力张量的行为是相似的。由此得出一个假设除了对称性的键路径之外QTAIM法和应力张量法对应的性质也是相似的。为了验证以上的结论紧接着检验了 C1-C2 BCP所对应的三种椭圆率,如图4-2(a),然而却得到了与上述结论不一致的结果。因此如果仅仅使用BCP椭圆率ε,εσ和εH是不足以用来比较QTAIM和应力张量的。因此,选用矢量应力张量轨迹Tσ(s)来作为确定应力张量法最优选方向的测量方法。当扭转角施加在C1-C2 BCP上时会使乙烯分子的能量偏离最小相对能的位置,在图4-3上可以清楚的看到逆时针方向和顺时针方向的变化。从表4-1中可以看到对于CW方向来说,与最优选方向相比而言,在最不优选方向上往往有着最大的投影值。即(e2σ.dr)max>(e1σ·dr)max,在补充材料中7中包含着CCW方向的详细数据。这足以表明e1σ特征向量指代的是最优选的电子运动方向。换句话说,e1σ特征向量表示最可压缩的方向。前面章节提供了一种基于新一代分子中的原子理论键路径框架集,B={p,q,r}。由此在应力张量理论中继续对其进行拓展,由于应力张量的椭圆率有两种表示方式εσ,和εσH,所以应力张量的键路径框架集与QTAIM不同的是其有两种不同形式的表示方式,BσH={pσH,qσH,r}和Bσ={pσ,qσ,r}。在上个篇章中提到的与BCP相关的三条基本路径p,q和r,在本章节把它们都做到了可视化,其中路径p,q和r分别用淡蓝色线条,洋红色线条以及黑色线条表示。可视化的操作可以直观的看到应力张量参数e1σ,e2σ以及QTAIM中的参数e1,e2的关系。如图4-6。通过对这三种形式的键路径框架集进行对比分析,可以找到最恰当的应力张量形式。特征向量对应的键路径长度(H*,H),(HσH*,HσH)和(Hσ*,Hσ)是直接与(p-,q-),(pσH-,qσH-)和(pσ-,qσ-)路径相关的,可以作为辨析应力张量与QTAIM是否具有相似性的方法,见图4-4。由图可知C1-C2 BCP键路径长度均大于由键径(r)确定的键径长度。但是对于C-H BCPs而言与(pσ-,qσ-)和(pσH-,qσH-)相关的(Hσ*,Hσ)和(HσH*,HσH)键径长度则有的比由键径(r)确定的键径长度要短。同时可以看出应力张量特征向量所对应的键路径长度(Hσ*,Hσ),(HσH*,HσH)伴随着角度的变化比QTAIM法中键路径长度(H*,H)要更加靠近键路径r。对于C1-C2 BCP来说应力张量特征向量的两种形式的键路径长度的数值在某种程度上难以区分。因此得出结论,对于这两种标量下的(Hσ*,Hσ),(HσH*,HσH)都不可作为与QTAIM中键路径长度(H*,H)有着相似变化的近似方案。从其中同样可得知的信息是C1-C2 BCP,C1-H3 BCP,C2-H6 BCP的键路径长度的值均为H*