#!/bin/bash #Creator: Stephan N. Steinmann #Help can be obtained from: # stephan.steinmann@ens-lyon.fr # or form # carine.michel@ens-lyon.fr #Purpose: Compute the electronic Free Energy as a function of the electrochemical potential #Requirements: vasp 5.x Output and LOCPOT files in folders that are named according to their NELECT tag. # workfunction.x utility, to be obtained from # https://github.com/WMD-Bath/Workfunction/blob/master/workfunction_v2.f # Or provided by this package # Gnuplot, awk #How to use it: #a) make sure you have at least 3 converged charge states, including the neutral system and the corresponding LOCPOT files. #b) you need to know the index of the point central point in "empty space" in z-direction (if you slab is symmetric around z=0.0, it will be close to the middle of the cell, i.e. NGZF/2). #c) Compute the factor z0/Z: Z= out-of-plane lattice vecotr of the unit cell; z0=z-direction not occupied by the metal slab. #d) Run the script in the "main" folder, i.e., in the folder where the subfolders with the names of the NELECT values are located. # Use as: ./Filhol-Neurock.sh NELECT_Neutral INDEX_central_point VALUE_of_z0_divided_by_Z #Output: #a) if no workfunction.x is found in your directory, it will be created, assuming that you have either ifort or gfortran in your path. #b) dat file: Data used for fitting in gnuplot. #c) fit.plot Inputfile for gnuplot: useful if you have to delete data points that deviate from the parabola because of charge-spilling into the background #d) fit.log: log file for the fit, containing the paramerters for the parabola #e) dat.dat: values of the parabola in a regular interval, given in the fit.plot file #f) charge.dat: the values of the charge as a function of the electrochemical potential #g) pot.png: image showing the fitted G(V) parabola and the corresponding points. #h) charge.png: same, but the electronic charge as a function of the electrochemical potential #Info: # We strongly recommend to use at least an implicit solvent model, e.g., # http://theory.mse.cornell.edu/software/vaspsol.html #Setup the necessary files: echo "3" > 3 #Gnuplot fitting file: change the interval and resolution to match your needs. exec 3> fit.plot cat >&3 <<-EOF set term png set xlabel "Potential(vs SHE)/V" set ylabel "G/eV" set output "pot.png" f(x)=a*x*x+b*x+c f2(x)=a2*x*x+b2*x+c2 fit f(x) 'dat' u 1:3 via a,b,c fit f2(x) 'dat' u 1:4 via a2,b2,c2 plot 'dat' u 1:3 ti "data", f(x) ti "fit" set output "charge.png" plot 'dat' u 1:4 ti "data", f2(x) ti "fit" #Change for resolution: set samples 12 set table "dat.dat" set format y "%+-14.6f" #Change for interval of your interest plot [-2.000:0.7500] f(x) set table "charge.dat" set format y "%+-14.6f" set samples 15 plot [-2.000:1.5000] f2(x) print f(-2) EOF exec 3>&- if [ -z `command -v workfunction.x` ] ; then # # Close the file descriptor for the workfunction file. # exec 3> workfunction.x cat >&3 <<-EOF C C Adapted from the program VTOTAV that was C distributed with VASP 4 and is of uncertain origin. C C V1 - Original distro version C V2 - Various updates from JLF Da Silva and A Walsh (2007-2009) C http://people.bath.ac.uk/aw558 (a.walsh@bath.ac.uk) C C Reads in the 3D electrostatic potential (LOCPOT) from VASP C which is generated with LVHAR= .TRUE. in the INCAR file. C C Outputs the 1D potential averaged along one of the lattice C vectors, which can be used, e.g., to compute the vacuum level. C C Note: to obtain the Hartree + Exc contributions, the flag C LVTOT = .TRUE. should be set instead. In VASP 4, this flag (confusingly) C provided the Hatree contributions only. C C To do: C 1. Automatically align values using the calculated vacuum level. C 2. Check if dipole corrections are required (slope of potential in the vacuum). C 3. Change axis units to Angstrom (for all IDIR). C 4. Update to read artbitrary number of ion types. C PROGRAM VTOTAV PARAMETER(NGXM=256,NOUTM=1024) CHARACTER*80 HEADER DIMENSION NATOM(10) DIMENSION VLOCAL(NGXM*NGXM*NGXM),VAV(NOUTM) WRITE(*,*) '========================' WRITE(*,*) 'Workfunction V2 (2009)' WRITE(*,*) '========================' WRITE(*,*) '' WRITE(*,*) 'Which direction to average along? (1=x;2=y;3=z)' READ(*,*) IDIR IDIR=MOD(IDIR+20,3)+1 C Read in LOCPOT OPEN(20,FILE='LOCPOT',STATUS='OLD',ERR=1000) 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 NIONS=0 Do i=1,10 NATOM(i)=0 ENDDO READ(HEADER,*,ERR=12,END=12) (NATOM(i),i=1,10) 12 DO i=1,10 NIONS=NIONS+NATOM(i) ENDDO READ(20,'(A)',ERR=1000,END=1000) HEADER WRITE(*,*) 'Number of ions:',NIONS DO 10 I=1,NIONS READ(20,*,ERR=1000,END=1000) RDUM1,RDUM2,RDUM3 10 CONTINUE WRITE(*,*) '(i) 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 READ(20,*,ERR=1000,END=1000) (VLOCAL(I),I=1,NPLWV) WRITE(*,*) '(ii) 3D potential read...' CLOSE(20) C Average 3D Potential VAVT=0.0 DO 20 I=1,NOUTM 20 VAV(I)=0. SCALE=1./FLOAT(NPLWV/NOUT) C 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 write(*,*) IX,VAVT,VAV(IX) VAVT=VAVT+VAV(IX)/NGX 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 VAVT=VAVT+VAV(IY)/NGY 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 VAVT=VAVT+VAV(IZ)/NGZ 350 CONTINUE ELSE WRITE(*,*) 'Wrong value of IDIR',IDIR STOP ENDIF Z0 = 0.d0 OPEN(20,FILE='vplanar.txt') WRITE(20,*) NOUT,IDIR C WRITE(20,*) NOUT,IDIR, C3 DO 500 I=1,NOUT WRITE(20,'(I6,2X,E18.11)') I,VAV(I) C WRITE(20,'(I6,2X,F10.6,E18.11)') I,(C3/(NOUT-1))*(I-1),VAV(I) 500 CONTINUE WRITE(*,*) '' WRITE(*,*) 'Done. Always check for convergence!' C VAVT=VAVT/SCALE write(*,*) 'Average potential ', VAVT write(20,*) '#Average potential ', VAVT CLOSE(20) STOP 1000 WRITE(*,*) 'Error opening or reading file LOCPOT!' STOP END EOF # # Close the file descriptor for the workfunction file. # exec 3>&- if [ ! -z `command -v ifort` ] ; then ifort workfunction.f -o workfunction.x elif [ ! -z `command -v gfortran` ] ; then gfortran workfunction.f -o workfunction.x else echo "No Fortran compiler identified" fi else echo "Found workfunction" fi #Define the function to extract useful data from OUTCAR and vplanar.txt files: get_el(){ line=$2 if [ -e $1/OUTCAR ] ; then echo ${1%%/} ` awk '/y=/{E=$7}/E-fermi/{F=$3} END{print E, F}' $1/OUTCAR` `awk -v n=$line '{if($1==n){E=$2}if($0~/Aver/){A=$3}}END{printf "%10.5f\n", (E-A)*1.0}' $1/vplanar.txt` elif [ -e $1/OUTCAR.gz ] ; then echo ${1%%/} `zgrep 'E-fermi\|y=' $1/OUTCAR.gz | awk '/y=/{E=$7}/E-fermi/{F=$3} END{print E, F}'` `awk -v n=$line '{if($1==n){E=$2}if($0~/Aver/){A=$3}}END{printf "%10.5f\n", (E-A)*1.0}' $1/vplanar.txt` else echo "No OUTCAR or OUTCAR.gz file" fi } #Store the arguments from the command line: ref=$1 point=$2 z0overZ=$3 #Create the 1D electrostatic potentials using the workfunction utility: for i in `ls */ -d` ; do if [ -e $i/LOCPOT ] ; then cd $i workfunction.x < ../3 #Zip the LOCPOT files, after all, they are quite big bzip2 LOCPOT cd .. fi done #Process data: for i in `ls */ -d` ; do if [ -e $i/vplanar.txt ] ; then get_el $i $point ; fi ; done | sort -k 1,1 > raw awk -v z=$z0overZ -v a=$ref '{ N[NR]= $1; E[NR]=$2; EF[NR]=$3; EV[NR]=$4; W[NR]=$4-$3; if($1==a){ref=NR}} END{Int=0; #Reference point print W[ref]-4.44,E[ref],E[ref],a #Positively charged for(i=ref-1;i>=1;i--){ if(iref){ Int-=(N[i-1]-N[i])*(EV[i]+EV[i-1])/2; print W[i]-4.44,E[ref]+z*(E[i]-E[ref]+Int),E[ref]+z*(E[i]-E[ref]-Int+W[i]*(N[i]-N[ref])),N[ref]+z*(N[i]-a),N[i],W[i]*(N[i]-N[ref])*z}} }' raw > dat gnuplot < fit.plot 2> fit.log #Clean up: rm raw 3