-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbm3_thermal_2.py
10074 lines (8433 loc) · 327 KB
/
bm3_thermal_2.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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Ab initio Elasticity and Thermodynamics of Minerals
#
# Version 2.8.1 18/10/2022
#
# Comment the following three lines to produce the documentation
# with readthedocs
# from IPython import get_ipython
# get_ipython().magic('cls')
# get_ipython().magic('reset -sf')
# Set to True to compile an executable version
exe_flag=False
import datetime
import os
import sys
import warnings
import scipy
import numpy as np
import matplotlib as mpl
if exe_flag:
mpl.use('Qt5Agg')
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
# from matplotlib import rc
import pandas as pd
import sympy as sym
import parame as pr
from scipy.optimize import curve_fit, fmin, minimize_scalar, minimize
from scipy.interpolate import UnivariateSpline, Rbf
from scipy import integrate
from plot import plot_class
from mineral_data import mineral, load_database, equilib, reaction,\
pressure_react, export, field, import_database, name_list
from mineral_data import ens, cor, py, coe, q, fo, ky, sill, andal, per, sp, \
mao, fmao, stv, cc, arag, jeff, jeff_fe, jeff_fe3p, jeff_feb, zrc, diam, \
tui, flo
warnings.filterwarnings('ignore')
import_database()
mpl.rcParams['figure.dpi']=80
class latex_class():
"""
Setup for the use of LaTeX for axis labels and titles; sets of parameters
for graphics output.
"""
def __init__(self):
self.flag=False
self.dpi=300
self.font_size=14
self.tick_size=12
self.ext='jpg'
mpl.rc('text', usetex=False)
def on(self):
self.flag=True
mpl.rc('text', usetex=True)
def off(self):
self.flag=False
mpl.rc('text', usetex=False)
def set_param(self, dpi=300, fsize=14, tsize=12, ext='jpg'):
"""
Args:
dpi: resolution of the graphics file (default 300)
fsize: size of the labels of the axes in points (default 14)
tsize: size of the ticks in points (default 12)
ext: extension of the graphics file (default 'jpg'); this argument
is only used in those routines where the name of the file is
automatically produced by the program (e.g. check_poly or
check_spline functions). In other cases, the extension is
directly part of the name of the file given as argument to
the function itself, and 'ext' is ignored.
"""
self.dpi=dpi
self.font_size=fsize
self.tick_size=tsize
self.ext=ext
def get_dpi(self):
return self.dpi
def get_fontsize(self):
return self.font_size
def get_ext(self):
return self.ext
def get_tsize(self):
return self.tick_size
class flag:
def __init__(self,value):
self.value=value
self.jwar=0
def on(self):
self.value=True
def off(self):
self.value=False
def inc(self):
self.jwar += 1
def reset(self):
self.jwar=0
class zero_point:
def __init__(self):
self.flag=True
def on(self):
self.flag=True
print("Zero point energy computation on")
def off(self):
self.flag=False
print("Zero point energy computation off")
class verbose_class():
def __init__(self,value):
self.flag=value
def on(self):
self.flag=True
print("Verbose mode on")
def off(self):
self.flag=False
print("Verbose mode off")
class BM3_error(Exception):
pass
class vol_corr_class:
def __init__(self):
self.flag=False
self.v0_init=None
def on(self):
self.flag=True
def off(self):
self.flag=False
def set_volume(self,vv):
self.v0_init=vv
class data_info():
"""
Stores information about the current settings
"""
def __init__(self):
self.min_static_vol=None
self.max_static_vol=None
self.static_points=None
self.min_freq_vol=None
self.max_freq_vol=None
self.freq_points=None
self.min_select_vol=None
self.max_select_vol=None
self.select_points=None
self.freq_sets=None
self.fit_type='No fit'
self.min_vol_fit=None
self.max_vol_fit=None
self.fit_points=None
self.fit_degree=None
self.fit_smooth=None
self.k0=None
self.kp=None
self.v0=None
self.temp=None
self.k0_static=None
self.kp_static=None
self.v0_static=None
self.popt=None
self.popt_orig=None
self.min_names=name_list.mineral_names
self.title=None
def show(self):
"""
Prints information about the current settings stored in the classes
"""
if self.title !=None:
print(self.title)
print("\nCurrent settings and results\n")
if self.min_static_vol != None:
print("Static data ** min, max volumes: %8.4f, %8.4f; points: %d"\
% (self.min_static_vol, self.max_static_vol, self.static_points))
if self.min_freq_vol != None:
print("Frequency volume range ** min, max volumes: %8.4f, %8.4f; points: %d"\
% (self.min_freq_vol, self.max_freq_vol, self.freq_points))
if self.min_select_vol != None:
print("Selected freq. sets ** min, max volumes: %8.4f, %8.4f; points: %d"\
% (self.min_select_vol, self.max_select_vol, self.select_points))
print("Frequency sets: %s" % str(self.freq_sets))
if self.fit_type != 'No fit':
if self.fit_type=='poly':
print("\nFit of frequencies ** type: %s, degree: %d" \
% (self.fit_type, self.fit_degree))
else:
print("\nFit of frequencies ** type: %s, degree: %d, smooth: %2.1f" \
% (self.fit_type, self.fit_degree, self.fit_smooth))
print(" min, max volumes: %8.4f, %8.4f; points %d" %\
(self.min_vol_fit, self.max_vol_fit, self.fit_points))
else:
print("No fit of frequencies")
if supercell.flag:
print("\n*** This is a computation performed on SUPERCELL data")
print(" (SCELPHONO and QHA keywords in CRYSTAL). Number of cells: %3i" % supercell.number)
if self.k0_static != None:
print("\n*** Static EoS (BM3) ***")
print("K0: %6.2f GPa, Kp: %4.2f, V0: %8.4f A^3" %\
(self.k0_static, self.kp_static, self.v0_static))
if static_range.flag:
print("\n*** Static EoS is from a restricted volume range:")
print("Minimum volume: %8.3f" % static_range.vmin)
print("Maximum volume: %8.3f" % static_range.vmax)
if p_stat.flag:
print("\n*** Static EoS from P(V) data ***")
print("Data points num: %3i" % p_stat.npoints)
print("Volume range: %8.4f, %8.4f (A^3)" % (p_stat.vmin, p_stat.vmax))
print("Pressure range: %5.2f, %5.2f (GPa)" % (p_stat.pmax, p_stat.pmin))
print("EoS -- K0: %6.2f (GPa), Kp: %4.2f, V0: %8.4f (A^3)" % (p_stat.k0,\
p_stat.kp, p_stat.v0))
print("Energy at V0: %12.9e (hartree)" % p_stat.e0)
if self.k0 != None:
print("\n** BM3 EoS from the last computation, at the temperature of %5.2f K **" % self.temp)
print("K0: %6.2f GPa, Kp: %4.2f, V0: %8.4f A^3" %\
(self.k0, self.kp, self.v0))
if not f_fix.flag:
print("Kp not fixed")
else:
print("Kp fixed")
if exclude.ex_mode != []:
uniq=np.unique(exclude.ex_mode)
print("\nZone center excluded modes: %s" % str(uniq))
else:
print("\nNo zone center excluded modes")
if zp.flag:
print("\nZero point energy is explicitly computed\n")
else:
print("\nZero point energy is not computed: it must be included")
print("in the static energy file\n")
if disp.ex_flag:
uniq=np.unique(disp.excluded_list)
print("Off center excluded modes: %s" % str(uniq))
else:
print("No off center excluded modes")
if kieffer.flag==True:
print("\nKieffer model on; frequencies %5.2f %5.2f %5.2f cm^-1" %\
(kieffer.kief_freq_inp[0], kieffer.kief_freq_inp[1], \
kieffer.kief_freq_inp[2]))
else:
print("\nKieffer model off")
if anharm.flag:
print("\nAnharmonic correction for mode(s) N. %s" % str(anharm.mode).strip('[]'))
print("Brillouin flag(s): %s" % str(anharm.brill).strip('[]'))
print("Anharmonic F(V,T) fit - degree V: %2i, degree T: %2i" % (anharm.vdeg, anharm.tdeg))
if disp.flag:
print("\n--------------- Phonon dispersion --------------------")
print("\nDispersion correction activated for the computation of entropy and")
print("specific heat:")
print("Number of frequency sets: %3i" % disp.nset)
if disp.nset > 1:
if disp.fit_type == 0:
print("Polynomial fit of the frequencies; degree: %3i " % disp.fit_degree)
else:
print("Spline fit of the frequencies; degree: %3i, smooth: %3.1f"\
% (disp.fit_degree, disp.fit_type))
print("Number of off-centered modes: %5i" % disp.f_size)
if disp.eos_flag:
print("\nThe phonon dispersion is used for the computation of the bulk modulus")
print("if the bulk_dir or the bulk_modulus_p functions are used, the latter")
print("in connection with the noeos option.")
if disp.fit_vt_flag:
print("The required V,T-fit of the free energy contribution from")
print("the off-centered modes is ready. Fit V,T-powers: %3i, %3i"
% (disp.fit_vt_deg_v, disp.fit_vt_deg_t))
else:
print("The required V,T-fit of the free energy contribution from")
print("the off-centered mode is NOT ready.")
else:
print("\nThe phonon dispersion correction is not used for the computation")
print("of the bulk modulus")
if disp.thermo_vt_flag & (disp.nset > 1):
print("\nVT-phonon dispersion correction to the thermodynamic properties")
elif (not disp.thermo_vt_flag) & (disp.nset > 1):
print("\nT-phonon dispersion correction to the thermodynamic properties")
print("Use disp.thermo_vt_on() to activate the V,T-correction")
print("\n --------------------------------------------------------")
if ac_approx.flag:
print("\nEstimation of the contribution to the vibrational pressure")
print("from acoustic phonon branches (no DISP computation) is activated\n")
print("Number of bins of the histogram of frequencies: %2i" % ac_approx.nbin)
print("\nLower frequency bin:")
print("Number of modes: %3i" % ac_approx.nmode)
print("Minimum frequency: %5.1f" % ac_approx.fmin)
print("Maximum frequency: %5.1f" % ac_approx.fmax)
print("Mean frequency: %5.1f" % ac_approx.fmean)
if lo.flag:
out_lo=(lo.mode, lo.split)
df_out=pd.DataFrame(out_lo, index=['Mode', 'Split'])
df_out=df_out.T
df_out['Mode']=np.array([int(x) for x in df_out['Mode']], dtype=object)
print("\nFrequencies corrected for LO-TO splitting.\n")
if verbose.flag:
print(df_out.to_string(index=False))
print("---------------------------------------------")
print("\n**** Volume driver for volume_dir function ****")
print("Delta: %3.1f; degree: %2i; left: %3.1f; right: %3.1f, Kp_fix: %s; t_max: %5.2f"\
% (volume_ctrl.delta, volume_ctrl.degree, volume_ctrl.left, volume_ctrl.right,\
volume_ctrl.kp_fix, volume_ctrl.t_max))
print("EoS shift: %3.1f; Quad_shrink: %2i; T_dump: %3.1f; Dump fact.: %2.1f, T_last %4.1f" % \
(volume_ctrl.shift, volume_ctrl.quad_shrink, volume_ctrl.t_dump, volume_ctrl.dump,\
volume_ctrl.t_last))
print("Upgrade shift: %r" % volume_ctrl.upgrade_shift)
print("\n**** Volume driver for volume_from_F function ****")
print("In addition to the attributes set in the parent volume_control_class:")
print("shift: %3.1f, flag: %r, upgrade_shift: %r" % (volume_F_ctrl.get_shift(), \
volume_F_ctrl.get_flag(), volume_F_ctrl.get_upgrade_status()))
print("\n**** Numerical T-derivatives driver class (delta_ctrl) ****")
if not delta_ctrl.adaptive:
print("Delta: %3.1f" % delta_ctrl.delta)
print("Degree: %3i" % delta_ctrl.degree)
print("N. of points %3i" % delta_ctrl.nump)
else:
print("Adaptive scheme active:")
print("T_min, T_max: %4.1f, %6.1f K" % (delta_ctrl.tmin, delta_ctrl.tmax))
print("Delta_min, Delta_max: %4.1f, %6.1f K" % (delta_ctrl.dmin, delta_ctrl.dmax))
print("Degree: %3i" % delta_ctrl.degree)
print("N. of points %3i" % delta_ctrl.nump)
if verbose.flag:
print("\n--------- Database section ---------")
print("Loaded phases:")
print(self.min_names)
class exclude_class():
"""
Contains the list of modes to be excluded from the
calculation of the Helmholtz free energy.
It can be constructed by using the keyword EXCLUDE
in the input.txt file.
"""
def __init__(self):
self.ex_mode=[]
self.ex_mode_keep=[]
self.flag=False
def __str__(self):
return "Excluded modes:" + str(self.ex_mode)
def add(self,modes):
"""
Args:
n : can be a scalar or a list of modes to be excluded
"""
if type(modes) is list:
self.ex_mode.extend(modes)
self.flag=True
elif type(modes) is int:
self.ex_mode.append(modes)
self.flag=True
else:
print("** Warning ** exclude.add(): invalid input type")
return
def restore(self):
"""
Restores all the excluded modes
"""
if self.flag:
self.ex_mode_keep=self.ex_mode
self.ex_mode=[]
self.flag=False
def on(self):
self.ex_mode=self.ex_mode_keep
self.flag=True
class fix_flag:
def __init__(self,value=0.):
self.value=value
self.flag=False
def on(self,value=4):
self.value=value
self.flag=True
def off(self):
self.value=0.
self.flag=False
class fit_flag:
def __init__(self):
pass
def on(self):
self.flag=True
def off(self):
self.flag=False
class spline_flag(fit_flag):
"""
Sets up the spline fit of the frequencies as functions of
the volume of the unit cell.
Several variables are defined:
1. flag: (boolean); if True, frequencies are fitted with splines
2. degree: degree of the spline
3. smooth: *smoothness* of the spline
4. flag_stack: (boolean) signals the presence of the spline stack
5. pol_stack: it is the stack containing parameters for the spline fit
Note:
The spline stack can be set up and initialized by using the keyword\
SPLINE under the keyword FITVOL in the *input.txt* file
Methods:
"""
def __init__(self,flag=False,degree=3,smooth=0):
super().__init__()
self.flag=False
self.flag_stack=False
self.degree=degree
self.smooth=smooth
self.pol_stack=np.array([])
def on(self):
super().on()
def off(self):
super().off()
def set_degree(self,degree):
self.degree=int(degree)
def set_smooth(self,smooth):
self.smooth=smooth
def stack(self):
self.pol_stack=freq_stack_spline()
self.flag_stack=True
def vol_range(self,v_ini, v_fin, npoint):
self.fit_vol=np.linspace(v_ini, v_fin, npoint)
class poly_flag(fit_flag):
def __init__(self,flag=False,degree=2):
super().__init__()
self.flag=flag
self.flag_stack=False
self.degree=degree
self.pol_stack=np.array([])
def on(self):
super().on()
def off(self):
super().off()
def set_degree(self,degree):
self.degree=int(degree)
def stack(self):
self.pol_stack=freq_stack_fit()
self.flag_stack=True
def vol_range(self,v_ini, v_fin, npoint):
self.fit_vol=np.linspace(v_ini, v_fin, npoint)
class kieffer_class():
def __str__(self):
return "Application of the Kieffer model for acoustic phonons"
def __init__(self,flag=False):
self.flag=False
self.stack_flag=False
self.kief_freq=None
self.kief_freq_inp=None
self.t_range=None
self.f_list=None
self.input=False
def stack(self, t_range, f_list):
self.t_range=t_range
self.f_list=f_list
def get_value(self,temperature):
free=scipy.interpolate.interp1d(self.t_range, self.f_list, kind='quadratic')
return free(temperature)*zu
def on(self):
self.flag=True
print("Kieffer correction on")
if disp.flag:
disp.flag=False
print("Phonon dispersion is deactivated")
if not self.stack_flag:
free_stack_t(pr.kt_init,pr.kt_fin,pr.kt_points)
def off(self):
self.flag=False
print("Kieffer correction off")
def freq(self,f1,f2,f3):
self.kief_freq_inp=np.array([f1, f2, f3])
self.kief_freq=self.kief_freq_inp*csl*h/kb
free_stack_t(pr.kt_init,pr.kt_fin,pr.kt_points)
def rescale(self, fact=1., save=False):
k0=self.kief_freq_inp[0]*fact
k1=self.kief_freq_inp[1]*fact
k2=self.kief_freq_inp[2]*fact
k_list=[k0, k1, k2]
k_print=list(round(ik,1) for ik in k_list)
print("Kieffer frequencies rescaled to ", k_print)
if save:
self.kief_freq_inp=np.array([k0, k1, k2])
self.kief_freq=np.array([k0, k1, k2])*csl*h/kb
free_stack_t(pr.kt_init,pr.kt_fin,pr.kt_points)
def plot(self):
plt.figure()
plt.plot(self.t_range, self.f_list, "k-")
plt.xlabel("Temperature (K)")
plt.ylabel("F free energy (J/mol apfu)")
plt.title("Free energy from acustic modes (Kieffer model)")
plt.show(block=False)
class acoustic_approx_class():
"""
Class used for the estimation of the contribution to the vibrational
pressure from acoustic modes, in case supercell data are not available.
Use the method 'on' to switch on the correction ('off' to switch it off)
The algorithm is used to compute such contribution if the total pressure
at given temperature and volume is evaluated by the pressure_dir function
Example:
>>> ac_approx.on()
>>> bulk_dir(298)
"""
def __init__(self):
self.flag=False
self.fit=None
self.deg=3
self.nt=12
self.tmax=1000.
self.tmin=20.
self.nbin=4
self.fit_flag=True
self.disp=False
def set(self, tmin=20., tmax=1000., deg=3, nt=12, nbin=4, on=True):
"""
Sets parameters for the estimation of the vibrational pressure
of low frequency modes, as a function of temperature
Args:
tmin, tmax, nt: minimum, maximum and number of points defining
the temperature range for the estimation of the
vibrational pressure as a function of T
deg: degree of the polynomial fitting the P(T) values
nbin: number of bins for the construction of the frequency
histogram
on: if True, the P(T) polynomial is computed after the set
method was invoked, and the acoustic correction is
switched on; if 'on' is False, the requested parameters are
changed but the correction is not switched on, or recomputed.
"""
self.tmin=tmin
self.tmax=tmax
self.deg=deg
self.nt=nt
self.nbin=nbin
if on:
self.on()
def on(self):
self.fit=acoustic_func()
if disp.flag:
print("\n**** Note ****")
disp.off()
self.disp=True
if self.fit_flag:
self.flag=True
else:
self.flag=False
print("The acoustic approximation requires the polynomial")
print("fits of the frequencies")
def off(self):
self.flag=False
if self.disp:
print("\n**** Note ****")
disp.on()
class bm4_class():
"""
Set up and information for a 4^ order Birch-Murnaghan EoS (BM4)
It provides:
1. energy: function; Volume integrated BM4 (V-BM4)
2. pressure: function; BM4
3. bm4_static_eos: BM4 parameters for the static energy
calculation as a function of V
4. en_ini: initial values for the BM4 fit
5. bm4_store: BM4 parameters from a fitting at a given
temperature
methods:
"""
def __init__(self):
self.flag=False
self.start=True
self.energy=None
self.pressure=None
self.en_ini=None
self.bm4_static_eos=None
self.bm4_store=None
def __str__(self):
return "BM4 setting: " + str(self.flag)
def on(self):
"""
Switches on the BM4 calculation
"""
self.flag=True
if self.start:
self.energy, self.pressure=bm4_def()
self.start=False
def estimates(self,v4,e4):
"""
Estimates initial values of BM4 parameters for the fit
"""
ini=init_bm4(v4,e4,4.0)
new_ini,dum=curve_fit(v_bm3, v4, e4, \
p0=ini,ftol=1e-15,xtol=1e-15)
kpp=(-1/new_ini[1])*((3.-new_ini[2])*\
(4.-new_ini[2])+35./9.)*1e-21/conv
self.en_ini=[new_ini[0], new_ini[1],\
new_ini[2], kpp, new_ini[3]]
k0_ini=new_ini[1]*conv/1e-21
print("\nBM4-EoS initial estimate:")
print("\nV0: %6.4f" % self.en_ini[0])
print("K0: %6.2f" % k0_ini)
print("Kp: %6.2f" % self.en_ini[2])
print("Kpp: %6.2f" % self.en_ini[3])
print("E0: %8.6e" % self.en_ini[4])
def store(self,bm4st):
"""
Stores BM4 parameters from a fit a given temperature
"""
self.bm4_store=bm4st
def upload(self,bm4_eos):
"""
Loads the parameters from the static calculation
(that are then stored in bm4_static_eos)
"""
self.bm4_static_eos=bm4_eos
def upgrade(self):
"""
Uses the stored values of parameters [from the application of
store()] to upgrade the initial estimation done with estimates()
"""
self.en_ini=self.bm4_store
def off(self):
"""
Switches off the BM4 calculation
"""
self.flag=False
def status(self):
"""
Informs on the status of BM4 (on, or off)
"""
print("\nBM4 setting: %s " % self.flag)
class gamma_class():
"""
Store coefficients of a gamma(T) fit
"""
def __init__(self):
self.flag=False
self.degree=1
self.pol=np.array([])
def upload(self,deg,pcoef):
self.flag=True
self.degree=deg
self.pol=pcoef
class super_class():
"""
Store supercell data: number of cells on which the frequencies
computation was done. To be used in connection with CRYSTAL
calculations performed with SCELPHONO and QHA keywords.
Default value: 1
"""
def __init__(self):
self.number=1
self.flag=False
def set(self,snum):
self.flag=True
self.number=snum
print("\n*** Supercell *** Number of cells: %3i" % snum)
def reset(self):
self.flag=False
self.number=1
print("\n*** Supercell deactivated *** Number of cells set to 1")
class lo_class():
"""
LO/TO splitting correction.
The class stores a copy of the original TO frequencies, the modes
affected by LO/TO splitting and the splitting values.
Modes are identified by their progressive number (starting from 0) stored
in the *mode* attribute.
When the correction is activated, new values of frequencies (*f_eff*)
are computed for the relevant modes, according to the formula:
f_eff = 2/3 f_TO + 1/3 f_LO
where f_LO = f_TO + split.
Correction is activated by the keyword LO in the input.txt file,
followed by the name of the file containing the splitting data (two
columns: mode number and the corresponding split in cm^-1).
Internally, the methods *on* and *off* switch respectively on and off
the correction. The method *apply* does the computation of the frequencies
*f_eff*.
"""
def __init__(self):
self.flag=False
self.mode=np.array([])
self.split=np.array([])
self.data_freq_orig=np.array([])
self.data_freq=np.array([])
def on(self):
self.apply()
if flag_spline.flag:
flag_spline.stack()
elif flag_poly.flag:
flag_poly.stack()
self.flag=True
print("Frequencies corrected for LO-TO splitting")
def off(self):
self.flag=False
self.data_freq=np.copy(self.data_freq_orig)
if flag_spline.flag:
flag_spline.stack()
elif flag_poly.flag:
flag_poly.stack()
print("LO-TO splitting not taken into account")
def apply(self):
for ifr in np.arange(lo.mode.size):
im=lo.mode[ifr]
for iv in int_set:
freq_lo=self.data_freq_orig[im,iv+1]+self.split[ifr]
self.data_freq[im,iv+1]=(2./3.)*self.data_freq_orig[im,iv+1]\
+(1./3.)*freq_lo
class anh_class():
def __init__(self):
self.flag=False
self.disp_off=0
def off(self):
self.flag=False
exclude.restore()
if disp.input_flag:
disp.free_exclude_restore()
print("Anharmonic correction is turned off")
print("Warning: all the excluded modes are restored")
def on(self):
self.flag=True
self.flag_brill=False
for im, ib in zip(anharm.mode, anharm.brill):
if ib == 0:
exclude.add([im])
elif disp.input_flag:
disp.free_exclude([im])
self.flag_brill=True
if self.flag_brill:
disp.free_fit_vt()
print("Anharmonic correction is turned on")
class static_class():
"""
Defines the volume range for the fit of the static EoS
If not specified (default) such range is defined from the
volumes found in the static energies file.
"""
def __init__(self):
self.flag=False
def set(self, vmin, vmax):
"""
Sets the minimum and maximum volumes for the V-range
Args:
vmin: minimum volume
vmax: maximum volume
"""
self.vmin=vmin
self.vmax=vmax
def off(self):
"""
Restores the original V-range (actually, it switches off the volume
selection for the fit of the static EoS)
"""
self.flag=False
def on(self):
"""
It switches on the volume selection for the fit of the static EoS
Note:
The minimum and maximum V-values are set by the 'set' method
of the class
"""
self.flag=True
class p_static_class():
def __init__(self):
self.flag=False
self.vmin=None
self.vmax=None
self.pmin=None
self.pmax=None
self.npoints=None
self.k0=None
self.kp=None
self.v0=None
self.e0=None
class volume_control_class():
"""
Defines suitable parameters for the volume_dir function
"""
def __init__(self):
self.degree=2
self.delta=2.
self.t_max=500.
self.shift=0.
self.t_dump=0.
self.dump=1.
self.quad_shrink=4
self.kp_fix=False
self.debug=False
self.upgrade_shift=False
self.skew=1.
self.t_last=0.
self.t_last_flag=False
self.v_last=None
def set_degree(self, degree):
"""
Sets the degree of polynomial used to fit the (P(V)-P0)^2 data.
The fitted curve is the minimized to get the equilibrium volume
at each T and P.
For each of the single parameter revelant in this class, there exist
a specific method to set its value. The method set_all can be used to
set the values of a number of that, at the same time, by using appropriate
keywords as argument. The arguments to set_all are:
Args:
degree: degree of the fitting polynomial (default=2)
delta: volume range where the minimum of the fitting function
is to be searched (default=2.)
skew: the Volume range is centered around the equilibrium
volume approximated by the EoS-based new_volume function
The symmetry around such point can be controlled by
the skew parameter (default=1.: symmetric interval)
shift: Systematic shift from the new_volume estimation (default=0.)
t_max: In the initial estimation of the volume at P/T with the EoS-based
new_volume function, the Kp is refined if T < t_max.
If T > t_max and kp_fix=True, Kp is fixed at the value
refined at t_max (default=500K)
kp_fix: See t_max (default=True)
quad_shrink: if degree=2, it restricts the volume range around the
approximated volume found. The new range is
delta/quad_shrink (default=4)
upgrade_shift: at the end of the computation, the difference between
the volume found and the initial one (from the EoS-
based new_volume function) is calculated. The shift
attribute is then upgraded if upgrade_shift is True
(default=False)
debug: if True, the (P(V)-P0)**2 function is plotted as a function
of V (default=False)
t_dump: temperature over which a dumping on the shift parameter is
applied (default=0.)
dump: dumping on the shift parameter (shift=shift/dump; default=1.)
t_last: if t_last > 10., the last volume computed is used as the
initial guess value (vini) for the next computation at a
new temperature.
"""
self.degree=degree
def set_delta(self, delta):
self.delta=delta
def set_tmax(self,tmax):
self.t_max=tmax
def set_skew(self, skew):
self.left=skew+1
self.right=(skew+1)/skew
def kp_on(self):
self.kp_fix=True
def kp_off(self):
self.kp_fix=False
def debug_on(self):
self.debug=True
def debug_off(self):
self.debug=False
def set_shift(self, shift):
self.shift=shift
def upgrade_shift_on(self):
self.upgrade_shift=True
def upgrade_shift_off(self):
self.ugrade_shift=False
def set_shrink(self, shrink):
self.quad_shrink=shrink
def shift_reset(self):
self.shift=0.
def set_t_dump(self,t_dump=0., dump=1.0):
self.t_dump=t_dump
self.dump=dump
def set_t_last(self, t_last):
self.t_last=t_last
def set_all(self,degree=2, delta=2., skew=1., shift=0., t_max=500.,\
quad_shrink=4, kp_fix=True, upgrade_shift=False, debug=False,\
t_dump=0., dump=1., t_last=0.):
self.degree=degree
self.delta=delta
self.t_max=t_max
self.kp_fix=kp_fix
self.debug=debug
self.left=skew+1
self.right=(skew+1)/skew
self.shift=shift
self.quad_shrink=quad_shrink
self.upgrade_shift=upgrade_shift
self.skew=skew
self.t_last=t_last
class volume_F_control_class():
"""
Class controlling some parameters relevant for the computation of
volume and thermal expansion by using the volume_from_F function.
Precisely, the initial volume (around which the refined volume vref
is to be searched) is set to vini+shift, where vini is the
output from the volume_dir, whereas shift is from this class.
Shift is computed as the difference vref-vini; it can be upgraded
provided the flag upgrade_shift is set to True.
"""
def __init__(self):
self.shift=0.
self.upgrade_shift=False
self.flag=False
def on(self):
self.flag=True
def off(self):
self.flag=False
def set_shift(self, sh):
self.shift=sh
def upgrade_on(self):
self.upgrade_shift=True
def upgrade_off(self):
self.upgrade_shift=False
def get_shift(self):
return self.shift
def get_upgrade_status(self):
return self.upgrade_shift
def get_flag(self):
return self.flag
class delta_class():
"""
Control parameters for the numerical evaluation of the first and second
derivatives of the Helmholtz free energy as a function of T. They are
relevant for the entropy_v function that computes both the entropy and
specific heat at a fixed volume, as well as the computation of thermal
expansion.
Initial values of delta, degree and number of points are read
from the parameters file 'parame.py'
New values can be set by the methods set_delta, set_degree and set_nump
of the class. values can be retrieved by the corresponding 'get' methods.