24小时热门版块排行榜

南方科技大学公共卫生及应急管理学院2025级博士研究生招生报考通知

返回列表

【奖励】本帖被评价11次，作者yjcmwgk增加金币 9.6 个

yjcmwgk

禁虫 (文坛精英)

[资源] 【转帖】极化连续介质模型（PCM）

（转载）极化连续介质模型（PCM）
作者：分子模拟论坛的lxq05先生

Tomasi和他的同事提出的极化连续介质模型（Polarizable Continuum Model, PCM）是一个经常用到的连续溶剂化方法，并且这些年来已经有很多的改进。PCM模型计算分子溶液中的分子自由能为三项的加和：
G(sol)=G(es)+G(dr)+G(cav)
这些成分代表了静电的（es）和散射-推斥对自由能（dr）的贡献，和空穴能（cav）。所有这三项都是由以原子位置为中心的连锁范德华球定义的空穴来计算的。反应场是通过位于分子空穴表面的点电荷（表观表面电荷模型）表示的。这里要讨论的PCM的特殊版本是一种用为Hartree-Fock的联合原子模型来构建空穴。在这种模型中，范德华表面是只由定位于重元素（非氢元素）的一些球构建的联合原子方法。每个原子的范德华半径是原子类型，连接性，分子总电荷，和连接氢原子数目的一个函数。在评估方程（1）中的三项时，这个空穴的用法稍微有些不同。(图1)
当计算空穴能Gcav时，采用范德华球定义的表面，溶剂可接近表面被用来计算散射-推斥对自由能（dr）的贡献。后一表面不同于前一表面，在后者中额外考虑了（理想化的）溶剂半径。在溶液中静电对自由能的贡献 Ges采用一个近似版本的溶剂排斥面，该排斥面通过用一个常数因子约化所有的半径，并在之后加入更多的不以原子为中心的球，以便得到稍微平滑些的表面。
定域化和表面电荷的计算是通过系统地将球表面分割成已知面积的镶嵌块计算每个表面元素的一个点电荷来达到的。
在Gaussian98中执行PCM/UAHF模型可以用SCRF关键词结合PCM专用修饰词。溶剂可以用对SCRF关键词给出Solvent= modifier来指定，可接受的溶剂名称是Water（水）, DMSO（二甲亚砜）, NitroMethane（硝基甲烷）, Methanol（甲醇）, Ethanol（乙醇）, Acetone（丙酮）, DiChloroEthane（二氯乙烷）, DiChloroMethane（二氯甲烷）, TetraHydroFuran（四氢呋喃）, Aniline（苯胺）, ChloroBenzene（氯苯）, Chloroform（氯仿）, Ether（乙醚）, Toluene（甲苯）, Benzene（苯）, CarbonTetrachloride（四氯化碳）, Cyclohexane（环己烷）, Hepaten, and Acetonitrile（乙腈）。附加的选项可以在输入文件的末尾指定并用给SCRF关键词指定Read修饰词来读入。PCM溶剂化模型可以用于HF和DFT水平下的能量计算和梯度计算。 PCM计算产生的输出可以用DUMP选项显著地扩展。
下面的例子输入解释Cs对称性的乙醇的水溶剂化自由能的单电能计算（没有几何优化,图2）：
#P B3LYP/6-31G(d) scf=tight int=finegrid SCRF=(PCM,Read,Solvent=Water)

pcm/b3lyp/6-31G(d) sp ethanol in water (Cs)

0 1
O1
C2  1  r2
C3  2  r3  1  a3
H4  3  r4  2  a4  1  180.0
H5  3  r5  2  a5  4  d5
H6  3  r5  2  a5  4  -d5
H7  2  r7  3  a7  1  d7
H8  2  r7  3  a7  1  -d7
H9  1  r9  2  a9  3  180.0

r2=1.42492915
r3=1.51965095
r4=1.09569807
r5=1.09496362
r7=1.10264669
r9=0.96904984
a3=107.81130783
a4=110.63999342
a5=110.37205263
a7=109.90077195
a9=107.87777748
d5=-120.23659087
d7=-121.12750852

DUMP

由PCM溶剂化模型引起的附加输出是由负责SCF计算的L502产生的：
-------------------------------------------------------------------
Solvent: WATER
Model  : PCM/UAHF, Icomp = 4
Version: MATRIX INVERSION
Cavity : PENTAKISDODECAHEDRA with  60 initial tesserae
-------------------------------------------------------------------
Nord Group  Hybr  Charge Alpha Radius          Bonded to
1 OH sp3 0.00 1.20  1.590    C2  [s]
2 CH2 sp3 0.00 1.20  1.860    O1  [s]  C3  [s]
3 CH3 sp3 0.00 1.20  1.950    C2  [s]
-------------------------------------------------------------------

