View source for EOS comp plot.py

This code plots the energy vs lattice parameter curve from DFT and EOS fitting for multiple equations of state and compares the difference between each equation of state to values found directly from DFT from [[Code: Quantum Espresso]] simulations. 

This code is part of a group of five codes all used for post-processing data from Quantum Espresso simulations.
* The [[EvA EvV plot.py]] code generates plots for the energy vs volume and energy vs lattice parameter curves for any set of EvsA an EvsV files generated by the ev_curve.bash script.
* The [[Convergence plots.py]] code generates plots for the convergency study from the <code> SUMMARY</code> files.
* The [[EOS plot.py]] code plots the difference in the energy values found from DFT and from the chosen EOS.
* The [[EOS comp plot.py]] code plots the energy vs lattice parameter curve from DFT and EOS fitting for multiple equations of state and compares the difference between each equation of state to values found directly from DFT.
* The [[Ecut conv.py]] code plots the values found from the ecut convergence studies.

{|border  ="0"
|<pre>

# For Parsing Esv files and evfit files into images
import numpy as np
import matplotlib.pyplot as plt
import os

# The values are stored in 2 Columns as:
# Lattice Parameter(angstrom)       energy(eV)
path = os.getcwd()


#setting up path to needed files##
##################################

#for Narrow run
esva_n1 = open(path+'/narrow/1/EvsA') #should point to narrow EvsA.txt file in 1
esvv_n1 = open(path+'/narrow/1/EvsV') #should point to narrow EvsV.txt file in 1 
evfit_n1 = open(path+'/narrow/1/evfit.1') #should point to narrow evfit file in 1
SUMn1   = open(path+'/narrow/1/SUMMARY')#should point to narrow summary file in 1 

esva_n4 = open(path+'/narrow/4/EvsA') #should point to narrow EvsA.txt file in 4
esvv_n4 = open(path+'/narrow/4/EvsV') #should point to narrow EvsV.txt file in 4 
evfit_n4 = open(path+'/narrow/4/evfit.4') #should point to narrow evfit file in 4
sum_n4  = open(path+'/narrow/4/SUMMARY')#should point to narrow summary file in 4

#for Medium run
esva_m1 = open(path+'/med/1/EvsA') #should point to med EvsA.txt file in 1
esvv_m1 = open(path+'/med/1/EvsV') #should point to med EvsV.txt file in 1 
evfit_m1 = open(path+'/med/1/evfit.1') #should point to med evfit file in 1
sum_m1  = open(path+'/med/1/SUMMARY')#should point to med summary file in 1 

esva_m4 = open(path+'/med/4/EvsA') #should point to med EvsA.txt file in 4
esvv_m4 = open(path+'/med/4/EvsV') #should point to med EvsV.txt file in 4 
evfit_m4 = open(path+'/med/4/evfit.4') #should point to med evfit file in 4
sum_m4  = open(path+'/med/4/SUMMARY')#should point to med summary file in 4

#for wide run
esva_w1 = open(path+'/wide/1/EvsA') #should point to wide EvsA.txt file in 1
esvv_w1 = open(path+'/wide/1/EvsV') #should point to wide EvsV.txt file in 1 
evfit_w1 = open(path+'/wide/1/evfit.1') #should point to wide evfit file in 1
sum_w1  = open(path+'/wide/1/SUMMARY')#should point to wide summary file in 1 

esva_w4 = open(path+'/wide/4/EvsA') #should point to wide EvsA.txt file in 4
esvv_w4 = open(path+'/wide/4/EvsV') #should point to wide EvsV.txt file in 4 
evfit_w4 = open(path+'/wide/4/evfit.4') #should point to wide evfit file in 4
sum_w4  = open(path+'/wide/4/SUMMARY')#should point to wide summary file in 4


#Setting Colors
color1 = '#00035b'
color2 = '#8b2e16'

#Setting font/figure sizes
fig_size_x  = 30
fig_size_y  = 20
title_font  = 60
subtitle_font = 54
mark_size   = 18
tick_font   = 20
legend_size = 30
y_label     = 30
x_label     = 42

#setting marker type
marker1 = 'o'
marker2 = '*'

#setting energy offset
energy_offset = 2858.8298734 # Set by running simulation with very large lattice parameter

#Legend Location
leg_loc = "lower center"


###Initializing arrays### 
#########################

# array of Lattice Parameter in Angrstrom
A_n1, A_n4, A_m1, A_m4, A_w1, A_w4 = ([] for i in range(6))


# array of Volume in Angrstrom^3
A3_n1, A3_n4, A3_m1, A3_m4, A3_w1, A3_w4 = ([] for i in range(6))

# array of Energy in eV from DFT 
E_n1, E_n4, E_m1, E_m4, E_w1, E_w4 = ([] for i in range(6))

#array of Energy in eV from EOS fitting
Ev_eos_n1, Ev_eos_n4, Ev_eos_m1, Ev_eos_m4, Ev_eos_w1, Ev_eos_w4 = ([] for i in range(6))

#array of Error in eV from EOS fitting
Err_n1, Err_n4, Err_m1, Err_m4, Err_w1, Err_w4 = ([] for i in range(6))

##Getting data for all files##
##############################

##########################
# Narrow run with EOS 1  #
##########################
# EvA/EvV data
line = True
while line: 
    line = esva_n1.readline()
    values = line.split()
    if not values:
        break
    A_n1.append(float(values[0]))
    E_n1.append(float(values[1]))
esva_n1.close()

line = True
while line:   
    line = esvv_n1.readline()
    values = line.split()
    if not values:
        break
    A3_n1.append(float(values[0]))
esvv_n1.close()


#evfit data
line = True
while line:  
    line  = evfit_n1.readline()
    lines = line.split()
    if lines != []:
        if lines[0] =="A0=":
            A_min_n1 = float(lines[1])
            K0_n1    = float(lines[3])
            DK0_n1   = float(lines[3])
            D2K0_n1  = float(lines[3])
        if lines[0] =="Emin=":
            Emin_n1  = float(lines[1])
        if (lines[0] != "A0=") and (lines[0] != "Equation") and (lines[0] != "Emin="):
            Ev_eos_n1.append(float(lines[2]))
            Err_n1.append(float(lines[3]))
evfit_n1.close()


##########################
# Narrow run with EOS 4  #
##########################
# EvA/EvV data
line = True
while line: 
    line = esva_n4.readline()
    values = line.split()
    if not values:
        break
    A_n4.append(float(values[0]))
    E_n4.append(float(values[1]))
esva_n4.close()

line = True
while line:   
    line = esvv_n4.readline()
    values = line.split()
    if not values:
        break
    A3_n4.append(float(values[0]))
esvv_n4.close()


#evfit data
line = True
while line:  
    line  = evfit_n4.readline()
    lines = line.split()
    if lines != []:
        if lines[0] =="A0=":
            A_min_n4 = float(lines[1])
            K0_n4    = float(lines[3])
            DK0_n4   = float(lines[3])
            D2K0_n4  = float(lines[3])
        if lines[0] =="Emin=":
            Emin_n4  = float(lines[1])
        if (lines[0] != "A0=") and (lines[0] != "Equation") and (lines[0] != "Emin="):
            Ev_eos_n4.append(float(lines[2]))
            Err_n4.append(float(lines[3]))
evfit_n4.close()

##########################
# Medium run with EOS 1  #
##########################
# EvA/EvV data
line = True
while line: 
    line = esva_m1.readline()
    values = line.split()
    if not values:
        break
    A_m1.append(float(values[0]))
    E_m1.append(float(values[1]))
esva_m1.close()

line = True
while line:   
    line = esvv_m1.readline()
    values = line.split()
    if not values:
        break
    A3_m1.append(float(values[0]))
esvv_m1.close()

#evfit data
line = True
while line:  
    line  = evfit_m1.readline()
    lines = line.split()
    if lines != []:
        if lines[0] =="A0=":
            A_min_m1 = float(lines[1])
            K0_m1    = float(lines[3])
            DK0_m1   = float(lines[3])
            D2K0_m1  = float(lines[3])
        if lines[0] =="Emin=":
            Emin_m1  = float(lines[1])
        if (lines[0] != "A0=") and (lines[0] != "Equation") and (lines[0] != "Emin="):
            Ev_eos_m1.append(float(lines[2]))
            Err_m1.append(float(lines[3]))
evfit_m1.close()


##########################
# Medium run with EOS 4  #
##########################
# EvA/EvV data
line = True
while line: 
    line = esva_m4.readline()
    values = line.split()
    if not values:
        break
    A_m4.append(float(values[0]))
    E_m4.append(float(values[1]))
esva_m4.close()

line = True
while line:   
    line = esvv_m4.readline()
    values = line.split()
    if not values:
        break
    A3_m4.append(float(values[0]))
esvv_m4.close()

#evfit data
line = True
while line:  
    line  = evfit_m4.readline()
    lines = line.split()
    if lines != []:
        if lines[0] =="A0=":
            A_min_m4 = float(lines[1])
            K0_m4    = float(lines[3])
            DK0_m4   = float(lines[3])
            D2K0_m4  = float(lines[3])
        if lines[0] =="Emin=":
            Emin_m4  = float(lines[1])
        if (lines[0] != "A0=") and (lines[0] != "Equation") and (lines[0] != "Emin="):
            Ev_eos_m4.append(float(lines[2]))
            Err_m4.append(float(lines[3]))
evfit_m4.close()



##########################
# Wide run with EOS 1    #
##########################
# EvA/EvV data
line = True
while line: 
    line = esva_w1.readline()
    values = line.split()
    if not values:
        break
    A_w1.append(float(values[0]))
    E_w1.append(float(values[1]))
esva_w1.close()

line = True
while line:   
    line = esvv_w1.readline()
    values = line.split()
    if not values:
        break
    A3_w1.append(float(values[0]))
esvv_w1.close()


#evfit data
line = True
while line:  
    line  = evfit_w1.readline()
    lines = line.split()
    if lines != []:
        if lines[0] =="A0=":
            A_min_w1 = float(lines[1])
            K0_w1    = float(lines[3])
            DK0_w1   = float(lines[3])
            D2K0_w1  = float(lines[3])
        if lines[0] =="Emin=":
            Emin_w1  = float(lines[1])
        if (lines[0] != "A0=") and (lines[0] != "Equation") and (lines[0] != "Emin="):
            Ev_eos_w1.append(float(lines[2]))
            Err_w1.append(float(lines[3]))
evfit_w1.close()


##########################
# Wide run with EOS 4  #
##########################
# EvA/EvV data
line = True
while line: 
    line = esva_w4.readline()
    values = line.split()
    if not values:
        break
    A_w4.append(float(values[0]))
    E_w4.append(float(values[1]))
esva_w4.close()

line = True
while line:   
    line = esvv_w4.readline()
    values = line.split()
    if not values:
        break
    A3_w4.append(float(values[0]))
esvv_w4.close()


#evfit data
line = True
while line:  
    line  = evfit_w4.readline()
    lines = line.split()
    if lines != []:
        if lines[0] =="A0=":
            A_min_w4 = float(lines[1])
            K0_w4    = float(lines[3])
            DK0_w4   = float(lines[3])
            D2K0_w4  = float(lines[3])
        if lines[0] =="Emin=":
            Emin_w4  = float(lines[1])
        if (lines[0] != "A0=") and (lines[0] != "Equation") and (lines[0] != "Emin="):
            Ev_eos_w4.append(float(lines[2]))
            Err_w4.append(float(lines[3]))
evfit_w4.close()





# Converting everything to numpy arrays for manpulation
latt_param_n1 = np.array(A_n1)
latt_param_n4 = np.array(A_n4) 
latt_param_m1 = np.array(A_m1) 
latt_param_m4 = np.array(A_m4) 
latt_param_w1 = np.array(A_w1) 
latt_param_w4 = np.array(A_w4) 

a3_n1 = np.array(A3_n1)
a3_n4 = np.array(A3_n4)
a3_m1 = np.array(A3_m1)
a3_m4 = np.array(A3_m4)
a3_w1 = np.array(A3_w1)
a3_w4 = np.array(A3_w4)

dft_eng_n1 = np.array(E_n1)
dft_eng_n4 = np.array(E_n4)
dft_eng_m1 = np.array(E_m1)
dft_eng_m4 = np.array(E_m4)
dft_eng_w1 = np.array(E_w1)
dft_eng_w4 = np.array(E_w4)

#shifting the energy values
dft_eng_n1 = dft_eng_n1 + energy_offset 
dft_eng_n4 = dft_eng_n4 + energy_offset
dft_eng_m1 = dft_eng_m1 + energy_offset
dft_eng_m4 = dft_eng_m4 + energy_offset
dft_eng_w1 = dft_eng_w1 + energy_offset
dft_eng_w4 = dft_eng_w4 + energy_offset



ev_eos_n1 = np.array(Ev_eos_n1) 
ev_eos_n4 = np.array(Ev_eos_n4)
ev_eos_m1 = np.array(Ev_eos_m1)
ev_eos_m4 = np.array(Ev_eos_m4)
ev_eos_w1 = np.array(Ev_eos_w1)
ev_eos_w4 = np.array(Ev_eos_w4)


ev_eos_n1 = ev_eos_n1 + energy_offset
ev_eos_n4 = ev_eos_n4 + energy_offset
ev_eos_m1 = ev_eos_m1 + energy_offset
ev_eos_m4 = ev_eos_m4 + energy_offset
ev_eos_w1 = ev_eos_w1 + energy_offset
ev_eos_w4 = ev_eos_w4 + energy_offset


err_n1 = np.array(Err_n1)
err_n4 = np.array(Err_n4)
err_m1 = np.array(Err_m1)
err_m4 = np.array(Err_m4)
err_w1 = np.array(Err_w1)
err_w4 = np.array(Err_w4)



#consolidating to a single plot
god_fig, axarr = plt.subplots(3,3, figsize = (fig_size_x,fig_size_y), sharex='col')


#Settign an overall title
god_fig.suptitle("Analysis of Birch 1st Order and Murnaghan Equation of State", fontsize=title_font)



#Top Row.  DFT Values vs EOS Values.  Narrow, Medium, Wide band. EOS = Birch 1st order.
axarr[0, 0].plot(latt_param_n1, dft_eng_n1, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[0, 0].plot(latt_param_n1, ev_eos_n1, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[0, 0].set_ylabel(' Birch 1st Order EOS \n Energy ($eV$)' , fontsize = y_label)
axarr[0, 0].set_title('\n Narrow Band', fontsize = subtitle_font)
axarr[0, 0].grid(True, linestyle='-.')
axarr[0, 0].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[0, 0].axhline(linewidth=2, color='k')

axarr[0, 1].plot(latt_param_m1, dft_eng_m1, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[0, 1].plot(latt_param_m1, ev_eos_m1, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[0, 1].set_title('\n Medium Band', fontsize = subtitle_font)
axarr[0, 1].grid(True, linestyle='-.')
axarr[0, 1].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[0, 1].axhline(linewidth=2, color='k')

axarr[0, 2].plot(latt_param_w1, dft_eng_w1, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[0, 2].plot(latt_param_w1, ev_eos_w1, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[0, 2].set_title('\n Wide Band', fontsize = subtitle_font)
axarr[0, 2].grid(True, linestyle='-.')
axarr[0, 2].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[0, 2].axhline(linewidth=2, color='k')

#Middle Row.  DFT Values vs EOS Values.  Narrow, Medium, Wide band, EOS = Murnaghan. 
axarr[1, 0].plot(latt_param_n4, dft_eng_n4, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[1, 0].plot(latt_param_n4, ev_eos_n4, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[1, 0].set_ylabel(' Mirnaghan EOS \n  Energy ($eV$) ' , fontsize = y_label)
axarr[1, 0].grid(True, linestyle='-.')
axarr[1, 0].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[1, 0].axhline(linewidth=2, color='k')

axarr[1, 1].plot(latt_param_m4, dft_eng_m4, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[1, 1].plot(latt_param_m4, ev_eos_m4, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[1, 1].grid(True, linestyle='-.')
axarr[1, 1].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[1, 1].axhline(linewidth=2, color='k')

axarr[1, 2].plot(latt_param_w4, dft_eng_w4, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[1, 2].plot(latt_param_w4, ev_eos_w4, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[1, 2].grid(True, linestyle='-.')
axarr[1, 2].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[1, 2].axhline(linewidth=2, color='k')


#Bottom Row.  Difference in energy values between the two EOSs.  Narrow, Medium, Wide band
axarr[2, 0].plot(latt_param_n1, err_n1, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[2, 0].plot(latt_param_n4, err_n4, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
axarr[2, 0].set_xlabel(' Lattice Parameter ($\AA$)' , fontsize = x_label)
axarr[2, 0].set_ylabel(' Difference between \n EOS and DFT ($eV$)' , fontsize = y_label)
axarr[2, 0].set_ylim(-1,1)
axarr[2, 0].grid(True, linestyle='-.')
axarr[2, 0].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[2, 0].axhline(linewidth=2, color='k')

axarr[2, 1].plot(latt_param_m1, err_m1, marker1, fillstyle= 'none', markersize = mark_size, color = color1)
axarr[2, 1].plot(latt_param_m4, err_m4, marker2, fillstyle = 'none', markersize = mark_size, color = color2)
#axarr[2, 1].set_xlabel(' Lattice Parameter ($\AA$)' , fontsize = x_label)
axarr[2, 1].set_ylim(-1,1)
axarr[2, 1].grid(True, linestyle='-.')
axarr[2, 1].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[2, 1].axhline(linewidth=2, color='k')

axarr[2, 2].plot(latt_param_w1, err_w1, marker1, fillstyle= 'none', markersize = mark_size, color = color1, label="Birch 1st Order")
axarr[2, 2].plot(latt_param_w4, err_w4, marker2, fillstyle = 'none', markersize = mark_size, color = color2, label="Murnaghan")
axarr[2, 2].set_xlabel(' Lattice Parameter ($\AA$)' , fontsize = x_label)
axarr[2, 2].set_ylim(-1,1)
axarr[2, 2].grid(True, linestyle='-.')
axarr[2, 2].tick_params(labelcolor='k', labelsize=tick_font, width=3)
axarr[2, 2].axhline(linewidth=2, color='k')


#Creating a legend for the overall plot
handles, labels = axarr[2, 2].get_legend_handles_labels()
god_fig.legend(handles, labels, loc=leg_loc, prop={'size': legend_size})
#bbox_to_anchor=(0.5, -0.01),

# # Printing plots to a .pdf file
god_fig.savefig("EOS_analysis.pdf", bbox_inches = 'tight')
god_fig.savefig("EOS_analysis.png", bbox_inches = 'tight')


plt.show()





</pre>
|}


[[Code: Quantum Espresso]]
[[Category:Repository of Codes]]
[[Category: Python]]