python_scripts/thermal_expansion_output_comparison.py at main · PascalSalzbrenner/python_scripts · GitHub

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
# short script to generate files suitable for comparing the lattice parameters/volumes/free energies of different structures as a function of temperature
# input: directories containing minimum_energy.txt files, as generated by postprocess_geometry_optimisation_cubic.sh

# written by Pascal Salzbrenner - pts28@cam.ac.uk

import sys
import math
import numpy as np

T_max = int(input("What is the maximum temperature [K] up to which you would like to plot the structures? "))
quantity = input("Which quantity would you like to plot? (lattice_parameter, volume, energy) ").lower()

if quantity.startswith("l"):
    # lattice parameter
    index = 1
    subtraction_string = ""
    y_label = "Lattice Parameter [A]"
elif quantity.startswith("v"):
    # volume
    index = 2
    subtraction_string = ""
    y_label = "Volume [A^3]"
else:
    # energy - default

    # add the actual number once the value has been found - the minimum energy at T_max
    index = 3
    subtraction_string="-"
    energy_list = []
    y_label = "Free Energy [meV/atom]"

for directory in sys.argv[1:]:
    # iterate over and read quantities

    infile = open("{}/minimum_energy.txt".format(directory.rstrip("/")), "r")
    outfile = open("{}/{}_temperature.plot".format(directory.rstrip("/"), quantity), "w")

    header_line = infile.readline().split()
    outfile.write("#T [K], {}\n".format(header_line[index].rstrip(",")))

    # read past next line - static lattice
    infile.readline()

    for line in infile:
        split_line = line.split()
        outfile.write("{} {}\n".format(split_line[0].split("_")[1].rstrip("K"), split_line[index]))

        if int(split_line[0].split("_")[1].rstrip("K")) >= T_max:
            if index == 3:
                energy_list.append(float(split_line[3]))

            infile.close()
            outfile.close()

            break

if index == 3:
    subtraction_string += "({})".format(min(energy_list))

# add writing of gnuplot file to plot fit
with open("{}_temperature.gnu".format(quantity), "w") as plotfile:

    plotfile.write("set terminal postscript eps colour font 'Helvetica,20'\n")
    plotfile.write("set style data points\n")
    plotfile.write("set key box lt -1 lw 2 width 2 height 1.5 opaque font 'Helvetica,15'\n")
    plotfile.write("set output '| epstopdf --filter --outfile={}_temperature.pdf'\n".format(quantity))

    # redefine gnuplot linetypes with nice colours
    plotfile.write("set linetype 1 lc rgb '#DC143C'\n")
    plotfile.write("set linetype 2 lc rgb '#66A61E'\n")
    plotfile.write("set linetype 3 lc rgb '#E6AB02'\n")
    plotfile.write("set linetype 4 lc rgb '#D95F02'\n")
    plotfile.write("set linetype 5 lc rgb '#191979'\n")
    plotfile.write("set linetype 6 lc rgb '#7570B3'\n")
    plotfile.write("set linetype 7 lc rgb '#E7298A'\n")
    plotfile.write("set linetype 8 lc rgb '#1E90FF'\n")
    plotfile.write("set linetype 9 lc rgb '#1B9E77'\n")
    plotfile.write("set linetype 10 lc rgb '#B8860B'\n")
    plotfile.write("set linetype cycle 10\n")

    plotfile.write("set xlabel 'T [K]'\n")
    plotfile.write("set ylabel '{}'\n".format(y_label))

    plot_string = "plot"

    for directory in sys.argv[1:]:
        plot_string += " '{0}/{1}_temperature.plot' u 1:($2{2}) w points pt 7 ps 0.5 title '{0}',".format(directory.rstrip("/"), quantity, subtraction_string)

    plot_string.rstrip(",")

    plotfile.write(plot_string)