|
| 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 | +import time |
| 10 | +import cellconstructor.Methods |
| 11 | + |
| 12 | +import sscha, sscha.SchaMinimizer |
| 13 | +import sscha.Ensemble |
| 14 | +import sscha.Relax |
| 15 | +import sscha.Utilities |
| 16 | +import argparse |
| 17 | + |
| 18 | + |
| 19 | +""" |
| 20 | +This is part of the program python-sscha |
| 21 | +Copyright (C) 2018 Lorenzo Monacelli |
| 22 | +
|
| 23 | +This program is free software: you can redistribute it and/or modify |
| 24 | +it under the terms of the GNU General Public License as published by |
| 25 | +the Free Software Foundation, either version 3 of the License, or |
| 26 | +(at your option) any later version. |
| 27 | +
|
| 28 | +This program is distributed in the hope that it will be useful, |
| 29 | +but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 30 | +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 31 | +GNU General Public License for more details. |
| 32 | +
|
| 33 | +You should have received a copy of the GNU General Public License |
| 34 | +along with this program. If not, see <https://www.gnu.org/licenses/>. |
| 35 | +""" |
| 36 | + |
| 37 | +DESCRIPTION = ''' |
| 38 | +This program plots the results of a SSCHA minimization |
| 39 | +saved with the Utilities.IOInfo class. |
| 40 | +
|
| 41 | +It plots both frequencies and the free energy minimization. |
| 42 | +
|
| 43 | +Many files can be specified, in this case, they are plotted one after the other, |
| 44 | +separated by vertical dashed lines. |
| 45 | +
|
| 46 | +Usage: |
| 47 | + |
| 48 | + >>> sscha-plot-data.x file1 ... |
| 49 | +
|
| 50 | +The code looks for data file called file.dat and file.freqs, |
| 51 | +containing, respectively, the minimization data and the auxiliary frequencies. |
| 52 | +
|
| 53 | +These files are generated only by python-sscha >= 1.2. |
| 54 | +
|
| 55 | +''' |
| 56 | + |
| 57 | + |
| 58 | +def plot_data(): |
| 59 | + print(DESCRIPTION) |
| 60 | + |
| 61 | + if len(sys.argv) < 2: |
| 62 | + ERROR_MSG = ''' |
| 63 | +Error, you need to specify at least one file. |
| 64 | +Exit failure! |
| 65 | +''' |
| 66 | + print(ERROR_MSG) |
| 67 | + return 1 |
| 68 | + #raise ValueError(ERROR_MSG) |
| 69 | + |
| 70 | + plt.rcParams["font.family"] = "Liberation Serif" |
| 71 | + LBL_FS = 12 |
| 72 | + DPI = 120 |
| 73 | + TITLE_FS = 15 |
| 74 | + |
| 75 | + freqs_files = [sys.argv[x] + '.freqs' for x in range(1, len(sys.argv))] |
| 76 | + minim_files = [sys.argv[x] + '.dat' for x in range(1, len(sys.argv))] |
| 77 | + |
| 78 | + # Check if all the files exist |
| 79 | + plot_frequencies = None |
| 80 | + plot_minim = None |
| 81 | + for f, m in zip(freqs_files, minim_files): |
| 82 | + if os.path.exists(f): |
| 83 | + if plot_frequencies is None: |
| 84 | + plot_frequencies = True |
| 85 | + if not plot_frequencies: |
| 86 | + raise IOError("Error, file {} found, but not others .freqs".format(f)) |
| 87 | + else: |
| 88 | + if plot_frequencies: |
| 89 | + raise IOError("Error, file {} not found.".format(f)) |
| 90 | + plot_frequencies = False |
| 91 | + |
| 92 | + |
| 93 | + if os.path.exists(m): |
| 94 | + if plot_minim is None: |
| 95 | + plot_minim = True |
| 96 | + if not plot_minim: |
| 97 | + raise IOError("Error, file {} found, but not others .freqs".format(m)) |
| 98 | + else: |
| 99 | + if plot_minim: |
| 100 | + raise IOError("Error, file {} not found.".format(m)) |
| 101 | + plot_minim = False |
| 102 | + |
| 103 | + if not plot_frequencies and not plot_minim: |
| 104 | + print("Nothing to plot, check if the .dat and .freqs files exist.") |
| 105 | + exit() |
| 106 | + |
| 107 | + if plot_minim: |
| 108 | + print("Preparing the minimization data...") |
| 109 | + minim_data = np.concatenate([np.loadtxt(f) for f in minim_files]) |
| 110 | + |
| 111 | + # Insert the x axis in the plotting data |
| 112 | + xsteps = np.arange(minim_data.shape[0]) |
| 113 | + new_data = np.zeros(( len(xsteps), 8), dtype = np.double) |
| 114 | + new_data[:,0] = xsteps |
| 115 | + new_data[:, 1:] = minim_data |
| 116 | + minim_data = new_data |
| 117 | + |
| 118 | + fig_data, axarr = plt.subplots(nrows=2, ncols = 2, sharex = True, dpi = DPI) |
| 119 | + |
| 120 | + # Plot the data |
| 121 | + axarr[0,0].fill_between(minim_data[:,0], minim_data[:,1] - minim_data[:, 2]*.5 , |
| 122 | + minim_data[:, 1] + minim_data[:, 2] * .5, color = "aquamarine") |
| 123 | + axarr[0,0].plot(minim_data[:,0], minim_data[:,1], color = "k") |
| 124 | + axarr[0,0].set_ylabel("Free energy / unit cell [meV]", fontsize = LBL_FS) |
| 125 | + |
| 126 | + |
| 127 | + axarr[0,1].fill_between(minim_data[:,0], minim_data[:,3] - minim_data[:, 4]*.5 , |
| 128 | + minim_data[:, 3] + minim_data[:, 4] * .5, color = "aquamarine") |
| 129 | + axarr[0,1].plot(minim_data[:,0], minim_data[:,3], color = "k") |
| 130 | + axarr[0,1].set_ylabel("FC gradient", fontsize = LBL_FS) |
| 131 | + |
| 132 | + axarr[1,1].fill_between(minim_data[:,0], minim_data[:,5] - minim_data[:, 6]*.5 , |
| 133 | + minim_data[:, 5] + minim_data[:, 6] * .5, color = "aquamarine") |
| 134 | + axarr[1,1].plot(minim_data[:,0], minim_data[:,5], color = "k") |
| 135 | + axarr[1,1].set_ylabel("Structure gradient", fontsize = LBL_FS) |
| 136 | + axarr[1,1].set_xlabel("Good minimization steps", fontsize = LBL_FS) |
| 137 | + |
| 138 | + |
| 139 | + axarr[1,0].plot(minim_data[:,0], minim_data[:,7], color = "k") |
| 140 | + axarr[1,0].set_ylabel("Effective sample size", fontsize = LBL_FS) |
| 141 | + axarr[1,0].set_xlabel("Good minimization steps", fontsize = LBL_FS) |
| 142 | + fig_data.tight_layout() |
| 143 | + |
| 144 | + |
| 145 | + |
| 146 | + |
| 147 | + if plot_frequencies: |
| 148 | + print("Plotting the frequencies") |
| 149 | + |
| 150 | + # Load all the data |
| 151 | + freqs_data = np.concatenate([np.loadtxt(f) for f in freqs_files]) |
| 152 | + # Now plot the frequencies |
| 153 | + fig_freqs = plt.figure(dpi = DPI) |
| 154 | + ax = plt.gca() |
| 155 | + |
| 156 | + N_points, Nw = np.shape(freqs_data) |
| 157 | + |
| 158 | + for i in range(Nw): |
| 159 | + ax.plot(freqs_data[:, i] * CC.Units.RY_TO_CM) |
| 160 | + |
| 161 | + ax.set_xlabel("Good minimization steps", fontsize = LBL_FS) |
| 162 | + ax.set_ylabel("Frequency [cm-1]", fontsize = LBL_FS) |
| 163 | + ax.set_title("Frequcency evolution", fontsize = TITLE_FS) |
| 164 | + fig_freqs.tight_layout() |
| 165 | + |
| 166 | + plt.show() |
| 167 | + return 0 |
| 168 | + |
| 169 | +def main(): |
| 170 | + PROG_NAME = "sscha" |
| 171 | + VERSION = "1.5" |
| 172 | + DESCRIPTION = r""" |
| 173 | +
|
| 174 | +
|
| 175 | +
|
| 176 | + * * * * * * * * * * * * * * * * * * * * * * * * * * * |
| 177 | + * * |
| 178 | + * STOCHASTIC SELF-CONSISTENT HARMONIC APPROXIMATION * |
| 179 | + * * |
| 180 | + * * * * * * * * * * * * * * * * * * * * * * * * * * * |
| 181 | +
|
| 182 | +
|
| 183 | +
|
| 184 | + Wellcome to the SSCHA code. |
| 185 | + This stand-alone script performs SSCHA minimization |
| 186 | + according to the options parsed from the input file. |
| 187 | + Note, you can also run the SSCHA with the python-APIs. |
| 188 | + The python-APIs give access to much more customizable options, |
| 189 | + not available with the stand-alone app. |
| 190 | +
|
| 191 | + Refer to the documentation for more details. |
| 192 | +
|
| 193 | +
|
| 194 | + """ |
| 195 | + |
| 196 | + |
| 197 | + # Define a custom function to raise an error if the file does not exist |
| 198 | + def is_valid_file(parser, arg): |
| 199 | + if not os.path.exists(arg): |
| 200 | + parser.error("The file %s does not exists!" % arg) |
| 201 | + else: |
| 202 | + return arg |
| 203 | + |
| 204 | + # Prepare the parser for the command line |
| 205 | + parser = argparse.ArgumentParser(prog = PROG_NAME, |
| 206 | + description = DESCRIPTION, |
| 207 | + formatter_class = argparse.RawTextHelpFormatter) |
| 208 | + |
| 209 | + # Add the arguments |
| 210 | + parser.add_argument("-i", help="Path to the input file", |
| 211 | + dest="filename", required = True, metavar = "FILE", |
| 212 | + type = lambda x : is_valid_file(parser, x)) |
| 213 | + parser.add_argument("--plot-results", |
| 214 | + help="Plot the results", |
| 215 | + dest="plot", action = "store_true") |
| 216 | + |
| 217 | + parser.add_argument("--save-data", dest="save_data", metavar = "FILE", |
| 218 | + help = "Save the minimization details") |
| 219 | + |
| 220 | + |
| 221 | + # Get the arguments from the command line |
| 222 | + args = parser.parse_args() |
| 223 | + |
| 224 | + # Get the input filename and initialize the sscha minimizer. |
| 225 | + minim = sscha.SchaMinimizer.SSCHA_Minimizer() |
| 226 | + |
| 227 | + # Open the input file |
| 228 | + fp = open(args.filename, "r") |
| 229 | + line_list = fp.readlines() |
| 230 | + fp.close() |
| 231 | + |
| 232 | + # Check if the mandatory arguments are defined |
| 233 | + namelist = CC.Methods.read_namelist(line_list) |
| 234 | + if not sscha.SchaMinimizer.__SCHA_NAMELIST__ in namelist.keys(): |
| 235 | + print("Error, %s namelist required." % sscha.SchaMinimizer.__SCHA_NAMELIST__) |
| 236 | + return 1 |
| 237 | + |
| 238 | + # Check if the relax is hypothized |
| 239 | + if sscha.Relax.__RELAX_NAMESPACE__ in namelist.keys(): |
| 240 | + relax = sscha.Relax.SSCHA() |
| 241 | + print ("Found a %s namespace!" % sscha.Relax.__RELAX_NAMESPACE__) |
| 242 | + print ("Performing a whole relaxation.") |
| 243 | + |
| 244 | + relax.setup_from_namelist(namelist) |
| 245 | + |
| 246 | + # Save the final dynamical matrix |
| 247 | + print ("Saving the final dynamical matrix in %s" % (relax.minim.dyn_path + "_popuation%d_" % relax.minim.population)) |
| 248 | + relax.minim.dyn.save_qe(relax.minim.dyn_path + "_population%d_" % minim.population) |
| 249 | + |
| 250 | + # Plot the results |
| 251 | + relax.minim.plot_results(args.save_data, args.plot) |
| 252 | + |
| 253 | + exit(0) |
| 254 | + |
| 255 | + |
| 256 | + # Parse the input file |
| 257 | + t1 = time.time() |
| 258 | + print ("Loading everything...") |
| 259 | + minim.setup_from_namelist(namelist) |
| 260 | + t2 = time.time() |
| 261 | + print ("Loaded in %.4f seconds." % (t2 - t1)) |
| 262 | + |
| 263 | + # Start the minimization |
| 264 | + print () |
| 265 | + print ("Initialization...") |
| 266 | + minim.init() |
| 267 | + print ("System initialized.") |
| 268 | + minim.print_info() |
| 269 | + print () |
| 270 | + |
| 271 | + # If present, setup the custom functions |
| 272 | + if sscha.Utilities.__UTILS_NAMESPACE__ in namelist: |
| 273 | + cfs = sscha.Utilities.get_custom_functions_from_namelist(namelist, minim.dyn) |
| 274 | + minim.run(custom_function_pre = cfs[0], |
| 275 | + custom_function_gradient = cfs[1], |
| 276 | + custom_function_post = cfs[2]) |
| 277 | + else: |
| 278 | + minim.run() |
| 279 | + |
| 280 | + print () |
| 281 | + print ("Minimization ended.") |
| 282 | + minim.finalize(2) |
| 283 | + |
| 284 | + # Save the final dynamical matrix |
| 285 | + print ("Saving the final dynamical matrix in %s" % ('final_sscha_dyn_pop%d' % minim.population)) |
| 286 | + print ('Saving the final centroids into {}'.format('final_sscha_structure_pop%d.scf' % minim.population)) |
| 287 | + print () |
| 288 | + minim.dyn.save_qe('final_sscha_dyn_pop%d' % minim.population) |
| 289 | + minim.dyn.structure.save_scf('final_sscha_structure_pop%d.scf' % minim.population) |
| 290 | + |
| 291 | + # Plot the results |
| 292 | + minim.plot_results(args.save_data, args.plot) |
| 293 | + return 0 |
0 commit comments