研究小分子与生物大分子的结合:量子化学方法

研究小分子和生物学大分子的结合或相互作用，不外乎两种方式：以共价方式结合，以配体-受体方式结合。前者涉及到了化学键的断裂和生成，需要考虑化学知识。量子化学可帮助理解这种结合方式，所以这段时间参加了几次workshop。逐渐建立自己的量子化学知识体系。至于配体-受体结合方式，涉及到氢键、范德华力、库仑力等弱相互作用，多用分子力学知识处理。

目录

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
1.量子化学基础知识
2.量化计算基础知识
 2.1. 理论theory
 2.2. 基组basis set
 2.3. 特殊效应
3.Gaussian计算
 3.1. 能量相关计算
 3.2. 过渡态TS、IRC等
 3.3. 结构相关计算
 3.4. 分子性质计算
 3.5. 光谱计算
4.波函数分析与Multiwfn
 4.1. 分析理论
 4.2. 电子密度electron density(和电子定域性)
 4.3. 轨道分析orbital analysis(和轨道成分、轨道定域化)
 4.4. 布居分析population analysis(和原子电荷)
 4.5. 键级分析bond order
 4.6. 属性分解分析
 4.7. 实空间函数及可视化
5.专题:定量分子表面分析(预测结合方向、活性位点等)
6.专题:弱相互作用分析

1. 量子化学基础知识

薛定谔方程(Schrodinger) (仅适用于速度不太大的非相对论粒子)与波函数

  1. 设描述微观粒子状态的波函数为Ψ(r,t)，质量为m的微观粒子在势场V(r,t)中运动的薛定谔方程。在给定初始条件和边界条件以及波函数所满足的单值、有限、连续的条件下，可解出波函数Ψ(r,t)。
  2. 波函数ψ，相应的能量E:
  3. 要完整描述电子状态，必须要四个量子数: 主量子数n、角量子数l、磁量子数m、自旋磁量子数
  4. 经典力学波函数也写成Ψ(x,t)，特别是光的波函数，是理解微观粒子波粒二象性的开始。
  5. 物质波的波函数是概率波，即没有实际意义，不描述具体的运动形式，而是反映粒子的状态；但其模的平方代表粒子在该处出现的概率密度。 既然是概率波，那么它当然具有归一性。
  6. 求解分子体系的Schrodinger方程：Born-Oppenheimer近似(视原子核固定，作为经验出现)，并且数值求解。所以计算精度在于理论和基组。

分子轨道理论(MO): 分子轨道是原子轨道的线性组合；离域性很强，没有严格的物理意义；半经验、HF、DFT的轨道都是MO，Post-HF没有MO的概念

分子体系总的的电子特征

  • 有无单电子
  • 开壳(有单电子)和闭壳(无单电子):
  • 自旋多重度: 2S+1; S是总自旋量子数
    • 只含成对电子 S=0, 自旋多重度＝1 （单线态）
    • 含一个不成对电子S=1/2, 自旋多重度＝2 （双线态）
    • 含两个不成对电子 S=1, 自旋多重度＝3 （三线态）

常见单位
1 a.u. = 1 Bohr = 0.529 A 1 nm = 10 A
1 a.u. = 1 Hartree = 627.51 kcal/mol = 2625.5 kJ/mol (1 kcal=4.184 kJ)



2. 量化计算基础知识

量子化学计算的关键在于理论方法和基组这两个部分

2.1. 理论theory
  • Semi-empirical:是HF方法的近似，e.g. Huckel分子轨道法, ZINDO, PM6D3, PM7；除了此法，剩下的是ab initio/第一性原理计算
  • Hartree-Fock(HF): 仅具有理论意义，实际不用
  • DFT(density functional theory): 计算花费和HF一样，但精度达到Post-HF，是主流. e.g. B3LYP, M06-2X(M062X), B2PLYP, PWPB95
    • [DFT]-D3是用于色散矫正，在研究弱相互作用!!!v(如氢键、卤键、VDW)时，如B3LYP-D3, M062X-D3, B2PLYP-D3
    • 计算激发态有TDDFT
  • Post-HF: 组态相互作用(CI), 微扰(MP), 耦合簇(CC)等多种方法，e.g. MP2, CCSD(T)(高精度首选，也不需要色散矫正)

