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