Merge branch 'hotfix/1.0.17'

MineralsCloud · Jun 13, 2019 · 28f9f0a · 28f9f0a
2 parents 973ccff + 324c7d4
commit 28f9f0a
Show file tree

Hide file tree

Showing 8 changed files with 129 additions and 28 deletions.
diff --git a/qha/multi_configurations/different_phonon_dos.py b/qha/multi_configurations/different_phonon_dos.py
@@ -34,7 +34,7 @@ class PartitionFunction:
 
     .. math::
 
-       Z_{\\text{all configs}}(T, V) = \sum_{j = 1}^{N_{c}} g_{j} Z_{j}(T, V),
+       Z_{\\text{all configs}}(T, V) = \\sum_{j = 1}^{N_{c}} g_{j} Z_{j}(T, V),
 
     where :math:`N_{c}` stands for the number of configurations and :math:`g_{j}` stands for degeneracy for the
     :math:`j` th configuration.
@@ -113,7 +113,7 @@ def partition_functions_for_each_configuration(self):
 
         .. math::
 
-           Z_{j}(T, V) = \exp \\bigg( -\\frac{ F_{j}(T, V) }{ k_B T } \\bigg).
+           Z_{j}(T, V) = \\exp \\bigg( -\\frac{ F_{j}(T, V) }{ k_B T } \\bigg).
 
         :return: A matrix, the partition function of each configuration of each volume.
         """
@@ -123,7 +123,29 @@ def partition_functions_for_each_configuration(self):
             raise ImportError("Install ``bigfloat`` package to use {0} object!".format(self.__class__.__name__))
 
         with bigfloat.precision(self.precision):
-            return np.array([bigfloat.exp(d) for d in  # shape = (# of volumes for each configuration, 1)
+            # shape = (# of configurations, # of volumes for each configuration)
+            exp = np.vectorize(bigfloat.exp)
+            return exp(-self.aligned_free_energies_for_each_configuration / (K * self.temperature))
+
+    @LazyProperty
+    def partition_functions_for_all_configurations(self):
+        """
+        Sum the partition functions for all configurations.
+
+        .. math::
+
+           Z_{\\text{all configs}}(T, V) = \\sum_{j} Z_{j}(T, V).
+
+        :return: A vector, the partition function of each volume.
+        """
+        try:
+            import bigfloat
+        except ImportError:
+            raise ImportError("Install ``bigfloat`` package to use {0} object!".format(self.__class__.__name__))
+
+        with bigfloat.precision(self.precision):
+            # shape = (# of volumes,)
+            return np.array([bigfloat.exp(d) for d in
                              logsumexp(-self.aligned_free_energies_for_each_configuration.T / (K * self.temperature),
                                        axis=1, b=self.degeneracies)])
 
@@ -133,7 +155,7 @@ def get_free_energies(self):
 
         .. math::
 
-           F_{\\text{all configs}}(T, V) = - k_B T \ln Z_{\\text{all configs}}(T, V).
+           F_{\\text{all configs}}(T, V) = - k_B T \\ln Z_{\\text{all configs}}(T, V).
 
         :return: The free energy on a temperature-volume grid.
         """
@@ -143,5 +165,5 @@ def get_free_energies(self):
             raise ImportError("Install ``bigfloat`` package to use {0} object!".format(self.__class__.__name__))
 
         with bigfloat.precision(self.precision):