总结：高精度用CCSD(T), 较高精度用双杂化泛函，一般精度用普通泛函，低精度用半经验
计算Tips: 单个结构能量：几个十几个用CCSD(T); 几十等原子体系用DFT；几百几千原子体系用半经验方法，如与分子动力学模拟结合等；动力学模拟时：十几个原子用DFT；然后用全原子力场(AA force filed)；更多的用粗粒化方法

2.2. 基组basis set

基函数: 每个基函数由多个高斯型函数(Gauss Type Function)组成(也会用到Contracted GTF)
弥散函数: 延伸到非常广的空间,对色散作用主导的弱相互作用等。不需要时绝对不能加!!!
极化函数: 6-31G, 6-31G** 使得原子价轨道在空间取向上更“柔软”，使之容易成键
扩展基: 2-zeta, 3-zeta, n-zeta: 每个价层原子轨道由n个基函数描述

基组类型:

  • Pople系列: 6-31G; 是所有计算的底限; 6-31G*, 6-311G，通用各种理论方法
  • def/def2系列: def2-TZVP def2-QZVP，适合DFT，当然post-HF也可以
  • Dunning系列: cc-pVTZ cc-pVQZ，专门适合Post-HF
  • 赝势基组：处理内层电子，对重金属元素，如对过渡金属用lanl系列基组
  • 混合基组：如3-21G; 有机大分子，不重要的部分可以用此，并做位置限制

弥散函数，极化函数，相对论效应(内层电子)
不同基组下计算的能量没有可比性！！！

2.3. 特殊效应

相对论效应, 溶剂效应, 周期性



3. Gaussian计算

量子化学计算只能设定单个体系总的电荷和自旋多重度
输入文件格式:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 %nproc=2
 %mem=1400MB
 %chk=C:\x.chk
 %rwf=C:\x.rwf
 # B3LYP/6-31G* geom=check guess=read out=wfn
 Title
 0   1
 O                 -0.64410480    0.28384279    0.00000000
 H                  0.31589520    0.28384279    0.00000000
 H                 -0.96455939    1.18877862    0.00000000
C:\x.wfn

Note:
geom=allcheck 输入结构从chk中读取(%chk指定的)
guess=read 从chk文件中读入波函数为初猜
out=wfn: 输出wfn文件
Post-HF输出时，需要额外写density
高斯输出文件: fch(Gaussian私有格式)、wfn、out、chk文件和rwf文件

可以通过linux命令监控log文件:
tail -f x.out
tail -f x.out | grep -E ‘keyword1|keyword2’ -A 2 -B 2 -C 4
tail -f x.out | grep -e ‘keyword1\|keyword2’ -A 2 -B 2 -C 4

3.1. 能量相关计算
  • 单点能single point: 是指电子的能量，不包括电子动能、电子间互斥能、电子与原子核的吸引能、核间互斥能。，与自由能不同，
  • 势能面potential energy surface: 可进行势能面扫描
3.2. 过渡态TS、IRC等

IRC是连接势能面上两个能量极小点的最低能量路径。而TS是其能量最高点。

3.3. 结构相关计算

几何优化geometric optimization: opt(tight|loose, maxcyc=n)命令

有两个条件：1. 各原子受力皆接近0；2. 振动分析时没有虚频

  1. 在优化过程中，观察每一步优化的结构: 在Open时选择Read intermediate Geometries选项。若不收敛，赶紧停止。
  2. 几何优化对基组不敏感，可以用一般，如B3LYP/6-31G*；而做单点能计算时用好的，如CCSD(T)/cc-pVTZ
  3. 同时做freq振动分析，因为没有虚频的几何优化结果才正确
  4. 连续两种方法进行几何优化时，除了用geom=check, guess=read从chk文件读取初始几何构象和波函数外，还要用opt(restart)。坐标部分可以删除
  5. 几何构型优化不收敛：
    首先要分清scf不收敛和几何构型优化不收敛：scf不收敛指的是自洽场叠代不收敛，可以认为是对指定结构的波函数不断优化的过程，是为了找到这个某个指定结构下能量最低的波函数，而几何构型优化是对结构的优化的过程，是为了找到某个指定的组分下能量极小结构（注意，不一定是能量最小结构）。在量子化学计算的几何构型优化中，每一步的几何构型优化都包含的很多次的scf计算

振动分析freq: 可以用geom=check, guess=read从chk文件读取初始几何构象和波函数外

3.4. 分子性质计算

热力学量
opt freq 计算热力学量默认的条件是temperature=293K pressure=1atm; 较大的分子可以采取预优化+正式优化的策略!

