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+|(?<!\s)(?=-)', line)))
+                freqs = list(filter(lambda s: s, re.split(r"\s+|(?<!\s)(?=-)", line)))
                 x = np.hstack((x, freqs))
 
             frequencies.append(x)
@@ -203,12 +217,16 @@ def read_frequency_file(inp: str) -> 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