EOS comp plot.py
From EVOCD
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
SUMMARY
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.
# 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() |