Skip to content

Make emc.py compatible with Python3, fix KPOINTS for CASTEP #25

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ganphys wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Make emc.py compatible with Python3, fix KPOINTS for CASTEP #25

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8a48152
Provide Python3 compatibility for emc.py
ganphys Oct 12, 2020
3e35c86
Remove kpoints weights from KPOINTS for CASTEP
ganphys Oct 12, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
108 changes: 57 additions & 51 deletions emc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#!/usr/bin/env python

from __future__ import print_function

EMC_VERSION = '1.51py'
STENCIL = 3 # or 5
#
Expand Down Expand Up @@ -180,7 +182,7 @@ def jacobi(ainput):
tol = max(tolmin,tol*0.99e-2)
#
if nrot != 0:
print 'jacobi: [WARNING] Jacobi iteration did not converge in 50 passes!'
print('jacobi: [WARNING] Jacobi iteration did not converge in 50 passes!')
#
# Sort eigenvectors and values into increasing order
e = [0.0 for i in range(n)] # zerovector
Expand Down Expand Up @@ -211,10 +213,10 @@ def fd_effmass_st3(e, h):
m[2][0] = m[0][2]
m[2][1] = m[1][2]
#
print '-> fd_effmass_st3: Effective mass tensor:\n'
print('-> fd_effmass_st3: Effective mass tensor:\n')
for i in range(len(m)):
print '%15.8f %15.8f %15.8f' % (m[i][0], m[i][1], m[i][2])
print ''
print('%15.8f %15.8f %15.8f' % (m[i][0], m[i][1], m[i][2]))
print('')
#
return m

Expand All @@ -237,10 +239,10 @@ def fd_effmass_st5(e, h):
m[2][0] = m[0][2]
m[2][1] = m[1][2]
#
print '-> fd_effmass_st5: Effective mass tensor:\n'
print('-> fd_effmass_st5: Effective mass tensor:\n')
for i in range(3):
print '%15.8f %15.8f %15.8f' % (m[i][0], m[i][1], m[i][2])
print ''
print('%15.8f %15.8f %15.8f' % (m[i][0], m[i][1], m[i][2]))
print('')
#
return m

Expand All @@ -252,7 +254,7 @@ def generate_kpoints(kpt_frac, st, h, prg, basis):
basis_r = [[ m[i][j]*2.0*pi for j in range(3) ] for i in range(3) ]
#
kpt_rec = MAT_m_VEC(T(basis_r), kpt_frac)
print '-> generate_kpoints: K-point in reciprocal coordinates: %5.3f %5.3f %5.3f' % (kpt_rec[0], kpt_rec[1], kpt_rec[2])
print('-> generate_kpoints: K-point in reciprocal coordinates: %5.3f %5.3f %5.3f' % (kpt_rec[0], kpt_rec[1], kpt_rec[2]))
#
if prg == 'V' or prg == 'P':
h = h*(1/Bohr) # [1/A]
Expand Down Expand Up @@ -310,9 +312,9 @@ def parse_EIGENVAL_VASP(eigenval_fh, band, diff2_size, debug=False):
eigenval_fh.readline()
#
nelec, nkpt, nband = [int(s) for s in eigenval_fh.readline().split()]
if debug: print 'From EIGENVAL: Number of the valence band is %d (NELECT/2)' % (nelec/2)
if debug: print('From EIGENVAL: Number of the valence band is %d (NELECT/2)' % (nelec/2))
if band > nband:
print 'Requested band (%d) is larger than total number of the calculated bands (%d)!' % (band, nband)
print('Requested band (%d) is larger than total number of the calculated bands (%d)!' % (band, nband))
sys.exit(1)

energies = []
Expand All @@ -324,7 +326,7 @@ def parse_EIGENVAL_VASP(eigenval_fh, band, diff2_size, debug=False):
if band == j:
energies.append(float(line.split()[1])*ev2h)

