Supercell Builder from QE Output

To download the new version of the QE.py code, click here: www.levilentz.com/Codes/QE.py

Hello All,

As I read about the new discovery of gravity waves at MIT (check it out: http://news.mit.edu/2016/ligo-first-detection-gravitational-waves-0211), I figured I could make myself feel a bit more impactful by adding into my QE module one other important feature.

Recently I have been studying machine learning options for material science applications. This involves generation large sets of data, and with all the computational power at my fingertips, why not automate this so that we can properly abuse large clusters such as Stampede, Edison, and our local computational resources. One necessary extension is the ability to generate a N1xN2xN3 supercell of a fully relaxed system.

  def to_Supercell(self,array):
    if isinstance(array,list): #Want to make sure that we have a list for the (Nx,Ny,Nz) array
      tmp = copy.copy(self) # I dont want to edit the existing structure, I want to return a supercell
      #The following is for conversion to crystal coordinates
      conversion = np.linalg.inv(tmp.From_Crystal())
      counter = 0
      for i in tmp.atoms:
        tmpa = np.array([i[1],i[2],i[3]])
        dot = np.dot(conversion,tmpa)
        tmp.atoms[counter] = [i[0],dot[0],dot[1],dot[2]]
        counter += 1
      COORDs = tmp.atoms[:]
      # Populating the atom positions
      for i in tmp.atoms:
        for j in range(0,array[0]):
          for k in range(0,array[1]):
            for z in range(0,array[2]):
              if (j == 0) and (k == 0) and (z == 0): #We already have the (0,0,0) structure
                next
              else:
                COORDs.append([i[0],i[1] + j, i[2] + k,i[3] + z])
      counter = 0
      #expanding the lattice
      for i in COORDs:
        tmpa = np.array([i[1],i[2],i[3]])
        dot = np.dot(tmp.From_Crystal(),tmpa)
        COORDs[counter] = [i[0],dot[0],dot[1],dot[2]]
        counter += 1
      for i in tmp.lattice:
        tmp.lattice[i] = np.array([tmp.lattice[i][0]*float(array[0]), tmp.lattice[i][1]*float(array[1]), tmp.lattice[i][2]*float(array[2])])
      tmp.atoms = COORDs
      tmp.natoms = tmp.natoms*(array[0]*array[1]*array[2])
      tmp.energy *= float(array[0])*float(array[1])*float(array[2])
      indexing = ['a','b','c']
      for i in indexing:
        tmp.norms[i] *= array[indexing.index(i)]
      
      return tmp
    else:
      print("Invalid supercell dimensions")

This portion of code simply takes in a list of number, [Nx,Ny,Nz] and returns a QE.Struct expanded to the proper size. The only trick is that we do all of this in crystallographic coordinates, but QE.Struct stores data internally in A and eV. So first we make a copy of the current QE.Struct, convert the coordinates to crystal, expand the system, then convert it back to Angstrom. I know it seems like a bear to have to do that, but I prefer to store the data internally in Angstrom and eV as it makes control QE a bit easier.

Using this, it allows you to do something such as:

#!/usr/bin/env python3.4

import QE as qe

f = qe.Struct()
f.File_Process('scf.out')

f.CIF("1.cif")
test = f.to_Supercell([5,5,2])
print("\n\n")
test.CIF("2.cif")

Which would give you two cif files, one that is the original calculation and the other which is a 5x5x2 supercell of the data in scf.out. This has been great for me to quickly generate large calculations with defects, deformations, etc. Hopefully my next module for QE will be to reduce it down to the primitive unit cell, but that is a bit more tricky and mathy.

Happy computing,

Levi