使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线

使用CP2K结合Multiwfn绘制密度差图、平面平均密度差曲线和电荷位移曲线

文/Sobereva@北京科音

First release: 2022-Mar-10  Last update: 2022-Mar-18

电子密度差（electron density difference, EDD）是分析成键、电子转移非常常用的手段。波函数分析程序Multiwfn（http://sobereva.com/multiwfn）在绘制和分析密度差方面非常灵活、方便、强大，笔者之前有不少文章与此相关，见《使用Multiwfn作电子密度差图》（http://sobereva.com/113）、《使用Multiwfn考察分子轨道、NBO等轨道对密度差、福井函数的贡献》（http://sobereva.com/502）、《使用Multiwfn计算激发态之间的密度差》（http://sobereva.com/429）、《使用Multiwfn做电子密度、ELF、静电势、密度差等函数的盆分析》（http://sobereva.com/179）等。之前笔者写的密度差有关的文章主要是对于孤立体系写的，本文的目的之一是介绍怎么将非常流行、免费强大的CP2K第一性原理程序与Multiwfn相结合，对周期性体系绘制密度差图。在此顺带特别强调，electron density difference的中文是密度差，“差分密度”是对密度差的严重不当称呼，绝对不要用这个词！

本文另一个目的是介绍怎么通过Multiwfn绘制平面平均密度差（plane-averaged EDD）曲线，它是密度差的衍生。密度差是个三维空间函数，经常通过绘制平面图、等值面图方式图形化考察。虽然这两种图很直观，但平面图会受到截面选取的影响，而等值面图会受到等值面数值（isovalue）选取的影响，因此都有一定任意性，也不容易定量讨论。对于电子转移有比较明确方向的情况，可以绘制平面平均密度差曲线，由此可以清楚地考察垂直于选取的方向上每个截面上的密度差积分值，在定量讨论、对比上比较便利。例如，固体表面上吸附一个小分子，吸附导致垂直于表面方向有明显的电子转移，就可以在垂直于表面的方向上绘制平面平均密度差曲线来准确、细致地考察不同截面处的电子净增、减情况。还有一种图叫做电荷位移曲线（charge displacement curve），它相当于平面平均密度差曲线的积分曲线，对于定量讨论电荷转移量情况很有帮助，本文也会介绍怎么通过Multiwfn绘制这种图。

本文使用的例子是NaCl板吸附一个CO分子，NaCl和CO将被定义为两个片段，通过密度差考察吸附导致的电子转移情况。本文文件包里的opt目录下的NaCl_CO-1.restart是CP2K在PBE-D3(BJ)/DZVP-MOLOPT-SR-GTH级别下对这个体系做优化产生的restart文件，我们将基于这个结构绘制密度差。为简化问题，在优化的时候整个NaCl板的坐标被冻结为了晶体坐标的状态，因此不涉及表面重构。

下文涉及到的相关文件可以在http://sobereva.com/attach/638/file.rar下载。cub文件由于体积太大，文件包里就不提供了。

本文的计算使用CP2K 9.1，安装方法见《CP2K第一性原理程序在CentOS中的简易安装方法》（http://sobereva.com/586）。VMD使用的是1.9.3版，可在http://www.ks.uiuc.edu/Research/vmd/免费下载。读者务必使用Multiwfn最新版，如果不知道是不是最新的就现在立刻去主页http://sobereva.com/multiwfn下载一个。如果对Multiwfn不了解，看《Multiwfn FAQ》（http://sobereva.com/452）和《Multiwfn入门tips》（http://sobereva.com/167）。不了解cub文件的话可以看看《Gaussian型cube文件简介及读、写方法和简单应用》（http://sobereva.com/125）。

这一节我们要得到整个体系、NaCl板和CO的各自的电子密度的cub文件，然后由Multiwfn求差。

我们先创建一个计算整个体系的电子密度cub文件的CP2K输入文件。启动Multiwfn，然后输入 NaCl_CO-1.restart  //在本文文件包里 cp2k  //创建CP2K输入文件。在《使用Multiwfn非常便利地创建CP2K程序的输入文件》（http://sobereva.com/587）一文里有详述 NaCl_CO.inp  //产生的输入文件的路径

现在进入了创建CP2K输入文件的界面。为了之后绘制的平面平均密度差曲线比较好观看，最好把当前体系挪到盒子中央附近。先看一眼当前的结构是什么样，输入 -11  //进入几何操作界面。此界面的详细介绍见《Multiwfn中非常实用的几何操作和坐标变换功能介绍》（http://sobereva.com/610） 0  //观看当前结构