if debug: print ''
if debug: print('')
return energies
#
def parse_nscf_PWSCF(eigenval_fh, band, diff2_size, debug=False):
Expand Down Expand Up @@ -357,11 +359,11 @@ def parse_nscf_PWSCF(eigenval_fh, band, diff2_size, debug=False):
#
a.extend(line.strip().split())
#
#print a
#print(a)
assert len(a) <= band, 'Length of the energies array at a k-point is smaller than band param'
energies.append(float(a[band-1])*ev2h)
#
#print engrs_at_k
#print(engrs_at_k)
return energies
#
def parse_inpcar(inpcar_fh, debug=False):
Expand All @@ -378,33 +380,33 @@ def parse_inpcar(inpcar_fh, debug=False):
p = re.search(r'^\s*(-*\d+\.\d+)\s+(-*\d+\.\d+)\s+(-*\d+\.\d+)', inpcar_fh.readline())
if p:
kpt = [float(p.group(1)), float(p.group(2)), float(p.group(3))]
if debug: print "Found k point in the reduced reciprocal space: %5.3f %5.3f %5.3f" % (kpt[0], kpt[1], kpt[2])
if debug: print("Found k point in the reduced reciprocal space: %5.3f %5.3f %5.3f" % (kpt[0], kpt[1], kpt[2]))
else:
print "Was expecting k point on the line 0 (3 floats), didn't get it, exiting..."
print("Was expecting k point on the line 0 (3 floats), didn't get it, exiting...")
sys.exit(1)

p = re.search(r'^\s*(\d+\.\d+)', inpcar_fh.readline())
if p:
stepsize = float(p.group(1))
if debug: print "Found stepsize of: %5.3f (1/Bohr)" % stepsize
if debug: print("Found stepsize of: %5.3f (1/Bohr)" % stepsize)
else:
print "Was expecting a stepsize on line 1 (1 float), didn't get it, exiting..."
print("Was expecting a stepsize on line 1 (1 float), didn't get it, exiting...")
sys.exit(1)

p = re.search(r'^\s*(\d+)', inpcar_fh.readline())
if p:
band = int(p.group(1))
if debug: print "Requested band is : %5d" % band
if debug: print("Requested band is : %5d" % band)
else:
print "Was expecting band number on line 2 (1 int), didn't get it, exiting..."
print("Was expecting band number on line 2 (1 int), didn't get it, exiting...")
sys.exit(1)

p = re.search(r'^\s*(\w)', inpcar_fh.readline())
if p:
prg = p.group(1)
if debug: print "Program identifier is: %5c" % prg
if debug: print("Program identifier is: %5c" % prg)
else:
print "Was expecting program identifier on line 3 (1 char), didn't get it, exiting..."
print("Was expecting program identifier on line 3 (1 char), didn't get it, exiting...")
sys.exit(1)

for i in range(3):
Expand All @@ -413,11 +415,11 @@ def parse_inpcar(inpcar_fh, debug=False):
basis.append([float(p.group(1)), float(p.group(2)), float(p.group(3))])

if debug:
print "Real space basis:"
print("Real space basis:")
for i in range(len(basis)):
print '%9.7f %9.7f %9.7f' % (basis[i][0], basis[i][1], basis[i][2])
print('%9.7f %9.7f %9.7f' % (basis[i][0], basis[i][1], basis[i][2]))

if debug: print ''
if debug: print('')

return kpt, stepsize, band, prg, basis

