wmaee.scopes.eos

  1"""
  2Convenient functions for fitting the energy-volume data.
  3"""
  4
  5import numpy as np
  6from scipy.optimize import curve_fit
  7from typing import Tuple, Optional, Iterable, Callable
  8from numpy.typing import ArrayLike, NDArray
  9
 10
 11def eos_birch_murnaghan(V: ArrayLike, E0: float, V0: float, B0: float, Bp: float) -> NDArray:
 12    """
 13    Birch-Murnaghan equation of state.
 14
 15    Parameters
 16    ----------
 17    V : ArrayLike
 18        Volume value(s).
 19    E0 : float
 20        Equilibrium energy.
 21    V0 : float
 22        Equilibrium volume.
 23    B0 : float
 24        Equilibrium bulk modulus.
 25    Bp : float
 26        Pressure derivative of the equilibrium bulk modulus.
 27
 28    Returns
 29    -------
 30    NDArray
 31        Energies corresponding to the given volumes.
 32
 33    Notes
 34    -----
 35    The Birch-Murnaghan equation of state is given by:
 36    E(V) = E0 + (9*B0*V0/16) * (((V0/V)**(2/3) - 1)**3 * Bp + ((V0/V)**(2/3) - 1)**2 * (6 - 4*(V0/V)**(2/3)))
 37    where E(V) is the energy as a function of volume (V).
 38    """
 39    return E0 + (9 * B0 * V0 / 16) * (((V0 / V)**(2/3) - 1)**3 * Bp + ((V0 / V)**(2/3) - 1)**2 * (6 - 4 * (V0 / V)**(2/3)))
 40
 41
 42def birch_murnaghan_fit(volumes: ArrayLike,
 43                        energies: ArrayLike,
 44                        p0: Optional[Iterable[float]] = None,
 45                        show_errors: bool = False,
 46                        eos_function: bool = False) -> Tuple[Iterable[float], Callable[[NDArray], NDArray]]:
 47    """
 48    Perform a Birch-Murnaghan fit on a given set of volumes and energies.
 49
 50    Parameters
 51    ----------
 52    volumes : ArrayLike
 53        Volumes.
 54    energies : ArrayLike
 55        Energies.
 56    p0 : Optional[Iterable[float]], optional
 57        Starting parameters in the order E0_guess, V0_guess, B0_guess, Bp_guess, by default None.
 58        If not provided, default values are calculated based on the input data.
 59    show_errors : bool, optional
 60        Whether to print the standard deviation of the fitted parameters, by default False.
 61    eos_function : bool, optional
 62        Whether to return the EOS function as a callable, by default False.
 63
 64    Returns
 65    -------
 66    Tuple[Iterable[float], Callable[[NDArray], NDArray]]
 67        Tuple containing the optimum parameters and, optionally, the EOS function.
 68
 69    optimum : Iterable[float]
 70        Optimal parameters E0, V0, B0, Bp.
 71    eos_function : Callable[[NDArray], NDArray], optional
 72        The EOS function as a callable.
 73
 74    Notes
 75    -----
 76    The Birch-Murnaghan equation of state is used for fitting.
 77
 78    If `eos_function` is True, the function returns the EOS function as a callable.
 79
 80    If `show_errors` is True, the function prints the standard deviation of the fitted parameters.
 81    
 82    If `p0` is not provided, default values are calculated as follows:
 83    - V0_guess: (np.amin(volumes) + np.amax(volumes))/2
 84    - E0_guess: np.amin(energies)
 85    - B0_guess: 1  # Presumably in eV/Ang^3; 160.2 in GPa
 86    - Bp_guess: 1
 87    """ 
 88    if p0 is None:
 89        V0_guess = (np.amin(volumes) + np.amax(volumes))/2
 90        E0_guess = np.amin(energies)
 91        B0_guess = 1  # Presumably in eV/Ang^3; 160.2 in GPa
 92        Bp_guess = 1
 93        p0 = (E0_guess, V0_guess, B0_guess, Bp_guess)
 94
 95    optimum, p_cov = curve_fit(eos_birch_murnaghan, np.array(volumes), np.array(energies), p0=p0)
 96    
 97    if show_errors:
 98        p_err = np.sqrt(np.diag(p_cov))
 99        print(f'Standard deviation of E0, V0, B0, Bp = {p_err}')
