1414plt .rcParams .update ({"font.size" : font_size })
1515
1616def flip (items , ncol ):
17- return list (itertools .chain (* [items [i ::ncol ] for i in range (ncol )]))
17+ return list (itertools .chain (* [items [i ::ncol ] for i in range (ncol )]))
1818
1919def inner_mod (a ,b ):
20- res = a % b
21- return res if not res else res - b if a < 0 else res
20+ res = a % b
21+ return res if not res else res - b if a < 0 else res
2222
2323colors = {
24- "steel_blue" : "#1F77B4" ,
25- "light_steel_blue" : "#AEC7E8" ,
26- "orange" : "#FF7F0E" ,
27- "light_orange" : "#FFBB78" ,
28- "forest_green" : "#2CA02C" ,
29- "light_green" : "#98DF8A" ,
30- "firebrick_red" : "#D62728" ,
31- "soft_red" : "#FF9896" ,
32- "lavender" : "#9467BD" ,
33- "light_lavender" : "#C5B0D5" ,
34- "brown" : "#8C564B" ,
35- "tan" : "#C49C94" ,
36- "orchid" : "#E377C2" ,
37- "light_orchid" : "#F7B6D2" ,
38- "gray" : "#7F7F7F" ,
39- "light_gray" : "#C7C7C7" ,
40- "yellow_green" : "#BCBD22" ,
41- "light_yellow_green" : "#DBDB8D" ,
42- "turquoise" : "#17BECF" ,
43- "light_turquoise" : "#9EDAE5"
24+ "steel_blue" : "#1F77B4" ,
25+ "light_steel_blue" : "#AEC7E8" ,
26+ "orange" : "#FF7F0E" ,
27+ "light_orange" : "#FFBB78" ,
28+ "forest_green" : "#2CA02C" ,
29+ "light_green" : "#98DF8A" ,
30+ "firebrick_red" : "#D62728" ,
31+ "soft_red" : "#FF9896" ,
32+ "lavender" : "#9467BD" ,
33+ "light_lavender" : "#C5B0D5" ,
34+ "brown" : "#8C564B" ,
35+ "tan" : "#C49C94" ,
36+ "orchid" : "#E377C2" ,
37+ "light_orchid" : "#F7B6D2" ,
38+ "gray" : "#7F7F7F" ,
39+ "light_gray" : "#C7C7C7" ,
40+ "yellow_green" : "#BCBD22" ,
41+ "light_yellow_green" : "#DBDB8D" ,
42+ "turquoise" : "#17BECF" ,
43+ "light_turquoise" : "#9EDAE5"
4444}
4545
4646custom_cmap = ListedColormap (colors .values ())
@@ -63,12 +63,12 @@ def inner_mod(a,b):
6363
6464elements = "" # Initialize an empty string to store species data
6565with open ("brawl.inp" , 'r' ) as file :
66- for line in file :
67- if 'species_name' in line : # Check if the line contains the species_name variable
68- parts = line .split ("=" ) # Split the line by '=' and get everything after it
69- if len (parts ) > 1 :
70- elements = parts [1 ].strip ().split () # Get the data after '=' and strip spaces
71- break # Exit the loop after finding species_name
66+ for line in file :
67+ if 'species_name' in line : # Check if the line contains the species_name variable
68+ parts = line .split ("=" ) # Split the line by '=' and get everything after it
69+ if len (parts ) > 1 :
70+ elements = parts [1 ].strip ().split () # Get the data after '=' and strip spaces
71+ break # Exit the loop after finding species_name
7272n_species = len (elements )
7373concentrations = 1.0 / n_species * np .ones (n_species )
7474
@@ -103,19 +103,19 @@ def inner_mod(a,b):
103103
104104wcs = np .zeros ((len (U_data ),n_species ,n_species ,2 ))
105105for i in range (len (U_data )):
106- wcoft = np .zeros ((n_species ,n_species ,2 ))
107- for j in range (len (elements )):
108- for k in range (len (elements )):
109- wcoft [j ,k ,0 ] = 1.0 - 1.0 / 8.0 * rho [i ,1 ,j ,k ]/ concentrations [j ]
110- wcoft [j ,k ,1 ] = 1.0 - 1.0 / 6.0 * rho [i ,2 ,j ,k ]/ concentrations [j ]
111- wcs [i ] = wcoft
106+ wcoft = np .zeros ((n_species ,n_species ,2 ))
107+ for j in range (len (elements )):
108+ for k in range (len (elements )):
109+ wcoft [j ,k ,0 ] = 1.0 - 1.0 / 8.0 * rho [i ,1 ,j ,k ]/ concentrations [j ]
110+ wcoft [j ,k ,1 ] = 1.0 - 1.0 / 6.0 * rho [i ,2 ,j ,k ]/ concentrations [j ]
111+ wcs [i ] = wcoft
112112
113113pairs = np .zeros ([int (len (elements )* (len (elements )+ 1 )/ 2 ),2 ], dtype = np .int16 )
114114k = 0
115115for i in range (len (elements )):
116- for j in range (i , len (elements )):
117- pairs [k ] = [i , j ]
118- k += 1
116+ for j in range (i , len (elements )):
117+ pairs [k ] = [i , j ]
118+ k += 1
119119
120120fig_height = 4 ; fig_width = fig_height * 1.68
121121
@@ -133,8 +133,8 @@ def inner_mod(a,b):
133133prob = np .zeros (len (bin_edges )- 1 )
134134# Reweight histogram
135135for ibin , edge in enumerate (bin_edges [:- 1 ]):
136- bin_energy = edge + 0.5 * bin_width
137- prob [ibin ] = wl_logdos [ibin ]- beta * bin_energy
136+ bin_energy = edge + 0.5 * bin_width
137+ prob [ibin ] = wl_logdos [ibin ]- beta * bin_energy
138138prob = prob - np .max (prob )/ 2
139139prob = np .exp (prob )
140140# Normalise
@@ -147,68 +147,68 @@ def inner_mod(a,b):
147147# Loop over temperatures of interest
148148for itemp , new_temp in enumerate (temperatures ):
149149
150- beta = 1.0 / (kb_ev * new_temp )
151-
152- # Reweighted histogram
153- prob = np .zeros (len (bin_edges )- 1 )
154-
155- # Reweight histogram
156- for ibin , edge in enumerate (bin_edges [:- 1 ]):
157- bin_energy = edge + 0.5 * bin_width
158- prob [ibin ] = wl_logdos [ibin ]- beta * bin_energy
159-
160- prob = prob - np .max (prob )
161- prob = np .exp (prob )
162-
163- # Normalise
164- prob = prob / (np .sum (prob ))
165- prob = np .nan_to_num (prob )
166-
167- # Only plot every 5th histogram to avoid crowding the axes
168- if np .isin (temperatures_plot , int (new_temp )).any () == True :
169- strlabel = "T={} K" .format (int (new_temp ))
170- index_to_zero = np .where (prob < 1e-8 )
171- hist_prob = copy .deepcopy (prob )
172- hist_prob [index_to_zero ] = 0
173- non_zero = np .nonzero (hist_prob )
174- hist_ax .stairs (prob , bin_edges / n_atoms * ev_to_mev - zero_energy , fill = True , alpha = 0.1 )
175- hist_ax .stairs (prob , bin_edges / n_atoms * ev_to_mev - zero_energy , label = strlabel , fill = False , alpha = 1 )
176- if (bin_edges [np .max (non_zero )+ 1 ]/ n_atoms * ev_to_mev - zero_energy > hist_max ):
177- hist_max = bin_edges [np .max (non_zero )+ 1 ]/ n_atoms * ev_to_mev - zero_energy
178- if (bin_edges [np .min (non_zero )]/ n_atoms * ev_to_mev - zero_energy < hist_min ):
179- hist_min = bin_edges [np .min (non_zero )]/ n_atoms * ev_to_mev - zero_energy
180- if (np .max (prob ) < 0.4 and np .max (prob ) > prob_max ):
181- prob_max = np .max (prob )
182-
183- # Mean energy
184- mean_energy = np .dot (bin_edges [:- 1 ]+ 0.5 * bin_width , prob )
185- mean_energies [itemp ] = mean_energy
186-
187- for ipair , pair in enumerate (pairs ):
188- asr_order_1 = np .dot (wcs [:,pair [0 ],pair [1 ],0 ], prob )
189- asr_orders_1 [ipair ,itemp ] = asr_order_1
190-
191- # Compute heat capacity using the histogram
192- msq_dev = np .zeros (len (bin_edges )- 1 )
193- for ibin , edge in enumerate (bin_edges [:- 1 ]):
194- bin_energy = edge + 0.5 * bin_width
195- msq_dev [ibin ] = (bin_energy - mean_energies [itemp ])** 2
196-
197- heat_caps [itemp ] = np .dot (msq_dev , prob )* bin_width / (kb_ev * new_temp ** 2 )
150+ beta = 1.0 / (kb_ev * new_temp )
151+
152+ # Reweighted histogram
153+ prob = np .zeros (len (bin_edges )- 1 )
154+
155+ # Reweight histogram
156+ for ibin , edge in enumerate (bin_edges [:- 1 ]):
157+ bin_energy = edge + 0.5 * bin_width
158+ prob [ibin ] = wl_logdos [ibin ]- beta * bin_energy
159+
160+ prob = prob - np .max (prob )
161+ prob = np .exp (prob )
162+
163+ # Normalise
164+ prob = prob / (np .sum (prob ))
165+ prob = np .nan_to_num (prob )
166+
167+ # Only plot every 5th histogram to avoid crowding the axes
168+ if np .isin (temperatures_plot , int (new_temp )).any () == True :
169+ strlabel = "T={} K" .format (int (new_temp ))
170+ index_to_zero = np .where (prob < 1e-8 )
171+ hist_prob = copy .deepcopy (prob )
172+ hist_prob [index_to_zero ] = 0
173+ non_zero = np .nonzero (hist_prob )
174+ hist_ax .stairs (prob , bin_edges / n_atoms * ev_to_mev - zero_energy , fill = True , alpha = 0.1 )
175+ hist_ax .stairs (prob , bin_edges / n_atoms * ev_to_mev - zero_energy , label = strlabel , fill = False , alpha = 1 )
176+ if (bin_edges [np .max (non_zero )+ 1 ]/ n_atoms * ev_to_mev - zero_energy > hist_max ):
177+ hist_max = bin_edges [np .max (non_zero )+ 1 ]/ n_atoms * ev_to_mev - zero_energy
178+ if (bin_edges [np .min (non_zero )]/ n_atoms * ev_to_mev - zero_energy < hist_min ):
179+ hist_min = bin_edges [np .min (non_zero )]/ n_atoms * ev_to_mev - zero_energy
180+ if (np .max (prob ) < 0.4 and np .max (prob ) > prob_max ):
181+ prob_max = np .max (prob )
182+
183+ # Mean energy
184+ mean_energy = np .dot (bin_edges [:- 1 ]+ 0.5 * bin_width , prob )
185+ mean_energies [itemp ] = mean_energy
186+
187+ for ipair , pair in enumerate (pairs ):
188+ asr_order_1 = np .dot (wcs [:,pair [0 ],pair [1 ],0 ], prob )
189+ asr_orders_1 [ipair ,itemp ] = asr_order_1
190+
191+ # Compute heat capacity using the histogram
192+ msq_dev = np .zeros (len (bin_edges )- 1 )
193+ for ibin , edge in enumerate (bin_edges [:- 1 ]):
194+ bin_energy = edge + 0.5 * bin_width
195+ msq_dev [ibin ] = (bin_energy - mean_energies [itemp ])** 2
196+
197+ heat_caps [itemp ] = np .dot (msq_dev , prob )* bin_width / (kb_ev * new_temp ** 2 )
198198
199199heat_caps = np .gradient (mean_energies , temperatures )
200200
201201gibbs_energies [0 ] = mean_energies [0 ]
202202gibbs_energies [1 ] = mean_energies [0 ]
203203for itemp in range (2 ,len (temperatures )):
204204
205- beta_i = 1.0 / (kb_ev * temperatures [itemp ])
206- beta_j = 1.0 / (kb_ev * temperatures [itemp - 1 ])
205+ beta_i = 1.0 / (kb_ev * temperatures [itemp ])
206+ beta_j = 1.0 / (kb_ev * temperatures [itemp - 1 ])
207207
208- gibbs_energies [itemp ] = beta_j * gibbs_energies [itemp - 1 ]/ beta_i + ((mean_energies [itemp - 1 ]+ mean_energies [itemp ])* (beta_i - beta_j ))/ (2 * beta_i )
208+ gibbs_energies [itemp ] = beta_j * gibbs_energies [itemp - 1 ]/ beta_i + ((mean_energies [itemp - 1 ]+ mean_energies [itemp ])* (beta_i - beta_j ))/ (2 * beta_i )
209209
210210for itemp , new_temp in enumerate (temperatures ):
211- entropies [itemp ] = (mean_energies [itemp ]- gibbs_energies [itemp ])/ new_temp
211+ entropies [itemp ] = (mean_energies [itemp ]- gibbs_energies [itemp ])/ new_temp
212212
213213
214214# Complete plots using data computed above
@@ -246,11 +246,11 @@ def inner_mod(a,b):
246246asro_min = 0
247247cv_ax1 .plot (temperatures , heat_caps , '-o' , markersize = 4 , label = "Heat Capacity" )
248248for ipair , pair in enumerate (pairs ):
249- cv_ax1_asro .plot (temperatures , asr_orders_1 [ipair ], label = elements [pair [0 ]] + '-' + elements [pair [1 ]])
250- if (np .max (asr_orders_1 [ipair ]) > asro_max ):
251- asro_max = np .max (asr_orders_1 [ipair ])
252- if (np .min (asr_orders_1 [ipair ]) < asro_min ):
253- asro_min = np .min (asr_orders_1 [ipair ])
249+ cv_ax1_asro .plot (temperatures , asr_orders_1 [ipair ], label = elements [pair [0 ]] + '-' + elements [pair [1 ]])
250+ if (np .max (asr_orders_1 [ipair ]) > asro_max ):
251+ asro_max = np .max (asr_orders_1 [ipair ])
252+ if (np .min (asr_orders_1 [ipair ]) < asro_min ):
253+ asro_min = np .min (asr_orders_1 [ipair ])
254254
255255cv_ticks = 400
256256x_ticks_subplots = np .arange (np .around (start_temp / cv_ticks , decimals = 0 )* cv_ticks , np .around (end_temp / cv_ticks , decimals = 0 )* cv_ticks + cv_ticks , cv_ticks )
@@ -267,33 +267,33 @@ def inner_mod(a,b):
267267
268268# Loop through the legend handles and change their linewidth
269269for handle in hist_ax_legend .legend_handles :
270- handle .set_linewidth (3 ) # Increase the line width for each legend line
270+ handle .set_linewidth (3 ) # Increase the line width for each legend line
271271
272272alternate_colors = list (colors .values ())[::2 ]
273273cv_ax1_asro .set_prop_cycle (cycler (color = alternate_colors ))
274274
275275margin_pct = 0.05
276276for ax in [hist_ax , cv_ax1 , cv_ax1_asro ]:
277- # Get axis limits
278- x_min , x_max = ax .get_xlim ()
279- y_min , y_max = ax .get_ylim ()
277+ # Get axis limits
278+ x_min , x_max = ax .get_xlim ()
279+ y_min , y_max = ax .get_ylim ()
280280
281- # Calculate the margin values (5% of the range)
282- x_margin = (x_max - x_min ) * margin_pct
283- y_margin = (y_max - y_min ) * margin_pct
281+ # Calculate the margin values (5% of the range)
282+ x_margin = (x_max - x_min ) * margin_pct
283+ y_margin = (y_max - y_min ) * margin_pct
284284
285- # Get the current tick positions
286- x_ticks = ax .get_xticks ()
287- y_ticks = ax .get_yticks ()
285+ # Get the current tick positions
286+ x_ticks = ax .get_xticks ()
287+ y_ticks = ax .get_yticks ()
288288
289- # Filter out ticks within the margin
290- x_ticks_filtered = [tick for tick in x_ticks if x_min + x_margin <= tick <= x_max - x_margin or tick == 0 ]
289+ # Filter out ticks within the margin
290+ x_ticks_filtered = [tick for tick in x_ticks if x_min + x_margin <= tick <= x_max - x_margin or tick == 0 ]
291291 y_ticks_filtered = [tick for tick in y_ticks if y_min + y_margin <= tick <= y_max - y_margin or tick == 0 ]
292292
293293
294- # Set the new ticks
295- #ax.set_xticks(x_ticks_filtered)
296- #ax.set_yticks(y_ticks_filtered)
294+ # Set the new ticks
295+ #ax.set_xticks(x_ticks_filtered)
296+ #ax.set_yticks(y_ticks_filtered)
297297
298298asro_diff = np .abs (asro_min - asro_max )
299299hist_ax .set_xlim (hist_min ,hist_max )
0 commit comments