------------------------------------------------------
         Dielectric Const  = 78.39000
         High.Fr.D.Const = 1.77600
         d(Diel.Const.)/dT = -0.35620
         Molar Volume    = 18.07000
         Therm.Exp.Coeff.  = 0.00026
         Radius          = 1.38500
         Absolute temper.  =  298.00000
         Number of spheres = 3
         OMEGA          = 40.00000
         RET             = 0.20000
         FRO             = 0.70000
         Accuracy       = 0.1D-05
------------------------------------------------------
头四行重复对水指定的设置或通常的PCM缺省设置。溶质空穴由范德华球构成，范德华球是由规则的pentakisdodecahedra（五个十二面体——六十面体？）表示的，将每个球表面分成相同尺寸的60个单元。
后面的四行列出了UAHF分析的结果，确定只有三个中心（联合原子）。对于每个中心，假设的杂化和它的形式电荷，最后半径及溶剂专用约化参数Alpha一起列出。后者通常用缺省值1.2，但是也可以用选项
ALPHA=x.x
直接指定。
输出的最后部分列出了溶剂专用参数如介电常数和有效溶剂半径，和一些更多的缺省PCM设置，如初始球的数目和参数OMEGA, RET,和FRO的当前值。后面的这三个参数控制着渐入更多球（非以原子位置为中心）的过程，以便平滑表面。新的OMEGA值可以用
OMEGA=n,n
选项设置。有意义的值在40.0-90.0之间（较高的值给出较少的加入球）。新的FRO值可以用
FRO=m，m
指定。有意义的值在0.7到0.2之间（较小的值给出较少的加入球）。RET指定加入新球的最小半径，新值可以用
RET=1.1
指定。增加的数值给出较少的附加球，非常大的值完全消除附加球。
PCM算法是首先进行一步气相能量计算，以便得到后面溶剂化自由能计算的参考点。在结束气相SCF循环后，列出了空穴产生的迭代过程的细节：
------------------------------------------------------
----------  CAVITY for ELECTROSTATIC term  -----------
------------------------------------------------------
-------  The SOLUTE is enclosed in ONE CAVITY  -------
      Total N.of Tesserae =  132
         Surface Area (Ang**2) = 97.71529
         Volume    (Ang**3) = 84.60281
------------------------------------------------------
Original Sphere On Atom Re0 Alpha    Surface
      1          O1    1.590  1.200    24.08112
      2          C2    1.860  1.200    27.26919
      3          C3    1.950  1.200    46.36498
------------------------------------------------------
----------------------------------------------------------------------
   AT CONVERGENCE
  132 Tesserae over a maximum of 1500
  Surface Area (Ang**2) = 97.71529
  Volume    (Ang**3) = 84.60281
Escaped Charge= 0.13334
Error on NUCLEAR pol.charges = 0.21898 Error on ELECTR. pol.charges =-0.33812
----------------------------------------------------------------------
------------------------------------------------------------------------------
dG(solv)/dEps       (kcal/mol) = 0.00000
------------------------------------------------------
IN VACUO Dipole moment (Debye):
X=  0.0176  Y=  1.5625  Z=  0.0000  Tot=  1.5626
IN SOLUTION Dipole moment (Debye):
X=  0.1053  Y=  1.9429  Z= -0.0019  Tot=  1.9457
------------------------------------------------------
Tessera    X       Y       Z       QTot    QSN    QSE
1    2.84136 0.54870 3.42913 0.00354  -0.17977 0.18331
.
.
.

  132 -2.70116  -2.26783  -4.13447  -0.00393  -0.25391 0.24997