100
101    if eos_function:
102        return optimum, lambda V: eos_birch_murnaghan(V, *optimum)
103    
104    return optimum
105
106
107
108def eos_murnaghan(V: ArrayLike, E0: float, V0: float, B0: float, Bp: float) -> NDArray:
109    """
110    Murnaghan equation of state.
111
112    Parameters
113    ----------
114    V : ArrayLike
115        Volume value(s).
116    E0 : float
117        Equilibrium energy.
118    V0 : float
119        Equilibrium volume.
120    B0 : float
121        Equilibrium bulk modulus.
122    Bp : float
123        Pressure derivative of the equilibrium bulk modulus.
124
125    Returns
126    -------
127    NDArray
128        Energies corresponding to the given volumes.
129
130    Notes
131    -----
132    The Murnaghan equation of state is given by:
133    E(V) = E0 + B0 * V0 * (1/(Bp*(Bp-1)) * (V/V0)**(1-Bp) + V/(Bp*V0) - 1/(Bp-1))
134    where E(V) is the energy as a function of volume (V).
135    """
136    return E0 + B0 * V0 * (1/(Bp * (Bp - 1)) * (V / V0)**(1 - Bp) + V / (Bp * V0) - 1 / (Bp - 1))
137
138
139
140def murnaghan_fit(volumes: ArrayLike,
141                  energies: ArrayLike,
142                  p0: Optional[Iterable[float]] = None,
143                  show_errors: bool = False,
144                  eos_function: bool = False) -> Tuple[Iterable[float], Callable[[NDArray], NDArray]]:
145    """
146    Perform a Murnaghan fit on a given set of volumes and energies.
147
148    Parameters
149    ----------
150    volumes : ArrayLike
151        Volumes.
152    energies : ArrayLike
153        Energies.
154    p0 : Optional[Iterable[float]], optional
155        Starting parameters in the order E0_guess, V0_guess, B0_guess, Bp_guess, by default None.
156        If not provided, default values are calculated based on the input data.
157    show_errors : bool, optional
158        Whether to print the standard deviation of the fitted parameters, by default False.
159    eos_function : bool, optional
160        Whether to return the EOS function as a callable, by default False.
161
162    Returns
163    -------
164    Tuple[Iterable[float], Callable[[NDArray], NDArray]]
165        Tuple containing the optimum parameters and, optionally, the EOS function.
166
167    optimum : Iterable[float]
168        Optimal parameters E0, V0, B0, Bp.
169    eos_function : Callable[[NDArray], NDArray], optional
170        The EOS function as a callable.
171
172    Notes
173    -----
174    The Murnaghan equation of state is used for fitting.
175
176    If `eos_function` is True, the function returns the EOS function as a callable.
177
178    If `show_errors` is True, the function prints the standard deviation of the fitted parameters.
179
180    If `p0` is not provided, default values are calculated as follows:
181    - V0_guess: (np.amin(volumes) + np.amax(volumes))/2
182    - E0_guess: np.amin(energies)
183    - B0_guess: 1  # Presumably in eV/Ang^3; 160.2 in GPa
184    - Bp_guess: 1
185    """
186    if p0 is None:
187        V0_guess = (np.amin(volumes) + np.amax(volumes))/2
188        E0_guess = np.amin(energies)
189        B0_guess = 1  # Presumably in eV/Ang^3; 160.2 in GPa
190        Bp_guess = 1
191        p0 = (E0_guess, V0_guess, B0_guess, Bp_guess)
192
193    optimum, p_cov = curve_fit(eos_murnaghan, np.array(volumes), np.array(energies), p0=p0)
194    
195    if show_errors:
196        p_err = np.sqrt(np.diag(p_cov))
197        print(f'Standard deviation of E0, V0, B0, Bp = {p_err}')
198
199    if eos_function:
200        return optimum, lambda V: eos_murnaghan(V, *optimum)
201    
202    return optimum
203
204
205def polynomial_fit(volumes: ArrayLike,
206                   energies: ArrayLike,
207                   order: int = 3,
208                   eos_function: bool = False) -> Tuple[float, float, float, Optional[Callable[[NDArray], NDArray]]]:
209    """
210    Perform a polynomial fit on a given set of volumes and energies.
211
212    Parameters
213    ----------
214    volumes : ArrayLike
215        Volumes.
216    energies : ArrayLike
217        Energies.
218    order : int, optional
219        Order of the polynomial fit, by default 3.
220    eos_function : bool, optional
221        Whether to return the EOS function as a callable, by default False.
222
223    Returns
224    -------
225    Tuple[float, float, float, Optional[Callable[[NDArray], NDArray]]]
226        Tuple containing the equilibrium energy, equilibrium volume, bulk modulus,
227        and, optionally, the EOS function.
228
229    E_eq : float
230        Equilibrium energy.
231    V_eq : float
232        Equilibrium volume.
233    B_eq : float
234        Bulk modulus.
235    eos_function : Callable[[NDArray], NDArray], optional
236        The EOS function as a callable.
237
238    Notes
239    -----
240    The function performs a polynomial fit of the data and calculates equilibrium properties.
241
242    If `eos_function` is True, the function returns the EOS function as a callable.
243    """
244    # Extract the smallest and the biggest volume
245    V_min = np.amin(volumes)
246    V_max = np.amax(volumes)
247
248    # Polynomial fit of the data
249    coefficients = np.polyfit(volumes, energies, order)
250
251    # Calculate the coefficients of the first derivative
252    first_derivative = np.polyder(coefficients, 1)
253
254    # Calculate at which volume the first derivative is zero
255    roots = np.roots(first_derivative)
256
257    # Now we need to extract indices of the correct
258    # minimum volume from roots
259    v_min_indices = np.where((V_min <= roots) & (roots <= V_max))
260
261    # With v_min_indices we can extract the correct minimum
262    # volume from roots
263    V_eq = float(roots[v_min_indices])
264
265    # Further, we can now evaluate the polynomial at V_eq
266    # to get the equilibrium energy
267    E_eq = np.polyval(coefficients, V_eq)
268
269    # Now, we will calculate the bulk modulus, B_eq
270    # In general, B_eq = - V_eq * (dp/dV)
271    # p(V) = first derivative, hence dp/dV is the second
272    # derivative of the energy polynomial
273    second_derivative = np.polyder(coefficients, 2)
274
275    # Calculating the bulk modulus
276    B_eq = -V_eq * (-np.polyval(second_derivative, V_eq))
277
278    fitted_params = (E_eq, V_eq, B_eq)
279    
280    if eos_function:
281        return fitted_params, lambda V: np.polyval(coefficients, V)
282    
283    return fitted_params