Quantum Espresso meets bash : Primer and tidying up
March 13, 2012 Posted by Emre S. Tasci
In my previous postdoc position, DFT calculations were governing my research field but with my current position at the Bilbao Crystallographic Server, we are more on the theoretical side of things then calculations. Nevertheless, recently, my old ‘flame’ have been kindled and I decided to scrape the rust.
As a change, I selected and installed the Quantum Espresso (QE) package instead of VASP that I was using in the past, main reason QE being open source.
As a newbie to QE and long-gone returner to DFT, I started to practice by mostly ye olde trial/error approach.
The first thing that surprised me the lack of a commenting option in regard to QE input files! (Actually, I’m almost sure that it surely is supported but I have yet to find!). Here’s my sample input file:
&control title='Diamond Relaxing', # Calculation Type ========== calculation='scf' #calculation='relax', #calculation='vc-relax' #calculation='cp' # Restart =================== restart_mode='from_scratch', #restart_mode='restart' # Read/Write units ========== #ndr = 51, #ndw = 51, #nstep = 100, #nstep = 1000, nstep = 2000, # max_seconds = 600 # jobs stops after max_seconds CPU time. iprint = 10, # band energies are written every iprint iterations #tstress = .TRUE., # calculate stress. It is set to .TRUE. automatically if # calculation='vc-md' or 'vc-relax' #tprnfor = .TRUE., # print forces. Set to .TRUE. if calculation='relax','md','vc-md' #dt = 5.0d0, # time step for molecular dynamics, in Rydberg atomic units # (1 a.u.=4.8378 * 10^-17 s : beware, the CP code use # Hartree atomic units, half that much!!!) prefix='diamond_rel', # prepended to input/output filenames pseudo_dir = '/usr/share/espresso-4.3.2/pseudo/', outdir='/home/sururi/shared/tmp/' etot_conv_thr = 1.d-6 # convergence threshold on total energy (a.u) for ionic minimization forc_conv_thr = 1.d-3 # convergence threshold on forces (a.u) for ionic minimization / &system ## ibrav ## 0 "free", see above not used ## 1 cubic P (sc) not used ## 2 cubic F (fcc) not used ## 3 cubic I (bcc) not used ## 4 Hexagonal and Trigonal P celldm(3)=c/a ## 5 Trigonal R, 3fold axis c celldm(4)=cos(alpha) ## -5 Trigonal R, 3fold axis <111> celldm(4)=cos(alpha) ## 6 Tetragonal P (st) celldm(3)=c/a ## 7 Tetragonal I (bct) celldm(3)=c/a ## 8 Orthorhombic P celldm(2)=b/a,celldm(3)=c/a ## 9 Orthorhombic base-centered(bco) celldm(2)=b/a,celldm(3)=c/a ## 10 Orthorhombic face-centered celldm(2)=b/a,celldm(3)=c/a ## 11 Orthorhombic body-centered celldm(2)=b/a,celldm(3)=c/a ## 12 Monoclinic P, unique axis c celldm(2)=b/a,celldm(3)=c/a, ## celldm(4)=cos(ab) ## -12 Monoclinic P, unique axis b celldm(2)=b/a,celldm(3)=c/a, ## celldm(5)=cos(ac) ## 13 Monoclinic base-centered celldm(2)=b/a,celldm(3)=c/a, ## celldm(4)=cos(ab) ## 14 Triclinic celldm(2)= b/a, ## celldm(3)= c/a, ## celldm(4)= cos(bc), ## celldm(5)= cos(ac), ## celldm(6)= cos(ab ibrav= 2, # Cell dimensions ======================== # a,b,c in Angstrom ======= A = 3.56712 B = 2.49 C = 2.49 cosAB=0.5 cosAC=0.5 cosBC=0.5 ## OR ## ## celldm(1:6) in Bohr ===== ## 1 Bohr = 0.529177249 A ## 1 A = 1.889725989 Bohr # celldm(1) = 6.7409 # alat # celldm(2) = 4.7054 # b/a # celldm(3) = 4.7054 # c/a # celldm(4) = 0.5 # cos(ab) # celldm(5) = 0.5 # cos(ac) # celldm(6) = 0.5 # cos(ab) nat=2, ntyp= 1, ecutwfc = 40, # kinetic energy cutoff (Ry) for wavefunctions # ecutrho = 160 # kinetic energy cutoff (Ry) for charge density and potential. # For norm-conserving pseudopotential you should stick to the # default value # nbnd = 6 # number of electronic states (bands) to be calculated. # Default: # for an insulator, nbnd = number of valence bands (nbnd = # of electrons /2); # for a metal, 20% more (minimum 4 more) / &electrons emass = 400.d0 emass_cutoff = 2.5d0 #electron_dynamics = 'sd' # CP specific #electron_dynamics = 'bfgs' # CP specific conv_thr = 1D-6 # Convergence threshold for selfconsistency / &ions #ion_dynamics = 'none' #ion_dynamics = 'bfgs' #ion_dynamics = 'sd' # === damping === 0 == # ion_dynamics = 'damp' # ion_damping = 0.2 # ion_velocities = 'zero' # CP specific -- initial ionic velocities # === damping === 1 == #ion_nstepe = 10 # CP specific -- number of electronic steps per ionic step. (Def: 1) / &cell cell_dynamics = 'none', #cell_dynamics = 'bfgs', # ion_dynamics must be 'bfgs' too #cell_dofree = 'xyz' #cell_factor = 3 press = 0.0d, # press_conv_thr = 0.05 # Convergence threshold on the pressure for variable cell relaxation / ATOMIC_SPECIES C 12.0107 C.pz-vbc.UPF ATOMIC_POSITIONS C 0.0 0.0 0.0 C 0.25 0.25 0.25 K_POINTS automatic 4 4 4 1 1 1
so, it’s not the tidiest parameter file around but nevertheless, I’ve included some explanations and it’s easily customizable by simple commenting out, right?…
And then the problem of (yet undiscovered) comments in the QE input files hits us. But with bash, it’s pretty easy. So, if you save the above contents to a file called, say, “diamond.scf.in” and filter it via grep and sed, you’re done:
sururi@husniya:~/shared/qe/diamond_tutorial$ cat diamond.scf.in | grep -v -E "^[\t ]*#" | sed "s:\([^#]*\)#.*:\1:" &control title='Diamond Relaxing', calculation='scf' restart_mode='from_scratch', nstep = 2000, iprint = 10, prefix='diamond_rel', pseudo_dir = '/usr/share/espresso-4.3.2/pseudo/', outdir='/home/sururi/shared/tmp/' etot_conv_thr = 1.d-6 forc_conv_thr = 1.d-3 / &system ibrav= 2, A = 3.56712 B = 2.49 C = 2.49 cosAB=0.5 cosAC=0.5 cosBC=0.5 nat=2, ntyp= 1, ecutwfc = 40, / &electrons emass = 400.d0 emass_cutoff = 2.5d0 conv_thr = 1D-6 / &ions / &cell cell_dynamics = 'none', press = 0.0d, / ATOMIC_SPECIES C 12.0107 C.pz-vbc.UPF ATOMIC_POSITIONS C 0.0 0.0 0.0 C 0.25 0.25 0.25 K_POINTS automatic 4 4 4 1 1 1
I’ve installed the QE to take advantage of the cluster nodes via MPI, so I submit my jobs via:
mpirun -n #_of_processes /usr/share/espresso-4.3.2/bin/pw.x < input_file > output_file
which I guess is more or less similar to what you are doing to run your jobs.
So, combining this “comment filtering” and job submission, I’ve written the following bash script (“mpi-pw.x“) to automatize (and keep log of) the job:
sururi@cowboy:~$ cat /home/sururi/bin/mpi-pw.x
#!/bin/bash
# Script to submit Quantum Espresso jobs via mpi
# EST, Fri Mar 9 17:02:10 CET 2012
function timer()
{
if [[ $# -eq 0 ]]; then
echo $(date '+%s')
else
local stime=$1
etime=$(date '+%s')
if [[ -z "$stime" ]]; then stime=$etime; fi
dt=$((etime - stime))
ds=$((dt % 60))
dm=$(((dt / 60) % 60))
dh=$((dt / 3600))
printf '%d:%02d:%02d' $dh $dm $ds
fi
}
if [ $# -ne 2 ]
then
echo "Usage: mpi-pw.x num_of_processes infile.in"
exit
fi
job=$2
job_wo_ext=${job%.*}
outext=".out"
logext=".log"
echo -e "Running: mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < $2 > $job_wo_ext$outext\n"
echo -e "Running: mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < $2 > $job_wo_ext$outext\n">$job_wo_ext$logext
echo "Start : "$(date)
echo "Start : "$(date)>>$job_wo_ext$logext
# Take out the comments from the input file
#grep -v -E "^[\t ]*#" $2 > in_wo_comments.in.tmp
grep -v -E "^[\t ]*#" $2 | sed "s:\([^#]*\)#.*:\1:" > in_wo_comments.in.tmp
t=$(timer)
mpirun -n $1 /usr/share/espresso-4.3.2/bin/pw.x < in_wo_comments.in.tmp > $job_wo_ext$outext
echo "Finish: "$(date)
echo "Finish: "$(date)>>$job_wo_ext$logext
printf 'Elapsed time:%s\n\n' $(timer $t)
printf 'Elapsed time:%s\n\n' $(timer $t) >>$job_wo_ext$logext
rm -f in_wo_comments.in.tmp(I’ve picked the “timer” function from Mitch Frazier’s entry in LinuxJournal)
What it does is, really simple: it strips out the extension of the input file ($2 – the second argument in calling); filters out the comments, constructing the temporary input file “in_wo_comments.in.tmp”; starts the timer; submits the job using the number of processes passed in the calling ($1 – first argument) and directs the output to a file whose name is the extension stripped input filename + “.out”; when the job is (this way or that) done, stops the timer and calculates the elapsed time via the ‘timer’ function. During all this time, it prints these information both on screen and to a file called extension stripped input filename + “.log”.
I have a very similar another script called “mpi-cp.x” that I’ve forked for Carr-Parinelli jobs (just search & replace all the occurrences of “mpi-pw.x” with “mpi-cp.x” and you’re done 8).
And here is the script in action:
sururi@cowboy:~/shared/qe/diamond_tutorial$ ll total 16 drwxr-xr-x 2 sururi sururi 4096 2012-03-13 15:59 ./ drwxr-xr-x 7 sururi sururi 4096 2012-03-13 15:10 ../ -rw-r--r-- 1 sururi sururi 4954 2012-03-13 15:59 diamond.scf.in sururi@cowboy:~/shared/qe/diamond_tutorial$ mpi-pw.x 4 diamond.scf.in Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < diamond.scf.in > diamond.scf.out Start : Tue Mar 13 16:18:20 CET 2012 Finish: Tue Mar 13 16:18:21 CET 2012 Elapsed time:0:00:01 sururi@cowboy:~/shared/qe/diamond_tutorial$ grep -e ! diamond.scf.out ! total energy = -22.70063050 Ry sururi@cowboy:~/shared/qe/diamond_tutorial$ cat diamond.scf.log Running: mpirun -n 4 /usr/share/espresso-4.3.2/bin/pw.x < diamond.scf.in > diamond.scf.out Start : Tue Mar 13 16:18:20 CET 2012 Finish: Tue Mar 13 16:18:21 CET 2012 Elapsed time:0:00:01
Now that we have the mpi-pw.x and mpi-cp.x, our next entry in the series will be the automated search for optimized values.
March 13, 2012 at 11:43 pm
[…] 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. […]
November 19, 2012 at 10:54 am
Comments in namelists can be introduced using “!”.
Comments in “cards” can be introduced by a “#” in the first column.
Paolo