在 VASP 的差分电荷密度计算及图像处理 中介绍了差分电荷密度的计算方法与三维图像处理，在文献中常常能见到二维的差分电荷图与平均到某个方向的差分电荷曲线（如图）。一般二维的差分电荷密度图在 VESTA 中 Utilitie / 2D Data Display 便可以进行处理，本文简单介绍一下平面平均差分电荷 (The plane-average electron difference) 的数据处理方式。
image
首先使用 chgdiff.pl 对计算得到的 CHGCAR 进行差分处理，这里对脚本进行了简单的修改（输出结果位置 % 前添加空格），使得输出结果与 VASP 的 CHGCAR 的格式保持一致，脚本如下：
#!/usr/bin/env perl
#;-*- Perl -*-
@args = @ARGV;
@args == 2 || die "usage: chgdiff.pl <reference CHGCAR> <CHGCAR2>\n";
open (IN1,$args[0]) || die ("Can't open file $!");
open (IN2,$args[1]) || die ("Can't open file $!");
open (OUT,">CHGCAR_diff");
for ($i=0; $i<5; $i++) {
$line1 = <IN1>;
$line2 = <IN2>;
$header1 .= $line1;
}
# check whether it is a vasp5 format
$line1 = <IN1>;
$header1 .= $line1;
$line1 =~ s/^\s+//;
@line1 = split(/\s+/,$line1);
if ($line1[0] =~ /^\d+$/) {
@atoms1 = @line1;
} else {
$atoms1 = <IN1>;
$header1 .= $atoms1;
@atoms1 = split(/\s+/,$atoms1);
}
$line2 = <IN2>;
$line2 =~ s/^\s+//;
@line2 = split(/\s+/,$line2);
if ($line2[0] =~ /^\d+$/) {
@atoms2 = @line2;
} else {
$atoms2 = <IN2>;
@atoms2 = split(/\s+/,$atoms2);
}
$sum1 += $_ for @atoms1;
$sum2 += $_ for @atoms2;
print "Atoms in file1: ".$sum1.", Atoms in file2: ".$sum2."\n";
for ($i=0; $i<$sum1+2; $i++) {
$header1 .= <IN1>;
}
for ($i=0; $i<$sum2+2; $i++) {
$dummy = <IN2>;
}
$points1 = <IN1>;
$header1 .= $points1;
$points2 = <IN2>;
@points1 = split(/\s+/,$points1);
@points2 = split(/\s+/,$points2);
$psum1 = 1;
$psum2 = 1;
for ($i=1; $i<@points1; $i++) {
$psum1 *= $points1[$i];
$psum2 *= $points2[$i];
}
print "Points in file1: ".$psum1.", Points in file2: ".$psum2."\n";
if ($psum1 != $psum2) {die ("Number of points not same in two files!");}
print OUT $header1;
for ($i=0; $i<$psum1/5; $i++) {
$line1 = <IN1>;
$line2 = <IN2>;
@line1 = split(/\s+/,$line1);
@line2 = split(/\s+/,$line2);
for ($j=1; $j<@line1; $j++) {
$line1[$j] = $line2[$j]-$line1[$j];
}
printf OUT " %18.11E %18.11E %18.11E %18.11E %18.11E\n",$line1[1],$line1[2],$line1[3],$line1[4],$line1[5];
}
close(OUT);
close(IN2);
close(IN1);
使用命令 chgdiff.pl file1 file2 处理数据时，file2 对应总的 charge， 而 file1 则是需要减去的。得到差分电荷密度后，使用修改后的 vtotav.f 处理即可得到 The plane-average electron difference。这里有两点需要注意，
- vtotav.f 是处理 LOCPOT 文件的程序（计算功函），需要简单修改源码从而读取 CHGCAR 处理数据。
- 处理得到的数据导入绘图软件作图时，需注意，横坐标为所选取方向的格点数，纵坐标对应电荷密度乘 cell 的体积，可做相应的单位转换。
附修改后源码
PROGRAM VTOTAV
PARAMETER(NGXM=256,NOUTM=1024)
CHARACTER*80 HEADER
DIMENSION VLOCAL(NGXM*NGXM*NGXM),VAV(NOUTM)
I=0
WRITE(*,*) 'Which direction to keep? (1-3 --- 1=X,2=Y,3=Z)'
READ(*,*) IDIR
IDIR=MOD(IDIR+20,3)+1
OPEN(20,FILE='CHGCAR',STATUS='OLD',ERR=1000)
C READ(20,*,ERR=1000,END=1000) NIONS,IDUM1,IDUM2
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,'(A)',ERR=1000,END=1000) HEADER
I=0; II=0; III=0; IIII=0
READ(HEADER,*,ERR=12,END=12) I,II,III,IIII
12 NIONS=I+II+III+IIII
C READ(20,*,ERR=1000,END=1000) NIONS
READ(20,'(A)',ERR=1000,END=1000) HEADER
WRITE(*,*) NIONS
DO 10 I=1,NIONS
READ(20,*,ERR=1000,END=1000) RDUM1,RDUM2,RDUM3
10 CONTINUE
WRITE(*,*) 'positions read'
READ(20,'(A)',ERR=1000,END=1000) HEADER
READ(20,*,ERR=1000,END=1000) NGX,NGY,NGZ
NPLWV=NGX*NGY*NGZ
IF (IDIR.EQ.1) NOUT=NGX
IF (IDIR.EQ.2) NOUT=NGY
IF (IDIR.EQ.3) NOUT=NGZ
IF (NPLWV.GT.(NGXM*NGXM*NGXM)) THEN
WRITE(*,*) 'NPLWV .GT. NGXM**3 (',NPLWV,').'
STOP
ENDIF
IF (NOUT.GT.NOUTM) THEN
WRITE(*,*) 'NOUT .GT. NOUTM (',NOUT,').'
STOP
ENDIF
C READ(20,'(10F8.3)',ERR=1000,END=1000) (VLOCAL(I),I=1,NPLWV)
READ(20,*,ERR=1000,END=1000) (VLOCAL(I),I=1,NPLWV)
WRITE(*,*) 'charge density read'
CLOSE(20)
DO 20 I=1,NOUTM
20 VAV(I)=0.
SCALE=1./FLOAT(NPLWV/NOUT)
WRITE(*,*) SCALE
IF (IDIR.EQ.1) THEN
DO 150 IX=1,NGX
DO 100 IZ=1,NGZ
DO 100 IY=1,NGY
IPL=IX+((IY-1)+(IZ-1)*NGY)*NGX
VAV(IX)=VAV(IX)+VLOCAL(IPL)*SCALE
100 CONTINUE
150 CONTINUE
ELSE IF (IDIR.EQ.2) THEN
DO 250 IY=1,NGY
DO 200 IZ=1,NGZ
DO 200 IX=1,NGX
IPL=IX+((IY-1)+(IZ-1)*NGY)*NGX
VAV(IY)=VAV(IY)+VLOCAL(IPL)*SCALE
200 CONTINUE
250 CONTINUE
ELSE IF (IDIR.EQ.3) THEN
DO 350 IZ=1,NGZ
DO 300 IY=1,NGY
DO 300 IX=1,NGX
IPL=IX+((IY-1)+(IZ-1)*NGY)*NGX
VAV(IZ)=VAV(IZ)+VLOCAL(IPL)*SCALE
300 CONTINUE
350 CONTINUE
ELSE
WRITE(*,*) 'Hmmm?? Wrong IDIR ',IDIR
STOP
ENDIF
OPEN(20,FILE='PACHG')
WRITE(20,*) NOUT,IDIR
DO 500 I=1,NOUT
WRITE(20,'(I6,2X,E18.11)') I,VAV(I)
500 CONTINUE
CLOSE(20)
STOP
1000 WRITE(*,*) 'Error opening or reading file CHGCAR.'
WRITE(*,*) 'item :',I
STOP
END
使用前需编译
ifort -o vtotav.x vtotav.f
使用 ./vtotav.x 命令，选取所需方向，即可处理。
同样的道理，可以使用 vtotav.x 分别处理 CHGCAR，再在 Origin 或是 Excel 里进行列相减，输出的结果是一致的。PS: 使用王老师的 vaspkit 也可以很方便的处理平面平均电荷 (Planar-Average CHG), 再手动处理一下就可以得到 The plane-average electron difference。
附朱全喜老师的一段话
数据后处理分析不能直接就想依赖于现成的脚本，应该是读懂后再进行相应处理后才能得到自己想要的结果.
对老师们慷慨分享的 code, 首先需要感谢老师们的贡献，比如可在论文中引用相关的文章，其次，搞不懂 code 的运行原理，也要搞明白里面是怎么操作的与正确的使用方法，切不可刻舟求剑，生搬硬套。
本文相关 code 与处理方法来自网络与个人经验，感谢侯柱峰，王伟和朱全喜老师的无私分享。
图片来源：DOI 10.1039/C7TA02109G