1
2
3
4
内能(U):   Utotal=Utrans+Urot+Uvib+Uele; 即体系的平动、转动、振动、电子
焓(H):     H=U+PV
熵(S):
吉布斯自由能(G): G=H-TS

  • 方法一：

    1
    2
    3
    4
    5
    6
    7
    8
    0K时:存在电子基态能量(单点能)、ZPE(零点振动能)
          U(0)=εele+ZPE
          H(0)=G(0)=U(0)
    >0K时:U(0)|H(0)|G(0)加上温度(298.15K、1atm)对其各自的贡献
          U(T)= εele+ZPE+ΔU0→r
          H(T)= εele+ZPE+ΔH0→r
          G(T)= εele+ZPE+ΔG0→r
    后两者成为热力学校正能量，在B3LYP/6-31G*下就可以；单点能在高级别条件下计算。二者相加即为所求。

  • 方法二：热力学组合方法，以CCSD(T)/CBS为目标，耗时长;有Gn系列和CBS系列

3.5. 光谱计算

红外光谱与拉曼光谱, 电子光谱, NMR



4. 波函数分析与Multiwfn

首先通过Gaussian获得波函数: 用Post-HF法时，Gauss中写上density关键词，才能输出波函数信息。
波函数分析=分析方法+分析程序+实际问题

4.1. 分析理论
  • AIM理论: Atoms in molecules。电子密度、拉普拉斯函数、动能|势能|能量密度、拓扑分析、盆分析、定域化|离域化指数、范德华表面
  • 概念密度泛函理论: 涵盖前线轨道FMO理论，电负性、软硬度、Fukui函数、双描述符、亲电|亲核指数、反应活性位点
  • NBO理论: Natural Bond Orbital，定义了很多轨道
4.2. 电子密度electron density(和电子定域性)
  • 薛定谔方程与Born-Oppenheimer近似
  • 波函数(wave function)和密度矩阵(对体系波函数的一种紧凑描述)
  • 电子密度: ρ(r), 直接衍生于波函数(模的平方)

电子密度、价层电子密度
电子定域性: 拉普拉斯值函数、ELF函数、LOL函数 (通常用于分析强相互作用;成键、孤对电子、原子壳层结构)

1
2
“化学键、孤对电子、sigma/pi/delta键、轨道杂化、lewis结构、VSEPR这些经典理论诞生于大半个世纪前，属于“概念化学”范畴。这些模型化的表述从当今量子化学角度看似乎缺乏根据，是一些很粗俗的观点，但实际上它们并非是完全人为臆断的，它们从量子化学的角度是可以找到充分依据的，定域化分析正是力图通过波函数或其衍生信息还原这些信息。”
“使用实空间坐标为变量的函数来描述电子定域性特征，这不需要明确地计算定域化轨道，因而避免了上述问题。这些函数将复杂的3N维波函数所携带的关键信息约化、抽象成三维函数，从这些函数在分子附近空间中各点的数值可以分析分子各个位置定域性特征的强弱，使化学键、孤对电子、原子壳层这样电子定域性强的区域充分凸现出来。”

1
2
3
4
5
6
拉普拉斯值函数的值被定义为r处电子密度Hessian矩阵的迹，也就是拉普拉斯算符作用到电子密度的结果，亦即电子密度在三个正交方向曲率和。拉普拉斯值函数在某处为正，说明此处电子密度以发散为主；若为负，说明以聚集为主。非极性共价键、孤对电子处由于电子呈聚集状态，因而它们都处于负值区域；极性键也有电子凝聚区域，相对于非极性键发生了偏移、变形。
[动能密度、势能密度、能量密度]
ELF(Electron Localization Function，电子定域化函数): ELF是个三维实空间函数，数值范围在0~1之间(一般关注ELF>0.5的区域)。简而言之，数值较高的ELF的等值面包围的区域，电子在里面定域性较强，不容易跑出去，相对地，电子能够容易地在这样的区域内随意运动（在此区域里随意离域）。在那些ELF数值较低的区域，电子定域性弱，即如果把电子放在那里，就很容易离域到其它区域去。能跑到以及容易跑到哪里，可直接通过考察费米穴的分布来衡量。
LOL(Localized orbital locator，定域化轨道指示函数): 数值范围在0~1之间0-0.75. 轨道定域化有随意性；但电子密度定域化(离域化)有确定的物理意义。完美的电子定域化域里，里面电子不会出去，外面电子不会出来。有些区域电子定域性很强(如内核区域)，有些地方相对较高(共价键区域、孤对电子区)
4.3. 轨道分析orbital analysis(和轨道成分、轨道定域化)
  • 原子轨道组合成分子轨道MO
  • 分子轨道(MO): 很强的离域性，而且与化学键、孤对电子等概念没有直接对应
    • Canonical MO (CMO, 半经验、HF和Korn-Sham DFT获得的)：分为内核MO、价层占据MO、价层非占据MO、高价非占据MO
    • Natural orbital (NO, Post-HF获得的)
  • 自然键轨道NBO: 比MO高度定域，但不能用于弱相互作用
    MO原始密度矩阵(CMO, NO) --- PNAO --- NAO(自然原子轨道)--- NAO下密度矩阵 ---[基于NAO的各种分析]--- PNBO --- NBO(自然分子轨道) --- [基于NBO的各种分析]---LMO --- MO (P是未正交化的)
    NAO(Cor, Val, Ryd)和NBO(CR, LP, BD, BD*, RY)有对应关系
    NBO已经高度定域化，但占据数是非整数，进一步LMO(localized MO)

