Merge pull request #65 from MineralsCloud/pull-64-restyled

Restyle feat: Add entropy output

Author: singularitti

qha/calculator.py

@@ -26,7 +26,7 @@
 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
+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? =====================
@@ -48,7 +48,8 @@ def __init__(self, user_settings: Dict[str, Any]):
             try:
                 runtime_settings.update({key: user_settings[key]})
             except KeyError:
-                continue  # If a key is not set in user settings, use the default.
+                # If a key is not set in user settings, use the default.
+                continue
 
         self._settings = runtime_settings
 
@@ -106,12 +107,15 @@ 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!")
+            raise RuntimeError(
+                "Check the input file to make sure the volume decreases!")
 
         self._formula_unit_number: int = formula_unit_number
         self._volumes = volumes
@@ -127,7 +131,8 @@ def where_negative_frequencies(self) -> Optional[Vector]:
         :return:
         """
         if self._frequencies is None:
-            print("Please invoke ``read_input`` method first!")  # ``None`` is returned
+            # ``None`` is returned
+            print("Please invoke ``read_input`` method first!")
         else:
             _ = np.transpose(np.where(self._frequencies < 0))
             if _.size == 0:
@@ -156,20 +161,23 @@ def refine_grid(self):
         try:
             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!")
+            raise KeyError(
+                "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'])
         else:
-            self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid(self.volumes, self.vib_ry)
+            self._finer_volumes_bohr3, self._f_tv_ry, self._v_ratio = r.refine_grid(
+                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):
@@ -218,7 +226,8 @@ def desired_pressure_status(self) -> None:
                 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'])
+            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("""\
                            !!!ATTENTION!!!
@@ -228,12 +237,17 @@ def desired_pressure_status(self) -> None:
                            Please reduce the NTV accordingly, for example, try to set NTV < {:4d}.
                            """.format(ntv_max)))
 
-            raise ValueError("DESIRED PRESSURE is too high (NTV is too large), qha results might not be right!")
+            raise ValueError(
+                "DESIRED PRESSURE is too high (NTV is too large), qha results might not be right!")
 
     @LazyProperty
     def p_tv_au(self):
         return pressure(self.finer_volumes_bohr3, self.f_tv_ry)
 
+    @LazyProperty
+    def s_tv_j(self):
+        return ry_to_j(entropy(self.temperature_array, self.f_tv_ry))
+
     @LazyProperty
     def f_tv_ev(self):
         return ry_to_ev(self.f_tv_ry)
@@ -373,12 +387,15 @@ 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")
-                raise ValueError("Check the input file to make sure the volumes are monotonic decreasing!")
+                save_to_output(
+                    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!")
 
             formula_unit_numbers.append(nm_tmp)
             volumes.append(volumes_tmp)
@@ -393,12 +410,15 @@ def read_input(self):
         q_weights = np.array(q_weights)
 
         if not len(set(formula_unit_numbers)) == 1:
-            raise RuntimeError("All the formula unit number in all inputs should be the same!")
+            raise RuntimeError(
+                "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!")
 
-        self._formula_unit_number = formula_unit_numbers[0]  # Choose any of them since they are all the same
+        # Choose any of them since they are all the same
+        self._formula_unit_number = formula_unit_numbers[0]
         self._volumes = volumes
         self._static_energies = static_energies
         self._frequencies = frequencies
@@ -408,11 +428,12 @@ def read_input(self):
     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']
+            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()
+            mat[i] = different_phonon_dos.PartitionFunction(
+                t, *(arg for arg in args)).get_free_energies()
 
         return mat
 
@@ -424,9 +445,10 @@ 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']
+            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()
+            mat[i] = same_phonon_dos.FreeEnergy(
+                t, *(arg for arg in args)).get_free_energies()
         return mat

qha/cli/runner.py

@@ -49,7 +49,8 @@ def run(self, namespace):
         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'])
@@ -64,10 +65,12 @@ def run(self, namespace):
             print("You have single-configuration calculation assumed.")
         elif calculation_type == 'same phonon dos':
             calc = SamePhDOSCalculator(user_settings)
-            print("You have multi-configuration calculation with the same phonon DOS assumed.")
+            print(
+                "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.")
+            print(
+                "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'!")
@@ -81,7 +84,8 @@ def run(self, namespace):
 
         print("Caution: If negative frequencies found, they are currently treated as 0!")
         tmp = calc.where_negative_frequencies
-        if tmp is not None and not (tmp.T[-1].max() <= 2):  # Don't delete this parenthesis!
+        # Don't delete this parenthesis!
+        if tmp is not None and not (tmp.T[-1].max() <= 2):
             if calc.frequencies.ndim == 4:  # Multiple configuration
                 for indices in tmp:
                     print(
@@ -122,20 +126,29 @@ def run(self, namespace):
                               }
 
         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)
+        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)
+        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)
+        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']
+                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)
+                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'
@@ -145,15 +158,19 @@ def run(self, namespace):
                 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)
+                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_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))