How to read QE's wfc*.dat files with python

Question

I'm working with QE 7.0 for the computation of some properties of solids. My working theory requires me to access all the information inside the wfc#c.dat/h5 files, created in the prefix.save directory for further post-processing. For the post-processing I created a python script, that works well with the .h5 files.
My issue is that for a few of my system the calculations were made with a version of QE that was compiled without hdf5, and since the calculation is very large for the size of the HPC I have access to, I prefer not to redo the calculation with another compilation of QE.

I'm now trying to read the data in the wfc#k.dat files, but can't figure out exactly how. I've tried using scipy.io.FortranFile but I can't manage to figure it out. I've also tried looking for converters for the .dat files to .h5 format, but couldn't find any.

I have close to no experience with Fortran, and no experience at all with Fortran unformatted binary files, so the information in the QE snippets https://gitlab.com/QEF/q-e/-/snippets/1869208 and https://gitlab.com/QEF/q-e/-/snippets/1869202 are of no help to me.

I'd highly appreciate any help with this issue.

edit: In the meanwhile I've managed to redo the calculation with a hdf5-enabled version, but out of curiosity I'd still appreciate insights into this.

Answer 1

With python, you can open a binary file using the 'rb' option in your open statement. Then using numpy.fromfile you can read the file provided you know the datatypes that are in the file.

First you have to know that, when you write an unformatted file in Fortran, each time you use a write statement, the program will also write, before the data, 4 bytes that tells the size of the data that is contained in the "data block", then the actual data, and then on 4 bytes again, the size of the data that is contained in the "data block". As this size doesn't matter in your case, you can just skip all of these blocks using python's file.seek method that allows you to move the cursor.

Looking at the wiki, you'll see that the format of the wfc file is the following :

ik, xk, ispin, gamma_only, scalef
ngw, igwx, npol, nbnd
b1, b2, b3
mill(1:3,1:igwx)
evc(1:npoligwx,1)
evc(1:npoligwx,2)
...
evc(1:npol*igwx,nbnd)

Where the datatypes are :

  INTEGER :: ik
  REAL(8) :: xk(3)
  INTEGER :: ispin
  LOGICAL :: gamma_only
  REAL(8) :: scalef
  INTEGER :: ngw
  INTEGER :: igwx
  INTEGER :: npol
  INTEGER :: nbnd
  REAL(8) :: b1(3), b2(3), b3(3)
  INTEGER :: mill(3,igwx)
  COMPLEX(8) :: evc(npol*igwx,nbnd)

Then, you'll need to know the conversion from the fortran types to the numpy ones

• INTEGER becomes int32
• REAL(8) becomes float64
• COMPLEX(8) becomes complex128
• LOGICAL is written on 4 bytes, you can just read it as an int32

Here's a quick and dirty way to read these files, you can group the np.fromfile if you define correctly the types, but I thought that it would be clearer doing it step by step.

import numpy as np

with open('wfc1.dat', 'rb') as f:
    # Moves the cursor 4 bytes to the right
    f.seek(4)

    ik = np.fromfile(f, dtype='int32', count=1)[0]
    xk = np.fromfile(f, dtype='float64', count=3)
    ispin = np.fromfile(f, dtype='int32', count=1)[0]
    gamma_only = bool(np.fromfile(f, dtype='int32', count=1)[0])
    scalef = np.fromfile(f, dtype='float64', count=1)[0]

    # Move the cursor 8 byte to the right
    f.seek(8, 1)

    ngw = np.fromfile(f, dtype='int32', count=1)[0]
    igwx = np.fromfile(f, dtype='int32', count=1)[0]
    npol = np.fromfile(f, dtype='int32', count=1)[0]
    nbnd = np.fromfile(f, dtype='int32', count=1)[0]

    # Move the cursor 8 byte to the right
    f.seek(8, 1)

    b1 = np.fromfile(f, dtype='float64', count=3)
    b2 = np.fromfile(f, dtype='float64', count=3)
    b3 = np.fromfile(f, dtype='float64', count=3)

    f.seek(8,1)
    
    mill = np.fromfile(f, dtype='int32', count=3*igwx)
    mill = mill.reshape( (igwx, 3) ) 

    evc = np.zeros( (nbnd, npol*igwx), dtype="complex128")

    f.seek(8,1)
    for i in range(nbnd):
        evc[i,:] = np.fromfile(f, dtype='complex128', count=npol*igwx)
        f.seek(8, 1)

Hope it helped.