MO可视化
Multiwfn-0 Show molecular structure and view orbitals
Multiwfn-6 Check & modify wavefunction
Multiwfn-200下有3.Generate cube file for multiple orbital wavefunctions

NBO分析:Natural bond orbital

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# B3LYP/6-31G* opt   pop=nbo6 external=’nbo6.bat’
Title
0   1
 O                 -0.64410480    0.28384279    0.00000000
 H                  0.31589520    0.28384279    0.00000000
 H                 -0.96455939    1.18877862    0.00000000
----------------------------------------------------------
# B3LYP/6-31G* opt    pop=nbo6read   external=’nbo6.bat’
Title
0   1
 O                 -0.64410480    0.28384279    0.00000000
 H                  0.31589520    0.28384279    0.00000000
 H                 -0.96455939    1.18877862    0.00000000
$NBO  plot file=filename $END

Note: NBO基于密度矩阵，若用Post-HF法，必需用density关键词，否则给出HF级别的密度矩阵; Opt和nbo关键词一起用时，先对初始结构做NBO分析，再对最终结构做NBO分析; NBO程序输出文件 .31包含基函数信息，.32~.40，分别是PNAO,NAO,PNHO,NHO,PNBO,NBO,PNLMO,NLMO,MO,NAO(Cor,Val,Ryd)和NBO(CR,LP,BD,BD*,RY)有对应关系

FMO分析: 前线轨道HOMO和LUMO;可以对它们进行轨道成分分析

轨道成分分析

  • 分析MO上任意一个原子或片段的贡献
  • 成键轨道是由哪些原子轨道组成
  • 预测反应活性(前线轨道理论): 分析HOMO和LUMO上原子的贡献。占HOMO(LUMO)成分越大的原子越有可能成为亲电(亲核)位点。
    计算原子或片段对MO的贡献，首选Hirshfeld或Becke;计算原子轨道对MO的贡献，首选NAO;计算基函数对MO贡献， 首选SCPA

轨道定域化和AdNDP分析

  • 定域化分子轨道 localized molecular orbit(LMO): 这类轨道多数定义在单、双键中心，可以和化学键、孤对电子等概念对应。
  • CMO得到LMO：LMO和正则MO轨道完全等价，而且LMO不同的方法形状可能不同，但物理性质完全一样。Pipek-Mezey LMO和Boys LMO
  • NLMO得到LMO: NAO—NBO—LMO。不同于MO，NBO已经不是离域化，而且更轻松得到LMO
  • LMO有附加功能: 氧化态计算，LOL，ELF
    CDA分析和FO
    一般是原子轨道组合成分子轨道，这儿用片段轨道(Fragment Orbital)组合成分子轨道
    可以分为donor和back-donation
    DOS图
4.4. 布居分析population analysis(和原子电荷)

