# Hex, Bugs and More Physics | Emre S. Tasci

a blog about physics, computation, computational physics and materials…

## Quantum Espresso meets bash : Automated running for various values

### March 13, 2012 Posted by Emre S. Tasci

Now that we can include comments in our input files and submit our jobs in a tidied up way, it’s time to automatize the “convergence search” process. For this, we’ll be using the same diamond example of the previous entry.

Suppose that we want to calculate the energies corresponding to various values of the lattice parameter a, ranging from 3.2 to 4, incrementing by 2: 3.2, 3.4, … , 4.0

Incrementation:
In bash, you can define a for loop to handle this:

```for (( i = 1; i &lt;= 9; i++ ))
do
echo \$i;
done

1
2
3
4
5
6
7
8
9```

But to work with rational numbers, things get a little bit complex:
```for (( i = 1; i &lt;= 5; i++ ))
do
val=\$(awk "BEGIN {print 3+\$i/5 }")
echo \$val;
done

3.2
3.4
3.6
3.8
4```

Alas, it is still not very practical to everytime calculate the necessary values to produce our intended range. This is where the mighty Octave once again comes to the rescue: [A:I:B] in Octave defines a range that starts from A, and continues by incrementing by I until B is reached (or exceeded). So, here is our example in Octave:
```octave:1&gt; [3.2:0.2:4]
ans =

3.2000   3.4000   3.6000   3.8000   4.0000```

You can generate temporary input scripts for Octave but better yet, you can also feed simple scripts directly to Octave:
```sururi@husniya:~/shared/qe/diamond_tutorial\$ echo "[3.2:0.2:4]"|octave -q
ans =

3.2000   3.4000   3.6000   3.8000   4.0000```

The “-q” parameter causes Octave to run in “quiet” mode, so we just feed the range via pipe to Octave and will harvest the output. In general, we’d like to keep the “ans =” and the blank lines out, and then, when there are more than 1 line’s worth of output, Octave puts the “Column” information as in:
```octave:4&gt; [3:0.1:4]
ans =

Columns 1 through 9:

3.0000   3.1000   3.2000   3.3000   3.4000   3.5000   3.6000   3.7000   3.8000

Columns 10 and 11:

3.9000   4.0000```

all these nuisances can go away via a simple “grep -v” (with the ‘v’ parameter inverting the filter criteria):
```sururi@vala:/tmp\$ echo "[3.2:0.1:4]"|octave -q|sed -n "2,1000p"| grep -v "Column"

3.2000   3.3000   3.4000   3.5000   3.6000   3.7000   3.8000   3.9000

4.0000

sururi@vala:/tmp\$ echo "[3.2:0.1:4]"|octave -q|sed -n "2,1000p"| grep -v "Column"|grep -v -e "^\$"
3.2000   3.3000   3.4000   3.5000   3.6000   3.7000   3.8000   3.9000
4.0000```

Replacement:
Now that we have the variable’s changing values, how are we gonna implement it to the input script? By search and replace. Suppose, we have a template in which the variable’s value is specified by “#REPLACE#”, all we need to do is to “sed” the file:
`cat diamond.scf.in | sed "s:#REPLACE#:\$a:g"&gt;tmp_\$a.in`

where “\$a” is the current value in the iteration.

Parsing the Energy:
We submit this job via mpirun (in our case via mpi-pw.x of the previous post) and parse the output for the final energy:

`grep -e ! tmp_\$a.out|tail -n1|awk '{print \$(NF-1)}'`

For each iteration we thus obtain the energy for our \$a value, so we put them together in two columns into a data file I’d like to call “E_vs_\$a.txt”.

The Script:
Putting all these together, here is the script in its full glory:

```#!/bin/bash

# Takes in the template Quantum Espresso input file, runs it iteratively
# changing the #REPLACE# parameters iteratively with the range elements,
# grepping the energy and putting it along with the current value of the
# #REPLACE# parameter.
#
# Emre S. Tasci, 03/2012

# Example: Suppose that the file "diamond.scf.in" contains the following line:
#   A = #REPLACE#
# Calling "qe-opti.sh diamond.scf.in a 3.2 0.2 4.2" would cause that line to be modified as
#   A = 3.2
# and submitted, followed by
#   A = 3.4, 3.6, .., 4.2 in the sequential runs.
#
# at the end, having constructed the E_vs_a.txt file.

if [ \$# -ne 5 ]
then
echo "Usage: qe-opti.sh infile_template.in variable_name initial_value increment     final_value"
echo "Example: qe-opti.sh diamond.scf.in a 3.2 0.2 4.2"
echo -e "\nNeeds octave to be accessible via \"octave\" command for the range to     work"
echo -e "\nEmre S. Tasci, 03/2012"
exit
fi

logext=".txt"
logfilename=E_vs_\$2\$logext
rm -f \$logfilename

range=\$(echo "[\$3:\$4:\$5]"|octave -q|sed -n "2,1000p"|grep -vr "^\$"|grep -v "Column")
for a in \$range
do
# echo \$a
cat \$1|sed "s:#REPLACE#:\$a:g"&gt;tmp_\$a.in
/home/sururi/bin/mpi-pw.x 4 tmp_\$a.in
energ=\$(grep -e ! tmp_\$a.out|tail -n1|awk '{print \$(NF-1)}')
# echo \$energ
echo -e \$a "    " \$energ &gt;&gt; \$logfilename
echo -e \$a "    " \$energ
done

echo "Results stored in: "\$logfilename```

