分子的静电势(Molecular electrostatic potential，MEP)广泛应用于计算化学以及化学反应性的理论中[1]，比如利用静电势分子解释氢键、卤键、硫键和磷键等由静电驱动的分子相互作用[2，3]；预测分子的化学反应位点[4]；预测分子的溶剂化自由能[5]、pKb值[6]、分子间相互作用能[7]等。主流的量子化学程序几乎都支持静电势的计算，但由于静电势的计算标度是O(N2)，故大体系静电势的计算将相当耗时。因此提出了一种密度拟合分子静电势(Density fitting molecular electrostatic potential，DF-MEP)的方法，用以快速计算大体系的静电势，相关文章发表于Journal of Computational Chemistry，2023，44(7): 806-813，密度的简单计算，链接为。，
1. 密度拟合近似
本部分我们将给出DF-MEP的推导过程。分子体系在空间中某一点r'的静电势为：
其中ZA为原子核A的核电荷数，RA是原子核的核坐标，ρ(r')为在位置r'的电子密度，其形式为：
其中D是密度矩阵，μ(r)和ν(r)为基函数，Nbas是基函数的数量。将(2)式代入到(1)式中可得：
上式中Naux为辅助基函数数目，C为展开系数，将(4)式代入到电子密度的表达式(2)可得：
(5)式中的dk被称为密度拟合系数，通过(5)式可以将电子密度的计算标度由二次方降低到一次方，而其中最关键的是求解出展开系数，为此我们定义一个残差函数：
以及残差函数的模：
这里的(μν|μν)，(μν|k)和(k|k´)分别是四中心、三中心和双中心双电子积分，将(7)式进行极小化可得：
在(8)式两端乘以密度矩阵元Dμν后求和可以得到密度拟合系数dk所满足的方程：
求得dk后，将(5)式代入到(2)式可获得密度拟合近似下的静电势：
2. 密度拟合静电势的精度和速度
2.1 精度测试
选取H2O···H2O、EtCOON、FeFC以及Valine共4个分子测试DF-MEP的精度与速度。图1给出了4个分子在ρ=0.001 e/Bohr3的电子密度等值面上的静电势。图1中每个体系左侧的图是Conv-MEP方法获得的静电势图，右侧则是使用DF-MEP方法获得的静电势图计算的，计算级别皆为B3LYP/def2-SVP，其中辅助基组为def2/J。从图1可以观察到，两种方法绘制的静电势图几乎是一致的，说明DF-MEP方法能够获得与Conv-MEP方法效果一致的静电势图。
图1. 四个体系在ρ=0.001 e/Bohr3的等值面上的静电势图(kcal/mol)
2.2 速度测试
我们选取了丙氨酸链(Ala)n测试Conv-MEP和DF-MEP随着n的增加的耗时的情况，相关结果列于表1中。从表1中可知：DF-MEP比Conv-MEP的速度快一个数量级以上。随着体系的增大，Conv-MEP与DF-MEP的耗时比率也越来越大，也就是说体系越大，DF-MEP的优势越明显。需要特别说明的是，密度拟合本身也需要一定的计算量，但通过表1可知密度拟合过程对DF-MEP的耗时增加不大，且密度拟合只需要计算一次，另外与SCF相比密度拟合的时间可忽略不计。
1、定义：单位体积某种物质的质量。2、公式：密度＝质量/体积 ρ＝m/v 3、单位：克/厘米3 或者 千克/米3 换算关系：1克/厘米3＝1000千克/米3 公式：ρ＝m/v中ρ、m、v都是对同一物体而言。对于同一物质，ρ一定。
另外我们还对比了蛋白质Trp-Cage在DF-MEP与Conv-MEP方法下的计算时间与静电势图，其中计算的格点数目为699200，相关图展示于图2。从图2中可知：两种方法的静电势图几乎是一致的，在16核并行的情况下，Conv-MEP的耗时为2060.2 s，而DF-MEP的耗时仅为16.6 s。最后选取1011个原子的Amylose分子，使用DF-MEP计算并绘制了图2所示的静电势图，在16核并行的情况下，计算1011个原子的体系静电势所需的时间仅为287.60 s。因此，DF-MEP方法不仅能够获得与Conv-MEP方法一致的静电势图，还能大大地缩短计算大体系静电势所需的时间，这也为更大体系静电势的获得提供了可能。
图2. Trp-Cage与Amylose在ρ=0.001 e/Bohr3的等值面上的静电势图(kcal/mol)
3. 密度拟合静电势的计算
利用pfdmdf和aesp两个程序可以获得密度拟合下的静电势，程序可到
下载。pfdmdf是将Gaussian计算得到的fchk文件进行密度拟合得到密度拟合系数，aesp是使用pfdmdf得到的密度拟合文件计算密度拟合的静电势。
pdfmdf的使用方式为：准备好Gaussian获得的fchk文件，比如名称为test.fchk，则使用命令：
计算公式：ρ=m/V 密度等于物体的质量除以体积，可以用符号ρ（读作rou）表示，国际单位制和中国法定计量单位中，密度的单位为千克每立方米，符号是kg/m3。
pfdmdf test.fchk def2-J -aesp 4
获得密度拟合文件test.df，文件中存储着辅助基函数和密度拟合系数。上述命令中的def2-J为def2/J辅助基组，4为并行核数，用户可以根据机器情况进行修改。使用aesp需要准备一个输入文件，比如test.in，具体为：
密度的计算公式：式中V为包含P点的体积元；M为该体积元的质量。在厘米·克·秒制中，密度的单位为克/厘米3。在国际单位制和中国法定计量单位中，密度的单位为千克/米3。
file test.df
task esp
unit angs
read type2
0.20 4.0
0.15 3.0
0.10 2.0
file字段后面为test.df文件的名称；task字段后面为任务类型，这里指的是计算静电势；unit字段后面为长度单位，这里指的是使用Å为单位；read字段后面为读取格点的类型，上面的例子的含义为：分子在x方向两端延拓4.0 Å，每个点的间隔为0.20 Å；分子在y方向两端延拓3.0 Å，每个点的间隔为0.15 Å；分子在z方向两端延拓2.0 Å，每个点的间隔为0.10 Å。将test.df文件与test.in放在同一目录下，然后输入命令：
aesp test.in
4. 结论
使用密度拟合静电势的DF-MEP方法，可以获得与传统静电势方法一致的静电势图，相同的精度，但计算的速度却可提升两个数量级以上。利用DF-MEP可以轻松地实现上千个原子的静电势计算，体系越大则DF-MEP方法计算静电势越具有优势。
[1] K. S. Murray，K. Sen，Molecular electrostatic potentials: concepts and applications，1996，Elsevier.
[2] P. Politzer，J. S. Murray，T. Clark，Phys. Chem. Chem. Phys. 2010，12，7748-7757.
[3] G. Cavallo，P. Metrangolo，T. Pilati，G. Resnati，G. Terraneo，Cryst. Growth Des. 2014，14，2697-2702.
密度的计算公式：密度＝质量÷体积，所以要计算密度，需要先知道质量和体积。
[4] C. H. Suresh，G. S. Remya，P. K. Anjalikrishna，WIREs Comput. Mol. Sci. 2022，12，e1601.
[7] N. Mohan，C. H. Suresh，J. Phys. Chem. A 2014，118，1697-1705.
