From 37b60c01497fce2396bff3e7a70b57535f6e2490 Mon Sep 17 00:00:00 2001
From: Sourcery AI <>
Date: Sun, 17 Jul 2022 14:42:00 +0000
Subject: [PATCH] 'Refactored by Sourcery'
---
Exciting_Xml2Vasp.py | 212 ++++++------
Extract_elastic_energy_VASP.py | 306 +++++++++---------
.../Extract_hessian.py | 62 ++--
.../parsing hessian matrix | 17 +-
get_forces_all_ionicSteps.py | 223 ++++++-------
src/Main_file.py | 37 ++-
src/contcar_poscar.py | 259 ++++++++-------
src/create_phonon_directories.py | 36 ++-
src/elastic_constants.py | 101 +++---
src/energy.py | 59 ++--
src/evaluate_manually_Elastic_constants.py | 130 ++++----
src/fit_energy_vs_vol.py | 154 +++++----
src/poscar.py | 157 +++++----
13 files changed, 866 insertions(+), 887 deletions(-)
diff --git a/Exciting_Xml2Vasp.py b/Exciting_Xml2Vasp.py
index d95afe4..95ccede 100644
--- a/Exciting_Xml2Vasp.py
+++ b/Exciting_Xml2Vasp.py
@@ -33,8 +33,8 @@
#*********************DEFINITION FOR INPUTFILE**********************************
if len(sys.argv) < 2:
- sys.exit('Usage: %s <Exciting xml inputfile>' % sys.argv[0])
-
+ sys.exit(f'Usage: {sys.argv[0]} <Exciting xml inputfile>')
+
tree = xml.parse(sys.argv[1])
root = tree.getroot()
@@ -42,122 +42,108 @@
lattice = struct.getchildren()
-#*********************CONTROL SECTION FOR INPUT FILE***************************
-
-f = open('POSCAR.vasp', 'w+')
-
-print (colored("\nExtracting information from input file... ", 'yellow'))
-print ('-'*80)
-print (colored(" | Title is: ", 'yellow'), root.findtext('title'))
-f.write(root.findtext('title'))
+with open('POSCAR.vasp', 'w+') as f:
+ print (colored("\nExtracting information from input file... ", 'yellow'))
+ print ('-'*80)
+ print (colored(" | Title is: ", 'yellow'), root.findtext('title'))
+ f.write(root.findtext('title'))
-#*************************DETECTING # OF SPECIES******************************
+ #*************************DETECTING # OF SPECIES******************************
-i = 0
-for child in struct:
- i = i+1
-print (colored(" | Species detected: ", 'yellow'), i-1 )
-f.write(" |\t Number of species detected: {:d}\n".format(i-1))
+ i = 0
+ for _ in struct:
+ i = i+1
+ print (colored(" | Species detected: ", 'yellow'), i-1 )
+ f.write(" |\t Number of species detected: {:d}\n".format(i-1))
-#***************************SCALING DEFINITION********************************
-#reading from a file
-for sca in struct.findall('crystal'):
- for key in sca.attrib:
- if key == 'scale': # SCALING HAS TO BE CHANGED
- sc = sca.get('scale')
+ #***************************SCALING DEFINITION********************************
+ #reading from a file
+ for sca in struct.findall('crystal'):
+ for key in sca.attrib:
+ sc = sca.get('scale') if key == 'scale' else 1.0
print (colored(" | Scaling factor is: ", 'yellow'), float(sc))
- else:
- sc = float(1.00)
- print (colored(" | Scaling factor is: ", 'yellow'), float(sc))
-
-
-print (colored(" | Converting bohr to angstrom ... ", 'yellow'))
-print (colored(" | Writing to a File in VASP format ... ", 'yellow'))
-f.write("%f\n" % 1.0)
-print ('-'*80)
-
-#***********************LATTICE VECTORS SECTION************************
-#reading from a file
-lv = []
-for vect in struct.findall('crystal'):
- #lat = vect.find('basevect')
- for k in vect.findall('basevect'):
- lv.append(k.text.split())
-###
-
-convertion = Bohr2angstrom*float(sc)
-for list in lv:
- print ("\t{:10.8f} {:10.8f} {:10.8f}".format(float(list[0])*convertion,
- float(list[1])*convertion, float(list[2])*convertion ) )
- f.write ("\t{:10.8f} {:10.8f} {:10.8f}\n".format(float(list[0])*convertion,
- float(list[1])*convertion, float(list[2])*convertion ) )
-
-#******************SPECIES CONTROL SECTION*********************************
-
-lab = []
-for coord in struct.findall('species'):
- spe = coord.get('speciesfile')
- jj = os.path.splitext(spe)
- lab.append(jj[0])
-
-sp = []
-for x in range(i-1):
- print("{:s}".format(lab[x]))
- f.write("\t%s" % lab[x])
- for s in struct[x+1].iter('atom'):
- sp.append(lab[x])
-f.write("\n")
-
-for h in range(i-1):
- print("\t{:d}".format(sp.count(lab[h])) )
- f.write("\t%d" % sp.count(lab[h]))
-f.write("\n")
-
-#*******************CONTROL SECTION FOR CARTESIAN COORDINATE***************
-
-lll = []
-ccl = []
-checker = False
-for cart in root.findall('structure'):
- ca = cart.get('cartesian')
- for key in cart.attrib:
- if key == 'cartesian':
- checker = True
- if ca == "true":
- print ("Cartesian")
- f.write("Cartesian\n")
- for i in range(i-1):
- for s in struct[i+1].iter('atom'):
- atom = s.get('coord')
- ccl.append(atom.split())
- lll.append(lab[i])
+ print (colored(" | Converting bohr to angstrom ... ", 'yellow'))
+ print (colored(" | Writing to a File in VASP format ... ", 'yellow'))
+ f.write("%f\n" % 1.0)
+ print ('-'*80)
+
+ #***********************LATTICE VECTORS SECTION************************
+ #reading from a file
+ lv = []
+ for vect in struct.findall('crystal'):
+ #lat = vect.find('basevect')
+ lv.extend(k.text.split() for k in vect.findall('basevect'))
+ ###
+
+ convertion = Bohr2angstrom*float(sc)
+ for list in lv:
+ print ("\t{:10.8f} {:10.8f} {:10.8f}".format(float(list[0])*convertion,
+ float(list[1])*convertion, float(list[2])*convertion ) )
+ f.write ("\t{:10.8f} {:10.8f} {:10.8f}\n".format(float(list[0])*convertion,
+ float(list[1])*convertion, float(list[2])*convertion ) )
+
+ #******************SPECIES CONTROL SECTION*********************************
+
+ lab = []
+ for coord in struct.findall('species'):
+ spe = coord.get('speciesfile')
+ jj = os.path.splitext(spe)
+ lab.append(jj[0])
+
+ sp = []
+ for x in range(i-1):
+ print("{:s}".format(lab[x]))
+ f.write("\t%s" % lab[x])
+ sp.extend(lab[x] for _ in struct[x+1].iter('atom'))
+ f.write("\n")
+
+ for h in range(i-1):
+ print("\t{:d}".format(sp.count(lab[h])) )
+ f.write("\t%d" % sp.count(lab[h]))
+ f.write("\n")
+
+ #*******************CONTROL SECTION FOR CARTESIAN COORDINATE***************
+
+ lll = []
+ ccl = []
+ checker = False
+ for cart in root.findall('structure'):
+ ca = cart.get('cartesian')
+ for key in cart.attrib:
+ if key == 'cartesian':
+ checker = True
+ if ca == "true":
+ print ("Cartesian")
+ f.write("Cartesian\n")
+ for i in range(i-1):
+ for s in struct[i+1].iter('atom'):
+ atom = s.get('coord')
+ ccl.append(atom.split())
+ lll.append(lab[i])
# print (atom, lab[i])
-
- for list in ccl:
- print ("{:5.8f} {:5.8f} {:5.8f}".format(float(list[0])* Bohr2angstrom,
- float(list[1])* Bohr2angstrom, float(list[2])* Bohr2angstrom ) )
- f.write ("{:5.8f} {:5.8f} {:5.8f}\n".format(float(list[0])* Bohr2angstrom,
- float(list[1])* Bohr2angstrom, float(list[2])* Bohr2angstrom ) )
-
+
+ for list in ccl:
+ print ("{:5.8f} {:5.8f} {:5.8f}".format(float(list[0])* Bohr2angstrom,
+ float(list[1])* Bohr2angstrom, float(list[2])* Bohr2angstrom ) )
+ f.write ("{:5.8f} {:5.8f} {:5.8f}\n".format(float(list[0])* Bohr2angstrom,
+ float(list[1])* Bohr2angstrom, float(list[2])* Bohr2angstrom ) )
#**************************DIRECT COORDINATE*********************
-
- else:
- print ("Direct")
- f.write ("Direct\n")
- for i in range(i-1):
- for s in struct[i+1].iter('atom'):
- atom = s.get('coord')
- ccl.append(atom.split())
- lll.append(lab[i])
+
+ else:
+ print ("Direct")
+ f.write ("Direct\n")
+ for i in range(i-1):
+ for s in struct[i+1].iter('atom'):
+ atom = s.get('coord')
+ ccl.append(atom.split())
+ lll.append(lab[i])
# print (atom, lab[i])
# print (ccl)
-
- for list in ccl:
- print ("{:5.8f} {:5.8f} {:5.8f}".format(float(list[0]),
- float(list[1]), float(list[2]) ) )
- f.write ("{:5.8f} {:5.8f} {:5.8f}\n".format(float(list[0]),
- float(list[1]), float(list[2]) ) )
-
-print (colored("File generated ... ", 'green'))
-
-f.close()
+
+ for list in ccl:
+ print ("{:5.8f} {:5.8f} {:5.8f}".format(float(list[0]),
+ float(list[1]), float(list[2]) ) )
+ f.write ("{:5.8f} {:5.8f} {:5.8f}\n".format(float(list[0]),
+ float(list[1]), float(list[2]) ) )
+
+ print (colored("File generated ... ", 'green'))
diff --git a/Extract_elastic_energy_VASP.py b/Extract_elastic_energy_VASP.py
index d781a80..ae7d586 100644
--- a/Extract_elastic_energy_VASP.py
+++ b/Extract_elastic_energy_VASP.py
@@ -42,9 +42,12 @@ def poscar():
print (' ERROR: POSCAR does not exist here.')
sys.exit(0)
print('Reading POSCAR/CONTCAR: \n')
- pos = []; kk = []; lattice = []; sum = 0
+ pos = []
+ kk = []
+ lattice = []
+ sum = 0
file = open('POSCAR','r') or open('CONTCAR','r')
-
+
firstline = file.readline() # IGNORE first line comment
secondfline = file.readline() # scale
Latvec1 = file.readline()
@@ -56,11 +59,10 @@ def poscar():
elementtype=file.readline().split()
if (str.isdigit(elementtype[0])):
sys.exit("VASP 4.X POSCAR detected. Please add the atom types")
- print ("Types of elements:", str(elementtype), end = '\n')
+ print("Types of elements:", elementtype, end = '\n')
numberofatoms=file.readline()
Coordtype=file.readline()
print ("Coordtype:", (Coordtype), end = '\n')
-
########################---------------------------------------------------------
print (">>>>>>>>>-------------------# of Atoms--------------------")
nat = numberofatoms.split()
@@ -72,37 +74,31 @@ def poscar():
print ("Number of atoms:", (numberofatoms), end = '\n')
########################---------------------------------------------------------
- #print (">>>>>>>>>---------------Atomic positions------------------")
- for x in range(int(numberofatoms)):
+ #print (">>>>>>>>>---------------Atomic positions------------------")
+ for _ in range(int(numberofatoms)):
coord = file.readline().split()
coord = [float(i) for i in coord]
pos = pos + [coord]
pos = np.array(pos)
#print (pos)
-
+
file.close()
-########################---------------------------------------------------------
- a=[]; b=[]; c=[];
Latvec1=Latvec1.split()
Latvec2=Latvec2.split()
- Latvec3=Latvec3.split()
- for ai in Latvec1:
- a.append(float(ai))
- for bi in Latvec2:
- b.append(float(bi))
- for ci in Latvec3:
- c.append(float(ci))
-
+ Latvec3=Latvec3.split()
+ a = [float(ai) for ai in Latvec1]
+ b = [float(bi) for bi in Latvec2]
+ c = [float(ci) for ci in Latvec3]
########################---------------------------------------------------------
- print (">>>>>>>>>---------------Lattice vectors distortions-----------------")
+ print (">>>>>>>>>---------------Lattice vectors distortions-----------------")
lattice = np.array([a] + [b] + [c])
#determinant = np.linalg.det(lattice)
lld = lattic_distortion.local_lattice_distortion(a,b,c)
- print ("lattice distortion parameter g: {}".format(lld) )
-
- print (">>>>>>>>>---------------Space group-----------------")
+ print(f"lattice distortion parameter g: {lld}")
+
+ print (">>>>>>>>>---------------Space group-----------------")
print (" ")
sp, symm = space_group_analyse(lattice, pos)
print ( sp, symm )
@@ -124,17 +120,19 @@ def poscar():
###
class lattic_distortion():
@classmethod
- def local_lattice_distortion(a1,b1,c1):
+ def local_lattice_distortion(cls, b1, c1):
#print ("The lattice distortion in paracrystals is measured by the lattice distortion parameter g")
#print (Back.YELLOW + "Wang, S. Atomic structure modeling of multi-principal-element alloys by the principle")
#print (Back.YELLOW + "of maximum entropy. Entropy 15, 5536–5548 (2013).")
#print ("")
- a=np.linalg.norm(a1); b=np.linalg.norm(b1); c=np.linalg.norm(c1)
+ a = np.linalg.norm(cls)
+ b=np.linalg.norm(b1)
+ c=np.linalg.norm(c1)
d = np.array([a,b,c])
- d_mean = np.mean(d); d_std = np.std(d)
+ d_mean = np.mean(d)
+ d_std = np.std(d)
d_square_mean = (a**2 + b**2 + c**2)/3
- g = np.sqrt( d_square_mean/(d_mean)**2 - 1 )
- return g
+ return np.sqrt( d_square_mean/(d_mean)**2 - 1 )
###
def local_lattice_distortion_DEF1():
@@ -142,7 +140,9 @@ def local_lattice_distortion_DEF1():
#print (Back.YELLOW + "Wang, S. Atomic structure modeling of multi-principal-element alloys by the principle")
#print (Back.YELLOW + "of maximum entropy. Entropy 15, 5536–5548 (2013).")
print ("+"*40,"HUME ROTHERY RULE","+"*40)
- C_i=C=0.2 ; r_avg = 0.0; del_sum=0.0
+ C_i=C=0.2
+ r_avg = 0.0
+ del_sum=0.0
elements = ["Nb", "Hf", "Ta", "Ti", "Zr"]
eta = {
"Nb" : 1.98,
@@ -150,17 +150,17 @@ def local_lattice_distortion_DEF1():
"Ta" : 2.00,
"Ti" : 1.76,
"Zr" : 2.06, }
-
+
print (" {element: atomic radius}")
print (eta)
-
+
for i in elements:
r_avg = r_avg + C * eta[i]
-
+
for j in elements:
del_sum = del_sum + C * ( 1 - float(eta[j]) / r_avg )**2
- del_sum = 100 * np.sqrt(del_sum)
- print("HEA_atomic_size_mismatch: \u03B4={}".format(del_sum))
+ del_sum = 100 * np.sqrt(del_sum)
+ print(f"HEA_atomic_size_mismatch: \u03B4={del_sum}")
###
def local_lattice_distortion_DEF2():
@@ -169,42 +169,43 @@ def local_lattice_distortion_DEF2():
print (" Phys. Rev. Mater. 1, 23404 (2017).")
print (" (***) Different definition of the atomic radius for the description ")
print (" of the local lattice distortion in HEAs")
-
+
if not os.path.exists('POSCAR' and 'CONTCAR'):
print ('>>> ERROR: POSCAR & CONTCAR does not exist (Both should be in the same directory)')
sys.exit(0)
print('Reading POSCAR and CONTCAR ... \n')
-
- x = []; y =[]; z=[]
- xp =[]; yp = []; zp = []; temp=0
-
- f = open('POSCAR','r')
- lines_poscar = f.readlines()
- f.close()
-
- f = open('CONTCAR','r')
- lines_contcar = f.readlines()
- f.close()
-
+
+ x = []
+ y =[]
+ z=[]
+ xp =[]
+ yp = []
+ zp = []
+ temp=0
+
+ with open('POSCAR','r') as f:
+ lines_poscar = f.readlines()
+ with open('CONTCAR','r') as f:
+ lines_contcar = f.readlines()
file_P = ase.io.read('POSCAR')
pos = file_P.get_cell_lengths_and_angles()
- print (CRED + "POSCAR=>Length&Angles->{}".format(pos) + CEND)
+ print(CRED + f"POSCAR=>Length&Angles->{pos}" + CEND)
file_C = ase.io.read('CONTCAR')
- con = file_C.get_cell_lengths_and_angles()
- print (CRED + "CONTCAR=>Length&Angles->{}".format(con) + CEND)
+ con = file_C.get_cell_lengths_and_angles()
+ print(CRED + f"CONTCAR=>Length&Angles->{con}" + CEND)
print ("Cell vectors difference:: ",con-pos)
-
+
sum_atoms = lines_poscar[6].split() ### reading 7th lines for reading # of atoms
sum_atoms = [int(i) for i in sum_atoms]
sum_atoms = sum(sum_atoms)
-
+
for i in lines_poscar:
if "Direct" in i:
lp=lines_poscar.index(i)
for j in lines_contcar:
if "Direct" in j:
lc=lines_contcar.index(j)
-
+
for i in range(sum_atoms):
x, y, z = lines_poscar[lp+1+i].split()
xc, yc, zc = lines_contcar[lp+1+i].split()
@@ -212,7 +213,7 @@ def local_lattice_distortion_DEF2():
xc = float(xc); yc = float(yc); zc = float(zc)
temp = temp + np.sqrt( (x-xc)**2 + (y-yc)**2 + (z-zc)**2 )
temp = temp/sum_atoms
- print("local lattice distortion (LLD): \u0394d={}".format(temp))
+ print(f"local lattice distortion (LLD): \u0394d={temp}")
###
def space_group_analyse(lattice, pos):
@@ -250,50 +251,44 @@ def space_group_analyse(lattice, pos):
def poscar_VASP42VASP5():
if not os.path.exists('POSCAR' and 'POTCAR'):
print (' ERROR: POSCAR does not exist here.')
- sys.exit(0)
- file1 = open("POSCAR",'r')
- line1 = file1.readlines()
- file1.close()
-
- file2 = open("POTCAR",'r')
- line2 = file2.readlines()
- file2.close()
-
+ sys.exit(0)
+ with open("POSCAR",'r') as file1:
+ line1 = file1.readlines()
+ with open("POTCAR",'r') as file2:
+ line2 = file2.readlines()
atom_number=[]
for i in line1:
if ("Direct" or "direct" or "d" or "D") in i:
PP=line1.index(i)
atom_number = line1[5].split()
print(atom_number)
-
- elementtype=[]; count=0
+
+ elementtype=[]
+ count=0
for i in line2:
if ("VRHFIN") in i:
count+=1
#print (i.split('=')[1].split(':')[0])
elementtype.append(i.split('=')[1].split(':')[0])
-
- test = open("POSCAR_W",'w')
-
- for i in range(5):
- test.write( line1[i] )
-
- for j in elementtype:
- test.write("\t" + j)
- test.write("\n" )
-
- for j in atom_number:
- test.write("\t" + j )
- test.write("\n" )
-
- test.write("Selective dynamics")
- test.write("\n" )
-
- for i in range(len(line1)-PP):
- test.write(line1[PP+i] )
-
- test.close()
-
+
+ with open("POSCAR_W",'w') as test:
+ for i in range(5):
+ test.write( line1[i] )
+
+ for j in elementtype:
+ test.write("\t" + j)
+ test.write("\n" )
+
+ for j in atom_number:
+ test.write("\t" + j )
+ test.write("\n" )
+
+ test.write("Selective dynamics")
+ test.write("\n" )
+
+ for i in range(len(line1)-PP):
+ test.write(line1[PP+i] )
+
print (" File is converted: POSCAR_W")
###
@@ -306,11 +301,14 @@ def poscar_VASP42VASP5():
def main_poscar():
count = 0
os.system("rm out_POSCARS.dat")
- VOL_P = []; pos = []; kk = []; lattice = [];
+ VOL_P = []
+ pos = []
+ kk = []
+ lattice = [];
mypath = os.getcwd()
- print ("-"*100)
+ print ("-"*100)
print (Back.YELLOW + "{:15s} {:15s} {:15.6s} {:15.6s} {:15.6s} {:15.15s}".format("Directory", "# of atoms", "||a||", \
- "||b||", "||c||", "VOL_POS[A^3]"),end="\n" )
+ "||b||", "||c||", "VOL_POS[A^3]"),end="\n" )
print ("-"*100)
for entry in os.listdir(mypath):
@@ -318,101 +316,93 @@ def main_poscar():
#print (entry)
for file in os.listdir(entry):
if file == "POSCAR":
- count+=1; sum = 0
+ count+=1
+ sum = 0
filepath = os.path.join(entry, file)
- #f = open(filepath, 'r')
- #print (f.read())
- #f.close()
- fo = open(filepath, 'r')
- ofile=open('out_POSCARS.dat','a')
-
- #print (colored('>>>>>>>> Name of the file: ','red'), fo.name, end = '\n', flush=True)
-
- ofile.write (fo.name + '\n')
- ofile.write ("")
- firstline = fo.readline()
- secondfline = fo.readline()
- Latvec1 = fo.readline()
- #print ("Lattice vector 1:", (Latvec1), end = '')
- #ofile.write (Latvec1)
- Latvec2 = fo.readline()
- #print ("Lattice vector 2:", (Latvec2), end = '')
- #ofile.write (Latvec2)
- Latvec3 = fo.readline()
- #print ("Lattice vector 3:", (Latvec3), end = '')
- #ofile.write (Latvec3)
- elementtype=fo.readline()
- elementtype = elementtype.split()
- #print ("Types of elements:", str(elementtype), end = '')
- #ofile.write (str(elementtype))
- numberofatoms=fo.readline()
- #print ("Number of atoms:", (numberofatoms), end = '')
- #ofile.write ((numberofatoms))
- Coordtype=fo.readline()
-
+ with open(filepath, 'r') as fo:
+ ofile=open('out_POSCARS.dat','a')
+
+ #print (colored('>>>>>>>> Name of the file: ','red'), fo.name, end = '\n', flush=True)
+
+ ofile.write (fo.name + '\n')
+ ofile.write ("")
+ firstline = fo.readline()
+ secondfline = fo.readline()
+ Latvec1 = fo.readline()
+ #print ("Lattice vector 1:", (Latvec1), end = '')
+ #ofile.write (Latvec1)
+ Latvec2 = fo.readline()
+ #print ("Lattice vector 2:", (Latvec2), end = '')
+ #ofile.write (Latvec2)
+ Latvec3 = fo.readline()
+ #print ("Lattice vector 3:", (Latvec3), end = '')
+ #ofile.write (Latvec3)
+ elementtype=fo.readline()
+ elementtype = elementtype.split()
+ #print ("Types of elements:", str(elementtype), end = '')
+ #ofile.write (str(elementtype))
+ numberofatoms=fo.readline()
+ #print ("Number of atoms:", (numberofatoms), end = '')
+ #ofile.write ((numberofatoms))
+ Coordtype=fo.readline()
##########################---------------------------------------------------------
- #print ("**********-------------------# of Atoms--------------------")
-
- nat = numberofatoms.split()
- nat = [int(i) for i in nat]
- for i in nat:
- sum = sum + i
- numberofatoms = sum
- #print ("{} :: Number of atoms: {}".format(nat, numberofatoms) )
+ #print ("**********-------------------# of Atoms--------------------")
+
+ nat = numberofatoms.split()
+ nat = [int(i) for i in nat]
+ for i in nat:
+ sum = sum + i
+ numberofatoms = sum
+ #print ("{} :: Number of atoms: {}".format(nat, numberofatoms) )
##########################---------------------------------------------------------
- #print ("-----------------------Atomic positions-----------------")
- #print ("Coordtype:", (Coordtype), end = '')
- for x in range(int(numberofatoms)):
- coord = fo.readline().split()
- coord = [float(i) for i in coord]
- pos = pos + [coord]
- pos = np.array(pos)
- #print (pos)
-
- ofile.write ("\n")
- fo.close()
-##########################---------------------------------------------------------
-
- a=[]; b=[]; c=[];
+ #print ("-----------------------Atomic positions-----------------")
+ #print ("Coordtype:", (Coordtype), end = '')
+ for _ in range(int(numberofatoms)):
+ coord = fo.readline().split()
+ coord = [float(i) for i in coord]
+ pos = pos + [coord]
+ pos = np.array(pos)
+ #print (pos)
+
+ ofile.write ("\n")
Latvec1=Latvec1.split()
Latvec2=Latvec2.split()
Latvec3=Latvec3.split()
-
-##########################---------------------------------------------------------
- for ai in Latvec1: a.append(float(ai))
- for bi in Latvec2: b.append(float(bi))
- for ci in Latvec3: c.append(float(ci))
+
+ a = [float(ai) for ai in Latvec1]
+ b = [float(bi) for bi in Latvec2]
+ c = [float(ci) for ci in Latvec3]
#print ('a=', a)
- ofile.write ("'a=' {}\n".format(a))
+ ofile.write(f"'a=' {a}\n")
#print ('b=', b)
- ofile.write ("'b=' {}\n".format(b))
+ ofile.write(f"'b=' {b}\n")
#print ('c=', c)
- ofile.write ("'c=' {}\n".format(c))
+ ofile.write(f"'c=' {c}\n")
lld = lattic_distortion.local_lattice_distortion(a,b,c)
##########################---------------------------------------------------------
-
+
alpha, beta, gamma = lattice_angles(a,b,c)
- VOL_POS = np.dot(a, np.cross(b,c))
+ VOL_POS = np.dot(a, np.cross(b,c))
VOL_P.append(VOL_POS)
- ofile.write ("'\u03B1=' {} '\u03B2=' {} '\u03B3=' {}\n".format(alpha,beta,gamma))
- ofile.write ("'||a||=' {}\n".format(np.linalg.norm(a)))
- ofile.write ("'||b||=' {}\n".format(np.linalg.norm(b)))
- ofile.write ("'||c||=' {}\n".format(np.linalg.norm(c)))
+ ofile.write(f"'\u03B1=' {alpha} '\u03B2=' {beta} '\u03B3=' {gamma}\n")
+ ofile.write(f"'||a||=' {np.linalg.norm(a)}\n")
+ ofile.write(f"'||b||=' {np.linalg.norm(b)}\n")
+ ofile.write(f"'||c||=' {np.linalg.norm(c)}\n")
#print ("#####------------------------------------------------")
-
+
#print ("a={} \t ||a||={:10.6f}".format(a, np.linalg.norm(a)) )
#print ("b={} \t ||b||={:10.6f}".format(b, np.linalg.norm(b)) )
#print ("c={} \t ||c||={:10.6f}".format(c, np.linalg.norm(c)) )
print ("{:15s} {:6d} {:15.6f} {:15.6f} {:15.6f} {:15.6f}".format(fo.name, numberofatoms, np.linalg.norm(a), \
np.linalg.norm(b), np.linalg.norm(c), VOL_POS) )
-
+
print ("'\u03B1=' {:6.6f} '\u03B2=' {:6.6f} '\u03B3=' {:6.6f} g={:6.6f}".format(alpha,beta,gamma,lld))
print ("."*5)
#print ('Vol= {:6.6f} A^3'.format(VOL_POS))
ofile.write ("***************************************************\n")
ofile.close()
- print ("_"*30)
+ print ("_"*30)
print ("Number of folders detected: ", count)
return VOL_P
diff --git a/Parsing-Hessian-matrix-from-vasprun/Extract_hessian.py b/Parsing-Hessian-matrix-from-vasprun/Extract_hessian.py
index e96b86d..ac0d825 100644
--- a/Parsing-Hessian-matrix-from-vasprun/Extract_hessian.py
+++ b/Parsing-Hessian-matrix-from-vasprun/Extract_hessian.py
@@ -68,31 +68,25 @@
# H = (3N,3N)
#-----------------------------------------------------------------------------------
-hessian = root.xpath('*/dynmat/varray/v/text()')
-hess_v = []
+hessian = root.xpath('*/dynmat/varray/v/text()')
+hess_v = [ind_hess.split() for ind_hess in hessian]
-for ind_hess in hessian:
- hess_v.append( ind_hess.split() )
hess_v = np.array(hess_v)
print("*"*80)
#print(hess_v) # for debugging ONLY
-row = len(hess_v); print("# of rows:", row)
-col = len(hess_v[0]); print("# of columns:", col)
+row = len(hess_v)
+print("# of rows:", row)
+col = len(hess_v[0])
+print("# of columns:", col)
print("*"*80)
-#-----------------------------------------------------------------------------------
-
-#*********************** writing Hessian to a file ********************************
-out = open('hessian.dat','w')
-
-for i in range( int(row/2) ):
- for j in range(col):
- #print(hess_v[i][j], end=' ' )
- out.write( "{:18.11s}".format( hess_v[i][j]) )
- out.write('\n')
- #print()
-
-out.close()
+with open('hessian.dat','w') as out:
+ for i in range(row // 2):
+ for j in range(col):
+ #print(hess_v[i][j], end=' ' )
+ out.write( "{:18.11s}".format( hess_v[i][j]) )
+ out.write('\n')
+ #print()
#******************************End Hessian******************************************
@@ -112,24 +106,22 @@
#******************************End basevectors********************************
'''
#---------------------------------'''METHOD 2'''---------------------------------
-
-a = [] ; l=0
-f = open('hessian.dat', 'r')
-
-for line in f.readlines():
- a.append([])
- #print (line[0])
- if (line[0] == '#'):
- continue
- for i in line.strip().split():
- a[-1].append(float(i))
- l+=1
-#print (a)
-
-f.close()
+
+a = []
+l=0
+with open('hessian.dat', 'r') as f:
+ for line in f:
+ a.append([])
+ #print (line[0])
+ if (line[0] == '#'):
+ continue
+ for i in line.strip().split():
+ a[-1].append(float(i))
+ l+=1
#------------------------------------------------------------------------
#a = np.random.random((900, 900))
-s = np.shape(a); print('size of a Matrix', s)
+s = np.shape(a)
+print('size of a Matrix', s)
a = np.matrix(a)
#print(a)
print ('--------------------------------EIGVAL AND EIGVECTORS ARE:')
diff --git a/Parsing-Hessian-matrix-from-vasprun/parsing hessian matrix b/Parsing-Hessian-matrix-from-vasprun/parsing hessian matrix
index 5372676..8bc9422 100644
--- a/Parsing-Hessian-matrix-from-vasprun/parsing hessian matrix
+++ b/Parsing-Hessian-matrix-from-vasprun/parsing hessian matrix
@@ -25,18 +25,15 @@ subprocess.call(['sed','-i','/^[[:space:]]*$/d',sys.argv[1]], shell = False)
#---------------------------------'''METHOD 2'''
a = []
-f = open(filename, 'r')
-
-for line in f.readlines():
- a.append([])
- for i in line.strip().split():
- a[-1].append(float(i))
-#print (a)
-
-f.close()
+with open(filename, 'r') as f:
+ for line in f:
+ a.append([])
+ for i in line.strip().split():
+ a[-1].append(float(i))
#------------------------------------------------------------------------
#a = np.random.random((900, 900))
-s = np.shape(a); print('size of a Matrix', s)
+s = np.shape(a)
+print('size of a Matrix', s)
a = np.matrix(a)
#print(a)
diff --git a/get_forces_all_ionicSteps.py b/get_forces_all_ionicSteps.py
index a67d817..e3e5800 100644
--- a/get_forces_all_ionicSteps.py
+++ b/get_forces_all_ionicSteps.py
@@ -32,115 +32,122 @@ def get_ediff(out_car):
# INITIALIZATION OF THE DATABASE
ini_path = os.path.join(os.getcwd(), 'Arg_2_metal')
-cwd = os.path.abspath(os.path.dirname(__file__))
+cwd = os.path.abspath(os.path.dirname(__file__))
outdbpath = os.path.join(cwd, '2M_forces_database.db')
# ITERATE OVER ALL 2M-dir's and get the forces at each ionic steps
for j, d in enumerate(tqdm(next(os.walk(ini_path))[1])):
- outcarfile = os.path.join(ini_path, d, '03_ISIF3snsb', 'OUTCAR')
- print(j, outcarfile)
- outcar = open(outcarfile, 'r')
- outcarlines = outcar.readlines()
-
- #Find max iterations
- nelmax = int(subprocess.Popen('grep "NELM " ' + outcarfile, stdout=subprocess.PIPE, shell=True).communicate()[0].split()[2][0:-1])
- nsw = int(subprocess.Popen('grep "NSW" ' + outcarfile + ' | head -1 | cut -d"=" -f2', stdout=subprocess.PIPE, shell=True).communicate()[0] )
- natoms = get_number_of_atoms(outcarfile)
- ediff = np.log10(float(get_ediff(outcarfile)))
-
- re_energy = re.compile("free energy")
- re_iteration = re.compile("Iteration")
- re_timing = re.compile("LOOP:")
- re_force = re.compile("TOTAL-FORCE")
- re_mag = re.compile("number of electron")
- re_volume = re.compile("volume of cell")
-
- lastenergy = 0.0
- energy = 0.0
- steps = 1
- iterations = 0
- cputime = 0.0
- totaltime = 0.0
- dE = 0.0
- magmom = 0.0
- volume = 0.0
- spinpolarized = False
- #average = 0.0
- #maxforce = 0.0
-
- i = 0
- for line in outcarlines:
-
- if re_iteration.search(line):
- iterations = iterations + 1
-
- if re_force.search(line):
- # Calculate forces here...
- forces = []
- magnitudes = []
- for j in range(0,natoms):
- parts = outcarlines[i+j+2].split()
- x = float(parts[3])
- y = float(parts[4])
- z = float(parts[5])
- forces.append([x,y,z])
- magnitudes.append(np.sqrt(x*x + y*y + z*z))
-
- average = sum(magnitudes)/natoms
- maxforce = max(magnitudes)
-
- if re_mag.search(line):
- parts = line.split()
- if len(parts) > 5 and parts[0].strip() != "NELECT":
- spinpolarized = True
- magmom = float(parts[5])
-
- if re_timing.search(line):
- cputime = cputime + float(line.split()[6])/60.0
-
- if re_volume.search(line):
- parts = line.split()
- if len(parts) > 4:
- volume = float(parts[4])
-
- if re_energy.search(line):
- lastenergy = energy
- energy = float(line.split()[4])
- dE = np.log10(abs(energy-lastenergy+1.0E-12))
-
- # Construct output string
- try:
- stepstr = str(steps).rjust(4)
- energystr = "Energy: " + ("%3.6f" % (energy)).rjust(12)
- logdestr = "Log|dE|: " + ("%1.3f" % (dE)).rjust(6)
- iterstr = "SCF: " + ("%3i" % (iterations))
- avgfstr="Avg|F|: " + ("%2.3f" % (average)).rjust(6)
- maxfstr="Max|F|: " + ("%2.3f" % (maxforce)).rjust(6)
- timestr="Time: " + ("%3.2fm" % (cputime)).rjust(6)
- volstr="Vol.: " + ("%3.1f" % (volume)).rjust(5)
- except NameError:
- print ("Cannot understand this OUTCAR file...try to read ahead")
- continue
-
- if iterations == nelmax:
- sys.stdout.write(FAIL)
- #print " ^--- SCF cycle reached NELMAX. Check convergence!"
-
- if (dE < ediff):
- sys.stdout.write(OKGREEN)
-
- if spinpolarized:
- magstr="Mag: " + ("%2.2f" % (magmom)).rjust(6)
- print ("%s %s %s %s %s %s %s %s %s" % (stepstr,energystr,logdestr,iterstr,avgfstr,maxfstr,volstr,magstr,timestr))
- else:
- print ("%s %s %s %s %s %s %s %s" % (stepstr,energystr,logdestr,iterstr,avgfstr,maxfstr,volstr,timestr))
- # sys.stdout.write(("Step %3i Energy: %+3.6f Log|dE|: %+1.3f Avg|F|: %.6f Max|F|: %.6f SCF: %3i Time: %03.2fm") % (steps,energy,dE,average,maxforce,iterations,cputime))
-
- sys.stdout.write(ENDC)
-
- steps = steps + 1
- iterations = 0
- totaltime = totaltime + cputime
- cputime = 0.0
-
- i = i + 1
+ outcarfile = os.path.join(ini_path, d, '03_ISIF3snsb', 'OUTCAR')
+ print(j, outcarfile)
+ outcar = open(outcarfile, 'r')
+ outcarlines = outcar.readlines()
+
+ #Find max iterations
+ nelmax = int(
+ subprocess.Popen(
+ 'grep "NELM " ' + outcarfile, stdout=subprocess.PIPE,
+ shell=True).communicate()[0].split()[2][:-1])
+ nsw = int(subprocess.Popen('grep "NSW" ' + outcarfile + ' | head -1 | cut -d"=" -f2', stdout=subprocess.PIPE, shell=True).communicate()[0] )
+ natoms = get_number_of_atoms(outcarfile)
+ ediff = np.log10(float(get_ediff(outcarfile)))
+
+ re_energy = re.compile("free energy")
+ re_iteration = re.compile("Iteration")
+ re_timing = re.compile("LOOP:")
+ re_force = re.compile("TOTAL-FORCE")
+ re_mag = re.compile("number of electron")
+ re_volume = re.compile("volume of cell")
+
+ lastenergy = 0.0
+ energy = 0.0
+ steps = 1
+ iterations = 0
+ cputime = 0.0
+ totaltime = 0.0
+ dE = 0.0
+ magmom = 0.0
+ volume = 0.0
+ spinpolarized = False
+ #average = 0.0
+ #maxforce = 0.0
+
+ i = 0
+ for line in outcarlines:
+
+ if re_iteration.search(line):
+ iterations = iterations + 1
+
+ if re_force.search(line):
+ # Calculate forces here...
+ forces = []
+ magnitudes = []
+ for j in range(natoms):
+ parts = outcarlines[i+j+2].split()
+ x = float(parts[3])
+ y = float(parts[4])
+ z = float(parts[5])
+ forces.append([x,y,z])
+ magnitudes.append(np.sqrt(x**2 + y**2 + z**2))
+
+ average = sum(magnitudes)/natoms
+ maxforce = max(magnitudes)
+
+ if re_mag.search(line):
+ parts = line.split()
+ if len(parts) > 5 and parts[0].strip() != "NELECT":
+ spinpolarized = True
+ magmom = float(parts[5])
+
+ if re_timing.search(line):
+ cputime = cputime + float(line.split()[6])/60.0
+
+ if re_volume.search(line):
+ parts = line.split()
+ if len(parts) > 4:
+ volume = float(parts[4])
+
+ if re_energy.search(line):
+ lastenergy = energy
+ energy = float(line.split()[4])
+ dE = np.log10(abs(energy-lastenergy+1.0E-12))
+
+ # Construct output string
+ try:
+ stepstr = str(steps).rjust(4)
+ energystr = "Energy: " + ("%3.6f" % (energy)).rjust(12)
+ logdestr = "Log|dE|: " + ("%1.3f" % (dE)).rjust(6)
+ iterstr = "SCF: " + ("%3i" % (iterations))
+ avgfstr="Avg|F|: " + ("%2.3f" % (average)).rjust(6)
+ maxfstr="Max|F|: " + ("%2.3f" % (maxforce)).rjust(6)
+ timestr="Time: " + ("%3.2fm" % (cputime)).rjust(6)
+ volstr="Vol.: " + ("%3.1f" % (volume)).rjust(5)
+ except NameError:
+ print ("Cannot understand this OUTCAR file...try to read ahead")
+ continue
+
+ if iterations == nelmax:
+ sys.stdout.write(FAIL)
+ #print " ^--- SCF cycle reached NELMAX. Check convergence!"
+
+ if (dE < ediff):
+ sys.stdout.write(OKGREEN)
+
+ if spinpolarized:
+ magstr="Mag: " + ("%2.2f" % (magmom)).rjust(6)
+ print(
+ f"{stepstr} {energystr} {logdestr} {iterstr} {avgfstr} {maxfstr} {volstr} {magstr} {timestr}"
+ )
+ else:
+ print(
+ f"{stepstr} {energystr} {logdestr} {iterstr} {avgfstr} {maxfstr} {volstr} {timestr}"
+ )
+ # sys.stdout.write(("Step %3i Energy: %+3.6f Log|dE|: %+1.3f Avg|F|: %.6f Max|F|: %.6f SCF: %3i Time: %03.2fm") % (steps,energy,dE,average,maxforce,iterations,cputime))
+
+ sys.stdout.write(ENDC)
+
+ steps = steps + 1
+ iterations = 0
+ totaltime = totaltime + cputime
+ cputime = 0.0
+
+ i = i + 1
diff --git a/src/Main_file.py b/src/Main_file.py
index 6267c3c..59ca25d 100644
--- a/src/Main_file.py
+++ b/src/Main_file.py
@@ -28,22 +28,27 @@ def Introduction():
from strain import *
from create_phonon_directories import *
from evaluate_manually_Elastic_constants import *
-
+
from colorama import Fore, Back, Style, init
- from termcolor import colored
+ from termcolor import colored
from pylab import *
import multiprocessing as mp
import compileall
-
+
print ("This may take a while!")
compileall.compile_dir(".", force=1)
-
+
Introduction()
init(autoreset=True)
print(colored('@'*80,'red'), end = '\n', flush=True)
print("Number of processors Detected: ", mp.cpu_count())
- print(Back.MAGENTA + ' NB: POSCAR should be in VASP 5 format & without selective dynamics', end = '\n', flush=True)
- print(Style.RESET_ALL)
+ print(
+ f'{Back.MAGENTA} NB: POSCAR should be in VASP 5 format & without selective dynamics',
+ end='\n',
+ flush=True,
+ )
+
+ print(Style.RESET_ALL)
print(colored('-'*80,'red'), end = '\n', flush=True)
print('>>> USAGE: execute by typing python3 sys.argv[0]')
print(colored('~'*80,'red'), end = '\n', flush=True)
@@ -68,33 +73,33 @@ def Introduction():
while True:
option = input("Enter the option as listed above: ")
option = int(option)
-
+
if (option == 1):
poscar()
-
+
elif (option == 2):
VOL_P = main_poscar()
VOL_C = main_contcar()
volume_diff(VOL_P, VOL_C)
-
+
elif (option == 3):
energy_vs_volume()
-
+
elif (option == 4):
print("OUTCAR should be in the same directory from which this script is run ")
pool = mp.Pool(mp.cpu_count())
elastic_matrix_VASP_STRESS()
pool.close()
-
+
elif (option == 5):
fitting_energy_vs_volume_curve()
-
+
elif (option == 6):
poscar_VASP42VASP5()
-
+
elif (option == 7):
Elastic_strain()
-
+
elif (option == 8):
create_phonon_directories()
@@ -103,10 +108,10 @@ def Introduction():
elif (option == 10):
local_lattice_distortion_DEF1()
-
+
elif (option == 11):
local_lattice_distortion_DEF2()
-
+
else:
print ("INVALID OPTION")
diff --git a/src/contcar_poscar.py b/src/contcar_poscar.py
index 6566187..cee51e5 100644
--- a/src/contcar_poscar.py
+++ b/src/contcar_poscar.py
@@ -14,11 +14,14 @@ def main_poscar():
count = 0
os.system("rm out_POSCARS.dat")
- VOL_P = []; pos = []; kk = []; lattice = [];
+ VOL_P = []
+ pos = []
+ kk = []
+ lattice = [];
mypath = os.getcwd()
- print ("-"*100)
+ print ("-"*100)
print (Back.YELLOW + "{:15s} {:15s} {:15.6s} {:15.6s} {:15.6s} {:15.15s}".format("Directory", "# of atoms", "||a||", \
- "||b||", "||c||", "VOL_POS[A^3]"),end="\n" )
+ "||b||", "||c||", "VOL_POS[A^3]"),end="\n" )
print ("-"*100)
for entry in os.listdir(mypath):
@@ -26,101 +29,93 @@ def main_poscar():
#print (entry)
for file in os.listdir(entry):
if file == "POSCAR":
- count+=1; sum = 0
+ count+=1
+ sum = 0
filepath = os.path.join(entry, file)
- #f = open(filepath, 'r')
- #print (f.read())
- #f.close()
- fo = open(filepath, 'r')
- ofile=open('out_POSCARS.dat','a')
-
- #print (colored('>>>>>>>> Name of the file: ','red'), fo.name, end = '\n', flush=True)
-
- ofile.write (fo.name + '\n')
- ofile.write ("")
- firstline = fo.readline()
- secondfline = fo.readline()
- Latvec1 = fo.readline()
- #print ("Lattice vector 1:", (Latvec1), end = '')
- #ofile.write (Latvec1)
- Latvec2 = fo.readline()
- #print ("Lattice vector 2:", (Latvec2), end = '')
- #ofile.write (Latvec2)
- Latvec3 = fo.readline()
- #print ("Lattice vector 3:", (Latvec3), end = '')
- #ofile.write (Latvec3)
- elementtype=fo.readline()
- elementtype = elementtype.split()
- #print ("Types of elements:", str(elementtype), end = '')
- #ofile.write (str(elementtype))
- numberofatoms=fo.readline()
- #print ("Number of atoms:", (numberofatoms), end = '')
- #ofile.write ((numberofatoms))
- Coordtype=fo.readline()
+ with open(filepath, 'r') as fo:
+ ofile=open('out_POSCARS.dat','a')
+ #print (colored('>>>>>>>> Name of the file: ','red'), fo.name, end = '\n', flush=True)
+
+ ofile.write (fo.name + '\n')
+ ofile.write ("")
+ firstline = fo.readline()
+ secondfline = fo.readline()
+ Latvec1 = fo.readline()
+ #print ("Lattice vector 1:", (Latvec1), end = '')
+ #ofile.write (Latvec1)
+ Latvec2 = fo.readline()
+ #print ("Lattice vector 2:", (Latvec2), end = '')
+ #ofile.write (Latvec2)
+ Latvec3 = fo.readline()
+ #print ("Lattice vector 3:", (Latvec3), end = '')
+ #ofile.write (Latvec3)
+ elementtype=fo.readline()
+ elementtype = elementtype.split()
+ #print ("Types of elements:", str(elementtype), end = '')
+ #ofile.write (str(elementtype))
+ numberofatoms=fo.readline()
+ #print ("Number of atoms:", (numberofatoms), end = '')
+ #ofile.write ((numberofatoms))
+ Coordtype=fo.readline()
##########################---------------------------------------------------------
- #print ("**********-------------------# of Atoms--------------------")
-
- nat = numberofatoms.split()
- nat = [int(i) for i in nat]
- for i in nat:
- sum = sum + i
- numberofatoms = sum
- #print ("{} :: Number of atoms: {}".format(nat, numberofatoms) )
+ #print ("**********-------------------# of Atoms--------------------")
+
+ nat = numberofatoms.split()
+ nat = [int(i) for i in nat]
+ for i in nat:
+ sum = sum + i
+ numberofatoms = sum
+ #print ("{} :: Number of atoms: {}".format(nat, numberofatoms) )
##########################---------------------------------------------------------
- #print ("-----------------------Atomic positions-----------------")
- #print ("Coordtype:", (Coordtype), end = '')
- for x in range(int(numberofatoms)):
- coord = fo.readline().split()
- coord = [float(i) for i in coord]
- pos = pos + [coord]
- pos = np.array(pos)
- #print (pos)
-
- ofile.write ("\n")
- fo.close()
-##########################---------------------------------------------------------
+ #print ("-----------------------Atomic positions-----------------")
+ #print ("Coordtype:", (Coordtype), end = '')
+ for _ in range(int(numberofatoms)):
+ coord = fo.readline().split()
+ coord = [float(i) for i in coord]
+ pos = pos + [coord]
+ pos = np.array(pos)
+ #print (pos)
- a=[]; b=[]; c=[];
+ ofile.write ("\n")
Latvec1=Latvec1.split()
Latvec2=Latvec2.split()
Latvec3=Latvec3.split()
-
-##########################---------------------------------------------------------
- for ai in Latvec1: a.append(float(ai))
- for bi in Latvec2: b.append(float(bi))
- for ci in Latvec3: c.append(float(ci))
+
+ a = [float(ai) for ai in Latvec1]
+ b = [float(bi) for bi in Latvec2]
+ c = [float(ci) for ci in Latvec3]
#print ('a=', a)
- ofile.write ("'a=' {}\n".format(a))
+ ofile.write(f"'a=' {a}\n")
#print ('b=', b)
- ofile.write ("'b=' {}\n".format(b))
+ ofile.write(f"'b=' {b}\n")
#print ('c=', c)
- ofile.write ("'c=' {}\n".format(c))
+ ofile.write(f"'c=' {c}\n")
lld = local_lattice_distortion(a,b,c)
##########################---------------------------------------------------------
-
+
alpha, beta, gamma = lattice_angles(a,b,c)
- VOL_POS = np.dot(a, np.cross(b,c))
+ VOL_POS = np.dot(a, np.cross(b,c))
VOL_P.append(VOL_POS)
- ofile.write ("'\u03B1=' {} '\u03B2=' {} '\u03B3=' {}\n".format(alpha,beta,gamma))
- ofile.write ("'||a||=' {}\n".format(np.linalg.norm(a)))
- ofile.write ("'||b||=' {}\n".format(np.linalg.norm(b)))
- ofile.write ("'||c||=' {}\n".format(np.linalg.norm(c)))
+ ofile.write(f"'\u03B1=' {alpha} '\u03B2=' {beta} '\u03B3=' {gamma}\n")
+ ofile.write(f"'||a||=' {np.linalg.norm(a)}\n")
+ ofile.write(f"'||b||=' {np.linalg.norm(b)}\n")
+ ofile.write(f"'||c||=' {np.linalg.norm(c)}\n")
#print ("#####------------------------------------------------")
-
+
#print ("a={} \t ||a||={:10.6f}".format(a, np.linalg.norm(a)) )
#print ("b={} \t ||b||={:10.6f}".format(b, np.linalg.norm(b)) )
#print ("c={} \t ||c||={:10.6f}".format(c, np.linalg.norm(c)) )
print ("{:15s} {:6d} {:15.6f} {:15.6f} {:15.6f} {:15.6f}".format(fo.name, numberofatoms, np.linalg.norm(a), \
np.linalg.norm(b), np.linalg.norm(c), VOL_POS) )
-
+
print ("'\u03B1=' {:6.6f} '\u03B2=' {:6.6f} '\u03B3=' {:6.6f} g={:6.6f}".format(alpha,beta,gamma,lld))
print ("."*5)
#print ('Vol= {:6.6f} A^3'.format(VOL_POS))
ofile.write ("***************************************************\n")
ofile.close()
- print ("_"*30)
+ print ("_"*30)
print ("Number of folders detected: ", count)
return VOL_P
@@ -128,85 +123,85 @@ def main_poscar():
def main_contcar():
count = 0
os.system("rm out_CONTCARS.dat")
- VOL_C = []; pos = []; kk = []; lattice = []; sum = 0
- mypath = os.getcwd()
- print ("-"*100)
+ VOL_C = []
+ pos = []
+ kk = []
+ lattice = []
+ sum = 0
+ mypath = os.getcwd()
+ print ("-"*100)
print (Back.GREEN + "{:15s} {:15s} {:15.6s} {:15.6s} {:15.6s} {:15.15s}".format("Directory", "# of atoms", "||a||", \
- "||b||", "||c||", "VOL_CON[A^3]"),end="\n" )
- print ("-"*100)
- print(Style.RESET_ALL)
+ "||b||", "||c||", "VOL_CON[A^3]"),end="\n" )
+ print ("-"*100)
+ print(Style.RESET_ALL)
for entry in os.listdir(mypath):
if os.path.isdir(os.path.join(mypath, entry)):
for file in os.listdir(entry):
if file == "CONTCAR":
- count+=1; sum = 0
+ count+=1
+ sum = 0
filepath = os.path.join(entry, file)
- fo = open(filepath, 'r')
-
- ofile=open('out_CONTCARS.dat','a')
-
- #print (colored('>>>>>>>> Name of the file: ','yellow'), fo.name, end = '\n', flush=True)
- ofile.write (fo.name + '\n')
- ofile.write ("")
- firstline = fo.readline()
- secondfline = fo.readline()
- Latvec1 = fo.readline()
- Latvec2 = fo.readline()
- Latvec3 = fo.readline()
- elementtype=fo.readline()
- elementtype = elementtype.split()
- numberofatoms=fo.readline()
- Coordtype=fo.readline()
- #print ("Coordtype:", (Coordtype), end = '')
+ with open(filepath, 'r') as fo:
+ ofile=open('out_CONTCARS.dat','a')
+
+ #print (colored('>>>>>>>> Name of the file: ','yellow'), fo.name, end = '\n', flush=True)
+ ofile.write (fo.name + '\n')
+ ofile.write ("")
+ firstline = fo.readline()
+ secondfline = fo.readline()
+ Latvec1 = fo.readline()
+ Latvec2 = fo.readline()
+ Latvec3 = fo.readline()
+ elementtype=fo.readline()
+ elementtype = elementtype.split()
+ numberofatoms=fo.readline()
+ Coordtype=fo.readline()
+ #print ("Coordtype:", (Coordtype), end = '')
##########################---------------------------------------------------------
- #print ("**********-------------------# of Atoms--------------------")
-
- nat = numberofatoms.split()
- nat = [int(i) for i in nat]
- for i in nat:
- sum = sum + i
- numberofatoms = sum
- #print ("{} :: Number of atoms: {}".format(nat, numberofatoms) )
+ #print ("**********-------------------# of Atoms--------------------")
+
+ nat = numberofatoms.split()
+ nat = [int(i) for i in nat]
+ for i in nat:
+ sum = sum + i
+ numberofatoms = sum
+ #print ("{} :: Number of atoms: {}".format(nat, numberofatoms) )
##########################---------------------------------------------------------
- #print ("//////---------------Atomic positions-----------------")
- for x in range(int(numberofatoms)):
- coord = fo.readline().split()
- coord = [float(i) for i in coord]
- pos = pos + [coord]
- pos = np.array(pos)
- #print (pos)
-
- ofile.write ("\n")
- fo.close()
-##########################---------------------------------------------------------
- a=[]; b=[]; c=[];
+ #print ("//////---------------Atomic positions-----------------")
+ for _ in range(int(numberofatoms)):
+ coord = fo.readline().split()
+ coord = [float(i) for i in coord]
+ pos = pos + [coord]
+ pos = np.array(pos)
+ #print (pos)
+
+ ofile.write ("\n")
Latvec1=Latvec1.split()
Latvec2=Latvec2.split()
- Latvec3=Latvec3.split()
-##########################---------------------------------------------------------
- for ai in Latvec1: a.append(float(ai))
- for bi in Latvec2: b.append(float(bi))
- for ci in Latvec3: c.append(float(ci))
+ Latvec3=Latvec3.split()
+ a = [float(ai) for ai in Latvec1]
+ b = [float(bi) for bi in Latvec2]
+ c = [float(ci) for ci in Latvec3]
#print ('a=', a)
- ofile.write ("'a=' {}\n".format(a))
+ ofile.write(f"'a=' {a}\n")
#print ('b=', b)
- ofile.write ("'b=' {}\n".format(b))
+ ofile.write(f"'b=' {b}\n")
#print ('c=', c)
- ofile.write ("'c=' {}\n".format(c))
- lld = local_lattice_distortion(a,b,c)
+ ofile.write(f"'c=' {c}\n")
+ lld = local_lattice_distortion(a,b,c)
##########################---------------------------------------------------------
alpha, beta, gamma = lattice_angles(a,b,c)
VOL_CON = np.dot(a, np.cross(b,c))
VOL_C.append(VOL_CON)
-
- ofile.write ("'\u03B1=' {} '\u03B2=' {} '\u03B3=' {}\n".format(alpha,beta,gamma))
- ofile.write ("'||a||=' {}\n".format(np.linalg.norm(a)))
- ofile.write ("'||b||=' {}\n".format(np.linalg.norm(b)))
- ofile.write ("'||c||=' {}\n".format(np.linalg.norm(c)))
+
+ ofile.write(f"'\u03B1=' {alpha} '\u03B2=' {beta} '\u03B3=' {gamma}\n")
+ ofile.write(f"'||a||=' {np.linalg.norm(a)}\n")
+ ofile.write(f"'||b||=' {np.linalg.norm(b)}\n")
+ ofile.write(f"'||c||=' {np.linalg.norm(c)}\n")
#print ("-"*80)
-
+
#print ("a={} \t ||a||={:10.6f}".format(a, np.linalg.norm(a)) )
#print ("b={} \t ||b||={:10.6f}".format(b, np.linalg.norm(b)) )
#print ("c={} \t ||c||={:10.6f}".format(c, np.linalg.norm(c)) )
@@ -218,8 +213,8 @@ def main_contcar():
#print ('Vol= {:6.6f} A^3'.format(VOL_CON), end="\n")
ofile.write ("***************************************************\n")
ofile.close()
- #print (VOL_C)
- print ("-"*80)
+ #print (VOL_C)
+ print ("-"*80)
print ("Number of folders detected: ", count)
return VOL_C
diff --git a/src/create_phonon_directories.py b/src/create_phonon_directories.py
index 97c6e23..f7378c6 100644
--- a/src/create_phonon_directories.py
+++ b/src/create_phonon_directories.py
@@ -26,8 +26,8 @@
def copy_fntn(k):
for i in range(1, int(k)+1):
- shutil.rmtree('disp-'+str(i).zfill(3), ignore_errors=True) #overwrite a directory
- os.mkdir('disp-'+str(i).zfill(3))
+ shutil.rmtree(f'disp-{str(i).zfill(3)}', ignore_errors=True)
+ os.mkdir(f'disp-{str(i).zfill(3)}')
#subprocess.check_call(['mkdir', 'disp-'+str(i).zfill(3)])
#subprocess.check_call('cd disp-'+str(i).zfill(3), shell=True)
# print(str(i).zfill(3))
@@ -35,14 +35,18 @@ def copy_fntn(k):
# The zfill()/rjust() is a function associated with the string object.
# 3 is the expected length of the string after padding.
#
- os.system('cp INCAR disp-'+str(i).zfill(3) )
- os.system('cp POTCAR disp-'+str(i).zfill(3) )
- os.system('cp KPOINTS disp-'+str(i).zfill(3) )
+ os.system(f'cp INCAR disp-{str(i).zfill(3)}')
+ os.system(f'cp POTCAR disp-{str(i).zfill(3)}')
+ os.system(f'cp KPOINTS disp-{str(i).zfill(3)}')
#os.system('cp job.sh disp-'+str(i).zfill(3) )
#os.system('cp WAVECAR disp-'+str(i).zfill(3) )
- subprocess.call(['cp','-r','POSCAR-'+str(i).zfill(3),'disp-'+str(i).zfill(3)], shell = False)
- os.chdir('disp-'+str(i).zfill(3))
- shutil.copyfile("POSCAR-"+str(i).zfill(3), "POSCAR")
+ subprocess.call(
+ ['cp', '-r', f'POSCAR-{str(i).zfill(3)}', f'disp-{str(i).zfill(3)}'],
+ shell=False,
+ )
+
+ os.chdir(f'disp-{str(i).zfill(3)}')
+ shutil.copyfile(f"POSCAR-{str(i).zfill(3)}", "POSCAR")
#os.system('ls')
os.chdir('../')
@@ -55,13 +59,15 @@ def create_phonon_directories():
print (' ERROR: POSCAR does not exist here.')
sys.exit(0)
else:
- counter=0
- mypath = os.getcwd()
- for entry in os.listdir(mypath):
- if os.path.isfile(os.path.join(mypath, entry)):
- if fnmatch.fnmatchcase(entry,'POSCAR-*'):
- counter+=1
- print ("# of POSCAR-* files generated with Phonopy code: ----> {}".format(counter) )
+ mypath = os.getcwd()
+ counter = sum(
+ 1
+ for entry in os.listdir(mypath)
+ if os.path.isfile(os.path.join(mypath, entry))
+ and fnmatch.fnmatchcase(entry, 'POSCAR-*')
+ )
+
+ print(f"# of POSCAR-* files generated with Phonopy code: ----> {counter}")
print("*"*80)
start_time = time.time()
diff --git a/src/elastic_constants.py b/src/elastic_constants.py
index 5b95e9c..0e2e79a 100644
--- a/src/elastic_constants.py
+++ b/src/elastic_constants.py
@@ -11,16 +11,18 @@
class Elastic_Matrix:
def print_Cij_Matrix(): ###EXERCISE
- Bij = []; C = "C"
- for i in range(0, 6, 1):
+ Bij = []
+ C = "C"
+ for i in range(6):
Bij.append([])
- for j in range(1,7, 1):
+ for j in range(1, 7):
Bij[i].append((C + str(i+1) + str(j)))
- l = np.matrix(Bij)
- return l
+ return np.matrix(Bij)
def langragian_strain():
- kk = ('x', 'y', 'z'); C = "E"; Eij=[]
+ kk = ('x', 'y', 'z')
+ C = "E"
+ Eij=[]
eta = {
"xx" : 11,
"yy" : 22,
@@ -30,42 +32,37 @@ def langragian_strain():
"xy" : 12 }
print ( "The Langragian strain in Voigt notation: ")
print ( eta.items())
-
+
for n in kk:
- for m in kk:
- Eij.append( (C + str(n) + str(m)) )
- ls = np.matrix(Eij).reshape(3,3)
- return ls
+ Eij.extend(C + str(n) + str(m) for m in kk)
+ return np.matrix(Eij).reshape(3,3)
def elastic_matrix_VASP_STRESS():
- while True:
+ while True:
if not os.path.exists('OUTCAR'):
print (' ERROR: OUTCAR does not exist here.')
sys.exit(0)
-
+
s=np.zeros((6,6))
- c=np.zeros((6,6))
- file = open("OUTCAR",'r')
- lines = file.readlines()
- file.close()
-
+ c=np.zeros((6,6))
+ with open("OUTCAR",'r') as file:
+ lines = file.readlines()
for i in lines:
if "TOTAL ELASTIC MODULI (kBar)" in i:
ll=lines.index(i)
if "LATTYP" in i:
crystaltype=str(i.split()[3])
- print ("DETECTED CRYSTAL FROM OUTCAR:", crystaltype)
- print (" ")
- for i in range(0,6):
+ print ("DETECTED CRYSTAL FROM OUTCAR:", crystaltype)
+ print (" ")
+ for i in range(6):
l=lines[ll+3+i] # indexing a line in huge file
word = l.split()
s[i][:] = word[1:7]
- for j in range(0,6):
+ for j in range(6):
c[i][j] = float(s[i][j])/10.0
-
##########-------------------stress tensor------------------
-
+
Cij = np.matrix(c)
np.set_printoptions(precision=4, suppress=True)
print (Cij)
@@ -76,114 +73,108 @@ def elastic_matrix_VASP_STRESS():
if evals.all() > 0:
print(evals)
print("All Eigen values are positive")
-
##########-------------------Compliance tensor------------------
##########------------------- s_{ij} = C_{ij}^{-1}
Sij = np.linalg.inv(Cij)
-
#---------------------------- ELASTIC PROPERTIES -----------------------------------
stability_test(Cij, crystaltype)
-
######## -----------------------------VOIGT-------------------------------
-
+
'''Voigt bulk modulus (GPa)'''
-
+
#9K_v = (C_{11}+C_{22}+C_{33}) + 2(C_{12} + C_{23} + C_{31})
Bv = ((Cij[0,0] + Cij[1,1] + Cij[2,2]) + 2 * (Cij[0,1] + Cij[1,2] + Cij[2,0])) / 9.0
-
+
'''Voigt shear modulus (GPa)'''
-
+
#15*G_v = (C_{11}+C_{22}+C_{33}) - (C_{12} + C_{23} + C_{31}) + 3(C_{44} + C_{55} + C_{66})$
Gv = ((Cij[0,0] + Cij[1,1] + Cij[2,2]) - (Cij[0,1] + Cij[1,2] + Cij[2,0])
+ 3 * (Cij[3,3] + Cij[4,4] + Cij[5,5]))/15.0
-
+
## Young's: Voigt
Ev = (9*Bv*Gv)/(3*Bv + Gv)
-
+
## Poisson's ratio: Voigt
NuV = (3*Bv - Ev)/(6*Bv)
-
+
######## -----------------------------REUSS-------------------------------
-
+
# Reuss bulk modulus K_R $(GPa)$
# 1/K_R = (s_{11}+s_{22}+s_{33}) + 2(s_{12} + s_{23} + s_{31})$
Br = 1/((Sij[0,0] + Sij[1,1] + Sij[2,2]) + 2*(Sij[0,1] + Sij[1,2] + Sij[2,0]))
-
+
# Reuss shear modulus G_v $(GPa)$
# 15/G_R = 4*(s_{11}+s_{22}+s_{33}) - 4*(s_{12} + s_{23} + s_{31}) + 3(s_{44} + s_{55} + s_{66})$
Gr = (4 * (Sij[0,0] + Sij[1,1] + Sij[2,2]) - 4*(Sij[0,1] + Sij[1,2] + Sij[2,0])
+ 3 * (Sij[3,3] + Sij[4,4] + Sij[5,5]))
Gr = 15.0/Gr
-
+
## Young's: Reuss
Er = (9*Br*Gr)/(3*Br + Gr)
-
+
## Poisson's ratio: Reuss
NuR = (3*Br - Er)/(6*Br)
##########################################################################
######## -----------------------------Averages-------------------------------
-
+
# #Hill bulk modulus K_{VRH}$ $(GPa)$
# K_{VRH} = (K_R + K_v)/2
B_H = (Bv + Br)/2
#print ("VRH bulk modulus (GPa): %20.8f " %(B_H) )
-
+
# Hill shear modulus G_{VRH}$ $(GPa)$
# G_{VRH} = (G_R + G_v)/2
G_H = (Gv + Gr)/2
#print ("VRH shear modulus (GPa): %20.8f " %(G_H) )
-
+
# Young modulus E = 9BG/(3B+G)
#E_H = (9 * B_H * G_H) / (3 * B_H + G_H)
E_H = (Ev + Er)/2
#print ("Young modulus E : {:1.8s} {:20.8f}".format(" ",E_H) )
-
+
# ### Isotropic Poisson ratio $\mu
# $\mu = (3K_{VRH} - 2G_{VRH})/(6K_{VRH} + 2G_{VRH})$
#nu_H = (3 * B_H - 2 * G_H) / (6 * B_H + 2 * G_H )
nu_H = (NuV + NuR) / 2
#print ("Isotropic Poisson ratio: {:15.8f} ".format(nu_H) )
-
+
## Elastic Anisotropy
## Zener anisotropy for cubic crystals only
A = 2*(c44)/(c11-c12)
-
+
# Universal Elastic Anisotropy AU
AU = (Bv/Br) + 5*(Gv/Gr) - 6.0
-
+
# C' tetragonal shear modulus
C = (c11-c12)/2
-
+
ratio_V = Bv/Gv
ratio_R = Br/Gr
- print ("{:_^80}".format("GPa"))
- print ("{:25.8s} {:15.8s} {:15.8s} {:15.8s}".format(" ","Voigt", "Reuss ", "Hill") )
+ print ("{:_^80}".format("GPa"))
+ print ("{:25.8s} {:15.8s} {:15.8s} {:15.8s}".format(" ","Voigt", "Reuss ", "Hill") )
print ("{:16.20s} {:15.3f} {:15.3f} {:15.3f}".format("Bulk Modulus",Bv, Br, B_H) )
print ("{:16.20s} {:15.3f} {:15.3f} {:15.3f}".format("Shear Modulus",Gv, Gr, G_H) )
print ("{:16.20s} {:15.3f} {:15.3f} {:15.3f}".format("Young Modulus",Ev, Er, E_H) )
- print ("{:-^80}".format("-"))
+ print ("{:-^80}".format("-"))
print ("{:16.20s} {:15.3f} {:15.3f} {:15.3f}".format("Poisson ratio ", NuV, NuR, nu_H) )
print ("{:16.20s} {:15.3f} {:15.3f} {:15.3f}({:5.3f})".format("B/G ratio ",Bv/Gv,Br/Gr, B_H/G_H, G_H/B_H) )
- print ("{:16.20s} {:15.3s} {:15.3s} {:15.3f}".format("Zener ratio Az",'','', A) )
+ print ("{:16.20s} {:15.3s} {:15.3s} {:15.3f}".format("Zener ratio Az",'','', A) )
print ("{:16.20s} {:15.3s} {:15.3s} {:15.3f}".format("Avr ratio ",'','', (Gv-Gr)/(Gv+Gr)) )
print ("{:16.20s} {:15.3s} {:15.3s} {:15.3f}".format("AU ",'','', AU) )
print ("{:16.20s} {:15.3s} {:15.3s} {:15.3f}".format("Cauchy pressure ",'','', (c12-c44)) )
print ("{:16.20s} {:15.3s} {:15.3s} {:15.3f}".format("C'tetra Shear ",'','', C) )
-
+
print ("-"*80)
return Sij
def ductile_test(ratio):
- if(ratio > 1.75):
- return "ductile"
- else:
- return "brittle"
+ return "ductile" if (ratio > 1.75) else "brittle"
###
def stability_test(matrix, crystaltype):
diff --git a/src/energy.py b/src/energy.py
index dccacdc..517e014 100644
--- a/src/energy.py
+++ b/src/energy.py
@@ -16,67 +16,66 @@ def energy_vs_volume():
os.system("rm energy-vs-volume energy-vs-strain")
eV2Hartree=0.036749309
Ang32Bohr3=6.74833304162
-
- E=[]; dir_list=[]; count = 0; dir_E=[];
- vol_cell=[]; strain_file=[]; strain_value=[] # strain_value is deformation
-
+
+ E=[]
+ dir_list=[]
+ count = 0
+ dir_E=[];
+ vol_cell=[]
+ strain_file=[]
+ strain_value=[] # strain_value is deformation
+
print (" >>>>> Extracting Energy from directories <<<<<<")
for entry in os.listdir(mypath):
if not os.path.exists('strain-01'):
print (' ERROR: strain-* does not exist here.')
- sys.exit(0)
- if fnmatch.fnmatchcase(entry,'strain-*'):
- f = open(entry,'r')
- lines = f.readline() #Read first line only
- strain_value.append( float(lines) )
- f.close()
+ sys.exit(0)
+ if fnmatch.fnmatchcase(entry,'strain-*'):
+ with open(entry,'r') as f:
+ lines = f.readline() #Read first line only
+ strain_value.append( float(lines) )
if os.path.isfile(os.path.join(mypath, entry)):
strain_file.append(entry)
if os.path.isdir(os.path.join(mypath, entry)):
dir_list.append(entry)
-
+
for file in os.listdir(entry):
if file == "OUTCAR":
filepath = os.path.join(entry, file)
if not os.path.exists(filepath):
print (' ERROR: OUTCAR does not exist here.')
sys.exit(0)
- f = open(filepath,'r')
- lines = f.readlines()
- f.close()
-
+ with open(filepath,'r') as f:
+ lines = f.readlines()
for i in lines:
if " free energy TOTEN =" in i:
m=float(i.split()[4])
if " volume of cell :" in i:
v=float(i.split()[4])
- vol_cell.append(v)
+ vol_cell.append(v)
E.append(m)
- count+=1
+ count+=1
print("# of folders detected: ", count)
print ("Directory :%10.6s %14s %18s %25.20s " % ("Folder", "Energy(eV)", "Vol_of_cell(A^3)", "strain_deformation" ))
-
+
for i in range(math.floor(count/2)): # 0 to 4
print ("Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f" %(dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
if (bool(math.floor(count/2))):
i = math.floor(count/2)
print(Back.GREEN + 'Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f <--Ref' % (dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
print(Style.RESET_ALL, end="")
- for i in range(math.ceil(count/2), count, 1):
- print ("Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f" %(dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
+ for i in range(math.ceil(count/2), count):
+ print ("Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f" %(dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
#rc = subprocess.Popen(['bash', 'extract_energy.sh'])
-
+
print (colored('ENERGIES & VOLUMES ARE WRITTEN IN ATOMIC UNITS TO A FILE <energy-vs-volume>','yellow'), end = '\n', flush=True)
print (colored('IT WILL BE READ BY ELASTIC SCRIPTS FOR POSTPROCESSING eV-->Ha; A^3-->Bohr^3','yellow'), end = '\n', flush=True)
-
- file = open("energy-vs-volume",'w')
- for i in range(count):
- file.write ("%14.6f %14.6f\n" %(vol_cell[i] * Ang32Bohr3, E[i] * eV2Hartree))
- file.close()
- file = open("energy-vs-strain",'w')
- for i in range(count):
- file.write ("%12.6f %14.6f\n" %(strain_value[i], E[i] * eV2Hartree))
- file.close()
+ with open("energy-vs-volume",'w') as file:
+ for i in range(count):
+ file.write ("%14.6f %14.6f\n" %(vol_cell[i] * Ang32Bohr3, E[i] * eV2Hartree))
+ with open("energy-vs-strain",'w') as file:
+ for i in range(count):
+ file.write ("%12.6f %14.6f\n" %(strain_value[i], E[i] * eV2Hartree))
###
\ No newline at end of file
diff --git a/src/evaluate_manually_Elastic_constants.py b/src/evaluate_manually_Elastic_constants.py
index 6248fec..a5335a3 100644
--- a/src/evaluate_manually_Elastic_constants.py
+++ b/src/evaluate_manually_Elastic_constants.py
@@ -28,11 +28,11 @@ def mechanical_properties():
# paper. For three distortions of the crystal three constants are extracted:
# C11, C44, C12. For cubic system the matrix is symmetric.
- print (CRED + "Values has to be supplied in the script::" + CEND)
+ print(f"{CRED}Values has to be supplied in the script::{CEND}")
########################## INPUT PARAMETERS assuming cubic #######################
-
- c11=c22=c33=123.22 ;
+
+ c11=c22=c33=123.22 ;
c44=c55=c66=47
B = 120
#c12=c21=c13=c31=c23=c32=(1/6) * (9 * B - 3 * c11)
@@ -50,107 +50,109 @@ def mechanical_properties():
[0 ,0 ,0 ,0 ,0 ,c66] ]
Cij=np.matrix(Cij).reshape((6,6))
print (Cij)
-
+
########------------Calculate variance and Mean of the Cij -----------
# ONLY for C11 , C22 , C33
Cij_stat = [ Cij[0,0] , Cij[1,1] , Cij[2,2] ]
print (">>>", (Cij_stat) )
- STD = np.std(Cij_stat); Var = np.var(Cij_stat) ; Mean = np.mean(Cij_stat)
+ STD = np.std(Cij_stat)
+ Var = np.var(Cij_stat)
+ Mean = np.mean(Cij_stat)
print ("{:6.3s} {:6.3s} {:6.3s}".format("C", "STD", "Var") )
print ("{:6.3f} {:6.3f} {:6.3f}".format(Mean, STD/np.sqrt(3), Var) )
######## ------------------------------------------------------------
-
+
evals, eigenvec = LA.eig(Cij)
print ("-"*40)
print("Eigenvalues are: ", evals>0)
print ("-"*40)
# ### Compliance tensor s_{ij}$ $(GPa^{-1})$
# s_{ij} = C_{ij}^{-1}$
-
+
print (CRED +"{:_^80s}".format("The COMPLIANCE MATRIX Sij is:")+ CEND)
Sij = np.linalg.inv(Cij)
- print ("{} ".format(Sij))
+ print(f"{Sij} ")
print ("-"*80)
-
+
######## -----------------------------VOIGT-------------------------------
-
+
'''Voigt bulk modulus (GPa)'''
-
+
#9K_v = (C_{11}+C_{22}+C_{33}) + 2(C_{12} + C_{23} + C_{31})
Bv = ((Cij[0,0] + Cij[1,1] + Cij[2,2]) + 2 * (Cij[0,1] + Cij[1,2] + Cij[2,0])) / 9.0
-
+
'''Voigt shear modulus (GPa)'''
-
+
#15*G_v = (C_{11}+C_{22}+C_{33}) - (C_{12} + C_{23} + C_{31}) + 3(C_{44} + C_{55} + C_{66})$
Gv = ((Cij[0,0] + Cij[1,1] + Cij[2,2]) - (Cij[0,1] + Cij[1,2] + Cij[2,0])
+ 3 * (Cij[3,3] + Cij[4,4] + Cij[5,5]))/15.0
-
+
## Young's: Voigt
Ev = (9*Bv*Gv)/(3*Bv + Gv)
-
+
## Poisson's ratio: Voigt
NuV = (3*Bv - Ev)/(6*Bv)
-
+
######## -----------------------------REUSS-------------------------------
-
+
# Reuss bulk modulus K_R $(GPa)$
# 1/K_R = (s_{11}+s_{22}+s_{33}) + 2(s_{12} + s_{23} + s_{31})$
Br = 1/((Sij[0,0] + Sij[1,1] + Sij[2,2]) + 2*(Sij[0,1] + Sij[1,2] + Sij[2,0]))
-
+
# Reuss shear modulus G_v $(GPa)$
# 15/G_R = 4*(s_{11}+s_{22}+s_{33}) - 4*(s_{12} + s_{23} + s_{31}) + 3(s_{44} + s_{55} + s_{66})$
Gr = (4 * (Sij[0,0] + Sij[1,1] + Sij[2,2]) - 4*(Sij[0,1] + Sij[1,2] + Sij[2,0])
+ 3 * (Sij[3,3] + Sij[4,4] + Sij[5,5]))
Gr = 15.0/Gr
-
+
## Young's: Reuss
Er = (9*Br*Gr)/(3*Br + Gr)
-
+
## Poisson's ratio: Reuss
NuR = (3*Br - Er)/(6*Br)
##########################################################################
######## -----------------------------Averages-------------------------------
-
+
# #Hill bulk modulus K_{VRH}$ $(GPa)$
# K_{VRH} = (K_R + K_v)/2
B_H = (Bv + Br)/2
#print ("VRH bulk modulus (GPa): %20.8f " %(B_H) )
-
+
# Hill shear modulus G_{VRH}$ $(GPa)$
# G_{VRH} = (G_R + G_v)/2
G_H = (Gv + Gr)/2
#print ("VRH shear modulus (GPa): %20.8f " %(G_H) )
-
+
# Young modulus E = 9BG/(3B+G)
#E_H = (9 * B_H * G_H) / (3 * B_H + G_H)
E_H = (Ev + Er)/2
#print ("Young modulus E : {:1.8s} {:20.8f}".format(" ",E_H) )
-
+
# ### Isotropic Poisson ratio $\mu
# $\mu = (3K_{VRH} - 2G_{VRH})/(6K_{VRH} + 2G_{VRH})$
#nu_H = (3 * B_H - 2 * G_H) / (6 * B_H + 2 * G_H )
nu_H = (NuV + NuR) / 2
#print ("Isotropic Poisson ratio: {:15.8f} ".format(nu_H) )
-
+
## Elastic Anisotropy
## Zener anisotropy for cubic crystals only
A = 2*(c44)/(c11-c12)
-
+
# Universal Elastic Anisotropy AU
AU = (Bv/Br) + 5*(Gv/Gr) - 6.0
-
+
# C' tetragonal shear modulus
C = (c11-c12)/2
-
+
ratio_V = Bv/Gv
ratio_R = Br/Gr
-
+
print ("{:30.8s} {:20.8s} {:20.8s} {:20.8s}".format(" ","Voigt", "Reuss ", "Hill") )
- print ("{:_^80}".format("GPa"))
+ print ("{:_^80}".format("GPa"))
print ("{:16.20s} {:20.3f} {:20.3f} {:20.3f}".format("Bulk Modulus",Bv, Br, B_H) )
print ("{:16.20s} {:20.3f} {:20.3f} {:20.3f}".format("Shear Modulus",Gv, Gr, G_H) )
print ("{:16.20s} {:20.3f} {:20.3f} {:20.3f}".format("Young Modulus",Ev, Er, E_H) )
@@ -161,7 +163,7 @@ def mechanical_properties():
print ("{:16.20s} {:20.3s} {:20.3s} {:20.3f}".format("AU ",'','', AU) )
print ("{:16.20s} {:20.3s} {:20.3s} {:20.3f}".format("Cauchy pressure ",'','', (c12-c44)) )
print ("{:16.20s} {:20.3s} {:20.3s} {:20.3f}".format("C'tetra Shear ",'','', C) )
-
+
print ("-"*80)
return Sij
@@ -169,18 +171,15 @@ def elastic_anisotropy(S):
#Elastic anisotropy of crystals is important since it correlates with the
#possibility to induce micro-cracks in materials
v=[]
- f = open('CONTCAR','r')
- lines = f.readlines()
- f.close()
-
+ with open('CONTCAR','r') as f:
+ lines = f.readlines()
s = float(lines[1])
- for i in range(2,5,1):
+ for i in range(2, 5):
l = s*float(lines[i].split()[0]), s*float(lines[i].split()[1]), s*float(lines[i].split()[2])
v.append( l )
v = np.mat(v)
print("Reading lattice vectors into matrix form")
print(v[:])
-
# l1, l2, l3 are the direction cosines
n1 = np.linalg.norm(v[0])
n2 = np.linalg.norm(v[1])
@@ -188,7 +187,7 @@ def elastic_anisotropy(S):
l1 = math.degrees(math.acos(np.dot(v[0],np.transpose(v[1]))/(n1*n2)))
l2 = math.degrees(math.acos(np.dot(v[1],np.transpose(v[2]))/(n2*n3)))
l3 = math.degrees(math.acos(np.dot(v[2],np.transpose(v[0]))/(n3*n1)))
- print ("l's are direction Cosines:: l1={:6.6f} l2={:6.6f} l3={:6.6f}".format(l1, l2, l3) )
+ print ("l's are direction Cosines:: l1={:6.6f} l2={:6.6f} l3={:6.6f}".format(l1, l2, l3) )
B = 1/(3*(S[0,0] + 2*S[0,1] ))
G = 1/( S[3,3] +4*(S[0,0] - S[0,1] -1/(2*S[3,3]) ) * (l1**2 * l2**2 + l2**2 * l3**2 + l1**2 * l3**2) )
E = 1/( S[0,0] -2*(S[0,0] - S[0,1] -1/(2*S[3,3]) ) * (l1**2 * l2**2 + l2**2 * l3**2 + l1**2 * l3**2) )
@@ -200,7 +199,9 @@ def local_lattice_distortion_DEF1():
#print (Back.YELLOW + "Wang, S. Atomic structure modeling of multi-principal-element alloys by the principle")
#print (Back.YELLOW + "of maximum entropy. Entropy 15, 5536–5548 (2013).")
print (">"*10,"HUME ROTHERY RULE")
- C_i=C=0.2 ; r_avg = 0.0; del_sum=0.0
+ C_i=C=0.2
+ r_avg = 0.0
+ del_sum=0.0
elements = ["Nb", "Hf", "Ta", "Ti", "Zr"]
eta = {
"Nb" : 1.98,
@@ -208,17 +209,17 @@ def local_lattice_distortion_DEF1():
"Ta" : 2.00,
"Ti" : 1.76,
"Zr" : 2.06, }
-
+
print (" {element: atomic radius}")
print (eta)
-
+
for i in elements:
r_avg = r_avg + C * eta[i]
-
+
for j in elements:
del_sum = del_sum + C * ( 1 - float(eta[j]) / r_avg )**2
- del_sum = 100 * np.sqrt(del_sum)
- print("HEA_atomic_size_mismatch: \u03B4={}".format(del_sum))
+ del_sum = 100 * np.sqrt(del_sum)
+ print(f"HEA_atomic_size_mismatch: \u03B4={del_sum}")
############
def local_lattice_distortion_DEF2():
@@ -227,42 +228,43 @@ def local_lattice_distortion_DEF2():
print (" Phys. Rev. Mater. 1, 23404 (2017).")
print (" (***) Different definition of the atomic radius for the description ")
print (" of the local lattice distortion in HEAs")
-
+
if not os.path.exists('POSCAR' and 'CONTCAR'):
print ('>>> ERROR: POSCAR & CONTCAR does not exist (Both should be in the same directory)')
sys.exit(0)
print('Reading POSCAR and CONTCAR ... \n')
-
- x = []; y =[]; z=[]
- xp =[]; yp = []; zp = []; temp=0
-
- f = open('POSCAR','r')
- lines_poscar = f.readlines()
- f.close()
-
- f = open('CONTCAR','r')
- lines_contcar = f.readlines()
- f.close()
-
+
+ x = []
+ y =[]
+ z=[]
+ xp =[]
+ yp = []
+ zp = []
+ temp=0
+
+ with open('POSCAR','r') as f:
+ lines_poscar = f.readlines()
+ with open('CONTCAR','r') as f:
+ lines_contcar = f.readlines()
file_P = ase.io.read('POSCAR')
pos = file_P.get_cell_lengths_and_angles()
- print (CRED + "POSCAR=>Length&Angles->{}".format(pos) + CEND)
+ print(CRED + f"POSCAR=>Length&Angles->{pos}" + CEND)
file_C = ase.io.read('CONTCAR')
- con = file_C.get_cell_lengths_and_angles()
- print (CRED + "CONTCAR=>Length&Angles->{}".format(con) + CEND)
+ con = file_C.get_cell_lengths_and_angles()
+ print(CRED + f"CONTCAR=>Length&Angles->{con}" + CEND)
print ("Cell vectors difference:: ",con-pos)
-
+
sum_atoms = lines_poscar[6].split() ### reading 7th lines for reading # of atoms
sum_atoms = [int(i) for i in sum_atoms]
sum_atoms = sum(sum_atoms)
-
+
for i in lines_poscar:
if "Direct" in i:
lp=lines_poscar.index(i)
for j in lines_contcar:
if "Direct" in j:
lc=lines_contcar.index(j)
-
+
for i in range(sum_atoms):
x, y, z = lines_poscar[lp+1+i].split()
xc, yc, zc = lines_contcar[lp+1+i].split()
@@ -270,7 +272,7 @@ def local_lattice_distortion_DEF2():
xc = float(xc); yc = float(yc); zc = float(zc)
temp = temp + np.sqrt( (x-xc)**2 + (y-yc)**2 + (z-zc)**2 )
temp = temp/sum_atoms
- print("local lattice distortion (LLD): \u0394d={}".format(temp))
+ print(f"local lattice distortion (LLD): \u0394d={temp}")
if __name__ == "__main__":
diff --git a/src/fit_energy_vs_vol.py b/src/fit_energy_vs_vol.py
index 0d065bb..e853950 100644
--- a/src/fit_energy_vs_vol.py
+++ b/src/fit_energy_vs_vol.py
@@ -16,69 +16,68 @@ def energy_vs_volume():
os.system("rm energy-vs-volume energy-vs-strain")
eV2Hartree=0.036749309
Ang32Bohr3=6.74833304162
-
- E=[]; dir_list=[]; count = 0; dir_E=[];
- vol_cell=[]; strain_file=[]; strain_value=[] # strain_value is deformation
-
+
+ E=[]
+ dir_list=[]
+ count = 0
+ dir_E=[];
+ vol_cell=[]
+ strain_file=[]
+ strain_value=[] # strain_value is deformation
+
print (" >>>>> Extracting Energies from directories <<<<<<")
for entry in os.listdir(mypath):
if not os.path.exists('strain-01'):
print (' ERROR: strain-* files does not exist. create a strain file for each deformation values.')
- sys.exit(0)
- if fnmatch.fnmatchcase(entry,'strain-*'):
- f = open(entry,'r')
- lines = f.readline() #Read first line only
- strain_value.append( float(lines) )
- f.close()
+ sys.exit(0)
+ if fnmatch.fnmatchcase(entry,'strain-*'):
+ with open(entry,'r') as f:
+ lines = f.readline() #Read first line only
+ strain_value.append( float(lines) )
if os.path.isfile(os.path.join(mypath, entry)):
strain_file.append(entry)
if os.path.isdir(os.path.join(mypath, entry)):
dir_list.append(entry)
-
+
for file in os.listdir(entry):
if file == "OUTCAR":
filepath = os.path.join(entry, file)
if not os.path.exists(filepath):
print (' ERROR: OUTCAR does not exist here.')
sys.exit(0)
- f = open(filepath,'r')
- lines = f.readlines()
- f.close()
-
+ with open(filepath,'r') as f:
+ lines = f.readlines()
for i in lines:
if " free energy TOTEN =" in i:
m=float(i.split()[4])
if " volume of cell :" in i:
v=float(i.split()[4])
- vol_cell.append(v)
+ vol_cell.append(v)
E.append(m)
- count+=1
+ count+=1
print("# of folders detected: ", count)
print ("Directory :%10.6s %14s %18s %25.20s " % ("Folder", "Energy(eV)", "Vol_of_cell(A^3)", "strain_deformation" ))
-
+
for i in range(math.floor(count/2)): # 0 to 4
print ("Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f" %(dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
if (bool(math.floor(count/2))):
i = math.floor(count/2)
print(Back.GREEN + 'Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f <--Ref' % (dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
print(Style.RESET_ALL, end="")
- for i in range(math.ceil(count/2), count, 1):
- print ("Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f" %(dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
+ for i in range(math.ceil(count/2), count):
+ print ("Folder name: %10.10s %16.8f %16.8f %16.12s %14.4f" %(dir_list[i], E[i], vol_cell[i], strain_file[i], strain_value[i] ))
#rc = subprocess.Popen(['bash', 'extract_energy.sh'])
-
+
print (colored('ENERGIES & VOLUMES ARE WRITTEN IN ATOMIC UNITS TO A FILE <energy-vs-volume>','yellow'), end = '\n', flush=True)
print (colored('IT WILL BE READ BY ELASTIC SCRIPTS FOR POSTPROCESSING eV-->Ha; A^3-->Bohr^3','yellow'), end = '\n', flush=True)
-
- file = open("energy-vs-volume",'w')
- for i in range(count):
- file.write ("%14.6f %14.6f\n" %(vol_cell[i] * Ang32Bohr3, E[i] * eV2Hartree))
- file.close()
-
- file = open("energy-vs-strain",'w')
- for i in range(count):
- file.write ("%12.6f %14.6f\n" %(strain_value[i], E[i] * eV2Hartree))
- file.close()
+
+ with open("energy-vs-volume",'w') as file:
+ for i in range(count):
+ file.write ("%14.6f %14.6f\n" %(vol_cell[i] * Ang32Bohr3, E[i] * eV2Hartree))
+ with open("energy-vs-strain",'w') as file:
+ for i in range(count):
+ file.write ("%12.6f %14.6f\n" %(strain_value[i], E[i] * eV2Hartree))
###
'''
@@ -90,21 +89,22 @@ def energy_vs_volume():
def fitting_energy_vs_volume_curve():
from sys import stdin
import matplotlib.pyplot as plt
- import matplotlib.ticker as ptk
+ import matplotlib.ticker as ptk
import pylab as pyl
import matplotlib.style
-
+
bohr_radius = 0.529177
bohr32ang3 = 0.14818474347
joule2hartree = 4.3597482
joule2rydberg = joule2hartree/2.
unitconv = joule2hartree/bohr_radius**3*10.**3
-
+
if (str(os.path.exists('energy-vs-volume'))=='False'):
sys.exit("ERROR: file energy-vs-volume not found!\n")
- energy = []; volume = []
+ energy = []
+ volume = []
read_energy = open('energy-vs-volume',"r")
-
+
while True:
line = read_energy.readline()
line = line.strip()
@@ -123,23 +123,33 @@ def fitting_energy_vs_volume_curve():
print ("0 --> Others")
print ("===============================\n")
scheck = input("Enter lattice symmetry code [default 0] >>>> ").replace(" ", "")
-
+
'''
These factors are for the conversion from the conventional cell to primitive cell
BCC: (a^3)/2 primitive cell volume
FCC: (a^3)/4 primitive cell volume
'''
-
- isym = 0; factor = 1
- if ( scheck == "1" ): isym = 1 ; factor=1 ; slabel = "(sc) "
- if ( scheck == "2" ): isym = 2 ; factor=2 ; slabel = "(bcc)"
- if ( scheck == "3" ): isym = 3 ; factor=4 ; slabel = "(fcc)"
+
+ isym = 0
+ factor = 1
+ if scheck == "1":
+ isym = 1
+ factor=1
+ slabel = "(sc) "
+ elif scheck == "2":
+ isym = 2
+ factor=2
+ slabel = "(bcc)"
+ elif scheck == "3":
+ isym = 3
+ factor=4
+ slabel = "(fcc)"
print ("Verification lattice symmetry code >>>> %d " %(isym) )
#-------------------------------------------------------------------------------
- print ('%s' %('-'*105) )
+ print(f"{'-' * 105}")
print ('%20.25s %29.30s %21.30s %11.12s %18.30s' %("Opt_vol Bohr^3 (Ang^3)", "Lattice_const Bohr (A)", "Bulk_modulus [GPa]", "Log(chi)", "Polynomial_order"))
- print ('%s' %('-'*105) )
+ print(f"{'-' * 105}")
for order_of_fit in range(2, 11): #order of polynomial fitting
if order_of_fit % 2 == 0:
order_of_fit = int(order_of_fit)
@@ -149,16 +159,18 @@ def fitting_energy_vs_volume_curve():
bulk = np.poly1d(np.polyder(fitr,2))
bpri = np.poly1d(np.polyder(fitr,3)) #
vmin = np.roots(np.polyder(fitr))
- dmin=[]
- for i in range(len(vmin)):
- if (abs(vmin[i].imag) < 1.e-10):
- if (volume[0] <= vmin[i] and vmin[i] <= volume[-1]):
- if(bulk(vmin[i]) > 0): dmin.append(vmin[i].real)
-
+ dmin = [
+ vmin[i].real
+ for i in range(len(vmin))
+ if (abs(vmin[i].imag) < 1.0e-10)
+ and volume[0] <= vmin[i] <= volume[-1]
+ and (bulk(vmin[i]) > 0)
+ ]
+
xvol = np.linspace(volume[0],volume[-1],100)
if (len(dmin) > 1): print ("WARNING: Multiple minima are found!\n")
- if (len(dmin) == 0): print ("WARNING: No minimum in the given xrange!\n")
-
+ if not dmin: print ("WARNING: No minimum in the given xrange!\n")
+
chi = 0
for i in range(len(energy)):
chi=chi+(energy[i]-curv(volume[i]))**2
@@ -174,14 +186,15 @@ def fitting_energy_vs_volume_curve():
print("%12.6f (%11.6f) %12.6f (%9.6f) %17.6f %13.2f %10d\n" %(v0, v0*bohr32ang3, a0, a0*bohr_radius, b0, math.log10(chi), order_of_fit), end="")
else:
print("%12.6f(%12.6f) %12.6f(%12.6f) %17.6f %13.2f %10d\n" %(v0, v0*bohr32ang3, a0, a0*bohr_radius, b0, math.log10(chi), order_of_fit), end="")
- print ('%s' %('-'*105) )
+ print(f"{'-' * 105}")
#-------------------------------------------------------------------------------
- xlabel = u'Volume [Bohr\u00B3]'; ylabel = r'Energy [Ha]'
+ xlabel = u'Volume [Bohr\u00B3]'
+ ylabel = r'Energy [Ha]'
fontlabel=20
fonttick=14
-
+
params = {'ytick.minor.size': 6,
'xtick.major.pad': 8,
'ytick.major.pad': 4,
@@ -190,15 +203,15 @@ def fitting_energy_vs_volume_curve():
'lines.linewidth': 1.8,
'lines.markersize': 8.0,
'axes.formatter.limits': (-4, 6)}
-
+
plt.rcParams.update(params)
plt.subplots_adjust(left=0.21, right=0.93,
bottom=0.18, top=0.88,
wspace=None, hspace=None)
-
+
yfmt = ptk.ScalarFormatter(useOffset=True,useMathText=True)
-
- figure = plt.figure(1, figsize=(9,9))
+
+ figure = plt.figure(1, figsize=(9,9))
ax = figure.add_subplot(111)
ax.text(0.5,-0.18,xlabel,size=fontlabel, transform=ax.transAxes,ha='center',va='center')
ax.text(-0.25,0.5,ylabel,size=fontlabel, transform=ax.transAxes,ha='center',va='center',rotation=90)
@@ -208,11 +221,11 @@ def fitting_energy_vs_volume_curve():
plt.xticks(size=fonttick)
plt.yticks(size=fonttick)
pyl.grid(True)
- plt.plot(xvol,curv(xvol),'b-',label='n='+str(order_of_fit)+' fit')
+ plt.plot(xvol, curv(xvol), 'b-', label=f'n={str(order_of_fit)} fit')
plt.plot(volume,energy,'go',label='calculated')
plt.plot(dmin,curv(dmin),'ro')
plt.legend(loc=9,borderaxespad=.8,numpoints=1)
-
+
ymax = max(max(curv(xvol)),max(energy))
ymin = min(min(curv(xvol)),min(energy))
dxx = abs(max(xvol)-min(xvol))/18
@@ -220,19 +233,20 @@ def fitting_energy_vs_volume_curve():
ax.yaxis.set_major_formatter(yfmt)
ax.set_xlim(min(xvol)-dxx,max(xvol)+dxx)
ax.set_ylim(ymin-dyy,ymax+dyy)
-
+
ax.xaxis.set_major_locator(MaxNLocator(7))
-
+
plt.savefig('PLOT.png',orientation='portrait',format='png')
###
def sortvolume(s,e):
- ss=[]; ee=[]; ww=[]
- for i in range(len(s)): ww.append(s[i])
- ww.sort()
- for i in range(len(s)):
- ss.append(s[s.index(ww[i])])
- ee.append(e[s.index(ww[i])])
- return ss, ee
+ ss=[]
+ ee=[]
+ ww = [s[i] for i in range(len(s))]
+ ww.sort()
+ for i in range(len(s)):
+ ss.append(s[s.index(ww[i])])
+ ee.append(e[s.index(ww[i])])
+ return ss, ee
###
\ No newline at end of file
diff --git a/src/poscar.py b/src/poscar.py
index 17fb3f6..95aa3da 100644
--- a/src/poscar.py
+++ b/src/poscar.py
@@ -14,9 +14,12 @@ def poscar():
print (' ERROR: POSCAR does not exist here.')
sys.exit(0)
print('Reading POSCAR/CONTCAR: \n')
- pos = []; kk = []; lattice = []; sum = 0
+ pos = []
+ kk = []
+ lattice = []
+ sum = 0
file = open('POSCAR','r') or open('CONTCAR','r')
-
+
firstline = file.readline() # IGNORE first line comment
secondfline = file.readline() # scale
Latvec1 = file.readline()
@@ -28,11 +31,10 @@ def poscar():
elementtype=file.readline().split()
if (str.isdigit(elementtype[0])):
sys.exit("VASP 4.X POSCAR detected. Please add the atom types")
- print ("Types of elements:", str(elementtype), end = '\n')
+ print("Types of elements:", elementtype, end = '\n')
numberofatoms=file.readline()
Coordtype=file.readline()
print ("Coordtype:", (Coordtype), end = '\n')
-
########################---------------------------------------------------------
print (">>>>>>>>>-------------------# of Atoms--------------------")
nat = numberofatoms.split()
@@ -44,37 +46,31 @@ def poscar():
print ("Number of atoms:", (numberofatoms), end = '\n')
########################---------------------------------------------------------
- #print (">>>>>>>>>---------------Atomic positions------------------")
- for x in range(int(numberofatoms)):
+ #print (">>>>>>>>>---------------Atomic positions------------------")
+ for _ in range(int(numberofatoms)):
coord = file.readline().split()
coord = [float(i) for i in coord]
pos = pos + [coord]
pos = np.array(pos)
#print (pos)
-
+
file.close()
-########################---------------------------------------------------------
- a=[]; b=[]; c=[];
Latvec1=Latvec1.split()
Latvec2=Latvec2.split()
- Latvec3=Latvec3.split()
- for ai in Latvec1:
- a.append(float(ai))
- for bi in Latvec2:
- b.append(float(bi))
- for ci in Latvec3:
- c.append(float(ci))
-
+ Latvec3=Latvec3.split()
+ a = [float(ai) for ai in Latvec1]
+ b = [float(bi) for bi in Latvec2]
+ c = [float(ci) for ci in Latvec3]
########################---------------------------------------------------------
- print (">>>>>>>>>---------------Lattice vectors distortions-----------------")
+ print (">>>>>>>>>---------------Lattice vectors distortions-----------------")
lattice = np.array([a] + [b] + [c])
#determinant = np.linalg.det(lattice)
lld = local_lattice_distortion(a,b,c)
- print ("lattice distortion parameter g: {}".format(lld) )
-
- print (">>>>>>>>>---------------Space group-----------------")
+ print(f"lattice distortion parameter g: {lld}")
+
+ print (">>>>>>>>>---------------Space group-----------------")
print (" ")
sp, symm = space_group_analyse(lattice, pos)
print ( sp, symm )
@@ -99,12 +95,14 @@ def local_lattice_distortion(a1,b1,c1):
#print (Back.YELLOW + "Wang, S. Atomic structure modeling of multi-principal-element alloys by the principle")
#print (Back.YELLOW + "of maximum entropy. Entropy 15, 5536–5548 (2013).")
#print ("")
- a=np.linalg.norm(a1); b=np.linalg.norm(b1); c=np.linalg.norm(c1)
+ a=np.linalg.norm(a1)
+ b=np.linalg.norm(b1)
+ c=np.linalg.norm(c1)
d = np.array([a,b,c])
- d_mean = np.mean(d); d_std = np.std(d)
+ d_mean = np.mean(d)
+ d_std = np.std(d)
d_square_mean = (a**2 + b**2 + c**2)/3
- g = np.sqrt( d_square_mean/(d_mean)**2 - 1 )
- return g
+ return np.sqrt( d_square_mean/(d_mean)**2 - 1 )
###
# Song, H. et al. Local lattice distortion in high-entropy alloys. Phys. Rev. Mater. 1, 23404 (2017).
# Senkov, O. N. & Miracle, D. B. Effect of the atomic size distribution on glass forming ability of amorphous metallic alloys. Mater. Res. Bull. 36, 2183–2198 (2001).
@@ -115,7 +113,9 @@ def local_lattice_distortion_DEF1():
#print (Back.YELLOW + "Wang, S. Atomic structure modeling of multi-principal-element alloys by the principle")
#print (Back.YELLOW + "of maximum entropy. Entropy 15, 5536–5548 (2013).")
print ("+"*40,"HUME ROTHERY RULE","+"*40)
- C_i=C=0.2 ; r_avg = 0.0; del_sum=0.0
+ C_i=C=0.2
+ r_avg = 0.0
+ del_sum=0.0
elements = ["Nb", "Hf", "Ta", "Ti", "Zr"]
eta = {
"Nb" : 1.98,
@@ -123,17 +123,17 @@ def local_lattice_distortion_DEF1():
"Ta" : 2.00,
"Ti" : 1.76,
"Zr" : 2.06, }
-
+
print (" {element: atomic radius}")
print (eta)
-
+
for i in elements:
r_avg = r_avg + C * eta[i]
-
+
for j in elements:
del_sum = del_sum + C * ( 1 - float(eta[j]) / r_avg )**2
- del_sum = 100 * np.sqrt(del_sum)
- print("HEA_atomic_size_mismatch: \u03B4={}".format(del_sum))
+ del_sum = 100 * np.sqrt(del_sum)
+ print(f"HEA_atomic_size_mismatch: \u03B4={del_sum}")
###
def local_lattice_distortion_DEF2():
@@ -141,34 +141,35 @@ def local_lattice_distortion_DEF2():
print ("Phys. Rev. Mater. 1, 23404 (2017).")
print ("_____| Different definition of the atomic radius for the description ")
print (" of the local lattice distortion in HEAs")
-
+
if not os.path.exists('POSCAR' and 'CONTCAR'):
print (' ERROR: POSCAR & CONTCAR does not exist')
sys.exit(0)
print('Reading POSCAR and CONTCAR ... \n')
-
- x = []; y =[]; z=[]
- xp =[]; yp = []; zp = []; temp=0
-
- f = open('POSCAR','r')
- lines_poscar = f.readlines()
- f.close()
-
- f = open('CONTCAR','r')
- lines_contcar = f.readlines()
- f.close()
-
+
+ x = []
+ y =[]
+ z=[]
+ xp =[]
+ yp = []
+ zp = []
+ temp=0
+
+ with open('POSCAR','r') as f:
+ lines_poscar = f.readlines()
+ with open('CONTCAR','r') as f:
+ lines_contcar = f.readlines()
sum_atoms = lines_poscar[6].split() ### reading 7th lines for reading # of atoms
sum_atoms = [int(i) for i in sum_atoms]
sum_atoms = sum(sum_atoms)
-
+
for i in lines_poscar:
if "Direct" in i:
lp=lines_poscar.index(i)
for j in lines_contcar:
if "Direct" in j:
lc=lines_contcar.index(j)
-
+
for i in range(sum_atoms):
x, y, z = lines_poscar[lp+1+i].split()
xp, yp, zp = lines_contcar[lp+1+i].split()
@@ -176,7 +177,7 @@ def local_lattice_distortion_DEF2():
xp = float(xp); yp = float(yp); zp = float(zp)
temp = temp + np.sqrt( (x-xp)**2 + (y-yp)**2 + (z-zp)**2 )
temp = temp/sum_atoms
- print("local lattice distortion: \u0394d={}".format(temp))
+ print(f"local lattice distortion: \u0394d={temp}")
###
def space_group_analyse(lattice, pos):
@@ -214,50 +215,44 @@ def space_group_analyse(lattice, pos):
def poscar_VASP42VASP5():
if not os.path.exists('POSCAR' and 'POTCAR'):
print (' ERROR: POSCAR does not exist here.')
- sys.exit(0)
- file1 = open("POSCAR",'r')
- line1 = file1.readlines()
- file1.close()
-
- file2 = open("POTCAR",'r')
- line2 = file2.readlines()
- file2.close()
-
+ sys.exit(0)
+ with open("POSCAR",'r') as file1:
+ line1 = file1.readlines()
+ with open("POTCAR",'r') as file2:
+ line2 = file2.readlines()
atom_number=[]
for i in line1:
if ("Direct" or "direct" or "d" or "D") in i:
PP=line1.index(i)
atom_number = line1[5].split()
print(atom_number)
-
- elementtype=[]; count=0
+
+ elementtype=[]
+ count=0
for i in line2:
if ("VRHFIN") in i:
count+=1
#print (i.split('=')[1].split(':')[0])
elementtype.append(i.split('=')[1].split(':')[0])
-
- test = open("POSCAR_W",'w')
-
- for i in range(5):
- test.write( line1[i] )
-
- for j in elementtype:
- test.write("\t" + j)
- test.write("\n" )
-
- for j in atom_number:
- test.write("\t" + j )
- test.write("\n" )
-
- test.write("Selective dynamics")
- test.write("\n" )
-
- for i in range(len(line1)-PP):
- test.write(line1[PP+i] )
-
- test.close()
-
+
+ with open("POSCAR_W",'w') as test:
+ for i in range(5):
+ test.write( line1[i] )
+
+ for j in elementtype:
+ test.write("\t" + j)
+ test.write("\n" )
+
+ for j in atom_number:
+ test.write("\t" + j )
+ test.write("\n" )
+
+ test.write("Selective dynamics")
+ test.write("\n" )
+
+ for i in range(len(line1)-PP):
+ test.write(line1[PP+i] )
+
print (" File is converted: POSCAR_W")
###
@@ -267,7 +262,7 @@ def poscar_VASP42VASP5():
def volume(a,b,c,alpha,beta,gamma):
ang2atomic = 1.889725988579 # 1 A = 1.889725988579 [a.u]
Ang32Bohr3 = 6.74833304162 # 1 A^3 = 6.7483330371 [a.u]^3
-
+
length = np.linalg.norm(a) * np.linalg.norm(b) * np.linalg.norm(c)
volume = length * ( np.sqrt(1 + 2 * math.cos(alpha) * math.cos(beta) * math.cos(gamma) - math.cos(alpha)**2 - math.cos(beta)**2 - math.cos(gamma)**2) )
vol_au = volume * Ang32Bohr3