现在出现了图形窗口。点击菜单上的Other settings - Toggle showing cell frame把盒子边框显示出来，再选Other settings - Toggle showing all boundary atoms把边界原子显示出来。现在看到的图如下

可见当前整个体系在盒子底部，NaCl有三层，由于最下层的原子正好在Z=0的位置，因此这些边界原子在盒子顶端也显示了出来。点击右上角的RETURN按钮返回。为了把体系挪到盒子中央去，输入 24  //平移体系使得选中的部分在晶胞中居中 1-110  //当前体系里所有原子。也可以直接按回车选中整个体系 0  //再次观看体系 可见此时体系确实在盒子中央了

关闭图形窗口，接着输入 -10  //返回CP2K输入文件创建界面 -7  //设置周期性 XY  //XY方向二维周期性 -3  //要求产生cube文件 1  //电子密度 0  //产生CP2K输入文件

现在当前目录下就有了NaCl_CO.inp。读者可根据要求修改里面的CUTOFF等设置，本例不那么讲究，就用默认的。运行之，算完后当前目录下就有了NaCl_CO-ELECTRON_DENSITY-1_0.cube，将之改名为NaCl_CO.cub。

将NaCl_CO.inp复制为NaCl.inp，用文本编辑器打开，把其中PROJECT NaCl_CO改为PROJECT NaCl，把&COORD里面C和O原子删掉。用CP2K运行NaCl.inp，将得到的cub文件改名为NaCl.cub。

将NaCl_CO.inp复制为CO.inp，用文本编辑器打开，把PROJECT NaCl_CO改为PROJECT CO，把&COORD里面Na和Cl原子删掉。用CP2K运行CO.inp，将得到的cub文件改名为CO.cub。

现在用Multiwfn对三个cub文件求差。启动Multiwfn，然后输入 NaCl_CO.cub  //输入此文件实际路径 13  //处理格点数据 11  //格点数据运算 4  //减去另一个格点数据 NaCl.cub  //输入此文件实际路径 11  //格点数据运算 4  //减去另一个格点数据 CO.cub  //输入此文件实际路径

现在内存里的格点数据就已经是密度差格点数据了。想观看一下等值面的话，可以进入选项-2 Visualize isosurface of present grid data。由于物理吸附对应的密度差的数量级往往很小，所以要用比较小的等值面数值才能清楚看到。在图形窗口右侧把等值面数值设为0.0005后，看到的图像如下（点save picture按钮保存的图像比直接在图形窗口里看到的更好）

绿色和蓝色等值面分别代表格点数据数值为正和为负的部分。当前密度差是按照ρ(NaCl...CO)-ρ(NaCl)-ρ(CO)算的，显然绿色体现吸附后电子密度增加的部分，蓝色体现电子密度减少的部分。点击右上角RETURN关闭图形窗口。

若想得到密度差的cub文件，就选0 Export present grid data to Gaussian-type cube file (.cub)，输入EDD.cub，然后内存里的格点数据就被导出为了当前目录下的EDD.cub。

使用《在VMD里将cube文件瞬间绘制成效果极佳的等值面图的方法》（http://sobereva.com/483）里提供的我写的脚本，在VMD程序里只需要输入一行命令cub EDD 0.0004，就可以把EDD.cub的0.0004等值面显示出来（如果没显示出来，一个字一个字看http://sobereva.com/483这篇文章）。在文本窗口里输入pbc box把盒子显示出来，再改成正交视角，经过Tachyon渲染看到的就是下面的图。图中蓝色小球是Na，棕色是Cl，等值面的蓝色和绿色和前面Multiwfn里的定义相同。

从此图可以明显看出，CO的碳结合到Na上后，C和O上的电子密度整体大幅减少，很大程度往Na的方向上转移了。这也容易理解，当前体系里Na带显著的正电荷，CO结合上去，特别是C还带着丰富的孤对电子，必然其电子会被往Na那边显著地极化。

再来看看把等值面数值设得更小的情况。在VMD文本窗口里输入cubiso 0.0001把等值面数值调为0.0001，稍微调整视图后看到的图如下

