Skip to content

Commit c7f0a76

Browse files
Merge pull request #72 from SSCHAcode/tar_file_cluster
Now the cluster automatic submission to export tar files
2 parents 85122c8 + db710d0 commit c7f0a76

10 files changed

Lines changed: 519 additions & 109 deletions

File tree

Modules/Cluster.py

Lines changed: 285 additions & 91 deletions
Large diffs are not rendered by default.

Modules/Ensemble.py

Lines changed: 41 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1639,6 +1639,39 @@ def get_free_energy(self, return_error = False):
16391639
return free_energy, error
16401640
return free_energy
16411641

1642+
def get_entropy(self, return_error = False):
1643+
r"""
1644+
GET THE ENTROPY
1645+
===============
1646+
1647+
Get the total anharmonic entropy.
1648+
1649+
The equation implemented is the analytical derivative of the free energy,
1650+
and assumes that the SSCHA free energy is minimized.
1651+
1652+
.. math::
1653+
1654+
S = - \frac{dF}{dT}
1655+
1656+
S = S_{harm} - \left<V - {\mathcal V}\right>\sum_\mu \frac{1}{1 + 2n_\mu} \frac{dn_\mu}{dT}
1657+
1658+
1659+
where :math:`S_{harm}` is the 'harmonic' entropy computed from the dynamucal matrix,
1660+
plus a correction accounting for the ensemble anharmonicity.
1661+
1662+
1663+
Parameters
1664+
----------
1665+
return_error : bool
1666+
If true, returns also the error
1667+
1668+
Results
1669+
-------
1670+
entropy, error : float
1671+
Returns the entropy and [optionally] the stochastic error.
1672+
"""
1673+
raise NotImplementedError("Error, to be implemented")
1674+
16421675

