From c3f754d4073241bb4516e654a99f4a5cae8a8d47 Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Wed, 12 Jun 2019 22:52:21 -0400 Subject: [PATCH 01/10] Fix a bug mentioned in issue #66 1. Deprecate the summation over configurations in `partition_functions_for_each_configuration` 2. Add an API `partition_functions_for_all_configurations` --- .../different_phonon_dos.py | 25 +++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/qha/multi_configurations/different_phonon_dos.py b/qha/multi_configurations/different_phonon_dos.py index 3aef618..0515708 100644 --- a/qha/multi_configurations/different_phonon_dos.py +++ b/qha/multi_configurations/different_phonon_dos.py @@ -123,7 +123,28 @@ 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) + return np.array( + [bigfloat.exp(d) for d in -self.aligned_free_energies_for_each_configuration / (K * self.temperature)]) + + 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)]) @@ -143,5 +164,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 From 18f4563e7bcf9b599c6ce14a7f52b08059c81617 Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Wed, 12 Jun 2019 22:53:21 -0400 Subject: [PATCH 02/10] doc: Fix invalid escape sequences --- qha/multi_configurations/different_phonon_dos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/qha/multi_configurations/different_phonon_dos.py b/qha/multi_configurations/different_phonon_dos.py index 0515708..1d0324e 100644 --- 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. """ @@ -154,7 +154,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. """ From 380a761eb0ee09d651cfd6e861c6f840e4ae298a Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Wed, 12 Jun 2019 23:14:05 -0400 Subject: [PATCH 03/10] test: Fix a bug mentioned in issue #67 --- qha/tests/test_read_input.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qha/tests/test_read_input.py b/qha/tests/test_read_input.py index 3016c17..95485e9 100644 --- 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) From a75cab0c44c0f33e9da1450e5e3279d655064d4d Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Wed, 12 Jun 2019 23:24:28 -0400 Subject: [PATCH 04/10] test: Fix a bug mentioned in issue #68 --- qha/tests/test_single_configuration.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/qha/tests/test_single_configuration.py b/qha/tests/test_single_configuration.py index f77ae9c..4c678d3 100644 --- a/qha/tests/test_single_configuration.py +++ b/qha/tests/test_single_configuration.py @@ -10,8 +10,8 @@ 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__': From abc7c67884ded9fcfe56fcdb227f5930a0a3ea6d Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Wed, 12 Jun 2019 23:28:34 -0400 Subject: [PATCH 05/10] Fix an import --- qha/tests/test_single_configuration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qha/tests/test_single_configuration.py b/qha/tests/test_single_configuration.py index 4c678d3..bc6abdc 100644 --- a/qha/tests/test_single_configuration.py +++ b/qha/tests/test_single_configuration.py @@ -3,7 +3,7 @@ import unittest -from qha.single_configuration import ho_free_energy +from qha.statmech import ho_free_energy class TestSingleConfiguration(unittest.TestCase): From 1458aec412cf91913cd98abe297c3433b0a8e5be Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Wed, 12 Jun 2019 23:34:38 -0400 Subject: [PATCH 06/10] docs: Fix invalid escape sequences --- qha/multi_configurations/same_phonon_dos.py | 28 ++++++++++----------- qha/tests/test_different_phonon_dos.py | 0 qha/tools.py | 6 ++--- qha/v2p.py | 2 +- 4 files changed, 18 insertions(+), 18 deletions(-) create mode 100644 qha/tests/test_different_phonon_dos.py diff --git a/qha/multi_configurations/same_phonon_dos.py b/qha/multi_configurations/same_phonon_dos.py index 4ba188d..373f0ce 100644 --- 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/tests/test_different_phonon_dos.py b/qha/tests/test_different_phonon_dos.py new file mode 100644 index 0000000..e69de29 diff --git a/qha/tools.py b/qha/tools.py index 315c405..1e1e377 100644 --- 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 `_. :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 `_. .. 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 `_. :param array: An array of monotonic increasing real numbers. diff --git a/qha/v2p.py b/qha/v2p.py index 460e554..7002417 100644 --- 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. From 0f902476f620ddbd1abbbbf48101a25c2dbb3dc9 Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Thu, 13 Jun 2019 00:40:10 -0400 Subject: [PATCH 07/10] Fix a warning mentioned in issue #69 --- qha/settings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qha/settings.py b/qha/settings.py index 1980035..093fbdc 100644 --- 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 From 76d9931bb2416294e1f267bd5469ffa7a58b28c2 Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Thu, 13 Jun 2019 02:15:19 -0400 Subject: [PATCH 08/10] test: Add tests for `multi_configurations.different_phonon_dos` --- qha/tests/test_different_phonon_dos.py | 79 ++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/qha/tests/test_different_phonon_dos.py b/qha/tests/test_different_phonon_dos.py index e69de29..5230d77 100644 --- 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]) From 38d1cdf348a1b4f735b57340ad9f1a24d8c5c925 Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Thu, 13 Jun 2019 02:15:42 -0400 Subject: [PATCH 09/10] Simplify code and fix a bug in `multi_configurations.different_phonon_dos` --- qha/multi_configurations/different_phonon_dos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/qha/multi_configurations/different_phonon_dos.py b/qha/multi_configurations/different_phonon_dos.py index 1d0324e..ce93d7f 100644 --- a/qha/multi_configurations/different_phonon_dos.py +++ b/qha/multi_configurations/different_phonon_dos.py @@ -124,8 +124,8 @@ def partition_functions_for_each_configuration(self): with bigfloat.precision(self.precision): # shape = (# of configurations, # of volumes for each configuration) - return np.array( - [bigfloat.exp(d) for d in -self.aligned_free_energies_for_each_configuration / (K * self.temperature)]) + exp = np.vectorize(bigfloat.exp) + return exp(-self.aligned_free_energies_for_each_configuration / (K * self.temperature)) def partition_functions_for_all_configurations(self): """ @@ -164,5 +164,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_all_configurations], 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 From 324c7d401de73bf38dc9bc68f2023718d6a2bee1 Mon Sep 17 00:00:00 2001 From: Qi Zhang Date: Thu, 13 Jun 2019 02:17:04 -0400 Subject: [PATCH 10/10] Make `partition_functions_for_all_configurations` a `LazyProperty` --- qha/multi_configurations/different_phonon_dos.py | 3 ++- qha/tests/test_different_phonon_dos.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/qha/multi_configurations/different_phonon_dos.py b/qha/multi_configurations/different_phonon_dos.py index ce93d7f..6aedfb1 100644 --- a/qha/multi_configurations/different_phonon_dos.py +++ b/qha/multi_configurations/different_phonon_dos.py @@ -127,6 +127,7 @@ def partition_functions_for_each_configuration(self): 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. @@ -164,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_all_configurations()], 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/tests/test_different_phonon_dos.py b/qha/tests/test_different_phonon_dos.py index 5230d77..b7c4c61 100644 --- a/qha/tests/test_different_phonon_dos.py +++ b/qha/tests/test_different_phonon_dos.py @@ -72,7 +72,7 @@ def test_partition_function(self): (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, + 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,