Merge branch 'release/1.0.9'

Author: singularitti

qha/__init__.py

@@ -10,4 +10,4 @@
 __credits__ = {'Renata M. M. Wentzcovitch': '[email protected]'}
 __date__ = 'Feb 17, 2018'
 __maintainer__ = 'Tian Qin, Qi Zhang'
-__version__ = '1.0.8'
+__version__ = '1.0.9'

qha/bmf.py

@@ -29,6 +29,9 @@
 from qha.unit_conversion import gpa_to_ry_b3, ry_b3_to_gpa, b3_to_a3, ry_to_j_mol, ry_to_ev
 from qha.v2p import v2p
 
+# ===================== What can be exported? =====================
+__all__ = ['Calculator', 'SamePhDOSCalculator', 'DifferentPhDOSCalculator']
+
 
 class Calculator:
     def __init__(self, user_settings: Dict[str, Any]):
@@ -350,8 +353,8 @@ def vib_ry(self):
         v = np.empty(self.temperature_array.shape)
 
         for i, t in enumerate(self.temperature_array):
-            v[i] = same_phonon_dos.FreeEnergy(t, self.volume_energy.as_matrix(), self.degeneracies,
-                                              self.q_weights, self.frequencies).total
+            v[i] = same_phonon_dos.FreeEnergy(t, self.degeneracies, self.q_weights, self.volume_energy.as_matrix(),
+                                              self.frequencies).total
         return v
 
 
@@ -415,7 +418,7 @@ 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.static_energies, self.degeneracies, self.q_weights, self.frequencies, self._volumes, \
+        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]))

qha/calculator.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+"""
+:mod: bmf_all_t
+================================
+
+.. module bmf_all_t
+   :platform: Unix, Windows, Mac, Linux
+   :synopsis:
+.. moduleauthor:: Tian Qin <[email protected]>
+"""
+
+from typing import Optional
+
+import numpy as np
+from numpy.linalg import inv
+
+# ===================== What can be exported? =====================
+__all__ = ['polynomial_least_square_fitting', 'birch_murnaghan_finite_strain_fitting']
+
+
+def polynomial_least_square_fitting(x, y, new_x, order: Optional[int] = 3):
+    """
+    Equation of calculate the coefficients are from
+    `Wolfram Mathworld <http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html>`_.
+
+    :param x: Eulerian strain of calculated volumes (sparse)
+    :param y: Free energy of these calculated volumes (sparse)
+    :param new_x: Eulerian strain at a greater dense vector
+    :param order: orders to fit Birch Murnaghan EOS
+    :return: Free energy at a denser strain vector (denser volumes vector)
+    """
+    order += 1  # The order needed is 1 more than ``numpy.vander`` default value.
+    X = np.vander(x, order, increasing=True)  # This will make a Vandermonde matrix that will be used in BM fitting.
+    X_T = X.T
+    a = inv(X_T @ X) @ X_T @ y
+    new_y = np.vander(new_x, order, increasing=True) @ a
+    return a, new_y
+
+
+def birch_murnaghan_finite_strain_fitting(eulerian_strain, free_energy, strain, order: Optional[int] = 3):
+    """
+    Calculate the F(T,V) for given strain.
+
+    :param eulerian_strain: Eulerian strain of calculated volumes (sparse)
+    :param free_energy: Free energy of these calculated volumes (sparse)
+    :param strain: Eulerian strain at a greater dense vector
+    :param order: orders to fit Birch Murnaghan EOS
+    :return: Free energy at a dense (T, V) grid.
+    """
+    temperature_amount, _ = free_energy.shape
+    dense_volume_amount = len(strain)
+    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(eulerian_strain, free_energy[i], strain, order)
+        f_v_t[i] = f_i
+    return f_v_t

qha/fitting.py

@@ -16,7 +16,7 @@
 from numba import vectorize, float64, jit, int64
 from numba.types import UniTuple
 
-from qha.bmf import bmf_all_t, bmf
+from qha.fitting import polynomial_least_square_fitting, birch_murnaghan_finite_strain_fitting
 from qha.type_aliases import Vector, Matrix
 from qha.unit_conversion import gpa_to_ry_b3
 
