注：本文转载自“” (彭先森)的个人文章。
VASP的差分电荷密度
一 差分电荷密度计算
of AB:∆ρ = ρAB − ρA − ρB
NOTE:
AB, A, B 需放在相同大小的空间格子。
In of the two , the are fixed as those they have in the AB .
Some NGX NGY NGZ
结构优化后的自洽计算
mkdir AB
mkdir A
mkdir BINCAR file
SYSTEM=X
ISTART=0
ENCUT=350
EDIFF=1E-5
IBRION=-1
POTIM=0.25
NSW=0 # Only electronic-SC loops are performed
EDIFFG=-1E-2
ISMEAR=0
SIGMA=0.05
PREC=ACCURATE
ISIF=2
NPAR=4
#LWAVE=FALSE
#LCHARG=FALSE # CHGCAR is written
LREAL=Auto
#IALGO=48
ISYM=0
NGX=x # 增加
NGY=y # 增加
NGZ=z # 增加vi ~/.bashrc
alias ngx="grep -A3 'NGX' OUTCAR"
source ~/.bashrcfile
cp opt/CONTCAR POSCARA, B 各自保留优化后的 部分，编辑删除 文件
vi POSCAR
:set nu # 获取编辑行数
:x,y d # 删除 x 行至 y 行
:wq与 信息相符即可，与自洽计算时相同电荷密度，提交作业，计算完成后，获得三个 文件电荷密度，使用 导出。
二 图像处理
使用软件 VESTA 打开 AB 的
Edit - Edit Data - Data - (A ) - from data - OK
Edit - Edit Data - Data - (B ) - from data - OK
Edit - Bonds - New - 和 Do not atoms… - Max.=3 - OK
取消勾选 date 下的 Show
Click … to the
the you like and save.
在VASP 的差分电荷密度计算及图像处理中介绍了差分电荷密度的计算方法与三维图像处理，在文献中常常能见到二维的差分电荷图与平均到某个方向的差分电荷曲线（如图）。一般二维的差分电荷密度图在 VESTA 中 / 2D Data 便可以进行处理，下面继续简单介绍一下平面平均差分电荷 (The plane- ) 的数据处理方式。
首先使用 .pl 对计算得到的 进行差分处理，这里对脚本进行了简单的修改（输出结果位置 % 前添加空格），使得输出结果与 VASP 的 的格式保持一致，脚本如下：
#!/usr/bin/env perl #;-*- Perl -*- @args = @ARGV; @args == 2 || die "usage: chgdiff.pl 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 = ; $line2 = ; $header1 .= $line1; } # check whether it is a vasp5 format $line1 = ; $header1 .= $line1; $line1 =~ s/^s+//; @line1 = split(/s+/,$line1); if ($line1[0] =~ /^d+$/) { @atoms1 = @line1; } else { $atoms1 = ; $header1 .= $atoms1; @atoms1 = split(/s+/,$atoms1); } $line2 = ; $line2 =~ s/^s+//; @line2 = split(/s+/,$line2); if ($line2[0] =~ /^d+$/) { @atoms2 = @line2; } else { $atoms2 = ; @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 .= ; } for ($i=0; $i<$sum2+2; $i++) { $dummy = ; } $points1 = ; $header1 .= $points1; $points2 = ; @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 = ; $line2 = ; @line1 = split(/s+/,$line1); @line2 = split(/s+/,$line2); for ($j=1; $j<@line1; $j++) { $line1[$j] = $line2[$j]-$line1[$j]; } printf OUT " .11E .11E .11E .11E .11En",$line1[1],$line1[2],$line1[3],$line1[4],$line1[5]; } close(OUT); close(IN2); close(IN1);
使用命令.pl file1 file2处理数据时，file2 对应总的 ， 而 file1 则是需要减去的。得到差分电荷密度后，使用修改后的 .f 处理即可得到 The plane- 。这里有两点需要注意，
.f 是处理 文件的程序（计算功函），需要简单修改源码从而读取 处理数据。
处理得到的数据导入绘图软件作图时，需注意，横坐标为所选取方向的格点数，纵坐标对应电荷密度乘 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
使用./.x命令，选取所需方向，即可处理。
同样的道理，可以使用 .x 分别处理 ，再在 或是 Excel 里进行列相减，输出的结果是一致的。PS: 使用王老师的 也可以很方便的处理平面平均电荷 (- CHG), 再手动处理一下就可以得到 The plane- 。
附朱全喜老师的一段话：
数据后处理分析不能直接就想依赖于现成的脚本，应该是读懂后再进行相应处理后才能得到自己想要的结果。
对老师们慷慨分享的 code, 首先需要感谢老师们的贡献，比如可在论文中引用相关的文章，其次，搞不懂 code 的运行原理，也要搞明白里面是怎么操作的与正确的使用方法，切不可刻舟求剑，生搬硬套。
本文相关 code 与处理方法来自网络与个人经验，感谢侯柱峰，王伟和朱全喜老师的无私分享。
图片来源：DOI 10.1039/
学术之友向广大读者们征稿，内容包括但不限于：科研文献解读，科研软件使用经验分享，科研基础知识总结归纳，科研招聘广告发布，科研生活闲聊杂谈等。我们会免费为大家编辑和推广。
投稿联系：。
限时特惠:本站持续每日更新海量各大内部创业课程，一年会员仅需要98元，全站资源免费下载
点击查看详情
站长微信:Jiucxh