Population analysis用来获得不同基函数、原子轨道、壳层、原子、片段所带有的电子数(布居数)
原子电荷=原子核电荷数-原子布居数
用途: 亲电试剂和亲核试剂、反应位点、静电势等。
计算: DFT理论/6-31G*就可以了
电荷计算charge

  • 直接基于波函数方法
    • Mulliken电荷: 快，不随基组增大而收敛，所以没法判断哪个基组下是较为准确的
    • NPA电荷: Natural population analysis被误认为NBO电荷
  • 拟合静电势方法
    对静电势(ESP)问题效果最好(用于力场较好)，多用于分子动力学模拟，但可移植性差，依赖构象(不适合柔性构象的分子模拟，所以分子对接用什么电荷???)
    范德华表面的静电势
    • MK|Kollman电荷:
    • CHELPG电荷: 推荐用于刚性分析的MD
    • RESP电荷: Amber力场的标准电荷模型，推荐用于柔性分子的MD。AM1-BCC电荷是快速估计RESP电荷; 类似AM1-BCC的是MMFF96电荷(用于量化算不动的巨大体系)，后两者都对静电势重复性很好。
  • 基于实空间下电子密度的方法
    • AIM电荷:
    • Hirshfeld电荷: 预测反应位点、反应活性较好，但对静电势重复很差，所以在MD中很差
      ADCH电荷: Hirshfeld的改进，有其各种优点，同时静电势很好重现
    • Becke电荷:
  • 基于电负均衡电荷
    • Gasteiger电荷:
4.5. 键级分析bond order

衡量原子间相互作用特征的最简单、最明确的指标；反映原子间共享电子对数或作用强度
Mayer键级(MBO): 最普世的键级计算方法
Laplacian键级(LBO):
多中心键级:

4.6. 属性分解分析

如能量分解、电荷转移分解，轨道成分分解

4.7. 实空间函数及可视化

各种实空间函数都是基于轨道的(可以计算各个轨道对实空间函数的贡献)
轨道波函数、电子密度(以及电子密度径向分布函数)、自旋密度、电子拉普拉斯函数(Laplacian of density)、电子定域化函数(ELF)、定域化轨道指示函数(LOL)、静电势(ESP)、价层电子密度(Valence density)、约化密度梯度(RDG)

  • 实空间函数: 轨道波函数
  • 实空间函数: 电子密度与定域性
    电子密度: 与价层电子密度()
    电子定域性: 拉普拉斯值函数、ELF函数、LOL函数 (通常用于分析强相互作用;成键、孤对电子、原子壳层结构)
  • 实空间函数: 分子静电势ESP
    见定量分子表面分析. ESP的可视化有两种:
    • 静电势等值面: 显示出来的面上的每个点的静电势值都与设定的静电势的isovalue值相同
    • 静电势映射的分子表面: 常用范德华表面

实函数的积分(用于定量分析)
对这些函数在各个原子空间内进行积分，这样每个原子上就有一个数值，便于不同原子之间进行定量比较。这样转化后，称呼也变了，比如：

  • 1
    2
    3
    4
    5
    电子密度 --> 原子的电子布居数，或简称原子布居数 (Atomic population)
    自旋密度 --> 原子的自旋布居数，或简称原子自旋布居 (Atomic spin population)
    福井函数 --> 简缩福井函数 (Condensed Fukui function)
    双描述符 --> 简缩双描述符 (Condensed dual descriptor)
    Parr函数 --> 简缩Parr函数 (Condensed Parr function)

实函数的可视化
输出图像格式以及像素数等，都是可以在settings.ini中定义

  • 点/线/面的性质：在一维或二维上考察各种实函数值的分布, 如轨道波函数、电子密度、电子密度拉普拉斯函数、自旋密度、定域化轨道指示函数(LOL)、电子定域化函数(ELF)、静电势(ESP)、平均局部离子化能(ALIE)、约化密度梯度(RDG)等。
    Muliwfn-1 Output all properties at a point
    Muliwfn-3 Output and plot specific property in a line
    Muliwfn-4 Output and plot specific property in a plane
    Muliwfn-5 Output and plot specific property within a spatial region (calc. grid data)

  • 分子表面：分子表面没有明确的定义。最常用的范德华表面根据电子密度等值面(Bader 所定义, 即电子密度为 0.001 a.u.的等值面)定义的，还有溶剂可及表面积(SAS)

  • 体数据：即格点数据，cube文件
    Multiwfn-5是Output and plot specific property within a spatial region (calc. grid data)，专门输出各种实函数的格点数据
    Multiwfn-13 Process grid data可以对cube文件进行各种后期处理。
    Multiwfn-100下面有个Draw scatter graph between two functions and generate their cube files，可以输出各种实函数的cube格式

  • 初属性和变形属性的绘制: 任何属性都有初属性(Promolecular property)和变形属性(Deformation property)
    实际属性=初属性+变形属性
    Muliwfn-3 Output and plot specific property in a line
    Muliwfn-4 Output and plot specific property in a plane
    Muliwfn-5 Output and plot specific property within a spatial region (calc. grid data)里面都有初属性和变形属性的计算

