irreducible-representations / irrep

GNU General Public License v3.0
59 stars 30 forks source link

Modifications for hdf5 support? #38

Closed MSteggles closed 2 years ago

MSteggles commented 2 years ago

Hi!

First, thanks so much for making such a great tool. I'm using BandUpPy for band unfolding, and I've got it working fine on a work desktop. But, our university cluster uses the more modern hdf5 format on our compilation of Quantum Espresso. This causes kpoint.py in Irrep to fail

"File "**/lib/python3.8/site-packages/irrep/kpoint.py", line 907, in __init_espresso"

because it expects a .dat for the wavefunctions.

I'm not familiar with the formatting of the hdf5 or dat files, or with fortran at all - could you perhaps outline what would need to be done to mod in hdf5 support please?

All the best, M

MSteggles commented 2 years ago

Fixed - if you use the h5py library it's actually super simple and you don't have to interact with the underlying code very much at all: my solution is a bit janky but works fine. If you modify the kpoint.py __init_espresso from line ~900 on to

        hdf5_flag = 0
        npw = int(kptxml.find("npw").text)
        #        kg= np.random.randint(100,size=(npw,3))-50
        npwtot = npw * (2 if self.spinor else 1)
        CG = np.zeros((IBend - IBstart, npwtot), dtype=complex)
        wfcname="wfc{}{}".format({None:"","dw":"dw","up":"up"}[spin_channel],ik+1)
        try:
            fWFC = h5py.File("{}.save/{}.hdf5".format(prefix,wfcname.upper()),'r')
            hdf5_flag = 1
        except FileNotFoundError:
            fWFC = h5py.File("{}.save/{}.hdf5".format(prefix,wfcname.lower()),'r')
            hdf5_flag = 1
        except FileNotFoundError:
            fWFC=FF("{}.save/{}.dat".format(prefix,wfcname.lower()),"r")
        except FileNotFoundError:
            fWFC=FF("{}.save/{}.dat".format(prefix,wfcname.upper()),"r")

        if hdf5_flag==0: #This is the old code
            rec = record_abinit(fWFC, "i4,3f8,i4,i4,f8")[0]
            ik, xk, ispin, gamma_only, scalef = rec
            #        xk/=bohr
            #        xk=xk.dot(np.linalg.inv(RecLattice))

            rec = record_abinit(fWFC, "4i4")
            #        print ('rec=',rec)
            ngw, igwx, npol, nbnd = rec

            rec = record_abinit(fWFC, "(3,3)f8")
            #        print ('rec=',rec)
            B = np.array(rec)
            #        print (np.mean(B/RecLattice))
            self.K = xk.dot(np.linalg.inv(B))

            rec = record_abinit(fWFC, "({},3)i4".format(igwx))
            #        print ('rec=',rec)
            kg = np.array(rec)
            #        print (np.mean(B/RecLattice))
            #        print ("k-point {0}: {1}/{2}={3}".format(ik, self.K,xk,self.K/xk))
            #        print ("k-point {0}: {1}".format(ik,self.K ))

            for ib in range(IBend):
                cg_tmp = record_abinit(fWFC, "{}f8".format(npwtot * 2))
                if ib >= IBstart:
                    CG[ib - IBstart] = cg_tmp[0::2] + 1.0j * cg_tmp[1::2]

            return sortIG(self.ik0, kg, self.K, CG, B, Ecut0, Ecut, self.spinor)

        elif hdf5_flag==1: #This is my addition

            #This is kind of jank
            attrnames = ['gamma_only', 'igwx', 'ik', 'ispin', 'nbnd', 'ngw', 'npol', 'scale_factor', 'xk']
            attributes =  []
            for atr in attrnames:
                attributes.append(fWFC.attrs[atr])

            gamma_only, igwx, ik, ispin, nbnd, ngw, npol, scale_factor, xk = attributes
            xk = np.array(xk)

            Miller_Indices =  fWFC['MillerIndices']

            B = np.array([Miller_Indices.attrs['bg1'],Miller_Indices.attrs['bg2'],Miller_Indices.attrs['bg3'] ])
            self.K = xk.dot(np.linalg.inv(B))
            kg = np.array(Miller_Indices[::])

            evc = fWFC['evc']
            for ib in range(IBend):
                cg_tmp = evc[ib,:]
                if ib >= IBstart:
                    CG[ib - IBstart] =cg_tmp[0::2] +  1.0j * cg_tmp[1::2]

            return sortIG(self.ik0, kg, self.K, CG, B, Ecut0, Ecut, self.spinor)

then this works for me. I also had to modify an import of irrep.__aux to irrep.utility, not really sure why. Stepan please feel free to add this, or I can try to clean it up myself and submit a merge request.

stepan-tsirkin commented 2 years ago

@MSteggles Thank you for your question and solution. Would you like to make a pull request out of it and become a contributor of IrRep ?

MSteggles commented 2 years ago

Hi Stepan,

Sure! I'll try to test it over the weekend and just make sure it's working.

Cheers, Matthew