使用Multiwfn做基于电子密度差描述符的电荷转移分析

sobereva

使用Multiwfn做基于电子密度差描述符的电荷转移分析

Using Multiwfn to analyze charge transfer based on electron density difference descriptors

文/Sobereva@北京科音   2026-Aug-11

1 前言

电子密度差（electron density difference, EDD。差分电荷密度之类都是其非常恶心的谬称）是量子化学和第一性原理领域被广为使用的考察电子转移情况的三维实空间函数，详细介绍见《谈谈量子化学研究中电子转移情况怎么分析考察》（http://sobereva.com/777）。EDD通常是作图观看，这往往不太容易定量讨论，而且等值面往往比较复杂、不好图形描述。J. Chem. Theory Comput., 7, 2498 (2011)中提出了基于EDD的一种分析思想，对EDD定义了一系列描述符，从而明显更便于分析讨论和横向对比EDD展现的信息，并且还将EDD的正值、负值部分分别以Gaussian函数平滑化，使之主体分布更便于考察。此方法的普适性很强，只要能算EDD就可以用这种分析方法。Multiwfn很早以前就支持了这个分析方法，并且我在原文的基础上对此方法做了诸多改进，使之明显更为普适、好用。本文就介绍一下此方法，第2节介绍其原理，第3节示例使用Multiwfn做这种分析考察基态的分子间电子转移，第4节示例用这种方法分析电子激发导致的电子转移。

本文使用Multiwfn 2026.7.15版，Gaussian 16 C.02版。强烈建议不了解Multiwfn者阅读《Multiwfn FAQ》（http://sobereva.com/452）、《Multiwfn入门tips》（http://sobereva.com/167）、《详谈Multiwfn支持的输入文件类型、产生方法以及相互转换》（http://sobereva.com/379）。Multiwfn可以在官网http://sobereva.com/multiwfn免费下载。使用Multiwfn做本文介绍的分析发文章时，应按照程序使用条款明确引用Multiwfn启动时提示的程序原文。

特别值得一提的是，Multiwfn中独家的空穴-电子分析方法在考察电子激发特征方面早已非常流行，这在《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（http://sobereva.com/434）中做了非常详尽的介绍。如果你用的是研究激发态问题最常用的TDDFT（spin-flip TDDFT、sTD-DFT/sTDA形式也支持），那么在考察电子激发导致的电子分布变化时没任何必要用本文介绍的功能。对于分析电子激发，本文介绍的方法的额外主要好处是支持更多的激发态计算方法，比如EOM-CCSD、delta-SCF等，只要你用的计算程序能产生Multiwfn支持的记录了基态和激发态波函数的文件，或者能产生它们的电子密度格点数据，就可以用Multiwfn的这个功能分析对应的EDD。


2 方法原理

这个方法定义了很多描述符，在Multiwfn手册3.21.3节有具体介绍并且有完整的计算公式，这一节只是用文字简要描述一下，细节请看手册。

此方法把EDD的正值和负值部分分别称为ρ+和ρ-函数。ρ+的全空间积分值，或者ρ-的全空间积分的绝对值，被称为q_CT，这体现了有多少电子的分布实际发生了改变。

ρ+和ρ-的质心坐标可以计算出来，分别记作R+和R-。R+和R-之间的距离称为D指数，可以衡量电子转移距离，其笛卡尔分量Dx、Dy、Dz可以衡量不同方向的转移距离。注意用这个衡量转移距离只适合单方向的电子转移，而如果是比如完全对称的电子转移，D指数会精确为0、无法说明任何问题。

由R+和R-的坐标还可以获得电子转移导致的偶极矩的变化。例如Y+和Y-分别是R+和R-的Y坐标，电子转移导致的偶极矩Y分量的变化为Y+减Y-再乘以q_CT。

对ρ+和ρ-可以分别计算衡量它们空间分布广度的σ+和σ-。例如ρ-在Z方向的分布广度可以通过σ(-,Z)表示，它对应于这个方向上ρ-函数分布的RMSD（根均方偏差）。Δσ指数定义为|σ+|-|σ-|，衡量ρ+和ρ-整体分布广度的差异程度。

