1.Poisson-Boltzmann, PB模型============================
在量子化学计算和使用经典分子力学进行模拟中,常会用到隐式溶剂模型,特别是作为一种快速估算溶剂化自由能的方法,在蛋白-底物分子作用计算,和蛋白质折叠等方面特别有用。量子化学gauss中,还会经常用来估算溶剂中某些分子的pKa。那么,隐式溶剂如何做到不使用溶剂而估算溶剂的影响呢?它的原理是什么?我们首先从最基本的问题谈起,即溶剂化效应是怎么回事。溶剂化效应是发生在溶剂中的很多化学现象的关键。它影响反应的自由能势能面,改变反应速率和热力学平衡。我们知道,一旦知道一个反应的势能面,我们就知道了它的全部热力学和动力学信息。举个例子,比如乙烷的二面角的旋转,如果我们把二面角作为反应坐标,则如果我们知道了这个旋转势能面,那么我们就知道了每个二面角取值的能量(势能),那么根据Boltzmann分布,我们可以计算任意温度下,乙烷对各二面角取值的分布概率: {P_{torsion}}=exp(-E_{torsion}/kT) / P_0} {P_{torsion}}是能量为{E_{torsion}}的构象的权重,而P0是配分函数,P0是把所有构象的权重加起来。这个概率就是所谓的热力学平衡分布,即体系平衡后某构象的浓度等于总浓度乘以该构象的概率: {C_{torsion} = C_{total} * P_{torsion}} 而动力学,就是翻越构象能垒的速率,可由能垒高度ΔE得到:反应速率常数:{ k= exp^{-\\Delta E/kT}} 反应速率由阿伦尼乌斯公式,即产物浓度变化:{ d/dt = A*k*C_{reactant}} (C_reactant为反应物浓度) 由上面可知反应势能面真的很厉害。有了反应势能面,别无他求。有些初学者往往觉得跑个Gauss没什么,做MD才给力。其实,如果量子化学计算能够计算势能面,干嘛要跑MD呢?MD不过是模拟体系在那个势能面上随时间的变化而已。不过,量子化学计算个单分子在气相中的反应,搜个过渡态啥的,很准确。但是对于溶剂中的反应呢?这时候,溶剂的作用,即溶剂分子在反应物分子周围来来往往的影响就不能不考虑了。回到上面乙烷的例子,如果在水里面旋转,肯定要受到水分子的阻力什么的影响。那么多的水分子,随便哪个分子动一下,阻力就不一样,这就是熵效应,即所有水分子的不同构象对乙烷旋转势能面的影响。在考虑了周围溶剂的影响后,乙烷的势能面已经不叫势能面了,而叫旋转自由能势能面。自由能包括内能和熵两个部分。当然,熵不仅包括溶剂的,也包括乙烷自身的构象熵,比
如它自身在水里面的振动,转动,平动,都影响最后的自由能势能面结果。那么,怎么考虑溶剂化效应对反应的势能面的影响呢?即,如果我们要计算自由能势能面,那有2个方法:第一,跑MD,随机模拟溶液体系的一些构象,然后根据体系中乙烷取各构象的概率,反推出自由能势能面;第二,不跑MD,还是用量子化学的方法来计算单个乙烷的旋转势能面,但是使用隐式溶剂模型,即使用一个经验的公式,来把溶剂的影响的平均结果考虑进来,这样,直接得出自由能势能面。怎么平均化溶剂作用影响呢?学过高中物理的都知道,设置一个等效介电常数。即,真空中两电荷作用力{F= (1/4\\pi\\epsilon_0)*q_1*q_2/r^2 }
那,放到溶剂中,{F= 1/4\\pi\\epsilon0*\\epsilon_r*q_1*q_2/r^2 } 比如,水的{\\epsilon_r=78 } 当然,我们在拉格朗日图景中,要使用能量而不是牛顿力学中的力。即,我们要求的双点电荷的Coulomb势能U(r12):
{U(r_{12})=(1/4\\pi\\epsilon_0)*q_1*q_2/r } 十分简单。But,wait a minute...实际的隐式溶剂模型,并不使用Coulomb公式,而是使用Poisson-Boltzmann公式!!!为什么呢?模型不一样。Coulomb公式用的是点电荷模型,q1, q2,距离r。但是,对于真实的分子,是由原子核加核外电子云构成,电子云是有体积的。量子化学计算就是核外电子云的密度变化。即使是使用经典分子力学力场, AMBER, CHARMM什么的,也是把原子看成是带电荷的,有半径r的球的,这样,必然每个原子附近都有不同的电荷密度。如果把原子看作点电荷,进行对加和计算,则收敛慢且截断误差很大。一个更聪明,更快速的办法是:把体系分成均匀格点,计算每个格点的静电势Φ(r),然后在该点上的电荷q(r)所具有的能量则为{U(r) = q(r)*\\phi(r) } 这样,最后把各点上的U(r)加起来,就得到总的静电势能U(注意静电势与静电势能的差别)。这就是Poisson方法的思想。Poisson与Coulomb对加和计算的差别,类似于量子化学计算中DFT方法与HF的差别, PB把Coulomb的多体作用转化为随空间变化的密度问题,势能等于密度乘以静电势。现在,最关键的,是静电势如何得到?它是由Poisson解决的。简单说,Poisson给出了一个微分方程,这个微分方程描述了介电常数,电荷密度分布和静电势的关系,给出了一个由一大片电子云或几个空间分布的带电原子球ρ(r)计算静电势Φ(r)的关系。Poisson公式: {\\nabla =
\\rho(r)} ?:拉普拉斯算符,可理解为一阶偏导r:表示空间格点,你可以把空间用线
切割开,把r想象成一个个格点位置(注意跟上面的表示两电荷距离的r定义不同),格点是使用有限差分方法求解微分方程的常用方法ε:介电常数Φ:所要求的静电势ρ:各格点上的电荷密度分布这里,Φ(r)是代表体系r处的静电势,是我们要求的量, 它的二阶导由该点的电荷密度ρ(r)及其介电常数ε(r)决定。虽然这个公式的物理意义我理解不了,但是我只需要知道通过它我们可以得到静电势Φ(r)就足够了,即由上式可得: {\\nabla^2\\phi(r)=-4\\pi*\\rho(r)/\\epsilon(r)}
这样,我们得到在各点r处静电势能的二阶导?^2Φ(r),最后对整个盒子所有格点进行二重数值积分即可得到最终静电势,乘以各点上的电荷得到静电势能U(r) = q(r)*Φ(r)。MD中一般使用的计算静电势能的PME方法,即Particle Mesh Ewald方法,原理就是利用Poisson公式进行格点积分。上面只说了Poisson公式,不是说隐式溶剂用的是Poisson-Boltamann公式吗? Boltzmann干了啥?问题是这样的,上面的模型,用来计算蛋白质在一定浓度电解质溶液中的静电势能,存在一个系统误差,原因是没有考虑反离子的分布的贡献。蛋白质往往是带电荷的,结果在真实溶液体系中,在一个蛋白的不同表面处,往往分布着较高浓
度的相反电荷。这叫离子氛模型。它造成的结果是,溶液中的静电作用的衰减要比由库仑公式计算的势能面的衰减快得多。这叫离子屏蔽效应。电解质溶液理论中使用Debye长度来屏蔽效应的强弱。但是Debye长度是一个经验量,无法很方便地定义。怎么能考虑这个屏蔽效应呢?一个直接的想法是,把溶质比如蛋白附近的反离子浓度计算出来。即一个改变了的电荷分布:ρ(r)+ Δρ(r),其中后面一项Δρ(r)就是溶液体系中反离子的浓度在r处的分布。即:{\\nabla = \\rho(r)+ \\Delta\\rho(r)} 通过这个改进的Poisson公式,可以更好的考虑蛋白质在电解质溶液即离子溶液中的静电势能。当然,由于反离子分布本身改变静电势能,而静电势能又改变反离子浓度分布,因此,这必然是个自洽计算。那么怎么计算反离子浓度呢? Poisson-Boltzmann公式中的Boltzmann部分: {\\Delta\\rho(r) = n_i(r) = n_i^0*exp(-q_i*\\phi/kT)} {n_i(r)}: 溶液中离子种类i在r处的数密度(可以是分数个) {n_i^0(r)}: 溶液中离子种类i的平均数密度{q_i}:该种类离子的电荷{\\phi(r)}:盒子空间中该区域的静电势上式的物理意义是,某离子在盒子内某处数
隐式溶剂模型