实函数的拓扑分析和盆分析
对电子密度、ELF、LOL等各种各样的实空间函数进行拓扑分析和盆分析

  • 拓扑分析: 电子密度及其拉普拉斯函数的拓扑分析是AIM分析中的核心，而ELF、LOL的拓扑分析在分析成键方式、芳香性等问题上也有重要应用
    电子密度拓扑分析(即AIM拓扑分析):
  • 盆分析: 最早在AIM理论中对电子密度进行的，可以划分原子空间，还可以对ELF等其它各种实函数进行盆分析。可以产生零通量面(zero-flux surface)|盆间面(inter-basin surface)和吸引子(attractor)。获得AIM盆后，可以(1)对盆里面的实函数进行积分，如积分电子密度就得到盆里面的电子布居数;(2)计算盆内的电多极矩;(3)计算盆内的定域化指数和盆间的离域化指数;

5. 专题:定量分子表面分析(预测结合方向、活性位点等)

定量分子表面分析对于预测反应位点、预测分子间结合模式、预测分子热力学性质有重要意义。

  • 主要是ESP和ALIE在分子范德华表面的分布。分子范德华表面的定义非常多，最常用的是Bader的定义，也就是对于气相分子，使用电子密度为0.001 e/bohr^3的等值面作为分子范德华表面。
  • 此外指标还有福井函数、双描述符，与ALIE一样，在亲核反应|亲电反应中起作用。
  • HOMO和LUMO原子成分分析也可以预测活性位点
  • [概念的区别: 静电势(整体分子的静电场对电荷的作用强度)、电负性(原子吸引电子的能力)和亲电性|亲核性(从分子获得电子的能力|带负电荷或孤对电子的试剂即亲核试剂对亲电子原子的进攻的能力)]

静电势(ESP): Electronic Static Potential, 大小、面积。该方法可以确定出静电势或平均局部离子化能等实空间函数在分子范德华表面上的极大点和极小点. 通常认为距离分子表面上静电势最小点(最大点)位置最近的原子将是最可能发生亲电(亲核)反应的位点.
最常见的有如下应用(引用卢天老师博文)：

1
2
3
4
(1)在化学反应初期分子间通过静电吸引而相互接近的过程与分子对其外部（范德华表面及更远区域）产生的静电效应密切相关，因此可以通过分析分子范德华表面上静电势的分布来预测在什么位点上最有可能发生何种化学反应，通常认为静电势越负（越正）的地方对应的原子越容易发生亲电（亲核）反应。实际上，通过原子电荷来预测反应位点在本质上是这种静电势分析方法的抽象、简化。
(2)若已知某复合物是以静电为主导的方式结合而成（如氢键、二氢键、卤键等形式），利用静电势可以预测、解释复合物中分子的相对朝向如何、结合的强度大概有多大，这对于研究分子晶体堆积方式、受体配体结合模式、分子吸附等问题很有用。其依据是分子之间容易以静电势互补方式接触，即分子表面静电势正值区域倾向于接触负值区域，而且正值越正、负值越负时这种倾向越强，这样能最大程度降低整体能量。
(3)分子的凝聚相性质，比如蒸发/升华焓、密度、沸点、粘度、表面张力、LogP等等，都与分子间相互作用密切相关，对于极性（或具有局部极性）的分子，静电相互作用通常比范德华相互作用大很多，因此通过分析在分子表面上的静电势，对于预测这些分子的这些性质十分有益。
Politzer等人提出的General interaction properties function(GIPF)就建立了这么一套关系。依靠基于分子表面上静电势分布定义的分子描述符，通过拟合QSPR方程，可以很好地预测上述性质。这些分子描述符包括分子整体/静电势为正/静电势为负区域的表面积，表面上静电势最大值、最小值，表面上静电势的整体/正值/负值部分的平均值/平均偏差/方差，以及它们进一步相互运算得到的几个描述符（具体公式见Multiwfn 2.5手册3.15.1节）。计算并利用这些描述符也是定量分子表面分析的主要组成部分。

