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