Parse Quantum Espresso Output File

Hello All,

Sorry for the large delays between updates to this blog. I have been really busy with school, and I will make it a habit to write a weekly blog post from here on it. I just added a calendar event to remind me next week.

I have been using python for a while now, and I have just come to fall in love with it (along with every other script kid out there). I am currently working on a project to store all our group’s simulation data in a search-able database. I want to make something easy for our group to upload to the web and will store the important information:

  • Relaxed Atomic Positions
  • Total Energy of the System
  • Bandgap
  • Lattice Constants
  • Band Diagram
  • Fermi Level

And other information such as what code version of the code ran, etc. I was playing around with this for a bit, and I came up with a nifty python script that extracts all the necessary information out of a Quantum Espresso output file. I eventually want to make this for other programs such as VASP, GPAW, etc. All the units are in Angstrom and electron volts.

I have tested this on Quantum Espresso 5.2. Download here: Code Download

#!/usr/bin/python3.4
#Created by Levi Lentz
#Covered under the MIT license
import sys
import numpy as np
import json
class Struct:
  def __init__(self):
    self.BOHRtoA = 0.529177249
    self.RYtoeV = 13.605698066
    self.program_version = ""
    self.volume = 0.0
    self.alat = 0.0
    self.natoms = 0
    self.nat = 0
    self.nelect = 0
    self.Ecut = 0.0
    self.RhoCut = 0.0
    self.Econv = 0.0
    self.beta = 0.0
    self.Exch = ""
    self.energy = 0.0
    self.natoms = 0
    self.bandgap = 0.0
    self.bands = 0
    self.lattice = {'a':np.zeros(3),'b':np.zeros(3),'c':np.zeros(3)}
    self.atoms = []
    self.norms = {'a':0.0,'b':0.0,'c':0.0}
    self.angles = {'alpha':0.0,'beta':0.0,'gamma':0.0}
    self.kpts = 0
    self.bnddiagram = np.zeros(0)
    self.FermiTest = False
    self.Fermi = 0.0
  def Process(self,filestring):
    f = open(filestring,'r')
    linenum = 0
    for i in f:
      if linenum < 1000:
        if "number of k points=" in i:
          self.kpts = int(i.split()[4])
          next
        if "Program PWSCF" in i:
          self.program_version = i.split()[2]
          next
        if "lattice parameter (alat)" in i:
          self.alat = float(i.split()[4])*self.BOHRtoA
          next
        if "number of Kohn-Sham states" in i:
          self.bands = int(i.split()[4]) 
        if "unit-cell volume" in i and "new" not in i:
          self.volume = float(i.split()[3])*(self.BOHRtoA**3.0)
          next
        if "number of atoms/cell" in i:
          self.natoms = int(i.split()[4])
          next
        if "number of atomic types" in i:
          self.nat = int(i.split()[5])
          next
        if "number of electrons" in i:
          self.nelect = float(i.split()[4])
          next
        if "kinetic-energy cutoff" in i:
          self.Ecut = float(i.split()[3])*self.RYtoeV
          next
        if "charge density cutoff" in i:
          self.RhoCut = float(i.split()[4])*self.RYtoeV
          next
        if "convergence threshold" in i:
          self.Econv = float(i.split()[3])
          next
        if "mixing beta" in i:
          self.beta = float(i.split()[3])
          next
        if "Exchange-correlation" in i:
          self.Exch = i.split()[2]
          next
        if "a(1) =" in i:
          tmp = i.split()
          for j in range(0,3):
            self.lattice['a'][j] = tmp[j+3]
          next
        if "a(2) =" in i:
          tmp = i.split()
          for j in range(0,3):
            self.lattice['b'][j] = tmp[j+3]
          next
        if "a(3) =" in i:
          tmp = i.split()
          for j in range(0,3):
            self.lattice['c'][j] = tmp[j+3] 
          next
        if "site n.     atom                  positions (alat units)" in i:
          for j in range(0,self.natoms):
            line = next(f).split()
            self.atoms.append([line[1],float(line[6])*self.alat,float(line[7])*self.alat,float(line[8])*self.alat])
          next
      if "!" in i:
        self.energy= float(i.split()[4])*self.RYtoeV
      if "new unit-cell volume" in i:
        self.volume = float(i.split()[4])*(self.BOHRtoA**3)
      if "Begin final coordinates" in i:
        #print(self.atoms)
        while "End final coordinates" not in line:
          line = next(f)
          if "CELL_PARAMETERS" in line:
            self.alat = float(line.replace(')','').split()[2])*self.BOHRtoA
            for j in ['a','b','c']:
              line = next(f)
              tmp = line.split()
              for k in range(0,3):
                self.lattice[j][k] = self.alat*float(tmp[k])
          if "ATOMIC_POSITIONS" in line:
            if "angstrom" in line:
              for j in range(0,self.natoms):
                line = next(f).split()
                self.atoms[j] = [line[0],float(line[1]),float(line[2]),float(line[3])]
      if "End of self-consistent calculation" in i: #The band diagram is stored in lines of 8 entries
        if np.floor(self.bands/8.)*8. <= self.bands:
          numlines = int(np.floor(self.bands/8.) + 1)
          remainder = int(self.bands - np.floor(self.bands/8.)*8.)
        else: 
          numlines = int(np.floor(self.bands/8.))
          remainder = 0
        self.bnddiagram = np.zeros((self.kpts,self.bands))
        counter = 0
        while counter < self.kpts:
          line = next(f)
          if "k =" in line:
            line = next(f)
            counter1 = 0
            for j in range(0,numlines):
              line = next(f)
              for k in range(0,len(line.split())):
                self.bnddiagram[counter][counter1 + k] = float(line.split()[k])
              counter1 += 8
            counter += 1  
        next
      if "highest occupied, lowest unoccupied level (ev)" in i:
        self.bandgap = float(i.split()[7]) - float(i.split()[6])
        next
      if "the Fermi energy is" in i:
        self.Fermi = float(i.split()[4])
        self.FermiTest = True
        next
      linenum += 1
    f.close()
    for i in ['a','b','c']:
      self.norms[i] = np.linalg.norm(self.lattice[i])
    self.angles['alpha'] = np.arccos(np.dot(self.lattice['b'],self.lattice['c'])/(self.norms['c']*self.norms['b'])) * 180./np.pi
    self.angles['gamma'] = np.arccos(np.dot(self.lattice['a'],self.lattice['b'])/(self.norms['a']*self.norms['b'])) * 180./np.pi
    self.angles['beta'] = np.arccos(np.dot(self.lattice['a'],self.lattice['c'])/(self.norms['a']*self.norms['c'])) * 180./np.pi
    if self.FermiTest == True: #The bandgap is now in the band diagram
      self.bnddiagram = np.subtract(self.bnddiagram,self.Fermi)
      emin = np.zeros(self.kpts)
      emax = np.zeros(self.kpts)
      counter = 0
      for j in self.bnddiagram:
        emin[counter] = j[np.searchsorted(j,  0.0,side='right')-1]
        emax[counter] = j[np.searchsorted(j,  0.0,side='right')]
        counter += 1
      self.bandgap = float(np.min(emax-emin))
  def to_JSON(self):
    for i in self.lattice:
      self.lattice[i] = self.lattice[i].tolist() 
    self.bnddiagram = self.bnddiagram.tolist() #JSON doesnt like numpy arrays
    return json.dumps(self, default=lambda o: o.__dict__, sort_keys=True, indent=4)