在这个例子（一切正常的）中不需要附加球。总的表面用132个镶嵌块("Tesserae"

表示。由于电子波函数的长尾（事实上，没有结尾），用当前表面定义的分子体积不包含体系所有电子密度，导致一些“逃逸电荷”。在每个表面元中心有一个表面电荷“QTot”，它含有一个从溶质核电荷而来的组分和一个从溶质电子电荷而来的组分。
在解电子薛定谔方程（包括附加的反应场效应）的求解过程得到的结果给出如下：
SCF Done:  E(RB+HF-LYP) =  -155.041616090    A.U. after 10 cycles
         Convg  = 0.1222D-08          -V/T =  2.0093
         S**2 = 0.0000
KE= 1.536131750668D+02 PE=-5.252140797608D+02 EE= 1.349508863800D+02
------------------------------------------------------
-------------- VARIATIONAL  PCM  RESULTS -------------
------------------------------------------------------
      (a.u.) = -155.033805
(a.u.) = -155.040760
(a.u.) = -155.041609
      (a.u.) = -155.032883
(a.u.) = -155.041616
Total free energy in sol.
(with non electrost.terms) (a.u.) = -155.041162
------------------------------------------------------
(Unpol.Solute)-Solvent (kcal/mol) =    -4.36
(Polar.Solute)-Solvent (kcal/mol) =    -5.48
Solute Polarization (kcal/mol) =    0.58
Total Electrostatic (kcal/mol) =    -4.90
------------------------------------------------------
Cavitation energy    (kcal/mol) =    8.92
Dispersion energy    (kcal/mol) = -11.39
Repulsion energy    (kcal/mol) =    2.75
Total non electr.    (kcal/mol) =    0.29
------------------------------------------------------
DeltaG (solv)       (kcal/mol) =    -4.62
------------------------------------------------------

这里列出的：
(a.u.) = -155.033805
是非微扰的气相SCF解，它被用作后面所有步骤的参考。后面描述为：
(a.u.) = -155.040760
的能量包括非极化的溶质和非极化的溶剂之间的相互作用。和气相参考能量对比得到相应的作用能：
(Unpol.Solute)-Solvent (kcal/mol) = -4.36
报告对应于非极化的溶质和极化的溶剂相互作用的总能量之后，下面的重要信息是关于极化的溶质的总能量。
(a.u.) = -155.032883
与非极化的气相总能量的能量差列出如下：
Solute Polarization (kcal/mol) = 0.58
而且永远应该是正的。极化的溶质和极化的溶剂的全部作用体系给出总能量：
(a.u.) = -155.041616
相应的静电作用能列出如下：
Total Electrostatic (kcal/mol) = -4.90
溶剂化能的非静电部分在同一个block里给出，以空穴化和散射-推斥能的加和结束：
Total non electr. (kcal/mol) = 0.29
非静电和静电贡献的加和给出了总的溶剂化自由能：
DeltaG (solv) (kcal/mol) = -4.62
然而，需要认识到，这里说的溶剂化自由能是指温度为0 K的气相静止体系。为了在特定温度下得到热力学有意义的溶剂化自由能，这些溶剂化自由能必须用气相热化学的标准处理进行修正。对于计算溶剂化自由能自身，作为气相和溶液自由能的差值，这一步经常是省略的，PCM得出的数值被直接用来和实验值相比。对乙醇来说，实验溶剂化自由能已经被测量为-5.0kcal/mol。和这个值相比PCM预测的-4.6kcal/mol可以被认为是相当精确的。
当溶剂化能作为气相和溶液相自由能的差计算的时候，必须注意到它们各自标准态的定义。因为气相和液相浓度都是以mol值（mol/l）给出的，气相和液相的数据可以直接对比。然而，气相值经常是指分压为1atm。假设是理想气体行为，在298.15K这对应于1/24.46 mol/l。
改进的溶剂化自由能预测应该涵盖溶液中的结构松弛效应。用PCM模型的几何优化是可能的，但是比气相优化消耗多许多的时间。这不仅是由于每次的能量和提速计算需要更多的CPU时间，而且是由于优化过程的收敛缓慢和时常出现的振荡。有两个选项对于减轻一些收敛问题是有用的，它们是TSNUM和TSARE：
TSNUM指定每个球的表面元的数目。 PCM算法选择规则的多面体，该多面体的表面元的数目尽可能地接近TSNUM。除了缺省的数目60，一些更大的值64，80，或100或许对于减轻一些几何优化中的共振行为是有帮助的。
TSARE以单位（Angstrom）2指定表面元的面积。有意义的值的范围从0.4到0.2，越小的值导致越大的表面元数目。将表面元的尺寸设为特定的值导致等尺寸的表面元，而不考虑球的半径（对于TSNUM设置改变的情况并非如此）。
然而，在这两种情况中，总能量依赖于表面元的实际选择，对不同体系或不同异构体的对比只有在选项选择相同的时候才有意义。对于这里用到的乙醇的例子，几何优化在用缺省设置时在23圈仍没收敛，但是在用TSARE=0.3时在7圈内收敛，得到最后总的溶剂化能-155.043097702 au.和用气相几何作的PCM的单点能计算相比（总能量为-155.041371 au），这暗示一个-1.1kcal/mol的结构松弛能，因此对溶剂化能“改进的”的预测是-5.7kcal/mol。
在应用当前的PCM/UAHF模型于溶液中反应路径时的一个问题是为范德华半径的导出使用杂化和连结性的直接后果。因为杂化和连结性都不是平滑变化的，而是沿着反应路径突然地从一点到另一点，UAHF方法一定会在溶剂化自由能图上产生突然的断点。这些问题可以，从原理上，通过平滑地约化从一组半径到另外一组，或者通过完全地避免UAHF方法，选择与连结性无关的半径。对后者，一个通常的选择是Pauling半径，用选项RADII=Pauling可以做到。

[ Last edited by yjcmwgk on 2010-1-11 at 19:05 ]

回复此楼