Проблема с pylab при построении графика ⇐ Python

[Anonymous]

Код: Выделить всё

I am trying to fit two curves (idential function) to my two sets of data and plot the two fiited curves and two sets of data together using python. Following is my script whcich uses two data files volume-1.dat and volume-pp-fe.dat.

====================================== ==============
#!/usr/bin/env python
# Сопоставьте уравнения состояния (EOS) с кривыми энергии и объема и, при необходимости, постройте график результатов

Код: Выделить всё

import sys, numpy, math
from scipy.optimize import leastsq

# Birch-Murnaghan equation of state-set1
def eos_birch_murnaghan_1(params, vol):
'From Phys. Rev. B 70, 224107'
E1, B1, Bp1, V1 = params
eta1 = (V1/vol_1)**(1.0/3.0)
E = E1 + 9.0*B1*V1/16.0 * (eta1**2-1.0)**2 * (6.0 + Bp1*(eta1**2-1.0) - 4.0*eta1**2)
return E

# Birch-Murnaghan equation of state-set2
def eos_birch_murnaghan_2(params, vol):
'From Phys. Rev.  B 70, 224107'
E2, B2, Bp2, V2 = params
eta2 = (V2/vol_1)**(1.0/3.0)
E = E2 + 9.0*B2*V2/16.0 * (eta2**2-1.0)**2 * (6.0 + Bp2*(eta2**2-1.0) - 4.0*eta2**2)
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

