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