1#!/bin/sh
2
3# run from directory where this script is
4cd `echo $0 | sed 's/\(.*\)\/.*/\1/'` # extract pathname
5EXAMPLE_DIR=`pwd`
6
7# check whether echo has the -e option
8if test "`echo -e`" = "-e" ; then ECHO=echo ; else ECHO="echo -e" ; fi
9
10$ECHO
11$ECHO "$EXAMPLE_DIR : starting"
12$ECHO
13$ECHO "This example shows a calculation of STM maps."
14
15# set the needed environment variables
16. ../../../environment_variables
17
18# required executables and pseudopotentials
19BIN_LIST="pw.x pp.x plotrho.x projwfc.x sumpdos.x"
20PSEUDO_LIST="Al.pz-vbc.UPF As.pz-bhs.UPF"
21
22$ECHO
23$ECHO "  executables directory: $BIN_DIR"
24$ECHO "  pseudo directory:      $PSEUDO_DIR"
25$ECHO "  temporary directory:   $TMP_DIR"
26$ECHO "  checking that needed directories and files exist...\c"
27
28# check for directories
29for DIR in "$BIN_DIR" "$PSEUDO_DIR" ; do
30    if test ! -d $DIR ; then
31        $ECHO
32        $ECHO "ERROR: $DIR not existent or not a directory"
33        $ECHO "Aborting"
34        exit 1
35    fi
36done
37for DIR in "$TMP_DIR" "$EXAMPLE_DIR/results" ; do
38    if test ! -d $DIR ; then
39        mkdir $DIR
40    fi
41done
42cd $EXAMPLE_DIR/results
43
44# check for executables
45for FILE in $BIN_LIST ; do
46    if test ! -x $BIN_DIR/$FILE ; then
47        $ECHO
48        $ECHO "ERROR: $BIN_DIR/$FILE not existent or not executable"
49        $ECHO "Aborting"
50        exit 1
51    fi
52done
53
54# check for gnuplot
55GP_COMMAND=`which gnuplot 2>/dev/null`
56if [ "$GP_COMMAND" = "" ]; then
57        $ECHO
58        $ECHO "gnuplot not in PATH"
59        $ECHO "Results will not be plotted"
60fi
61
62# check for pseudopotentials
63for FILE in $PSEUDO_LIST ; do
64    if test ! -r $PSEUDO_DIR/$FILE ; then
65       $ECHO
66       $ECHO "Downloading $FILE to $PSEUDO_DIR...\c"
67            $WGET $PSEUDO_DIR/$FILE $NETWORK_PSEUDO/$FILE 2> /dev/null
68    fi
69    if test $? != 0; then
70        $ECHO
71        $ECHO "ERROR: $PSEUDO_DIR/$FILE not existent or not readable"
72        $ECHO "Aborting"
73        exit 1
74    fi
75done
76$ECHO " done"
77
78# how to run executables
79PW_COMMAND="$PARA_PREFIX $BIN_DIR/pw.x $PARA_POSTFIX"
80PP_COMMAND="$PARA_PREFIX $BIN_DIR/pp.x $PARA_POSTFIX"
81PLOTRHO_COMMAND="$BIN_DIR/plotrho.x"
82PROJWFC_COMMAND="$PARA_PREFIX $BIN_DIR/projwfc.x $PARA_POSTFIX"
83SUMPDOS_COMMAND="$BIN_DIR/sumpdos.x"
84$ECHO
85$ECHO "  running pw.x as:      $PW_COMMAND"
86$ECHO "  running pp.x as:      $PP_COMMAND"
87$ECHO "  running plotrho.x as: $PLOTRHO_COMMAND"
88$ECHO "  running projwfc.x as: $PROJWFC_COMMAND"
89$ECHO "  running sumpdos.x as: $SUMPDOS_COMMAND"
90$ECHO "  running gnuplot as:   $GP_COMMAND"
91$ECHO
92
93# self-consistent calculation
94cat > AlAs110re.scf.in << EOF
95 &control
96    calculation = 'scf'
97    restart_mode='from_scratch',
98    title='AlAs 110 surface slab, relaxed (central plane fixed)'
99    pseudo_dir = '$PSEUDO_DIR/',
100    outdir='$TMP_DIR/',
101    prefix='AlAs110'
102 /
103 &system
104    ibrav=  8, celldm(1) =7.424621202, celldm(2)=1.414213576,
105    celldm(3)= 6.00000,
106    nat= 14, ntyp= 2,
107    ecutwfc =14.0,
108 /
109 &electrons
110    mixing_mode = 'plain'
111    mixing_beta = 0.7
112    conv_thr =  1.0d-6
113 /
114 ATOMIC_SPECIES
115 Al  26.98  Al.pz-vbc.UPF
116 As  74.92  As.pz-bhs.UPF
117ATOMIC_POSITIONS alat
118 As  0.000000000  -0.044777195  -0.058158722
119 Al  0.500000000   0.251460333   0.113525468
120 As  0.500000000   0.712279009   0.504183885
121 Al  0.000000000   1.067633546   0.480460620
122 As  0.000000000  -0.003937059   0.995826731
123 Al  0.500000000   0.351400965   1.004220212
124 As  0.000000000  -0.003937059   2.004173269
125 Al  0.500000000   0.351400965   1.995779788
126 As  0.500000000   0.712279009   2.495816115
127 Al  0.000000000   1.067633546   2.519539380
128 As  0.000000000  -0.044777195   3.058158722
129 Al  0.500000000   0.251460333   2.886474532
130 As  0.500000000   0.707106800   1.500000000
131 Al  0.000000000   1.060660200   1.500000000
132K_POINTS {automatic}
133 6 2 1 0 0 0
134EOF
135$ECHO "  running the scf calculation...\c"
136$PW_COMMAND < AlAs110re.scf.in > AlAs110re.scf.out
137check_failure $?
138$ECHO " done"
139
140cat > AlAs110re.nonscf.in << EOF
141 &control
142    calculation = 'nscf'
143    restart_mode='from_scratch',
144    title='AlAs 110 surface slab, relaxed (central plane fixed)'
145    pseudo_dir = '$PSEUDO_DIR/',
146    outdir='$TMP_DIR/',
147    prefix='AlAs110'
148 /
149 &system
150    ibrav=  8, celldm(1) =7.424621202, celldm(2)=1.414213576,
151    celldm(3)= 6.00000,
152    nat= 14, ntyp= 2,
153    occupations='smearing', smearing='gaussian', degauss=0.01, nbnd=34,
154    ecutwfc =14.0,
155 /
156 &electrons
157    mixing_mode = 'plain'
158    mixing_beta = 0.7
159    conv_thr =  1.0d-6
160 /
161 ATOMIC_SPECIES
162 Al  26.98  Al.pz-vbc.UPF
163 As  74.92  As.pz-bhs.UPF
164ATOMIC_POSITIONS alat
165 As  0.000000000  -0.044777195  -0.058158722
166 Al  0.500000000   0.251460333   0.113525468
167 As  0.500000000   0.712279009   0.504183885
168 Al  0.000000000   1.067633546   0.480460620
169 As  0.000000000  -0.003937059   0.995826731
170 Al  0.500000000   0.351400965   1.004220212
171 As  0.000000000  -0.003937059   2.004173269
172 Al  0.500000000   0.351400965   1.995779788
173 As  0.500000000   0.712279009   2.495816115
174 Al  0.000000000   1.067633546   2.519539380
175 As  0.000000000  -0.044777195   3.058158722
176 Al  0.500000000   0.251460333   2.886474532
177 As  0.500000000   0.707106800   1.500000000
178 Al  0.000000000   1.060660200   1.500000000
179K_POINTS {automatic}
180 12 4 1 0 0 0
181EOF
182$ECHO "  running the non-scf calculation...\c"
183$PW_COMMAND < AlAs110re.nonscf.in > AlAs110re.nonscf.out
184check_failure $?
185$ECHO " done"
186
187# post-processing for stm images (sample bias given in Ry!)
188cat > AlAs110.pp_stm-.in << EOF
189 &inputpp
190    prefix  = 'AlAs110'
191    outdir='$TMP_DIR/',
192    filplot = 'AlAsresm-1.0'
193    sample_bias=-0.0735d0,
194    plot_num= 5
195 /
196 &plot
197   nfile=1
198   filepp(1)='AlAsresm-1.0'
199   weight(1)=1.0
200   iflag=2
201   output_format=2
202   e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
203   e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
204   x0(1)=0.0, x0(2)=-0.18,   x0(3)=3.25
205   nx=36 ,ny=56
206   fileout='AlAs110-1.0'
207 /
208EOF
209$ECHO
210$ECHO "  running the post-processing phase, negative bias...\c"
211$PP_COMMAND < AlAs110.pp_stm-.in > AlAs110.pp_stm-.out
212check_failure $?
213$ECHO " done"
214
215# run plotrho to do the figure
216cat > AlAs110.plotrho-.in << EOF
217AlAs110-1.0
218AlAs110-1.0eV.ps
219n
2200.00005 0.0078  8
221EOF
222$ECHO "  running plotrho on negative bias data...\c"
223$PLOTRHO_COMMAND < AlAs110.plotrho-.in > AlAs110.plotrho-.out
224check_failure $?
225$ECHO " done"
226
227# post-processing for stm images (negative bias, constant current)
228cat > AlAs110.pp_isostm-.in << EOF
229 &inputpp
230 /
231 &plot
232    nfile=1
233    filepp(1)='AlAsresm-1.0'
234    weight(1)=1.0
235    iflag=2
236    output_format=7    
237    fileout='AlAs110.pp_isostm-.dat'
238    e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
239    e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
240    nx=150, ny=150
241    isostm_flag=.true.
242    isovalue=0.00005
243    heightmin=0.4
244    heightmax=0.75
245    direction=1
246 /
247EOF
248$ECHO
249$ECHO "  STM image, negative bias and constant current...\c"
250$PP_COMMAND < AlAs110.pp_isostm-.in > AlAs110.pp_isostm-.out
251check_failure $?
252$ECHO " done"
253
254# run gnuplot to do the figure
255if [ "$GP_COMMAND" = "" ]; then
256    break
257else
258cat > gnuplot.tmp <<EOF
259set term postscript enhanced color solid lw 3 24
260set output 'AlAs110-1.0eV.isoplot.ps'
261set xlabel "x (bohr)"
262set ylabel "y (bohr)"
263set pm3d map
264set size ratio -1
265set palette rgb 21,22,23
266set tics out
267unset key
268splot [0:51.972][0:52.500] 'AlAs110.pp_isostm-.dat'
269EOF
270$ECHO
271$ECHO "  plotting results ...\c"
272$GP_COMMAND < gnuplot.tmp
273$ECHO " done"
274rm gnuplot.tmp
275fi
276
277# post-processing for stm images (as before, but positive bias)
278cat > AlAs110.pp_stm+.in << EOF
279 &inputpp
280    prefix  = 'AlAs110'
281    outdir='$TMP_DIR/',
282    filplot = 'AlAsresm+1.0'
283    sample_bias=0.0735d0,
284    plot_num= 5
285 /
286 &plot
287   nfile=1
288   filepp(1)='AlAsresm+1.0'
289   weight(1)=1.0
290   iflag=2
291   output_format=2
292   e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
293   e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
294   x0(1)=0.0, x0(2)=-0.18,   x0(3)=3.25
295   nx=36 ,ny=56
296   fileout='AlAs110+1.0'
297 /
298EOF
299$ECHO "  running the post-processing phase, positive bias...\c"
300$PP_COMMAND < AlAs110.pp_stm+.in > AlAs110.pp_stm+.out
301check_failure $?
302$ECHO " done"
303
304# plotrho
305cat > AlAs110.plotrho+.in << EOF
306AlAs110+1.0
307AlAs110+1.0eV.ps
308n
3090.00002 0.0021  8
310EOF
311$ECHO "  running plotrho on positive bias data...\c"
312$PLOTRHO_COMMAND < AlAs110.plotrho+.in > AlAs110.plotrho+.out
313check_failure $?
314$ECHO " done"
315
316# post-processing for stm images (positive bias, constant current)
317cat > AlAs110.pp_isostm+.in << EOF
318 &inputpp
319 /
320 &plot
321    nfile=1
322    filepp(1)='AlAsresm+1.0'
323    weight(1)=1.0
324    iflag=2
325    output_format=7    
326    fileout='AlAs110.pp_isostm+.dat'
327    e1(1)=7.0, e1(2)=0.0,     e1(3)=0.0
328    e2(1)=0.0, e2(2)=7.07107, e2(3)=0.0
329    nx=150, ny=150
330    isostm_flag=.true.
331    isovalue=0.00005
332    heightmin=0.4
333    heightmax=0.75
334    direction=1
335 /
336EOF
337$ECHO
338$ECHO "  STM image, positive bias and constant current...\c"
339$PP_COMMAND < AlAs110.pp_isostm+.in > AlAs110.pp_isostm+.out
340check_failure $?
341$ECHO " done"
342
343# run gnuplot to do the figure
344if [ "$GP_COMMAND" = "" ]; then
345    break
346else
347cat > gnuplot.tmp <<EOF
348set term postscript enhanced color solid lw 3 24
349set output 'AlAs110+1.0eV.isoplot.ps'
350set xlabel "x (bohr)"
351set ylabel "y (bohr)"
352set pm3d map
353set size ratio -1
354set palette rgb 21,22,23
355set tics out
356unset key
357splot [0:51.972][0:52.500] 'AlAs110.pp_isostm+.dat'
358EOF
359$ECHO
360$ECHO "  plotting results ...\c"
361$GP_COMMAND < gnuplot.tmp
362$ECHO " done"
363rm gnuplot.tmp
364fi
365
366# Projection of the DOS on volumes (boxes)
367cat > AlAs110.box.projwfc.in << EOF
368 &projwfc
369    prefix  = 'AlAs110'
370    outdir='$TMP_DIR/',
371    ngauss=0
372    degauss=0.01
373    DeltaE=0.02
374    tdosinboxes=.true.
375    plotboxes=.true.
376    n_proj_boxes=8
377
378!! Boxes centered on the first vacuum layer:
379  !! 1) above the surface Al
380    irmin(1,1)= 0, irmax(1,1)= 2, irmin(2,1)= 0, irmax(2,1)= 2, irmin(3,1)=63, irmax(3,1)=65,
381  !! 2) above the surface As
382    irmin(1,2)= 9, irmax(1,2)=11, irmin(2,2)= 5, irmax(2,2)= 7, irmin(3,2)=63, irmax(3,2)=65,
383  !! 3) above the 2nd layer Al
384    irmin(1,3)= 9, irmax(1,3)=11, irmin(2,3)=14, irmax(2,3)=16, irmin(3,3)=63, irmax(3,3)=65,
385  !! 4) as large as the surface unit cell
386    irmin(1,4)= 1, irmax(1,4)=18, irmin(2,4)= 1, irmax(2,4)=27, irmin(3,4)=63, irmax(3,4)=65,
387
388!! Same as above, centered on the second vacuum layer:
389    irmin(1,5)= 0, irmax(1,5)= 2, irmin(2,5)= 0, irmax(2,5)= 2, irmin(3,5)=72, irmax(3,5)=74,
390    irmin(1,6)= 9, irmax(1,6)=11, irmin(2,6)= 5, irmax(2,6)= 7, irmin(3,6)=72, irmax(3,6)=74,
391    irmin(1,7)= 9, irmax(1,7)=11, irmin(2,7)=14, irmax(2,7)=16, irmin(3,7)=72, irmax(3,7)=74,
392    irmin(1,8)= 1, irmax(1,8)=18, irmin(2,8)= 1, irmax(2,8)=27, irmin(3,8)=72, irmax(3,8)=74,
393 /
394EOF
395$ECHO
396$ECHO "  running local DOS calculation...\c"
397$PROJWFC_COMMAND < AlAs110.box.projwfc.in > AlAs110.box.projwfc.out
398check_failure $?
399$ECHO " done"
400
401# Projection of the DOS on atomic wavefunctions
402cat > AlAs110.projwfc.in << EOF
403 &projwfc
404    prefix  = 'AlAs110'
405    outdir='$TMP_DIR/',
406    ngauss=0
407    degauss=0.01
408    DeltaE=0.02
409    tdosinboxes=.false.
410 /
411EOF
412$ECHO
413$ECHO "  running projected DOS calculation...\c"
414$PROJWFC_COMMAND < AlAs110.projwfc.in > AlAs110.projwfc.out
415check_failure $?
416$ECHO " done"
417
418$ECHO
419$ECHO "  summing the atomic PDOS...\c"
420$SUMPDOS_COMMAND AlAs110.pdos_atm\#10\(Al\)_wfc* > "AlAs110.pdos_atm#10(Al)" 2> /dev/null
421$SUMPDOS_COMMAND AlAs110.pdos_atm\#11\(As\)_wfc* > "AlAs110.pdos_atm#11(As)" 2> /dev/null
422$ECHO " done"
423
424#
425#  if gnuplot was found, the results are plotted
426#
427if [ "$GP_COMMAND" = "" ]; then
428    break
429else
430eFermi=`grep "Fermi" AlAs110re.nonscf.out | cut -d  \  -f 14`
431cat > gnuplot.tmp <<EOF
432set term postscript enhanced color solid lw 3 24
433set output 'AlAs110.box.projwfc.ps'
434ef=$eFermi
435set xlabel "Energy - E_F (eV)"
436set ylabel "Local DOS (states/eV)"
437set style data lines
438set key top left Left reverse
439set border 31 lw 0.2
440set title "Projected DOS"
441plot \\
442 './AlAs110.pdos_atm#10(Al)' u (\$1-ef):2 t "Surface Al" ,\\
443 './AlAs110.pdos_atm#11(As)' u (\$1-ef):2 t "Surface As"
444set title "Local DOS centered in the first vacuum layer"
445plot \\
446 './AlAs110.ldos_boxes' u (\$1-ef):4 t "Above Al" ,\\
447 './AlAs110.ldos_boxes' u (\$1-ef):5 t "Above As" ,\\
448 './AlAs110.ldos_boxes' u (\$1-ef):(\$7/54) t "Surface average"
449set title "Local DOS centered in the second vacuum layer"
450plot \\
451 './AlAs110.ldos_boxes' u (\$1-ef):8 t "Above Al" ,\\
452 './AlAs110.ldos_boxes' u (\$1-ef):9 t "Above As" ,\\
453 './AlAs110.ldos_boxes' u (\$1-ef):(\$11/54) t "Surface average"
454EOF
455$ECHO
456$ECHO "  plotting DOS results ...\c"
457$GP_COMMAND < gnuplot.tmp
458$ECHO " done"
459rm gnuplot.tmp
460fi
461
462$ECHO
463$ECHO "  To visualize a volume in which the DOS is integrated, execute:"
464$ECHO "    xcrysden --xsf results/AlAs110.box#1.xsf"
465$ECHO "  and plot the isosurface corresponding to isovalue 0.5"
466
467$ECHO
468# clean TMP_DIR
469$ECHO "  cleaning $TMP_DIR...\c"
470rm -rf $TMP_DIR/AlAs110.*
471$ECHO
472$ECHO " $EXAMPLE_DIR: done"
473