-            log_z = np.array([bigfloat.log(d) for d in self.partition_functions_for_each_configuration], dtype=float)
+            log_z = np.array([bigfloat.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
@@ -33,13 +33,13 @@ class PartitionFunction:
 
     .. math::
 
-       Z_{\\text{all configs}}(T, V) = \sum_{j = 1}^{N_{c}} g_{j}
-       \\bigg\{
-         \exp \\Big( -\\frac{ E_j(V) }{ k_B T } \\Big)
-         \prod_{\mathbf{q}s} \\bigg(
-           \\tfrac{\exp \\big(-\\tfrac{ \hbar \omega_{\mathbf{ q }s}^j(V) }{ 2 k_B T }\\big)}{1 - \exp \\big(-\\tfrac{ \hbar \omega_{\mathbf{ q }s}^j(V) }{ k_B T }\\big)}
-         \\bigg)^{w_{\mathbf{ q }}^j}
-       \\bigg\}.
+       Z_{\\text{all configs}}(T, V) = \\sum_{j = 1}^{N_{c}} g_{j}
+       \\bigg\\{
+         \\exp \\Big( -\\frac{ E_j(V) }{ k_B T } \\Big)
+         \\prod_{\\mathbf{q}s} \\bigg(
+           \\tfrac{\\exp \\big(-\\tfrac{ \\hbar \\omega_{\\mathbf{ q }s}^j(V) }{ 2 k_B T }\\big)}{1 - \\exp \\big(-\\tfrac{ \\hbar \\omega_{\\mathbf{ q }s}^j(V) }{ k_B T }\\big)}
+         \\bigg)^{w_{\\mathbf{ q }}^j}
+       \\bigg\\}.
 
     :param temperature: The temperature at which the partition function is calculated.
     :param degeneracies: An array of degeneracies of each configuration, which will not be normalized in the
@@ -119,7 +119,7 @@ def get_free_energies(self):
 
         .. math::
 
-           F_{\\text{all configs}}(T, V) = - k_B T \ln Z_{\\text{all configs}}(T, V).
+           F_{\\text{all configs}}(T, V) = - k_B T \\ln Z_{\\text{all configs}}(T, V).
 
         :return: The free energy on the temperature-volume grid.
         """
@@ -140,12 +140,12 @@ class FreeEnergy:
 
     .. math::
 
-       F_{\\text{all configs}}(T, V) = - k_B T \ln Z_{\\text{all configs}}(T, V)
-       = - k_B T \ln \\bigg( \sum_{j = 1}^{N_{c}} g_{j} \exp \\Big( -\\frac{ E_j(V) }{ k_B T } \\Big) \\bigg)
-         + \sum_{\mathbf{ q }s} w_\mathbf{ q }
-             \\bigg\{ \\frac{ \hbar \omega_{\mathbf{ q }s}(V) }{ 2 }
-             + k_B \ln \\bigg( 1 - \exp \\Big( -\\frac{ \hbar \omega_{\mathbf{ q }s}(V) }{ k_B T } \\Big) \\bigg)
-         \\bigg\}.
+       F_{\\text{all configs}}(T, V) = - k_B T \\ln Z_{\\text{all configs}}(T, V)
+       = - k_B T \\ln \\bigg( \\sum_{j = 1}^{N_{c}} g_{j} \\exp \\Big( -\\frac{ E_j(V) }{ k_B T } \\Big) \\bigg)
+         + \\sum_{\\mathbf{ q }s} w_\\mathbf{ q }
+             \\bigg\\{ \\frac{ \\hbar \\omega_{\\mathbf{ q }s}(V) }{ 2 }
+             + k_B \\ln \\bigg( 1 - \\exp \\Big( -\\frac{ \\hbar \\omega_{\\mathbf{ q }s}(V) }{ k_B T } \\Big) \\bigg)
+         \\bigg\\}.
 
     :param temperature: The temperature at which the partition function is calculated.
     :param degeneracies: An array of degeneracies of each configuration, which will not be normalized in the

diff --git a/qha/settings.py b/qha/settings.py
@@ -80,7 +80,7 @@ def from_yaml(filename: str) -> Settings:
     :return: A ``Settings`` class.
     """
     with open(filename, 'r') as f:
-        return Settings(yaml.load(f))
+        return Settings(yaml.load(f, Loader=yaml.CLoader))
 
 
 global energy_unit

diff --git a/qha/tests/test_different_phonon_dos.py b/qha/tests/test_different_phonon_dos.py
@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+
+import os
+import unittest
+
+import numpy as np
+
+from qha.calculator import DifferentPhDOSCalculator
+from qha.multi_configurations.different_phonon_dos import PartitionFunction
+from qha.settings import from_yaml
+
+
+class TestPartitionFunction(unittest.TestCase):
+    def setUp(self) -> None:
+        os.chdir('../../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'):
+            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)
+
+    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})
+
+    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._volumes.shape, (52, 6))
+        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])
diff --git a/qha/tests/test_read_input.py b/qha/tests/test_read_input.py
@@ -23,7 +23,7 @@ def test_read_si_input(self):
 
     def test_read_ice_input(self):
         for i in range(1, 53):
