[Wien] "Plane average potential" for surface

[Oleg Rubel]

Dear Salman,

I attached a set of files used to produce a 3D slices of potential (or 
e-density) along Z-axis using LAPW5 and MATLAB.

"t-Se.in5c_tmpl" is a template input file for LAPW5.

"run_lapw5.sh" is a script that will generate a set of input files, will 
run LAPW5 using them and store the results as a set of "case.rho_.###" 
files, where ### is a slice number [1...N]. You will need to edit the 
INPUT PARAMETERS section in order to provide your geometry, case name 
and identify complex calculation.

"avrg_potential.m" is a MATLAB script that will read "case.rho_.###" 
files (again the Load data section needs your input) and shape them into 
a 3D array.

After that you can do with this array what you want.

Sorry in advance, if you find the script not very user-friendly and well 
commented. It was never meant for external users. Feel free to modify 
it, if needed.


Best regards
Oleg


P.S. I realized that email systems block shell scripts. You need to 
rename run_lapw5.hs -> run_lapw5.sh


On 06/10/2013 1:17 PM, Salman Zarrini wrote:
>
> Dear Oleg,
>
> Thank you for your answer, I know how to get case.vcoul by adjusting
> case.in0 and then lapw0 and lapw5, however what I am interested in is
> variationof coulomb potential respect to just one dimension where vacuum
> and slab are elongated(V-Z). An example of such a plot as well as my
> case.vcoul, case.rho and case.struct have been enclosed.
>
> Cheers,
>
> Salman

-- 
Oleg Rubel, PhD
Scientist, Thunder Bay Regional Research Institute
Adjunct Professor, Dept Physics, Lakehead University
290 Munro St, Thunder Bay, P7A 7T1, Ontario, Canada
Phone: +1-807-7663350
Fax: +1-807-3441948
E-mail: orubel at lakeheadu.ca
Homepage: http://www.tbrri.com/~orubel/

-------------- next part --------------
% This programm reads the potential from *.rho_* files and performs
% its averaging
%
% (c) 02.12.2008 by Oleg Rubel

clear all;

%
%   Load data
%

fprefix = 't-Se.rho_.';  % name of the WIEN2k 'case'
nZ = 200; % number of slices alonz Z axis
nX = 17;
nY = 17;
nZstep = 19; % thickness of a pice for averaging

V = zeros( nX , nY , nZ ); % allocate V

%   Read potential from different slices
for j = 0 : nZ-1
  fname = sprintf('%s%s%s', fprefix, num2str(j));
  V(:,:,j+1) = reshape( load(fname) , nX , nY ); % check! can be '...nY,nX)'
end;

mean(mean(mean(V)))

%
%   Perform averaging along Z
%

Vavg = zeros( nZ-nZstep , 1 );
z = zeros( nZ-nZstep , 1 );
for j = 1 : nZ-nZstep
  z(j) = j;
  Vsub = V(:,:,j:j+nZstep);
  Vsub2 = Vsub(:);
  Vavg(j) = mean( Vsub2 );
end;
plot(z,Vavg)
stop

-------------- next part --------------
#!/bin/sh
#
# This script is intended to build a 3D map of potential (or electron
# density) of a supercell by running 'lapw5' programm of WIEN2k
#
# (c) 21.03.2008 by Oleg Rubel

####################
# INPUT PARAMETERS #
####################

fprefix='t-Se'
fsuffix='c'
dev=10000000   # devisor
a=`echo | awk '{print 8.25}'`   # lattice constant along X [Bohr]
b=`echo $a`                         #                        Y
c=`echo | awk '{print 94.1584}'`  #                        Z
meshStep=`echo | awk '{print 0.47}'` # mesh step in Bohr


##########################
# PRINT INPUT PARAMETERS #
##########################
echo INPUT PARAMETERS:
echo fprefix = $fprefix
echo fsuffix = $fsuffix
echo dev = $dev
echo a = $a [Bohr]
echo b = $b [Bohr]
echo c = $c [Bohr]
echo meshStep = $meshStep [Bohr]


#################################
# DETERMINE NUMBER OF DIVISIONS #
#################################
echo CALCULATED VALUES:
nX=`echo $a $meshStep | awk '{printf "%i", $1/$2}'`
nY=`echo $b $meshStep | awk '{printf "%i", $1/$2}'`
nZ=`echo $c $meshStep | awk '{printf "%i", $1/$2}'`
echo nX = $nX
echo nY = $nY
echo nZ = $nZ


#############
# MAIN LOOP #
#############

for (( i = 0; i <= nZ-1; i++ ))
do
   echo i = $i of `expr $nZ - 1`
   origin1=0
   origin2=0
   origin3=`echo $nZ $dev $i | awk '{printf "%i", (1.0/$1)*$2*$3}'`
   echo "  origin =" $origin1 , $origin2 , $origin3
   xend1=`echo $nX $dev | awk '{printf "%i", (1.0-1.0/$1)*$2}'`
   xend2=0
   xend3=$origin3
   echo "  xend =" $xend1 , $xend2 , $xend3
   yend1=0
   yend2=`echo $nY $dev | awk '{printf "%i", (1.0-1.0/$1)*$2}'`
   yend3=$origin3
   echo "  yend =" $yend1 , $yend2 , $yend3

   sed "s/origin1/${origin1}/g;\
      s/origin2/${origin2}/g;\
      s/origin3/${origin3}/g;\
      s/dev/${dev}/g;\
      s/xend1/${xend1}/g;\
      s/xend2/${xend2}/g;\
      s/xend3/${xend3}/g;\
      s/yend1/${yend1}/g;\
      s/yend2/${yend2}/g;\
      s/yend3/${yend3}/g;\
      s/nX/${nX}/g;\
      s/nY/${nY}/g" < ${fprefix}.in5${fsuffix}_tmpl > ${fprefix}.in5${fsuffix}
   sleep 1
   # run 'lapw5(c)'
   echo "  run lapw5${fsuffix} lapw5.def"
   lapw5${fsuffix} lapw5.def
   sleep 3
   echo "  cp ${fprefix}.rho ${fprefix}.rho_${i}"
   cp ${fprefix}.rho ${fprefix}.rho.${i}
   echo "  cp ${fprefix}.in5${fsuffix} ${fprefix}.in5${fsuffix}.${i}"
   cp ${fprefix}.in5${fsuffix} ${fprefix}.in5${fsuffix}.${i}
   echo "  reformat < ${fprefix}.rho.${i} | sed '/^$/d' > ${fprefix}.rho_.${i}"
   reformat < ${fprefix}.rho.${i} | sed '/^$/d' > ${fprefix}.rho_.${i}
done

echo end script
exit

-------------- next part --------------
    origin1    origin2    origin3      dev
    xend1      xend2      xend3        dev
    yend1      yend2      yend3        dev
3 2 3
nX nY
RHO NO 
ATU VAL NODEBUG
NONORTHO