Created
May 5, 2020 07:16
-
-
Save sriharijayaram5/62f3ba4fcf1f3a1585776b78324df45a to your computer and use it in GitHub Desktop.
EOS Poly fit
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# Fit polynomial equations of state (EOS) to energy vs volume curves | |
# and optionally plot the results | |
# Davide Ceresoli <dceresoli@gmail.com>, 10/03/2012 | |
# Inspired by ase.utils.eosase2 | |
import sys, numpy, math | |
from scipy.optimize import leastsq | |
import numpy.polynomial.polynomial as poly | |
# Customized input with default and accepted values | |
def myinput(prompt, default, accepted): | |
while True: | |
res = input(prompt + " [default=%s]: " % (default)) | |
if res == '': res = default | |
if res in accepted: | |
break | |
else: | |
print("accepted values:", accepted) | |
return res | |
print("Welcome to eos-fit-poly.py") | |
print() | |
fname = input("Filename containing energy vs volume [volume.dat]: ") | |
if fname == '': fname = "volume.dat" | |
try: | |
f = open(fname, 'rt') | |
except IOError: | |
sys.stderr.write("Error opening or reading file %s\n" % (fname)) | |
sys.exit(1) | |
# read data from file | |
print() | |
print("Data read from file:") | |
vol = [] | |
ene = [] | |
while True: | |
line = f.readline().strip() | |
if line == '': break | |
if line[0] == '#' or line[0] == '!': continue | |
v, e = [float(x) for x in line.split()[:2]] | |
vol.append(v) | |
ene.append(e) | |
print(v, e) | |
print() | |
f.close() | |
# transform to numpy arrays | |
vol = numpy.array(vol) | |
ene = numpy.array(ene) | |
# further questions | |
is_volume = True | |
ans = myinput("Volume or lattice spacing (vol/latt)", "vol", ["vol", "latt"]) | |
is_volume = (ans == 'vol') | |
if is_volume: | |
vol_unit = myinput("Volume units (ang3/bohr3)", "ang3", ["ang3", "bohr3"]) | |
if vol_unit == "bohr3": vol *= 0.14818471 # convert to angstrom^3 | |
else: | |
vol_unit = myinput("Lattice units (ang/bohr)", "ang", ["ang", "bohr"]) | |
if vol_unit == "bohr": vol *= 0.52917721 # convert to angstrom | |
latt = myinput("Lattice type (sc/fcc/bcc)", "sc", ["sc", "bcc", "fcc"]) | |
fact = 1.0 | |
if latt == "fcc": fact = 0.25 | |
if latt == "bcc": fact = 0.5 | |
# lattice to volume | |
vol = fact * vol**3 | |
ene_unit = myinput("Energy units (eV/Ry/Ha)", "eV", ["eV", "Ry", "Ha"]) | |
if ene_unit == "Ry": ene *= 13.605692 # convert to eV | |
if ene_unit == "Ha": ene *= 27.211383 # convert to eV | |
order = input("Order of polynomial [4]: ") | |
try: | |
order = int(order) | |
assert 2 < order < len(vol) | |
except: | |
print("wrong input, using 4-th order") | |
order = 4 | |
print() | |
# fit a sixth-order polynomial to the data | |
p = poly.polyfit(vol, ene, order) | |
p1 = poly.polyder(p) # first derivative | |
roots = poly.polyroots(p1) | |
p2 = poly.polyder(p1) # second derivative | |
p3 = poly.polyder(p2) # third derivative | |
# check for roots in the interval | |
label = "%ith-order" % (order) | |
print() | |
for r in roots: | |
if r.imag == 0.0 and vol[0] <= r.real <= vol[-1]: | |
V0 = r.real | |
E0 = poly.polyval(V0, p) | |
B0 = poly.polyval(V0, p2) * V0 | |
Bp = -V0*V0/B0 * poly.polyval(V0, p3) | |
print(label, ": E0 = %f eV" % (E0)) | |
print(label, ": B0 = %f GPa" % (B0*160.21765)) | |
print(label, ": Bp = %f" % (Bp)) | |
print(label, ": V0 = %f angstrom^3" % (V0)) | |
print() | |
try: | |
import pylab | |
except ImportError: | |
sys.stderr.write("pylab module non available, skipping plot") | |
sys.exit(0) | |
# plotting | |
ans = myinput("Do you want to plot the result (yes/no)", "yes", ["yes", "no"]) | |
if ans == "no": sys.exit(0) | |
#import pylab | |
vfit = numpy.linspace(min(vol),max(vol),100) | |
pylab.plot(vol, ene-E0, 'ro') | |
pylab.plot(vfit, poly.polyval(vfit, p)-E0, label=label) | |
pylab.xlabel('Volume ($\AA^3$)') | |
pylab.ylabel('Energy (eV)') | |
pylab.legend(loc='best') | |
pylab.show() | |
quit() | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment