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