test = Struct()
test.Process("vc.out")

print(test.to_JSON())

The code itself is pretty self explanatory, the output of which looks like (bnd diagram removed):

{
"BOHRtoA": 0.529177249,
"Econv": 1e-08,
"Ecut": 816.34188396,
"Exch": "BLYP",
"Fermi": 0.0,
"FermiTest": false,
"RYtoeV": 13.605698066,
"RhoCut": 9796.10260752,
"alat": 5.430681517862499,
"angles": {
"alpha": 90.0,
"beta": 90.0,
"gamma": 90.0
},
"atoms": [
[
"Si",
0.0,
0.0,
0.0
],
[
"Si",
0.0,
2.752464825,
2.752464825
],
[
"Si",
2.752464825,
2.752464825,
0.0
],
[
"Si",
2.752464825,
0.0,
2.752464825
],
[
"Si",
4.12869281,
1.376226702,
4.12869281
],
[
"Si",
1.376226702,
1.376226702,
1.376226702
],
[
"Si",
1.376226702,
4.12869281,
4.12869281
],
[
"Si",
4.12869281,
4.12869281,
1.376226702
],
[
"Si",
0.0,
0.0,
0.0
],
[
"Si",
0.0,
2.7524556656287773,
2.7524556656287773
],
[
"Si",
2.7524556656287773,
2.7524556656287773,
0.0
],
[
"Si",
2.7524556656287773,
0.0,
2.7524556656287773
],
[
"Si",
4.128678882363876,
1.3762221305987947,
4.128678882363876
],
[
"Si",
1.3762221305987947,
1.3762221305987947,
1.3762221305987947
],
[
"Si",
1.3762221305987947,
4.128678882363876,
4.128678882363876
],
[
"Si",
4.128678882363876,
4.128678882363876,
1.3762221305987947
]
],
"bandgap": 1.3244999999999996,
"bands": 100,
"beta": 0.7,
"bnddiagram": [
[
]
],
"energy": -848.6786362150125,
"kpts": 10,
"lattice": {
"a": [
5.504913503530162,
0.0,
0.0
],
"b": [
0.0,
5.504913503530162,
0.0
],
"c": [
0.0,
0.0,
5.504913503530162
]
},
"nat": 1,
"natoms": 8,
"nelect": 32.0,
"norms": {
"a": 5.504913503530162,
"b": 5.504913503530162,
"c": 5.504913503530162
},
"program_version": "v.5.1",
"volume": 166.8227981355072
}

You could easily combine this with Excel, or SQL, or something else to easily catalog your calculation. If you are like me, you will have a huge need for this. The main issue left is that it needs to check for different coordinate systems for the atoms (currently it only checks for angstrom). This is relatively trivial to add; I will update this to the code with that once it is completed.

Happy simulations,

Levi

3 comments Add yours

Leave a Reply

Your email address will not be published. Required fields are marked *