mirror of https://github.com/QMCPACK/qmcpack.git
163 lines
4.3 KiB
Python
Executable File
163 lines
4.3 KiB
Python
Executable File
#! /usr/bin/env python3
|
|
|
|
import os
|
|
import sys
|
|
from optparse import OptionParser
|
|
|
|
|
|
|
|
def find_nexus_modules():
|
|
import sys
|
|
nexus_lib = os.path.abspath(os.path.join(__file__,'..','..','lib'))
|
|
assert(os.path.exists(nexus_lib))
|
|
sys.path.append(nexus_lib)
|
|
#end def find_nexus_modules
|
|
|
|
|
|
def import_nexus_module(module_name):
|
|
import importlib
|
|
return importlib.import_module(module_name)
|
|
#end def import_nexus_module
|
|
|
|
|
|
# Load Nexus modules
|
|
try:
|
|
# Attempt specialized path-based imports.
|
|
# (The executable should still work even if Nexus is not installed)
|
|
find_nexus_modules()
|
|
|
|
generic = import_nexus_module('generic')
|
|
obj = generic.obj
|
|
del generic
|
|
|
|
developer = import_nexus_module('developer')
|
|
DevBase = developer.DevBase
|
|
error = developer.error
|
|
ci = developer.ci
|
|
del developer
|
|
|
|
hdfreader = import_nexus_module('hdfreader')
|
|
read_hdf = hdfreader.read_hdf
|
|
del hdfreader
|
|
except:
|
|
from generic import obj
|
|
from developer import DevBase,error,ci
|
|
from hdfreader import read_hdf
|
|
#end try
|
|
|
|
|
|
|
|
|
|
def read_eshdf_nofk_data(filename,Ef):
|
|
from numpy import array,pi,dot,sqrt,abs,zeros
|
|
from numpy.linalg import inv,det
|
|
|
|
def h5int(i):
|
|
return array(i,dtype=int)[0]
|
|
#end def h5int
|
|
|
|
E_fermi = Ef + 1e-8
|
|
h = read_hdf(filename,view=True)
|
|
gvu = array(h.electrons.kpoint_0.gvectors)
|
|
axes = array(h.supercell.primitive_vectors)
|
|
kaxes = 2*pi*inv(axes).T
|
|
gv = dot(gvu,kaxes)
|
|
Ngv = len(gv[:,0])
|
|
kmag = sqrt((gv**2).sum(1))
|
|
nk = h5int(h.electrons.number_of_kpoints)
|
|
ns = h5int(h.electrons.number_of_spins)
|
|
occpaths = obj()
|
|
data = obj()
|
|
for k in range(nk):
|
|
kin_k = obj()
|
|
eig_k = obj()
|
|
k_k = obj()
|
|
nk_k = obj()
|
|
nelec_k = zeros((ns,),dtype=float)
|
|
kp = h.electrons['kpoint_'+str(k)]
|
|
gvs = dot(array(kp.reduced_k),kaxes)
|
|
gvk = gv.copy()
|
|
for d in range(3):
|
|
gvk[:,d] += gvs[d]
|
|
#end for
|
|
kinetic=(gvk**2).sum(1)/2 # Hartree units
|
|
for s in range(ns):
|
|
#print ' ',(k,s),(nk,ns)
|
|
kin_s = []
|
|
eig_s = []
|
|
k_s = gvk
|
|
nk_s = 0*kmag
|
|
nelec_s = 0
|
|
path = 'electrons/kpoint_{0}/spin_{1}'.format(k,s)
|
|
spin = h.get_path(path)
|
|
eig = convert(array(spin.eigenvalues),'Ha','eV')
|
|
nst = h5int(spin.number_of_states)
|
|
for st in range(nst):
|
|
e = eig[st]
|
|
if e<E_fermi:
|
|
stpath = path+'/state_{0}/psi_g'.format(st)
|
|
occpaths.append(stpath)
|
|
psi = array(h.get_path(stpath))
|
|
nk_orb = (psi**2).sum(1)
|
|
kin_orb = (kinetic*nk_orb).sum()
|
|
nelec_s += nk_orb.sum()
|
|
nk_s += nk_orb
|
|
kin_s.append(kin_orb)
|
|
eig_s.append(e)
|
|
#end if
|
|
#end for
|
|
data[k,s] = obj(
|
|
kpoint = array(kp.reduced_k),
|
|
kin = array(kin_s),
|
|
eig = array(eig_s),
|
|
k = k_s,
|
|
nk = nk_s,
|
|
ne = nelec_s,
|
|
)
|
|
#end for
|
|
#end for
|
|
res = obj(
|
|
orbfile = filename,
|
|
E_fermi = E_fermi,
|
|
axes = axes,
|
|
kaxes = kaxes,
|
|
nkpoints = nk,
|
|
nspin = ns,
|
|
data = data,
|
|
)
|
|
return res
|
|
#end def read_eshdf_nofk_data
|
|
|
|
|
|
|
|
def kinetic():
|
|
# read command line inputs
|
|
usage = '''usage: %prog kinetic [options] [eshdf_file]'''
|
|
parser = OptionParser(usage=usage,add_help_option=False,version='%prog 0.1')
|
|
|
|
options,args = parser.parse_args()
|
|
scalar_files = list(sorted(args[1:]))
|
|
#end def kinetic
|
|
|
|
|
|
|
|
operations = obj(
|
|
kinetic = kinetic,
|
|
)
|
|
|
|
|
|
|
|
if __name__=='__main__':
|
|
|
|
if len(sys.argv)<2:
|
|
error('First argument must be operation to perform on ESHDF file.\ne.g. to examine kinetic energies, type "eshdf kinetic ..."\nValid operations are: {0}'.format(sorted(operations.keys())))
|
|
#end if
|
|
op_type = sys.argv[1]
|
|
if op_type in operations:
|
|
operations[op_type]()
|
|
else:
|
|
error('Unknown operation: {0}\nValid options are: {1}'.format(op_type,sorted(operations.keys())))
|
|
#end if
|
|
|
|
#end if
|