还可以把ρ+和ρ-分别转化为三维高斯函数C+和C-描述，中心分别在R+和R-，往四周衰减的速度分别由σ+和σ-决定。C+的全空间积分值和ρ+相同，C-的全空间积分值和ρ-相同。C+和C-的图像特征比ρ+和ρ-明显更为简单、更易于肉眼考察，但失去了EDD的细节。

进一步还可以计算一些量，由于牵扯到一些公式，我这里就直接用量子化学波函数分析与Multiwfn程序培训班（http://www.keinsci.com/WFN）全面系统讲解各种电子激发分析一节里的幻灯片直接说明了：

如果大家看过《使用Multiwfn做空穴-电子分析全面考察电子激发特征》（http://sobereva.com/434）一文的话，会发现上面这些描述符和Multiwfn的空穴-电子分析里的非常类似，但有以下不同：
(1)空穴-电子分析只能用于考察电子激发导致的电子分布变化，而本文介绍的方法可以用于考察任意因素导致电子分布变化的情况，对于《谈谈量子化学研究中电子转移情况怎么分析考察》（http://sobereva.com/777）里提到的各种机制导致的电子转移原理上都可以用（至于实际效果怎么样是另一回事）。
(2)空穴-电子分析里不是基于ρ+和ρ-计算那些描述符，而是基于空穴和电子分布计算的。对于研究电子激发问题，ρ+和ρ-体现的是电子净增加和减少的部分，无法把被激发的电子离开的区域和去的区域的重叠的部分描述出来（即存在相互抵消因此损失了内在信息），而空穴和电子的分布则可以完整全面展现电子是从哪里离开的，又去了哪里，且能展现出重叠区域。
(3)空穴-电子分析对应的是非弛豫密度，而本文的分析方法可以基于弛豫也可以基于非弛豫密度。弛豫和非弛豫密度的差别在http://sobereva.com/434里专门做过说明。用非弛豫密度原理上更加准确，但需要花很多额外耗时去计算，而且一次只能算一个激发态的密度，从这点来说明显不如空穴-电子分析用着方便。

由于上述原因，如果你考察的是电子激发造成的电子密度变化，除非用的是空穴-电子分析不支持的激发态计算方法，否则不建议用本文的方法。而且空穴-电子分析框架下还能做广阔得多的额外分析，如考察原子/基函数对空穴/电子的贡献量、片段间电子转移（IFCT）定理分析（http://sobereva.com/433）、电荷转移光谱CTS（http://sobereva.com/628）、电荷转移矩阵（http://sobereva.com/436）等，没看过《Multiwfn支持的电子激发分析方法一览》（http://sobereva.com/437）的话建议看看，这些分析都提到了。


3 基于密度差分析基态分子间电子转移实例：TFF-TCNE复合物

四硫富瓦烯（tetrathiafulvalene, TTF）是个很强的电子给体，四氰乙烯（tetracyanoethylene, TCNE）是很强的电子受体，二者形成的复合物不需要通过吸收光子，光是在基态下就可以出现很显著的TTF到TCNE的显著的电子转移，这是个很经典的有机电荷转移复合物。这一节通过基于EDD的描述符对其电子转移情况进行考察，使用wB97XD/6-31+G*级别。

这一节涉及的文件都在http://sobereva.com/attach/776/file.7z里的TTF-TCNE目录下。其中optfreq.gjf是对TFF-TCNE复合物做优化和振动分析的输入文件，计算后产生的optfreq.chk文件用formchk转换后得到的optfreq.fchk提供在文件包里了，其中记录了复合物的极小点结构下的波函数信息。把优化完的复合物中的TTF和TCNE部分分别直接抠出来创建单点任务的输入文件TTF.gjf和TCNE.gjf，注意一定要写上nosymm以避免Gaussian将片段自动摆到标准朝向下使得坐标与整体不符而令EDD没有意义，这一点在《使用Multiwfn作电子密度差图》（http://sobereva.com/113）中明确强调了。计算完毕后得到了记录两个片段的波函数文件TTF.fchk和TCNE.fchk，这俩文件也在文件包里。