@@ -83,7 +83,7 @@ def approaching_to_best_ratio(self, volumes: Vector, free_energies: Vector, init
         """
         strains, finer_volumes = interpolate_volumes(volumes, self.ntv, initial_ratio)
         eulerian_strain = calc_eulerian_strain(volumes[0], volumes)
-        f_v_tmax = bmf(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.p_desire)
         # find the index of the first pressure value that slightly smaller than p_desire
@@ -97,8 +97,8 @@ def refine_grids(self, volumes: Vector, free_energies: Matrix,
         Get the appropriate volume grid for interpolation.
         Avoid to use a too large volume grid to obtain data, which might lose accuracy.
         :param free_energies: Calculated Helmholtz Free Energies for input volumes (sparse).
-        :param volumes: olumes of these calculations were perform (sparse).
-        :param ratio:  This ratio is used to get a larger volume grid
+        :param volumes: Volumes of these calculations were perform (sparse).
+        :param ratio: This ratio is used to get a larger volume grid
         :return: volume, Helmholtz free energy at a denser vector, and the `ratio` used in this calculation
         """
         if ratio is not None:
@@ -112,6 +112,5 @@ def refine_grids(self, volumes: Vector, free_energies: Matrix,
 
         eulerian_strain = calc_eulerian_strain(volumes[0], volumes)
         strains, finer_volumes = interpolate_volumes(volumes, self.ntv, new_ratio)
-        f_tv = bmf_all_t(eulerian_strain, free_energies, strains, self.option)
-
-        return finer_volumes, f_tv, new_ratio
+        f_tv_bfm = birch_murnaghan_finite_strain_fitting(eulerian_strain, free_energies, strains, self.option)
+        return finer_volumes, f_tv_bfm, new_ratio

qha/grid_interpolation.py

@@ -10,13 +10,15 @@
 .. moduleauthor:: Qi Zhang <[email protected]>
 """
 
+from typing import Optional
+
 import numpy as np
 from lazy_property import LazyProperty
 from scipy.constants import Boltzmann
 from scipy.special import logsumexp
 
 import qha.settings
-from qha.bmf import bmf
+from qha.fitting import polynomial_least_square_fitting
 from qha.grid_interpolation import calc_eulerian_strain
 from qha.single_configuration import free_energy
 from qha.type_aliases import Array4D, Scalar, Vector, Matrix
@@ -33,8 +35,9 @@
 
 
 class PartitionFunction:
-    def __init__(self, temperature: Scalar, static_energies: Matrix, degeneracies: Vector, q_weights: Matrix,
-                 frequencies: Array4D, volumes: Matrix, static_only: bool, precision: int = 500):
+    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(
@@ -59,7 +62,7 @@ def __init__(self, temperature: Scalar, static_energies: Matrix, degeneracies: V
         self.volumes = volumes
         self.static_only = static_only
         self.precision = int(precision)
-        self.__ntv = 400
+        self.order = int(order)
 
     @LazyProperty
     def helmoholtz_configs(self):
@@ -79,7 +82,8 @@ def helmholtz_at_ref_v(self):
             # strains, finer_volumes[i, :] = interpolate_volumes(self.volumes[i], self.__ntv, 1.05)
             eulerian_strain = calc_eulerian_strain(self.volumes[i][0], self.volumes[i])
             strains = calc_eulerian_strain(self.volumes[i][0], self.volumes[0])
-            helmholtz_fitted[i, :] = bmf(eulerian_strain, self.helmoholtz_configs[i], strains)
+            _, helmholtz_fitted[i, :] = polynomial_least_square_fitting(eulerian_strain, self.helmoholtz_configs[i],
+                                                                        strains, order=self.order)
         return helmholtz_fitted
 
     @LazyProperty

qha/multi_configurations/different_phonon_dos.py

@@ -9,6 +9,8 @@
 .. moduleauthor:: Qi Zhang <[email protected]>
 """
 
+from typing import Optional
+
 import numpy as np
 from lazy_property import LazyProperty
 from scipy.constants import Boltzmann
@@ -30,8 +32,8 @@
 
 
 class PartitionFunction:
-    def __init__(self, temperature: Scalar, static_energies: Matrix, degeneracies: Vector, q_weights: Vector,
-                 frequencies: Array3D, precision: 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!')
@@ -88,7 +90,7 @@ def derive_free_energy(self):
 
 
 class FreeEnergy:
-    def __init__(self, temperature: Scalar, static_energies: Matrix, degeneracies: Vector, q_weights: Vector,
+    def __init__(self, temperature: Scalar, degeneracies: Vector, q_weights: Vector, static_energies: Matrix,
                  frequencies: Array3D):
         if not np.all(np.greater_equal(degeneracies, 0)):
             raise ValueError('Degeneracies should all be integers greater equal than 0!')