dft_eddp_dev_enthalpies.py

# script to plot the DFT and EDDP enthalpies contained in DFT/EDDP_enthalpy.agr files (generated by cryan and renamed)
# need to ensure both DFT and EDDP data are referenced to the same structure
# EDDP deviations are added as line broadenings to the EDDP enthalpy curves

import os
import sys
import numpy as np
import matplotlib.pyplot as plt

from copy import deepcopy
from scipy import interpolate

# note that the default labels in the _enthalpy.agr files are chemical_formula-structure_name
# for consistency with the dev_press_vol files, we use only the structure_name part wherever the name of a structure is concerned

# we convert to int, there is no good reason really to start at fractional pressures, and int makes things cleaner
low_press = int(sys.argv[1])
high_press = int(sys.argv[2])
reference_structure = sys.argv[3]

press_range = high_press - low_press

# check if user-supplied labels are present
# labels.txt format:
# structure_name_in_file label_for plot
# ...

labels = {}
if "labels.txt" in os.listdir():
	with open("labels.txt", "r") as labels_file:
		for line in labels_file:
			label = line.split()
			labels[label[0]] = label[1]
else:
	# for unified access, we fill the labels dictionary anyways, even though the label is just identical to the name
	labels[reference_structure] = reference_structure

# read input data
dft_data = {reference_structure: [np.linspace(low_press, high_press, high_press-low_press+1), np.zeros(high_press-low_press+1)]}
eddp_data = {reference_structure: [np.linspace(low_press, high_press, high_press-low_press+1), np.zeros(high_press-low_press+1)]}
structure_list = [reference_structure]

for theory in ["DFT", "EDDP"]:

	infile = open("{}_enthalpy.agr".format(theory), "r")

	for line in infile:

		# the data we want is stored in the blocks at the end of the file; the beginning of these blocks is indicated by a few header lines setting the style. This is what we detect.
		if "legend" in line and "s" in line:
			# extra check for "s" to avoid detecting the lines including "legend" in the general file preamble setting global style
			
			if reference_structure in line:
				# we have already set the data for the reference structure to 0; ignore it
				continue
			else:
				# this we read
				# extract the structure name
				# here, we must remove the chemical_formula before the first hyphen to make the naming consistent
				structure_name_parts = line.split()[-1].strip('"').split("-")
				structure_name = "-".join(structure_name_parts[1:])

				if not structure_name in structure_list:
					structure_list.append(structure_name)

				if not structure_name in labels.keys():
					labels[structure_name] = structure_name

				# skip the next 6 lines which are headers
				for i in range(6):
					infile.readline()

				# now we get to the data
				pressure = []
				enthalpy = []

				for data_line in infile:

					if data_line.startswith("&"):
						# "&" marks the end of a data block
						if theory == "DFT":
							dft_data[structure_name] = [np.array(deepcopy(pressure)), np.array(deepcopy(enthalpy))]
						else:
							eddp_data[structure_name] = [np.array(deepcopy(pressure)), np.array(deepcopy(enthalpy))]
						break

					else:
						data = data_line.split()

						current_pressure = float(data[0])

						if (current_pressure > low_press or np.isclose(current_pressure, low_press)) and (current_pressure < high_press or np.isclose(current_pressure, high_press)):

							pressure.append(current_pressure)
							enthalpy.append(1000*float(data[1])) # convert from eV/atom to meV/atom

		else:
			continue

	infile.close()

# read dev as a function of pressure
deviations = {}
for structure, pressure_enthalpy in eddp_data.items():
	dev_pressures = []
	
	# we add and subtract the deviations to get the width of the curve
	deviation_list = []
	with open("dev_press_vol_{}.txt".format(structure), "r") as dev_file:
		dev_file.readline()
		for line in dev_file:
			data = line.split()
			dev_pressures.append(float(data[0]))
			deviation_list.append(float(data[2]))

	# interpolate the deviations to ensure we have can access them at the pressures where we have the enthalpy data
	# then add and subtract it from the energy itself to get the upper and lower bounds for the curve width

	deviation = interpolate.splrep(dev_pressures, deviation_list, s=0)
	deviations_p = pressure_enthalpy[1]+interpolate.splev(pressure_enthalpy[0], deviation, der=0)
	deviations_n = pressure_enthalpy[1]-interpolate.splev(pressure_enthalpy[0], deviation, der=0)
	deviations[structure] = [pressure_enthalpy[0], deviations_p, deviations_n]

# plotting
# define colour list
colours = ["#E6AB02", "#66A61E", "#8000C4", "#7570B3", "#E7298A", "#1E90FF", "#1B9E77", "#20C2C2", "#D95F02", "#DC143C"]

plt.xlabel("Pressure [GPa]")
plt.ylabel("Enthalpy [meV/formula unit]")
plt.xlim(low_press, high_press)

# determine sensible tick spacing on the x-axis
max = round(press_range/50)*5

for i in range(max, 0, -1):

	if press_range%i == 0:
		tick_step = i
		break

plt.xticks(np.arange(low_press, high_press+1, tick_step))

# to set a sensible enthalpy window, determine the largest and lowest enthalpies
max_enthalpy = 0
min_enthalpy = 0

# plot EDDP data
for structure, pressure_enthalpy in eddp_data.items():

	if structure == reference_structure:
		struc_label = labels[structure]
	else:
		struc_label = "{} - EDDP".format(labels[structure])

	plt.plot(pressure_enthalpy[0], pressure_enthalpy[1], color=colours[structure_list.index(structure)], linestyle="solid", label=struc_label)
	plt.fill_between(pressure_enthalpy[0], deviations[structure][1], deviations[structure][2], color=colours[structure_list.index(structure)], alpha=0.5)

	# to determine the maximum (minimum) value, we need only consider the energies with positive (negative) deviation; this is guaranteed to have be higher (lower) than the actual enthalpy
	if np.amax(deviations[structure][1]) > max_enthalpy:
		max_enthalpy = np.amax(deviations[structure][1])
	if np.amin(deviations[structure][2]) < min_enthalpy:
		min_enthalpy = np.amin(deviations[structure][2])

# plot DFT data
for structure, pressure_enthalpy in dft_data.items():

	if structure == reference_structure:
		continue

	plt.plot(pressure_enthalpy[0], pressure_enthalpy[1], color=colours[structure_list.index(structure)], linestyle="dashed", label="{} - DFT".format(labels[structure]))

	if np.amax(pressure_enthalpy[1]) > max_enthalpy:
		max_enthalpy = np.amax(pressure_enthalpy[1])
	if np.amin(pressure_enthalpy[1]) < min_enthalpy:
		min_enthalpy = np.amin(pressure_enthalpy[1])

# round maximum and minimum enthalpies to the next 10 to get y range
y_max = round(max_enthalpy, -1)
y_min = round(min_enthalpy, -1)

# make sure the entire range is included
if y_max < max_enthalpy:
	y_max += 10
if y_min > min_enthalpy:
	y_min -= 10

plt.ylim(y_min, y_max)

plt.legend()

plt.savefig("dft_eddp_dev_enthalpies.png", dpi=300)
#plt.savefig("dft_eddp_dev_enthalpies.pdf")
plt.close()