通常B3LYP/6-31G**级别的波函数已经能很好地重现分子表面的静电势分布了，继续提高基组和理论级别对静电势影响不大。
惯绘制分子表面的ESP填色图(ALIE图也同样方法):

  1. Multiwfn — Quantitative Analysis of Molecular Surface — ESP (OR ALIE)
  2. Generate: 1. pdb file(represented as dot) OR 2 cube file. Both are Ok!
  3. VMD—可视化

平均局部离子化能(ALIE): Averaged Localized Ionizing Energy.数值越低, 表明此处的电子被束缚得越弱, 因此活性越高, 越容易参与亲电或自由基反应. 若某原子处在数值越小的极小点附近, 则该原子的反应活性越高. 若原子附近的分子表面上没有极小点出现,则表明此原子参与反应的可能性很小.

1
在ALIE的诸多应用中，最为重要的就是它预测反应位点的能力。在分子表面上，ALIE越小处，由于电子束缚得越弱，或者说电子活性越强，就越容易发生亲电及自由基反应。通过定量分子表面分析算法，找出分子表面上ALIE最小的几个点，看这几个点离哪些原子比较近，就能判断哪些原子最可能是反应活性位点。ALIE的极小点容易出现在pi电子、孤对电子附近，因为这些区域的电子被束缚得较弱。可以认为ALIE越小处，局部极化率越大，这与密度泛函理论中的局部软度/硬度概念也是密切相关的。

通过静电势和ALIE在预测反应位点上是互补的(引用卢天老师博文):

1
2
3
4
单靠静电势很多情况都不能正确预测反应位点；而单靠ALIE，在一些情况下也是失败的（而且ALIE也不能用于预测亲核位点），此时必须将静电势综合考虑。这是因为，两个分子相互反应的完整过程分为以下两个步骤。注意这两个步骤没有明确界限，是平滑过渡的。
1. (弱相互作用主导时)两个分子从相距很远逐渐接近：这个过程中分子的电子分布并不发生明显改变，顶多是稍微被另一个分子的电场所极化。长程静电相互作用在这个过程中是关键，通过分析分子表面的静电势，可以很好地估计反应物容易被拉到哪里。也因此，光靠分析ALIE来预测反应位点往往不行，因为即便靠ALIE预测出某些原子容易发生反应，但是若反应物不容易靠静电作用被拉到这里，则在这里发生反应的几率还是不大。
2. (强相互作用时)范德华表面快要接触的两个分子进一步靠近，通过成键、断键这样的电子结构重排过程完成化学反应：这一步在哪里容易进行，关键看哪个原子附近电子活性大，用ALIE最适合研究这个问题。这也表现出单纯靠静电势预测反应位点的不可靠性，因为比如某原子附近静电势越负，只能说明亲核试剂越容易被拉到这原子附近，却不代表这原子附近的电子越活泼因而越容易参与反应。再者，当分子间接近到快要反应时，分子的实际电子分布与它在孤立时的电子分布已有很大差异，前面介绍静电势时提到的“对分子的电荷分布不产生任何影响”的前提条件已经不能满足，所以用分子孤立状态的静电势已经不能很好说明问题。
3. 综上所述，由于分子间反应涉及到了两个步骤，而两个步骤本质不同，分别适合用静电势和ALIE来分析，因此需要将它们在分子表面上的分布综合考虑以判断反应位点。通常可以先用ALIE判断，如果遇到困难，比如好几个原子都有对应的ALIE极小点（而且数值和最小点数值差不多），或者极小点介于两个原子之间，那么可以再考察静电势的分布，看看反应物更容易被拉到哪个原子附近。另外，有时需要根据分子三维结构，将位阻效应对试剂接近过程的影响考虑进去。

福井函数、双描述符、相对软度, 实际上和ALIE一样，着眼的也是反应过程的第二步.因此用它们分析活性位点时其实也应该将静电势考虑进去。此外，Hirshfeld电荷等也可以用于预测活性反应位点。
HOMO|LOMO的原子成分分析, 体系的最高占据轨道(HOMO)和亲电反应有关, 体系的最低空轨道(LUMO)和亲核反应有关. 根据HOMO或LUMO 的轨道成分多少可以判断发生亲电反应或亲核反应的位点. 分析 HOMO(LUMO)的轨道成分等同于考察这个轨道(设为第i个轨道)的电子密度在各原子附近分布的总量. 有多种方法都可以计算分子轨道中原子的成分, 例如Mulliken方法、AOIM方法、自然原子轨道法(NAOMO, 但不适合分析LUMO)
参考:《Multiwfn的定量分子表面分析功能预测反应位点等》http://blog.sciencenet.cn/blog-482864-599667.html