启动Multiwfn，载入optfreq.fchk，然后输入
5  //计算格点数据
0  //自定义运算
2  //有两个波函数文件将与当前波函数文件运算
-,TTF.fchk  //假定此文件在当前目录下，所以不用写目录名，后同
-,TCNE.fchk
1  //电子密度
3  //高质量格点（计算约1728000个点）

马上就算完了。从计算过程中屏幕上提示的信息可以看到格点间距为0.202641 Bohr。格点间距越小，之后基于EDD计算的描述符越准确，因为前述的那些描述符在Multiwfn里是基于均匀格点积分的。当前这样0.2 Bohr的格点间距就足够精细了。对于尺寸明显更大的体系，“高质量格点”对应的格点间距也可能也偏大，此时在设定格点的界面中建议选择4 Input the number of points or grid spacing in X,Y,Z, covering whole system并且自己直接把格点间距设为0.2 Bohr。关于Multiwfn中格点设定的更多信息看《Multiwfn FAQ》（http://sobereva.com/452）里的Q39。

现在进入了后处理菜单，直接选-1看一下EDD的原本的等值面是什么样。把等值面数值调到0.001后看到的图像如下所示，绿色和蓝色分别是正值和负值等值面，显然分别对应形成复合物后电子密度增加和减少的区域，明显可看到电子是整体从下方（TTF）往上方（TCNE）转移。

点Multiwfn图形窗口右上角的RETURN按钮关闭之。现在EDD的格点数据已经在内存中了，之后可以直接做下面说的基于EDD的描述符计算。如果你想把EDD格点数据导出成cube格式的文件的话就选择2，这样以后再次分析EDD时就不需要像上面这样重新算一遍了，直接在Multiwfn启动后载入EDD的cube文件即可。

接着输入
0  //返回主菜单
18  //电子激发分析
3  //基于EDD格点数据分析电子转移

立马看到以下结果，各个部分和第2节介绍的各种指数的对应关系不言自明。

q_CT (positive and negative parts):   0.356  -0.356 a.u.
Barycenter of positive part in x,y,z (Angstrom):   0.680   1.060   0.408
Barycenter of negative part in x,y,z (Angstrom):  -0.110  -0.733  -0.323
Distance of CT in x,y,z (Angstrom):   0.790   1.793   0.731  D index:   2.091
Dipole moment variation (a.u.) :  -0.532  -1.208  -0.492 Norm:   1.408
Dipole moment variation (Debye):  -1.352  -3.069  -1.251 Norm:   3.580
RMSD of positive part in x,y,z (Angstrom):  1.998  1.789  1.560 Total:   3.103
RMSD of negative part in x,y,z (Angstrom):  2.415  1.786  1.529 Total:   3.370
Difference between RMSD of positive and negative parts (Angstrom):
X:  -0.417  Y:   0.003  Z:   0.031  delta_sigma index:  -0.267
H_x:  2.206  H_y:  1.788  H_z:  1.544  H_CT:  1.826  H index:  3.236 Angstrom
t index:   0.265 Angstrom
Overlap integral between C+ and C- (i.e. S+- index):  0.837311

可见，形成复合物导致分布发生了变化的电子的数目，即q_CT，是0.356，对于有机分子间的基态电子转移来说这已经是很大的数值了。这里特别要强调的是绝对不要简单地把q_CT直接当做分子间的电子转移量，因为在TTF和TCNE区域内EDD都是同时有正有负、ρ+和ρ-和并非分别完全分布在TCNE和TTF区域中，显然不能把分子内的电子密度的分布调整也算到电子转移量里。若想得到真正的电子转移量，应该算片段电荷来考察，见《谈谈量子化学研究中电子转移情况怎么分析考察》（http://sobereva.com/777）。顺带一提，以上输出中的0.356和-0.356分别是Multiwfn对ρ+和ρ-的积分值，当格点定义得合理时（盒子足够大、格点间距足够小），二者的绝对值是相同的，如果绝对值明显不同就暗示给出的各种描述符都可能不太准确，需要调整格点定义。