Action:
Modify the sample input from the previous post so that it’s lattice parameter A is designated as something like “A = #REPLACE#”, as in:
```sururi@bebop:~/shared/qe/diamond_tutorial\$ grep "REPLACE" -B2 -A2 diamond.scf.in
# a,b,c in Angstrom =======
#A = 3.56712
A = #REPLACE#
B = 2.49
C = 2.49```

Then call the qe-opti.sh script with the related parameters:
```sururi@bebop:~/shared/qe/diamond_tutorial\$ qe-opti.sh diamond.scf.in a 3.2 0.2 4.2
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x &lt; tmp_3.2000.in &gt; tmp_3.2000.out

Start : Tue Mar 13 23:21:07 CET 2012
Finish: Tue Mar 13 23:21:08 CET 2012
Elapsed time:0:00:01

3.2000      -22.55859979
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x &lt; tmp_3.4000.in &gt; tmp_3.4000.out

Start : Tue Mar 13 23:21:08 CET 2012
Finish: Tue Mar 13 23:21:09 CET 2012
Elapsed time:0:00:01

3.4000      -22.67717405
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x &lt; tmp_3.6000.in &gt; tmp_3.6000.out

Start : Tue Mar 13 23:21:09 CET 2012
Finish: Tue Mar 13 23:21:10 CET 2012
Elapsed time:0:00:01

3.6000      -22.69824702
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x &lt; tmp_3.8000.in &gt; tmp_3.8000.out

Start : Tue Mar 13 23:21:10 CET 2012
Finish: Tue Mar 13 23:21:11 CET 2012
Elapsed time:0:00:01

3.8000      -22.65805473
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x &lt; tmp_4.0000.in &gt; tmp_4.0000.out

Start : Tue Mar 13 23:21:11 CET 2012
Finish: Tue Mar 13 23:21:13 CET 2012
Elapsed time:0:00:02

4.0000      -22.57776354
Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x &lt; tmp_4.2000.in &gt; tmp_4.2000.out

Start : Tue Mar 13 23:21:13 CET 2012
Finish: Tue Mar 13 23:21:14 CET 2012
Elapsed time:0:00:01

4.2000      -22.47538891
Results stored in: E_vs_a.txt```

Where the “E_vs_a.txt” contains the … E vs a values!:
```sururi@bebop:~/shared/qe/diamond_tutorial\$ cat E_vs_a.txt
3.2000      -22.55859979
3.4000      -22.67717405
3.6000      -22.69824702
3.8000      -22.65805473
4.0000      -22.57776354
4.2000      -22.47538891```

You can easily plot it using GNUPlot or whatever you want. In my case, I prefer GNUPlot:
`sururi@husniya:~/shared/qe/diamond_tutorial\$ plot E_vs_a.txt`

giving me:

Other than the ‘plot’, I’ve also written another one, fit for data sets just like these where you have a parabol-like data and you’d like to try fitting it and here is what it does (automatically):

```sururi@husniya:~/shared/qe/diamond_tutorial\$ plotfitparabol.sh E_vs_a.txt
func    f(x)=0.674197*x*x-4.881275*x-13.855234
min, f(min)    3.620066 -22.690502```

Here are the source codes of the ‘plot’ and ‘plotfitparabol.sh’ scripts. Doei!

```sururi@mois:/vala/bin\$ cat plot
#!/bin/bash
gnuplot -persist <<END
set term postscript enhanced color
set output "\$1.ps"
plot "\$1" \$2 \$3

set term wxt
replot
END
sururi@mois:/vala/bin\$ cat plotfitparabol.sh
#!/bin/bash

# Usage : plotfitparabol.sh <datafile>

# Reads the data points from the given file
# then uses Octave's "polyfit" function to calculate the coefficients
# for the quadratic fit and also the minimum (maximum) value.
# Afterwards, prepares a GNUPlot script to display:
#    the data set (scatter)
#    fit function
#    the extremum point's value
#
# Written with Quantum Espresso optimizing procedures in mind.

if [ \$# -ne 1 ]
then
echo "Usage: plotfitparabol.sh <datafile>"
echo "Example: plotfitparabol.sh E_vs_A.txt"
echo -e "\nNeeds octave to be accessible via \"octave\" command to work"
echo -e "\nEmre S. Tasci <emre.tasci@ehu.es>, 03/2012"
exit
fi

oct=\$(echo \$octave_script|octave -q)
#echo \$oct

func=\${oct/|*/ }
minfmin=\${oct/*|/ }
echo -e "func \t" \$func
echo -e "min, f(min) \t" \$minfmin

echo '\$func
set term postscript enhanced color
set label "\$1" at graph 0.03, graph 0.94
set output "\$1.ps"
plot "\$1" title "data", f(x) title "\$func", "<echo '\$minfmin'" title "\$minfmin"

set term wxt
replot'

gnuplot -persist <<END
\$func
set term postscript enhanced color
set label "\$1" at graph 0.03, graph 0.94
set output "\$1.ps"
plot "\$1" title "data", f(x) title "\$func", "<echo '\$minfmin'" title "\$minfmin"

set term wxt
replot
END
```

### One Response to “Quantum Espresso meets bash : Automated running for various values”

1. Mustafa Kurban Says:

Thank you for the information !