6. 专题:弱相互作用分析

  • 强相互作用的化学键(>200KJ/mol)>静电主导的弱相互作用(氢键、卤键)>色散主导的弱相互作用(π-π堆积、VDW)
  • 本质都是电磁力，或称为静电力。强弱相互作用，相对于kT而言(k为玻尔兹曼常数，T为温度)
  • 可以用来分析分子间作用力，大分子的机构，受体-配体的特异性识别等。
  • 所有静电为主导的弱相互作用都可以写成Lewis酸碱
  • 从类型上讲有氢键、卤键、π-π堆积、范德华力、静电作用等
    • 氢键: 与Donor和Acceptor之间的距离(<3.5A)以及角度(<30°)有关;
    • π-氢键: (π电子区显负电性，可以代替ONF作为H受体0); F…H…F 为167JK/mol。
    • 卤键: 卤族元素和C|N|O|S等之间 17-40KJ/mol
    • π-π堆积:
    • 范德华力: <10KJ/mol (除非大体系)

前面如LOL、ELF等通常是用来分析强相互作用范畴的问题。分析弱相互作用可以用：

  • 静电势分析，特别是静电势填色的范德华表面，以及范德华表面的静电势极大点、极小点位置和数值: 特别适合分子静电主导的弱相互作用。正因为如此，对静电势重现性好的原子电荷，也可以用于分析弱相互作用。
  • RDG分析(约化密度梯度, Reduced density gradient,): 可以基于波函数的RDG，也可以基于promolecular近似。凸显出体系中涉及弱相互作用的区域，以便直观地了解到分子中哪些区域与弱相互作用有关;
    • RDG=0.6等值面可以显示弱相互作用区域和原子核附近
    • 结合使用RDG函数和ρ(r)函数(电子密度)，就可以区别分子中成键区域和弱相互作用区域(左图)
      [ρ在原子核附近最大，化学键区域中等，在弱相互作用区小]
    • 结合使用RDG函数和sign(λ_2)ρ函数函数，就可以区别分子中成键区域和弱相互作用区域(右图)
    • RDG函数还可以与ELF函数联合作图:
    • 对于的大体系，可以不用真实的电子密度，而使用Promolecule密度(不涉及波函数)，再做RDG函数和sign(λ_2)ρ函数函数图，可以寻找弱相互作用。这种方式，可以直接用pdb、xyz文件
    • aRDG: averaged RDG (aRDG), 使得RDG分析方法可以用于分析多帧结构，尤为适合结合分子动力学模拟技术来研究平衡的动态环境中的弱相互作用。
  • DORI分析: Density overlap regions indicator, 只需要一个函数就可以把弱相互作用和共价键作用区域同时展现。

Multiwfn实现RDG分析:

1
2
3
4
Multiwfn: 100-1可以绘制x,y散点图和cube文件, 考察弱相互作用。还同时考察大分子局部的弱相互作用: 可以指定某些参数，使得只分析指定区域的RDG，并增加该区域内格点数，可以有更清晰的结果。
默认设置是RDG_maxrho=0.05改为0, http://sobereva.com/291
其中，vmd可视化cube文件，通过RDG函数的等值面来展现弱相互作用区域，并且将sign(λ_2) ρ函数用不同的颜色投影到RDG等值面上。
使原子坐标保持在形成分子时的状态，将自由原子密度叠加得到的密度称为Promolecule密度，可以视为在形成分子前，电子密度尚未驰豫的电子密度。构建这种密度不需要量子化学计算，只要有分子结构文件和自由原子密度就能十分容易地构建。由于在弱相互作用区域Promolecule密度和实际密度差异并不像在成键区域那么显著，Promolecule密度在此分析方法中可以近似代替实际电子密度，等值面的基本形状不会有太大差异，但是定量上，也即等值面细部特征会有一定差异，因为电子密度驰豫后会在吸引作用区域聚集，尤其是会在位阻效应区域疏散以减小斥力，位阻效应越强，改变量越大。

赏

谢谢你请我吃糖果

支付宝
微信

扫一扫，分享到微信

微信分享二维码
这不是科学的乐园<br>而是科学的森林<br>无处不充满着引人入胜的奥秘<br>Email:sunnysean@qq.com<br>