From c73d6cbcf0b4d5fd59a03a019fcd47d9d435b799 Mon Sep 17 00:00:00 2001 From: singularitti Date: Sat, 21 Oct 2023 03:49:44 -0400 Subject: [PATCH 1/2] Reformat all files uing blackformatter --- qha/__init__.py | 15 +- qha/basic_io/input_maker.py | 129 ++++++--- qha/basic_io/out.py | 44 +-- qha/basic_io/read_input.py | 57 +++- qha/calculator.py | 207 ++++++++++---- qha/cli/__init__.py | 8 +- qha/cli/converter.py | 6 +- qha/cli/parser.py | 27 +- qha/cli/plotter.py | 100 ++++--- qha/cli/runner.py | 264 +++++++++++------ qha/fitting.py | 22 +- qha/grid_interpolation.py | 43 ++- .../different_phonon_dos.py | 103 +++++-- qha/multi_configurations/same_phonon_dos.py | 100 +++++-- qha/plotting.py | 220 ++++++++------- qha/settings.py | 42 +-- qha/single_configuration.py | 33 ++- qha/statmech.py | 48 ++-- qha/tests/test_different_phonon_dos.py | 265 ++++++++++++++---- qha/tests/test_overall_run.py | 42 ++- qha/tests/test_read_input.py | 27 +- qha/tests/test_samevdos_overall.py | 54 ++-- qha/tests/test_single_configuration.py | 2 +- qha/tests/test_unit_conversion.py | 10 +- qha/thermodynamics.py | 46 +-- qha/tools.py | 73 +++-- qha/type_aliases.py | 2 +- qha/unit_conversion.py | 118 ++++---- qha/v2p.py | 30 +- 29 files changed, 1417 insertions(+), 720 deletions(-) diff --git a/qha/__init__.py b/qha/__init__.py index cb5ab6f..3c2c73e 100644 --- a/qha/__init__.py +++ b/qha/__init__.py @@ -4,12 +4,11 @@ # In case of user who does not have ``pip`` above version 9.0.0 if sys.version_info < (3, 5): - raise EnvironmentError('Please use Python version higher than 3.5!') + raise EnvironmentError("Please use Python version higher than 3.5!") -__author__ = {'Tian Qin': 'qinxx197@umn.edu', - 'Qi Zhang': 'qz2280@columbia.edu'} -__copyright__ = 'Copyright (c) 2018, Renata group' -__credits__ = {'Renata M. M. Wentzcovitch': 'rmw2150@columbia.edu'} -__date__ = 'Feb 17, 2018' -__maintainer__ = 'Tian Qin, Qi Zhang' -__version__ = '1.0.22' +__author__ = {"Tian Qin": "qinxx197@umn.edu", "Qi Zhang": "qz2280@columbia.edu"} +__copyright__ = "Copyright (c) 2018, Renata group" +__credits__ = {"Renata M. M. Wentzcovitch": "rmw2150@columbia.edu"} +__date__ = "Feb 17, 2018" +__maintainer__ = "Tian Qin, Qi Zhang" +__version__ = "1.0.22" diff --git a/qha/basic_io/input_maker.py b/qha/basic_io/input_maker.py index fc3ecde..ecd5380 100644 --- a/qha/basic_io/input_maker.py +++ b/qha/basic_io/input_maker.py @@ -25,7 +25,7 @@ from yaml import Loader # ===================== What can be exported? ===================== -__all__ = ['FromQEOutput'] +__all__ = ["FromQEOutput"] class FromQEOutput: @@ -68,12 +68,12 @@ def read_file_list(self) -> None: """ Read all the files' names for frequency files given by Quantum ESPRESSO program ``matdyn.x``. """ - with open(self._inp_file_list, 'r') as f: + with open(self._inp_file_list, "r") as f: d = load(f, Loader=Loader) - self.formula_unit_number = int(d['formula_unit_number']) - self.comment: str = d['comment'] - self._frequency_files: List[str] = d['frequency_files'] + self.formula_unit_number = int(d["formula_unit_number"]) + self.comment: str = d["comment"] + self._frequency_files: List[str] = d["frequency_files"] def read_static(self) -> None: """ @@ -84,16 +84,22 @@ def read_static(self) -> None: volumes = [] energies = [] - with open(self._inp_static, 'r') as f: + with open(self._inp_static, "r") as f: print( - "Reading static data: emtpy lines or lines starting with '#' will be ignored!") + "Reading static data: emtpy lines or lines starting with '#' will be ignored!" + ) for line in f: - if not line.strip() or line.startswith('#'): # Ignore empty line or comment line + if not line.strip() or line.startswith( + "#" + ): # Ignore empty line or comment line continue - match = re.match(r"p\s*=\s*(-?\d*\.?\d*)\s*v\s*=\s*(-?\d*\.?\d*)\s*e\s*=\s*(-?\d*\.?\d*)", line, - flags=re.IGNORECASE) + match = re.match( + r"p\s*=\s*(-?\d*\.?\d*)\s*v\s*=\s*(-?\d*\.?\d*)\s*e\s*=\s*(-?\d*\.?\d*)", + line, + flags=re.IGNORECASE, + ) if match is None: continue @@ -113,29 +119,35 @@ def read_q_points(self) -> None: q_coordinates = [] q_weights = [] - with open(self._inp_q_points, 'r') as f: + with open(self._inp_q_points, "r") as f: print( - "Reading q-points file: emtpy lines or lines starting with '#' will be ignored!") + "Reading q-points file: emtpy lines or lines starting with '#' will be ignored!" + ) regex = re.compile( - r"\s*(-?\d*\.?\d*)\s*(-?\d*\.?\d*)\s*(-?\d*\.?\d*)\s*(-?\d*\.?\d*)") + r"\s*(-?\d*\.?\d*)\s*(-?\d*\.?\d*)\s*(-?\d*\.?\d*)\s*(-?\d*\.?\d*)" + ) for line in f: - if not line.strip() or line.startswith('#'): # Ignore empty line or comment line + if not line.strip() or line.startswith( + "#" + ): # Ignore empty line or comment line continue match = regex.match(line) if not regex.match(line): raise RuntimeError( - "Unknown line! Should be 3 coordinates and 1 weight, other lines should be commented with '#'.") + "Unknown line! Should be 3 coordinates and 1 weight, other lines should be commented with '#'." + ) else: g = match.groups() q_coordinates.append(g[0:3]) # TODO: Check possible bug, what if the regex match fails q_weights.append(g[3]) - self.q_coordinates = np.array(q_coordinates, - dtype=float) # TODO: Possible bug, ``np.array([])`` is regarded as ``False`` + self.q_coordinates = np.array( + q_coordinates, dtype=float + ) # TODO: Possible bug, ``np.array([])`` is regarded as ``False`` self.q_weights = np.array(q_weights, dtype=float) @staticmethod @@ -165,15 +177,17 @@ def read_frequency_file(inp: str) -> Tuple[Vector, Matrix]: if not line.strip(): continue - if 'nbnd' or 'nks' in line: + if "nbnd" or "nks" in line: match = regex.search(line) if not match: raise RuntimeError( - "The head line '{0}' is not complete! Here 'nbnd' and 'nks' are not found!".format(line)) + "The head line '{0}' is not complete! Here 'nbnd' and 'nks' are not found!".format( + line + ) + ) else: - bands_amount, q_points_amount = strings_to_integers( - match.groups()) + bands_amount, q_points_amount = strings_to_integers(match.groups()) break gen: Iterator[str] = text_stream.generator_starts_from(offset) @@ -192,7 +206,7 @@ def read_frequency_file(inp: str) -> Tuple[Vector, Matrix]: # Sometimes QE prints negative numbers that coalesce with the previous one, # so we need to split not only spaces but also minus signs # Regex from https://stackoverflow.com/a/30858977/3260253 - freqs = list(filter(lambda s: s, re.split(r'\s+|(? Tuple[Vector, Matrix]: if q_coordinates.shape[0] != q_points_amount: raise RuntimeError( "The number of q-points detected, {0}, is not the same as what specified in head line!".format( - q_coordinates.shape[0])) + q_coordinates.shape[0] + ) + ) if frequencies.shape != (q_points_amount, bands_amount): raise RuntimeError( "The frequencies array shape '{0}' is not the same as '{1}'!".format( - frequencies.shape, (q_points_amount, bands_amount))) + frequencies.shape, (q_points_amount, bands_amount) + ) + ) return q_coordinates, frequencies @@ -219,25 +237,30 @@ def read_frequency_files(self) -> None: """ frequencies_for_all_files = [] - if any(_ is None for _ in (self.q_coordinates, self.q_weights)): # If any of them is ``None`` + if any( + _ is None for _ in (self.q_coordinates, self.q_weights) + ): # If any of them is ``None`` self.read_q_points() # Fill these 2 properties for i in range(len(self._frequency_files)): q_coordinates, frequencies = self.read_frequency_file( - self._frequency_files[i]) + self._frequency_files[i] + ) # Here I use `allclose` rather than `array_equal` since they may have very little # differences even they are supposed to be the same, because of the digits QE gave. if not np.allclose(q_coordinates, self.q_coordinates): - warnings.warn("The q-points' coordinates are different from what specified in the q-point file!", - stacklevel=1) + warnings.warn( + "The q-points' coordinates are different from what specified in the q-point file!", + stacklevel=1, + ) frequencies_for_all_files.append(frequencies) # Shape: (# volumes, # q-points, # bands on each point) self.frequencies = np.array(frequencies_for_all_files) - def write_to_file(self, outfile='input') -> None: + def write_to_file(self, outfile="input") -> None: """ Write all data to a file *outfile*, which will be regarded as standard input file for ``qha``. @@ -246,31 +269,47 @@ def write_to_file(self, outfile='input') -> None: path = pathlib.Path(outfile) if path.is_file(): print( - "Old '{0}' file found, I will backup it before continue.".format(outfile)) - path.rename(outfile + '.backup') + "Old '{0}' file found, I will backup it before continue.".format( + outfile + ) + ) + path.rename(outfile + ".backup") - with open(outfile, 'w') as f: + with open(outfile, "w") as f: f.write("# {0}\n".format(self.comment)) - f.write('# The file contains frequencies and weights at the END!\n') + f.write("# The file contains frequencies and weights at the END!\n") f.write( - '# Number of volumes (nv), q-vectors (nq), normal mode (np), formula units(nm)\n') + "# Number of volumes (nv), q-vectors (nq), normal mode (np), formula units(nm)\n" + ) # TODO: Possible bug introduced in formatting - f.write("{0} {1} {2} {3}\n\n".format(len(self.volumes), len(self.q_weights), - self.frequencies.shape[-1], self.formula_unit_number)) + f.write( + "{0} {1} {2} {3}\n\n".format( + len(self.volumes), + len(self.q_weights), + self.frequencies.shape[-1], + self.formula_unit_number, + ) + ) for i in range(len(self.volumes)): - f.write("P={0:20.10f} V={1:20.10f} E={2:20.10f}\n".format(self.pressures[i], self.volumes[i], - self.static_energies[i])) + f.write( + "P={0:20.10f} V={1:20.10f} E={2:20.10f}\n".format( + self.pressures[i], self.volumes[i], self.static_energies[i] + ) + ) for j in range(len(self.q_weights)): - f.write("{0:12.8f} {1:12.8f} {2:12.8f}\n".format( - *self.q_coordinates[j])) + f.write( + "{0:12.8f} {1:12.8f} {2:12.8f}\n".format(*self.q_coordinates[j]) + ) for k in range(self.frequencies.shape[-1]): - f.write("{0:20.10f}\n".format( - self.frequencies[i, j, k])) + f.write("{0:20.10f}\n".format(self.frequencies[i, j, k])) - f.write('\nweight\n') + f.write("\nweight\n") for j in range(len(self.q_weights)): - f.write("{0:12.8f} {1:12.8f} {2:12.8f} {3:12.8f}\n".format( - *self.q_coordinates[j], self.q_weights[j])) + f.write( + "{0:12.8f} {1:12.8f} {2:12.8f} {3:12.8f}\n".format( + *self.q_coordinates[j], self.q_weights[j] + ) + ) diff --git a/qha/basic_io/out.py b/qha/basic_io/out.py index 95f9f6e..1323368 100644 --- a/qha/basic_io/out.py +++ b/qha/basic_io/out.py @@ -14,17 +14,17 @@ def save_to_output(fn_output, text): - with open(fn_output, 'a') as f: - f.write(text + '\n') + with open(fn_output, "a") as f: + f.write(text + "\n") def save_x_tp(df, t, desired_pressures_gpa, p_sample_gpa, outfile_name): # To fix the last 2 temperature points of Cp is not accurate, we added 4 temperature points, # Now to get rid of the added temperature points before saving calculated properties into files. df = pd.DataFrame(df, index=t, columns=desired_pressures_gpa).iloc[:-4, :] - df.columns.name = 'T(K)\P(GPa)' + df.columns.name = "T(K)\P(GPa)" sample = df.loc[:, df.columns.isin(p_sample_gpa)] - with open(outfile_name, 'w') as f: + with open(outfile_name, "w") as f: f.write(sample.to_string()) @@ -32,48 +32,60 @@ def save_x_pt(df, t, desired_pressures_gpa, t_sample, outfile_name): # To fix the last 2 temperature points of Cp is not accurate, we added 4 temperature points, # Now to get rid of the added temperature points before saving calculated properties into files. df = pd.DataFrame(df[:-4].T, index=desired_pressures_gpa, columns=t[:-4]) - df.columns.name = 'P(GPa)\T(K)' + df.columns.name = "P(GPa)\T(K)" sample = df.loc[:, df.columns.isin(t_sample)] - with open(outfile_name, 'w') as f: + with open(outfile_name, "w") as f: f.write(sample.to_string()) def save_x_vt(x, t, volume_grid, t_sample, outfile_name): df = pd.DataFrame(x.T, index=volume_grid, columns=t) - df.columns.name = 'V(A^3)\T(K)' + df.columns.name = "V(A^3)\T(K)" sample = df.loc[:, df.columns.isin(t_sample)] - with open(outfile_name, 'w') as f: + with open(outfile_name, "w") as f: f.write(sample.to_string()) def save_x_tv(x, t, volume_grid, t_sample, outfile_name): df = pd.DataFrame(x, index=t, columns=volume_grid).iloc[:-4, :] - df.columns.name = 'T(K)\V(A^3)' + df.columns.name = "T(K)\V(A^3)" sample = df.loc[df.index.isin(t_sample[:-4]), :] - with open(outfile_name, 'w') as f: + with open(outfile_name, "w") as f: f.write(sample.to_string()) def make_starting_string() -> str: - return textwrap.dedent("""\ + return textwrap.dedent( + """\ ============================================================ Current time: {0:%Y-%m-%d %H:%M:%S} - """.format(datetime.utcnow())) + """.format( + datetime.utcnow() + ) + ) def make_tp_info(min_temperature, max_temperature, min_pressure, max_pressure): - return textwrap.dedent("""\ + return textwrap.dedent( + """\ ------------------------------------------------------------ Desired T range: {0:6.2f} to {1:6.2f} K Desired P range: {2:6.2f} to {3:6.2f} GPa ------------------------------------------------------------ - """.format(min_temperature, max_temperature, min_pressure, max_pressure)) + """.format( + min_temperature, max_temperature, min_pressure, max_pressure + ) + ) def make_ending_string(time_elapsed) -> str: - return textwrap.dedent("""\ + return textwrap.dedent( + """\ ------------------------------------------------------------ Total elapsed time is: {0:8.2f} seconds Thanks for using QHA code, have a nice one :) ============================================================ - """.format(time_elapsed)) + """.format( + time_elapsed + ) + ) diff --git a/qha/basic_io/read_input.py b/qha/basic_io/read_input.py index b965a1c..d7cfefa 100755 --- a/qha/basic_io/read_input.py +++ b/qha/basic_io/read_input.py @@ -18,10 +18,12 @@ from qha.type_aliases import Vector, Array3D # ===================== What can be exported? ===================== -__all__ = ['read_input'] +__all__ = ["read_input"] -def read_input(inp: Union[str, pathlib.PurePath]) -> Tuple[int, Vector, Vector, Array3D, Vector]: +def read_input( + inp: Union[str, pathlib.PurePath] +) -> Tuple[int, Vector, Vector, Array3D, Vector]: """ Read the standard "input" file for ``qha``. @@ -47,24 +49,36 @@ def read_input(inp: Union[str, pathlib.PurePath]) -> Tuple[int, Vector, Vector, regex0 = re.compile("\s*(\d+)[\s,]*(\d+)[\s,]*(\d+)[\s,]*(\d+)") for line, offset in gen: - if not line.strip() or line.startswith('#'): + if not line.strip() or line.startswith("#"): continue match = regex0.search(line) if match is None: continue else: - volumes_amount, q_points_amount, modes_per_q_point_amount, formula_unit_number = strings_to_integers( - match.groups()) + ( + volumes_amount, + q_points_amount, + modes_per_q_point_amount, + formula_unit_number, + ) = strings_to_integers(match.groups()) break # If the metadata is not found, check the *inp*! - if any(_ is None for _ in (volumes_amount, q_points_amount, modes_per_q_point_amount)): - raise ValueError("At least one of the desired values 'nv', 'nq', 'np' is not found in file {0}!".format(inp)) + if any( + _ is None for _ in (volumes_amount, q_points_amount, modes_per_q_point_amount) + ): + raise ValueError( + "At least one of the desired values 'nv', 'nq', 'np' is not found in file {0}!".format( + inp + ) + ) # Generate containers for storing the following data. volumes = np.empty(volumes_amount, dtype=float) static_energies = np.empty(volumes_amount, dtype=float) - frequencies = np.empty((volumes_amount, q_points_amount, modes_per_q_point_amount), dtype=float) + frequencies = np.empty( + (volumes_amount, q_points_amount, modes_per_q_point_amount), dtype=float + ) q_weights = np.empty(q_points_amount, dtype=float) # We start a new iterator from where we stopped. @@ -74,16 +88,23 @@ def read_input(inp: Union[str, pathlib.PurePath]) -> Tuple[int, Vector, Vector, j = 0 # q-point index, note it is not count like `i`! # Now we start reading the energies, volumes, and frequencies. - regex1 = re.compile("P\s*=\s*-?\d*\.?\d*\s*V\s*=(\s*\d*\.?\d*)\s*E\s*=\s*(-?\d*\.?\d*)", re.IGNORECASE) + regex1 = re.compile( + "P\s*=\s*-?\d*\.?\d*\s*V\s*=(\s*\d*\.?\d*)\s*E\s*=\s*(-?\d*\.?\d*)", + re.IGNORECASE, + ) for line in gen: if not line.strip(): continue - if '=' in line: + if "=" in line: match = regex1.search(line) if match is None: - raise ValueError("Search of pattern {0} failed in line '{1}!".format(regex1.pattern, line)) + raise ValueError( + "Search of pattern {0} failed in line '{1}!".format( + regex1.pattern, line + ) + ) else: volumes[i], static_energies[i] = match.groups() i += 1 @@ -92,14 +113,16 @@ def read_input(inp: Union[str, pathlib.PurePath]) -> Tuple[int, Vector, Vector, sp = line.split() if len(sp) == 3: - for k in range(modes_per_q_point_amount): # Note `k` is the index of mode, like `j`, not count like `i`. + for k in range( + modes_per_q_point_amount + ): # Note `k` is the index of mode, like `j`, not count like `i`. line = next(gen) frequencies[i - 1, j, k] = line j += 1 continue - if 'weight' in line.lower(): + if "weight" in line.lower(): break # Now we start reading q-point-weights. We start the exact iterator from where we stopped. @@ -109,10 +132,14 @@ def read_input(inp: Union[str, pathlib.PurePath]) -> Tuple[int, Vector, Vector, if not line.strip(): continue - q_weights[j] = line.split()[-1] # This line has format: q_x q_y q_z weight, and only weight is taken + q_weights[j] = line.split()[ + -1 + ] # This line has format: q_x q_y q_z weight, and only weight is taken j += 1 if i != volumes_amount: - raise ValueError('The number of volumes detected is not equal to what specified in head! Check your file!') + raise ValueError( + "The number of volumes detected is not equal to what specified in head! Check your file!" + ) return formula_unit_number, volumes, static_energies, frequencies, q_weights diff --git a/qha/calculator.py b/qha/calculator.py index a00fe78..a81ab47 100644 --- a/qha/calculator.py +++ b/qha/calculator.py @@ -26,23 +26,46 @@ from qha.single_configuration import free_energy from qha.thermodynamics import * from qha.type_aliases import Vector -from qha.unit_conversion import gpa_to_ry_b3, ry_b3_to_gpa, b3_to_a3, ry_to_j_mol, ry_to_ev, ry_to_j +from qha.unit_conversion import ( + gpa_to_ry_b3, + ry_b3_to_gpa, + b3_to_a3, + ry_to_j_mol, + ry_to_ev, + ry_to_j, +) from qha.v2p import v2p # ===================== What can be exported? ===================== -__all__ = ['Calculator', 'SamePhDOSCalculator', 'DifferentPhDOSCalculator'] +__all__ = ["Calculator", "SamePhDOSCalculator", "DifferentPhDOSCalculator"] class Calculator: def __init__(self, user_settings: Dict[str, Any]): runtime_settings = dict() - allowed_keys = ('input', 'calculation', - 'thermodynamic_properties', 'static_only', 'energy_unit', - 'T_MIN', 'NT', 'DT', 'DT_SAMPLE', - 'P_MIN', 'NTV', 'DELTA_P', 'DELTA_P_SAMPLE', - 'volume_ratio', 'order', 'p_min_modifier', - 'T4FV', 'output_directory', 'high_verbosity', 'qha_output') + allowed_keys = ( + "input", + "calculation", + "thermodynamic_properties", + "static_only", + "energy_unit", + "T_MIN", + "NT", + "DT", + "DT_SAMPLE", + "P_MIN", + "NTV", + "DELTA_P", + "DELTA_P_SAMPLE", + "volume_ratio", + "order", + "p_min_modifier", + "T4FV", + "output_directory", + "high_verbosity", + "qha_output", + ) for key in allowed_keys: try: @@ -107,15 +130,20 @@ def v_ratio(self) -> Optional[float]: def read_input(self): try: - formula_unit_number, volumes, static_energies, frequencies, q_weights = read_input( - self.settings['input']) + ( + formula_unit_number, + volumes, + static_energies, + frequencies, + q_weights, + ) = read_input(self.settings["input"]) except KeyError: - raise KeyError( - "The 'input' option must be given in your settings!") + raise KeyError("The 'input' option must be given in your settings!") if not qha.tools.is_monotonic_decreasing(volumes): raise RuntimeError( - "Check the input file to make sure the volume decreases!") + "Check the input file to make sure the volume decreases!" + ) self._formula_unit_number: int = formula_unit_number self._volumes = volumes @@ -149,35 +177,49 @@ def temperature_array(self) -> Vector: """ # Normally, the last 2 temperature points in Cp are not accurate. # Here 4 more points are added for calculation, but they will be removed at the output files. - minimum_temperature = self.settings['T_MIN'] + minimum_temperature = self.settings["T_MIN"] if minimum_temperature < 0: raise ValueError("Minimum temperature should be no less than 0!") - return qha.tools.arange(minimum_temperature, self.settings['NT'] + 4, self.settings['DT']) + return qha.tools.arange( + minimum_temperature, self.settings["NT"] + 4, self.settings["DT"] + ) def refine_grid(self): d = self.settings try: - p_min, p_min_modifier, ntv, order = d['P_MIN'], d['p_min_modifier'], d['NTV'], d['order'] + p_min, p_min_modifier, ntv, order = ( + d["P_MIN"], + d["p_min_modifier"], + d["NTV"], + d["order"], + ) except KeyError: raise KeyError( - "All the 'P_MIN', 'p_min_modifier', 'NTV', 'order' options must be given in your settings!") + "All the 'P_MIN', 'p_min_modifier', 'NTV', 'order' options must be given in your settings!" + ) r = FinerGrid(p_min - p_min_modifier, ntv, order=order) - if 'volume_ratio' in d: - self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid(self.volumes, self.vib_ry, - ratio=d['volume_ratio']) + if "volume_ratio" in d: + self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid( + self.volumes, self.vib_ry, ratio=d["volume_ratio"] + ) else: self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid( - self.volumes, self.vib_ry) + self.volumes, self.vib_ry + ) @LazyProperty def vib_ry(self): # We grep all the arguments once since they are being invoked for thousands of times, and will be an overhead. - args = self.q_weights, self.static_energies, self.frequencies, self.settings[ - 'static_only'] + args = ( + self.q_weights, + self.static_energies, + self.frequencies, + self.settings["static_only"], + ) mat = np.empty((self.temperature_array.size, self.volumes.size)) for i, t in enumerate(self.temperature_array): @@ -187,16 +229,20 @@ def vib_ry(self): @LazyProperty def thermodynamic_potentials(self) -> Dict[str, Any]: - return thermodynamic_potentials(self.temperature_array, self.finer_volumes_bohr3, self.f_tv_ry, self.p_tv_au) + return thermodynamic_potentials( + self.temperature_array, self.finer_volumes_bohr3, self.f_tv_ry, self.p_tv_au + ) @LazyProperty def temperature_sample_array(self): - return self.temperature_array[0::int(self.settings['DT_SAMPLE'] / self.settings['DT'])] + return self.temperature_array[ + 0 :: int(self.settings["DT_SAMPLE"] / self.settings["DT"]) + ] @LazyProperty def desired_pressures_gpa(self): d = self.settings - return qha.tools.arange(d['P_MIN'], d['NTV'], d['DELTA_P']) + return qha.tools.arange(d["P_MIN"], d["NTV"], d["DELTA_P"]) @LazyProperty def desired_pressures(self): @@ -204,7 +250,9 @@ def desired_pressures(self): @LazyProperty def pressure_sample_array(self): - return self.desired_pressures_gpa[0::int(self.settings['DELTA_P_SAMPLE'] / self.settings['DELTA_P'])] + return self.desired_pressures_gpa[ + 0 :: int(self.settings["DELTA_P_SAMPLE"] / self.settings["DELTA_P"]) + ] @LazyProperty def finer_volumes_ang3(self): @@ -221,24 +269,38 @@ def volumes_ang3(self): def desired_pressure_status(self) -> None: d = self.settings - if d['high_verbosity']: - save_to_output(d['qha_output'], "The pressure range can be dealt with: [{0:6.2f} to {1:6.2f}] GPa".format( - self.p_tv_gpa[:, 0].max(), self.p_tv_gpa[:, -1].min())) + if d["high_verbosity"]: + save_to_output( + d["qha_output"], + "The pressure range can be dealt with: [{0:6.2f} to {1:6.2f}] GPa".format( + self.p_tv_gpa[:, 0].max(), self.p_tv_gpa[:, -1].min() + ), + ) if self.p_tv_gpa[:, -1].min() < self.desired_pressures_gpa.max(): ntv_max = int( - (self.p_tv_gpa[:, -1].min() - self.desired_pressures_gpa.min()) / d['DELTA_P']) - - save_to_output(d['qha_output'], textwrap.dedent("""\ + (self.p_tv_gpa[:, -1].min() - self.desired_pressures_gpa.min()) + / d["DELTA_P"] + ) + + save_to_output( + d["qha_output"], + textwrap.dedent( + """\ !!!ATTENTION!!! DESIRED PRESSURE is too high (NTV is too large)! QHA results might not be right! Please reduce the NTV accordingly, for example, try to set NTV < {:4d}. - """.format(ntv_max))) + """.format( + ntv_max + ) + ), + ) raise ValueError( - "DESIRED PRESSURE is too high (NTV is too large), qha results might not be right!") + "DESIRED PRESSURE is too high (NTV is too large), qha results might not be right!" + ) @LazyProperty def p_tv_au(self): @@ -266,7 +328,7 @@ def f_tp_ev(self): @LazyProperty def u_tv_ry(self): - return self.thermodynamic_potentials['U'] + return self.thermodynamic_potentials["U"] @LazyProperty def u_tp_ry(self): @@ -278,7 +340,7 @@ def u_tp_ev(self): @LazyProperty def h_tv_ry(self): - return self.thermodynamic_potentials['H'] + return self.thermodynamic_potentials["H"] @LazyProperty def h_tp_ry(self): @@ -290,7 +352,7 @@ def h_tp_ev(self): @LazyProperty def g_tv_ry(self): - return self.thermodynamic_potentials['G'] + return self.thermodynamic_potentials["G"] @LazyProperty def g_tp_ry(self): @@ -342,11 +404,15 @@ def cv_tp_jmolk(self): @LazyProperty def gamma_tp(self): - return gruneisen_parameter(self.v_tp_bohr3, self.bt_tp_au, self.alpha_tp, self.cv_tp_au) + return gruneisen_parameter( + self.v_tp_bohr3, self.bt_tp_au, self.alpha_tp, self.cv_tp_au + ) @LazyProperty def bs_tp_au(self): - return adiabatic_bulk_modulus(self.bt_tp_au, self.alpha_tp, self.gamma_tp, self.temperature_array) + return adiabatic_bulk_modulus( + self.bt_tp_au, self.alpha_tp, self.gamma_tp, self.temperature_array + ) @LazyProperty def bs_tp_gpa(self): @@ -354,11 +420,15 @@ def bs_tp_gpa(self): @LazyProperty def cp_tp_au(self): - return isobaric_heat_capacity(self.cv_tp_au, self.alpha_tp, self.gamma_tp, self.temperature_array) + return isobaric_heat_capacity( + self.cv_tp_au, self.alpha_tp, self.gamma_tp, self.temperature_array + ) @LazyProperty def cp_tp_jmolk(self): - return isobaric_heat_capacity(self.cv_tp_jmolk, self.alpha_tp, self.gamma_tp, self.temperature_array) + return isobaric_heat_capacity( + self.cv_tp_jmolk, self.alpha_tp, self.gamma_tp, self.temperature_array + ) class DifferentPhDOSCalculator(Calculator): @@ -377,8 +447,8 @@ def volumes(self): return self._volumes[0] def read_input(self): - self._degeneracies = tuple(self.settings['input'].values()) - input_data_files = tuple(self.settings['input'].keys()) + self._degeneracies = tuple(self.settings["input"].values()) + input_data_files = tuple(self.settings["input"].keys()) formula_unit_numbers = [] volumes = [] @@ -387,15 +457,23 @@ def read_input(self): q_weights = [] for inp in input_data_files: - nm_tmp, volumes_tmp, static_energies_tmp, freq_tmp, weights_tmp = read_input( - inp) + ( + nm_tmp, + volumes_tmp, + static_energies_tmp, + freq_tmp, + weights_tmp, + ) = read_input(inp) if not qha.tools.is_monotonic_decreasing(volumes_tmp): # TODO: Clean this sentence save_to_output( - self.settings['qha_output'], "Check the input file to make sure the volume decreases") + self.settings["qha_output"], + "Check the input file to make sure the volume decreases", + ) raise ValueError( - "Check the input file to make sure the volumes are monotonic decreasing!") + "Check the input file to make sure the volumes are monotonic decreasing!" + ) formula_unit_numbers.append(nm_tmp) volumes.append(volumes_tmp) @@ -411,11 +489,11 @@ def read_input(self): if not len(set(formula_unit_numbers)) == 1: raise RuntimeError( - "All the formula unit number in all inputs should be the same!") + "All the formula unit number in all inputs should be the same!" + ) if len(volumes.shape) == 1: - raise RuntimeError( - "All configurations should have same number of volumes!") + raise RuntimeError("All configurations should have same number of volumes!") # Choose any of them since they are all the same self._formula_unit_number = formula_unit_numbers[0] @@ -427,13 +505,20 @@ def read_input(self): @LazyProperty def vib_ry(self): # We grep all the arguments once since they are being invoked for thousands of times, and will be an overhead. - args = self.degeneracies, self.q_weights, self.static_energies, self._volumes, self.frequencies, \ - self.settings['static_only'] + args = ( + self.degeneracies, + self.q_weights, + self.static_energies, + self._volumes, + self.frequencies, + self.settings["static_only"], + ) mat = np.empty((self.temperature_array.size, self._volumes.shape[1])) for i, t in enumerate(self.temperature_array): mat[i] = different_phonon_dos.PartitionFunction( - t, *(arg for arg in args)).get_free_energies() + t, *(arg for arg in args) + ).get_free_energies() return mat @@ -444,11 +529,19 @@ def __init__(self, user_settings): @LazyProperty def vib_ry(self): - args = self.degeneracies, self.q_weights[0], self.static_energies, self._volumes, self.frequencies[0], \ - self.settings['static_only'], self.settings['order'] + args = ( + self.degeneracies, + self.q_weights[0], + self.static_energies, + self._volumes, + self.frequencies[0], + self.settings["static_only"], + self.settings["order"], + ) mat = np.empty((self.temperature_array.size, self._volumes.shape[1])) for i, t in enumerate(self.temperature_array): mat[i] = same_phonon_dos.FreeEnergy( - t, *(arg for arg in args)).get_free_energies() + t, *(arg for arg in args) + ).get_free_energies() return mat diff --git a/qha/cli/__init__.py b/qha/cli/__init__.py index dc8cd7a..790282b 100644 --- a/qha/cli/__init__.py +++ b/qha/cli/__init__.py @@ -16,13 +16,13 @@ def main(): parser = QHAArgumentParser() qha_converter = QHAConverter() - parser.register_handler('convert', qha_converter, 'conv') + parser.register_handler("convert", qha_converter, "conv") qha_runner = QHARunner() - parser.register_handler('run', qha_runner) + parser.register_handler("run", qha_runner) qha_plotter = QHAPlotter() - parser.register_handler('plot', qha_plotter) + parser.register_handler("plot", qha_plotter) parser.load_plugins() @@ -30,5 +30,5 @@ def main(): parser.invoke_handler(namespace) -if __name__ == '__main__': +if __name__ == "__main__": main() diff --git a/qha/cli/converter.py b/qha/cli/converter.py index f2436fd..ca7a035 100644 --- a/qha/cli/converter.py +++ b/qha/cli/converter.py @@ -16,9 +16,9 @@ def __init__(self): def init_parser(self, parser): super().init_parser(parser) - parser.add_argument('inp_file_list', type=str) - parser.add_argument('inp_static', type=str) - parser.add_argument('inp_q_points', type=str) + parser.add_argument("inp_file_list", type=str) + parser.add_argument("inp_static", type=str) + parser.add_argument("inp_q_points", type=str) def run(self, namespace): inp_file_list = namespace.inp_file_list diff --git a/qha/cli/parser.py b/qha/cli/parser.py index f93b48b..abc1c92 100644 --- a/qha/cli/parser.py +++ b/qha/cli/parser.py @@ -31,9 +31,7 @@ class QHAArgumentParser: def __init__(self): self.parser = argparse.ArgumentParser() self.init_parser() - self.subparsers = self.parser.add_subparsers( - dest='command' - ) + self.subparsers = self.parser.add_subparsers(dest="command") self.handlers = dict() def init_parser(self): @@ -45,10 +43,11 @@ def init_parser(self): 2. short version: ``-h``, long version: ``--help``, the help document of ``qha``. """ self.parser.add_argument( - '-V', '--version', - action='version', + "-V", + "--version", + action="version", version="current qha version: {0}".format(__version__), - help='The current release version of ``qha``.' + help="The current release version of ``qha``.", ) def register_handler(self, command: str, handler: QHACommandHandler, *aliases): @@ -70,7 +69,9 @@ def register_handler(self, command: str, handler: QHACommandHandler, *aliases): raise TypeError("Argument *command* should be a string!") if not isinstance(handler, QHACommandHandler): - raise TypeError("Argument *handler* should be a ``QHACommandHandler`` instance!") + raise TypeError( + "Argument *handler* should be a ``QHACommandHandler`` instance!" + ) new_parser = self.subparsers.add_parser(command, aliases=list(aliases)) self.handlers.update( @@ -79,7 +80,7 @@ def register_handler(self, command: str, handler: QHACommandHandler, *aliases): SubCommandResolvers( handler=handler, parser=new_parser, - ) + ), ) ) handler.init_parser(new_parser) @@ -104,9 +105,11 @@ def invoke_handler(self, namespace): :param namespace: The namespace returned by ``parse_args`` method. """ try: - command: str = getattr(namespace, 'command') + command: str = getattr(namespace, "command") except AttributeError: - raise AttributeError("Argument *namespace* does not have an ``command`` attribute!") + raise AttributeError( + "Argument *namespace* does not have an ``command`` attribute!" + ) try: handler: QHACommandHandler = self.handlers[command].handler @@ -118,7 +121,7 @@ def load_plugins(self): """ Load plugins. Leave for the future plugin system. """ - for plugin in pkg_resources.iter_entry_points(group='qha.plugins'): + for plugin in pkg_resources.iter_entry_points(group="qha.plugins"): klass = plugin.load() - aliases = klass.aliases if 'aliases' in dir(klass) else None + aliases = klass.aliases if "aliases" in dir(klass) else None self.register_handler(plugin.name, klass(), *aliases) diff --git a/qha/cli/plotter.py b/qha/cli/plotter.py index 8052448..e710e93 100644 --- a/qha/cli/plotter.py +++ b/qha/cli/plotter.py @@ -22,74 +22,90 @@ def __init__(self): def init_parser(self, parser): super().init_parser(parser) - parser.add_argument('settings', type=str) - parser.add_argument('--outdir', type=str, help='output directory') + parser.add_argument("settings", type=str) + parser.add_argument("--outdir", type=str, help="output directory") def run(self, namespace): user_settings = {} # save necessary info for plotting later file_settings = namespace.settings settings = from_yaml(file_settings) - for key in ('energy_unit', 'NT', 'DT', 'DT_SAMPLE', 'P_MIN', 'NTV', 'DELTA_P', 'DELTA_P_SAMPLE', - 'thermodynamic_properties', 'T4FV', 'output_directory'): + for key in ( + "energy_unit", + "NT", + "DT", + "DT_SAMPLE", + "P_MIN", + "NTV", + "DELTA_P", + "DELTA_P_SAMPLE", + "thermodynamic_properties", + "T4FV", + "output_directory", + ): try: user_settings.update({key: settings[key]}) except KeyError: raise KeyError("Key '{0}' is not set in your settings!") - if not os.path.exists(user_settings['output_directory']): - raise FileNotFoundError("There is no results folder, please run: `qha-run` first! ") + if not os.path.exists(user_settings["output_directory"]): + raise FileNotFoundError( + "There is no results folder, please run: `qha-run` first! " + ) plotter = Plotter(user_settings) - desired_pressures_gpa = qha.tools.arange(user_settings['P_MIN'], user_settings['NTV'], user_settings['DELTA_P']) - user_settings.update({'DESIRED_PRESSURES_GPa': desired_pressures_gpa}) - - results_folder = pathlib.Path(user_settings['output_directory']) - - calculation_option = {'F': 'f_tp', - 'G': 'g_tp', - 'H': 'h_tp', - 'U': 'u_tp', - 'V': 'v_tp', - 'Cv': 'cv_tp_jmolk', - 'Cp': 'cp_tp_jmolk', - 'Bt': 'bt_tp_gpa', - 'Btp': 'btp_tp', - 'Bs': 'bs_tp_gpa', - 'alpha': 'alpha_tp', - 'gamma': 'gamma_tp', - } - - file_ftv_fitted = results_folder / 'f_tv_fitted_ev_ang3.txt' - user_settings.update({'f_tv_fitted': file_ftv_fitted}) - - file_ftv_non_fitted = results_folder / 'f_tv_nonfitted_ev_ang3.txt' - user_settings.update({'f_tv_non_fitted': file_ftv_non_fitted}) - - file_ptv_gpa = results_folder / 'p_tv_gpa.txt' - user_settings.update({'p_tv_gpa': file_ptv_gpa}) + desired_pressures_gpa = qha.tools.arange( + user_settings["P_MIN"], user_settings["NTV"], user_settings["DELTA_P"] + ) + user_settings.update({"DESIRED_PRESSURES_GPa": desired_pressures_gpa}) + + results_folder = pathlib.Path(user_settings["output_directory"]) + + calculation_option = { + "F": "f_tp", + "G": "g_tp", + "H": "h_tp", + "U": "u_tp", + "V": "v_tp", + "Cv": "cv_tp_jmolk", + "Cp": "cp_tp_jmolk", + "Bt": "bt_tp_gpa", + "Btp": "btp_tp", + "Bs": "bs_tp_gpa", + "alpha": "alpha_tp", + "gamma": "gamma_tp", + } + + file_ftv_fitted = results_folder / "f_tv_fitted_ev_ang3.txt" + user_settings.update({"f_tv_fitted": file_ftv_fitted}) + + file_ftv_non_fitted = results_folder / "f_tv_nonfitted_ev_ang3.txt" + user_settings.update({"f_tv_non_fitted": file_ftv_non_fitted}) + + file_ptv_gpa = results_folder / "p_tv_gpa.txt" + user_settings.update({"p_tv_gpa": file_ptv_gpa}) plotter.fv_pv() - for idx in user_settings['thermodynamic_properties']: - if idx in ['F', 'G', 'H', 'U']: - attr_name = calculation_option[idx] + '_' + user_settings['energy_unit'] - file_name = attr_name + '.txt' + for idx in user_settings["thermodynamic_properties"]: + if idx in ["F", "G", "H", "U"]: + attr_name = calculation_option[idx] + "_" + user_settings["energy_unit"] + file_name = attr_name + ".txt" file_dir = results_folder / file_name user_settings.update({idx: file_dir}) plotter.plot2file(idx) - if idx == 'V': - v_ang3 = calculation_option[idx] + '_' + 'ang3' - file_name_ang3 = v_ang3 + '.txt' + if idx == "V": + v_ang3 = calculation_option[idx] + "_" + "ang3" + file_name_ang3 = v_ang3 + ".txt" file_dir_ang3 = results_folder / file_name_ang3 user_settings.update({idx: file_dir_ang3}) plotter.plot2file(idx) - if idx in ['Cv', 'Cp', 'Bt', 'Btp', 'Bs', 'alpha', 'gamma']: + if idx in ["Cv", "Cp", "Bt", "Btp", "Bs", "alpha", "gamma"]: attr_name = calculation_option[idx] - file_name = attr_name + '.txt' + file_name = attr_name + ".txt" file_dir = results_folder / file_name user_settings.update({idx: file_dir}) plotter.plot2file(idx) diff --git a/qha/cli/runner.py b/qha/cli/runner.py index d59a036..1d9d4cd 100644 --- a/qha/cli/runner.py +++ b/qha/cli/runner.py @@ -12,7 +12,14 @@ import time from qha.calculator import Calculator, SamePhDOSCalculator, DifferentPhDOSCalculator -from qha.basic_io.out import save_x_tp, save_x_tv, save_to_output, make_starting_string, make_tp_info, make_ending_string +from qha.basic_io.out import ( + save_x_tp, + save_x_tv, + save_to_output, + make_starting_string, + make_tp_info, + make_ending_string, +) from qha.settings import from_yaml from .handler import QHACommandHandler @@ -26,7 +33,7 @@ def __init__(self): def init_parser(self, parser): super().init_parser(parser) - parser.add_argument('settings', type=str) + parser.add_argument("settings", type=str) def run(self, namespace): start_time_total = time.time() @@ -35,54 +42,85 @@ def run(self, namespace): file_settings = namespace.settings settings = from_yaml(file_settings) - for key in ('input', 'calculation', - 'thermodynamic_properties', 'static_only', 'energy_unit', - 'T_MIN', 'NT', 'DT', 'DT_SAMPLE', - 'P_MIN', 'NTV', 'DELTA_P', 'DELTA_P_SAMPLE', - 'volume_ratio', 'order', 'p_min_modifier', - 'T4FV', 'output_directory', 'high_verbosity'): + for key in ( + "input", + "calculation", + "thermodynamic_properties", + "static_only", + "energy_unit", + "T_MIN", + "NT", + "DT", + "DT_SAMPLE", + "P_MIN", + "NTV", + "DELTA_P", + "DELTA_P_SAMPLE", + "volume_ratio", + "order", + "p_min_modifier", + "T4FV", + "output_directory", + "high_verbosity", + ): try: user_settings.update({key: settings[key]}) except KeyError: continue - if not os.path.exists(user_settings['output_directory']): - os.makedirs(user_settings['output_directory']) + if not os.path.exists(user_settings["output_directory"]): + os.makedirs(user_settings["output_directory"]) - user_settings.update({'qha_output': os.path.join( - user_settings['output_directory'], 'output.txt')}) + user_settings.update( + { + "qha_output": os.path.join( + user_settings["output_directory"], "output.txt" + ) + } + ) try: - os.remove(user_settings['qha_output']) + os.remove(user_settings["qha_output"]) except OSError: pass - save_to_output(user_settings['qha_output'], make_starting_string()) + save_to_output(user_settings["qha_output"], make_starting_string()) - calculation_type = user_settings['calculation'].lower() - if calculation_type == 'single': + calculation_type = user_settings["calculation"].lower() + if calculation_type == "single": calc = Calculator(user_settings) print("You have single-configuration calculation assumed.") - elif calculation_type == 'same phonon dos': + elif calculation_type == "same phonon dos": calc = SamePhDOSCalculator(user_settings) print( - "You have multi-configuration calculation with the same phonon DOS assumed.") - elif calculation_type == 'different phonon dos': + "You have multi-configuration calculation with the same phonon DOS assumed." + ) + elif calculation_type == "different phonon dos": calc = DifferentPhDOSCalculator(user_settings) print( - "You have multi-configuration calculation with different phonon DOS assumed.") + "You have multi-configuration calculation with different phonon DOS assumed." + ) else: - raise ValueError("The 'calculation' in your settings in not recognized! It must be one of:" - "'single', 'same phonon dos', 'different phonon dos'!") - - save_to_output(user_settings['qha_output'], - make_tp_info(calc.temperature_array[0], calc.temperature_array[-1 - 4], - calc.desired_pressures_gpa[0], - calc.desired_pressures_gpa[-1])) + raise ValueError( + "The 'calculation' in your settings in not recognized! It must be one of:" + "'single', 'same phonon dos', 'different phonon dos'!" + ) + + save_to_output( + user_settings["qha_output"], + make_tp_info( + calc.temperature_array[0], + calc.temperature_array[-1 - 4], + calc.desired_pressures_gpa[0], + calc.desired_pressures_gpa[-1], + ), + ) calc.read_input() - print("Caution: If negative frequencies found, they are currently treated as 0!") + print( + "Caution: If negative frequencies found, they are currently treated as 0!" + ) tmp = calc.where_negative_frequencies # Don't delete this parenthesis! if tmp is not None and not (tmp.T[-1].max() <= 2): @@ -90,17 +128,26 @@ def run(self, namespace): for indices in tmp: print( "Found negative frequency in {0}th configuration {1}th volume {2}th q-point {3}th band".format( - *tuple(indices + 1))) + *tuple(indices + 1) + ) + ) elif calc.frequencies.ndim == 3: # Single configuration for indices in tmp: print( - "Found negative frequency in {0}th volume {1}th q-point {2}th band".format(*tuple(indices + 1))) + "Found negative frequency in {0}th volume {1}th q-point {2}th band".format( + *tuple(indices + 1) + ) + ) calc.refine_grid() - if user_settings['high_verbosity']: - save_to_output(user_settings['qha_output'], - 'The volume range used in this calculation expanded x {0:6.4f}'.format(calc.v_ratio)) + if user_settings["high_verbosity"]: + save_to_output( + user_settings["qha_output"], + "The volume range used in this calculation expanded x {0:6.4f}".format( + calc.v_ratio + ), + ) calc.desired_pressure_status() @@ -109,68 +156,107 @@ def run(self, namespace): temperature_sample = calc.temperature_sample_array p_sample_gpa = calc.pressure_sample_array - results_folder = pathlib.Path(user_settings['output_directory']) - - calculation_option = {'F': 'f_tp', - 'G': 'g_tp', - 'H': 'h_tp', - 'U': 'u_tp', - 'V': 'v_tp', - 'Cv': 'cv_tp_jmolk', - 'Cp': 'cp_tp_jmolk', - 'Bt': 'bt_tp_gpa', - 'Btp': 'btp_tp', - 'Bs': 'bs_tp_gpa', - 'alpha': 'alpha_tp', - 'gamma': 'gamma_tp', - } - - file_ftv_fitted = results_folder / 'f_tv_fitted_ev_ang3.txt' - save_x_tv(calc.f_tv_ev, temperature_array, - calc.finer_volumes_ang3, temperature_sample, file_ftv_fitted) - - file_ftv_non_fitted = results_folder / 'f_tv_nonfitted_ev_ang3.txt' - save_x_tv(calc.vib_ev, temperature_array, calc.volumes_ang3, - temperature_sample, file_ftv_non_fitted) - - file_ptv_gpa = results_folder / 'p_tv_gpa.txt' - save_x_tv(calc.p_tv_gpa, temperature_array, - calc.finer_volumes_ang3, temperature_sample, file_ptv_gpa) - - file_stv_j = results_folder / 's_tv_j.txt' - save_x_tv(calc.s_tv_j, temperature_array, - calc.finer_volumes_ang3, temperature_sample, file_stv_j) - - for idx in calc.settings['thermodynamic_properties']: - if idx in ['F', 'G', 'H', 'U']: - attr_name = calculation_option[idx] + \ - '_' + calc.settings['energy_unit'] - file_name = attr_name + '.txt' + results_folder = pathlib.Path(user_settings["output_directory"]) + + calculation_option = { + "F": "f_tp", + "G": "g_tp", + "H": "h_tp", + "U": "u_tp", + "V": "v_tp", + "Cv": "cv_tp_jmolk", + "Cp": "cp_tp_jmolk", + "Bt": "bt_tp_gpa", + "Btp": "btp_tp", + "Bs": "bs_tp_gpa", + "alpha": "alpha_tp", + "gamma": "gamma_tp", + } + + file_ftv_fitted = results_folder / "f_tv_fitted_ev_ang3.txt" + save_x_tv( + calc.f_tv_ev, + temperature_array, + calc.finer_volumes_ang3, + temperature_sample, + file_ftv_fitted, + ) + + file_ftv_non_fitted = results_folder / "f_tv_nonfitted_ev_ang3.txt" + save_x_tv( + calc.vib_ev, + temperature_array, + calc.volumes_ang3, + temperature_sample, + file_ftv_non_fitted, + ) + + file_ptv_gpa = results_folder / "p_tv_gpa.txt" + save_x_tv( + calc.p_tv_gpa, + temperature_array, + calc.finer_volumes_ang3, + temperature_sample, + file_ptv_gpa, + ) + + file_stv_j = results_folder / "s_tv_j.txt" + save_x_tv( + calc.s_tv_j, + temperature_array, + calc.finer_volumes_ang3, + temperature_sample, + file_stv_j, + ) + + for idx in calc.settings["thermodynamic_properties"]: + if idx in ["F", "G", "H", "U"]: + attr_name = calculation_option[idx] + "_" + calc.settings["energy_unit"] + file_name = attr_name + ".txt" file_dir = results_folder / file_name - save_x_tp(getattr(calc, attr_name), temperature_array, - desired_pressures_gpa, p_sample_gpa, file_dir) - - if idx == 'V': - v_bohr3 = calculation_option[idx] + '_' + 'bohr3' - file_name_bohr3 = v_bohr3 + '.txt' + save_x_tp( + getattr(calc, attr_name), + temperature_array, + desired_pressures_gpa, + p_sample_gpa, + file_dir, + ) + + if idx == "V": + v_bohr3 = calculation_option[idx] + "_" + "bohr3" + file_name_bohr3 = v_bohr3 + ".txt" file_dir_au = results_folder / file_name_bohr3 - v_ang3 = calculation_option[idx] + '_' + 'ang3' - file_name_ang3 = v_ang3 + '.txt' + v_ang3 = calculation_option[idx] + "_" + "ang3" + file_name_ang3 = v_ang3 + ".txt" file_dir_ang3 = results_folder / file_name_ang3 - save_x_tp(getattr(calc, v_bohr3), temperature_array, - desired_pressures_gpa, p_sample_gpa, file_dir_au) - save_x_tp(getattr(calc, v_ang3), temperature_array, - desired_pressures_gpa, p_sample_gpa, file_dir_ang3) - - if idx in ['Cv', 'Cp', 'Bt', 'Btp', 'Bs', 'alpha', 'gamma']: + save_x_tp( + getattr(calc, v_bohr3), + temperature_array, + desired_pressures_gpa, + p_sample_gpa, + file_dir_au, + ) + save_x_tp( + getattr(calc, v_ang3), + temperature_array, + desired_pressures_gpa, + p_sample_gpa, + file_dir_ang3, + ) + + if idx in ["Cv", "Cp", "Bt", "Btp", "Bs", "alpha", "gamma"]: attr_name = calculation_option[idx] - file_name = attr_name + '.txt' + file_name = attr_name + ".txt" file_dir = results_folder / file_name - save_x_tp(getattr(calc, attr_name), temperature_array, - desired_pressures_gpa, p_sample_gpa, file_dir) + save_x_tp( + getattr(calc, attr_name), + temperature_array, + desired_pressures_gpa, + p_sample_gpa, + file_dir, + ) end_time_total = time.time() time_elapsed = end_time_total - start_time_total - save_to_output(user_settings['qha_output'], - make_ending_string(time_elapsed)) + save_to_output(user_settings["qha_output"], make_ending_string(time_elapsed)) diff --git a/qha/fitting.py b/qha/fitting.py index 8c04300..5d351a9 100644 --- a/qha/fitting.py +++ b/qha/fitting.py @@ -14,7 +14,7 @@ from qha.type_aliases import Matrix, Vector # ===================== What can be exported? ===================== -__all__ = ['polynomial_least_square_fitting', 'apply_finite_strain_fitting'] +__all__ = ["polynomial_least_square_fitting", "apply_finite_strain_fitting"] def polynomial_least_square_fitting(xs, ys, new_xs, order: Optional[int] = 3): @@ -30,14 +30,20 @@ def polynomial_least_square_fitting(xs, ys, new_xs, order: Optional[int] = 3): :return: A tuple, the polynomial-fitting coefficients and the new vector of y-coordinates. """ order += 1 # The definition of order in ``numpy.vander`` is different from the order in finite strain by one. - xx = np.vander(xs, order, increasing=True) # This will make a Vandermonde matrix that will be used in EoS fitting. + xx = np.vander( + xs, order, increasing=True + ) # This will make a Vandermonde matrix that will be used in EoS fitting. a, _, _, _ = np.linalg.lstsq(xx, ys, rcond=None) new_y = np.vander(new_xs, order, increasing=True) @ a return a, new_y -def apply_finite_strain_fitting(strains_sparse: Vector, free_energies: Matrix, strains_dense: Vector, - order: Optional[int] = 3): +def apply_finite_strain_fitting( + strains_sparse: Vector, + free_energies: Matrix, + strains_dense: Vector, + order: Optional[int] = 3, +): """ Calculate the free energies :math:`F(T, V)` for some strains (*strains_dense*), with the free energies (*free_energies*) on some other strains (*strains_sparse*) known already. @@ -53,9 +59,13 @@ def apply_finite_strain_fitting(strains_sparse: Vector, free_energies: Matrix, s """ temperature_amount, _ = free_energies.shape dense_volume_amount = len(strains_dense) - f_v_t = np.empty((temperature_amount, dense_volume_amount)) # Initialize the F(T,V) array + f_v_t = np.empty( + (temperature_amount, dense_volume_amount) + ) # Initialize the F(T,V) array for i in range(temperature_amount): - _, f_i = polynomial_least_square_fitting(strains_sparse, free_energies[i], strains_dense, order) + _, f_i = polynomial_least_square_fitting( + strains_sparse, free_energies[i], strains_dense, order + ) f_v_t[i] = f_i return f_v_t diff --git a/qha/grid_interpolation.py b/qha/grid_interpolation.py index 7e77cb5..70787c7 100644 --- a/qha/grid_interpolation.py +++ b/qha/grid_interpolation.py @@ -19,10 +19,10 @@ # ===================== What can be exported? ===================== __all__ = [ - 'calculate_eulerian_strain', - 'from_eulerian_strain', - 'VolumeExpander', - 'FinerGrid' + "calculate_eulerian_strain", + "from_eulerian_strain", + "VolumeExpander", + "FinerGrid", ] @@ -123,7 +123,9 @@ def out_volumes_num(self, value: int): if not isinstance(value, int): raise TypeError("The argument *out_volumes_num* must be an integer!") if value <= 0: - raise ValueError("The argument *out_volumes_num* must be an integer larger than 0!") + raise ValueError( + "The argument *out_volumes_num* must be an integer larger than 0!" + ) self._out_volumes_num = value @property @@ -156,7 +158,9 @@ class is created, unless user is clear what is being done. # r = v_upper / v_max = v_min / v_lower v_lower, v_upper = v_min / self._ratio, v_max * self._ratio # The *v_max* is a reference value here. - s_upper, s_lower = calculate_eulerian_strain(v_max, v_lower), calculate_eulerian_strain(v_max, v_upper) + s_upper, s_lower = calculate_eulerian_strain( + v_max, v_lower + ), calculate_eulerian_strain(v_max, v_upper) self._strains = np.linspace(s_lower, s_upper, self._out_volumes_num) self._out_volumes = from_eulerian_strain(v_max, self._strains) @@ -171,7 +175,9 @@ class FinerGrid: :param order: The order of the Birch--Murnaghan finite-strain equation of state fitting. """ - def __init__(self, desired_p_min: float, dense_volumes_amount: int, order: Optional[int] = 3): + def __init__( + self, desired_p_min: float, dense_volumes_amount: int, order: Optional[int] = 3 + ): self.desired_p_min = float(desired_p_min) self.dense_volumes_amount = int(dense_volumes_amount) self.option = int(order) @@ -181,7 +187,9 @@ def __init__(self, desired_p_min: float, dense_volumes_amount: int, order: Optio def ratio(self) -> float: return self._ratio - def approach_to_best_ratio(self, volumes: Vector, free_energies: Vector, initial_ratio: float) -> float: + def approach_to_best_ratio( + self, volumes: Vector, free_energies: Vector, initial_ratio: float + ) -> float: """ Trying to find the best volume grids based on an a very large volume grids. @@ -194,7 +202,9 @@ def approach_to_best_ratio(self, volumes: Vector, free_energies: Vector, initial vr.interpolate_volumes() strains, finer_volumes = vr.strains, vr.out_volumes eulerian_strain = calculate_eulerian_strain(volumes[0], volumes) - _, f_v_tmax = polynomial_least_square_fitting(eulerian_strain, free_energies, strains, self.option) + _, f_v_tmax = polynomial_least_square_fitting( + eulerian_strain, free_energies, strains, self.option + ) p_v_tmax = -np.gradient(f_v_tmax) / np.gradient(finer_volumes) p_desire = gpa_to_ry_b3(self.desired_p_min) # Find the index of the first pressure value that slightly smaller than p_desire. @@ -202,8 +212,9 @@ def approach_to_best_ratio(self, volumes: Vector, free_energies: Vector, initial final_ratio = finer_volumes[idx] / max(volumes) return final_ratio - def refine_grid(self, volumes: Vector, free_energies: Matrix, - ratio: Optional[float] = None) -> Tuple[Vector, Matrix, float]: + def refine_grid( + self, volumes: Vector, free_energies: Matrix, ratio: Optional[float] = None + ) -> Tuple[Vector, Matrix, float]: """ Get the appropriate volume grid for interpolation. Avoid to use a too large volume grid to obtain data, which might lose accuracy. @@ -225,8 +236,14 @@ def refine_grid(self, volumes: Vector, free_energies: Matrix, self._ratio = new_ratio eulerian_strain = calculate_eulerian_strain(volumes[0], volumes) - vr = VolumeExpander(in_volumes=volumes, out_volumes_num=self.dense_volumes_amount, ratio=new_ratio) + vr = VolumeExpander( + in_volumes=volumes, + out_volumes_num=self.dense_volumes_amount, + ratio=new_ratio, + ) vr.interpolate_volumes() # As mentioned in ``VolumeExpander`` doc, call this method immediately. strains, dense_volumes = vr.strains, vr.out_volumes - dense_free_energies = apply_finite_strain_fitting(eulerian_strain, free_energies, strains, self.option) + dense_free_energies = apply_finite_strain_fitting( + eulerian_strain, free_energies, strains, self.option + ) return dense_volumes, dense_free_energies, new_ratio diff --git a/qha/multi_configurations/different_phonon_dos.py b/qha/multi_configurations/different_phonon_dos.py index f397380..45a9fd9 100644 --- a/qha/multi_configurations/different_phonon_dos.py +++ b/qha/multi_configurations/different_phonon_dos.py @@ -19,12 +19,15 @@ from qha.type_aliases import Array4D, Scalar, Vector, Matrix # ===================== What can be exported? ===================== -__all__ = ['PartitionFunction'] +__all__ = ["PartitionFunction"] -K = {'ha': pc['Boltzmann constant in eV/K'][0] / pc['Hartree energy in eV'][0], - 'ry': pc['Boltzmann constant in eV/K'][0] / pc['Rydberg constant times hc in eV'][0], - 'ev': pc['Boltzmann constant in eV/K'][0], - 'SI': pc['Boltzmann constant'][0]}[qha.settings.energy_unit] +K = { + "ha": pc["Boltzmann constant in eV/K"][0] / pc["Hartree energy in eV"][0], + "ry": pc["Boltzmann constant in eV/K"][0] + / pc["Rydberg constant times hc in eV"][0], + "ev": pc["Boltzmann constant in eV/K"][0], + "SI": pc["Boltzmann constant"][0], +}[qha.settings.energy_unit] class PartitionFunction: @@ -55,14 +58,24 @@ class PartitionFunction: :param order: The order of Birch--Murnaghan equation of state fitting, by default, is ``3``. """ - def __init__(self, temperature: Scalar, degeneracies: Vector, q_weights: Matrix, static_energies: Matrix, - volumes: Matrix, frequencies: Array4D, static_only: Optional[bool] = False, - precision: Optional[int] = 500, order: Optional[int] = 3): + def __init__( + self, + temperature: Scalar, + degeneracies: Vector, + q_weights: Matrix, + static_energies: Matrix, + volumes: Matrix, + frequencies: Array4D, + static_only: Optional[bool] = False, + precision: Optional[int] = 500, + order: Optional[int] = 3, + ): if not np.all(np.greater_equal(degeneracies, 0)): - raise ValueError('Degeneracies should all be greater equal than 0!') - if not np.all(np.greater_equal( - q_weights, 0)): # Weights should all be greater equal than 0, otherwise sum will be wrong. - raise ValueError('Weights should all be greater equal than 0!') + raise ValueError("Degeneracies should all be greater equal than 0!") + if not np.all( + np.greater_equal(q_weights, 0) + ): # Weights should all be greater equal than 0, otherwise sum will be wrong. + raise ValueError("Weights should all be greater equal than 0!") self.frequencies = np.array(frequencies) if self.frequencies.ndim != 4: @@ -93,8 +106,18 @@ def unaligned_free_energies_for_each_configuration(self): :return: A matrix, the "raw" free energy of each configuration of each volume. """ configurations_amount, _ = self.volumes.shape - return np.array([free_energy(self.temperature, self.q_weights[i], self.static_energies[i], self.frequencies[i], - self.static_only) for i in range(configurations_amount)]) + return np.array( + [ + free_energy( + self.temperature, + self.q_weights[i], + self.static_energies[i], + self.frequencies[i], + self.static_only, + ) + for i in range(configurations_amount) + ] + ) @LazyProperty def aligned_free_energies_for_each_configuration(self): @@ -103,8 +126,11 @@ def aligned_free_energies_for_each_configuration(self): :return: A matrix, the aligned free energy of each configuration of each volume. """ - return calibrate_energy_on_reference(self.volumes, self.unaligned_free_energies_for_each_configuration, - self.order) + return calibrate_energy_on_reference( + self.volumes, + self.unaligned_free_energies_for_each_configuration, + self.order, + ) @LazyProperty def partition_functions_for_each_configuration(self): @@ -120,12 +146,19 @@ def partition_functions_for_each_configuration(self): try: import mpmath except ImportError: - raise ImportError("Install ``mpmath`` package to use {0} object!".format(self.__class__.__name__)) + raise ImportError( + "Install ``mpmath`` package to use {0} object!".format( + self.__class__.__name__ + ) + ) with mpmath.workprec(self.precision): # shape = (# of configurations, # of volumes for each configuration) exp = np.vectorize(mpmath.exp) - return exp(-self.aligned_free_energies_for_each_configuration / (K * self.temperature)) + return exp( + -self.aligned_free_energies_for_each_configuration + / (K * self.temperature) + ) @LazyProperty def partition_functions_for_all_configurations(self): @@ -141,13 +174,25 @@ def partition_functions_for_all_configurations(self): try: import mpmath except ImportError: - raise ImportError("Install ``mpmath`` package to use {0} object!".format(self.__class__.__name__)) + raise ImportError( + "Install ``mpmath`` package to use {0} object!".format( + self.__class__.__name__ + ) + ) with mpmath.workprec(self.precision): # shape = (# of volumes,) - return np.array([mpmath.exp(d) for d in - logsumexp(-self.aligned_free_energies_for_each_configuration.T / (K * self.temperature), - axis=1, b=self.degeneracies)]) + return np.array( + [ + mpmath.exp(d) + for d in logsumexp( + -self.aligned_free_energies_for_each_configuration.T + / (K * self.temperature), + axis=1, + b=self.degeneracies, + ) + ] + ) def get_free_energies(self): """ @@ -162,8 +207,18 @@ def get_free_energies(self): try: import mpmath except ImportError: - raise ImportError("Install ``mpmath`` package to use {0} object!".format(self.__class__.__name__)) + raise ImportError( + "Install ``mpmath`` package to use {0} object!".format( + self.__class__.__name__ + ) + ) with mpmath.workprec(self.precision): - log_z = np.array([mpmath.log(d) for d in self.partition_functions_for_all_configurations], dtype=float) + log_z = np.array( + [ + mpmath.log(d) + for d in self.partition_functions_for_all_configurations + ], + dtype=float, + ) return -K * self.temperature * log_z diff --git a/qha/multi_configurations/same_phonon_dos.py b/qha/multi_configurations/same_phonon_dos.py index b292f9a..5332946 100644 --- a/qha/multi_configurations/same_phonon_dos.py +++ b/qha/multi_configurations/same_phonon_dos.py @@ -18,12 +18,15 @@ from qha.type_aliases import Array3D, Scalar, Vector, Matrix # ===================== What can be exported? ===================== -__all__ = ['PartitionFunction', 'FreeEnergy'] +__all__ = ["PartitionFunction", "FreeEnergy"] -K = {'ha': pc['Boltzmann constant in eV/K'][0] / pc['Hartree energy in eV'][0], - 'ry': pc['Boltzmann constant in eV/K'][0] / pc['Rydberg constant times hc in eV'][0], - 'ev': pc['Boltzmann constant in eV/K'][0], - 'SI': pc['Boltzmann constant'][0]}[qha.settings.energy_unit] +K = { + "ha": pc["Boltzmann constant in eV/K"][0] / pc["Hartree energy in eV"][0], + "ry": pc["Boltzmann constant in eV/K"][0] + / pc["Rydberg constant times hc in eV"][0], + "ev": pc["Boltzmann constant in eV/K"][0], + "SI": pc["Boltzmann constant"][0], +}[qha.settings.energy_unit] class PartitionFunction: @@ -53,14 +56,23 @@ class PartitionFunction: value, by default, is ``500``. """ - def __init__(self, temperature: Scalar, degeneracies: Vector, q_weights: Vector, static_energies: Matrix, - frequencies: Array3D, precision: Optional[int] = 500): - + def __init__( + self, + temperature: Scalar, + degeneracies: Vector, + q_weights: Vector, + static_energies: Matrix, + frequencies: Array3D, + precision: Optional[int] = 500, + ): if not np.all(np.greater_equal(degeneracies, 0)): - raise ValueError('Degeneracies should all be integers greater equal than 0!') - if not np.all(np.greater_equal(q_weights, - 0)): # Weights should all be greater equal than 0, otherwise sum will be wrong. - raise ValueError('Weights should all be greater equal than 0!') + raise ValueError( + "Degeneracies should all be integers greater equal than 0!" + ) + if not np.all( + np.greater_equal(q_weights, 0) + ): # Weights should all be greater equal than 0, otherwise sum will be wrong. + raise ValueError("Weights should all be greater equal than 0!") self.frequencies = np.array(frequencies) if self.frequencies.ndim != 3: @@ -87,11 +99,22 @@ def _static_part(self) -> Vector: import mpmath except ImportError: raise ImportError( - "You need to install ``mpmath`` package to use {0} object!".format(self.__class__.__name__)) + "You need to install ``mpmath`` package to use {0} object!".format( + self.__class__.__name__ + ) + ) with mpmath.workprec(self.precision): - return np.array([mpmath.exp(d) for d in # shape = (# of volumes for each configuration, 1) - logsumexp(-self.static_energies / (K * self.temperature), axis=1, b=self.degeneracies)]) + return np.array( + [ + mpmath.exp(d) + for d in logsumexp( # shape = (# of volumes for each configuration, 1) + -self.static_energies / (K * self.temperature), + axis=1, + b=self.degeneracies, + ) + ] + ) @property def _harmonic_part(self) -> Vector: @@ -101,7 +124,10 @@ def _harmonic_part(self) -> Vector: :return: The harmonic contribution on the temperature-volume grid. """ log_product_modes: Matrix = np.sum( - log_subsystem_partition_function(self.temperature, self.frequencies), axis=2, dtype=float) + log_subsystem_partition_function(self.temperature, self.frequencies), + axis=2, + dtype=float, + ) return np.exp(np.dot(log_product_modes, self._scaled_q_weights)) # (vol, 1) @LazyProperty @@ -111,7 +137,9 @@ def total(self) -> Vector: :return: The partition function on the temperature-volume grid. """ - return np.multiply(self._static_part, self._harmonic_part) # (vol, 1), product element-wise + return np.multiply( + self._static_part, self._harmonic_part + ) # (vol, 1), product element-wise def get_free_energies(self): """ @@ -126,7 +154,11 @@ def get_free_energies(self): try: import mpmath except ImportError: - raise ImportError("Install ``mpmath`` package to use {0} object!".format(self.__class__.__name__)) + raise ImportError( + "Install ``mpmath`` package to use {0} object!".format( + self.__class__.__name__ + ) + ) with mpmath.workprec(self.precision): log_z = np.array([mpmath.log(d) for d in self.total], dtype=float) @@ -162,13 +194,25 @@ class FreeEnergy: :param order: The order of Birch--Murnaghan equation of state fitting, by default, is ``3``. """ - def __init__(self, temperature: Scalar, degeneracies: Vector, q_weights: Vector, static_energies: Matrix, - volumes: Matrix, frequencies: Array3D, static_only: Optional[bool] = False, order: Optional[int] = 3): + def __init__( + self, + temperature: Scalar, + degeneracies: Vector, + q_weights: Vector, + static_energies: Matrix, + volumes: Matrix, + frequencies: Array3D, + static_only: Optional[bool] = False, + order: Optional[int] = 3, + ): if not np.all(np.greater_equal(degeneracies, 0)): - raise ValueError('Degeneracies should all be integers greater equal than 0!') - if not np.all(np.greater_equal(q_weights, - 0)): # Weights should all be greater equal than 0, otherwise sum will be wrong. - raise ValueError('Weights should all be greater equal than 0!') + raise ValueError( + "Degeneracies should all be integers greater equal than 0!" + ) + if not np.all( + np.greater_equal(q_weights, 0) + ): # Weights should all be greater equal than 0, otherwise sum will be wrong. + raise ValueError("Weights should all be greater equal than 0!") self.frequencies = np.array(frequencies) if self.frequencies.ndim != 3: @@ -194,7 +238,9 @@ def aligned_static_energies_for_each_configuration(self): :return: A matrix of aligned static energy of each configuration of each volume. """ - return calibrate_energy_on_reference(self.volumes, self.static_energies, self.order) + return calibrate_energy_on_reference( + self.volumes, self.static_energies, self.order + ) @LazyProperty def static_part(self) -> Vector: @@ -204,7 +250,9 @@ def static_part(self) -> Vector: :return: The static contribution on the temperature-volume grid. """ kt: float = K * self.temperature # k_B T - inside_exp: Matrix = -self.aligned_static_energies_for_each_configuration.T / kt # exp( E_n(V) / k_B / T ) + inside_exp: Matrix = ( + -self.aligned_static_energies_for_each_configuration.T / kt + ) # exp( E_n(V) / k_B / T ) return -kt * logsumexp(inside_exp, axis=1, b=self.degeneracies) @LazyProperty diff --git a/qha/plotting.py b/qha/plotting.py index cd058e8..17bc909 100644 --- a/qha/plotting.py +++ b/qha/plotting.py @@ -13,30 +13,29 @@ import numpy as np import pandas as pd -matplotlib.use('Agg') +matplotlib.use("Agg") # ===================== What can be exported? ===================== -__all__ = ['Plotter'] +__all__ = ["Plotter"] class Plotter: def __init__(self, user_settings): self.user_settings = user_settings - self.results_folder = pathlib.Path( - self.user_settings['output_directory']) + self.results_folder = pathlib.Path(self.user_settings["output_directory"]) def plot2file(self, plot_what): - if plot_what == 'alpha': + if plot_what == "alpha": self.plot_thermal_expansion() - if plot_what == 'gamma': + if plot_what == "gamma": self.plot_gruneisen_parameter() - if plot_what == 'Cp': + if plot_what == "Cp": self.plot_pressure_specific_heat_capacity() - if plot_what == 'Bt': + if plot_what == "Bt": self.plot_isothermal_bulk_modulus() - if plot_what == 'G': + if plot_what == "G": self.plot_gibbs_free_energy() - if plot_what == 'V': + if plot_what == "V": self.plot_volume() @staticmethod @@ -44,191 +43,194 @@ def plot_to_file(): size_small = 12 size_medium = 14 size_big = 16 - fontweight = 'normal' # normal, bold, bolder, lighter + fontweight = "normal" # normal, bold, bolder, lighter # Defind the font size - plt.rcParams['font.size'] = size_small - plt.rcParams['font.weight'] = fontweight - plt.rcParams['axes.titlesize'] = size_medium - plt.rcParams['axes.labelsize'] = size_medium - plt.rcParams['legend.fontsize'] = size_small - plt.rcParams['figure.titlesize'] = size_big - - plt.rcParams['axes.linewidth'] = 1.5 - plt.rcParams['xtick.major.size'] = 6 - plt.rcParams['xtick.minor.size'] = 3 - plt.rcParams['ytick.major.size'] = 6 - plt.rcParams['ytick.minor.size'] = 3 - plt.rcParams['lines.linewidth'] = 2 + plt.rcParams["font.size"] = size_small + plt.rcParams["font.weight"] = fontweight + plt.rcParams["axes.titlesize"] = size_medium + plt.rcParams["axes.labelsize"] = size_medium + plt.rcParams["legend.fontsize"] = size_small + plt.rcParams["figure.titlesize"] = size_big + + plt.rcParams["axes.linewidth"] = 1.5 + plt.rcParams["xtick.major.size"] = 6 + plt.rcParams["xtick.minor.size"] = 3 + plt.rcParams["ytick.major.size"] = 6 + plt.rcParams["ytick.minor.size"] = 3 + plt.rcParams["lines.linewidth"] = 2 # Define the x and y tick features - plt.rcParams['xtick.labelsize'] = size_small - plt.rcParams['ytick.labelsize'] = size_small - plt.rcParams['xtick.direction'] = 'in' - plt.rcParams['ytick.direction'] = 'in' - plt.rcParams['xtick.minor.visible'] = True - plt.rcParams['ytick.minor.visible'] = True + plt.rcParams["xtick.labelsize"] = size_small + plt.rcParams["ytick.labelsize"] = size_small + plt.rcParams["xtick.direction"] = "in" + plt.rcParams["ytick.direction"] = "in" + plt.rcParams["xtick.minor.visible"] = True + plt.rcParams["ytick.minor.visible"] = True # Define the default figure size in inches - plt.rcParams['figure.figsize'] = [8, 6] + plt.rcParams["figure.figsize"] = [8, 6] # Define Figures overall properties - plt.rcParams['figure.dpi'] = 300 + plt.rcParams["figure.dpi"] = 300 plt.figure() @staticmethod def save_fig(ax, fig_name): fig = ax.get_figure() - fig_name_pdf = str(fig_name) + '.pdf' - fig.savefig(fig_name_pdf, format='pdf', dpi=300, bbox_inches='tight') + fig_name_pdf = str(fig_name) + ".pdf" + fig.savefig(fig_name_pdf, format="pdf", dpi=300, bbox_inches="tight") def plot_thermal_expansion(self): - f_alpha = self.user_settings['alpha'] - alpha = pd.read_csv(f_alpha, sep='\s+', index_col='T(K)\P(GPa)') + f_alpha = self.user_settings["alpha"] + alpha = pd.read_csv(f_alpha, sep="\s+", index_col="T(K)\P(GPa)") self.plot_to_file() ax = alpha.plot(figsize=(6, 5)) - ax.set_xlabel('Temperature (K)') - ax.set_ylabel('Thermal Expansion') - ax.legend(frameon=False, bbox_to_anchor=(1.02, 1), - loc='upper left', title="P (GPa)") + ax.set_xlabel("Temperature (K)") + ax.set_ylabel("Thermal Expansion") + ax.legend( + frameon=False, bbox_to_anchor=(1.02, 1), loc="upper left", title="P (GPa)" + ) ax.set_xlim(alpha.index.min(), alpha.index.max()) - ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) - fig_name = self.results_folder / 'alpha_T' + ax.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) + fig_name = self.results_folder / "alpha_T" self.save_fig(ax, fig_name) def plot_isothermal_bulk_modulus(self): - f_isothermal_bulk_modulus = self.user_settings['Bt'] - btp = pd.read_csv(f_isothermal_bulk_modulus, - sep='\s+', index_col='T(K)\P(GPa)') + f_isothermal_bulk_modulus = self.user_settings["Bt"] + btp = pd.read_csv(f_isothermal_bulk_modulus, sep="\s+", index_col="T(K)\P(GPa)") self.plot_to_file() - ax = btp.plot(y=btp.columns[0], use_index=True, - figsize=(6, 5), color='k') - ax.set_xlabel('Temperature (K)') - ax.set_ylabel('Isothermal Bulk Modulus (GPa)') - ax.legend(frameon=True, loc='best', title="P (GPa)") + ax = btp.plot(y=btp.columns[0], use_index=True, figsize=(6, 5), color="k") + ax.set_xlabel("Temperature (K)") + ax.set_ylabel("Isothermal Bulk Modulus (GPa)") + ax.legend(frameon=True, loc="best", title="P (GPa)") ax.set_xlim(btp.index.min(), btp.index.max()) - fig_name = self.results_folder / 'Bt_T' + fig_name = self.results_folder / "Bt_T" self.save_fig(ax, fig_name) def plot_pressure_specific_heat_capacity(self): - f_cp = self.user_settings['Cp'] - cp = pd.read_csv(f_cp, sep='\s+', index_col='T(K)\P(GPa)') + f_cp = self.user_settings["Cp"] + cp = pd.read_csv(f_cp, sep="\s+", index_col="T(K)\P(GPa)") self.plot_to_file() ax = cp.plot(figsize=(6, 5)) - ax.set_xlabel('Temperature (K)') - ax.set_ylabel('Cp (J/mol K)') - ax.legend(frameon=False, bbox_to_anchor=(1.02, 1), - loc='upper left', title="P (GPa)") + ax.set_xlabel("Temperature (K)") + ax.set_ylabel("Cp (J/mol K)") + ax.legend( + frameon=False, bbox_to_anchor=(1.02, 1), loc="upper left", title="P (GPa)" + ) ax.set_xlim(cp.index.min(), cp.index[:-2].max()) - ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) - fig_name = self.results_folder / 'Cp_T' + ax.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) + fig_name = self.results_folder / "Cp_T" self.save_fig(ax, fig_name) def plot_gruneisen_parameter(self): - f_gamma = self.user_settings['gamma'] - gamma = pd.read_csv(f_gamma, sep='\s+', index_col='T(K)\P(GPa)') + f_gamma = self.user_settings["gamma"] + gamma = pd.read_csv(f_gamma, sep="\s+", index_col="T(K)\P(GPa)") self.plot_to_file() ax = gamma.plot(figsize=(6, 5)) - ax.set_xlabel('Temperature (K)') - ax.set_ylabel('Gruneisen Parameter') - ax.legend(frameon=False, bbox_to_anchor=(1.02, 1), - loc='upper left', title="P (GPa)") + ax.set_xlabel("Temperature (K)") + ax.set_ylabel("Gruneisen Parameter") + ax.legend( + frameon=False, bbox_to_anchor=(1.02, 1), loc="upper left", title="P (GPa)" + ) ax.set_xlim(gamma.index.min(), gamma.index[:-2].max()) - ax.ticklabel_format(style='sci', axis='y', scilimits=(0, 0)) - fig_name = self.results_folder / 'gamma_T' + ax.ticklabel_format(style="sci", axis="y", scilimits=(0, 0)) + fig_name = self.results_folder / "gamma_T" self.save_fig(ax, fig_name) def plot_gibbs_free_energy(self): - f_gibbs = self.user_settings['G'] - gibbs = pd.read_csv(f_gibbs, sep='\s+', index_col='T(K)\P(GPa)') + f_gibbs = self.user_settings["G"] + gibbs = pd.read_csv(f_gibbs, sep="\s+", index_col="T(K)\P(GPa)") self.plot_to_file() ax = gibbs.plot(figsize=(6, 5)) - ax.set_xlabel('Temperature (K)') - ax.set_ylabel('Gibbs Free Energy (ev)') - ax.legend(frameon=False, bbox_to_anchor=(1.02, 1), - loc='upper left', title="P (GPa)") + ax.set_xlabel("Temperature (K)") + ax.set_ylabel("Gibbs Free Energy (ev)") + ax.legend( + frameon=False, bbox_to_anchor=(1.02, 1), loc="upper left", title="P (GPa)" + ) ax.set_xlim(gibbs.index.min(), gibbs.index.max()) - fig_name = self.results_folder / 'G_T' + fig_name = self.results_folder / "G_T" self.save_fig(ax, fig_name) def plot_volume(self): - f_volume = self.user_settings['V'] - volume = pd.read_csv(f_volume, sep='\s+', index_col='T(K)\P(GPa)') + f_volume = self.user_settings["V"] + volume = pd.read_csv(f_volume, sep="\s+", index_col="T(K)\P(GPa)") self.plot_to_file() ax = volume.plot(figsize=(6, 5)) - ax.set_xlabel('Temperature (K)') - ax.set_ylabel(r'Volume ($\AA^3$)') - ax.legend(frameon=False, bbox_to_anchor=(1.02, 1), - loc='upper left', title="P (GPa)") + ax.set_xlabel("Temperature (K)") + ax.set_ylabel(r"Volume ($\AA^3$)") + ax.legend( + frameon=False, bbox_to_anchor=(1.02, 1), loc="upper left", title="P (GPa)" + ) ax.set_xlim(volume.index.min(), volume.index.max()) - fig_name = self.results_folder / 'V_T' + fig_name = self.results_folder / "V_T" self.save_fig(ax, fig_name) def fv_pv(self): - temperature_sample: list = list(map(float, self.user_settings['T4FV'])) + temperature_sample: list = list(map(float, self.user_settings["T4FV"])) f_min = [] v_f_min = [] - f_fitted_volume = self.user_settings['f_tv_fitted'] - volume_tv = pd.read_csv( - f_fitted_volume, sep='\s+', index_col='T(K)\V(A^3)') + f_fitted_volume = self.user_settings["f_tv_fitted"] + volume_tv = pd.read_csv(f_fitted_volume, sep="\s+", index_col="T(K)\V(A^3)") volume_tv.index = volume_tv.index.map(float) volume_vt = volume_tv.T volume_vt.index = volume_vt.index.map(float) volume_t = volume_vt.loc[:, volume_vt.columns.isin(temperature_sample)] volume_v_grid = np.asarray(volume_vt.index, float) - f_nonfitted_volume = self.user_settings['f_tv_non_fitted'] + f_nonfitted_volume = self.user_settings["f_tv_non_fitted"] volume_nonfitted_tv = pd.read_csv( - f_nonfitted_volume, sep='\s+', index_col='T(K)\V(A^3)') + f_nonfitted_volume, sep="\s+", index_col="T(K)\V(A^3)" + ) volume_nonfitted_tv.index = volume_nonfitted_tv.index.map(float) volume_nonfitted_vt = volume_nonfitted_tv.T volume_nonfitted_vt.index = volume_nonfitted_vt.index.map(float) - volume_nonfitted_t = volume_nonfitted_vt.loc[:, volume_nonfitted_vt.columns.isin( - temperature_sample)] + volume_nonfitted_t = volume_nonfitted_vt.loc[ + :, volume_nonfitted_vt.columns.isin(temperature_sample) + ] - f_p_volume = self.user_settings['p_tv_gpa'] - p_volume_tv = pd.read_csv( - f_p_volume, sep='\s+', index_col='T(K)\V(A^3)') + f_p_volume = self.user_settings["p_tv_gpa"] + p_volume_tv = pd.read_csv(f_p_volume, sep="\s+", index_col="T(K)\V(A^3)") p_volume_tv.index = p_volume_tv.index.map(float) p_volume_vt = p_volume_tv.T p_volume_vt.index = p_volume_vt.index.map(float) - p_volume_t = p_volume_vt.loc[:, - p_volume_vt.columns.isin(temperature_sample)] + p_volume_t = p_volume_vt.loc[:, p_volume_vt.columns.isin(temperature_sample)] # Get the 'V0' and 'E0' for selected Temperatures for t in temperature_sample: f_min.append(volume_vt.loc[:, volume_vt.columns == t].min().values) - v_f_min.append( - volume_vt.loc[:, volume_vt.columns == t].idxmin().values) + v_f_min.append(volume_vt.loc[:, volume_vt.columns == t].idxmin().values) f_min, v_f_min = np.asfarray(f_min), np.asfarray(v_f_min) self.plot_to_file() fig, axes = plt.subplots(nrows=1, ncols=2) - volume_t.plot(ax=axes[0], style='b-') - volume_nonfitted_t.plot(ax=axes[0], style='ko') - axes[0].plot(v_f_min, f_min, 'ro-', label=None) + volume_t.plot(ax=axes[0], style="b-") + volume_nonfitted_t.plot(ax=axes[0], style="ko") + axes[0].plot(v_f_min, f_min, "ro-", label=None) - axes[0].set_xlabel(r'Volume ($\AA$)') - axes[0].set_ylabel('Helmholtz Free Energy (eV)') + axes[0].set_xlabel(r"Volume ($\AA$)") + axes[0].set_ylabel("Helmholtz Free Energy (eV)") axes[0].legend().set_visible(False) axes[0].set_xlim(volume_v_grid.min(), volume_v_grid.max()) for t in temperature_sample: x, y = 1.002 * float(volume_v_grid.max()), volume_vt[t].tolist()[0] - axes[0].text(x, y, str(t) + ' K', fontsize=8) + axes[0].text(x, y, str(t) + " K", fontsize=8) axes[1].axhline( - y=self.user_settings['DESIRED_PRESSURES_GPa'][0], color='whitesmoke') + y=self.user_settings["DESIRED_PRESSURES_GPa"][0], color="whitesmoke" + ) axes[1].axhline( - y=self.user_settings['DESIRED_PRESSURES_GPa'][-1], color='whitesmoke') - p_volume_t.plot(ax=axes[1], style='k-') + y=self.user_settings["DESIRED_PRESSURES_GPa"][-1], color="whitesmoke" + ) + p_volume_t.plot(ax=axes[1], style="k-") - axes[1].set_xlabel(r'Volume ($\AA$)') - axes[1].set_ylabel('Pressure (GPa)') + axes[1].set_xlabel(r"Volume ($\AA$)") + axes[1].set_ylabel("Pressure (GPa)") axes[1].legend().set_visible(False) axes[1].set_xlim(volume_v_grid.min(), volume_v_grid.max()) fig.tight_layout(w_pad=3.2, h_pad=0) - fig_name_pdf = str(self.results_folder / 'FVT_PVT.pdf') - fig.savefig(fig_name_pdf, format='pdf', dpi=300, bbox_inches='tight') + fig_name_pdf = str(self.results_folder / "FVT_PVT.pdf") + fig.savefig(fig_name_pdf, format="pdf", dpi=300, bbox_inches="tight") diff --git a/qha/settings.py b/qha/settings.py index 1d5f93e..3cde671 100644 --- a/qha/settings.py +++ b/qha/settings.py @@ -21,23 +21,23 @@ # ===================== What can be exported? ===================== -__all__ = ['DEFAULT_SETTINGS', 'Settings', 'from_yaml'] +__all__ = ["DEFAULT_SETTINGS", "Settings", "from_yaml"] DEFAULT_SETTINGS: Dict[str, Any] = { - 'energy_unit': 'ry', - 'length_unit': 'angstrom', - 'order': 3, # BM fitting order, can be 3, 4 or 5, normally, 3rd order is sufficient. - 'p_min_modifier': 1.0, - 'target': 'parallel', - 'T_MIN': 0, - 'DT_SAMPLE': 10, - 'DELTA_P': 0.1, - 'DELTA_P_SAMPLE': 1, - 'static_only': False, + "energy_unit": "ry", + "length_unit": "angstrom", + "order": 3, # BM fitting order, can be 3, 4 or 5, normally, 3rd order is sufficient. + "p_min_modifier": 1.0, + "target": "parallel", + "T_MIN": 0, + "DT_SAMPLE": 10, + "DELTA_P": 0.1, + "DELTA_P_SAMPLE": 1, + "static_only": False, # output setting - 'output_directory': './results/', - 'T4FV': ['0', '300'], - 'high_verbosity': False + "output_directory": "./results/", + "T4FV": ["0", "300"], + "high_verbosity": False, } @@ -68,13 +68,15 @@ class Settings(collections.ChainMap): :param user_settings: It can be a single dictionary or a tuple of dictionaries, with strings as their keys. """ - def __init__(self, *user_settings: Union[Dict[str, Any], Tuple[Dict[str, Any], ...]]): + def __init__( + self, *user_settings: Union[Dict[str, Any], Tuple[Dict[str, Any], ...]] + ): super().__init__(*user_settings, DEFAULT_SETTINGS) def to_yaml_file(self, filename: str): - if not filename.endswith('.yaml'): - filename += '.yaml' - with open(filename, 'w') as f: + if not filename.endswith(".yaml"): + filename += ".yaml" + with open(filename, "w") as f: yaml.dump(self.maps, f) @@ -85,9 +87,9 @@ def from_yaml(filename: str) -> Settings: :param filename: The name of the YAML file. :return: A ``Settings`` class. """ - with open(filename, 'r') as f: + with open(filename, "r") as f: return Settings(load(f, Loader=Loader)) global energy_unit -energy_unit = DEFAULT_SETTINGS['energy_unit'] +energy_unit = DEFAULT_SETTINGS["energy_unit"] diff --git a/qha/single_configuration.py b/qha/single_configuration.py index 9fa3f1a..f4cfb0c 100644 --- a/qha/single_configuration.py +++ b/qha/single_configuration.py @@ -14,12 +14,17 @@ from qha.type_aliases import Scalar, Vector, Matrix, Array3D # ===================== What can be exported? ===================== -__all__ = ['free_energy', 'HOFreeEnergySampler'] +__all__ = ["free_energy", "HOFreeEnergySampler"] @jit(float64[:](float64, float64[:], float64[:], float64[:, :, :], boolean), cache=True) -def free_energy(temperature: Scalar, q_weights: Vector, static_energies: Vector, frequencies: Array3D, - static_only: bool = False) -> Vector: +def free_energy( + temperature: Scalar, + q_weights: Vector, + static_energies: Vector, + frequencies: Array3D, + static_only: bool = False, +) -> Vector: """ The total free energy at a certain temperature. @@ -36,13 +41,15 @@ def free_energy(temperature: Scalar, q_weights: Vector, static_energies: Vector, volumes. The default unit is the same as in function ``ho_free_energy``. """ if not np.all(np.greater_equal(q_weights, 0)): - raise ValueError('Weights should all be greater equal than 0!') + raise ValueError("Weights should all be greater equal than 0!") if static_only: return static_energies scaled_q_weights: Vector = q_weights / np.sum(q_weights) - vibrational_energies: Vector = np.dot(ho_free_energy(temperature, frequencies).sum(axis=2), scaled_q_weights) + vibrational_energies: Vector = np.dot( + ho_free_energy(temperature, frequencies).sum(axis=2), scaled_q_weights + ) return static_energies + vibrational_energies @@ -65,11 +72,13 @@ class HOFreeEnergySampler: """ def __init__(self, temperature: float, q_weights: Vector, frequencies: Array3D): - self.temperature = temperature # Fix temperature and volume, just sample q-points and bands + self.temperature = ( + temperature # Fix temperature and volume, just sample q-points and bands + ) self.q_weights = np.array(q_weights, dtype=float) self.omegas = np.array(frequencies, dtype=float) if not np.all(np.greater_equal(q_weights, 0)): - raise ValueError('Weights should all be greater equal than 0!') + raise ValueError("Weights should all be greater equal than 0!") else: self._scaled_q_weights = q_weights / np.sum(q_weights) @@ -98,7 +107,10 @@ def on_volume(self, i: int) -> Scalar: :param i: An integer labeling :math:`i` th volume. :return: The accumulated free energy on the :math:`i` th volume. """ - return np.vdot(ho_free_energy(self.temperature, self.omegas[i]).sum(axis=1), self._scaled_q_weights) + return np.vdot( + ho_free_energy(self.temperature, self.omegas[i]).sum(axis=1), + self._scaled_q_weights, + ) @LazyProperty def on_all_volumes(self) -> Vector: @@ -110,4 +122,7 @@ def on_all_volumes(self) -> Vector: """ # First calculate free energies on a 3D array, then sum along the third axis (bands), # at last contract weights and free energies on all q-points. - return np.dot(ho_free_energy(self.temperature, self.omegas).sum(axis=2), self._scaled_q_weights) + return np.dot( + ho_free_energy(self.temperature, self.omegas).sum(axis=2), + self._scaled_q_weights, + ) diff --git a/qha/statmech.py b/qha/statmech.py index 659a834..f42e97b 100644 --- a/qha/statmech.py +++ b/qha/statmech.py @@ -13,21 +13,35 @@ import qha.settings # ===================== What can be exported? ===================== -__all__ = ['ho_free_energy', 'subsystem_partition_function', 'log_subsystem_partition_function'] - -K = {'ha': pc['Boltzmann constant in eV/K'][0] / pc['Hartree energy in eV'][0], - 'ry': pc['Boltzmann constant in eV/K'][0] / pc['Rydberg constant times hc in eV'][0], - 'ev': pc['Boltzmann constant in eV/K'][0], - 'SI': pc['Boltzmann constant'][0]}[qha.settings.energy_unit] - -HBAR = {'ha': 100 / pc['electron volt-inverse meter relationship'][0] / pc['Hartree energy in eV'][0], - 'ry': 100 / pc['electron volt-inverse meter relationship'][0] / pc['Rydberg constant times hc in eV'][0], - 'ev': 100 / pc['electron volt-inverse meter relationship'][0], - 'SI': 100 / pc['electron volt-inverse meter relationship'][0] / pc['joule-electron volt relationship'][0]}[ - qha.settings.energy_unit] - - -@vectorize([float64(float64, float64)], nopython=True, target='parallel', cache=True) +__all__ = [ + "ho_free_energy", + "subsystem_partition_function", + "log_subsystem_partition_function", +] + +K = { + "ha": pc["Boltzmann constant in eV/K"][0] / pc["Hartree energy in eV"][0], + "ry": pc["Boltzmann constant in eV/K"][0] + / pc["Rydberg constant times hc in eV"][0], + "ev": pc["Boltzmann constant in eV/K"][0], + "SI": pc["Boltzmann constant"][0], +}[qha.settings.energy_unit] + +HBAR = { + "ha": 100 + / pc["electron volt-inverse meter relationship"][0] + / pc["Hartree energy in eV"][0], + "ry": 100 + / pc["electron volt-inverse meter relationship"][0] + / pc["Rydberg constant times hc in eV"][0], + "ev": 100 / pc["electron volt-inverse meter relationship"][0], + "SI": 100 + / pc["electron volt-inverse meter relationship"][0] + / pc["joule-electron volt relationship"][0], +}[qha.settings.energy_unit] + + +@vectorize([float64(float64, float64)], nopython=True, target="parallel", cache=True) def ho_free_energy(temperature, frequency): """ Calculate Helmholtz free energy of a single harmonic oscillator at a specific temperature. @@ -46,7 +60,7 @@ def ho_free_energy(temperature, frequency): return 1 / 2 * hw + kt * np.log(1 - np.exp(-hw / kt)) -@vectorize([float64(float64, float64)], nopython=True, target='parallel', cache=True) +@vectorize([float64(float64, float64)], nopython=True, target="parallel", cache=True) def subsystem_partition_function(temperature, frequency): """ Calculate the subsystem partition function of a single harmonic oscillator at a specific temperature. @@ -64,7 +78,7 @@ def subsystem_partition_function(temperature, frequency): return np.exp(x / 2) / (1 - np.exp(x)) -@vectorize([float64(float64, float64)], nopython=True, target='parallel', cache=True) +@vectorize([float64(float64, float64)], nopython=True, target="parallel", cache=True) def log_subsystem_partition_function(temperature, frequency): """ Calculate the natural logarithm of the subsystem partition function of a single harmonic oscillator diff --git a/qha/tests/test_different_phonon_dos.py b/qha/tests/test_different_phonon_dos.py index 56624cb..e69355f 100644 --- a/qha/tests/test_different_phonon_dos.py +++ b/qha/tests/test_different_phonon_dos.py @@ -12,68 +12,235 @@ class TestPartitionFunction(unittest.TestCase): def setUp(self) -> None: - os.chdir(Path(__file__).parent.parent.parent / 'examples/ice VII/') + os.chdir(Path(__file__).parent.parent.parent / "examples/ice VII/") self.user_settings = {} settings = from_yaml("settings.yaml") - for key in ('input', 'calculation', - 'thermodynamic_properties', 'static_only', 'energy_unit', - 'T_MIN', 'NT', 'DT', 'DT_SAMPLE', - 'P_MIN', 'NTV', 'DELTA_P', 'DELTA_P_SAMPLE', - 'volume_ratio', 'order', 'p_min_modifier', - 'T4FV', 'output_directory', 'high_verbosity'): + for key in ( + "input", + "calculation", + "thermodynamic_properties", + "static_only", + "energy_unit", + "T_MIN", + "NT", + "DT", + "DT_SAMPLE", + "P_MIN", + "NTV", + "DELTA_P", + "DELTA_P_SAMPLE", + "volume_ratio", + "order", + "p_min_modifier", + "T4FV", + "output_directory", + "high_verbosity", + ): try: self.user_settings.update({key: settings[key]}) except KeyError: continue self.calculator = DifferentPhDOSCalculator(self.user_settings) self.calculator.read_input() - self.partition_function = PartitionFunction(1000, self.calculator.degeneracies, self.calculator.q_weights, - self.calculator.static_energies, self.calculator._volumes, - self.calculator.frequencies) + self.partition_function = PartitionFunction( + 1000, + self.calculator.degeneracies, + self.calculator.q_weights, + self.calculator.static_energies, + self.calculator._volumes, + self.calculator.frequencies, + ) def test_parse_settings(self): - self.assertDictEqual(self.user_settings, { - 'input': {'input_01': 6, 'input_02': 96, 'input_03': 96, 'input_04': 24, 'input_05': 24, 'input_06': 192, - 'input_07': 384, 'input_08': 24, 'input_09': 384, 'input_10': 96, 'input_11': 192, - 'input_12': 384, 'input_13': 48, 'input_14': 96, 'input_15': 48, 'input_16': 96, 'input_17': 96, - 'input_18': 384, 'input_19': 384, 'input_20': 48, 'input_21': 384, 'input_22': 96, - 'input_23': 384, 'input_24': 96, 'input_25': 192, 'input_26': 192, 'input_27': 192, - 'input_28': 384, 'input_29': 192, 'input_30': 192, - 'input_31': 96, 'input_32': 192, 'input_33': 192, 'input_34': 384, 'input_35': 384, - 'input_36': 192, 'input_37': 24, 'input_38': 48, 'input_39': 96, 'input_40': 48, 'input_41': 384, - 'input_42': 12, 'input_43': 96, 'input_44': 48, 'input_45': 192, 'input_46': 12, 'input_47': 96, - 'input_48': 48, 'input_49': 6, 'input_50': 96, 'input_51': 24, 'input_52': 24}, - 'calculation': 'different phonon dos', - 'thermodynamic_properties': ['F', 'G', 'U', 'H', 'V', 'alpha', 'gamma', 'Cp', 'Cv', 'Bt', 'Btp', 'Bs'], - 'static_only': False, 'energy_unit': 'ry', 'T_MIN': 0, 'NT': 401, 'DT': 1, 'DT_SAMPLE': 1, 'P_MIN': 0, - 'NTV': 701, 'DELTA_P': 0.1, 'DELTA_P_SAMPLE': 1, 'order': 3, 'p_min_modifier': 1.0, 'T4FV': ['0', '300'], - 'output_directory': './results/', 'high_verbosity': True}) + self.assertDictEqual( + self.user_settings, + { + "input": { + "input_01": 6, + "input_02": 96, + "input_03": 96, + "input_04": 24, + "input_05": 24, + "input_06": 192, + "input_07": 384, + "input_08": 24, + "input_09": 384, + "input_10": 96, + "input_11": 192, + "input_12": 384, + "input_13": 48, + "input_14": 96, + "input_15": 48, + "input_16": 96, + "input_17": 96, + "input_18": 384, + "input_19": 384, + "input_20": 48, + "input_21": 384, + "input_22": 96, + "input_23": 384, + "input_24": 96, + "input_25": 192, + "input_26": 192, + "input_27": 192, + "input_28": 384, + "input_29": 192, + "input_30": 192, + "input_31": 96, + "input_32": 192, + "input_33": 192, + "input_34": 384, + "input_35": 384, + "input_36": 192, + "input_37": 24, + "input_38": 48, + "input_39": 96, + "input_40": 48, + "input_41": 384, + "input_42": 12, + "input_43": 96, + "input_44": 48, + "input_45": 192, + "input_46": 12, + "input_47": 96, + "input_48": 48, + "input_49": 6, + "input_50": 96, + "input_51": 24, + "input_52": 24, + }, + "calculation": "different phonon dos", + "thermodynamic_properties": [ + "F", + "G", + "U", + "H", + "V", + "alpha", + "gamma", + "Cp", + "Cv", + "Bt", + "Btp", + "Bs", + ], + "static_only": False, + "energy_unit": "ry", + "T_MIN": 0, + "NT": 401, + "DT": 1, + "DT_SAMPLE": 1, + "P_MIN": 0, + "NTV": 701, + "DELTA_P": 0.1, + "DELTA_P_SAMPLE": 1, + "order": 3, + "p_min_modifier": 1.0, + "T4FV": ["0", "300"], + "output_directory": "./results/", + "high_verbosity": True, + }, + ) def test_parameters(self): - self.assertEqual(self.calculator.degeneracies, ( - 6, 96, 96, 24, 24, 192, 384, 24, 384, 96, 192, 384, 48, 96, 48, 96, 96, 384, 384, 48, 384, 96, 384, 96, 192, - 192, 192, 384, 192, 192, 96, 192, 192, 384, 384, 192, 24, 48, 96, 48, 384, 12, 96, 48, 192, 12, 96, 48, 6, - 96, 24, 24)) - self.assertEqual(self.calculator.frequencies.shape, - (52, 6, 1, 144)) # frequency on each (configuration, volume, q-point and mode) - np.testing.assert_array_equal(self.calculator.q_weights, - np.ones((52, 1))) # We only have Γ points, whose weight is 1. - np.testing.assert_array_equal(self.calculator.volumes, - [2290.7412, 2179.0756, 1666.2973, 1362.8346, 1215.6528, 1120.4173]) + self.assertEqual( + self.calculator.degeneracies, + ( + 6, + 96, + 96, + 24, + 24, + 192, + 384, + 24, + 384, + 96, + 192, + 384, + 48, + 96, + 48, + 96, + 96, + 384, + 384, + 48, + 384, + 96, + 384, + 96, + 192, + 192, + 192, + 384, + 192, + 192, + 96, + 192, + 192, + 384, + 384, + 192, + 24, + 48, + 96, + 48, + 384, + 12, + 96, + 48, + 192, + 12, + 96, + 48, + 6, + 96, + 24, + 24, + ), + ) + self.assertEqual( + self.calculator.frequencies.shape, (52, 6, 1, 144) + ) # frequency on each (configuration, volume, q-point and mode) + np.testing.assert_array_equal( + self.calculator.q_weights, np.ones((52, 1)) + ) # We only have Γ points, whose weight is 1. + np.testing.assert_array_equal( + self.calculator.volumes, + [2290.7412, 2179.0756, 1666.2973, 1362.8346, 1215.6528, 1120.4173], + ) self.assertEqual(self.calculator._volumes.shape, (52, 6)) - self.assertEqual(self.calculator.static_energies.shape, - (52, 6)) # static energy on each (configuration, volume) + self.assertEqual( + self.calculator.static_energies.shape, (52, 6) + ) # static energy on each (configuration, volume) def test_partition_function(self): - self.assertEqual(self.partition_function.unaligned_free_energies_for_each_configuration.shape, - (52, 6)) # (# of configurations, # of volumes) - self.assertEqual(self.partition_function.aligned_free_energies_for_each_configuration.shape, - (52, 6)) # (# of configurations, # of volumes) - self.assertEqual(self.partition_function.partition_functions_for_each_configuration.shape, - (52, 6)) # (# of configurations, # of volumes) - self.assertEqual(self.partition_function.partition_functions_for_all_configurations.shape, - (6,)) # (# of volumes,) - np.testing.assert_array_almost_equal(self.partition_function.get_free_energies(), - [-550.74580132, -550.70964062, -550.37436235, -549.87365787, -549.43586034, - -549.03029969]) + self.assertEqual( + self.partition_function.unaligned_free_energies_for_each_configuration.shape, + (52, 6), + ) # (# of configurations, # of volumes) + self.assertEqual( + self.partition_function.aligned_free_energies_for_each_configuration.shape, + (52, 6), + ) # (# of configurations, # of volumes) + self.assertEqual( + self.partition_function.partition_functions_for_each_configuration.shape, + (52, 6), + ) # (# of configurations, # of volumes) + self.assertEqual( + self.partition_function.partition_functions_for_all_configurations.shape, + (6,), + ) # (# of volumes,) + np.testing.assert_array_almost_equal( + self.partition_function.get_free_energies(), + [ + -550.74580132, + -550.70964062, + -550.37436235, + -549.87365787, + -549.43586034, + -549.03029969, + ], + ) diff --git a/qha/tests/test_overall_run.py b/qha/tests/test_overall_run.py index be51502..e3c6cf2 100644 --- a/qha/tests/test_overall_run.py +++ b/qha/tests/test_overall_run.py @@ -11,28 +11,44 @@ class TestOverallRun(unittest.TestCase): def setUp(self): - self.root_directory = pathlib.Path(__file__).parent.parent.parent / 'examples' - self.command = 'qha run' - self.fixed_directory = 'results.benchmark' - self.new_results_directory = 'results.plot' + self.root_directory = pathlib.Path(__file__).parent.parent.parent / "examples" + self.command = "qha run" + self.fixed_directory = "results.benchmark" + self.new_results_directory = "results.plot" @staticmethod def listdir_nohidden(txt_path): for f in os.listdir(txt_path): - if not f.startswith('.') and f.find('_tp') > 0: + if not f.startswith(".") and f.find("_tp") > 0: yield f def compare_results(self, path_results_benchmark, path_results_new): d = dict() for f in self.listdir_nohidden(path_results_benchmark): - d.update({f: pd.read_csv(str(path_results_benchmark) + '/' + f, sep='\s+', index_col='T(K)\P(GPa)')}) + d.update( + { + f: pd.read_csv( + str(path_results_benchmark) + "/" + f, + sep="\s+", + index_col="T(K)\P(GPa)", + ) + } + ) d0 = dict() for f in self.listdir_nohidden(path_results_new): - d0.update({f: pd.read_csv(str(path_results_new) + '/' + f, sep='\s+', index_col='T(K)\P(GPa)')}) + d0.update( + { + f: pd.read_csv( + str(path_results_new) + "/" + f, + sep="\s+", + index_col="T(K)\P(GPa)", + ) + } + ) for k, v in d.items(): - print(k + ':', np.max(np.abs(v.as_matrix() - d0[k].as_matrix()))) + print(k + ":", np.max(np.abs(v.as_matrix() - d0[k].as_matrix()))) def prepare_results_new(self, path_results_new, path_results, path_run_command): os.makedirs(path_results_new, exist_ok=False) @@ -40,11 +56,11 @@ def prepare_results_new(self, path_results_new, path_results, path_run_command): command = "mv *.txt " + self.new_results_directory subprocess.run(command, shell=True, cwd=path_results) - def test_silicon(self, test_directory='silicon'): + def test_silicon(self, test_directory="silicon"): print("testing the examples/silicon") path_run_command = self.root_directory / test_directory - path_results = path_run_command / 'results' + path_results = path_run_command / "results" path_results_fixed = path_results / self.fixed_directory path_results_new = path_results / self.new_results_directory @@ -52,11 +68,11 @@ def test_silicon(self, test_directory='silicon'): self.compare_results(path_results_fixed, path_results_new) - def test_ice(self, test_directory='ice VII'): + def test_ice(self, test_directory="ice VII"): print("testing the examples/ice VII") path_run_command = self.root_directory / test_directory - path_results = path_run_command / 'results' + path_results = path_run_command / "results" path_results_benchmark = path_results / self.fixed_directory path_results_new = path_results / self.new_results_directory @@ -74,5 +90,5 @@ def test_ice(self, test_directory='ice VII'): # self.compare_results(path_results / 'results.bmf5', path_results / 'results.bfm5') -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/qha/tests/test_read_input.py b/qha/tests/test_read_input.py index 976b30a..deda758 100644 --- a/qha/tests/test_read_input.py +++ b/qha/tests/test_read_input.py @@ -11,22 +11,35 @@ class TestReadInput(unittest.TestCase): def setUp(self): - self.dir = pathlib.Path(__file__).parent.parent.parent / 'examples' + self.dir = pathlib.Path(__file__).parent.parent.parent / "examples" def test_read_si_input(self): - file_path = self.dir / 'silicon/input' + file_path = self.dir / "silicon/input" nm, volumes, static_energies, frequencies, q_weights = read_input(file_path) self.assertEqual(nm, 2) - np.testing.assert_array_equal(volumes, - [320.5259, 311.4549, 302.5568, 293.8297, 285.2721, 276.8823, 268.6586, 260.5994, - 252.703, 244.9677, 237.392]) + np.testing.assert_array_equal( + volumes, + [ + 320.5259, + 311.4549, + 302.5568, + 293.8297, + 285.2721, + 276.8823, + 268.6586, + 260.5994, + 252.703, + 244.9677, + 237.392, + ], + ) def test_read_ice_input(self): for i in range(1, 53): - file_path = self.dir / "{0}{1:02d}".format('ice VII/input_', i) + file_path = self.dir / "{0}{1:02d}".format("ice VII/input_", i) nm, volumes, static_energies, frequencies, q_weights = read_input(file_path) self.assertEqual(nm, 16) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/qha/tests/test_samevdos_overall.py b/qha/tests/test_samevdos_overall.py index 8ba3f21..d137cb4 100644 --- a/qha/tests/test_samevdos_overall.py +++ b/qha/tests/test_samevdos_overall.py @@ -11,28 +11,44 @@ class TestSamePhononDOS(unittest.TestCase): def setUp(self): - self.root_directory = pathlib.Path(__file__).parent.parent.parent / 'examples' - self.command = 'qha run' - self.fixed_directory = 'results.benchmark' - self.new_results_directory = 'results.same_vdos_new01' + self.root_directory = pathlib.Path(__file__).parent.parent.parent / "examples" + self.command = "qha run" + self.fixed_directory = "results.benchmark" + self.new_results_directory = "results.same_vdos_new01" @staticmethod def listdir_nohidden(txt_path): for f in os.listdir(txt_path): - if not f.startswith('.') and f.find('_tp') > 0: + if not f.startswith(".") and f.find("_tp") > 0: yield f def compare_results(self, path_results_benchmark, path_results_new): d = dict() for f in self.listdir_nohidden(path_results_benchmark): - d.update({f: pd.read_csv(str(path_results_benchmark) + '/' + f, sep='\s+', index_col='T(K)\P(GPa)')}) + d.update( + { + f: pd.read_csv( + str(path_results_benchmark) + "/" + f, + sep="\s+", + index_col="T(K)\P(GPa)", + ) + } + ) d0 = dict() for f in self.listdir_nohidden(path_results_new): - d0.update({f: pd.read_csv(str(path_results_new) + '/' + f, sep='\s+', index_col='T(K)\P(GPa)')}) + d0.update( + { + f: pd.read_csv( + str(path_results_new) + "/" + f, + sep="\s+", + index_col="T(K)\P(GPa)", + ) + } + ) for k, v in d.items(): - print(k + ':', np.max(np.abs(v.as_matrix() - d0[k].as_matrix()))) + print(k + ":", np.max(np.abs(v.as_matrix() - d0[k].as_matrix()))) def prepare_results_new(self, path_results_new, path_results, path_run_command): os.makedirs(path_results_new, exist_ok=False) @@ -40,19 +56,25 @@ def prepare_results_new(self, path_results_new, path_results, path_run_command): command = "mv *.txt " + self.new_results_directory subprocess.run(command, shell=True, cwd=path_results) - def test_samevdos(self, test_directory='results/same_vdos_example_ol'): - + def test_samevdos(self, test_directory="results/same_vdos_example_ol"): print("testing the results/same_vdos_example_ol same vdos ") path_run_command = self.root_directory / test_directory - path_results = path_run_command / 'results' + path_results = path_run_command / "results" path_results_fixed = path_results / self.fixed_directory path_run_command = self.root_directory / test_directory - path_results = path_run_command / 'results' - self.compare_results(path_results_fixed, path_results / 'results.same_vdos_new01') - self.compare_results(path_results / 'results.same_vdos_new', path_results / 'results.same_vdos_new01') - self.compare_results(path_results / 'results.diff_vdos', path_results / 'results.diff_vdos_new') + path_results = path_run_command / "results" + self.compare_results( + path_results_fixed, path_results / "results.same_vdos_new01" + ) + self.compare_results( + path_results / "results.same_vdos_new", + path_results / "results.same_vdos_new01", + ) + self.compare_results( + path_results / "results.diff_vdos", path_results / "results.diff_vdos_new" + ) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/qha/tests/test_single_configuration.py b/qha/tests/test_single_configuration.py index bc6abdc..45445c6 100644 --- a/qha/tests/test_single_configuration.py +++ b/qha/tests/test_single_configuration.py @@ -14,5 +14,5 @@ def test_ho_free_energy(self): self.assertAlmostEqual(ho_free_energy(100, 1000), 0.0045562989049199466) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/qha/tests/test_unit_conversion.py b/qha/tests/test_unit_conversion.py index 79e9d72..6d2e2a8 100644 --- a/qha/tests/test_unit_conversion.py +++ b/qha/tests/test_unit_conversion.py @@ -13,8 +13,8 @@ class TestUnitConversion(unittest.TestCase): def test_j_to_ev(self): - self.assertEqual(j_to_ev(1), 6.241509125883258e+18) - self.assertEqual(j_to_ev(900), 5.617358213294932e+21) + self.assertEqual(j_to_ev(1), 6.241509125883258e18) + self.assertEqual(j_to_ev(900), 5.617358213294932e21) def test_gpa_to_megabar(self): self.assertEqual(gpa_to_megabar(1), 0.01) @@ -69,11 +69,13 @@ def isufunc(obj): return True return False - all_ufuncs: List[np.ufunc] = [obj[1] for obj in getmembers(qha.unit_conversion, isufunc)] + all_ufuncs: List[np.ufunc] = [ + obj[1] for obj in getmembers(qha.unit_conversion, isufunc) + ] for ufunc in all_ufuncs: self.assertEqual(ufunc(vector).size, 50) self.assertEqual(ufunc(grid).shape, (3, 3)) -if __name__ == '__main__': +if __name__ == "__main__": unittest.main() diff --git a/qha/thermodynamics.py b/qha/thermodynamics.py index c616fb4..fbda966 100644 --- a/qha/thermodynamics.py +++ b/qha/thermodynamics.py @@ -14,17 +14,17 @@ # ===================== What can be exported? ===================== __all__ = [ - 'pressure', - 'entropy', - 'thermodynamic_potentials', - 'volume', - 'thermal_expansion_coefficient', - 'gruneisen_parameter', - 'isothermal_bulk_modulus', - 'adiabatic_bulk_modulus', - 'bulk_modulus_derivative', - 'isobaric_heat_capacity', - 'volumetric_heat_capacity' + "pressure", + "entropy", + "thermodynamic_potentials", + "volume", + "thermal_expansion_coefficient", + "gruneisen_parameter", + "isothermal_bulk_modulus", + "adiabatic_bulk_modulus", + "bulk_modulus_derivative", + "isobaric_heat_capacity", + "volumetric_heat_capacity", ] @@ -37,7 +37,9 @@ def calculate_derivatives(xs: Vector, fs: Matrix) -> Matrix: :return: A matrix, with shape :math:`(N_x, _)`. """ if xs.ndim > 1 or fs.ndim < 2: - raise ValueError('The argument *xs* should be a 1D array and *ys* should be a 2D array!') + raise ValueError( + "The argument *xs* should be a 1D array and *ys* should be a 2D array!" + ) return np.gradient(fs, axis=0) / np.gradient(xs)[:, None] # df(x)/dx. @@ -72,7 +74,9 @@ def entropy(temperature: Vector, free_energies: Matrix) -> Matrix: return -calculate_derivatives(temperature, free_energies) -def thermodynamic_potentials(temperature: Vector, vs: Vector, free_energies: Matrix, ps: Matrix): +def thermodynamic_potentials( + temperature: Vector, vs: Vector, free_energies: Matrix, ps: Matrix +): """ Calculate the enthalpy :math:`H(T, V)`, the internal energy :math:`U(T, V)`, and the Gibbs free energy :math:`G` on a :math:`(T, V)` grid from Helmholtz free energy :math:`F(T, V)` by @@ -93,11 +97,15 @@ def thermodynamic_potentials(temperature: Vector, vs: Vector, free_energies: Mat """ g: Matrix = free_energies + ps * vs # G(T,V) = F(T,V) + V * P(T,V) - u: Matrix = free_energies + entropy(temperature, free_energies) * temperature.reshape(-1, 1) # U(T,V) = F(T,V) + T * S(T,V) + u: Matrix = free_energies + entropy( + temperature, free_energies + ) * temperature.reshape( + -1, 1 + ) # U(T,V) = F(T,V) + T * S(T,V) h: Matrix = u + ps * vs # H(T,V) = U(T,V) + V * P(T,V) - return {'U': u, 'H': h, 'G': g} + return {"U": u, "H": h, "G": g} def volume(vs: Vector, desired_ps: Vector, ps: Matrix) -> Matrix: @@ -171,7 +179,9 @@ def isothermal_bulk_modulus(vs: Vector, ps: Matrix) -> Matrix: return -np.gradient(ps, axis=1) / np.gradient(vs) * vs -def adiabatic_bulk_modulus(bt: Matrix, alpha: Matrix, gamma: Matrix, temperature: Vector) -> Matrix: +def adiabatic_bulk_modulus( + bt: Matrix, alpha: Matrix, gamma: Matrix, temperature: Vector +) -> Matrix: """ Calculate the adiabatic bulk modulus by @@ -209,7 +219,9 @@ def bulk_modulus_derivative(ps: Vector, bt: Matrix) -> Matrix: return calculate_derivatives(ps, bt.T).T -def isobaric_heat_capacity(cv: Matrix, alpha: Matrix, gamma: Matrix, temperature: Vector) -> Matrix: +def isobaric_heat_capacity( + cv: Matrix, alpha: Matrix, gamma: Matrix, temperature: Vector +) -> Matrix: """ Calculate the isobaric heat capacity by diff --git a/qha/tools.py b/qha/tools.py index 1e1e377..017caea 100644 --- a/qha/tools.py +++ b/qha/tools.py @@ -17,8 +17,16 @@ from qha.type_aliases import Matrix, Scalar, Vector # ===================== What can be exported? ===================== -__all__ = ['find_nearest', 'vectorized_find_nearest', 'lagrange3', 'lagrange4', 'is_monotonic_decreasing', - 'is_monotonic_increasing', 'arange', 'calibrate_energy_on_reference'] +__all__ = [ + "find_nearest", + "vectorized_find_nearest", + "lagrange3", + "lagrange4", + "is_monotonic_decreasing", + "is_monotonic_increasing", + "arange", + "calibrate_energy_on_reference", +] def lagrange4(xs: Vector, ys: Vector) -> Callable[[float], float]: @@ -33,12 +41,12 @@ def lagrange4(xs: Vector, ys: Vector) -> Callable[[float], float]: :math:`[\\text{min}(xs), \\text{max}(xs)]`, where :math:`xs` denotes the argument *xs*. """ if not len(xs) == len(ys) == 4: - raise ValueError('The *xs* and *ys* must both have length 4!') + raise ValueError("The *xs* and *ys* must both have length 4!") x0, x1, x2, x3 = xs y0, y1, y2, y3 = ys - @vectorize(["float64(float64)"], target='parallel') + @vectorize(["float64(float64)"], target="parallel") def f(x: float) -> float: """ A helper function that only does the evaluation. @@ -46,10 +54,12 @@ def f(x: float) -> float: :param x: The variable on which the Lagrange polynomial is going to be applied. :return: The value of the Lagrange polynomial on :math:`x`, i.e., :math:`L(x)`. """ - return (x - x1) * (x - x2) * (x - x3) / (x0 - x1) / (x0 - x2) / (x0 - x3) * y0 + \ - (x - x0) * (x - x2) * (x - x3) / (x1 - x0) / (x1 - x2) / (x1 - x3) * y1 + \ - (x - x0) * (x - x1) * (x - x3) / (x2 - x0) / (x2 - x1) / (x2 - x3) * y2 + \ - (x - x0) * (x - x1) * (x - x2) / (x3 - x0) / (x3 - x1) / (x3 - x2) * y3 + return ( + (x - x1) * (x - x2) * (x - x3) / (x0 - x1) / (x0 - x2) / (x0 - x3) * y0 + + (x - x0) * (x - x2) * (x - x3) / (x1 - x0) / (x1 - x2) / (x1 - x3) * y1 + + (x - x0) * (x - x1) * (x - x3) / (x2 - x0) / (x2 - x1) / (x2 - x3) * y2 + + (x - x0) * (x - x1) * (x - x2) / (x3 - x0) / (x3 - x1) / (x3 - x2) * y3 + ) return f @@ -74,15 +84,15 @@ def lagrange3(xs: Vector, ys: Vector) -> Callable[[float], float]: :math:`[\\text{min}(xs), \\text{max}(xs)]`, where :math:`xs` denotes the argument *xs*. """ if not len(xs) == len(ys) == 3: - raise ValueError('The *xs* and *ys* must both have length 3!') + raise ValueError("The *xs* and *ys* must both have length 3!") if len(set(xs)) < 3: # The ``set`` will remove duplicated items. - raise ValueError('Some elements of *xs* are duplicated!') + raise ValueError("Some elements of *xs* are duplicated!") x0, x1, x2 = xs y0, y1, y2 = ys - @vectorize(["float64(float64)"], target='parallel') + @vectorize(["float64(float64)"], target="parallel") def f(x: float) -> float: """ A helper function that only does the evaluation. @@ -90,9 +100,11 @@ def f(x: float) -> float: :param x: The variable on which the Lagrange polynomial is going to be applied. :return: The value of the Lagrange polynomial on :math:`x`, i.e., :math:`L(x)`. """ - return (x - x1) * (x - x2) / (x0 - x1) / (x0 - x2) * y0 + \ - (x - x0) * (x - x2) / (x1 - x0) / (x1 - x2) * y1 + \ - (x - x0) * (x - x1) / (x2 - x0) / (x2 - x1) * y2 + return ( + (x - x1) * (x - x2) / (x0 - x1) / (x0 - x2) * y0 + + (x - x0) * (x - x2) / (x1 - x0) / (x1 - x2) * y1 + + (x - x0) * (x - x1) / (x2 - x0) / (x2 - x1) * y2 + ) return f @@ -140,7 +152,7 @@ def find_nearest(array: Vector, value: Scalar) -> int: return j_low -@guvectorize([(float64[:], float64[:], int64[:])], '(m),(n)->(n)') +@guvectorize([(float64[:], float64[:], int64[:])], "(m),(n)->(n)") def vectorized_find_nearest(array: Vector, values: Vector, result: Vector): """ A vectorized version of function ``find_nearest``. @@ -153,10 +165,9 @@ def vectorized_find_nearest(array: Vector, values: Vector, result: Vector): n: int = len(array) if len(values) != len(result): - raise ValueError('The *values* and *result* arguments should have same length!') + raise ValueError("The *values* and *result* arguments should have same length!") for i in range(len(values)): - if values[i] <= array[0]: result[i] = 0 elif values[i] >= array[-1]: @@ -230,8 +241,11 @@ def is_monotonic_increasing(array: Vector) -> bool: return np.all(dx >= 0) -def calibrate_energy_on_reference(volumes_before_calibration: Matrix, energies_before_calibration: Matrix, - order: Optional[int] = 3): +def calibrate_energy_on_reference( + volumes_before_calibration: Matrix, + energies_before_calibration: Matrix, + order: Optional[int] = 3, +): """ In multi-configuration system calculation, the volume set of each calculation may vary a little, This function would make the volume set of the first configuration (normally, the most populated configuration) @@ -248,17 +262,22 @@ def calibrate_energy_on_reference(volumes_before_calibration: Matrix, energies_b energies_after_calibration = np.empty(volumes_before_calibration.shape) for i in range(configurations_amount): - strains_before_calibration = calculate_eulerian_strain(volumes_before_calibration[i, 0], - volumes_before_calibration[i]) - strains_after_calibration = calculate_eulerian_strain(volumes_before_calibration[i, 0], volumes_for_reference) - _, energies_after_calibration[i, :] = polynomial_least_square_fitting(strains_before_calibration, - energies_before_calibration[i], - strains_after_calibration, - order=order) + strains_before_calibration = calculate_eulerian_strain( + volumes_before_calibration[i, 0], volumes_before_calibration[i] + ) + strains_after_calibration = calculate_eulerian_strain( + volumes_before_calibration[i, 0], volumes_for_reference + ) + _, energies_after_calibration[i, :] = polynomial_least_square_fitting( + strains_before_calibration, + energies_before_calibration[i], + strains_after_calibration, + order=order, + ) return energies_after_calibration -if __name__ == '__main__': +if __name__ == "__main__": import doctest doctest.testmod() diff --git a/qha/type_aliases.py b/qha/type_aliases.py index cb83849..244001a 100644 --- a/qha/type_aliases.py +++ b/qha/type_aliases.py @@ -17,7 +17,7 @@ from numba import float64 # ===================== What can be exported? ===================== -__all__ = ['Scalar', 'Vector', 'Matrix', 'Array3D', 'Array4D'] +__all__ = ["Scalar", "Vector", "Matrix", "Array3D", "Array4D"] Scalar = float64 # 0-dimensional float Vector = float64[:] # 1-dimensional floats diff --git a/qha/unit_conversion.py b/qha/unit_conversion.py index 3ca5270..7c59ae2 100644 --- a/qha/unit_conversion.py +++ b/qha/unit_conversion.py @@ -13,62 +13,62 @@ # ===================== What can be exported? ===================== __all__ = [ - 'j_to_ev', - 'ev_to_j', - 'gpa_to_megabar', - 'megabar_to_gpa', - 'b3_to_a3', - 'a3_to_b3', - 'eh_to_ev', - 'ev_to_eh', - 'ry_to_ev', - 'ev_to_ry', - 'j_to_eh', - 'eh_to_j', - 'eh_to_hz', - 'hz_to_eh', - 'eh_to_k', - 'k_to_eh', - 'eh_to_m_inverse', - 'm_inverse_to_eh', - 'eh_to_cm_inverse', - 'cm_inverse_to_eh', - 'ev_to_m_inverse', - 'm_inverse_to_ev', - 'ev_to_cm_inverse', - 'cm_inverse_to_ev', - 'ev_to_k', - 'k_to_ev', - 'ry_to_j', - 'j_to_ry', - 'gpa_to_ev_a3', - 'ev_a3_to_gpa', - 'gpa_to_ry_b3', - 'ry_b3_to_gpa', - 'gpa_to_ha_b3', - 'ha_b3_to_gpa', - 'ev_b3_to_gpa', - 'gpa_to_ev_b3', - 'ry_b_to_ev_a', - 'ha_b_to_ev_a', - 'ry_to_kj_mol', - 'ry_to_j_mol' + "j_to_ev", + "ev_to_j", + "gpa_to_megabar", + "megabar_to_gpa", + "b3_to_a3", + "a3_to_b3", + "eh_to_ev", + "ev_to_eh", + "ry_to_ev", + "ev_to_ry", + "j_to_eh", + "eh_to_j", + "eh_to_hz", + "hz_to_eh", + "eh_to_k", + "k_to_eh", + "eh_to_m_inverse", + "m_inverse_to_eh", + "eh_to_cm_inverse", + "cm_inverse_to_eh", + "ev_to_m_inverse", + "m_inverse_to_ev", + "ev_to_cm_inverse", + "cm_inverse_to_ev", + "ev_to_k", + "k_to_ev", + "ry_to_j", + "j_to_ry", + "gpa_to_ev_a3", + "ev_a3_to_gpa", + "gpa_to_ry_b3", + "ry_b3_to_gpa", + "gpa_to_ha_b3", + "ha_b3_to_gpa", + "ev_b3_to_gpa", + "gpa_to_ev_b3", + "ry_b_to_ev_a", + "ha_b_to_ev_a", + "ry_to_kj_mol", + "ry_to_j_mol", ] # ===================== Constants ===================== -BOHR = physical_constants['Bohr radius'][0] -EH_EV = physical_constants['Hartree energy in eV'][0] -RY_EV = physical_constants['Rydberg constant times hc in eV'][0] -EH_J = physical_constants['hartree-joule relationship'][0] -EH_HZ = physical_constants['hartree-hertz relationship'][0] -EH_K = physical_constants['hartree-kelvin relationship'][0] -EH_M_INVERSE = physical_constants['hartree-inverse meter relationship'][0] -EV_M_INVERSE = physical_constants['electron volt-inverse meter relationship'][0] -EV_K = physical_constants['electron volt-kelvin relationship'][0] -RY_J = physical_constants['Rydberg constant times hc in J'][0] +BOHR = physical_constants["Bohr radius"][0] +EH_EV = physical_constants["Hartree energy in eV"][0] +RY_EV = physical_constants["Rydberg constant times hc in eV"][0] +EH_J = physical_constants["hartree-joule relationship"][0] +EH_HZ = physical_constants["hartree-hertz relationship"][0] +EH_K = physical_constants["hartree-kelvin relationship"][0] +EH_M_INVERSE = physical_constants["hartree-inverse meter relationship"][0] +EV_M_INVERSE = physical_constants["electron volt-inverse meter relationship"][0] +EV_K = physical_constants["electron volt-kelvin relationship"][0] +RY_J = physical_constants["Rydberg constant times hc in J"][0] # ===================== Settings ===================== -_target = DEFAULT_SETTINGS['target'] +_target = DEFAULT_SETTINGS["target"] # ===================== Functions ===================== @@ -388,7 +388,7 @@ def gpa_to_ev_a3(value): :param value: The value to be converted. :return: The converted value. """ - return value * 1e9 / electron_volt * angstrom ** 3 + return value * 1e9 / electron_volt * angstrom**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -399,7 +399,7 @@ def ev_a3_to_gpa(value): :param value: The value to be converted. :return: The converted value. """ - return value / 1e9 * electron_volt / angstrom ** 3 + return value / 1e9 * electron_volt / angstrom**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -410,7 +410,7 @@ def gpa_to_ev_b3(value): :param value: The value to be converted. :return: The converted value. """ - return value * 1e9 / electron_volt * BOHR ** 3 + return value * 1e9 / electron_volt * BOHR**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -421,7 +421,7 @@ def ev_b3_to_gpa(value): :param value: The value to be converted. :return: The converted value. """ - return value / 1e9 * electron_volt / BOHR ** 3 + return value / 1e9 * electron_volt / BOHR**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -432,7 +432,7 @@ def gpa_to_ry_b3(value): :param value: The value to be converted. :return: The converted value. """ - return value * 1e9 / RY_J * BOHR ** 3 + return value * 1e9 / RY_J * BOHR**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -443,7 +443,7 @@ def ry_b3_to_gpa(value): :param value: The value to be converted. :return: The converted value. """ - return value / 1e9 * RY_J / BOHR ** 3 + return value / 1e9 * RY_J / BOHR**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -454,7 +454,7 @@ def gpa_to_ha_b3(value): :param value: The value to be converted. :return: The converted value. """ - return value * 1e9 / RY_J * BOHR ** 3 / 2 + return value * 1e9 / RY_J * BOHR**3 / 2 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) @@ -465,7 +465,7 @@ def ha_b3_to_gpa(value): :param value: The value to be converted. :return: The converted value. """ - return 2 * value / 1e9 * RY_J / BOHR ** 3 + return 2 * value / 1e9 * RY_J / BOHR**3 @vectorize([float64(float64)], nopython=True, cache=True, target=_target) diff --git a/qha/v2p.py b/qha/v2p.py index a5f3c29..fae055f 100644 --- a/qha/v2p.py +++ b/qha/v2p.py @@ -14,7 +14,7 @@ from qha.type_aliases import Matrix, Vector # ===================== What can be exported? ===================== -__all__ = ['v2p'] +__all__ = ["v2p"] @numba.jit(nopython=True, parallel=True) @@ -34,11 +34,12 @@ def _lagrange4(x: float, x0, x1, x2, x3, y0, y1, y2, y3) -> float: :param y3: The y-coordinate of the 4th point. :return: The y-coordinate of the point to be evaluated. """ - return (x - x1) * (x - x2) * (x - x3) / (x0 - x1) / (x0 - x2) / (x0 - x3) * y0 + \ - (x - x0) * (x - x2) * (x - x3) / (x1 - x0) / (x1 - x2) / (x1 - x3) * y1 + \ - (x - x0) * (x - x1) * (x - x3) / (x2 - x0) / (x2 - x1) / (x2 - x3) * y2 + \ - (x - x0) * (x - x1) * (x - x2) / \ - (x3 - x0) / (x3 - x1) / (x3 - x2) * y3 + return ( + (x - x1) * (x - x2) * (x - x3) / (x0 - x1) / (x0 - x2) / (x0 - x3) * y0 + + (x - x0) * (x - x2) * (x - x3) / (x1 - x0) / (x1 - x2) / (x1 - x3) * y1 + + (x - x0) * (x - x1) * (x - x3) / (x2 - x0) / (x2 - x1) / (x2 - x3) * y2 + + (x - x0) * (x - x1) * (x - x2) / (x3 - x0) / (x3 - x1) / (x3 - x2) * y3 + ) def v2p(func_of_t_v: Matrix, p_of_t_v: Matrix, desired_pressures: Vector) -> Matrix: @@ -58,9 +59,15 @@ def v2p(func_of_t_v: Matrix, p_of_t_v: Matrix, desired_pressures: Vector) -> Mat result = np.empty((t_amount, desired_pressures_amount)) extended_f = np.hstack( - (func_of_t_v[:, 3].reshape(-1, 1), func_of_t_v, func_of_t_v[:, -4].reshape(-1, 1))) + ( + func_of_t_v[:, 3].reshape(-1, 1), + func_of_t_v, + func_of_t_v[:, -4].reshape(-1, 1), + ) + ) extended_p = np.hstack( - (p_of_t_v[:, 3].reshape(-1, 1), p_of_t_v, p_of_t_v[:, -4].reshape(-1, 1))) + (p_of_t_v[:, 3].reshape(-1, 1), p_of_t_v, p_of_t_v[:, -4].reshape(-1, 1)) + ) for i in range(t_amount): rs = np.zeros(desired_pressures_amount) @@ -69,10 +76,11 @@ def v2p(func_of_t_v: Matrix, p_of_t_v: Matrix, desired_pressures: Vector) -> Mat for j in range(desired_pressures_amount): k = int(rs[j]) # Interpolate pressures around desired pressure, x1-x4 in `_lagrange4` - x1, x2, x3, x4 = extended_p[i, k - 1:k + 3] + x1, x2, x3, x4 = extended_p[i, k - 1 : k + 3] # Interpolate func, f1-f4 in `_lagrange4` - f1, f2, f3, f4 = extended_f[i, k - 1:k + 3] + f1, f2, f3, f4 = extended_f[i, k - 1 : k + 3] # Evaluate func at desired pressure result[i, j] = _lagrange4( - desired_pressures[j], x1, x2, x3, x4, f1, f2, f3, f4) + desired_pressures[j], x1, x2, x3, x4, f1, f2, f3, f4 + ) return result From 26c00d2cc9d67c9388876393fde1b8ffebd9fcca Mon Sep 17 00:00:00 2001 From: singularitti Date: Sat, 21 Oct 2023 03:53:57 -0400 Subject: [PATCH 2/2] Replace `numba.vectorize` with `numpy.vectorize` --- qha/grid_interpolation.py | 4 +- qha/statmech.py | 6 +-- qha/tools.py | 4 +- qha/unit_conversion.py | 80 +++++++++++++++++++-------------------- 4 files changed, 47 insertions(+), 47 deletions(-) diff --git a/qha/grid_interpolation.py b/qha/grid_interpolation.py index 70787c7..9e410cc 100644 --- a/qha/grid_interpolation.py +++ b/qha/grid_interpolation.py @@ -26,7 +26,7 @@ ] -@vectorize([float64(float64, float64)], nopython=True) +@np.vectorize def calculate_eulerian_strain(v0, vs): """ Calculate the Eulerian strain (:math:`f`s) of a given volume vector *vs* with respect to a reference volume *v0*, @@ -43,7 +43,7 @@ def calculate_eulerian_strain(v0, vs): return 1 / 2 * ((v0 / vs) ** (2 / 3) - 1) -@vectorize([float64(float64, float64)], nopython=True) +@np.vectorize def from_eulerian_strain(v0, fs): """ Calculate the corresponding volumes :math:`V`s from a vector of given Eulerian strains (*fs*) diff --git a/qha/statmech.py b/qha/statmech.py index f42e97b..073112b 100644 --- a/qha/statmech.py +++ b/qha/statmech.py @@ -41,7 +41,7 @@ }[qha.settings.energy_unit] -@vectorize([float64(float64, float64)], nopython=True, target="parallel", cache=True) +@np.vectorize def ho_free_energy(temperature, frequency): """ Calculate Helmholtz free energy of a single harmonic oscillator at a specific temperature. @@ -60,7 +60,7 @@ def ho_free_energy(temperature, frequency): return 1 / 2 * hw + kt * np.log(1 - np.exp(-hw / kt)) -@vectorize([float64(float64, float64)], nopython=True, target="parallel", cache=True) +@np.vectorize def subsystem_partition_function(temperature, frequency): """ Calculate the subsystem partition function of a single harmonic oscillator at a specific temperature. @@ -78,7 +78,7 @@ def subsystem_partition_function(temperature, frequency): return np.exp(x / 2) / (1 - np.exp(x)) -@vectorize([float64(float64, float64)], nopython=True, target="parallel", cache=True) +@np.vectorize def log_subsystem_partition_function(temperature, frequency): """ Calculate the natural logarithm of the subsystem partition function of a single harmonic oscillator diff --git a/qha/tools.py b/qha/tools.py index 017caea..4711787 100644 --- a/qha/tools.py +++ b/qha/tools.py @@ -46,7 +46,7 @@ def lagrange4(xs: Vector, ys: Vector) -> Callable[[float], float]: x0, x1, x2, x3 = xs y0, y1, y2, y3 = ys - @vectorize(["float64(float64)"], target="parallel") + @np.vectorize def f(x: float) -> float: """ A helper function that only does the evaluation. @@ -92,7 +92,7 @@ def lagrange3(xs: Vector, ys: Vector) -> Callable[[float], float]: x0, x1, x2 = xs y0, y1, y2 = ys - @vectorize(["float64(float64)"], target="parallel") + @np.vectorize def f(x: float) -> float: """ A helper function that only does the evaluation. diff --git a/qha/unit_conversion.py b/qha/unit_conversion.py index 7c59ae2..36c4b2e 100644 --- a/qha/unit_conversion.py +++ b/qha/unit_conversion.py @@ -72,7 +72,7 @@ # ===================== Functions ===================== -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def j_to_ev(value): """ Convert the *value* in unit joule to electronvolt. @@ -83,7 +83,7 @@ def j_to_ev(value): return value / electron_volt -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_to_j(value): """ Convert the *value* in unit electronvolt to joule. @@ -94,7 +94,7 @@ def ev_to_j(value): return value * electron_volt -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def gpa_to_megabar(value): """ Convert the *value* in unit gigapascal to megabar. @@ -105,7 +105,7 @@ def gpa_to_megabar(value): return value * 0.01 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def megabar_to_gpa(value): """ Convert the *value* in unit megabar to gigapascal. @@ -116,7 +116,7 @@ def megabar_to_gpa(value): return value / 0.01 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def b3_to_a3(value): """ Convert the *value* in unit cubic bohr radius to what in cubic angstrom. @@ -127,7 +127,7 @@ def b3_to_a3(value): return value * (BOHR / angstrom) ** 3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def a3_to_b3(value): """ Convert the *value* in unit cubic angstrom to what in cubic bohr radius. @@ -138,7 +138,7 @@ def a3_to_b3(value): return value * (angstrom / BOHR) ** 3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def eh_to_ev(value): """ Convert the *value* in unit hartree to electronvolt. @@ -149,7 +149,7 @@ def eh_to_ev(value): return value * EH_EV -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_to_eh(value): """ Convert the *value* in unit electronvolt to hartree. @@ -160,7 +160,7 @@ def ev_to_eh(value): return value / EH_EV -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ry_to_ev(value): """ Convert the *value* in unit rydberg to electronvolt. @@ -171,7 +171,7 @@ def ry_to_ev(value): return value * RY_EV -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_to_ry(value): """ Convert the *value* in unit electronvolt to rydberg. @@ -182,7 +182,7 @@ def ev_to_ry(value): return value / RY_EV -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def j_to_eh(value): """ Convert the *value* in unit joule to hartree. @@ -193,7 +193,7 @@ def j_to_eh(value): return value / EH_J -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def eh_to_j(value): """ Convert the *value* in unit hartree to joule. @@ -204,7 +204,7 @@ def eh_to_j(value): return value * EH_J -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def eh_to_hz(value): """ Convert the *value* in unit hartree to hertz. @@ -215,7 +215,7 @@ def eh_to_hz(value): return value * EH_HZ -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def hz_to_eh(value): """ Convert the *value* in unit hertz to hartree. @@ -226,7 +226,7 @@ def hz_to_eh(value): return value / EH_HZ -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def eh_to_k(value): """ Convert the *value* in unit hartree to kelvin. @@ -237,7 +237,7 @@ def eh_to_k(value): return value * EH_K -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def k_to_eh(value): """ Convert the *value* in unit kelvin to hartree. @@ -248,7 +248,7 @@ def k_to_eh(value): return value / EH_K -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def eh_to_m_inverse(value): """ Convert the *value* in unit hartree to :math:`\\text{m}^{-1}`. @@ -259,7 +259,7 @@ def eh_to_m_inverse(value): return value * EH_M_INVERSE -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def m_inverse_to_eh(value): """ Convert the *value* in unit :math:`\\text{m}^{-1}` to hartree. @@ -270,7 +270,7 @@ def m_inverse_to_eh(value): return value / EH_M_INVERSE -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def eh_to_cm_inverse(value): """ Convert the *value* in unit hartree to :math:`\\text{cm}^{-1}`. @@ -281,7 +281,7 @@ def eh_to_cm_inverse(value): return value * EH_M_INVERSE / 100 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def cm_inverse_to_eh(value): """ Convert the *value* in unit :math:`\\text{cm}^{-1}` to hartree. @@ -292,7 +292,7 @@ def cm_inverse_to_eh(value): return value / EH_M_INVERSE * 100 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_to_m_inverse(value): """ Convert the *value* in unit electronvolt to :math:`\\text{m}^{-1}`. @@ -303,7 +303,7 @@ def ev_to_m_inverse(value): return value * EV_M_INVERSE -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def m_inverse_to_ev(value): """ Convert the *value* in unit :math:`\\text{m}^{-1}` to electronvolt. @@ -314,7 +314,7 @@ def m_inverse_to_ev(value): return value / EV_M_INVERSE -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_to_cm_inverse(value): """ Convert the *value* in unit electronvolt to :math:`\\text{cm}^{-1}`. @@ -325,7 +325,7 @@ def ev_to_cm_inverse(value): return value * EV_M_INVERSE / 100 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def cm_inverse_to_ev(value): """ Convert the *value* in unit :math:`\\text{cm}^{-1}` to electronvolt. @@ -336,7 +336,7 @@ def cm_inverse_to_ev(value): return value / EV_M_INVERSE * 100 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_to_k(value): """ Convert the *value* in unit electronvolt to kelvin. @@ -347,7 +347,7 @@ def ev_to_k(value): return value * EV_K -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def k_to_ev(value): """ Convert the *value* in unit kelvin to electronvolt. @@ -358,7 +358,7 @@ def k_to_ev(value): return value / EV_K -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ry_to_j(value): """ Convert the *value* in unit rydberg to joule. @@ -369,7 +369,7 @@ def ry_to_j(value): return value * RY_J -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def j_to_ry(value): """ Convert the *value* in unit joule to rydberg. @@ -380,7 +380,7 @@ def j_to_ry(value): return value / RY_J -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def gpa_to_ev_a3(value): """ Convert the *value* in unit gigapascal to :math:`\\frac{ \\text{electronvolt} }{ \\text{angstrom}^3 }`. @@ -391,7 +391,7 @@ def gpa_to_ev_a3(value): return value * 1e9 / electron_volt * angstrom**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_a3_to_gpa(value): """ Convert the *value* in unit :math:`\\frac{ \\text{electronvolt} }{ \\text{angstrom}^3 }` to gigapascal. @@ -402,7 +402,7 @@ def ev_a3_to_gpa(value): return value / 1e9 * electron_volt / angstrom**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def gpa_to_ev_b3(value): """ Convert the *value* in unit gigapascal to :math:`\\frac{ \\text{electronvolt} }{ \\text{bohr radius}^3 }`. @@ -413,7 +413,7 @@ def gpa_to_ev_b3(value): return value * 1e9 / electron_volt * BOHR**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ev_b3_to_gpa(value): """ Convert the *value* in unit :math:`\\frac{ \\text{electronvolt} }{ \\text{bohr radius}^3 }` to gigapascal. @@ -424,7 +424,7 @@ def ev_b3_to_gpa(value): return value / 1e9 * electron_volt / BOHR**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def gpa_to_ry_b3(value): """ Convert the *value* in unit gigapascal to :math:`\\frac{ \\text{rydberg} }{ \\text{bohr radius}^3 }`. @@ -435,7 +435,7 @@ def gpa_to_ry_b3(value): return value * 1e9 / RY_J * BOHR**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ry_b3_to_gpa(value): """ Convert the *value* in unit :math:`\\frac{ \\text{rydberg} }{ \\text{bohr radius}^3 }` to gigapascal. @@ -446,7 +446,7 @@ def ry_b3_to_gpa(value): return value / 1e9 * RY_J / BOHR**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def gpa_to_ha_b3(value): """ Convert the *value* in unit gigapascal to :math:`\\frac{ \\text{hartree} }{ \\text{bohr radius}^3 }`. @@ -457,7 +457,7 @@ def gpa_to_ha_b3(value): return value * 1e9 / RY_J * BOHR**3 / 2 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ha_b3_to_gpa(value): """ Convert the *value* in unit :math:`\\frac{ \\text{hartree} }{ \\text{bohr radius}^3 }` to gigapascal. @@ -468,7 +468,7 @@ def ha_b3_to_gpa(value): return 2 * value / 1e9 * RY_J / BOHR**3 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ry_b_to_ev_a(value): """ The atomic force. @@ -481,7 +481,7 @@ def ry_b_to_ev_a(value): return value * RY_EV / (BOHR / angstrom) -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ha_b_to_ev_a(value): """ The atomic force. @@ -494,7 +494,7 @@ def ha_b_to_ev_a(value): return value * EH_EV / (BOHR / angstrom) -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ry_to_kj_mol(value): """ Convert the *value* is in unit Rydberg, the converted value is in unit @@ -506,7 +506,7 @@ def ry_to_kj_mol(value): return value * Avogadro * RY_J / 1000 -@vectorize([float64(float64)], nopython=True, cache=True, target=_target) +@np.vectorize def ry_to_j_mol(value): """ Convert the *value* is in unit Rydberg, the converted value is in unit