-
Notifications
You must be signed in to change notification settings - Fork 0
/
bulk_modulus_fit.py
executable file
·178 lines (145 loc) · 5.61 KB
/
bulk_modulus_fit.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
# non-linear fitting/regression
import sys
import math
import numpy as np
from scipy.optimize import leastsq
#####################################
# Birch-Murnaghan equation of state #
#####################################
def eos_birch_murnaghan(params, vol):
'From Phys. Rev. B 70, 224107'
'It is third-order Birch-Murnaghan EOS'
E0, B0, Bp, V0 = params
eta = (vol/V0)**(1.0/3.0)
E = E0 + 9.0*B0*V0/16.0 * (eta**2-1.0)**2 * (6.0 + Bp*(eta**2-1.0) - 4.0*eta**2)
return E
# Murnaghan equation of state
def eos_murnaghan(params, vol):
'From Phys. Rev. B 28, 5480 (1983)'
E0, B0, Bp, V0 = params
E = E0 + B0/Bp * vol * ((V0/vol)**Bp/(Bp-1.0)+1.0) - V0*B0/(Bp-1.0)
return E
# Birch equation of state
def eos_birch(params, vol):
'From Intermetallic compounds: Principles and Practice, Vol. I: Princples'
'Chapter 9 pages 195-210 by M. Mehl. B. Klein, D. Papaconstantopoulos'
E0, B0, Bp, V0 = params
E = (E0 + 9.0/8.0*B0*V0*((V0/vol)**(2.0/3.0) - 1.0)**2
+ 9.0/16.0*B0*V0*(Bp-4.)*((V0/vol)**(2.0/3.0) - 1.0)**3)
return E
# Vinet equation of state
def eos_vinet(params, vol):
'From Phys. Rev. B 70, 224107'
E0, B0, Bp, V0 = params
eta = (vol/V0)**(1.0/3.0)
E = E0 + 2.0*B0*V0/(Bp-1.0)**2 * \
(2.0 - (5.0 + 3.0*Bp*(eta-1.0) - 3.0*eta)*np.exp(-3.0*(Bp-1.0)*(eta-1.0)/2.0))
return E
# 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
# Note!
# If you want to run this code in your own server please delete the comment (''') in below
# and eliminate the prepared numpy array (right below of the #further questions)
print ("This is EOS-fit.py")
'''
#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 = np.array(vol)
#ene = np.array(ene)
'''
# further questions
'AlF3'
vol = np.array([88.650323, 88.597737, 88.545163, 88.808215, 88.860884])
ene = np.array([-29486.984323034, -29486.984302681, -29486.984253874, -29486.984230923, -29486.984149340])
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
print
# fit a parabola to the data and get inital guess for equilibirum volume
# and bulk modulus
a, b, c = np.polyfit(vol, ene, 2)
V0 = -b/(2*a)
E0 = a*V0**2 + b*V0 + c
B0 = 2*a*V0
Bp = 4.0
# initial guesses in the same order used in the Murnaghan function
x0 = [E0, B0, Bp, V0]
def print_params(label, params):
E0, B0, Bp, V0 = params
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
# fit the equations of state
target = lambda params, y, x: y - eos_murnaghan(params, x)
murn, ier = leastsq(target, x0, args=(ene,vol))
print_params("Murnaghan", murn)
target = lambda params, y, x: y - eos_birch_murnaghan(params, x)
birch_murn, ier = leastsq(target, x0, args=(ene,vol))
print_params("Birch-Murnaghan", birch_murn)
target = lambda params, y, x: y - eos_birch(params, x)
birch, ier = leastsq(target, x0, args=(ene,vol))
print_params("Birch", birch)
target = lambda params, y, x: y - eos_vinet(params, x)
vinet, ier = leastsq(target, x0, args=(ene,vol))
print_params("Vinet", vinet)
# plotting
ans = myinput("Do you want to plot the result (yes/no)", "yes", ["yes", "no"])
if ans == "no": sys.exit(0)
import pylab
vfit = np.linspace(min(vol),max(vol),100)
pylab.plot(vol, ene-E0, 'ro')
pylab.plot(vfit, eos_murnaghan(murn,vfit)-E0, label='Murnaghan')
pylab.plot(vfit, eos_birch_murnaghan(birch_murn,vfit)-E0, label='Birch-Murnaghan')
pylab.plot(vfit, eos_birch(birch,vfit)-E0, label='Birch')
pylab.plot(vfit, eos_vinet(vinet,vfit)-E0, label='Vinet')
pylab.xlabel('Volume ($\AA^3$)')
pylab.ylabel('Energy (eV)')
pylab.legend(loc='best')
pylab.show()
quit()