-            file_path = self.dir / "{0}{1:02d}".format('ice VII/input_conf', 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)
 

diff --git a/qha/tests/test_single_configuration.py b/qha/tests/test_single_configuration.py
@@ -3,15 +3,15 @@
 
 import unittest
 
-from qha.single_configuration import ho_free_energy
+from qha.statmech import ho_free_energy
 
 
 class TestSingleConfiguration(unittest.TestCase):
     def test_ho_free_energy(self):
         self.assertEqual(ho_free_energy(0, 0), 0)
         self.assertEqual(ho_free_energy(1, -2), 0)
-        self.assertEqual(ho_free_energy(0, 1000), 0.004556299262079407)
-        self.assertEqual(ho_free_energy(100, 1000), 0.0045562989049199466)
+        self.assertAlmostEqual(ho_free_energy(0, 1000), 0.004556299262079407)
+        self.assertAlmostEqual(ho_free_energy(100, 1000), 0.0045562989049199466)
 
 
 if __name__ == '__main__':

diff --git a/qha/tools.py b/qha/tools.py
@@ -24,7 +24,7 @@
 def lagrange4(xs: Vector, ys: Vector) -> Callable[[float], float]:
     """
     A third-order Lagrange polynomial function. Given 4 points for interpolation:
-    :math:`(x_0, y_0), \ldots, (x_3, y_3)`, evaluate the Lagrange polynomial on :math:`x`, referenced from
+    :math:`(x_0, y_0), \\ldots, (x_3, y_3)`, evaluate the Lagrange polynomial on :math:`x`, referenced from
     `Wolfram MathWorld <http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html>`_.
 
     :param xs: A vector of the x-coordinates' of the 4 points.
@@ -57,7 +57,7 @@ def f(x: float) -> float:
 def lagrange3(xs: Vector, ys: Vector) -> Callable[[float], float]:
     """
     A second-order Lagrange polynomial function. Given 3 points for interpolation:
-    :math:`(x_0, y_0), \ldots, (x_2, y_2)`, evaluate the Lagrange polynomial on :math:`x`, referenced from
+    :math:`(x_0, y_0), \\ldots, (x_2, y_2)`, evaluate the Lagrange polynomial on :math:`x`, referenced from
     `Wolfram MathWorld also <http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html>`_.
 
     .. doctest::
@@ -103,7 +103,7 @@ def find_nearest(array: Vector, value: Scalar) -> int:
     Given an *array* , and given a *value* , returns an index ``j`` such that *value* is between ``array[j]``
     and ``array[j+1]``. The *array* must be monotonic increasing. ``j=-1`` or ``j=len(array)`` is returned
     to indicate that *value* is out of range below and above respectively.
-    If *array* is unsorted, consider first using an :math:`O(n \log n)` sort and then use this function.
+    If *array* is unsorted, consider first using an :math:`O(n \\log n)` sort and then use this function.
     Referenced from `Stack Overflow <https://stackoverflow.com/questions/2566412/find-nearest-value-in-numpy-array>`_.
 
     :param array: An array of monotonic increasing real numbers.

diff --git a/qha/v2p.py b/qha/v2p.py
@@ -21,7 +21,7 @@
 def _lagrange4(x: float, x0, x1, x2, x3, y0, y1, y2, y3) -> float:
     """
     A third-order Lagrange polynomial function. Given 4 points for interpolation:
-    :math:`(x_0, y_0), \ldots, (x_3, y_3)`, evaluate the Lagrange polynomial on :math:`x`.
+    :math:`(x_0, y_0), \\ldots, (x_3, y_3)`, evaluate the Lagrange polynomial on :math:`x`.
 
     :param x: The x-coordinate of the point to be evaluated.
     :param x0: The x-coordinate of the 1st point.