从以上输出还可见，当前D指数是2.091埃，算是电子转移的距离较明显。电子转移导致偶极矩变化了3.580 Debye。ρ+和ρ-的分布广度相仿佛，RMSD都是三点几埃，其平均值H指数也是如此。S+-指数不小、t指数只是略正，都说明ρ+和ρ-的分布不算很充分分离开。

在后处理菜单中选择1 Show isosurface of C+ and C- functions simultaneously，在图形窗口中把等值面数值设为0.0003，然后在顶端的菜单中选Isosurface style - Use transparent face，就会看到下图。绿色椭圆和蓝色椭圆分别对应C+和C-，即分别把ρ+和ρ-用Gaussian函数平滑化后的样子。可见C+和C-虽然分别主要分布在TCNE和TTF上，但在俩分子之间重叠得还是很显著的，确实EDD的正、负部分分离得不充分。

也可以在关闭图形窗口后选择2把C+和C-的格点数据分别导出成当前目录下的Cpos.cub和Cneg.cub文件，以便在VMD里绘制效果更完美的等值面图。


4 基于密度差分析电子激发实例：DMABN的荧光发射

4-(N,N-二甲基氨基)苯腈（DMABN）在极性足够大的溶剂下具有双荧光特性，低波长的S1→S0发射是局域激发，并且分子整体是准平面的。高波长的S1→S0发射属于twisted intramolecular charge transfer (TICT)激发，其中二甲氨基相对于其它部分呈近乎90度垂直。这个例子我们用CAM-B3LYP/def-TZVP通过TDDFT方式产生其TICT垂直发射时的激发态密度，用CAM-B3LYP/def-TZVP产生相应结构下的基态密度，然后用Multiwfn计算EDD，并计算相关描述符。基态和激发态计算都用IEFPCM模型描述乙腈溶剂环境。这里假定读者已经有了北京科音初级量子化学培训班（http://www.keinsci.com/KEQC）或《Gaussian中用TDDFT计算激发态和吸收、荧光、磷光光谱的方法》（http://sobereva.com/314）中介绍的TDDFT计算常识。

这一节涉及的文件都在http://sobereva.com/attach/776/file.7z里的DMABN目录下。S1_opt.gjf是优化S1态极小点的输入文件，初猜结构就很接近TICT结构，因此收敛到了相应的极小点。这个输入文件关键词为#p CAM-B3LYP/TZVP TD opt scrf(solvent=acetonitrile) out=wfn，并且在坐标后面空一行处写了S1.wfn，因此任务完成后当前目录下就出现了S1.wfn，其中记录了优化出来的对应TICT的S1极小点结构下的S1态的弛豫密度的自然轨道（用《在Multiwfn中基于fch产生自然轨道的方法与激发态波函数、自旋自然轨道分析实例》http://sobereva.com/403里介绍的方法基于此任务的fch文件获得记录自然轨道的mwfn文件代替S1.wfn亦可，只不过需要多做一步）。S0.gjf是在这个极小点结构下做DFT基态计算的输入文件，关键词为#p CAM-B3LYP/TZVP scrf(solvent=acetonitrile)，并且写了%chk=S0.chk，因此算完后用formchk把S0.chk转成S0.fchk，此文件里就包含了此结构下的基态波函数信息（若用out=wfn把基态波函数信息导出成wfn文件用于后面的EDD计算亦可）。

启动Multiwfn，然后输入以下命令计算EDD格点数据并绘图。注意是S0减S1，因为是荧光垂直发射，如果是S1减S0那就成了研究垂直吸收了
S0.fchk
5  //计算格点数据
0  //自定义运算
1  //有一个波函数文件将与当前波函数文件运算
-,S1.wfn
1  //电子密度
3  //高质量格点
-1  //观看等值面图

把等值面数值改为0.005，看到的EDD等值面如下所示。绿色和蓝色分别对应正值和负值部分，分别对应荧光垂直发射导致电子密度增加和减少的部分。可以看出电子主要是从氰基向二甲氨基方向转移，但是正负等值面交错混杂在一起不太好看，特别是在二甲氨基与苯环连接的地方

