Coverage for mchammer/free_energy_tools.py: 88%
88 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-05-06 04:14 +0000
1from typing import List
2from icet import ClusterSpace
3from icet.core.sublattices import Sublattices
5import numpy as np
6from scipy.integrate import cumulative_trapezoid
8from ase import Atoms
9from ase.units import kB
11from mchammer import DataContainer
13import operator
14from collections import Counter
15from math import factorial
16from functools import reduce
19def _lambda_function_forward(n_steps: int, step: int) -> float:
20 """
21 Returns the current lambda value for a backward calculation.
23 Parameters
24 ----------
25 n_steps
26 Total number of steps..
27 step
28 Current step.
29 """
30 x = (step + 1) / (n_steps)
31 lam = x**5*(70*x**4 - 315*x**3 + 540*x**2 - 420*x + 126)
32 # Due to numerical precision lambda may be slightly outside the interval [0.0, 1.0]
33 lam = np.clip(lam, 0.0, 1.0)
34 return lam
37def _lambda_function_backward(n_steps: int, step: int) -> float:
38 """
39 Returns the current lambda value for a backward calculation.
41 Parameters
42 ----------
43 n_steps
44 Total number of steps.
45 step
46 Current step.
47 """
48 return 1 - _lambda_function_forward(n_steps, step)
51def _stirling(n: int) -> float:
52 """ Stirling approximation of log(n!). """
53 return n * np.log(n) - n + 0.5 * np.log(2*np.pi*n)
56def _lognpermutations(array: List[int]) -> float:
57 """If _npermutations becomes too large we resort to calculating the
58 logarithm directly using Stirling's approximation, i.e., we
59 calculate
61 .. math::
63 log(n!) - (log(a1!) + log(a2!) + log(a3!) + ... + log(ak!)
65 where we denote in the code
67 .. math::
69 v n = log(n!)
70 den = (log(a1!) + log(a2!) + log(a3!) + ... + log(ak!)
72 """
74 n = _stirling(len(array))
75 # Gets the number of each unique atom type
76 mults = Counter(array).values()
77 den = reduce(operator.add,
78 (_stirling(v) for v in mults),
79 1)
80 return n - den
83def _npermutations(array: List[int]) -> float:
84 """
85 Returns the total number of ways we can permutate the array which is given by
87 .. math::
89 n! / (a1! * a2! * a3! * ... * ak!)
91 where we denote in the code
93 .. math::
95 den = (a1! * a2! * a3! * ... * ak!)
96 """
97 n = factorial(len(array))
98 # Gets the number of each unique atom type
99 a = Counter(array).values()
100 den = reduce(operator.mul,
101 (factorial(v) for v in a),
102 1)
103 return n / den
106def _ideal_mixing_entropy(swap_sublattice_probabilities: List[int],
107 atoms_on_sublattices: np.ndarray,
108 boltzmann_constant: float) -> float:
109 """
110 Calculates the ideal mixing entropy for the supercell.
112 Parameters
113 ----------
114 swap_sublattice_probabilites
115 Sublattices that are active during the MC simulation.
116 atoms_on_sublattices
117 Sorted atomic numbers in sublattices lists.
118 boltzmann_constant
119 Value of Boltzmann constant in the units used in the MC simulation.
120 """
121 log_multiplicity = []
122 sublattices = np.array(swap_sublattice_probabilities) > 0
123 for aos in atoms_on_sublattices[sublattices]:
124 try:
125 log_multiplicity.append(np.log(_npermutations(aos)))
126 except Exception:
127 log_multiplicity.append(_lognpermutations(aos))
128 return boltzmann_constant * np.sum(log_multiplicity)
131def _get_atoms_on_sublattice(structure: Atoms,
132 sublattices: Sublattices) -> np.ndarray:
133 """Sorts the atom symbols in the structure to lists involving
134 specific sublattices, for example,
136 [[5, 5, 3, 3, 3], [11, 11, 10, 10, 10 10, 10]]
138 Parameters
139 ----------
140 structure
141 Atomic structure.
142 sublattices
143 Sublattices for the supercell obtained from the cluster space.
144 """
145 _atoms_on_sublattices = []
146 for sl in sublattices:
147 atoms_on_sublattice = []
148 for an in sl.atomic_numbers:
149 natoms = np.where(structure.numbers == an)[0].size
150 atoms_on_sublattice.extend([an] * natoms)
151 _atoms_on_sublattices.append(atoms_on_sublattice)
152 _atoms_on_sublattices = np.array(_atoms_on_sublattices,
153 dtype=object)
154 return _atoms_on_sublattices
157def get_free_energy_thermodynamic_integration(
158 dc: DataContainer,
159 cluster_space: ClusterSpace,
160 forward: bool,
161 max_temperature: float = np.inf,
162 sublattice_probabilities: List[float] = None,
163 boltzmann_constant: float = kB,
164) -> (np.ndarray, np.ndarray):
165 r"""Returns the free energy calculated via thermodynamic integration
166 using the :class:`ThermodynamicIntegrationEnsemble
167 <mchammer.ensembles.ThermodynamicIntegrationEnsemble>`.
169 The temperature dependence of the free energy can be extracted
170 from the thermodynamic integration as
172 .. math::
174 F(T) = \frac{F_0(\lambda)}{\lambda} + \frac{T_0}{\lambda} S_\text{B}
176 where :math:`S_\text{B}` is the Boltzmann entropy,
178 .. math::
180 T = \frac{T_0}{\lambda}
182 and
184 .. math::
186 F_0(\lambda) = \int_0^\lambda \left\langle\frac{\mathrm{d}H(\lambda)}
187 {\mathrm{d}\lambda}\right\rangle_{H} \mathrm{d}\lambda
189 Parameters
190 ----------
191 dc
192 Data container from the thermodynamic integration simulation.
193 cluster_space
194 Cluster space used to construct the cluster expansion used for the simulation.
195 forward
196 Whether or not the thermodynamic integration was carried out forward or backward.
197 max_temperature
198 Largest temperature to extract from the thermodynamic integration.
199 sublattice_probababilites
200 Sublattice probabilties that were provided to the thermodynamic integration
201 simulation.
202 boltzmann_constant
203 Boltzmann constant in the units used for the thermodynamic integration
204 simulation.
206 """
208 lambdas, potentials = dc.get('lambda', 'potential')
209 if not forward:
210 # We want to integrate from high to low temperature even though
211 # the simulation was from low to high.
212 potentials = potentials[::-1]
213 lambdas = lambdas[::-1]
215 sublattices = cluster_space.get_sublattices(dc.structure)
216 if sublattice_probabilities is None:
217 sublattice_probabilities = [True] * len(sublattices)
219 atoms_on_sublattices = _get_atoms_on_sublattice(dc.structure,
220 sublattices)
221 ideal_mixing_entropy = _ideal_mixing_entropy(sublattice_probabilities,
222 atoms_on_sublattices,
223 boltzmann_constant)
225 temperature = dc._ensemble_parameters['temperature']
226 with np.errstate(divide='ignore'):
227 # First value of lambdas will be zero (we ignore this warning)
228 temperatures = (1 / lambdas) * temperature
229 temp_inds = np.where(temperatures <= max_temperature)[0]
231 free_energy_change = cumulative_trapezoid(potentials, x=lambdas, initial=0)
232 with np.errstate(invalid='ignore', divide='ignore'):
233 # First value of lambdas will be zero (we ignore this warning)
234 free_energy = (free_energy_change / lambdas -
235 temperatures * ideal_mixing_entropy)
236 free_energy[np.isnan(free_energy)] = 0
237 return (temperatures[temp_inds], free_energy[temp_inds])
240def get_free_energy_temperature_integration(
241 dc: DataContainer,
242 cluster_space: ClusterSpace,
243 forward: bool,
244 temperature_reference: float,
245 free_energy_reference: float = None,
246 sublattice_probabilities: List[float] = None,
247 max_temperature: float = np.inf,
248 boltzmann_constant: float = kB,
249) -> (np.ndarray, np.ndarray):
250 r""" Returns the free energy calculated using temperature integration and
251 the corresponding temperature.
253 .. math::
255 \frac{A(T_{2})}{T_{2}} = \frac{A(T_{1})}{T_{1}}
256 - \int_{T_{1}}^{T_{2}}\frac{U(T)}{T^2}\mathrm{d}T
258 Parameters
259 ----------
260 dc
261 Data container from a canonical annealing simulation.
262 The first (last for :attr:`forward=False`) temperature in the
263 data container has to be at the same temperature as
264 :attr:`temperature_reference`.
265 cluster_space
266 Cluster space used to construct the cluster expansion
267 that was used in the simulations.
268 forward
269 If ``True`` the canonical annealing simulation was carried out
270 from high to low temperature, otherwise the opposite is assumed.
271 temperature_reference
272 Temperature at which :attr:`free_energy_reference` was calculated
273 free_energy_reference
274 Reference free energy. If set to ``None`` it will be assumeed that the
275 free energy at :attr:`temperature_reference` can be approximated by
276 :math:`T S_B` where :math:`S_B` is the ideal mixing entropy.
277 sublattice_probababilites
278 Sublattice probabilties that were provided to the canonical annealing
279 simulation.
280 max_temperature
281 Largest temperature to extract from the temperature integration.
282 boltzmann_constant
283 Boltzmann constant in the units used for the thermodynamic integration
284 simulation.
285 """
286 temperatures, potentials = dc.get('temperature', 'potential')
287 if not forward:
288 # We want to integrate from high to low temperature even though
289 # the simulation was from low to high.
290 potentials = potentials[::-1]
291 temperatures = temperatures[::-1]
293 if not np.isclose(temperatures[0], temperature_reference): 293 ↛ 294line 293 didn't jump to line 294, because the condition on line 293 was never true
294 raise ValueError(
295 'The first or the last temperature in the temperature list'
296 ' has to equal the reference temperature.'
297 f' reference_temperature: {temperature_reference}'
298 f' first canonical temperature: {temperatures[0]}')
300 temp_inds = np.where(temperatures <= max_temperature)[0]
302 if free_energy_reference is None: 302 ↛ 303line 302 didn't jump to line 303, because the condition on line 302 was never true
303 sublattices = cluster_space.get_sublattices(dc.structure)
304 if sublattice_probabilities is None:
305 sublattice_probabilities = [True] * len(sublattices)
307 atoms_on_sublattices = _get_atoms_on_sublattice(dc.structure,
308 sublattices)
309 ideal_mixing_entropy = _ideal_mixing_entropy(sublattice_probabilities,
310 atoms_on_sublattices,
311 boltzmann_constant)
312 free_energy_reference = -ideal_mixing_entropy * temperature_reference
314 reference = free_energy_reference / temperature_reference
315 integral = cumulative_trapezoid(potentials / temperatures**2,
316 x=temperatures, initial=0)
317 free_energy = (reference - integral) * temperatures
318 return (temperatures[temp_inds], free_energy[temp_inds])