实验同行做计算--电荷密度读取和处理

       在之前视频中，我们介绍了用VESTA操作显示差分电荷的内容。有朋友在评论区留言说，如果沿着一个方向去求和，怎么把电荷密度累加起来？这里给大家提供一个MATLAB的程序（讲解视频在这里）。

       当我们在INCAR里边设置LCHARG=T以后，一般就会看到两个跟电荷有关的文件，CHG和CHGCAR。后者的有效数字多一些，每行是5列，我们一般是处理这个文件。以下面分子为例，我们用VESTA做完差分的（也可以根据今天的程序自己做差分）。这里重点讲解Matlab程序的导入。可以看到，文件的前半部分是结构信息，对应POSCAR部分，这里4个碳6个氢，有晶胞参数，每个原子的具体坐标。这里要注意的是196，168，140这三个整数，对应电荷密度沿不同方向的划分，下面的数据就是每一个网格点上的电荷密度啊，这里有正有负的话呢，是因为做了差分。我们需要用程序把这些数据导成一个196×168×140的矩阵，然后我们就可以沿着任意一个方向去进行累加。这里要注意我们选的晶胞是长方体，就才可以用现在的程序（附录A）。如果比较复杂的晶胞的话，有通用的程序（@五班最强路人）。

      首先，我们要把文件前面的结构信息和三个整数那一行全删掉，保证下面是一个矩阵，这里用load的指令导入数据。这里a是921984行，5列，累计有4609920个数据，就是196x168x140，这里VASP就是按一定顺序把数据写出来的。这里的for循环就是把这些数据变成1行，4609920列，然后继续转化。这里nn=5是因为CHGCAR文件每行是5列数据。

       在程序中，不同文件，每次大家修改n2,n1,n3三个整数即可。这里有三个for循环，就是 把对应的数据放到矩阵all的相应的位置，最核心的就是这里，这里主要是根据VASP的写文件顺序确定的，z方向呢，就是140个分割点。

       生成完all矩阵后，就可以取其中的一层来画图，也可以把它加起来。变量tp的初值是第一层，从第二层一直到N3，全部把它累加，那就把电荷全加起来了，是一个二维矩阵，行列分别是168×19。可以用contour做等高线，也可以用image函数做彩图。这里用了flipud函数是为了显示和分子结构的方向一致。注意，matlab不同版本的colormap是不一样的，旧的版本image在1~60区分度很高，但其实比较新的版本是1~200多才有区分度。大家可以看到这些蓝的地方就是失电荷比较多，红的地方得电荷比较明显。在二维材料、表面等体系的电荷转移，这个显示效果就会好一点。

   附录：

clear

a=load('CHGDIFF.vasp');

sa=size(a);

b=zeros(1,sa(1)*sa(2));

nn=5;

for ii=1:length(a)

    b(1,nn*(ii-1)+1:nn*ii)=a(ii,:);

end

%  196  168  140

n2=196；n1=168；n3=140；

all=zeros(n1,n2,n3);

for ii=1:n3

    for jj=1:n1

        for kk=1:n2

            all(jj,kk,ii)=b(1,n1*n2*(ii-1)+n2*(jj-1)+kk);

        end

    end

end

tp=all(:,:,1);

for ii=2:n3

    tp=tp+all(:,:,ii);

end

contour(tp)

figure

tp=flipud(tp);

image((tp-min(min(tp)))*200/(max(max(tp))-min(min(tp))))

colormap(jet)