关闭图形窗口，接着在Multiwfn里输入
0  //返回主菜单
18  //电子激发分析
3  //基于EDD格点数据分析电子转移

马上看到以下统计信息。可见这个DMABN的TICT荧光发射造成了电子分布的十分显著的变化，q_CT达到1.042。正值和负值部分的质心有明显的分离，D指数达到2.419埃。荧光发射过程造成了体系偶极矩变化高达12.1 Debye。由于如之前的EDD图所看到的，电子密度增加和减少部分交错分布显著，因此t指数虽然明显为正但不算特别正，S+-也不小

q_CT (positive and negative parts):   1.042  -1.042 a.u.
Barycenter of positive part in x,y,z (Angstrom):  -0.098  -1.978  -0.000
Barycenter of negative part in x,y,z (Angstrom):   0.091   0.433   0.000
Distance of CT in x,y,z (Angstrom):   0.189   2.412   0.000  D index:   2.419
Dipole moment variation (a.u.) :   0.372   4.748   0.000 Norm:   4.762
Dipole moment variation (Debye):   0.945  12.068   0.000 Norm:  12.105
RMSD of positive part in x,y,z (Angstrom):  1.248  1.617  0.938 Total:   2.248
RMSD of negative part in x,y,z (Angstrom):  0.903  2.294  1.132 Total:   2.713
Difference between RMSD of positive and negative parts (Angstrom):
X:   0.345  Y:  -0.677  Z:  -0.194  delta_sigma index:  -0.465
H_x:  1.076  H_y:  1.956  H_z:  1.035  H_CT:  1.952  H index:  2.481 Angstrom
t index:   0.467 Angstrom
Overlap integral between C+ and C- (i.e. S+- index):  0.776950

之后选1观看C+和C-的等值面图，把等值面数值设为0.003后看到的图如下所示。这个图比EDD图看起来清楚直观不少，可见电子密度增加部分的中心差不多是二甲氨基的氮的位置，电子密度减少部分的中心差不多是在苯环中央。

之后选2可以把C+和C-分别导出成Cpos.cub和Cneg.cub，载入VMD 1.9.3（用其它版本后果自负），在Graphics - Representation界面里把二者的Drawing method都设为Isosurface，对Cpos.cub和Cneg.cub的等值面数值分别设0.003和-0.003，Material都设为Transparent，Coloring method选为ColorID并对Cpos.cub和Cneg.cub分别设为绿色和蓝色。然后在VMD文本窗口里输入以下命令把C+和C-的中心位置画成小球
draw color yellow
draw sphere { -0.098  -1.978  -0.000 } radius 0.2 resolution 15
draw sphere {  0.091   0.433   0.000 } radius 0.2 resolution 15
之后再用powerpoint之类标注箭头。现在看到的图如下所示，可见把荧光垂直发射造成的密度分布变化展现得很易于理解。

5 总结&其它

本文介绍了Multiwfn具有的基于EDD考察电子增加、减少部分的分布特征以及衡量电子转移情况描述符的功能，还介绍了能把EDD的正值、负值部分原本复杂的分布转化为简单直观的C+和C-函数的功能。这个功能非常普适，只要有EDD格点数据就能做。EDD格点数据的来源可以是：
(1)像本文这样由Multiwfn基于波函数文件直接算出来。
(2)你提供不同状态的电子密度格点数据，用Multiwfn主功能13的子功能11做格点数据求差运算得到EDD格点数据，示例见Multiwfn手册4.13.2节，《使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线》（http://sobereva.com/638）文中也示例了用Multiwfn对CP2K产生的不同电子密度格点数据之间求差得到EDD的过程。
(3)任意程序产生的EDD格点数据。
顺带一提，Multiwfn支持的格点数据格式非常丰富，绝不仅限于cube格式，全面介绍见Multiwfn手册2.5节。比如VA$P的记录电子密度的CHGCAR文件也可以直接被Multiwfn载入，可以由Multiwfn的主功能13的子功能0导出成cube格式，之后再用Multiwfn求差产生EDD格点数据。

最后再次提醒，使用Multiwfn做本文介绍的分析发文章时，请恰当引用Multiwfn启动时提示的程序原文。