-
Notifications
You must be signed in to change notification settings - Fork 0
/
InTerSect_TS_level.py
109 lines (77 loc) · 3.55 KB
/
InTerSect_TS_level.py
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
import sys
from sympy import *
import sympy as sym
import os
from itertools import chain
import pickle as pl
# Intial candidates for fit, per FU: - thus, the E vs V input data has to be per FU
E0_init = -941.510817926696
V0_init = 63.54960592453
B0_init = 76.3746233515232
B0_prime_init = 4.05340727164527
def BM(x, a, b, c, d):
return a + b*x + c*x**2 + d*x**3
def P(x, b, c, d):
return 4.3597482E+3 * (-b - 2*c*x - 3 *d*x**2)
def H(x, a, b, c, d):
return a + b*x + c*x**2 + d*x**3
filefolder_Calcite_I_SG_167 = '../../Files_Outputs/Calcite_I/TS_vs_V'
filefolder_Calcite_II_SG_14 = '../../Files_Outputs/Calcite_II/TS_vs_V'
filefolder_energetics = 'EL_vs_V'
# Calcite I (Red triangles):
V_C_I, E_C_I = np.loadtxt(os.path.join(filefolder_Calcite_I_SG_167, filefolder_energetics, 'EL_vs_V.dat'), skiprows = 1).T
# 14 (Empty grey triangles):
V_14, E_14 = np.loadtxt(os.path.join(filefolder_Calcite_II_SG_14, filefolder_energetics, 'EL_vs_V.dat'), skiprows = 1).T
init_vals = [E0_init, V0_init, B0_init, B0_prime_init]
popt_C_I, pcov_C_I = curve_fit(BM, V_C_I, E_C_I, p0=init_vals)
popt_14, pcov_14 = curve_fit(BM, V_14, E_14, p0=init_vals)
# Linspace for plotting the fitting curves:
V_C_I_lin = np.linspace(V_C_I[0], V_C_I[-1], 10000)
V_14_lin = np.linspace(V_14[0], V_14[-1], 10000)
fig_handle = plt.figure()
# Plotting the fitting curves:
p2, = plt.plot(V_C_I_lin, BM(V_C_I_lin, *popt_C_I), color='black', label='Cubic fit Calcite I' )
p6, = plt.plot(V_14_lin, BM(V_14_lin, *popt_14), 'b', label='Cubic fit Calcite II')
# Plotting the scattered points:
p1 = plt.scatter(V_C_I, E_C_I, color='red', marker="^", label='Calcite I', s=100)
p5 = plt.scatter(V_14, E_14, color='grey', marker="^", facecolors='none', label='Calcite II', s=100)
fontP = FontProperties()
fontP.set_size('15')
plt.legend((p1, p2, p5, p6), ("Calcite I", "Cubic fit Calcite I", "Calcite II", 'Cubic fit Calcite II'), prop=fontP)
global V0, B0, B0_prime
E0 = popt_C_I[0]
V0 = popt_C_I[1]
B0 = popt_C_I[2]
B0_prime = popt_C_I[3]
pressures_per_F_unit_C_I = P(V_C_I, V0, B0, B0_prime)
output_array_2 = np.vstack((E_C_I, V_C_I, pressures_per_F_unit_C_I)).T
np.savetxt('Volumes_and_pressures_C_I.dat', output_array_2, header="Energy / FU (a.u.) \t Volume / FU (A^3) \t Pressures (GPa)", fmt="%0.13f")
global V0_14, B0_14, B0_prime_14
E0_14 = popt_14[0]
V0_14 = popt_14[1]
B0_14 = popt_14[2]
B0_prime_14 = popt_14[3]
pressures_per_F_unit_14 = P(V_14, V0_14, B0_14, B0_prime_14)
output_array_2 = np.vstack((E_14, V_14, pressures_per_F_unit_14)).T
np.savetxt('Volumes_and_pressures_14.dat', output_array_2, header="Energy / FU (a.u.) \t Volume / FU (A^3) \t Pressures (GPa)", fmt="%0.13f")
# obtain the Temperature we are working at:
import subprocess
T = subprocess.check_output("basename `pwd`", shell=True)
print 'T = ', T
import string
new_T = string.replace(T, '_', ' ')
print 'new_T = ', new_T
T = T.rstrip()
new_T = new_T.rstrip()
plt.xlabel('$V$ / F.U. (Angstrom$^{3}$)', fontsize=20)
plt.ylabel(r'TS / F.U. (a.u.)', fontsize=20)
plt.suptitle("PBE-D3, pob-TZVP, SHRINK 8 8, Bipolar 18 18, TOLINTEG 8 18, XXLGRID, TOLDEE 8")
plt.title("(0.87 - 0.98)$V_{eq}$ and (0.98 - 1.08)$V_{eq}$. T = %s" %new_T, fontsize=10)
plt.ticklabel_format(useOffset=False)
plt.savefig('calcite_I_and_II_all_2_summary_better_plot.pdf', bbox_inches='tight')
pl.dump(fig_handle,file('sinus.pickle_calcite_I_and_II_all_2_summary_better_plot','w'))