16431676
def get_free_energy_interpolating(self, target_supercell, support_dyn_coarse = None, support_dyn_fine = None, error_on_imaginary_frequency = True, return_error = False):
16441677
"""
@@ -2694,7 +2727,7 @@ def get_dynamical_bubble(self, q, w, smearing = 1e-5):
26942727

26952728

26962729
def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, verbose = False, \
2697-
use_symmetries = True):
2730+
use_symmetries = True, return_d3 = False):
26982731
"""
26992732
GET THE FREE ENERGY ODD CORRECTION
27002733
==================================
@@ -2719,11 +2752,16 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v
27192752
use_symmetries : bool
27202753
If true, the d3 and d4 are symmetrized in real space.
27212754
It requires that spglib is installed to detect symmetries in the supercell correctly.
2755+
return_d3 : bool
2756+
If true, returns also the tensor of three phonon scattering.
27222757
27232758
Returns
27242759
-------
27252760
phi_sc : Phonons()
27262761
The dynamical matrix of the free energy hessian in (Ry/bohr^2)
2762+
d3 : ndarray (size = (3*nat_sc, 3*nat_sc, 3*nat_sc), Optional
2763+
Return the three-phonon-scattering tensor (in Ry atomic units).
2764+
Only if return_d3 is True.
27272765
"""
27282766
# For now the v4 is not implemented
27292767
# if include_v4:
@@ -2859,6 +2897,8 @@ def get_free_energy_hessian(self, include_v4 = False, get_full_hessian = True, v
28592897
dyn_hessian.dynmats[iq] = dynq_odd[iq, :, :]
28602898

28612899

2900+
if return_d3:
2901+
return dyn_hessian, d3* 2.0 # Ha to Ry
28622902
return dyn_hessian
28632903

28642904

Modules/Minimizer.py

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -187,6 +187,16 @@ def get_dyn_struct(self):
187187
struct = self.current_x[self.nq * self.n_modes * self.n_modes :].reshape(self.dyn.structure.coords.shape)
188188
return dynq, struct
189189

190+
def is_new_direction(self):
191+
"""
192+
Return if a new direction must be chosen.
193+
"""
194+
195+
if self.fixed_step:
196+
return True
197+
return self.new_direction
198+
199+
190200
def run_step(self, gradient, kl_new):
191201
"""
192202
Perform the minimization step with the line minimization

Modules/Optimizer.py

Lines changed: 23 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,9 @@ def __init__(self, starting_unit_cell):
4646
self.uc_0_inv = np.linalg.inv(self.uc_0)
4747
self.x_start = None
4848
self.reset_strain = False
49+
50+
# A control parameter for the line minimization (avoid failures)
51+
self.last_factor = 0
4952

5053
self.uc_old = None
5154

@@ -105,7 +108,8 @@ def get_line_step(self, grad):
105108
sys.stdout.flush()
106109

107110
print("[CELL] GRADIENT DOT DIRECTION = {}".format(factor))
108-
111+
failure_lmin = False
112+
109113
sys.stdout.flush()
110114
# Regularization (avoid run away)
111115
if factor > 2:
@@ -119,20 +123,36 @@ def get_line_step(self, grad):
119123
print(" incrementing by a 15 %")
120124

121125
elif factor < self.min_alpha_factor:
126+
# Check if the last time was good
127+
if self.last_factor != 0:
128+
if factor < self.last_factor:
129+
print("[CELL] Warning: reduced step and got worse; restarting line minimization with a good step.")
130+
failure_lmin = True
131+
self.last_factor = factor
122132
factor = self.min_alpha_factor
123133
good_step = False
124134

125135
if np.isnan(factor):
136+
self.last_factor = 0
126137
factor = self.min_alpha_factor
127138
good_step = False
128-
129-
self.alpha *= factor
139+
140+
# Check if the line minimization got contradditing result
141+
# In this case stop the line minimization and proceed with a standard gradient descend.
142+
if failure_lmin:
143+
good_step = True
144+
else:
145+
self.alpha *= factor
130146

131147

132148
#if self.alpha < self.min_alpha_step * self.alpha0:
133149
# self.alpha = self.min_alpha_step * self.alpha0
134150
#self.alpha *= factor
135151

152+
# Reset the check on the line minimization if this is a good step
153+
if good_step:
154+
self.last_factor = 0
155+
136156
return good_step
137157

138158
def reset(self, unit_cell):

Modules/Relax.py

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -536,6 +536,10 @@ def vc_relax(self, target_press = 0, static_bulk_modulus = 100,
536536
self.minim.ensemble.dyn_0 = self.minim.dyn.Copy()
537537
if pop != start_pop or not restart_from_ens:
538538
self.minim.ensemble.generate(self.N_configs)
539+
540+
# Save also the generation
541+
#if ensemble_loc is not None and self.save_ensemble:
542+
# self.minim.ensemble.save_bin(ensemble_loc, pop)
539543

540544
# Compute energies and forces
541545
self.minim.ensemble.compute_ensemble(self.calc, True, stress_numerical,

Modules/SchaMinimizer.py

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -680,7 +680,7 @@ def setup_from_namelist(self, input_file):
680680

681681
# Symmetrize the dynmat if requested
682682
if not self.neglect_symmetries:
683-
self.dyn.Symmetrize()
683+
self.dyn.Symmetrize(use_spglib = self.use_spglib)
684684

685685
if not __SCHA_DATADIR__ in keys:
686686
self.ensemble = Ensemble.Ensemble(self.dyn, 0)
@@ -1426,8 +1426,8 @@ def check_stop(self):
14261426

14271427
# Check the KL
14281428
kl = self.ensemble.get_effective_sample_size()
1429-
1430-
if kl / float(self.ensemble.N) < self.kong_liu_ratio and self.minimizer.new_direction:
1429+
1430+
if kl / float(self.ensemble.N) < self.kong_liu_ratio and self.minimizer.is_new_direction():
14311431
self.__converged__ = False
14321432
print ("KL:", kl, "KL/N:", kl / float(self.ensemble.N), "KL RAT:", self.kong_liu_ratio)
14331433
print (" According to your input criteria")

Modules/Utilities.py

Lines changed: 36 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -456,6 +456,7 @@ def __init__(self):
456456
self.save_atomic_positions = False
457457
self.__save_each_step = False
458458
self.current_struct = None
459+
self.minim_data = []
459460

460461
def Reset(self):
461462
"""
@@ -492,10 +493,15 @@ def SetupSaving(self, fname, save_each_step = True):
492493
"""
493494
Setup the system to save the data each time the function is called.
494495
496+
By default, the system will save frequencies and minimization info (like free energy, gradients and kong liu)
497+
The files are
498+
.freqs for the frequencies
499+
.dat for the minimization info
500+
495501
Parameters
496502
----------
497503
fname : string
498-
path to the file to save the frequencies vs time
504+
path to the file to save the files
499505
save_each_step : bool
500506
If true the file is saved (and updated) each time step.
501507
@@ -504,7 +510,7 @@ def SetupSaving(self, fname, save_each_step = True):
504510
self.__save_fname = fname
505511
self.__save_each_step = save_each_step
506512

507-
def Save(self, fname= None):
513+
def Save(self, fname = None):
508514
"""
509515
Save the data on a file
510516
@@ -518,12 +524,22 @@ def Save(self, fname= None):
518524
if not sscha.Parallel.am_i_the_master():
519525
return
520526

521-
if fname is None:
522-
if self.__save_fname is None:
523-
raise IOError("Error, a filename must be specified to save the frequencies.")
524-
np.savetxt(self.__save_fname, self.total_freqs, header = "Time vs Frequencies")
525-
else:
526-
np.savetxt(fname, self.total_freqs, header = "Time vs Frequencies")
527+
root_name = fname
528+
529+
if root_name is None:
530+
root_name = self.__save_fname
531+
532+
if root_name is None:
533+
raise IOError("Error, a filename must be specified to save the data.")
534+
535+
print("ROOT NAME:", root_name)
536+
print("SAVE NAME:", self.__save_fname)
537+
print("FNAME NAME:", fname)
538+
freq_name = root_name + '.freqs'
539+
data_name = root_name + '.dat'
540+
541+
np.savetxt(freq_name, self.total_freqs, header = "Time vs Frequencies")
542+
np.savetxt(data_name, self.minim_data, header = '# Free energy [meV] +- error; FC gradient +- error; Structure gradient +- error; Kong-Liu effective sample size')
527543

528544

529545
if self.save_weights:
@@ -547,7 +563,7 @@ def CFP_SaveAll(self, minim):
547563
if not sscha.Parallel.am_i_the_master():
548564
return
549565

550-
if minim.minimizer.new_direction:
566+
if minim.minimizer.is_new_direction():
551567

552568
# Get the weights if required
553569
if self.save_weights:
@@ -558,6 +574,17 @@ def CFP_SaveAll(self, minim):
558574

559575
if self.save_atomic_positions:
560576
self.current_struct = minim.dyn.structure.copy()
577+
578+
# Get the minimization data
579+
fe = minim.__fe__[-1] * sscha.SchaMinimizer.__RyTomev__
580+
fe_err = minim.__fe_err__[-1] * sscha.SchaMinimizer.__RyTomev__
581+
gc = minim.__gc__[-1] * sscha.SchaMinimizer.__RyTomev__
582+
gc_err = minim.__gc_err__[-1] * sscha.SchaMinimizer.__RyTomev__
583+
gw = minim.__gw__[-1] * sscha.SchaMinimizer.__RyTomev__
584+
gw_err = minim.__gw_err__[-1] * sscha.SchaMinimizer.__RyTomev__
585+
kl = minim.__KL__[-1]
586+
587+
self.minim_data.append([fe, fe_err, gc, gc_err, gw, gw_err, kl])
561588

562589
# This perform also the saving
563590
self.CFP_SaveFrequencies(minim)

scripts/plot_frequencies.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77

88
RY_TO_CM = 109691.40235
99

10-
print()
10+
print('Obsolete since v1.2, please use sscha-plot-data.py')
1111
print()
1212
print("PLOTTING THE FREQIENCIES")
1313
print("========================")

scripts/sscha-plot-data.py

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#!python
2+
3+
import numpy as np
4+
import matplotlib.pyplot as plt
5+
import sys, os
6+
7+
import cellconstructor as CC, cellconstructor.Units
8+
9+
10+
DESCRIPTION = '''
11+
This program plots the results of a SSCHA minimization
12+
saved with the Utilities.IOInfo class.
13+
14+
It plots both frequencies and the free energy minimization.
15+
16+
Many files can be specified, in this case, they are plotted one after the other,
17+
separated by vertical dashed lines.
18+
19+
Usage:
20+
21+
>>> sscha-plot-data.py file1 ...
22+
23+
The code looks for data file called file.dat and file.freqs,
24+
containing, respectively, the minimization data and the auxiliary frequencies.
25+
26+
These files are generated only by python-sscha >= 1.2.
27+
28+
'''
29+
30+
31+
def main():
32+
print(DESCRIPTION)
33+
34+
if len(sys.argv) < 2:
35+
ERROR_MSG = '''
36+
Error, you need to specify at least one file.
37+
Exit failure!
38+
'''
39+
print(ERROR_MSG)
40+
raise ValueError(ERROR_MSG)
41+
42+
plt.rcParams["font.family"] = "Liberation Serif"
43+
LBL_FS = 12
44+
DPI = 120
45+
TITLE_FS = 15
46+
47+
freqs_files = [sys.argv[x] + '.freqs' for x in range(1, len(sys.argv))]
48+
minim_files = [sys.argv[x] + '.dat' for x in range(1, len(sys.argv))]
49+
50+
# Check if all the files exist
51+
for f, m in zip(freqs_files, minim_files):
52+
assert os.path.exists(f), 'Error, file {} not found.'.format(f)
53+
assert os.path.exists(m), 'Error, file {} not found.'.format(m)
54+
55+
print("Reading the input...")
56+
57+
# Load all the data
58+
freqs_data = np.concatenate([np.loadtxt(f) for f in freqs_files])
59+
minim_data = np.concatenate([np.loadtxt(f) for f in minim_files])
60+
61+
print("Plotting...")
62+
63+
# Insert the x axis in the plotting data
64+
xsteps = np.arange(minim_data.shape[0])
65+
new_data = np.zeros(( len(xsteps), 8), dtype = np.double)
66+
new_data[:,0] = xsteps
67+
new_data[:, 1:] = minim_data
68+
minim_data = new_data
69+
70+
fig_data, axarr = plt.subplots(nrows=2, ncols = 2, sharex = True, dpi = DPI)
71+
72+
# Plot the data
73+
axarr[0,0].fill_between(minim_data[:,0], minim_data[:,1] - minim_data[:, 2]*.5 ,
74+
minim_data[:, 1] + minim_data[:, 2] * .5, color = "aquamarine")
75+
axarr[0,0].plot(minim_data[:,0], minim_data[:,1], color = "k")
76+
axarr[0,0].set_ylabel("Free energy / unit cell [meV]", fontsize = LBL_FS)
77+
78+
79+
axarr[0,1].fill_between(minim_data[:,0], minim_data[:,3] - minim_data[:, 4]*.5 ,
80+
minim_data[:, 3] + minim_data[:, 4] * .5, color = "aquamarine")
81+
axarr[0,1].plot(minim_data[:,0], minim_data[:,3], color = "k")
82+
axarr[0,1].set_ylabel("FC gradient", fontsize = LBL_FS)
83+
84+
axarr[1,1].fill_between(minim_data[:,0], minim_data[:,5] - minim_data[:, 6]*.5 ,
85+
minim_data[:, 5] + minim_data[:, 6] * .5, color = "aquamarine")
86+
axarr[1,1].plot(minim_data[:,0], minim_data[:,5], color = "k")
87+
axarr[1,1].set_ylabel("Structure gradient", fontsize = LBL_FS)
88+
axarr[1,1].set_xlabel("Good minimization steps", fontsize = LBL_FS)
89+
90+
91+
axarr[1,0].plot(minim_data[:,0], minim_data[:,7], color = "k")
92+
axarr[1,0].set_ylabel("Effective sample size", fontsize = LBL_FS)
93+
axarr[1,0].set_xlabel("Good minimization steps", fontsize = LBL_FS)
94+
fig_data.tight_layout()
95+
96+
97+
# Now plot the frequencies
98+
fig_freqs = plt.figure(dpi = DPI)
99+
ax = plt.gca()
100+
101+
N_points, Nw = np.shape(freqs_data)
102+
103+
for i in range(Nw):
104+
ax.plot(freqs_data[:, i] * CC.Units.RY_TO_CM)
105+
106+
ax.set_xlabel("Good minimization steps", fontsize = LBL_FS)
107+
ax.set_ylabel("Frequency [cm-1]", fontsize = LBL_FS)
108+
ax.set_title("Frequcency evolution", fontsize = TITLE_FS)
109+
fig_freqs.tight_layout()
110+
111+
plt.show()
112+
113+
if __name__ == "__main__":
114+
main()

0 commit comments

Comments
 (0)