Expand All @@ -443,7 +445,7 @@ def get_eff_masses(m, basis):
import datetime
import time
filename = 'emcpy.out_'+str(int(time.time()))
print 'Redirecting output to '+filename
print('Redirecting output to '+filename)
sys.stdout = open(filename, 'w')
#
if STENCIL == 3:
Expand All @@ -453,19 +455,19 @@ def get_eff_masses(m, basis):
fd_effmass = fd_effmass_st5
st = st5
else:
print 'main: [ERROR] Wrong value for STENCIL, should be 3 or 5.'
print('main: [ERROR] Wrong value for STENCIL, should be 3 or 5.')
sys.exit(1)
#
print 'Effective mass calculator '+EMC_VERSION
print 'Stencil: '+str(STENCIL)
print 'License: MIT'
print 'Developed by: Alexandr Fonari and Chris Sutton'
print 'Started at: '+datetime.datetime.now().strftime("%Y-%m-%d %H:%M")+'\n'
print('Effective mass calculator '+EMC_VERSION)
print('Stencil: '+str(STENCIL))
print('License: MIT')
print('Developed by: Alexandr Fonari and Chris Sutton')
print('Started at: '+datetime.datetime.now().strftime("%Y-%m-%d %H:%M")+'\n')
#
if len(sys.argv) == 1:
print "Run as:"
print " %s input.in [output.out]" % sys.argv[0]
print ""
print("Run as:")
print(" %s input.in [output.out]" % sys.argv[0])
print("")
sys.exit(1)
inpcar_fn = sys.argv[1]
#
Expand All @@ -474,11 +476,11 @@ def get_eff_masses(m, basis):
except IOError:
sys.exit("Couldn't open input file "+inpcar_fn+", exiting...\n")
#
print "Contents of the "+inpcar_fn+" file:\n"
print inpcar_fh.read()
print ""
print "=========="
print ""
print("Contents of the "+inpcar_fn+" file:\n")
print(inpcar_fh.read())
print("")
print("==========")
print("")
#
kpt, stepsize, band, prg, basis = parse_inpcar(inpcar_fh)
#
Expand All @@ -491,7 +493,7 @@ def get_eff_masses(m, basis):
sys.exit("Couldn't open input file "+output_fn+", exiting...\n")
#
if output_fn:
print 'Successfully opened '+output_fn+', preparing to parse it...\n'
print('Successfully opened '+output_fn+', preparing to parse it...\n')
#
energies = []
if prg.upper() == 'V' or prg.upper() == 'C':
Expand All @@ -507,19 +509,19 @@ def get_eff_masses(m, basis):
m = fd_effmass(energies, stepsize)
#
masses, vecs_cart, vecs_frac, vecs_n = get_eff_masses(m, basis)
print 'Principle effective masses and directions:\n'
print('Principle effective masses and directions:\n')
for i in range(3):
print 'Effective mass (%d): %12.3f' % (i, masses[i])
print 'Original eigenvectors: %7.5f %7.5f %7.5f' % (vecs_cart[i][0], vecs_cart[i][1], vecs_cart[i][2])
print 'Normal fractional coordinates: %7.5f %7.5f %7.5f\n' % (vecs_n[i][0], vecs_n[i][1], vecs_n[i][2])
print('Effective mass (%d): %12.3f' % (i, masses[i]))
print('Original eigenvectors: %7.5f %7.5f %7.5f' % (vecs_cart[i][0], vecs_cart[i][1], vecs_cart[i][2]))
print('Normal fractional coordinates: %7.5f %7.5f %7.5f\n' % (vecs_n[i][0], vecs_n[i][1], vecs_n[i][2]))
#
else:
print 'No output file provided, entering the Generation regime...\n'
print('No output file provided, entering the Generation regime...\n')
#
if prg.upper() == "C" and band != 1:
print " Band should be set to 1 for CRYSTAL calculations,"
print " desired band number is set as a parameter (-b) for cry-getE.pl script."
print ""
print(" Band should be set to 1 for CRYSTAL calculations,")
print(" desired band number is set as a parameter (-b) for cry-getE.pl script.")
print("")
sys.exit(1)
#
kpoints = generate_kpoints(kpt, st, stepsize, prg, basis)
Expand All @@ -528,8 +530,12 @@ def get_eff_masses(m, basis):
kpoints_fh.write("%d\n" % len(st))
kpoints_fh.write("Reciprocal\n")
#
for i, kpt in enumerate(kpoints):
kpoints_fh.write( '%15.10f %15.10f %15.10f 0.01\n' % (kpt[0], kpt[1], kpt[2]) )
if prg.upper() == 'P':
for i, kpt in enumerate(kpoints):
kpoints_fh.write( '%15.10f %15.10f %15.10f \n' % (kpt[0], kpt[1], kpt[2]) )
else:
for i, kpt in enumerate(kpoints):
kpoints_fh.write( '%15.10f %15.10f %15.10f 0.01\n' % (kpt[0], kpt[1], kpt[2]) )
#
kpoints_fh.close()
print 'KPOINTS file has been generated in the current directory...\n'
print('KPOINTS file has been generated in the current directory...\n')