The plane-average electron difference in VASP

chempeng / 2017-11-29

在VASP的差分电荷密度计算及图像处理 VASP 的差分电荷密度计算及图像处理中介绍了差分电荷密度的计算方法与三维图像处理，在文献中常常能见到二维的差分电荷图与平均到某个方向的差分电荷曲线（如图）。一般二维的差分电荷密度图在 VESTA 中 Utilitie / 2D Data Display 便可以进行处理，本文简单介绍一下平面平均差分电荷 (The plane-average electron difference) 的数据处理方式。

首先使用 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