分子的静电势(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)式中可得：

从(3)式可以看出，计算某个点处的静电势标度为O(N2)，因此大体系静电势计算的时间会随着体系尺寸的增加而快速增长。利用密度拟合近似，可以降低大体系分子静电势的计算时间和标度。在密度拟合近似中，(2)式中的基函数乘积μ(r)ν(r)可以使用一组辅助基函数ηk(r)展开：

上式中Naux为辅助基函数数目，C为展开系数，将(4)式代入到电子密度的表达式(2)可得：

(5)式中的dk被称为密度拟合系数，通过(5)式可以将电子密度的计算标度由二次方降低到一次方，而其中最关键的是求解出展开系数，为此我们定义一个残差函数：

以及残差函数的模：

这里的(μν|μν)，(μν|k)和(k|k´)分别是四中心、三中心和双中心双电子积分，将(7)式进行极小化可得：

在(8)式两端乘以密度矩阵元Dμν后求和可以得到密度拟合系数dk所满足的方程：

求得dk后，将(5)式代入到(2)式可获得密度拟合近似下的静电势：

我们将(10)式的计算静电势的方式称为DF-MEP，其计算标度为线性的，相比(3)式的Conv-MEP，DF-MEP可以大幅度缩短大体系静电势的计算时间。

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. DF-MEP(t1)，Conv-MEP(t2)和密度拟合(t3)随着丙氨酸链(Ala)n增大的耗时(s)，以及t2/t1的比率

另外我们还对比了蛋白质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，则使用命令：

pfdmdf test.fchk def2-J -aesp 4

获得密度拟合文件test.df，文件中存储着辅助基函数和密度拟合系数。上述命令中的def2-J为def2/J辅助基组，4为并行核数，用户可以根据机器情况进行修改。使用aesp需要准备一个输入文件，比如test.in，具体为：

file test.dftask espunit angsread type20.20 4.00.15 3.00.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

可获得密度拟合后静电势test.cub文件。更多关于aesp的功能和使用方法可参考链接中的手册。

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.

[5] J. S. Murray, F. Abu-Awwad, P. Politzer, J. Phys. Chem. A 1999, 103, 1853-1856.

[6] J. Sandoval-Lira, G. Mondragon-Solorzano, L. I. Lugo-Fuentes, J. Barroso-Flores, J. Chem. Inf. Model. 2020, 60, 1445-1452.

[7] N. Mohan, C. H. Suresh, J. Phys. Chem. A 2014, 118, 1697-1705.