SCF files the Python Way

Hello All,

Week two of my personal post-a-week challenge. I have been practicing automating my workflow and I have another quick script.

Using my QE module that I talked about last week, I made a quick script that will develop the SCF, NSCF, DOS, and PDOS files. These are the files needed for a density of states calculation. What is nice about this work flow is it allows you, in one run script, to run all these calculations. This is important for running on Super Computers as wait times can sometimes outweigh the actual calculations.

#!/usr/bin/python3.4

import QE as qe
import sys
import numpy as np

class Files:
  def __init__(self):
    self.prefix = ''
    self.pdos = """&projwfc
    outdir='.'
    prefix='PREFIX'
    io_choice='both'
    Emin=-100, Emax=100.0, DeltaE=0.1
    ngauss=1, degauss=0.02
 /
"""
    self.dos = """&dos
  outdir='.'
  prefix='PREFIX'
  fildos='PREFIX.dos'
  Emin=-100.0, Emax=100.0, DeltaE=0.1
/
"""
    self.scf = """"""
    self.nscf = """"""
    self.round = 5
  def Process(self,filestring,QE):
    try:
      f = open(filestring,'r')
    except:
      print("Cannot open %s" % filestring)
      sys.exit(1)
    self.scf = ''
    self.nscf = ''
    for i in f:
      if 'calculation' in i:
        self.scf += i.replace('vc-relax','scf').replace('relax','scf')
        self.nscf += i.replace('vc-relax','nscf').replace('relax','nscf')
      elif '&system' in i:
        tmpstr = '''&system\n  ibrav = 14\n'''
        for j in ['a','b','c']:
          tmpstr += "  " + j + " = " + str(round(QE.norms[j],self.round)) + "\n"
        tmpstr += "  cosBC = " + str(round(np.cos(QE.angles['alpha'] * np.pi/180.),self.round)) + "\n"
        tmpstr += "  cosAC = " + str(round(np.cos(QE.angles['beta'] * np.pi/180.),self.round)) + "\n" 
        tmpstr += "  cosAB = " + str(round(np.cos(QE.angles['gamma'] * np.pi/180.),self.round)) + "\n"
        self.scf += tmpstr
        self.nscf += tmpstr
        for j in range(0,7):
          next(f)
      elif 'ATOMIC_POSITIONS' in i.upper():
        tmpstr = 'ATOMIC_POSITIONS {angstrom}\n'
        for j in QE.atoms:
          tmpstr += j[0] + "\t" + str(j[1]) + "\t" + str(j[2]) + "\t" + str(j[3]) + "\n"
          next(f)
        self.scf += tmpstr
        self.nscf += tmpstr
      elif 'K_POINTS' in i.upper():
        self.scf += i 
        self.nscf += i 
        line = next(f)
        self.scf += line
        ks = line.split()
        for j in range(0,6):
          if j < 3:
            self.nscf += str(2*int(ks[j])) + " " 
          else:
            self.nscf += ks[j] + " "
        self.nscf += "\n"
      elif "prefix" in i.lower():
        self.prefix=i.replace('prefix','').replace('=','').replace("'","").replace('"','').strip()
      else:
        self.scf += i
        self.nscf += i
      self.pdos = self.pdos.replace('PREFIX',self.prefix)
      self.dos = self.dos.replace('PREFIX',self.prefix)
        
def main(command):
  QEStruct= qe.Struct()
  QEStruct.File_Process(command[1])
  outputs = Files()
  outputs.Process(command[0],QEStruct)
  f = open('scf.in','w')
  f.write(outputs.scf)
  f.close()
  f = open('nscf.in','w')
  f.write(outputs.nscf)
  f.close()
  f = open('dos.in','w')
  f.write(outputs.dos)
  f.close()
  f = open('pdos.in','w')
  f.write(outputs.pdos)
  f.close()

if __name__ == "__main__":
  if len(sys.argv) != 3:
    print("Incorrect number of arguments, run as ./Scf.py QEINPUT QEOUTPUT")
    sys.exit(6)
  command = [sys.argv[1],sys.argv[2]]
  main(command)

This generates a scf.in, nscf.in, dos.in, pdos.in, which allows you to run one input script for all files as follows:

pw.x < relax.in > relax.out
./ScfCreator.py relax.in relax.out
pw.x < scf.in > scf.out
pw.x < nscf.in > nscf.out
dos.x < dos.in > dos.out
projwfc.x < pdos.in > pdos.out

Eventually I hope to develop this into a python wrapper for the QE program.

Happy computing,

Levi