Last active
February 22, 2024 15:00
-
-
Save dceresoli/a83846043c970e5dfc385a02eb1b8f19 to your computer and use it in GitHub Desktop.
Fit polynomial Equations of State (EoS) to (lattice,enery) or (volume,energy) data
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