print("Welcome to eos-fit.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_1 = []
ene_1 = []
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_1.append(v)
ene_1.append(e)
print(v, e)
print()
f.close()

# transform to numpy arrays
vol_1 = numpy.array(vol_1)
ene_1 = numpy.array(ene_1)

# 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.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_2 = []
ene_2 = []
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_2.append(v)
ene_2.append(e)
print(v, e)
print()
f.close()

# transform to numpy arrays
vol_2 = numpy.array(vol_2)
ene_2 = numpy.array(ene_2)

# fit a parabola to the data and get inital guess for equilibirum volume
# and bulk modulus
a1, b1, c1 = numpy.polyfit(vol_1, ene_1, 2)
V1 = -b1/(2*a1)
E1 = a1*V1**2 + b1*V1 + c1
B1 = 2*a1*V1
Bp1 = 4.0

# initial guesses in the same order used in the Murnaghan function
x1 = [E1, B1, Bp1, V1]

def print_params(label, params):
E1, B1, B1, V1 = params
print(label, ": E1 = %f eV" % (E1))
print(label, ": B1 = %f GPa" % (B1*160.21765))
print(label, ": Bp1 = %f" % (Bp1))
if is_volume:
print(label, ": V1 = %f angstrom^3" % (V1))
else:
print(label, ": V1 = %f angstrom^3, a0 = %f angstrom" % (V1, (V1/fact)**(1./3.)))
print()

# fit the equations of state
target = lambda params, y, x: y - eos_birch_murnaghan_1(params, x)
birch_murn_1, ier = leastsq(target, x1, args=(ene_1,vol_1))
print_params("Birch-Murnaghan_1", birch_murn_1)

# fit a parabola to the data and get inital guess for equilibirum volume
# and bulk modulus
a2, b2, c2 = numpy.polyfit(vol_2, ene_2, 2)
V2 = -b2/(2*a2)
E2 = a2*V2**2 + b2*V2 + c2
B2 = 2*a2*V2
Bp2 = 4.0

# initial guesses in the same order used in the Murnaghan function
x2 = [E2, B2, Bp2, V2]

def print_params(label, params):
E2, B2, Bp2, V2 = params
print(label, ": E2 = %f eV" % (E2))
print(label, ": B2 = %f GPa" % (B2*160.21765))
print(label, ": Bp2 = %f" % (Bp2))
if is_volume:
print(label, ": V2 = %f angstrom^3" % (V2))
else:
print(label, ": V2 = %f angstrom^3, a0 = %f angstrom"  % (V2, (V2/fact)**(1./3.)))
print()

# fit the equations of state
target = lambda params, y, x: y - eos_birch_murnaghan_1(params, x)
birch_murn_2, ier = leastsq(target, x2, args=(ene_2,vol_2))
print_params("Birch-Murnaghan_1", birch_murn_2)

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_1 = numpy.linspace(min(vol_1),max(vol_1),100)
vfit_2 = numpy.linspace(min(vol_2),max(vol_2),100)

pylab.plot(vol_1, ene_1, 'ro')
pylab.plot(vfit_1, eos_birch_murnaghan_1(birch_murn_1,vfit_1), label='Birch-Murnaghan-1')
pylab.plot(vol_2, ene_2, 'bo')
pylab.plot(vfit_2, eos_birch_murnaghan_2(birch_murn_2,vfit_2), label='Birch-Murnaghan-2')
pylab.xlabel('Volume ($\AA^3$)')
pylab.ylabel('Energy (eV)')
pylab.legend(loc='best')
pylab.show()
quit()

Код: Выделить всё

**On runnning it gives value error. Here is its copy:**
=============================================
File "/lfs/sudipkm/fcc-fe/PP-Fe_sv/relax/alat/step-2/bmeos-poly-ev.py", line 113, in 
plot(vfit_1, eos_birch_murnaghan_1(birch_murn_1,vfit_1), 'red', label='BMEOS-zval-16')
File "/home/sudipkm/.conda/envs/phonopy/lib/python3.12/site-packages/matplotlib/pyplot.py", line 3794, in plot
return gca().plot(
^^^^^^^^^^^
File "/home/sudipkm/.conda/envs/phonopy/lib/python3.12/site-packages/matplotlib/axes/_axes.py", line 1779, in plot
lines = [*self._get_lines(self, *args, data=data, **kwargs)]
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/sudipkm/.conda/envs/phonopy/lib/python3.12/site-packages/matplotlib/axes/_base.py", line 296, in __call__
yield from self._plot_args(
^^^^^^^^^^^^^^^^
File "/home/sudipkm/.conda/envs/phonopy/lib/python3.12/site-packages/matplotlib/axes/_base.py", line 486, in _plot_args
raise ValueError(f"x and y must have same first dimension, but "
ValueError: x and y must have same first dimension, but have shapes (100,) and (9,)

Код: Выделить всё

**However if I run them singly, I get no error. Meaning if I run the foloowing code separately on the two .dat files it gives no error.**

====================================== =============================
#!/usr/bin/env python
# Подходящие уравнения состояния (EOS) для кривых зависимости энергии от объема и, при необходимости, построения графика результаты

Код: Выделить всё

import sys, numpy, math
from scipy.optimize import leastsq

# Birch-Murnaghan equation of state
def eos_birch_murnaghan(params, vol):
'From Phys. Rev.  B 70, 224107'
E0, B0, Bp, V0 = params
eta = (V0/vol)**(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

# 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.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)

# fit a parabola to the data and get inital guess for equilibirum volume
# and bulk modulus
a, b, c = numpy.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))
if is_volume:
print(label, ": V0 = %f angstrom^3" % (V0))
else:
print(label, ": V0 = %f angstrom^3, a0 = %f angstrom" % (V0, (V0/fact)**(1./3.)))
print()

# fit the equations of state
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)

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, 'ro')
pylab.plot(vfit, eos_birch_murnaghan(birch_murn,vfit), label='Birch-Murnaghan')
pylab.xlabel('Volume ($\AA^3$)')
pylab.ylabel('Energy (eV)')
pylab.legend(loc='best')
pylab.show()
quit()

Кажется, я не могу понять проблему, поскольку оба кода идентичны. Есть идеи, что является причиной этого? Я понимаю, что значения x и y не имеют длины smae. Но почему первый код выдает ошибку, а второй работает нормально? vfit-linspace не работает должным образом в первом случае, а во втором — работает. Спасибо большое.

Код: Выделить всё

Folllowing are the data files I have.
volume-1.dat
=============================================
#V  #E
30.96   -29.19558780
31.55   -29.72428941
32.16   -30.20134821
32.77   -30.62963030
33.39   -31.01182951
34.01   -31.35042683
34.65   -31.64781786
35.29   -31.90625727
35.94   -32.12787409

Код: Выделить всё

and volume-pp-fe.dat
==============================================
#V  #E
30.96   -28.67254616
31.55   -29.20333191
32.16   -29.68196612
32.77   -30.11133108
33.39   -30.49416679
34.01   -30.83309362
34.65   -31.13061068
35.29   -31.38909654
35.94   -31.61080325

Подробнее здесь: https://stackoverflow.com/questions/793 ... e-plotting