Coverage for mchammer/ensembles/thermodynamic_integration_ensemble.py: 93%
49 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-09-14 04:08 +0000
1from mchammer.free_energy_tools \
2 import (_lambda_function_forward, _lambda_function_backward)
4from ase import Atoms
5from ase.units import kB
7from .thermodynamic_base_ensemble import ThermodynamicBaseEnsemble
8from ..calculators.base_calculator import BaseCalculator
9from .. import DataContainer
10from icet.input_output.logging_tools import logger
12logger = logger.getChild('thermodynamic_integration_ensemble')
15class ThermodynamicIntegrationEnsemble(ThermodynamicBaseEnsemble):
16 r"""Instances of this class allow one to find the free energy of the
17 system. To this end, we use the :class:`canonncal ensemble
18 <mchammer.ensembles.CanonicalEnsemble>` with a modified
19 Hamiltonian,
21 .. math::
22 H(\lambda) = (1 - \lambda) H_{A} + \lambda H_{B}
24 The Hamiltonian is then sampled continuously from :math:`\lambda=0`
25 to :math:`\lambda=1`. :math:`H_{B}` is your cluster expansion
26 and :math:`H_{A}=0`, is a completely disordered system, with free
27 energy given by the ideal mixing entropy.
29 The free energy, A, of system B is then given by:
31 .. math::
32 A_{B} = A_{A} + \int_{0}^{1} \left\langle\frac{\mathrm{d}H(\lambda)}
33 {\mathrm{d}\lambda}\right\rangle_{H} \mathrm{d}\lambda
35 and since :math:`A_{A}` is known it is easy to compute :math:`A_{B}`
37 :math:`\lambda` is parametrized as,
39 .. math::
40 \lambda(x) = x^5(70x^4 - 315x^3 + 540x^2 - 420x + 126)
42 where :math:`x = \mathrm{step} / (\mathrm{n\_steps} - 1)`.
44 Parameters
45 ----------
46 structure
47 Atomic configuration to be used in the Monte Carlo simulation;
48 also defines the initial occupation vector.
49 calculator
50 Calculator to be used for calculating the potential changes
51 that enter the evaluation of the Metropolis criterion.
52 temperature
53 Temperature :math:`T` in appropriate units, commonly Kelvin.
54 n_lambdas
55 Number of :math:`\lambda` values to be sampled between 0 and 1.
56 forward
57 If this is set to ``True`` the simulation runs from :math:`H_A` to :math:`H_B`,
58 otherwise it runs from :math:`H_B` to :math:`H_A`.
59 :math:`H_B` is the cluster expansion and :math:`H_A = 0`, is the fully disordered system.
60 boltzmann_constant
61 Boltzmann constant :math:`k_B` in appropriate
62 units, i.e., units that are consistent
63 with the underlying cluster expansion
64 and the temperature units. Default: eV/K.
65 user_tag
66 Human-readable tag for ensemble. Default: ``None``.
67 random_seed
68 Seed for the random number generator used in the Monte Carlo simulation.
69 dc_filename
70 Name of file the data container associated with the ensemble
71 will be written to. If the file exists it will be read, the
72 data container will be appended, and the file will be
73 updated/overwritten.
74 data_container_write_period
75 Period in units of seconds at which the data container is
76 written to file. Writing periodically to file provides both
77 a way to examine the progress of the simulation and to back up
78 the data. Default: 600 s.
79 ensemble_data_write_interval
80 Interval at which data is written to the data container. This
81 includes for example the current value of the calculator
82 (i.e., usually the energy) as well as ensembles specific fields
83 such as temperature or the number of atoms of different species.
84 Default: Number of sites in the :attr:`structure`.
85 trajectory_write_interval
86 Interval at which the current occupation vector of the atomic
87 configuration is written to the data container.
88 Default: Number of sites in the :attr:`structure`.
89 sublattice_probabilities
90 Probability for picking a sublattice when doing a random swap.
91 This should be as long as the number of sublattices and should
92 sum up to 1.
95 Example
96 -------
97 The following snippet illustrate how to carry out a simple thermodynamic
98 integration. Here, the parameters of the cluster expansion are set to
99 emulate a simple Ising model in order to obtain an
100 example that can be run without modification. In practice, one should of
101 course use a proper cluster expansion::
103 >>> from ase.build import bulk
104 >>> from icet import ClusterExpansion, ClusterSpace
105 >>> from mchammer.calculators import ClusterExpansionCalculator
107 >>> # prepare cluster expansion
108 >>> # the setup emulates a second nearest-neighbor (NN) Ising model
109 >>> # (zerolet and singlet ECIs are zero; only first and second neighbor
110 >>> # pairs are included)
111 >>> prim = bulk('Au')
112 >>> cs = ClusterSpace(prim, cutoffs=[4.3], chemical_symbols=['Ag', 'Au'])
113 >>> ce = ClusterExpansion(cs, [0, 0, 0.1, -0.02])
115 >>> # prepare initial configuration
116 >>> structure = prim.repeat(3)
117 >>> for k in range(5):
118 >>> structure[k].symbol = 'Ag'
120 >>> # set up and run MC simulation
121 >>> calc = ClusterExpansionCalculator(structure, ce)
122 >>> mc = ThermodynamicIntegrationEnsemble(structure=structure, calculator=calc,
123 ... temperature=600,
124 ... n_steps=100000,
125 ... forward=True,
126 ... dc_filename='myrun_thermodynamic_integration.dc')
127 >>> mc.run()
129 """
131 def __init__(self,
132 structure: Atoms,
133 calculator: BaseCalculator,
134 temperature: float,
135 n_steps: int,
136 forward: bool,
137 user_tag: str = None,
138 boltzmann_constant: float = kB,
139 random_seed: int = None,
140 dc_filename: str = None,
141 data_container: str = None,
142 data_container_write_period: float = 600,
143 ensemble_data_write_interval: int = None,
144 trajectory_write_interval: int = None,
145 sublattice_probabilities: list[float] = None,
146 ) -> None:
148 self._ensemble_parameters = dict(temperature=temperature,
149 n_steps=n_steps)
150 self._last_state = dict()
152 super().__init__(
153 structure=structure,
154 calculator=calculator,
155 user_tag=user_tag,
156 random_seed=random_seed,
157 data_container=data_container,
158 dc_filename=dc_filename,
159 data_container_class=DataContainer,
160 data_container_write_period=data_container_write_period,
161 ensemble_data_write_interval=ensemble_data_write_interval,
162 trajectory_write_interval=trajectory_write_interval,
163 boltzmann_constant=boltzmann_constant)
165 if sublattice_probabilities is None: 165 ↛ 169line 165 didn't jump to line 169 because the condition on line 165 was always true
166 self._swap_sublattice_probabilities = \
167 self._get_swap_sublattice_probabilities()
168 else:
169 self._swap_sublattice_probabilities = sublattice_probabilities
171 sublattices = []
172 for sl in self.sublattices:
173 sublattices.append(sl.atomic_numbers)
175 # add species count to ensemble parameters
176 symbols = set([symbol for sub in calculator.sublattices
177 for symbol in sub.chemical_symbols])
178 for symbol in symbols:
179 key = 'n_atoms_{}'.format(symbol)
180 count = structure.get_chemical_symbols().count(symbol)
181 self._ensemble_parameters[key] = count
183 self._n_steps = n_steps
185 if forward:
186 self._lambda_function = _lambda_function_forward
187 self._lambda = 0
188 else:
189 self._lambda_function = _lambda_function_backward
190 self._lambda = 1
192 @property
193 def temperature(self) -> float:
194 """ Current temperature. """
195 return self._ensemble_parameters['temperature']
197 @property
198 def n_steps(self) -> int:
199 return self._n_steps
201 def _do_trial_step(self):
202 """ Carries out one Monte Carlo trial step. """
203 self._lambda = self._lambda_function(self.n_steps, self.step)
204 sublattice_index = self.get_random_sublattice_index(self._swap_sublattice_probabilities)
205 swap = self.do_thermodynamic_swap(sublattice_index=sublattice_index,
206 lambda_val=self._lambda)
207 return swap
209 def run(self):
210 """ Runs the thermodynamic integration. """
211 if self.step >= self.n_steps: 211 ↛ 212line 211 didn't jump to line 212 because the condition on line 211 was never true
212 logger.warning('The simulation is already done')
213 else:
214 super().run(self.n_steps - self.step)
216 def _get_ensemble_data(self):
217 data = super()._get_ensemble_data()
218 data['lambda'] = self._lambda
219 return data