可见此时的密度差图展现了比之前更多的信息。粉色虚线圈住的地方电子密度有一定降低，由于仅在等值面数值很小的时候才能显现出来，所以这部分降低量很少（想通过积分得到确切值的话需要用Multiwfn的盆分析功能，前面http://sobereva.com/179文中介绍过）。这部分电子密度降低是因为CO向Na这边转移了电子，由于交换互斥作用，使得原本这部分的电子被排斥开，向更远处转移。

这一节绘制密度差局部积分曲线和平面平均密度差曲线，首先了解一下定义。某个三维空间函数的局部积分曲线（local integral curve）可以对任意方向绘制，如果对比如Z方向绘制，则表达式为

其中p是被考察的函数。如果p取为密度差的话，那么这个曲线就是前述的密度差局部积分曲线了。

还有个概念叫积分曲线（integral curve），就是对局部积分曲线再进一步积分。例如对Z方向如下所示

其中积分的起点z_ini由自己决定。如果被考察的函数是密度差，则这个积分曲线对应的就是前述的电荷位移曲线。

平面平均曲线（plane-averaged curve）定义如下，其中A_XY是格点数据对应的盒子在XY方向的面积。这种曲线体现的是被考察的函数在特定截面上的平均值。

在Multiwfn里可以基于格点数据对任意函数绘制上述三种曲线，需要先把相应的格点数据存到内存里。比如可以通过主功能5格点数据运算功能算出来后自动存到内存里，也可以启动时从cub等格点数据文件里直接载入，等等。

下面绘制密度差局部积分曲线。如果你之前已经关了Multiwfn，就再启动Multiwfn，载入密度差格点数据文件EDD.cub，然后进入主功能13，再进入子功能18，这是用来绘制（局部）积分曲线和平面平均曲线的功能。如果你按照上一节操作之后还没关Multiwfn，则内存里的格点数据正是密度差格点数据，因此直接进入子功能18就行了。

在子功能18里，Multiwfn先问你对哪个方向绘制，我们想垂直于NaCl界面绘制，因此输入Z，因为NaCl界面平行于XY方向。之后Multiwfn问你绘制范围的下限和上限位置是什么，此例想对整个Z范围进行绘制，因此按照屏幕上的提示直接输入a。之后会看到一个菜单，通过相应选项可以直接绘制积分曲线、局部积分曲线、平面平均曲线，可以保存它们的图像，还可以把曲线数据导出为文本文件便于用Origin、gnuplot、qtgrace等第三方程序绘图。

现在选择屏幕上的选项2 Plot graph of local integral curve，由于当前内存里的格点数据是密度差，所以下图就是密度差局部积分曲线

读者如果想让等值面图和密度差局部积分曲线便于准确对照，可以用屏幕上相应选项让Multiwfn导出局部积分曲线数据为当前目录下的locintcurve.txt（每一列的含义看导出时屏幕上的提示），然后放到Origin等程序里作图。之后在VMD里用正交视角，把盒子边框显示出来，把保存的等值面图和局部积分曲线图放到一起，把前者的白色背景设为透明色或者抠掉，再调节位置和大小让两张图边框位置正好对上，如下所示。不会用Photoshop的话用版本不很老的PowerPoint也可以完成

再改一下透明度，裁切一下边缘去掉盒子边框，最终得到下面的图，效果很好。由这个图可以比较准确地捕捉各个截面上电子的净增减，可见在C与表层的NaCl相接的地方电子密度由于相互作用而增加明显（但绝对的数量级并不算大）。

如果想要绘制平面平均密度差曲线，前面选2 Plot graph of local integral curve的地方改成选3 Plot graph of plane-averaged curve即可，其它不变。 

在Multiwfn当前界面里选1 Plot graph of integral curve，由于当前内存里的格点数据是密度差，因此得到的积分曲线的就是电荷位移曲线了：

上图横轴左端对应于NaCl...CO体系Z=0的位置，横轴右端对应于NaCl...CO体系Z=25埃（盒子Z尺寸）的位置。这个图是对上一节的密度差局部积分曲线从左向右积分产生的，从肉眼可以清楚看出来对应关系。曲线最右边数值为0，体现了密度差全空间积分为0，也即CO与NaCl结合没有导致总电子数有净变化。

为了便于观看，将电荷位移曲线和结构图+等值面图叠加起来。我也把曲线最高的峰的横、纵坐标位置标注上了。电荷位移曲线纵坐标是无量纲的，因为对应的是电子数。

这张图体现出，在Z