From a372bf890d0bee3cb1bb029dc3dcb6c477990488 Mon Sep 17 00:00:00 2001 From: Riccardo Farris Date: Mon, 27 Jul 2026 15:33:19 +0200 Subject: [PATCH 1/5] Fix the batched replica-exchange mu-ladder swap rule, and six review findings BatchedReplicaExchange._accept_swap implemented only the temperature-ladder criterion (beta_j - beta_i)(Phi_j - Phi_i), which is identically zero when the replicas share a temperature. Every mu-ladder swap was therefore accepted with p = 1, so configurations random-walked freely across the ladder and no replica sampled its own mu. Replace it with the general form beta_i Phi_X^(i) + beta_j Phi_Y^(j) - beta_i Phi_Y^(i) - beta_j Phi_X^(j) where Phi_Z^(k) is a configuration scored with slot k's chemical potentials. This reduces to the old expression for a shared mu, to beta (mu_i - mu_j) (N_j - N_i) for a shared temperature, and also covers a joint (T, mu) ladder. An EMT Ag(111)/Au ladder gives N_Au = [0, 0, 0, 14] with the fix and [0, 0, 0, 0] without it: unconditional swapping erased the mu dependence rather than merely degrading the statistics. Existing mu-ladder results are invalid and need re-running; temperature ladders were unaffected, as was the MPI ReplicaExchange, which selected _exchange_prob_mu all along. A degenerate ladder now warns at the end of a run. A 100% swap acceptance sat unremarked in the log of six consecutive production runs, so the all-accept and never-accept tallies say so explicitly. Checked on the whole-run tally rather than live, because a ladder whose replicas have not differentiated yet legitimately accepts every early swap. Also fixed: - PermutationMove with n_swaps > 1 could apply k swaps and then return the "could not propose" sentinel when a later iteration drew a species absent from the system. The ensembles read that sentinel as "the atoms were not touched" and skipped the rollback, leaving the configuration out of step with E_old. Species counts are swap-invariant, so the usable species are now resolved once, up front. Both ensembles also restore their snapshot on the sentinel path rather than trusting the contract. - Molecules whose centre of mass drifted above a CustomCell's top stopped being deletion candidates and dropped out of the per-species de Broglie count, inflating V/((N+1)Lambda^3) into the runaway insertion mode that docs/gcmc_acceptance_convention.rst describes. Cells grow a second predicate, is_point_exchangeable, which CustomCell overrides with the same dropped z upper bound that get_atoms_specie_inside_cell already applied to single atoms. find_molecules and the displacement one-way-door guard both use it, so candidacy and displacement regions always agree. - AlchemiCalculator(energy_only=True) discarded 'forces' from the model config, which a pre-loaded MACEWrapper shares with every other calculator built from it, silently disabling FIRE relaxation in an AlchemiFCalculator depending on construction order. The combination now raises. - GrandCanonicalEnsemble/CanonicalEnsemble.set_state restored the step count and exchange statistics from the incoming state. ReplicaExchange passes a full get_state() dict on every accepted swap, so the two ranks traded their swap tallies. Only the configuration travels now. - GrandCanonicalEnsemble.write_outfile returned before reset_counters() when the outfile was disabled, silently turning interval_ratios() into total_ratios() for the rest of the run. - MoveSelector abbreviated move names with a three-character slice, collapsing MoleculeInsertionMove, MoleculeDeletionMove and MoleculeDisplacementMove to one indistinguishable 'Mol'. Single-word moves keep their existing labels so old outfiles stay comparable. - A missing entry in species_radii raised a bare KeyError from inside the free-volume sampler; it now names the species. CI gains the flake8 job CLAUDE.md already documented. --- .github/workflows/tests.yml | 20 +++ docs/gcmc_acceptance_convention.rst | 15 ++ docs/reference/calculators.rst | 6 +- docs/reference/cells.rst | 11 +- examples/re_gcmc_co_cupd_batched.py | 9 +- mcpy/calculators/alchemi_calculator.py | 14 ++ mcpy/cell/cell.py | 36 +++++ mcpy/cell/custom_cell.py | 22 ++- mcpy/cell/dome_cell.py | 6 +- mcpy/cell/spherical_cell.py | 6 +- mcpy/ensembles/batched_replica_exchange.py | 80 ++++++++-- mcpy/ensembles/canonical_ensemble.py | 8 +- mcpy/ensembles/grand_canonical_ensemble.py | 38 +++-- mcpy/moves/molecule_displacement_move.py | 7 +- mcpy/moves/molecule_utils.py | 9 +- mcpy/moves/move_selector.py | 19 ++- mcpy/moves/permutation_move.py | 12 +- tests/test_audit_regressions.py | 134 ++++++++++++++++ tests/test_batched_re_acceptance.py | 175 +++++++++++++++++++++ tests/test_compound_moves.py | 43 +++++ tests/test_custom_cell_region.py | 81 ++++++++++ tests/test_molecule_moves.py | 6 +- 22 files changed, 692 insertions(+), 65 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 21dcb90..fddfb1e 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -34,3 +34,23 @@ jobs: - name: Run tests run: pytest tests/ -v + + lint: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + + - uses: actions/setup-python@v5 + with: + python-version: '3.12' + cache: 'pip' + + - name: Install flake8 + run: | + python -m pip install --upgrade pip + pip install flake8 + + # Same command CLAUDE.md documents; config lives in .flake8. + - name: Lint + run: flake8 mcpy/ diff --git a/docs/gcmc_acceptance_convention.rst b/docs/gcmc_acceptance_convention.rst index e9d8516..0798eb9 100644 --- a/docs/gcmc_acceptance_convention.rst +++ b/docs/gcmc_acceptance_convention.rst @@ -152,6 +152,21 @@ Atomic moves are unaffected: for them ``get_exchange_count()`` returns ``None`` and the ensembles fall back to the total-atom count documented above. +Because molecular moves *do* use the per-species count, they are exposed to +the runaway mode described in the previous section: any escape route out of +the counted region undercounts ``N`` and inflates +:math:`V/((N+1)\Lambda^3)`. Molecule candidacy therefore goes through +``cell.is_point_exchangeable(com)``, not ``cell.is_point_inside(com)``. The +two differ for :class:`mcpy.cell.CustomCell`, whose exchangeable region drops +the z upper bound exactly as ``get_atoms_specie_inside_cell`` already does for +single atoms, so a molecule that desorbs above the cell top stays both +countable and deletable. The floor exclusion is kept on both paths: molecules +that sink below ``bottom_z`` remain excluded, and the buried-species +irreversibility noted above applies to them unchanged. +:class:`mcpy.moves.MoleculeDisplacementMove` tests the same predicate for its +one-way-door guard, so the candidacy region and the region a molecule may be +displaced within always agree. + The textbook form above holds for ``min_insert=None``. When ``min_insert`` is set, ``MoleculeInsertionMove`` retries the random position/orientation draw (up to 1000 times) against the cell's ``species_radii`` atoms until it diff --git a/docs/reference/calculators.rst b/docs/reference/calculators.rst index b1e2855..03c3740 100644 --- a/docs/reference/calculators.rst +++ b/docs/reference/calculators.rst @@ -89,7 +89,11 @@ chunk_size=None)``, and ``run_md(...)``. ``get_potential_energies``; caps peak GPU memory at one chunk. ``None`` evaluates the whole batch in one pass. See :doc:`../calculators`. - ``energy_only`` (bool): drop force computation (no autograd graph) for a - ~12% memory saving. Energy is unchanged; forces are unavailable. + ~12% memory saving. Energy is unchanged; forces are unavailable. The switch + lives on the model's own config, so it cannot be scoped to one calculator: + combining it with a pre-loaded ``MACEWrapper`` (which other calculators may + share, including a relaxing ``AlchemiFCalculator``) raises ``ValueError``. + Pass a checkpoint path instead so this calculator loads its own model. AlchemiFCalculator diff --git a/docs/reference/cells.rst b/docs/reference/cells.rst index 6534f83..4e966e1 100644 --- a/docs/reference/cells.rst +++ b/docs/reference/cells.rst @@ -32,9 +32,18 @@ The full ASE simulation box as the active region. Random points are sampled uniformly in fractional coordinates. Its volume is the fixed box volume, so ``calculate_volume`` does no free-volume sampling. -- ``species_radii`` (dict, optional): per-element exclusion radii. +- ``species_radii`` (dict, optional): per-element exclusion radii. Every + species the system can hold needs one, including species inserted during the + run; a missing entry raises ``ValueError`` from ``calculate_volume``. - ``seed`` (int, optional): cell-local RNG seed. +Every cell exposes two point-membership predicates. ``is_point_inside(point)`` +bounds the *proposal* region -- the one ``get_random_point`` samples. +``is_point_exchangeable(point)`` bounds the region the reservoir may take +molecules back from, and is what the molecule moves use to pick candidates. +They coincide for every cell except ``CustomCell``; see +:doc:`../gcmc_acceptance_convention`. + CustomCell ---------- diff --git a/examples/re_gcmc_co_cupd_batched.py b/examples/re_gcmc_co_cupd_batched.py index d604b42..d531105 100644 --- a/examples/re_gcmc_co_cupd_batched.py +++ b/examples/re_gcmc_co_cupd_batched.py @@ -46,9 +46,12 @@ def parse_args(): p.add_argument('--gcmc-steps', type=int, default=80, help='GCMC steps per replica') p.add_argument('--exchange-interval', type=int, default=8, - help='Steps between replica-exchange attempts; use >= 50 ' - 'for per-mu coverage curves (small values over-mix ' - 'the ladder)') + help='Steps between replica-exchange attempts. Frequent ' + 'exchange is fine: a correct swap rule leaves each ' + "rung's distribution intact. (The old advice to use " + '>= 50 here was working around the mu-ladder swap ' + 'bug, which accepted every swap and so did flatten ' + 'the coverage curve the more often it ran.)') p.add_argument('--checkpoint', default='medium-mpa-0') p.add_argument('--rel-steps', type=int, default=30) p.add_argument('--rel-fmax', type=float, default=0.1) diff --git a/mcpy/calculators/alchemi_calculator.py b/mcpy/calculators/alchemi_calculator.py index 2b6c4ee..e9cebd2 100644 --- a/mcpy/calculators/alchemi_calculator.py +++ b/mcpy/calculators/alchemi_calculator.py @@ -54,6 +54,20 @@ def __init__( energy_only: bool = False, head: Union[str, int, None] = None, ) -> None: + if energy_only and isinstance(checkpoint, MACEWrapper): + # energy_only mutates the wrapper's own model_config (below), and a + # pre-loaded wrapper is shared by every calculator built from it — + # including an AlchemiFCalculator whose FIRE relaxation needs the + # forces this would switch off. There is no way to scope the change + # to one calculator, so refuse instead of breaking the other one + # silently and order-dependently. + raise ValueError( + "energy_only=True cannot be combined with a pre-loaded " + "MACEWrapper: dropping 'forces' from its active outputs would " + "also disable forces for every other calculator sharing that " + "wrapper. Pass the checkpoint path so this calculator loads " + "its own model, or set energy_only on all of the sharers." + ) self.device = device self.dtype = dtype self.max_neighbors = max_neighbors diff --git a/mcpy/cell/cell.py b/mcpy/cell/cell.py index 7d172d6..581a785 100644 --- a/mcpy/cell/cell.py +++ b/mcpy/cell/cell.py @@ -66,6 +66,28 @@ def get_species(self): """ return list(self.species_radii.keys()) + def _radii_for(self, atoms): + """Per-atom exclusion radii for the free-volume samplers. + + Raises a message that names the missing symbols instead of the bare + ``KeyError`` a dict lookup would throw from inside the sampler. The + check lives here rather than in ``__init__`` because GCMC inserts + species that are absent when the cell is built. + + :return: ndarray of radii, one per atom, in ``atoms`` order. + """ + symbols = atoms.get_chemical_symbols() + missing = sorted(set(symbols) - set(self.species_radii)) + if missing: + raise ValueError( + f'{type(self).__name__}.species_radii has no radius for ' + f'{missing}. Every species the system can contain needs one, ' + f'including species inserted during the run; got ' + f'{sorted(self.species_radii)}.' + ) + return np.fromiter((self.species_radii[s] for s in symbols), + dtype=float, count=len(symbols)) + def is_point_inside(self, point): """The box cell spans the whole periodic cell: every point is inside. @@ -73,3 +95,17 @@ def is_point_inside(self, point): any cell type (the region cells implement a real test). """ return True + + def is_point_exchangeable(self, point): + """Whether a molecule whose center of mass sits at ``point`` may be + exchanged with the reservoir (counted, deleted, displaced). + + The point counterpart of :meth:`get_atoms_specie_inside_cell`, and the + predicate the molecule moves use. It is separate from + :meth:`is_point_inside` because a region may deliberately accept + molecules it would never *propose* — :class:`CustomCell` drops the z + upper bound so a desorbed molecule stays deletable instead of + accumulating forever. Cells whose two regions coincide (the box, the + sphere, the dome) inherit this delegation unchanged. + """ + return self.is_point_inside(point) diff --git a/mcpy/cell/custom_cell.py b/mcpy/cell/custom_cell.py index a17374a..d796e88 100644 --- a/mcpy/cell/custom_cell.py +++ b/mcpy/cell/custom_cell.py @@ -75,12 +75,7 @@ def calculate_volume(self, atoms) -> float: cart_coords = frac_coords @ self.dimensions # in cell frame positions = self._periodic_images(atoms) - symbols = atoms.get_chemical_symbols() - radii = np.fromiter( - (self.species_radii[s] for s in symbols), - dtype=float, count=n_atoms, - ) - radii = np.tile(radii, len(positions) // n_atoms) + radii = np.tile(self._radii_for(atoms), len(positions) // n_atoms) covered = np.zeros(self.mc_sample_points, dtype=bool) for r in np.unique(radii): @@ -144,6 +139,21 @@ def is_point_inside(self, point): frac = (point - self.offset) @ self._dim_inv return bool(np.all((frac >= 0.0) & (frac < 1.0))) + def is_point_exchangeable(self, point): + """Same asymmetry as :meth:`get_atoms_specie_inside_cell`: inside the + xy footprint and at or above the cell floor, with the z upper bound + dropped. + + Without this a molecule that desorbs above the cell top would stop + being a deletion candidate *and* would drop out of the per-species + count in the acceptance factor ``V/((N+1)Λ³)`` -- the runaway + insertion mode described in docs/gcmc_acceptance_convention.rst. + Molecules below the floor stay excluded, like atoms. + """ + frac = (point - self.offset) @ self._dim_inv + return bool(np.all(frac[:2] >= 0.0) and np.all(frac[:2] < 1.0) + and frac[2] >= 0.0) + def get_species(self): """ Get the species present in the custom cell. diff --git a/mcpy/cell/dome_cell.py b/mcpy/cell/dome_cell.py index 1a5b71c..2a3990a 100644 --- a/mcpy/cell/dome_cell.py +++ b/mcpy/cell/dome_cell.py @@ -121,11 +121,7 @@ def calculate_volume(self, atoms): self.volume = dome_fraction * self.sphere_volume return - symbols = atoms.get_chemical_symbols() - radii = np.fromiter( - (self.species_radii[s] for s in symbols), - dtype=float, count=n_atoms, - ) + radii = self._radii_for(atoms) positions = atoms.positions covered = np.zeros(len(pts), dtype=bool) diff --git a/mcpy/cell/spherical_cell.py b/mcpy/cell/spherical_cell.py index 75c919f..e1bbb03 100644 --- a/mcpy/cell/spherical_cell.py +++ b/mcpy/cell/spherical_cell.py @@ -121,11 +121,7 @@ def calculate_volume(self, atoms): pts = self._sample_sphere_points(self.mc_sample_points) - symbols = atoms.get_chemical_symbols() - radii = np.fromiter( - (self.species_radii[s] for s in symbols), - dtype=float, count=n_atoms, - ) + radii = self._radii_for(atoms) positions = atoms.positions covered = np.zeros(self.mc_sample_points, dtype=bool) diff --git a/mcpy/ensembles/batched_replica_exchange.py b/mcpy/ensembles/batched_replica_exchange.py index 3c2d117..b5e4ae0 100644 --- a/mcpy/ensembles/batched_replica_exchange.py +++ b/mcpy/ensembles/batched_replica_exchange.py @@ -36,6 +36,9 @@ logger = logging.getLogger(__name__) +# Below this many attempts an all-accept / all-reject tally says nothing. +_SWAP_DEGENERACY_MIN_ATTEMPTS = 20 + class BatchedReplicaExchange: def __init__( @@ -224,9 +227,12 @@ def _batched_single_move(self, active: List[int]) -> None: result if isinstance(result, tuple) else (result, 0, None) ) if atoms_new is False or atoms_new is None: - # Move couldn't propose — atoms unchanged, snapshot harmless. - # Identity check, not truthiness: an empty Atoms (last atom - # deleted) is falsy but is a real proposal. + # Move couldn't propose. Identity check, not truthiness: an + # empty Atoms (last atom deleted) is falsy but is a real + # proposal. Restore rather than trust the sentinel to mean + # "nothing was touched" (see GrandCanonicalEnsemble). + r.atoms.arrays = snapshots[i] + r.atoms.set_constraint(constraint_snapshots[i]) continue if atoms_new is not r.atoms: raise RuntimeError( @@ -284,21 +290,37 @@ def _attempt_exchanges(self) -> None: self.exchange_successes[j] += 1 def _accept_swap(self, i: int, j: int) -> bool: - """Temperature-RE Metropolis: P = exp((β_j - β_i)(Φ_j - Φ_i)). + """Replica-exchange Metropolis for a temperature *or* a μ ladder. Replicas here are grand-canonical (fixed μ, fluctuating N), so configs - at different temperatures must be compared through the grand potential - Φ = E - Σ_s μ_s N_s rather than bare energy. With no chemical potential - this reduces to the standard energy-only swap. + must be compared through the grand potential Φ = E - Σ_s μ_s N_s rather + than bare energy. Exchanging config X (slot i) with config Y (slot j) + changes the joint weight by + + ln(w'/w) = β_i Φ_X^(i) + β_j Φ_Y^(j) - β_i Φ_Y^(i) - β_j Φ_X^(j) + + where Φ_Z^(k) is configuration Z scored with slot k's chemical + potentials. The cross terms are what makes this work on both ladders: + with a shared μ it collapses to the familiar (β_j - β_i)(Φ_j - Φ_i), + but on a μ ladder (one shared temperature) that shortened form is + identically zero and would accept every swap, destroying the ladder. + With shared β it reduces to β Σ_s (μ_i,s - μ_j,s)(N_Y,s - N_X,s), the + rule ``ReplicaExchange._exchange_prob_mu`` implements for MPI. """ ri, rj = self.replicas[i], self.replicas[j] beta_i, beta_j = ri.units.beta, rj.units.beta - phi_i, phi_j = self._grand_potential(ri), self._grand_potential(rj) - delta = (beta_j - beta_i) * (phi_j - phi_i) - p = min(1.0, float(np.exp(delta))) + phi_ii = self._grand_potential(ri) # config X, μ of i + phi_jj = self._grand_potential(rj) # config Y, μ of j + phi_ji = ri._minimum_score(rj.atoms, rj.E_old) # config Y, μ of i + phi_ij = rj._minimum_score(ri.atoms, ri.E_old) # config X, μ of j + delta = (beta_i * phi_ii + beta_j * phi_jj + - beta_i * phi_ji - beta_j * phi_ij) + # exp() overflows to inf for a large positive delta; short-circuit. + p = 1.0 if delta >= 0.0 else float(np.exp(delta)) self.logger.debug( - "swap %d<->%d: beta_i=%.3e beta_j=%.3e Phi_i=%.3f Phi_j=%.3f delta=%.3f p=%.3f", - i, j, beta_i, beta_j, phi_i, phi_j, delta, p, + "swap %d<->%d: beta_i=%.3e beta_j=%.3e Phi_ii=%.3f Phi_jj=%.3f " + "Phi_ji=%.3f Phi_ij=%.3f delta=%.3f p=%.3f", + i, j, beta_i, beta_j, phi_ii, phi_jj, phi_ji, phi_ij, delta, p, ) return self.rng.get_uniform() < p @@ -364,8 +386,42 @@ def _log_status(self, step: int) -> None: self.logger.info('RE %d/%d | N: %s | E(eV): %s | swap acc %s', step, self.gcmc_steps, n, e, swap) + def _warn_if_ladder_degenerate(self) -> None: + """Flag a ladder that stopped behaving like one. + + A ladder accepting *every* swap is not exchanging information between + rungs, it is averaging them together; one accepting *none* leaves each + replica sampling in isolation. Neither shows up as an error, and the + all-accept case in particular yields a plausible-looking but + μ-independent isotherm, so it needs saying out loud. + + Checked once on the whole-run tally rather than live: a ladder whose + replicas have not differentiated yet (identical starting configs, an + adsorbate that has not begun to adsorb) legitimately accepts every + early swap, and a warning that fires during equilibration would be + noise. The trade is that this is a post-mortem. + """ + attempts, successes = sum(self.exchange_attempts), sum(self.exchange_successes) + if attempts < _SWAP_DEGENERACY_MIN_ATTEMPTS: + return + if successes == attempts: + reason = ('every swap was accepted -- neighbouring rungs are ' + 'indistinguishable, so the ladder averages over them ' + 'and per-rung results lose their T/mu dependence') + elif successes == 0: + reason = ('no swap was ever accepted -- rungs are too far apart ' + 'to exchange, so each replica sampled in isolation') + else: + return + self.logger.warning( + 'replica-exchange ladder looks degenerate: %d/%d swaps accepted. ' + 'Widen or tighten the %s spacing.', + successes, attempts, 'mu' if self.mus is not None else 'temperature') + self.logger.warning(' %s', reason) + def _log_summary(self) -> None: """Final consolidated summary: per-replica move acceptances.""" + self._warn_if_ladder_degenerate() for i, r in enumerate(self.replicas): ratios = ' '.join( f'{x:.0%}' if x == x else 'n/a' diff --git a/mcpy/ensembles/canonical_ensemble.py b/mcpy/ensembles/canonical_ensemble.py index 77f1a9a..4465bac 100644 --- a/mcpy/ensembles/canonical_ensemble.py +++ b/mcpy/ensembles/canonical_ensemble.py @@ -95,12 +95,8 @@ def set_state(self, state): # config swap, without relying on the pickled .info surviving MPI. self.atoms.info.setdefault("key_value_pairs", {}) self.atoms.info["key_value_pairs"]["potential_energy"] = state["energy"] - if "step" in state: - self._step = state["step"] - if "exchange_attempts" in state: - self.exchange_attempts = state["exchange_attempts"] - if "exchange_successes" in state: - self.exchange_successes = state["exchange_successes"] + # Step count and exchange statistics stay with the slot, not the + # config -- see GrandCanonicalEnsemble.set_state. def _acceptance_condition(self, potential_diff: float) -> bool: if potential_diff <= 0: diff --git a/mcpy/ensembles/grand_canonical_ensemble.py b/mcpy/ensembles/grand_canonical_ensemble.py index ca8345b..fac1b6d 100644 --- a/mcpy/ensembles/grand_canonical_ensemble.py +++ b/mcpy/ensembles/grand_canonical_ensemble.py @@ -98,13 +98,11 @@ def set_state(self, state: Dict[str, any]) -> None: self.atoms = state["atoms"] self.E_old = state["energy"] self.n_atoms = state["n_atoms"] - # Restore optional bookkeeping if present (used by restart, not by RE). - if "step" in state: - self._step = state["step"] - if "exchange_attempts" in state: - self.exchange_attempts = state["exchange_attempts"] - if "exchange_successes" in state: - self.exchange_successes = state["exchange_successes"] + # Only the configuration travels. Step count and exchange statistics + # describe the *slot*, not the config: ``ReplicaExchange`` passes a + # full ``get_state()`` dict here on every accepted swap, so restoring + # them would trade the two ranks' swap tallies back and forth and make + # the per-rank "Accepted Exchange (%)" column meaningless. # The configuration changed under the cells: refresh their free # volumes now, or the next insertion/deletion acceptance would use # the previous configuration's volume (GCMC only recalculates on @@ -143,12 +141,17 @@ def write_outfile(self, step: int = None, energy: float = None) -> None: :class:`BaseEnsemble` signature but ignored — GCMC always logs its own ``_step``/``E_old`` so the row matches the sampler state. """ - if self._outfile is None or self._outfile_handle is None: - return if self._last_logged_step == self._step: return # already wrote this step (e.g. finalize_run after a triggered write) acceptance_ratios = self.move_selector.interval_ratios() + # Reset before the disabled-outfile bail-out: the counters are + # per-write-interval, so skipping the reset when there is nothing to + # write turns interval_ratios() into total_ratios() for the rest of + # the run. self.move_selector.reset_counters() + if self._outfile is None or self._outfile_handle is None: + self._last_logged_step = self._step + return ratio_str = ", ".join( f"{r * 100:.1f}%" if not np.isnan(r) else "N/A" for r in acceptance_ratios @@ -246,11 +249,18 @@ def do_gcmc_step(self) -> None: atoms_new, delta_particles, species = self.move_selector.do_trial_move(atoms) if atoms_new is False or atoms_new is None: - # Move couldn't be proposed (e.g. empty cell). The move did - # not mutate ``atoms``; MoveSelector already recorded the - # failure so it won't depress the acceptance ratio. Identity - # check, not truthiness: an empty Atoms (last atom deleted) - # is falsy but is a real proposal that must be scored. + # Move couldn't be proposed (e.g. empty cell). MoveSelector + # already recorded the failure so it won't depress the + # acceptance ratio. Identity check, not truthiness: an empty + # Atoms (last atom deleted) is falsy but is a real proposal + # that must be scored. + # + # Restore the snapshot rather than trusting the sentinel to + # mean "nothing was touched": a move that mutates before + # bailing would otherwise leave the configuration changed + # while ``E_old`` still describes the previous one. + atoms.arrays = saved_arrays + atoms.set_constraint(saved_constraints) continue if atoms_new is not atoms: diff --git a/mcpy/moves/molecule_displacement_move.py b/mcpy/moves/molecule_displacement_move.py index c154421..500e7d0 100644 --- a/mcpy/moves/molecule_displacement_move.py +++ b/mcpy/moves/molecule_displacement_move.py @@ -61,11 +61,14 @@ def do_trial_move(self, atoms) -> Atoms: # Reject displacements that carry the COM out of the region: the # reverse proposal would have zero probability (the molecule stops # being a candidate), so crossing would be a one-way door that - # strands molecules outside the grand-canonical bookkeeping. + # strands molecules outside the grand-canonical bookkeeping. The + # region tested is the exchangeable one -- the same predicate that + # decides candidacy in ``find_molecules`` -- so the door stays + # two-way for every cell type. new_com = com + shift if np.any(np.asarray(atoms.pbc)): new_com = wrap_positions([new_com], atoms.cell, pbc=atoms.pbc)[0] - if not self.cell.is_point_inside(new_com): + if not self.cell.is_point_exchangeable(new_com): return False, 0, self.name if self.max_angle is None: rotation = random_rotation_matrix(self.rng) diff --git a/mcpy/moves/molecule_utils.py b/mcpy/moves/molecule_utils.py index 8f172a2..2f2a176 100644 --- a/mcpy/moves/molecule_utils.py +++ b/mcpy/moves/molecule_utils.py @@ -42,8 +42,10 @@ def find_molecules(atoms, template_symbols, cell=None): ``template_symbols`` is the sorted list of the template's chemical symbols. When ``cell`` is given, only molecules whose center of mass - lies inside it (``cell.is_point_inside``) are returned; ``None`` skips - the spatial filter (used by the grand-potential bookkeeping). + the cell counts as exchangeable (``cell.is_point_exchangeable``) are + returned; ``None`` skips the spatial filter (used by the grand-potential + bookkeeping). The exchangeable region is not always the proposal region: + see :meth:`mcpy.cell.Cell.is_point_exchangeable`. """ ids = atoms.arrays.get('molecule_id') if ids is None: @@ -56,7 +58,8 @@ def find_molecules(atoms, template_symbols, cell=None): members = np.where(ids == mid)[0] if sorted(symbols[members]) != list(template_symbols): continue - if cell is not None and not cell.is_point_inside(molecule_com(atoms, members)): + if cell is not None and not cell.is_point_exchangeable( + molecule_com(atoms, members)): continue groups.append(members) return groups diff --git a/mcpy/moves/move_selector.py b/mcpy/moves/move_selector.py index 290502d..b4db54e 100644 --- a/mcpy/moves/move_selector.py +++ b/mcpy/moves/move_selector.py @@ -1,9 +1,26 @@ +import re import warnings import numpy as np from mcpy.utils import RandomNumberGenerator +def _abbreviate(class_name): + """Short but still distinct label for a move class. + + A plain three-character slice collapses ``MoleculeInsertionMove``, + ``MoleculeDeletionMove`` and ``MoleculeDisplacementMove`` to the same + 'Mol', which makes the outfile's acceptance-ratio header unreadable for + molecular runs. Split the CamelCase name, drop the trailing 'Move', and + keep three letters per remaining word: 'MolIns', 'MolDel', 'MolDis'. + Single-word moves are unchanged ('Ins', 'Del', 'Dis', 'Per'). + """ + words = re.findall('[A-Z][a-z0-9]*', class_name) or [class_name] + if len(words) > 1 and words[-1] == 'Move': + words = words[:-1] + return ''.join(word[:3] for word in words) + + class MoveSelector: """Class used to randomly select a Monte Carlo Move from a list of Moves. @@ -35,7 +52,7 @@ class MoveSelector: def __init__(self, probabilities, move_list, seed=None, n_moves=None): assert len(probabilities) == len(move_list) self.move_list = move_list - self.move_list_names = [move.__class__.__name__[:3] for move in move_list] + self.move_list_names = [_abbreviate(move.__class__.__name__) for move in move_list] if n_moves is None: n_moves = max(1, round(sum(probabilities))) if n_moves < 1: diff --git a/mcpy/moves/permutation_move.py b/mcpy/moves/permutation_move.py index 6bfda42..61ebb7b 100644 --- a/mcpy/moves/permutation_move.py +++ b/mcpy/moves/permutation_move.py @@ -33,12 +33,18 @@ def do_trial_move(self, atoms): """ nums = atoms.arrays['numbers'] symbols = np.asarray(atoms.get_chemical_symbols()) + # Resolve the usable species once, up front. A swap preserves every + # species count, so a species absent now stays absent for the whole + # trial: bailing out from inside the loop would leave the earlier + # swaps applied while returning the "couldn't propose" sentinel, which + # the ensembles read as "the atoms were not touched". + present = [s for s in self.species if np.any(symbols == s)] + if len(present) < 2: + return False, 0, 'X' for _ in range(self.n_swaps): - species_pair = self.rng.random.sample(self.species, 2) + species_pair = self.rng.random.sample(present, 2) indices_a = np.where(symbols == species_pair[0])[0] indices_b = np.where(symbols == species_pair[1])[0] - if len(indices_a) == 0 or len(indices_b) == 0: - return False, 0, 'X' i = int(self.rng.random.choice(indices_a)) j = int(self.rng.random.choice(indices_b)) nums[i], nums[j] = nums[j], nums[i] diff --git a/tests/test_audit_regressions.py b/tests/test_audit_regressions.py index 20ab53a..ac2bfa5 100644 --- a/tests/test_audit_regressions.py +++ b/tests/test_audit_regressions.py @@ -398,3 +398,137 @@ def test_import_does_not_mutate_matplotlib_rcparams(): "matplotlib.rcParams['figure.dpi']\n" ) subprocess.run([sys.executable, '-c', code], check=True) + + +# -------------------------------------------------------------------------- +# 2026-07-27 review +# -------------------------------------------------------------------------- + +# A move that mutates and then reports failure must not corrupt the sampler +# (bug: the ensembles trusted the sentinel to mean "atoms untouched" and +# skipped the rollback, leaving the config out of step with E_old) + +class _DirtyBailingMove(BaseMove): + """Contract violator: mutates the atoms, then returns the ``False`` + 'could not propose' sentinel.""" + + def __init__(self, seed=1): + super().__init__(NullCell(), ['X'], seed) + + def do_trial_move(self, atoms): + atoms.positions += 1.0 + return False, 0, 'X' + + +def _h2(): + return Atoms('H2', positions=[[0.0, 0.0, 0.0], [2.0, 0.0, 0.0]], + cell=[10.0, 10.0, 10.0]) + + +def test_gcmc_rolls_back_a_move_that_mutates_then_reports_failure(): + atoms = _h2() + before = atoms.positions.copy() + g = _gcmc(atoms, [NullCell()], MoveSelector([1], [_DirtyBailingMove()]), + {'H': 0.0}) + energy_before = g.E_old + g.do_gcmc_step() + np.testing.assert_allclose(g.atoms.positions, before) + assert g.E_old == energy_before + + +def test_batched_re_rolls_back_a_move_that_mutates_then_reports_failure(): + from mcpy.ensembles.batched_replica_exchange import BatchedReplicaExchange + + atoms = _h2() + before = atoms.positions.copy() + replica = _gcmc(atoms, [NullCell()], + MoveSelector([1], [_DirtyBailingMove()]), {'H': 0.0}) + + class _BatchCalc(StubCalc): + def get_potential_energies(self, atoms_list): + return np.array([self.get_potential_energy(a) for a in atoms_list]) + + pt = BatchedReplicaExchange.__new__(BatchedReplicaExchange) + pt.replicas = [replica] + pt.calculator = _BatchCalc() + pt._batched_single_move([0]) + np.testing.assert_allclose(replica.atoms.positions, before) + + +# Replica exchange swaps configurations, not slot bookkeeping (bug: get_state +# emitted step/exchange counters and set_state restored them, so accepted +# swaps traded the two ranks' tallies) + +def test_gcmc_set_state_keeps_local_step_and_exchange_counters(): + g = _gcmc(_h2(), [NullCell()], MoveSelector([1], [_DirtyBailingMove()]), + {'H': 0.0}) + g._step, g.exchange_attempts, g.exchange_successes = 40, 8, 3 + + partner = _gcmc(_h2(), [NullCell()], + MoveSelector([1], [_DirtyBailingMove()]), {'H': 0.0}) + partner._step, partner.exchange_attempts, partner.exchange_successes = 40, 8, 7 + g.set_state(partner.get_state()) + + assert (g._step, g.exchange_attempts, g.exchange_successes) == (40, 8, 3) + assert g.atoms is partner.atoms # the configuration still travels + + +def test_canonical_set_state_keeps_local_step_and_exchange_counters(tmp_path): + mc = _canonical(tmp_path) + mc._step, mc.exchange_attempts, mc.exchange_successes = 12, 4, 1 + mc.set_state({'atoms': _h2(), 'energy': -1.0, 'step': 99, + 'exchange_attempts': 40, 'exchange_successes': 39}) + assert (mc._step, mc.exchange_attempts, mc.exchange_successes) == (12, 4, 1) + assert mc._current_energy == -1.0 + + +# Per-interval acceptance counters are cleared even with the outfile disabled +# (bug: write_outfile returned before reset_counters, so interval_ratios() +# silently degraded into total_ratios()) + +def test_interval_ratios_reset_with_outfile_disabled(): + atoms = _h2() + g = _gcmc(atoms, [NullCell()], MoveSelector([1], [DisplacementMove( + species=['H'], seed=5, max_displacement=0.1)]), {'H': 0.0}) + assert g._outfile is None + g.do_gcmc_step() + assert g.move_selector.move_counter == [1] + g.write_outfile() + assert g.move_selector.move_counter == [0] + assert g.move_selector.move_counter_total == [1] # cumulative untouched + + +# Distinct outfile labels for the molecule moves (bug: a 3-character slice +# collapsed every Molecule*Move to 'Mol') + +def test_move_labels_distinguish_the_molecule_moves(): + from mcpy.moves.move_selector import _abbreviate + + assert _abbreviate('MoleculeInsertionMove') == 'MolIns' + assert _abbreviate('MoleculeDeletionMove') == 'MolDel' + assert _abbreviate('MoleculeDisplacementMove') == 'MolDis' + labels = [_abbreviate(n) for n in ('MoleculeInsertionMove', + 'MoleculeDeletionMove', + 'MoleculeDisplacementMove')] + assert len(set(labels)) == 3 + + +def test_move_labels_unchanged_for_single_word_moves(): + """Atomic runs keep the labels their existing outfiles were written with.""" + from mcpy.moves.move_selector import _abbreviate + + for name in ('InsertionMove', 'DeletionMove', 'DisplacementMove', + 'PermutationMove', 'ShakeMove', 'BrownianMove'): + assert _abbreviate(name) == name[:3] + + +# A missing exclusion radius names the species (bug: bare KeyError from inside +# the free-volume sampler) + +def test_missing_species_radius_raises_a_named_error(): + atoms = Atoms('Cu2O', positions=[[0, 0, 0], [2, 0, 0], [0, 2, 0]], + cell=[10.0, 10.0, 10.0]) + cell = CustomCell(atoms, custom_height=5.0, bottom_z=0.0, + species_radii={'Cu': 2.4}, mc_sample_points=64, seed=0) + with pytest.raises(ValueError, match=r"no radius for \['O'\]"): + cell.calculate_volume(atoms) diff --git a/tests/test_batched_re_acceptance.py b/tests/test_batched_re_acceptance.py index 8824dba..479b93d 100644 --- a/tests/test_batched_re_acceptance.py +++ b/tests/test_batched_re_acceptance.py @@ -7,6 +7,9 @@ """ import logging +import numpy as np +import pytest + from mcpy.ensembles.batched_replica_exchange import BatchedReplicaExchange @@ -85,6 +88,123 @@ def test_accept_swap_uses_grand_potential_not_bare_energy(): assert re._accept_swap(0, 1) is True +# -------------------------------------------------------------------------- +# mu-ladder swaps (bug: the temperature-only criterion collapses to delta=0 +# when both replicas share a temperature, accepting every swap unconditionally) +# -------------------------------------------------------------------------- + +_BETA_300K = 1.0 / (8.617333e-5 * 300.0) + + +def _mu_ladder(u): + """Two same-temperature replicas differing only in mu and N.""" + re = _bare_re() + ri = _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -5.0}, + symbols=['O'] * 10) + rj = _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -3.0}, + symbols=['O'] * 60) + re.replicas = [ri, rj] + re.rng = _FakeRng(u) + return re + + +def test_mu_ladder_swap_is_not_always_accepted(): + # beta_i == beta_j, so (beta_j - beta_i)(Phi_j - Phi_i) is identically 0 + # and the old rule returned p = 1 for every draw. The cross terms give + # exp(beta (mu_i - mu_j)(N_j - N_i)) = exp(-3868) here: always reject. + assert not _mu_ladder(u=1e-12)._accept_swap(0, 1) + + +def _partial_mu_ladder(): + """mu-ladder tuned to an acceptance well inside (0, 1), so that both the + 'just below p' and the 'just above p' draw are legal uniforms -- a p + pinned at 1 would let the buggy always-accept rule pass either way. + Energies differ too: they must cancel out of a same-temperature swap.""" + beta, mu_i, n_i, n_j = _BETA_300K, -5.0, 10, 12 + mu_j = mu_i + 0.5 / (beta * (n_j - n_i)) # -> delta = -0.5 + re = _bare_re() + re.replicas = [ + _FakeReplica(beta=beta, energy=-100.0, mu={'O': mu_i}, + symbols=['O'] * n_i), + _FakeReplica(beta=beta, energy=-42.0, mu={'O': mu_j}, + symbols=['O'] * n_j), + ] + p = float(np.exp(beta * (mu_i - mu_j) * (n_j - n_i))) + assert p == pytest.approx(np.exp(-0.5)) + return re, p + + +def test_mu_ladder_swap_matches_analytic_probability(): + re, p = _partial_mu_ladder() + re.rng = _FakeRng(p * 0.99) + assert re._accept_swap(0, 1) + re.rng = _FakeRng(p * 1.01) + assert not re._accept_swap(0, 1) + + +def test_mu_ladder_swap_ignores_the_energy_difference(): + """At one temperature the two configs' energies cancel exactly, so the + swap probability must not move when either replica's energy changes.""" + re, p = _partial_mu_ladder() + re.rng = _FakeRng(p * 0.99) + re.replicas[0].E_old += 37.0 + re.replicas[1].E_old -= 12.5 + assert re._accept_swap(0, 1) + re.rng = _FakeRng(p * 1.01) + assert not re._accept_swap(0, 1) + + +def test_mu_ladder_swap_toward_favoured_replica_always_accepted(): + # Reverse the population imbalance: the same expression is now positive, + # so the swap is downhill in the joint weight and must always be taken. + re = _bare_re() + re.replicas = [ + _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -3.0}, + symbols=['O'] * 10), + _FakeReplica(beta=_BETA_300K, energy=-100.0, mu={'O': -5.0}, + symbols=['O'] * 60), + ] + re.rng = _FakeRng(1.0 - 1e-12) + assert re._accept_swap(0, 1) # p clamped to 1 without an exp() overflow + + +def test_swap_probability_is_symmetric_in_slot_order(): + """Both orderings must agree, or the pairing loop's (i, j) choice would + silently bias the ladder.""" + re, p = _partial_mu_ladder() + for u in (p * 0.99, p * 1.01): + re.rng = _FakeRng(u) + assert re._accept_swap(0, 1) == re._accept_swap(1, 0) + + +def test_temperature_ladder_swap_unchanged_by_the_mu_fix(): + """With a shared mu the general form must reduce exactly to the old + (beta_j - beta_i)(Phi_j - Phi_i).""" + ri = _FakeReplica(beta=1.0, energy=10.0, mu={'Ag': 1.0}, symbols=['Ag'] * 5) + rj = _FakeReplica(beta=0.5, energy=20.0, mu={'Ag': 1.0}, symbols=['Ag'] * 18) + re = _bare_re() + re.replicas = [ri, rj] + # Phi_i = 5, Phi_j = 2 -> delta = (0.5 - 1.0)(2 - 5) = 1.5 -> p = 1. + legacy_delta = (rj.units.beta - ri.units.beta) * ( + re._grand_potential(rj) - re._grand_potential(ri)) + assert legacy_delta == pytest.approx(1.5) + re.rng = _FakeRng(1.0 - 1e-12) + assert re._accept_swap(0, 1) + + +def test_temperature_ladder_uphill_swap_matches_legacy_probability(): + ri = _FakeReplica(beta=1.0, energy=10.0, mu={'Ag': 1.0}, symbols=['Ag'] * 18) + rj = _FakeReplica(beta=0.5, energy=20.0, mu={'Ag': 1.0}, symbols=['Ag'] * 5) + re = _bare_re() + re.replicas = [ri, rj] + # Phi_i = -8, Phi_j = 15 -> delta = (0.5 - 1.0)(15 - -8) = -11.5. + p = float(np.exp(-11.5)) + re.rng = _FakeRng(p * 0.99) + assert re._accept_swap(0, 1) + re.rng = _FakeRng(p * 1.01) + assert not re._accept_swap(0, 1) + + def test_consolidated_status_line(caplog): """One console line covers all replicas; per-replica detail stays in files.""" import logging @@ -107,3 +227,58 @@ def test_consolidated_status_line(caplog): assert len(caplog.records) == 1 msg = caplog.records[0].getMessage() assert 'RE 30/100' in msg and '42' in msg and '44' in msg and '25%' in msg + + +# -------------------------------------------------------------------------- +# Degenerate-ladder warning. A 100%-accept ladder sat unremarked in the +# "Accepted Exchange (%)" column of six consecutive production runs before the +# mu-ladder bug above was found; an all-accept tally must not stay silent. +# -------------------------------------------------------------------------- + +def _tally_re(attempts, successes, mus=None): + pt = BatchedReplicaExchange.__new__(BatchedReplicaExchange) + pt.logger = logging.getLogger('mcpy.ensembles.batched_replica_exchange') + pt.exchange_attempts = attempts + pt.exchange_successes = successes + pt.mus = mus + return pt + + +@pytest.mark.parametrize('successes, expect', [ + ([40, 40], True), # every swap accepted + ([0, 0], True), # no swap ever accepted + ([12, 12], False), # healthy partial acceptance +]) +def test_degenerate_ladder_warning(caplog, successes, expect): + pt = _tally_re([40, 40], successes) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + warned = any(r.levelno == logging.WARNING for r in caplog.records) + assert warned is expect + + +def test_degenerate_ladder_warning_needs_enough_attempts(caplog): + """Two all-accept swaps prove nothing; don't cry wolf.""" + pt = _tally_re([1, 1], [1, 1]) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + assert not caplog.records + + +def test_degenerate_ladder_warning_names_the_ladder_kind(caplog): + pt = _tally_re([40, 40], [40, 40], mus=[{'CO': -1.0}, {'CO': -0.5}]) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + text = ' '.join(r.getMessage() for r in caplog.records) + assert '80/80' in text and 'mu spacing' in text + + caplog.clear() + pt = _tally_re([40, 40], [40, 40], mus=None) + with caplog.at_level(logging.WARNING, + logger='mcpy.ensembles.batched_replica_exchange'): + pt._warn_if_ladder_degenerate() + assert 'temperature spacing' in ' '.join( + r.getMessage() for r in caplog.records) diff --git a/tests/test_compound_moves.py b/tests/test_compound_moves.py index 44a9f44..e1447ce 100644 --- a/tests/test_compound_moves.py +++ b/tests/test_compound_moves.py @@ -78,3 +78,46 @@ def test_displacement_n_steps_exceeds_movable_raises(): move = DisplacementMove(species=['Au'], seed=1, n_steps=5) with pytest.raises(ValueError): move.do_trial_move(atoms) + + +# -------------------------------------------------------------------------- +# A trial must be all-or-nothing (bug: a multi-swap trial could apply k swaps +# and then report the "couldn't propose" sentinel, which the ensembles read as +# "the atoms were not touched" -- so the config changed while E_old did not) +# -------------------------------------------------------------------------- + +def test_permutation_absent_species_does_not_abort_a_multi_swap_trial(): + """Declaring a species that is not in the system must not turn a later + iteration into a mid-trial bail-out. Species counts are swap-invariant, so + the usable pair is resolved once, up front.""" + for seed in range(200): + atoms = _balanced_alloy() + before = atoms.get_atomic_numbers().copy() + move = PermutationMove(species=['Au', 'Pt', 'Ag'], seed=seed, n_swaps=3) + result, delta, _ = move.do_trial_move(atoms) + assert result is atoms, f'seed={seed} bailed out mid-trial' + assert delta == 0 + assert not np.array_equal(atoms.get_atomic_numbers(), before) + + +def test_permutation_without_two_present_species_is_a_clean_no_op(): + """The one case that genuinely cannot propose must leave the atoms + byte-identical, since the ensembles skip the rollback on that path.""" + atoms = Atoms('Au4', positions=[[0, 0, 0], [2, 0, 0], [0, 2, 0], [0, 0, 2]]) + before = atoms.get_atomic_numbers().copy() + move = PermutationMove(species=['Au', 'Pt'], seed=1, n_swaps=3) + result, delta, name = move.do_trial_move(atoms) + assert result is False + assert (delta, name) == (0, 'X') + np.testing.assert_array_equal(atoms.get_atomic_numbers(), before) + + +def test_permutation_stream_unchanged_when_every_species_is_present(): + """Filtering to the present species must be a no-op for a healthy setup, + or previously published runs would no longer reproduce.""" + results = [] + for species in (['Au', 'Pt'], ['Au', 'Pt']): + atoms = _balanced_alloy() + PermutationMove(species=species, seed=11, n_swaps=5).do_trial_move(atoms) + results.append(atoms.get_atomic_numbers().copy()) + np.testing.assert_array_equal(*results) diff --git a/tests/test_custom_cell_region.py b/tests/test_custom_cell_region.py index dfca434..f751e37 100644 --- a/tests/test_custom_cell_region.py +++ b/tests/test_custom_cell_region.py @@ -64,3 +64,84 @@ def test_subsurface_oxygen_is_excluded_from_deletion_candidates(): assert total_o == 2 # both O are present in the structure assert len(counted) == 1 # only the above-floor O is a deletion candidate + + +# -------------------------------------------------------------------------- +# Exchangeable region vs proposal region (bug: molecules that desorbed above +# the cell top stopped being deletion candidates AND dropped out of the +# per-species de Broglie count, driving runaway insertion) +# -------------------------------------------------------------------------- + +def _point_at(cell, frac_z): + """A point in the middle of the footprint at fractional height ``frac_z`` + (``> 1`` is above the cell top).""" + return np.array([0.5, 0.5, frac_z]) @ cell.dimensions + cell.offset + + +def test_point_above_the_cell_top_is_exchangeable_but_not_inside(): + """The two predicates are deliberately different: ``is_point_inside`` + bounds the *proposal* region, ``is_point_exchangeable`` the region the + reservoir can take molecules back from -- the same asymmetry + ``get_atoms_specie_inside_cell`` already applies to single atoms.""" + _, cell = make_cell() + escaped = _point_at(cell, 1.4) + assert not cell.is_point_inside(escaped) + assert cell.is_point_exchangeable(escaped) + + +def test_point_below_the_cell_floor_is_neither(): + _, cell = make_cell() + buried = _point_at(cell, -0.2) + assert not cell.is_point_inside(buried) + assert not cell.is_point_exchangeable(buried) + + +def test_point_outside_the_xy_footprint_is_neither(): + _, cell = make_cell() + outside = np.array([1.4, 0.5, 0.5]) @ cell.dimensions + cell.offset + assert not cell.is_point_inside(outside) + assert not cell.is_point_exchangeable(outside) + + +def test_desorbed_molecule_stays_a_deletion_candidate(): + """A CO whose center of mass drifts above the cell top must remain + findable, exactly as its individual C and O atoms already do. Otherwise it + can never be deleted and its absence from ``last_exchange_count`` inflates + V/((N+1)Lambda^3) -- the runaway insertion mode documented in + docs/gcmc_acceptance_convention.rst. + """ + from ase import Atoms + + from mcpy.moves.molecule_utils import find_molecules + + atoms, cell = make_cell() + atoms.new_array('molecule_id', np.full(len(atoms), -1, dtype=int)) + for mol_id, frac_z in enumerate((0.5, 1.4)): # inside, and above the top + frag = Atoms('CO', positions=[[0, 0, 0], [0, 0, 1.13]]) + frag.positions += _point_at(cell, frac_z) + frag.new_array('molecule_id', np.full(len(frag), mol_id, dtype=int)) + atoms += frag + + template = sorted(['C', 'O']) + assert len(find_molecules(atoms, template)) == 2 + assert len(find_molecules(atoms, template, cell)) == 2 + # The atomic path has always seen all four member atoms; the molecular + # path must not disagree with it. + assert len(cell.get_atoms_specie_inside_cell(atoms, ['C', 'O'])) == 4 + + +def test_buried_molecule_is_still_excluded(): + """The floor exclusion is intentional (buried species are kept), so the + fix above must not open that side of the region too.""" + from ase import Atoms + + from mcpy.moves.molecule_utils import find_molecules + + atoms, cell = make_cell() + atoms.new_array('molecule_id', np.full(len(atoms), -1, dtype=int)) + frag = Atoms('CO', positions=[[0, 0, 0], [0, 0, 1.13]]) + frag.positions += _point_at(cell, -0.5) + frag.new_array('molecule_id', np.zeros(len(frag), dtype=int)) + atoms += frag + + assert find_molecules(atoms, sorted(['C', 'O']), cell) == [] diff --git a/tests/test_molecule_moves.py b/tests/test_molecule_moves.py index 3eaf3dc..633a0c5 100644 --- a/tests/test_molecule_moves.py +++ b/tests/test_molecule_moves.py @@ -86,7 +86,7 @@ def test_find_molecules_missing_array_and_cell_filter(): assert find_molecules(plain, ['H', 'H']) == [] class _NowhereCell: - def is_point_inside(self, point): + def is_point_exchangeable(self, point): return False atoms = _water_box() @@ -448,13 +448,13 @@ def test_write_xyz_molecule_id_roundtrip(tmp_path): def test_find_molecules_queries_com_point(): - """The point handed to cell.is_point_inside must be the molecule's + """The point handed to cell.is_point_exchangeable must be the molecule's center of mass, not e.g. the first member's position.""" atoms = _water_box() recorded = {} class _RecordingCell: - def is_point_inside(self, point): + def is_point_exchangeable(self, point): recorded['point'] = point return True From ea0328516e7a1a3bb776129e22a4b1aa330dd01e Mon Sep 17 00:00:00 2001 From: Riccardo Farris Date: Mon, 27 Jul 2026 16:14:26 +0200 Subject: [PATCH 2/5] Expose chunk_size in the CO/CuPd replica-exchange example An acceptance-equalized mu ladder for this system needs ~29 rungs (spacing 0.064 eV at the bare end down to 0.020 eV at high coverage, from dmu = 1/sqrt(beta dN/dmu) on the measured isotherm). At ~460 atoms per replica that is ~13k atoms in the relax batch, several times the whole-batch ceiling, so the example could not run a correctly spaced ladder at all. chunk_size ties peak memory to the largest chunk instead of the replica count. --- examples/re_gcmc_co_cupd_batched.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/examples/re_gcmc_co_cupd_batched.py b/examples/re_gcmc_co_cupd_batched.py index d531105..f444a6f 100644 --- a/examples/re_gcmc_co_cupd_batched.py +++ b/examples/re_gcmc_co_cupd_batched.py @@ -64,6 +64,13 @@ def parse_args(): p.add_argument('--mol-disp-angle', type=float, default=None, help='Max rotation angle of the rigid move (rad); ' 'None = full uniform rotation') + p.add_argument('--chunk-size', type=int, default=None, + help='Relax at most this many replicas per forward pass. ' + 'Peak GPU memory follows the largest chunk rather ' + 'than the replica count, which is what makes a wide ' + 'ladder fit: acceptance-equalized spacing needs ~29 ' + 'rungs here, and 29 x ~460 atoms is several times the ' + 'whole-batch ceiling. None relaxes the whole batch.') p.add_argument('--no-compile', action='store_true') p.add_argument('--seed', type=int, default=7, help='Master seed for moves, RE, and the Pd placement') @@ -111,6 +118,7 @@ def main(): steps=args.rel_steps, fmax=args.rel_fmax, compile_model=not args.no_compile, + chunk_size=args.chunk_size, ) e_co = calculator.get_potential_energy( From e20b1242d25238917f917ee70360377c9b5129d7 Mon Sep 17 00:00:00 2001 From: Riccardo Farris Date: Tue, 28 Jul 2026 09:34:03 +0200 Subject: [PATCH 3/5] Document how to space a replica-exchange chemical-potential ladder A mis-spaced mu ladder fails silently: no error, no warning, and an output that reads as plausible physics. Two things make it hard to catch, and both are now written down. First, the cumulative acceptance column cannot be used to judge a ladder. Every replica starts from the same configuration, so early swaps are free and inflate the tally permanently. A five-rung CO/CuPd run reported cumulative per-slot acceptances of 28.6 / 19.6 / 5.9 / 2.0 / 0.0 %, which reads as a ladder that merely mixes poorly at one end; in the run's second half three of its four pairs accepted exactly zero swaps and four replicas were independent single-mu chains. The halving test recovers the real rate from the log without extra instrumentation. Second, uniform spacing cannot work across a coverage range at all, because dN/dmu grows with coverage while dmu stays fixed. The fix is dmu = 1/sqrt(beta dN/dmu) from fluctuation-dissipation, with the rung count as the integral of sqrt(beta dN/dmu). For CO on Cu375Pd30 at 400 K over -1.8..-1.0 eV that is 29 rungs rather than 5, spacing 0.064 eV at the bare end down to 0.020 eV at high coverage. Measured second-half acceptance of that ladder: min 18 %, median 40 %, max 65 %, no dead pair, and a monotonic isotherm from 0.1 to 57 CO. The page carries a reference implementation, verified to reproduce that ladder, for whoever lands the mu counterpart of utils.ladder.geometric_temperatures. It also records the calibration caveat (the rule targets the mean-dN exponent, but the realised acceptance runs higher, so the rung count is a safe upper bound to trim against) and the two practical consequences of going wide: chunk_size becomes mandatory, and the run gets faster in absolute terms because a narrow ladder leaves the GPU idle. --- docs/index.rst | 1 + docs/replica_exchange_ladder_spacing.rst | 192 +++++++++++++++++++++++ 2 files changed, 193 insertions(+) create mode 100644 docs/replica_exchange_ladder_spacing.rst diff --git a/docs/index.rst b/docs/index.rst index 137c59e..7bdaf81 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -61,6 +61,7 @@ the hybrid GCMC method it implements (see :doc:`bibliography`). molecular_adsorbates phase_diagrams gcmc_acceptance_convention + replica_exchange_ladder_spacing .. toctree:: :caption: Tutorials diff --git a/docs/replica_exchange_ladder_spacing.rst b/docs/replica_exchange_ladder_spacing.rst new file mode 100644 index 0000000..54ea9d4 --- /dev/null +++ b/docs/replica_exchange_ladder_spacing.rst @@ -0,0 +1,192 @@ +Spacing a replica-exchange ladder +================================= + +A replica-exchange run can look perfectly healthy and still be doing nothing. +This note records how to tell, and how to space a chemical-potential ladder so +that it works. It is the ladder-design companion to +:doc:`gcmc_acceptance_convention`, and it exists because the failure modes here +are silent: no error, no warning, and output that reads as plausible physics. + +.. contents:: + :local: + :depth: 1 + + +Read second-half acceptance, never the cumulative column +-------------------------------------------------------- + +Every replica starts from the same configuration, so early swaps are free: two +rungs holding identical states always exchange. Those free successes are +recorded permanently in the cumulative tally, which therefore overstates the +ladder's health for the rest of the run -- badly, and for a long time. + +A five-rung CO/CuPd run reported cumulative per-slot acceptances of 28.6, 19.6, +5.9, 2.0 and 0.0 %. That reads as a ladder that mixes well at one end and +poorly at the other. In the run's *second half*, three of the four pairs +accepted **exactly zero** swaps: only pair 0-1 was still functioning, and the +other four replicas were independent single-:math:`\mu` chains. + +The test needs no extra instrumentation. Attempts accumulate linearly, so +between the run's midpoint and its end the attempt count doubles. If a +cumulative percentage *halves exactly* over that span, the numerator never +changed and there were no successes at all in the second half: + +.. code-block:: text + + second-half rate = 2 * cum_end - cum_mid + +Slots ``0`` and ``n-1`` each belong to exactly one pair, so their columns read +those two pairs directly; interior slots average the pair on either side. + +``BatchedReplicaExchange`` warns at the end of a run when the whole-run tally is +all-accept or all-reject, which catches a catastrophically mis-spaced ladder. +It cannot catch the case above, where a ladder starts healthy and dies as the +replicas differentiate. For that, run the halving test. + + +Why uniform spacing cannot work +------------------------------- + +For a :math:`\mu` ladder at a single temperature the swap exponent is + +.. math:: + + \ln \frac{w'}{w} = \beta\,(\mu_i - \mu_j)(N_j - N_i) + = -\beta\,\Delta\mu\,\Delta N . + +:math:`\Delta N` between neighbouring rungs is not a constant: it grows with +coverage, because :math:`dN/d\mu` does. A spacing chosen where the surface is +bare is therefore far too coarse once the surface fills. In the run above, at +400 K with :math:`\Delta\mu = 0.2` eV and :math:`\Delta N \approx 16`, the +exponent reaches :math:`29 \times 0.2 \times 16 \approx 93`, i.e. +:math:`p \sim 10^{-40}`. No amount of sampling recovers that pair. + + +The spacing rule +---------------- + +Grand-canonical fluctuation-dissipation gives the width of each rung's +:math:`N` distribution, + +.. math:: + + \sigma_N^2 = \frac{1}{\beta}\frac{dN}{d\mu} , + +and neighbouring rungs exchange when their distributions overlap, i.e. when +:math:`\beta\,\Delta\mu\,\sigma_N \sim 1`. Substituting: + +.. math:: + + \Delta\mu(\mu) = \frac{1}{\sqrt{\beta\,dN/d\mu}} , + \qquad + n_\text{rungs} = \int \sqrt{\beta\,\frac{dN}{d\mu}}\; d\mu . + +Use :math:`dN/d\mu` from a measured isotherm, **not** the per-rung +:math:`\sigma_N` observed in an unconverged run: a stuck rung reports a +:math:`\sigma_N` far below equilibrium and a drifting one far above, because the +"fluctuation" is drift. + +Histogram reweighting is not an alternative route to this. Reweighting requires +overlap between adjacent rungs, and a mis-spaced ladder has none -- that absence +is the disease being diagnosed. + +.. code-block:: python + + import numpy as np + + def acceptance_equalized_mu_ladder(mu, n_ads, temperature, mu_min, mu_max): + """Rung positions for a chemical-potential ladder of even acceptance. + + Places rungs where the cumulative integral of sqrt(beta dN/dmu) crosses + successive integers, so every neighbouring pair has comparable + distribution overlap and therefore comparable swap acceptance. + + Parameters + ---------- + mu, n_ads : array_like + A measured isotherm: chemical potentials and mean adsorbate counts. + temperature : float + Ladder temperature in K (all rungs share it). + mu_min, mu_max : float + Range to span. + + Returns + ------- + ndarray + Rung chemical potentials, ascending. Spacing widens where the + isotherm is flat and tightens where it is steep. + """ + beta = 1.0 / (8.617333e-5 * temperature) + mu, n_ads = np.asarray(mu, float), np.asarray(n_ads, float) + slope = np.diff(n_ads) / np.diff(mu) # dN/dmu per interval + midpoints = 0.5 * (mu[1:] + mu[:-1]) + + fine = np.linspace(mu_min, mu_max, 2001) + density = np.sqrt(beta * np.interp(fine, midpoints, slope)) + cumulative = np.concatenate([[0.0], np.cumsum( + 0.5 * (density[1:] + density[:-1]) * np.diff(fine))]) + return np.interp(np.arange(0.0, cumulative[-1], 1.0), cumulative, fine) + + +Calibration +----------- + +The rule above is derived from the exponent at the *mean* :math:`\Delta N`, but +the acceptance that matters is :math:`\langle \min(1, \cdot) \rangle` over both +rungs' fluctuations. A favourable fluctuation accepts outright, so the realised +acceptance is higher than the point estimate suggests -- for the CO/CuPd system +the rule aimed at ~25 % and delivered ~40 %. + +Treat the rule as a safe upper bound on rung count, then trim once a run has +measured the real acceptance. A ladder derived from an unconverged isotherm is +also a lower bound on the rungs eventually needed, since :math:`dN/d\mu` +steepens as coverage converges. Those two biases push in opposite directions; +re-derive from each run rather than trusting one calculation. + + +Worked example: CO on a CuPd nanoparticle +----------------------------------------- + +Cu\ :sub:`375`\ Pd\ :sub:`30` at 400 K, :math:`\Delta\mu` spanning -1.8 to -1.0 +eV, ``AlchemiFCalculator``. + +.. list-table:: + :header-rows: 1 + + * - Interval (eV) + - :math:`dN/d\mu` + - Spacing + - Rungs + * - -1.8 .. -1.6 + - 8.5 + - 0.064 eV + - 3 + * - -1.6 .. -1.4 + - 45.1 + - 0.028 eV + - 7 + * - -1.4 .. -1.2 + - 54.5 + - 0.025 eV + - 8 + * - -1.2 .. -1.0 + - 88.9 + - 0.020 eV + - 10 + +29 rungs, against the 5 that a uniform 0.2 eV ladder would use. Measured +second-half acceptance: minimum 18 %, median 40 %, maximum 65 %, with no dead +pair -- and a smooth, strictly monotonic isotherm from 0.1 to 57 CO. + +Two practical consequences of going wide: + +- **Chunking becomes mandatory.** 29 replicas of ~460 atoms is ~13k atoms, + several times the whole-batch relaxation ceiling. Pass ``chunk_size`` so peak + memory follows the largest chunk rather than the replica count. +- **It is faster, not slower, in absolute terms.** A narrow ladder leaves the + GPU idle. Five replicas ran at ~43 % utilisation; 29 replicas at + ``chunk_size=8`` reached 94-98 % and 15.6 GB of 32.6 GB. + +Once the ladder is healthy, convergence becomes the binding constraint rather +than the sampler. Chain runs with ``--init-dir`` pointing at the previous +output and watch the tail-half drift per rung, not just the acceptance. From 2ee6ab9f786650d8a7378b85a918790415d7f142 Mon Sep 17 00:00:00 2001 From: Riccardo Farris Date: Tue, 28 Jul 2026 13:20:37 +0200 Subject: [PATCH 4/5] Release v1.4.0 Minor rather than patch: the release adds a public cell predicate (`is_point_exchangeable`), adds a degenerate-ladder warning, and changes the outfile acceptance-ratio header for molecule moves, which any downstream parser of that column will notice. Also restores the `[1.3.0]` compare link that was omitted when 1.3.0 was cut. --- CHANGELOG.md | 23 +++++++++++++++++++++++ pyproject.toml | 2 +- 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 69d96a3..cdf226d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,27 @@ All notable changes to this project are documented here. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [1.4.0] - 2026-07-28 + +### Fixed +- **Batched replica exchange accepted every chemical-potential swap.** `BatchedReplicaExchange._accept_swap` implemented only the temperature-ladder criterion `(beta_j - beta_i)(Phi_j - Phi_i)`, which is identically zero when the replicas share a temperature: every mu-ladder swap was accepted with p = 1, configurations random-walked freely across the ladder, and no replica sampled its own mu. **Results produced by `BatchedReplicaExchange(..., mus=[...])` in 1.3.0 or earlier are invalid and need re-running**; temperature ladders were always correct, as was the MPI `ReplicaExchange`, which selected its mu criterion all along. The replacement carries the cross terms, `beta_i Phi_X^(i) + beta_j Phi_Y^(j) - beta_i Phi_Y^(i) - beta_j Phi_X^(j)`, where `Phi_Z^(k)` is a configuration scored with slot `k`'s chemical potentials; it reduces to the previous expression for a shared mu, to `beta (mu_i - mu_j)(N_j - N_i)` for a shared temperature, and additionally covers a joint (T, mu) ladder. The failure was not a statistical degradation: an EMT Ag(111)/Au ladder gives `N_Au = [0, 0, 0, 14]` with the fix and `[0, 0, 0, 0]` without it, because unconditional swapping averages the rungs together and erases the mu dependence entirely. +- **`PermutationMove` could mutate the configuration and then report that it had not.** With `n_swaps > 1`, drawing a species absent from the system in a later iteration returned the "could not propose" sentinel after earlier swaps had already been applied; both GCMC loops read that sentinel as "the atoms were not touched" and skipped the rollback, leaving the configuration out of step with the stored energy for the rest of the run. The usable species are now resolved once, before any mutation (species counts are swap-invariant, so a species absent at the start is absent for the whole trial), and both ensembles restore their pre-trial snapshot on the sentinel path rather than trusting the contract. +- Molecules whose center of mass drifted above a `CustomCell`'s top stopped being deletion candidates and dropped out of the per-species de Broglie count, inflating `V/((N+1)Lambda^3)` into the runaway insertion mode documented in `docs/gcmc_acceptance_convention.rst`. Molecule candidacy now goes through the new `is_point_exchangeable` predicate, which applies the same dropped z upper bound that `get_atoms_specie_inside_cell` already applied to single atoms; molecules below the cell floor stay excluded, as atoms do. +- `AlchemiCalculator(energy_only=True)` discarded `'forces'` from a `model_config` that a pre-loaded `MACEWrapper` shares with every calculator built from it, silently disabling FIRE relaxation in an `AlchemiFCalculator` depending on construction order. The combination now raises with the workaround in the message. +- `GrandCanonicalEnsemble.set_state` and `CanonicalEnsemble.set_state` restored the step count and exchange statistics from the incoming state. `ReplicaExchange` passes a full `get_state()` dict on every accepted swap, so the two ranks traded their swap tallies and the per-rank "Accepted Exchange (%)" column was meaningless. Only the configuration travels now. +- Per-interval acceptance ratios were never cleared when the outfile was disabled, silently degrading `interval_ratios()` into `total_ratios()` for the rest of the run. +- A missing `species_radii` entry raised a bare `KeyError` from inside the free-volume sampler; the error now names the cell and the missing species, and notes that species inserted during the run need radii too. + +### Added +- `Cell.is_point_exchangeable(point)`: the point counterpart of `get_atoms_specie_inside_cell`, used by the molecule moves to decide which molecules the reservoir may take back. It defaults to `is_point_inside`, so the box, spherical and dome cells are unchanged; `CustomCell` overrides it. `MoleculeDisplacementMove` tests the same predicate for its region guard, so candidacy and displacement always agree on the region. +- `BatchedReplicaExchange` warns at the end of a run when the whole-run swap tally is all-accept or all-reject, the two ways a ladder stops being a ladder. Checked on the whole-run tally rather than live, because replicas that have not differentiated yet legitimately accept every early swap. +- `--chunk-size` on `examples/re_gcmc_co_cupd_batched.py`: peak GPU memory follows the largest chunk rather than the replica count, which is what lets a correctly spaced ladder run at all (an acceptance-equalized ladder for that system needs ~29 rungs, i.e. ~13k atoms in the relax batch, several times the whole-batch ceiling). +- `docs/replica_exchange_ladder_spacing.rst`: how to detect a dead ladder (read second-half swap acceptance, never the cumulative column, which early free swaps inflate permanently), why uniform mu spacing cannot work across a coverage range (`dN/dmu` grows with coverage while `dmu` does not), and the `dmu = 1/sqrt(beta dN/dmu)` spacing rule with a reference implementation. Worked example on CO/Cu375Pd30 at 400 K: 29 rungs instead of 5 gives second-half acceptance of 18-65% (median 40%) with no dead pair, against three of four pairs at exactly zero accepted swaps on the uniform ladder, and GPU utilisation of 94-98% instead of 43%. +- CI runs `flake8 mcpy/`, the lint command the contributor docs already specified but the workflow never invoked. + +### Changed +- The outfile acceptance-ratio header abbreviates the molecule moves distinctly (`MolIns`, `MolDel`, `MolDis`) instead of collapsing all three to an indistinguishable `Mol`. Single-word move labels (`Ins`, `Del`, `Dis`, `Per`, `Sha`, `Bro`) are unchanged, so outfiles from atomic runs stay comparable across versions. + ## [1.3.0] - 2026-07-08 ### Added @@ -87,6 +108,8 @@ adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). Initial public release. +[1.4.0]: https://github.com/farrisric/mcpy/compare/v1.3.0...v1.4.0 +[1.3.0]: https://github.com/farrisric/mcpy/compare/v1.2.0...v1.3.0 [1.2.0]: https://github.com/farrisric/mcpy/compare/v1.1.0...v1.2.0 [1.1.0]: https://github.com/farrisric/mcpy/compare/v1.0.0...v1.1.0 [1.0.0]: https://github.com/farrisric/mcpy/releases/tag/v1.0.0 diff --git a/pyproject.toml b/pyproject.toml index dff563c..2d41b0f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "mcpy" -version = "1.3.0" +version = "1.4.0" description = "A python package to run atomistic Monte Carlo simulations" readme = "README.md" requires-python = ">=3.11" From f7fc5a80131ea941bea2a9b2a0f954d4acb93749 Mon Sep 17 00:00:00 2001 From: Riccardo Farris Date: Thu, 30 Jul 2026 09:48:54 +0200 Subject: [PATCH 5/5] Cite two conference abstracts in the research impact statement Add the CSI 2026 and AI4AM 2026 abstracts from the group's grand canonical modelling work as evidence of use in the research impact statement, and rebuild paper.pdf. --- paper/paper.bib | 21 +++++++++++++++++++++ paper/paper.md | 6 +++++- paper/paper.pdf | Bin 1053664 -> 1054415 bytes 3 files changed, 26 insertions(+), 1 deletion(-) diff --git a/paper/paper.bib b/paper/paper.bib index b044d56..07f52f3 100644 --- a/paper/paper.bib +++ b/paper/paper.bib @@ -111,3 +111,24 @@ @article{reuter2001 year = {2001}, doi = {10.1103/PhysRevB.65.035406} } + +@inproceedings{telari2026csi, + author = {Telari, Emanuele and Farris, Riccardo and Quinlivan Dom\'{i}nguez, Jon and Neyman, Konstantin M. and Bruix, Albert}, + title = {Modelling Nanoparticle Transformations in Reactive Environments}, + booktitle = {Book of Abstracts, CSI 2026 -- Cluster-Surface Interactions}, + address = {Paris, France}, + pages = {26}, + year = {2026}, + note = {Conference abstract, 1--4 June 2026}, + url = {https://csi2026.sciencesconf.org/data/pages/Book_Abstracts.pdf} +} + +@inproceedings{bruix2026ai4am, + author = {Bruix, Albert}, + title = {Accelerating the Structural and Chemical Characterization of Nanostructured Materials under Reaction Conditions with {ML}-guided Grand Canonical Global Optimization}, + booktitle = {Book of Abstracts, Artificial Intelligence for Advanced Materials ({AI4AM} 2026)}, + address = {Madrid, Spain}, + year = {2026}, + note = {Conference abstract, 19--21 May 2026}, + url = {https://phantomsfoundation.com/AI4AM/2026/Abstracts/AI4AM2026_Bruix_Albert_58.pdf} +} diff --git a/paper/paper.md b/paper/paper.md index 860bdcf..9a9fd64 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -180,7 +180,11 @@ and surfaces under reactive atmospheres at near-DFT accuracy. The replica-exchange GCMC workflow was built during doctoral research at the Universitat de Barcelona [@farris2025thesis] and underpins a forthcoming manuscript on the oxidation thermodynamics of silver nanoparticles, with further application to additional -catalytic metal/oxide systems. \autoref{fig:phase} shows +catalytic metal/oxide systems. Simulations run with `mcpy` have contributed to +the group's grand canonical modelling of nanostructured catalysts under reaction +conditions, presented at recent international conferences on cluster science and +on machine learning for materials [@telari2026csi; @bruix2026ai4am]. +\autoref{fig:phase} shows a representative result from this workflow: the hydrogenation phase diagram of a 201-atom PdAg nanoparticle, in which `mcpy` recovers the progressive hydrogen loading of the particle as the chemical potential increases. By coupling diff --git a/paper/paper.pdf b/paper/paper.pdf index dd93b40944cb2e083bd201ae77327de232f3e309..14a254b73a0db6d7ef2a992701ff34d41b9167c9 100644 GIT binary patch delta 26302 zcmb@tbyOTp_bnXU8Qe*RKyaIZ!QCA~aDoL3?kQ8iK<& zdGfyZ_uTK^KQC+5T6IpHbKQZoR0Jez0H zLKvK1`_-?7GplJLuZHHGOIkGd+KuYnf`_ORbm`?(i;e9jVw|C+8 zDOGH&j~Yxnsb9u&Bv{Poa{?bveF;qgWsIue7)!Z+a&&o_DHAB)g?e5U;yzd+~ptld_Lw&Uy$2C7;afG)6 zKCk3(oHs2T2=Jv!)MBSH?l5mNomdXf#QCs~mv{0XZlrTLA0=0NR|iSMc#nSVc9U>V zvqbNI$W&pZ3r^pY`A}sXdNRhRMClFXVJtOKU}D>Ve%yA7X>L-GCO$Qv_PFcWvi>-U zCmt+UQl(LA5c%wt{Qg<9HWrbX!%ll&`xF!}^zf4HLC<9gSya5e{*=<|naQaUPZ&Su z9(cC(069A_D~FQt8&*NZVeKsNeSJxT!=OxSw3JixQ75C{ne$O0JSi{^KBh`NA zY=*P(RfjUE*6QfX)r8Twf|U99dmPbEbkdez`V{l>W>ufAT@@$GHg2)&$M0IPioBBV zU{!q=pt$b-rqOZBg=2Q~FydmhhOVtcY^kE23>rM}EI=v%9|;fHyiHb_WFlJdSzBJv zTXDjZme57>hB-OrW-gu>^IZS*bC{(Is#0Tl@ADX7czBqbD?_z5-%uJ?hZ2uhg-%Y->FoZv@#L1DT@8z z`91LIj4&3nfjU(=<;SHI|EJlu@5P_dd&6aRD%{iRe7_#_@qTJ%vb;R4R5jcICVwKz zA(>zHNiWE3$>AY`H6O(!R6B|RFT znNMY;buT1`-O7@>TIL2|I>b3qo)G$W3$NT^9DhyIX&^Sq?_*px->{a_-n!R#t7D!w z&oWA;bP8}kCpyfLVJjlTGD#B2D-$7_Dr)d*V)%=mPp+-j9HtGwE!T2n$_`5X`aCC> zqRxTupKhzIoyA#L(OiC2=JC3=EtfNdg3tWN#>>nW5}Uu%4?0;zqgE!xK@{t8olJL~ z>9cAUNpT8C-gdZo3F^=K&`d47G74~G`YBX_B*u}ML|IX&{Siq+RXg}h%2ku5;06Ch zqiL)P!^e|J3WrY&{fyE(7U=7#jNIhrPlrFV+6F5yGRsTv@lk51m6MVjo$=IULj5?S zX{>G7bTgjwD^hFCPAJhl9b=}PSMEkpK_1PN<@D%mqBhzxc+>H2*{CeZnhd`U^t1Fy z*~_AVq0BJPMw@TyCS+qt)(okeZ1ZC>T@)h*9HMWxHX9-A*Xm0}p}Sw244c`zs;-Aw zS-uDz2mvpi9DmUKC1!}0jsh8FbA?KLz8XSK`EX{Y^T4dGzVvpmBklf|zW7)mHm;rT zlBHjOCoXq3XBO>Ky>Le=>BmeAtPnl_9v(=Ug11huPtc4Q`QRME76JN9!pVH+;BXFSh>q=7}m32iAMM;XGKHj&{S ztQ7ghgpM)DoB1Bjq-?Vfx17*$qT1}XcGRQI<|l7NQ51Si&TaXot6xWyAHln?I%&_O)NLn-X-F{`9M>>AC@nJ%I&GM7$ zxoNPxdX4SrJ6swBy}of;*1L}k zkYSvMrU7p*pf5J{j&G#qVXjf$F1@kQcvRMWW0@vk8E@KmOOcxG2^97^3yr%gzQlhd zu=L>_X}LJ2rW4JV5*O3^lY5;Ip%i|;v}qL}M=25~aw&?(p{3_7RhLA+gzT$Vp^nB? z3N$W_g!X1w$jNhpv{qk%xry(1>7DgSwK!`?(7Do%?;Ni`_rC5qj~w&Li==}+o&Lh9 zjYZseuzV<6E<02<9sqT)c>jilj!&*YRb89&*ND$X=ygf$??BO3_^WpF5GL||RXRo^XjVzj^05TBGtGR6_>tDRSU z;wt;{%QmB*yo8o>^IBtzLrn=OXKVvtgNizCtUEUy$v3v#-8aG*KD4g#u6SOUF3-{+ z-kQ`6(l7uE&*O@YqSGKZ*^inq^4tm8ZTxwe4zp3C2#+4pk&Hgj=O^i?eejqIYP&@A77woQ><$=gpY#o%m{3vV&HmcE*|82N)&mFi&!{nbZ#VF`3kSEry5B_OJ#TBMW3l; zVB>y$xvE4c5(Qk_Cp9=94PYxyY^R7<1wIv0A$^|+#!Kl73FtGcjLMq`^Z8t3G+#as zwU|xXg=XY@nsC6-?!m66w?caVDRZShr)(lbZ3lUK&yp5|!%}?7y)E{ZH2E;(TM%Ox zpgLsU6~)@1?r|W&FLpuCEaf&!rf!~-joRGFM0wPvy<<*^&pu<`wvFX_{7Gm@(*;!c zB2$azPAigXYnElM#>>t^O!W4(%AsbN2pe{zC3HF00j=h9O1plVnzpKCfU@OzS_6Ca zQ@)VmvLKG;c%Q|26u)fx53QW7{9V`)kq!dRPZREb%qZdG$<(WcbC_oz(^T=>%wo8w z^I;#Z;&=jS-Rm>%Rr{Q)S5sN@ZRG+d2~2sP?Je`_u}EQxwg^qXwAAM^Kdx?cAl7b@ zP=(&Qaa|HNj~U!HtcNcCI6~^d=(`oa$q{LA^JmSz8`&*r?MT0=fKmPaUgpa&457jE zi7+DPX(nB)m?}M{T*FHG%^ny3T*^srC%pEmZ1!=>>b?#hyZC!T8fQ8|ZLJB}a(?)E z3z>1H2d*doEa*wOAkPK<90s=((UbDIWpk(kO}tejqusA)6Tv2&HjRLW=7fOBa*8_C zTCa$HD(>H91Yg99(*Zhp~oW-jrQw5rSS*na%d#S#R) zBjzsT`sbQBP=_7 zn7YG2n?|I=xd9&%{O&NtvmNQ?;PK(J)iKdR`AVU0%n*sCR1Uwls-U}J zr3KedIibl;$Xvn0%oj!)M}d&t^0OnUm~+mA=u{=|Up&O3cYDDWli2xV*~G-(iYx$O z?_9!9B+?^W#ZWT{kKk~!6jO(s6r{bfl#y^_H)yiHPOj92BGDI3qE9uY?#J;^D$5dg zvvDrL$F?t8m$uu*=&Wx#Ce>bbf2+I^M4zXBRj2+b24MSZ9t-8XI{v4B8#9d+&Oy1- z{=B|8CfgGaA#c$cI$6F%)^_ut&EHMrI$|ycpH7slDs(b;3JY@X+}m{5=Fj;I-_HG9 zGcMcS+gN4Z%B<+u@5v3lZ3(`VgudY(x2S*q(tCNVgr2An@*+;(ty7Bk))XV=ht@C! zRWin82gt?pqYwttz?_ZGI2Dg_(DlFRzJuO-~B3Bj{S_gWR4&GEUD?1tSRT@EU2 zjT>s++3~zAsWrNn&mijP5gGT_uVfA^vimHUe`?NvmQUS_JX|l`xBcPzOHkhnN^1ie zv^qy^ph|r8tBV%$*%o)|?58Q2icFTT6f%uFQVb#Pp=b*S3i0~+S!!u}A!XgkMRyp6 zOfz=)iCcDMt}zOsf#!tZUuA~%8m2w#LjtF<-#&<+uOka93zd2pc=jJtltj>nZ;v~d zE6K~7U0A=9xh|}Eh(I+5J*;1E;<$GHWkN58Os37-{g7)VTL*`lbufUQCQ2yRbp2*KU+-z}h@J zTG>T`Tv=XFN9MMETbyj6ew}W`55A2Y2}izT50!6YUo(tGnFuIT<$U=OZfWRbZ`0#j z5vGh(8SZ!E2${JoKHM3bxLH*pkM>R!gZwUAC^51FCw0;_@75=`96z28T~EZ`Q00Or zgI%t#15Bh&s=Td*>q($hm|;%-iuxUupQ|16tn^6Y>gN}cjK%J_U1Ly9eM-+ez3k0x z%Rtz!VopML58Eqsy_*-%uX}qbB&J&j(Cg7>OfgABa+_DT8)d!F+mZ{*vd@!mkqk^) zo9QiLHp_qF!5_N6n_stPwVUwdBp<98&&G|6t8B*z>FJssU^ze|yOv*mF6jv<{AK@C zNH_=N!`kkUC!-Tt*aEZ=5Z>O8&onl0t+18~>IHD&lR0sU;F+lg5bLiot_GCp!K2O} zVmod%b9;7-O-!U$U-$NuUi3|WUb}xUHh8^***ck3V08ed^FSNc?-DCLcZO@1q=VCb zq|Rd;Mif_$)hj|pbJpA(pLWqV zF!Fgs_fZOiQyzBk3^NNp9F3u&3?nhe1i=;NF$xv`a0-GMM`PH)LxGryQn0tCIN(TG zb4*0Tg3Fp?k9EOTN^r8_l|Y1CIeV&W&dP|5%R)e~twHb7;r}FPOm_rhzT60PtX4W&m0gK|wxPg8>o=^waB6 zM~3dCp37IWI8nPqt4|~)l(fM^aA>JhggT!U=YyFy=lRVpktBMivUbTVN)esACF3RL z(3x!F+4Gs+$NT)V{I0vbXV8&3YVE-C7*{03r8oco=KC3?Gu@MWkOZP()+^fXf4@`tP>fn>P%_4P(4%FXrd&SagYEn5lb0Gj*Y zf8((6aCWXvP271elC|PQ4f##JCmT$-93ImmDM0J&@h0|33FfW77|Y>UXk4GL?|tOs zm61{Hqe8%==b|nF(dy2iAIq;cwU;tww5gUvn(zZb!f8XN=995$CCS8TXlYtVv*t-= z;0m0a*!daf~n5b$%K0P1xep7H8bJ?II$*^q~Q*px7y z*cqc&7DFE4OW!(p`CbUQSj7CH^}3=>ay;U;m2y2`$@{B>y=hohOCg6%leV2)?R;N^ z*#P^TI{~{rs~m++ly?fh?GwZ`Rd$_ic0y_#n4q`X>N;wo=uX0iv4#DODAHe7iBz(Y z;89~dt!kzPnZSg3n- zTU%A1&xe`|4`+d_Wb>@glT2{kr6d77gTHvIG`Wr3(cqr<+&tE2)bCwWEO@#Z6djQI(gNI~II3r70dhkyvme9IlEV{Z3Qw{?cMV`Yj}ltVhL z`*M1A`P$;}GhT1|vZc4-=;|C)?^42UdQ0{I&BM_%T-bA!c9(TU~_bv$U?$i=T0 zk-r<3COB0=lL8BDB3(URwlM`O35L=6vaOiIZ0dsP`YqgjsQH=rH{N@!ygs?>J*;iH z)EElL@{l9qn}4THVI9X2LL_MgWWYXLcC3^F@Zaneo3=|+ym0*GD4#^r);d}AO4SLH zxvmobWLx6vr)T3yIz&RMyL>otWbf_tE4#Sso4A2RQ0Mv2IOc7s90A|GtUgk5!YgD_ zW?r7sei(}4M6EI>StiDr(GIs;8M#!jb&wlg6*)WSTc6d~hVKMzn7(?+WThNu#O}R0 zJMy8&hy2!8%(t*Qs|r?=o*mUTm{s{6J~uopK9?4+@H4uYScX>I&*Cb#RZ_ux*BHlz zooAWY1*-KXeAWK}<&6W0b-3)87&gmgy_@3>%&WCcej*^w8?Wdwv1;ZcV(-wWM^EZ* z_>78Yp5i1_K4;r5j;{8X`C!gMbh2)bd1Wy4MXqyBeBQ2Lxjp^in}dyB3q2(+Ka0G* zQ#r9>d-LzJ#c>g_D$F2DC@-6Kzs9PDVHc@Kks=hyI=pgSJtj2RHctN+@l%tGmV~pu zL-dUABPwR`tkmMv6P^zSd6>=yBicDC_VEP^2K>8hPoyoQZ6;i(k3M^Ow=_$pKBeKW z+~@&Ii)`6HN8x9?JIGv)t?R!Np^kK@0Yvx-ILlY>G*R7tpu)ecEKfO%nqR>E{ca@A zWugdbE0b;|v?uX~#wOegFK|REzJ}WKdK6Y>WD}J#bL^(jr;3BL(^$a2cgD&0BK`Z_}+?Ijm?Xca$6|#8k&ehtU(HWo403|6677m7b%| z7)F=yo{gTv=ikEGe!I~E2c9Q7VFm53T%ExKt_;#~x6axFH6&Tb*mz2Mnne;!NEj2+ zIh&`P#cS9vumBZ9}@Su??Ti^f%=@VlxYd zzAKN*Ud=2ZbfB@8nXPrmkO^CIJxbpeA`BA4geQo3DFPS*grJzUysCVCSpvoc|Ks?19oD>KHOfN zXw`PkLK+^TMv17ibxrw{r6uXHE7vs0O7u0C3$ufDepNb3U;9e+Y}IzfPUuZ9l$f;Q zhiKr2gPSdvQyj-2ZZaAlUZ{Y>z0SPQ`|DEY+Knjr2?pTVkkb%-YJMD9O*pRA%-d^; zo1)ZON$trv*qMd@<&V=_Rwuiwn}^@;jrWACN>(a|qaKcxA0-pns_d@4qdlos@&w8> z=F|lVwkgDQ9GFwzuI5}*gwD1!L;IgDMVYQGXDISn@uH&Y>*M%3O5V)8w_xzuZrNeT z$a_(qEOmfSi9z;ebqptXg!=_bP>yasUno`<8tsDfS%cf$_fDm*|d{sKVcQEb;Ec)AC1#eL%{fd2p^mx+Cfm zBVzcNN4{EQCS@UBma`h(veO1Vhejf)pIoUZZn@mA0fTaX*$@Rtp{QRtW}-xBfFiia z@O(xbu{Za#*Y1ulsv=mj_n0$Y?y&LiqHgiLT%xWL9LZfV{@%B7g(=-wmpzwi(0>N z%7fw*I};+G{4|I3^MTLOvjyGL+h5x82BE*t^*+njltq*R&)=1| zG;j)zgbMLqBF_u`Ol-G;z6v+57E=4JE?W3tP%ZRHv}R~=qc?Zg?UiAcY|C~qs$f4V z{ii65A6)qMOoMbP_M=A2bty+wi^?_%>4yZ&-W!rcQ{BJPif4(Ks^!lHB;$n8Ik9fD zeb2(^-qhpiFHnUo22Z`FItH38;`A$cYyL{VP#>Y4FV=HZU;Y`d1Whm0`(YS=qjI<+ zFcE-+t*X}a_hQB0G%TDmACKI1wvjW=7Q zN0=M~RT=E7&(lQK5Rc3Pd4)F?%8rYlp8~)fl@a+SxDgU$@?DePxbx zih4Flec+xEKm`pn&v}If54Yu3v5|1p;wJ>ycV!((v# ztGqkWO1nidxZbsJ2eg<>9wh}qGB{@5HNq1EV#+tS-gU>$k*iu(`-tP5+!u538Ie- zJ+)8<+63hvDYi>@F66bvK-P6F@i$f41=nfbunR4i8N8$ISvJ!;_8(Y%Zfe(xWm9YY z_9i-?)Sf<*+V8;nEYg(Bu^=4IVepzAo3{02*Ex*eF7JggQ`cE9hqOK0ya-)u1cBi+ z3F)i(n)@24Kt{a0OFfnK8gAP`n}=9ZZt#Y4ex&rg84Aq>##d(i-SOdZreuhiH@2|x zX>iiPW1fpEM*ZMjyQQ};S-j=awwK#{>kNi<3Ou*O)@S05ecwRJC5ZgmC(<>769<2zBp=% z&>=18t5R&{Bgk4SSS39d@g_0k5T*T7_ph6d1kY+Ir~KPX_m05rGeYV3`teanb6s-k zi4@vMq>C{1H#hXLe%V1_?C>m6$^9^A&Jw1LCCSID%tU`>GRFL)(E~Rz4$)zJ>hX?O zm)-@a+H?@dVd+732hZGe2aFaAlT6tVa7!6|@>Xj|g$U%XJ(F9Xrp&x|u3;diSt)));A;OBy#mxX}g8uu85_8hW=- z$IK5b(hq)UQGH?BbsM$lypb`JNAW(EHeOf90)oXvJH!b;Ks44Is+(HkEX|~{J?fcI z)qXp-$GInMN<(o!U@qM59i^f5&XQR9th=|!FDl~U0(yM~*^YWAa_wi@JRa!2;a;%g zos8PwtDUbYtc+K_NK*aAp3Cp6|?W>tyc=2Ltqo69#RIVYL%Qs{A<#Fz*;e7xySSG4>qdo-6v(9iYqO{ zsLOkp~EA!MD#&JY8YKgPYha=tJ9-L7lvBD1RHn)XrDk#R$h(w=aM{?I&64o z$o5Icl8c6L>h&Z7)0HqG#$neAz9*sA8p{%jw7I^N0s7Dd*hMJm_A~JFyZA_`;F|HY z()y<|FkRA6hYu&<+hbP7)0|8@J_oAf-+0qlD4uitj94BqeXpL+O>>siU}^l~%o|7I z+#}hJHtLyd6X<5*o9FfPzwe?;BA!iE^DLm@rx6BeSX~)gi)}{_fLsQ0}vwciJJ^)C=kF6Z2T2x!JIWm`UfeFs>o$k zQ0-&QdlFxz`P0m(r>><9JMYg17wyMqhn52@ISvX^`mIVOf4&dsdPDGe7+UdM_{MwI zj>-`H?S9oR=QBwfas62uokcX3)q|tYgM0ZYo#+SN2$REKYc3Ux;B-Ek2$*`)i41V- ztGr{@V%POUKCPIlV^cFZYzjGzE9LRrWcs$Zq*vU?VwCTpSM%&LMq7>n#&rL1yVDzZ z7XklV=5lL6@iwfQNPe-w5^6X8!cf#BEwRM2xu>l5%+K=Bk~eyda!IncpPIJdg~NP| z5Zg_w=@YyN=7Tl*GuocnY>cD1{+iKlXWE#RMtV$+*KiNqU^98K+5o0vvR=FJp+r_WZyWvUA2aHOQVFMA!S}!CakEpP;K~itZH@*1v^# zizCO+Ss6ctx<;Yw)|Ray>9ld1Vn{Q6I2M5P9fY4|;p86>8?~hPh5_?KX^N zqTUbW`Uptl<xCNEa4eLJ%5!19WH`PHR3e&?= z(2BGTerBi{mnR%NUC1N>+=@?So*I&BlH^BDTzadIX@Nf~LWh;^69s*@pBHICNJ2Ts zs@EqSivIf7fB}aW@_Hvw@561-(bl{iZ%BA<7%d&PPEeeNt-E?(pU_sjRD!7{i>Z%% zd?<(MZ5#2_aN>G_4X%Z8n{;0up)DyVc$2}5=mVzc!=(kP_ju&KJ_E3xTLg_^>ZgJp zZ{lUi6FwXgB69{kZ+&Q?yUm9W*!a82*bXQMQ(sv2({$Q`cV0xuFpLl#MYsJ>xMH&X zS^y%bDItua-m5^l(*9CR+#AiErKLq@L&5qh&o?p{(|wjYZ9=k<1P^Szt=4f`tWiRT z$E7U)6>Z|hrpO|)ak4X$Bg5JzIfY4^ukRtvT2civP+1upO)(rAT1X}8{g6Y7QGHc< zU12t;{99XB^|nLb`8qb>{^GN4ukG>4!Xs;#K*6}6k=n4bA}pN~Dz?(DlhGlkCmyZf zHOI`i#xX?Ri%nl{zwhdPI{L;wm$QjlVxQ}YOLyPw?l-0W+0v)aCY?~sYmu5XGf3_m z9Je89O{a?dtm0}=zIMsF?#j!b9XYl}JMKS@HJwC#5A5xQy0)S^hLg^H@0!|%H zI~JwU3U5XAp{C14bz{s&XrFBI%tJKDv@>S(KNf|k=-u4aM}FDr^gT6fm{bAzJcv}V zFpicCN}eXb>kf`}nF65m_ub4{o(7ZNGf&A+G8P9-TBr%WAfyyWl=dx0MV&=+^CzZq z&C7++SOGC_w;em)lCC*6pRO^!RgHaZ`kwZSK`<*#J2YU$z}scPV2@w2s`)FB+f*=J zni-^e^g%kuljhFB&-dNNcS9NTiO!hgyd<~89{W6>^mCkFtKfkv%owl_Z+X)lw#a!mQ_(G_jXP##W2;Z@=c9WzHtHh z3@7>1bIGumg^3j$J;l=eCO5N;jBbP~-bfQnGl4fnDD%xqp=^gFBXt?f$idFvd{k(D z1BAx%u71H*8S>^SJxZ_z9ae_NeiRL=Nsa6*-jSr^NC@$3n_|X&bf&Jgb(x(U4c|4; z^Lj^_(fYoR;JE=#1T*rR@WPC}Sg3@iGYHS;vJTo$#)r4_M2Ga_Ir-lCH7kGnvJk)@ zu!tY_+p(sG`1OWO?i zNm*l#K>zlH291F46@;1{-G*mfrkr~|k`AeR;W$2dlYIf7JN=JhCVvJv?zs_>pLHfe zH^OtK(R_@0JSUyG2cGXS*A;~8joZk6@*zlYUimSdzNY`K&T-|V#K!`BS}#UmK%$EB zD96}~G2wz~!-!Oh$wX0X8=e+fqV6znp5mpuR5|jk_SM#ks70WpyVUbEqz$B^L5xEN z!seCT3?Cm~-+9pw>_zXDuJZ?-HLvn`p^TQ2?{J9NWoCynV(~&c$-$Jdt(@s1N}D)Q zaH24bthnJTY>gy^alK)7+CKVF@WAfuNHGa>2tmS-VWN+W!I2k*B=GkbgCiPiU)q<) zTU5GzE;`J^8jO$v2Il-ObGx`v4%6s|-i$Mw{9n^K>QQ`|ziil?axTn|@Jc2sps7B5 z30CD}g%e~^`fd}*TPM_X5!3cCQ68h(Q_~>cR@DntK9sDuwQMXWRA0DU=@36 z%X{MP-Y(WNaMVn{j2bx@2_76Ef))j`>j^30vI?amJ&%0pl!!qsJ5CWX6L@PPyKhlh z?I15BGy?kY9wgg&%HJncnJ~E$tv=p_csk;v%$5Ei)mm9BEpb9|95V9`t%RFH{z7kc z;cW^!(R<{^s^$k06^pMd2J1un@nfwuSgAkND{GG)3`enZ@vo6vT`R7WgP>of@7Rs5 z?cEN<=1mLTx9j?zV2vQjL`{ zC0sMJyBGbgzt2_ms(!vH4TZ)JHfPSonCslV=3te-v(b~4R-D;q)#)+7RAQ6 zmbxkxTLhHUTG9r}3>y3VcF=y7kITi>`-}-V(uZ`Sdd;#>T+~|pdk?z07^F{j9AU!h zUMH|$H3OYKh?xB$`;?)R>}T?G78g3wW4F+IstQtvhHCT(rP#~5J97}H!*a;-@|Y58 zHuePCDBc;0s$T@fFusd+=SnH8|996zn8eZ8a%$NMBUW<fxMrtMF50rC96QS z(|4R#Qr=mlQ_8AYiU@k(BtVlt#?1^Hp<~MP)Mx8J|tNmmHX9;>P#?j~;{PcAp zeuP@@P=j2V^B`>$DlJ`lQzrcVQ~YIXwra}OA&hmOzNt#amSoU)1{kW{a*7vQ zKQktnyW!4!rsMGV5j~JkczK8d%=q4^%0pg=F3IyIPj2G@1$Bs;*mn^8fCG88YE>{r_ArCVbHU?t;%T z;ZqS|+QSCeU?K~G{t*ozgZNW`IAqWoy-`PceuC~>wi}Fr zP0iIO(h@707?~>!9eBH`!1e*mWW<+7vW?N%49aAV0)@b}7p$a*7J^gs{A9cW0)#>S#Ecpncj`+|h5+MzEiqH{c8{k8WXy6Ps{-KFq({nF1_oo8? zC3rbi&jM^ZZi_li`Sv-#A!;){4a_o=VI6M(K z+6*AONyCC}8;9N)=+5`ab>~kZ+7J1 z#fDy(4caY-ny}gwCh(vg9Maf~B(Dg3nG&U?#1g3tR;aZ^3u02!Ir3WaHWg-~5$b>4 zzv>VcmsgAFwt%IUL|2>_Ds4;HY@MYq(P({^q>Y9bJY}j#2rQxb?wE39!)PAvRYoed zBJDTy%qOP8g5+C^N(n{9G_Qq}k-&UFo*)!M?nbVNZB#LpM?Mc8w#Bi+rjJpphtW(b z)UhN)$uZiyMDuRYjuaB$L&k!Hp<>2Qfi46hup)4i_zL?dU#H>_w$$%USyJci6Zus; zraWjza)a~QQFSxEB;FH?3(#AqqJ5|xbTk~$BT%1fdrA% zDnf5{Q-|W_PjuHdXn{l}m6YP`ALK=$+TME>WCc3j$kcNh=+vP_bY1gxcb2na-@fZE z=(X;%=_ib9Gwb|h9?YB8co|nVp=ZwLo@By;ODkR2=hML!-_m1mw7U^)v8E=xj|MAd z(GFtO-RDSY=G*KLi*KY{NSaM>*5Jz((VrEA_=SLk|eua){KGiy6le%!_v7QMk>Z7q8frK}^D|C85x_h(i9@eP1t!ktE- zmGkl?L!9xXf>N5CXcWUq43Dvb6LdF4c4G%=#nNGr>5$p%Xssmo?YG57T7QO1!tK1- z8griX!X-^!0MXB*4+cTMuWt^9Ng&T4J=5kLBH|;h#gwIO~FYyZ&0kU4t6I(9hx*rL!od0pag=7 zoce-3q^7-8oV(4F;?s$6YN)iKzTx@i{#?8fNq#ucsaZ=^`PuTWVSr_lVRZFDY^7t& zJr%iXcrEyHj$^NgVLRz^9$ka~d{v-M3!qmL7Gw20D565uxtrtJ(C3ljH^Zv|wHYo& zu^P26sATazFNfpR!#@@4;z7%6q+t&ql7hpgkiw+t7|odJP3Dk<={J9NvF)X4Y#*?` zn$X$jZ~KtS`)Npfmh>&7?{Q&xd-1qznKdWBIy(*O;m!>Dc=~Q-PL|`sw3yQ4l2>%d z+wA1Jyx#!dz&Gb&*O4F`uO@-KX;nr;UTk;Tr&wWptVS|{kHOo$+viYH+Je)?YZK$^ zh$I`-#pdI>orN+dck_3BzJxXN5%f#(sT>g*HS(-KFWXf=*k?h9V+4<8T(#f>(l*7f z)hn{e^&kzFPA>s4m0!tL*Lm-IX|m|?Z9^u=?OW)msUlc60%k3KqTY@N3!%)Q?GieCdhJmY+IDZ_TWi;8ubgMP^jqpmSSQP8 zzlgw%D|*MNa5lA={PEEsML1FftV5K>Tmv6H`fxeF9@d;$S`gN}$HzWo5qjcgcycdj zVXDamRpHMLu)1EHuPrh>pHxz#lc_Eud)YVOJF?kr>77E3l)Yi)lCrU^$ikA7Ep;vC zRg%UXxcd~R!M)dVzjCt8^)at8|MhKH{?Zm9`(+_&B8iOU1x9jDbOwG)R@n7KPJY&! zvU*9si;2Il;8d&9e)x=a>&kjvK3UgEpc+iqNdOf22#CO}RsMUowt9)UTeC3%2uKD0 zqd(NNbaV1>HMevFf&>I%8cYB&sNi4mqL|p9GdZU>w%DM5Jb!wIK#<^Hogzr^uRIYX z_*aAo68x(@1PT1t{;QG%3H}vFf&~6M{?~8xKgs_j{FD5TCR5AD8Kw>fa1e3>`S}Dn z1^z(!+kq8A0Cuo$EyM)sSOKiW+(0l0#0mKa4j&IUCs>gGZ^*Ri-w-xf4=(_NoS%=M z6C%J1gh0TY{186GX`&#fATQrPu}FfjU|xXeQvndC0Ji{;p9jpzij|G)@5FQ{O7{tlP z^UoAzVO4wp2815;A3|{X0sbU{f7G!5-8M47KRc5fmctKV;{-weFGGLTwtD&o2$d`d zksh(yi1p#;;fGZV0T{Ri|D(_ScK{zhm{X8P5a{S(Zx1y3_reGWrV9h~V3xuFGT1v| zKnV(%n;Qla0q`I?V3>;)Kpb8w0%$}*m>`&o1%MRhECFbMp$7thfApy|a3e_o3<;Ts z`_CW{cCQFvf{984m=N*2u+3foF$wp7CCmT)nfw_+1PlC$BZaL<0|5Alw$GV9p&`uwA%F~)@)D3oiir3J_Wx{jFo@@`jMqUE07pX1;lD_jVD?%77DSSNSRo<) zKS)6g{6h*Nj2!mC;twgH{~-l|=`S`2EK(aVj>yCRcX3j60HW9k^WUUAngViR5{m#J z9I--F$VmTW{?FPxwfG0cKctYL{pEyk=`Sa`h?4(92jbt4KPB10fA-VgE2bb~Z}K26 zWPh3wLK}`S!TbV19)v0#Ve$hFVWvv}I+%kc02Ah83CO|VLqzZ*F2)pA098gX1jLK5 zA&A!@MnE1AFCrfJ=jvkw{9jj4F!(R!sE|Lq@1OT00eVmg})gn()LPkOZ)k2P~T)@tu8(!nJ&WmQLwdN=l>!rMtWSfOJYor*z4) zD*yZTxt;GGt7xR(Ax^+52bs^EDI2|U zO1dpl>+yD+>(qB{WzOpF>pniOdDl|&#*Xk#K>AqqN!+FY<`WJzzL$pLOYBi9_wr}` zZ~>!hROKs$CsG!Z=hid zNL9fIA_gK$2PPMOXLAyKKl~?pm(|jnu7$h=<^lXdR)C0%jJxV?yOFt2D-+ z<@VTd^7~9hs{OA0#X`QvO6p1*syBsbEHMp3^IguXceY9s8R@qaw`27I^^#R|6fauUzp9X> zTX2@!E!pOUv`7`GOy2sOU=lski>=AFZy)MlKC9m;+bI|Me3h-bmfK?=DWaX%+UGg; z9@u)x>2)j{=}J^<3D>*kkw9~<8yF@dx@i+miAIM=Rk z{6n{5bta2WxVn234Z9r#MHcYMMhrve#fG3^j1zq=!L26mHTGV(wL({(m@+M){S1pa zTTN~1)ur(ml>WRvG4D6^E$t{gy$MIJzt{84B%aZtIU;( zR5YFQ*xF?9A40(L^VudUep{{H;o4O*P@pCFGm2v_{ZlmaZMPKRfE$O#vOA^DsEQku zRvNmdbcEWYS~8N{ijM4!lrs zahu%SoYF9~tD9i4R~}Iu3MRJnt#QzS>uFeb@fb(N!QuJi@pVw|&^WxLp4je4PGDK6!i5eT~8|rzQVoVbr zQW)kkuNh@3>O7#!&QHED?;QcfY9+P97U-nY2$5-UNuYx4EPLO=JH3nmDOJ~4)XRY` zYyKj&vJM1go3(sZe6RF|>*5kqW6?=C-+Z;pDO5`&dYqh-wz0ry834GMYeC6Vp>Fq* znb}vROb@!?gfF!CCeY}^<{--A9 z94@+|yeb`|-mPLIf8!~?+9TELYa`tq26zXKg!F{8(|oqVIFNaQ>M|UcesGaPWv91i zxouh0kU0y}i!JnfUKmr1!*XHe+GL5^hoVQffcTPtx7(?Ya?_xl^7}1%c8MT|p>-%P za}{84mhN5z8lVaMVIaCudA0{;^h!DJb9&+GZ~J)~Vqe>jx4AeT1zNK#k%Fz;J1ZyN zB;RNZlA=XsmwaxnKCyPSKi|5F&of~TO4}L;Q|(u;rM2OFwDt(h`s(VO(jaV8{PC-- zA{@yK9*drZ(1peEGjYKwGQI#N zTURTAGqMU(DXIX83nPo{uLP>pB)It#b!!V33*HCxMe?WepRjD)joS(HZwyZt>nHD8 zpN-wF2q?XM!pUF+ZJWcYGpi%=)D*0ew~xJb8{VkBHM|-QqM5HrK(WmG+{)^Og}%gK zPH*N%4lSk|=!P{}@0|>_VUJWt6WaHK<^d`Z|<#~v5A6y0mq^DLsdlaGHcuGID{ zk_G1}q2}B=IRYIg2CtJek;}FX_WKCCI%FbKLX}yIaS3rO92!bEh&WZ`yK%POdbg{e zExjPZc>kRdZBEKYkLK*PP}>3BnMUm?5%k*@$=zA@?FR79iaN2fdr}N@Hc26X9^$e@ zGU))(;~!;9Q0HUktybRh!ER=%Dxe5daG6kwMTSnY2Y)KFNAYf_>T&w)sE zid`#KXI(F2u+A^UHPZRwejsbMcgeKaaJ<8$^Bvu1%9;l6TKjvNbc%TUk+*EvS6LWhTw24p=>8u~Fn8!kE?mcW%lvwQz6P3EEM6*&mt^C{L zSMdv>u5`k}g#`t*E`DeOXW6Xma%+6VMgV{ZnAPCqi_+aC z9ENKL6r}SJ``@)BUXW(3n11-^@%`OT!EYM)l0lOGbl}0frlMt%n0*&mp-?G+Gr-Fo z!&lsca=t3_l*gQnf{2^}8<_aEE+nlT-QnqZxj1Q~v3H0vFPn#*! zXT*YU1vMXN$I`IlMOqA<2fqz}{u*0`_t0xY0jX~BT@-UyDjWg`WbxGk)A98cgBIGf4(X~_wd(o%@vnEAL;=t$Y&?S^y`Nn~;oY?+sN)qoF-ow$a+GUO8-|K*QgiQr_ur0l?GqwOQKGn z@uyl2RN4IFbU=K;$Jp#yipHTR9hYB~%`@K7e)x#rubK7}o|SjekH9-KLRop9Spwy7 zs!wrP8lhBQUeiGHM~Ijn$KmA&Hb`iY&TJT1V&_w1iL69iJDw6VYp=J7~hiMDS;_dYR?+QK-`E!o+5JWLN^s6d>=_2qOT z7bV?Pj9;q6-UW#@DGv;L%?aJ1mPx$O3glxUEBG|XR>=kxqlyNk=xWb-V050D29BU4 zh>waNsBUQ+FESngCNTM&?I*HNy{f`1tzyypD|Xe1e_l{SJy_68VxDzFQ^%IkL-gqO zb!=fTLhO)r`DL1D8L5@9@Ess0Lr?IQBW_n-zG&{NRiz^>!+M|X%;p^T5B1AV|)!!rB5wN9qAbzMR(p!2vo*m znXz#x{tmS|)~-$hiEYcaf_$+JMAy1Hb8}%sW+7k3oyBpn&Wq7hKAt^3W^?MTk1L2c zVz*s2W8tAoblVEx3R1fta9FXHem`vwjHAgNVo8P$%`c$kRBTZ5(nJaS-2oDAkoHqI z(#g@tj1G;L4^8PlWOSa5`8uhk_}M3=eLeD<>>{iOJ9snqaQozZ_^Ap!jvVp?nGCP+ z{=D_=wSB|meOi-$>pLsl4*eoK!=S)eMQ;&DA|%?fM~;6R&n=Vgw`~Ou&l7WLM^W;Z zEKGoILcE}u_-U6-+``Old_zn+LAlZp1gG$DSPYVupL?{Sw9)%R9&IX@b{65n`kbblK!WVvL0o-vXxQ^ez#RH2$HO5XWS zpORGHmjZ*_tOt~nUd=^PBrOfAJ|bF>&Oe8~aeGp5t>d(khanclCHZ?EZ@992E=NnE_`=28ECY*TO~(rRBjXnfmTs!n72!4$kup;+ zXv!_$u2@6zxyCtXx%OB6&*E0{4conkL>mMuo`I1hdtt7c2t@7LcSC7=nZnzYiEA7N zip^`@1|435Ea-3jkgl{WYQWGN3yYqcuum~E9$?-V00$wT7Ne;fmCc=Ivv{EudvnKz z&W~DDK}M98wviRr?wzC`X9>2t99^bfLQmHDQ;i~N)w(tV#qsGF9jk5&4!5-XL`YM* zF?&B}fl_PTI(3rTaVGOu^1TQS{l8?PFs?SY`>xJpeU2_R;_;^&X^|oi<9<%D5+&MQ z{na}Iw;!u5sv~+TTEoeUlzLaawqPUdtz(jGx?&tnYFW%px;8DQjdIT3R@Su<47#Ga zu5LU3#3QmW0k@Wm8ms!W%C++OcYCGNu5dKu^;qkXf4p3FS7P$|tL94BOAH5Oj!x_3 z2^F(2WsV<{1QPvhWY&<@luDGFU+8b{Csqg?){Dz#|Kv&Jkad@Hg*@lTKi1U~%j>lB zRr9-IwkbP?3CZX z?9OzXJt?^8r8?sROM3!aG-)cW(<6A2SupkfDK-n{GpE67g$q#*2sQ*05;-Db=r>NvLa0BAz*Yf$E1EQv z$x;=ki?1pqGk5j&Mboo(e5EOY+F;E=Uqns-*ak}AG`t3a7#!gII0CF+*9Tn6ahnPx ze(UyXhv_y&F%NK~)nf^2NT}jsY!Kb0^>59k&ibHtzZ<1Qpv23RGNr2wtJ9O~lN*p= zvBph0xJdroqpIjn@w6&SYt0sF)WtIv5a1wrqSg?UwL-M6Z;2&3hYF}?a_Jkv4;A7> zZuX?$kuxGPS1y(*5Ht4+UBMwRB-VxI6U>@k@usMf-lJ)?lO8^bMkBFpp4TNl2>V)K`6Q_ znjcrtDRdZ*-A|C9M=y>C#Z+##SGRKOXzLI)cjgahZo@t;AaKh%%L1Ku5PG`itVeu< zY(4vJ-VXI1&YUSM4G!Z%o5GiJ)>Dj$N60=&$E{VxRjR~(+2JdAUKm^wQT!u5yuE}Y z;z%!1<3-+57|#5z$E*9G523A{R=agamkzP*WjNED^dI_oMe6(f*LTto~rz&Erl0B<#byOUMnZ+D6@S?e2iE=>U7H=&(&Jnq8MZT;fQ)X+) zLq{oXk)@98ysatDDz^14qc(xpw`ON3*)5;##oZA-=^5jou4hD;WP=0=@%RLNS=WGj zBr>q>>jxD#4qV3lpAXldzrFwiv=xf?$66IB;kD936~B;-YZ(vM(dwjeHVk%Y2v&BTzX^?lib9ktI*1oB z*nDz7I4@l=$P#5`!d84#>=jZBWLLXg&T@to`91ac+Pk}}Mu1*wXPBlKY8lFnLdOeJ zSL3IvPN1xjl8V-f?h}?}3JmA6Yp=1HCTP;6S9m%dyXFSZ?RA}2vo8+O#^$(#1%R0Hsd_#r#Bv)6a)m`R@JpfRhtLS3HjR;F= zy!+dm*SZ?gUgI3RBZ|08xy&a&rGDDC)LzvO4w_D<-!;umwL37Qa-4tOo&Rv*@ED4f zIoDBTe#Xcynr!_77Medx!t_F_-I5-ovYOzHZ{*ohnE#-e!OYhR(Vuh~7bPX8rc$c{ zGQI{x`8HGaR@0N|=pnrEmh>Uu8!(1~`mH0SqfI4>ZTpqs$A4t6YF&597UyJC5*OY* zUzDsJB}@|oj!vkNd&$!=`v zQF9b>pNXeRX$OG9E0E3c&CvJaLvRxwxQt>^Yhuga~-jZ1%({`VJO&1eHr_Su0`T!pzf*#lm=>K2{c*`#g@xP9tf-hgNw{ww#!#^i0 z2N!!MY7pFu=C>;i9Gie1&;*F5cFz|mhX6M+fg5Hs*XsKL3z-o?e}f575ipq$^oN@M z#Y{kO7-kGV0=XFIf`|zE18mSie`qV*0Y(59mmD_w3P=hYDg{EYLH}hGpa*0F=B)Mo z1gC0<;0I*0RtBWSMuZbL_{IE2Rz$F>RvoY* z5V;>hp*#;vAVuW<7Ww-Aw^9~Os*ONYoUP+NnDn~%uI2j=AD zXkDB7Dz8M~NMGfN5n3JvFf8est_xW=qpaC-jVA;6^#=F}X@Ce~^$w$vO- zaG3Y+0w)kJ`1b-QBmsWVpC{n3=DZylAFj`|=8);lOM{S38f{`!vjS`Eyh7tzvAH(R zsc}I!78?{R3V!IaWRY8dNEPM_6ndlX_p>w?4(CI`B2NW=e8C0FxYDy5of0~C^d*{5 zASGoJ2RFH=dI!Wiw__tiT&+uCWLL_!Jm#s0mxASVTF3xu5mYTtSriPQiyz*Mx)}-I7{NDW_Ed}>>bCC$b{OXY zsrVSi`_$U923&L@oPvV1H@lnWFWj26T!d#6?09T64=FO%%J!7chcYFO*30`0s=j6} zO)sWk8Oq5qI^{9XG6SIa5E-ZZ61>^L-gL^rx!Q|J(fy#V=$L3;j*2xa=rY^h$GiR2 zQiJ7%vpc7|=IX)5*JJL$-T4DYE}P=GR5zfTH{)6sn~;W3i}ppy>$&pUf#OCHZQ&(R z{g_R@{o>Z>*=H8D%al|xXwx?5PjbNsN~|W_tlh*45mA6MEYwXTlsotHbBMxs zj>{peNOFZ^z>Z~+Rve?ZB2l3+$6l=8KxNix&yuemM%xW^W-pk~KWYl}la|hu76XW} z3nrRLT}K*_V_**>C%-KP1oPwtNhOS(rVM?@M9JWLI+*EsQ?u_ykYMa;MP*j_G3qRm zIM(m1jca;_WpWHlmi1|b%yX@wBg;=Ch^3Zmfh810xY?`?2a8;yw!UsP)^<1LjxDp7 zC68uNX5Cw%_bSSWWSizvFY!KhuX!tJ-kY?XxjLbwI(_z(N-C8Nx(5-p7D?jO>pK{b zSVyZ!l3TM2Ns^IV)$l23|b8~V@@_;2d`9%JoE~vlz=>M6S_~1H@-|yxW6S{JV;)H-Bfrhf_Mf+yM z5;h*z!6p~y=H>FD600Re#Rkh!T|5=5RB52N87`qX@a$Ehnc|eGxeP}fp$mc}4H{KM z*j7%>JzV(OJCu4gGQ_u>sxFUx*~Ri+^&84nwp0=G_ku+Ia4_{^)p(cry_Jr{ZJ%IQ zS!#;tmpJHtSW>#yqb~^*s~B3mU&QdGtZ^IPn6%RexzHxto8-`AZ8N8^>4^$ZO%N#`XdhQDi8zoa<(qN~Ir2ezVVlw&Ajpp5IJ zQY>as6~lmVJ*Am6p~&bVak)zO^jttt*6ae{(3p;BOB8CV}W-uZF+w3pkI6e|cdzx{PW!?En2OgRwS} zJUPiQqu~dBg^SF8T++id$xkFVPNdc*qEnNLaty|0QHf2fZ8RERN7+o(y*9VesCm6_ ztyb_%*0#7vmt0gwQADF4+UN~%n%`ybht5F`!m25ZU3G={^ zS2b_NJ|HeV86^kKvCazO2&!Jzmj@oue2#&V2n!h9OLH!u$(XcjiSWB%WNZ*vaQKId z%RAg7^Sc8Sxb9rY+fgbr9lxp2dMsfYWFk9zw<`nNOUS;Rz*neZx9B0`##_I-TTLZv-dgsw|{HxwbtIO#py2;|DG-GTaGp3gat*Lv;yu|;W}eg08j_+ zFIhM?Lh&y=aOVnu>CJ{J&4ZCLE=HAkUpq}9>3&3M;%AL%^I%Gs9*#8YY7@#ZW~#(K z_ItZ48%>TAwKPW?w=3sNP4GxeHA$Jjex*X9xz@rg1&q!~Ow_6diE1QJ-`QRnWd3cq2WFhK=saDKY4lS)%D6e z?C8pH=b}Dt#=<`0`nr7$g}MkzwBDnmr1ekC^{Y1bWqe^%sZ}{aFzHNrF_i)}H z!!M)ww6@*($d=AamuE)JKGv4h?EDUah)hAmXP@|%fs+GkP4%+((mo%n`U+csg)cBqf9Ti1)Ue!^e4D-ADdcI4gMKvEt&BF!>%p0K!>c`$=vI4M<*~hj8 z)G;+zUD3u}iX~n9LgW?tm*kI4Ev+dlYE!rlq;%D_f@FzI4CwFoMY@f7o*pzD$vpV* z%t6E3u%xeh*xSq1x<0P*$uQ~md1aQc8#XATd7;;4+dl}w?iF6(w@1SB( z`-e}i_oVo>HBQr&yAS)#V&q6f!S zN9S!PJ_U7LB+W9w5`r8;CWzodv<<#?&<(7WOx;ub#=L|AtAS@ z>^Q54_09R?7T4R=#`ht2-e-NS&x8iLjATQ-WURnO3Oo;3Q6G0;^WBp!d|Bb0V?)=r zXNg1i(*^><;f1~4wS-4`u4|@~lr>?M+)`Z`Ti*W9GoB#)H;$(+N_x$C4#E39De#s9 z`!BYycJcP=#cZ6b8{<-FEQ4@8&&(&|70rVnePuj!Hur;={peIx(bwi-tmpo%P3U8fg;?m8DXwKnq?Gz`)=70=x91) zYWKoISk{*_mBB${MQ>@abWy*ZVW0)z@l<6=^tRAA`qX6_Z!cG@E1TdHSe8qCS5D9h z8cNS=kIgv7xDu4(1%2U*VyQ?yjY12w6x|!n>ko$uTBj0cunNN}OegBnGTa}FtZ(Ku z-DJm~MM;ozZ8A>Q7>@(;K&*b~4cvshr&>nmn{2T1yyrI{y+&!9Il`!iGgO+lWOAuS z&cUW3((G-(eokO7qQy_RIv(-WC5<5IRu}VAtYnF`$ECS<*fV;@s4L{6l~lzjwIuas zwI?qtq{a=2IqwhO#BC+Cw}*>yb~u-ymkBA&P)p?!(!IQ2=kJ8y-Ml*)#kJn_v%Se< zyKp=5)kmtPZLv`4BbHjdA8lFmZ-PQ?&WQ-OJ`VLkE@ZM2tz?y~8`265fOo$Z6_T|F z?cuoK&h!iJ`B6KFRx*>wj34Bk8wC!7slK-p4WqdA8vWKjJ z`CBs0S!Pm($43jzfS&td8PO5gYXieR-k;7&@~39g*1t84d`Kj3Yn=Y5W|VkpdEed* zwG<8VY~O+(PRdyw^4Pe4<}t6Gc`2 zBLev>#s;yoDBi*u*B)Znn2%yNvBp*+iXz1+{bP!*jQ5>-4!*K^89=h+1WUjNQeI~A zENeXfk}Q}L&U~8-fo#IqVyeukOVOK*C$EjhU&hc>TYq+q);#o`8caLZ1~opN^)Zp? zeCCP%VJc}x?;#gm$9Sb;#%6?gkO8GOy1^4e=nuE_N;K`%Z|y)2$~i0RWWC5 zG#j<$S_^o;nLLNSD^2kjlpq*UAVqanhgGCDe))0MeRCDyH`EQNiDMaUIwqp9XL2kj zn(~r)Buk&uk3rnS6 zzt#z2LrE&ViYEPvJLOFXTdk-~lD#EcBc>WY+se+k%Ft_1yGoOFt*(F9j7u-Xs-gLw z!2-+{oQd-LP+889DtLSUhT$z#+JayFQ(=fdUDrtsn%r!p^RNZy7S9^YTry><*SJ6IBD2)P9x7n zt#O+_4D3{EJjpzJ<*#0icA-ygj}K7p|F&E$BRld=FP}A8MmOx??dKdNsI$12ssSb$ z^!7ywgtSNjc5J|8AgC5celwaj{GIk}_d)Mwx9GNLnxQXA_bl5R((O^7fJD?lazZro z`Dw27irTa=&Ue#Va*z0XuZrU%EN#^T_Ri%a+uWbu)VJlb`;M%FmEptn{SNj(h;&^W2??%&i z@xd}3%(D5uq4$p~A_%#J1{`kgeKXD@3ye7h-S^sC@vYYY(UmSyp&lu=aR_8c9+-s) z4bnXJdtbenTB053D$OHae$Q#K)>DCN=@#XKo?4;w+vU$XHIuW%awtk?{Oo+PtXXmb z>=@&?4jeAE&px*g(ev2B3@N8763sJAHtlVZz(VaRy^bJB>u$dU`dVV34(HU%soaK- zy-n=Jh35Br?A;Tz#ZaL>T`9drQQEw89hKJZT?BgdDtQ+2@AIhJyogP^V?j&>%rM_2 zwMN)vvgddowzmwq@x6CjlG;i0Cg&t|g1Fu^+Rv;Q{p@+|xF(;i=oHBl5tVRa`~TSH z;xn!whsnI86l7F{PJB!l5jhnVjnFS!NrVJ^R4Q0A8HrTp?nZghdOwV9UC6ai{_t~0 zbmZYTdh4wfJ@4_&bu6Rv?>;O(kCHo(Vcch8rO=#0hkMT+IrJ`%s2u5pc_Y@hiASgf z%PGfmn}VS5o!W=qAKJ20xqtrjPF{M8QrH%7_;KZxAKVH8K=EC(rblkNjAX>|pJ8_jQ{0fx^BDvv<*`VYWXT;$exeu&ELZa>(I~0T)(Fw=RLG|z zwXGX}B6)*6${Q%AcVDZJ*_t!MUAv5nYpLdaV$}q)n!KQ@S+4fum30eVTN&rq ziYLpuc2jZMvo((f;BcI`sRs82Zs;56UDa@J#0gBO>uJ6F?marGscK$L;;9%m*b*VD zqx7Xc2)S-OAUI!dB)<8SG4>0~3yNVUgM7~;v9S^0<8<9U?>&?C3-RYKv{ld()S`wa z2@W6l9Cd|!QD)e+8W8z6KPk=xPvgp^=gU`KhSo53Rag}te9io9JGsWR@}y8&sJ|96 z!g@87(ZSevFk}^0y3L=O;_exfs?mI3Sa7}t#T0xLweI@zj^B&h8?S&MeoE%xeV5yJ z8b#(`>f34(jJ6?ZonG#Cl0A9;Bnm0OBx%nMvU>D7E(pa^J@IBj9-@nP zykQY#pELaDe`e6ciFBT}5*UT@j!@__KkMpLe(n9-#$HbNgoiZ_D*;H@{YdJ(5Edz3&qwXs852w$hPrzlVP;bg@J*o z>}_YZX_+uVWrdYh0PfftEvX>iFyyW?v&fh}>eZS^xF_`Ry3BD)!ebh?d2qn9Mv%pM(2Dg2cJyZ_bB!E&lGypGH42G!x>7qiieAG@7Y-^ z+#C(IueKdqJVJY#va{_4owpIF|ZeY)bmxq z$0ZGn|K>K>w_n<>MA|yw+lQ7e!iMXG%l$ycC!z6fois6f=RX3z!|??2z8g@neS6Un z%uI(G$a5d%Q{gE{mo0-GV;PAZ03+Y&G@@SP#@yEXxj#cKVQNQl5@lq*xGP+FqvkCw zyZ!rH@(D1v>q`#%3#Z_97_`%He66>!mcqli)M);7M#1blG+2N3!6(QaBF$CBmh0l| zKCw^W7CFSzpxt$~AliWKmp_W`SH(m!OVb>qNQvwlc4P$Hq|5B?SDYB>#I4lOY1=)# zzjDXV?d+>u6SDOLBcUh}s|* zr5~T|cguc9cg%5Xq76IVqvvdTHbmV!50JDWAn2*^&hg7Vai*Z4rsl@HuIS*WPmkv|5^-^*b8Nlj&0Aq_5Yg~p zo5=a}a94(6xW|xVspt|>Hn!rTbyTVG7$s-%Zs#+z7ZU|b=<84_{)mpXUhERiU5~17 zJ?f>md>CGO(Lc`1z^}+<3+C*FK)ubIg6Up#yd!HljA)3R)P0Y1{9I=uHJf7iiej?A ztevm`!JFk{@GiuhvWe(!x>CNFE2-GR7Eh)JmXwJ<9~qEuw9QDaW^~g>*(95AwTb8hO`t zuWv=&bG}CJ=>cwI!Q<`_u~ltl)K0@H^4t&POAp?_4;$fUKj# zA2DSygdk9~uOxuI>Ld?O4;u@NQ53MszX)$u4omR&ebg}?w-Q>t2v4wTN9q>eL$t{; z9=J;R-L21x=sKxe;3~w}tyFn5*RLsHHv zB3aTEERz5B1bzLGRNz0lgyl#<|GlfqQ;w7ti`~o=2{aKz0{MVID3T8XghKdCz(6RU z1?C3?gfN3cA+V|ndC~za96=a{(_2MSZY*e3g(B&@NwhHm0I%u`7GMCYN}L1m92@hX zoC{!$e-j4(EeZ|eL zSnbKxR``?a14WND^Sbf}0*j{GXgzsy{*zAM2{!C&zuhO}H+C9Om#DR@ka|x1yE$Y6 z4cDeLX>%xpt26G*N_G0^tM$XZwOXI^shtOoN}FfrRt~6%)shB;0Aa{CeuDy%=CyBH z_uKGQKNNpjjx8+`8Rx8PNx44kN9A;{pt4I?7%uTsSlK^!zdKg9s*D2{W|FIUeQK|8 zxr_Mh!4^2^Yn$ufNrIX)3OQKo%+jpjF0?wM8szqGf=7(fzwjYCE&JGS>BC=dW?ld2 zbx3TrZ?F-!Fl%fua|?YgB;M~(_a_!d+?%^XNm(&igw-E%T>MyVa#nv>>x^-%3@*9I zIhrdyjlp|%hdZ%8bc`O;^O@i`dmPO>vej;;7spXA)Ek^1=-^Bi2C?w^uZq)FX#++lJr^lqA2ZcN>*DRE7{Yx>T!)B9bt9e5yyR^~w$CzxPY zH|TRJuJi(>GTijZl$rB98g@@qz~14#nd!`P-2&S9>YU0qsE}QAUat&ckOI%CEchdS zzDhHvD$!>5N0}y2xi0ZZVF_2yce-4pzh*~5bK#w7a^8h|%H^^q?WZ}s-EDWPWC8(8 z95x$MewyWd&-%p8o-!U32fTmJ!DVI}OlfVp#Fby!?tD$5rlg(gnuMsOP2`DWLNRq( z7`+#NttO3LL77|0K9zNsC7vr=xn4==56t!C1ARG#}7=@P9B2enxx0L9%*-O4UA@7zn*iPBbT)m^b(AvS_HuG5rRbR{yXciDw|I}hPe)!KnZQ>#_-4N-J4 zFs~T<43s+$;d0~s`k&LEvz|ad%qWH3SQP*MJ$xCzB0xEs&LVBV&xKDhmIlv_8qe%@ z-=wDrF}l=T>9TJV-m(m2&)gX7^t`PngL*;7qSUeqSt5y-^1!!X=ButreCR>PqYJeKYVvH4!0 zHgc}(35*OKzb7`o7bt(tzQ#2uc2_I7hnC}B!f{2)BMadMJ{4LAy>HryA$B2jI_s|H zvp0+5l|W&h!$GwNw+nNC5ySo~L6>WMf~}?rC=%p8%cCSskK~5s1`6i5R6v$D!bz@R zxs1NuO4NWMzb(-q+|sAlEG1Yeq)a4F&J1e^d5dC^eH`pYmpFpW__i?TaJ&^sgba7gm^_?GsB*Hx~^rHQuxHWXTLPmqjJWo!{?QWy=<}UgUaD+L}vm8pO#Nn`>E`H z7)G$~J_eh7-803N8maV;I-^~zuc)wP1BV!s6cXX>e{vAHj1Z;9ljtY=_@KT)Hj47o z`{-ors4Fd@GvtFE zznyhhf%Ofd*r^9tlGS#jbpr7v+xcaScw2tx4E`unX0pG>)MNH-N5%7rX*yuFuDCNZ z*~yy!Xn?mIE?H>VG@ak#jLkkb^r;PH&znf8L|8MH&$i26q|7Y+o#L8+YK{zujQg&I ze-V7ZpjKP9lDRbi^n}AJU%T2rEM^?b*s47f^AUPX@%m!4(VpW~(jk`z+RDTegx5Mx zq8HzzCBVJ^c8^Sje|MLtn5Ce;oOJdLc9+Qwt zS6%#YQ8g2@LMLh_t)}B!UYA|w)e0b%q(&ym2TguHawQoab*nJ)TbA=~=rAZBTsgq5 z&+brnPr?2Ix1jCHa#_3$Ow~d0saCZgXgo}3PNeEvvmRj?)80Ej)Fj37|4GTc!q&)^ z$Zy_{=i6k1MU7`~wT+9a_*T5Gkmp2v(I$VHdw-Y(_vO&i&@(0MtCn0bc+()SEF-{f zN6-!!NpkuuKaNlc6o5)l1wZi$ueJJN zTI6>~4Q~H@LdeFY)pv8sIdc;GK;9TagQKCf_z1?ECvgOe!_td0OnCb(s`r7Wmd^5%?8F6S%|~0^ zYcpqg@-%6d_&wIylKiafDEa}Oh>EO)O6-K&6K=~E{Z@~sC;r#{o~ z6y?3KZ`ISt6S1Yq1Zx1Ym@NCPk>fs=RuJV#^@jD=ny6s;AN~jRo6FB9 z1-=EkswxNMpht zS!n@WuW{mgGm)uI@EYYcQ#{US*--fo#gQ^%arIG^Hx5l#B6#+luq_Gy(2Xa#32JvoioQ zbV*C8^86z>HUjD6wbe`Gt=0N`WYvP{KP5kCYDe8PS(Z*a#_B<30Q>VQzsn3Cq(1l| zE%TFERjgb*cV|L+Xp*W(#?9)3)7TZF#gp=Vn`Ac2uv@1gLC}6;(7hgRpL0j zPR3kzI1eu+oXYPUEc@%RwFVORi51cyCfm5mipC(l2YOduO!XW-ciz zd*XX$vl%zR!Y!U6DP2BIPtmpHDhC&8wsKWyT*@Q!$FHG}B9x|wJ&Z1!Na6Fz1a0I( z-Rf*NzdpVy@BFxL)sp)ygWk%tZ&>_saUZ*s8RK4&iFrnGlS(}Yiz!MP&`b{vr}LiJ zlQk)|2<+4#Iz$Cnn}<-P9zs^<>iH*9p+4F1JH#{&l#8`3y@QW_M2CI(I8%MP6hG`< z_LRle>ux5plYRB#>U17<$DrB-l%oJ0$?;v79dymFG@C%u82}PyNsRVKHw4Qx5^}1- z`#nWIRyzBQrDAVB(nZiTl1Ioc1Egk&< zOAMvECoaQv{ig-fiFoXBq1mJ@u!&rMD)8#OUD(E}E7!eTrrS|I=2>!-*%hU9Q!TPO zjJ6uLipI4~0?H#55E+&qkIg*MCviLFd%VxR^Lrf@$^bj{-z+LDZp8}K+XCxJ?ICY z?!J#{QD2NaDzrEwKIHypIEx>^KbVWYTx zEe$Q|1~dirX=`5jl(cP4AM+RB0WVF+Z8m=36HbyLSpp_s%>O#5(_a{K#lTd%WV6*| zAB@f&nO+*VpS3#M6wIL68y9j04;`|bQ zAJr-@5iK>oRkn(S`xku4`PaCPGQm7SsKexOouO}R8!M0A7$j$Dg7XeP&yQQb z0=FCs9h_Q$sF{4PWQREV%lx=$kY{aW+SBt1 zvG0CgKoyAW?U;mB_oLUg{WQZbmwNA@$_R~DZn;0z)RXbB=OSEtpN?Oa9-(4Z!@$cu zo{svW1@d2+3<>UEdy#v-?%Dmu_bXzXXhefs8mpfAut+>es)3Y|!)-e8flrRH{H>%H zFhSD_hOI?Wsa%eDuTrnJ3};HWyCu6e(}*OlE*{ml16)d3@lHP)_`+qcX(_}$-inVr z(Vg1~vxG*LHahXuzW5MsNq>jsd}o&-_V^~sk3oG+I%kJ+;Q)wuz1!O^yTIG38D>|8 z?NTScJc%?5p%gN@)JQI+#J{l0Y&y8T+*4!z47e-e_$hHAn?H;4UQni8Gb8iXq?xaE zyvY|J!IxtF(eYm2f^4QAxccMerc+yqVo&I%NbhX--ZIqeVyPF5gkugvmg--M1!MQ;ckyt)4BrK3&fU*ethum{Ju7q^LF8(78IW)_vs2L=H4SZL%`8ipERM@I5-u%= zk|K{Xhw~S3XA}=c3EG)wM)L^+N`BlBI7=0VqX&MPsBYdcN0^%jDZJHGDbq9B{&7}9 z{FZMhg5P*Pwk#2%o`kSQcFIR+WRx_cu%d+>D`dYZnXfz7HI?+JHQ$DS2bF=vUn0jJ zHosu3IqY=H3L||Kt71`=8;U*Xn%G8F!2)*j5vsAvG$M_2HLkvSI2i^&C(z|wQ#F|L zl)n%bo4B1)0@t?{2BIpyPNtgrwCK9gSt}m%in<~T67Ge=nsK$VQ4JsKICb&Vc1ZemDEg_@~nvJb24#+ zrG?pc9FsIxU6Qxj@^VtCKYJksYR1@9)yOi{>*Q8{i>S6`gg~&n7xVH%LJ}a`@LgVt z$BA3ApXPTN&UCBlJzr+*i48n_;8VGMCOoI4@=797&fOWQnY)j>kGj-9yPO+Z>3?(M zG-mNSHV*YM$r@*TdZ&T^`fNkMW2mm$gS^n;biX`$3m=(HGH$MSMKUy5fx2Og5NLM{ z+V&%9mHULy;Szf=#vXI2Auv_=!O|L&aEh^*)k!+$evuzE`R23|%}>kX&vE<~N4JDi z_(Vy@Pycm_Zdzl#K2N?4qNXc%(4=UDrF_(waVGvtixN%-M{lx5YysCL`eps+B!1`v zq3@gBTB?HX68Ms>nB24pF$(qhMj7?jUA;!FfQZ-0Znt!bIEa<}Z!}bmi@Y2H#U+O5 z=a4W}1XECsMRs( z$ji&q!a-A@<w`oTk~IJaMEXUdx}7+MJlxI+c|E9a|j{@l^^nW_16RuBw|_Co^1 zQ1Cz5UokO4a|biWr#6J3-#x$X^jSgBKY3jc^iMh$1pku>2EqR%fI-kdSzQqHPu>{> z{gZA6LH|q$g8o4VLI0$bLGZuu;D6!49!0^9V zL(xDXfCvj13g-u6lHO1Rj30tw4*|mXkuV4=6p7%60O9CsVSqFX5GE)9hArgscLj?JOQ0Q-l^#5;$Xh#VEU9=G!9{Sq@01qP+$S(+I zg&=|a2oUJs(gMN!aM&+)GCV*8PymQPg5f|gm>&!TA^9=ujfBI1|0yqK?%yf|L*STK z1X&>%5s@$iD+G#BB?QU}fdMgu1%EHbKTx6kARq_`M1cQ+3j2-93>E}o6!>dtG3b~D z!X&!+5tu><1Vadj_yvmKvof>ew=!`svoq&+Fn19EVWzOAb-KZ zp!{G2h!p|{W3&ll6@h$({zf(7|81%Thu|CcF(VPFXIA4teQ693YlAOgvXp~VkHLjP-V zz(_&U^??zXQJ8!_Mr$fe|38es zyy9Or2Zn<%gt05|Z&7Wr zV854|Nf6_B2+S)OyT?4j7z`MRfybx^4CMc#TeP7#;2!!_`5(sI=mv4Xi`y6%1+xB% z{J%<&(yW;HkJ&&}Bmg|9-(|}Gc=D%6sH3a*Us{Jl!2B>wq;_z%vtu>}83{k0$bzX2d1jC=@I7#NCKEXY3re(zWR<8UYf6F@OL&p+qKc+CG| zdZ-|fAB5SN{xKfoyZ?*f5Cq1}!T%f$0-=8Gg8w6Z2n4g|{omC8d)NM-hXXOu8DmfX zV2=c%{x{PL!Z2F}Mp6G555km;{}=QHF~*8S{$J#W3;wU>2g5M?0R;Ka`F|fe{zvh_ zAVGdXjG6s|K86yyUIG9BVjSvs<^H$H|7CHQlMp`?vj_iU)UN~9g#a`T06dTh7qF;Rl@S?3A0zOb- z9{!&7ZzmK?xc>b?kC+?)jdc@@K))&hP@#bp03!6PD&RhvSpjgL8-&?$ez7%V70__h zbg(vcG&9G9ia*Dvzk><>*KlfdV-SD@ovi?1C&iHXBdj*-FP@RuWV#q)5HwT~z=^qm zRbdFX-b4dM0k_eregJy(V`V@BhU^~|{T=^L)wqgV1rUveDV_f}uF4pie@=-1LG#aP z@9!}Gmhvx}VQRl={`C(J7@8P%f6+vu8Px&C*dXLzdi|mf2*yRrBmc$KMF)^fg?aj) zibejVSY=&66&BXN74w%~AL#?Q7(tj_=AU}S*a)T*JuL)aLHlCVkH)-9gWk~x07yZ| ze`^x)N0VqF2>>D5(E#ujZ7l*IMcXO@AXTo008v78y70e9ezyc9-2!0_AisdGL;-YY zyH)^WRkao13JVMT-}QJOeQoob+y70OnBDak#JJeMxV^OnaA0WtTXV=inqvkb{}re7 zS&^9V|Le^0>oUYBmsKCFY6oCK&&>ny(e`$LL~J;wi_9MI2+I%+dI+HYEd@3N0Y?|> zVNszc?Eyq6jp}m0hNM8RPn9;UZ-i_$rLm~U?y};T0YmN(++h#FQp1}d2GjnK5dT4k zBTa};jr$&_P`!0xbo6|3s&-TO!Etu|4vzPjzI*+VgYiAxPOZpW zxTVKx@q7Qqr=hWO=3Za>3@OS@E+>7|a~mpNkHdjFsE^PRv(RxBad*0N^46C2A#3 z@T#1^7%_8kzBF1Jt|R^<<|DBqtfM>T(tO!hO04=4g;}Gjwp^w*OclwK;`ZrEIRZI+ z!(v@AlPd4kZsqvLAk8Io)gN;5bbTmI3sZNC@fDJtcuX3@wk$2PM{u@O@@Dd5Vap2} z0V)|3^uTSmZH2Dq-J-x^ZU42?%b^S5L$pgsSekNC-wx626ILCv&&GC54^+9&ru5tA zyh@$cn*E8|X(sqo)H`hfYJt8dJ!|@_ z+Rq=XI+NQ?D-oOKwG=K;7cRgEk)9t<>3E7!xM9{hVkIP@8z&X*Nznk0i%UjWZlfaq zVrup{?Wg7))>7K*u;}-R1BVPwB*4YXjjLO$zRFu~o;#bv3Q=zsl|D{W&?q}J7qV^7 zNkj{o&7`N4$lr>)BT8!#%htlp%*N^Z+;I$gCgZjo!Fkfn6;c0{5x05=T)ViP6N*pBQZjcS;!5D z$1gLz<=&)KycI!u=RMs_>2*-Q)CqkNwJX#es{N$1CN3%|hC~;*KDewiuA`%g(o=u4 zKDDo(Mzy`OHH{_rN#~7wXo9cmN7D2ae`tM&-W~i-SJ#S2yI>K4s<3CfJ zZiKj_tCloVFQ*@G-}vT};S)W!-M0*dg^iK>YBwyE?qG;I+-#rpS>8YHx}GPT9*us z8k9e-FtVrINZgFAVY}9eEl|+PcEB6HRe#TrNmP3s?^%0)Y<7$U!EpPiv?hx4(@6~B zJa8OmLOA2YlkEi=53FIynY9}HuLMZ%ANn10xnNQddC?QV6U}j8L zdEIp9ryWsNt7Pp`l`IdB5%(8A#m*x-tD&|x=WC1cS9?A}Q&NoJL^`;+-INq9h5FA{km5$aUWvU<$H<80Qw%gtWU4HIFME&yFMkXJ(q_X zfiuV&2#3>TlteYnPV0AkcU90*r?yJv&r@(#GuSaQ2K@ZaxV?V=FtR2^h`FR{rsO`Y zJ(c4uK{t7=pHIrG_^R$ z`LYLDfIT)7Q&Li^+j@(0X@IyY>#jVdrTn|34S%6|4&s=fZt|m>E3;QPb>fLkI$;KO ztFjG>^GUlmIP@LGRvB#LQQRIuN+=78*peo-{YKM1QtD#t<{hfb_z`W7+=9}qlC)I1 zEG$@(p3315Dv|cUe$P+45@KTEpvSw@IFj(e2}TRT&Z&sda^@;Rqv9=7kiiE0nRKSHC6;LCn>Br-@XyzpiOq90Q zLY!26AI7JAdJ&yDBa6)TSN5L6(rex_e`{pt8$nV&Yd+{$B(-TY=oqSqx*3*S0tuSv zN0lYwNg}lkZk78QHnDR(W?zxcNmiV{`^up~+j~p$d3%+nPxXl04w!hp4)nrCE`{mE z%Z1z~MHu#Vw5RGB>hpaRH&Sp-#7t63C-h0?!+O>^VlsEa&;2jUzT%3jG^(1QD8VOF zC(AYYPS{QLKB5+$9mF((y{f1z*-&C5BSJn>oCi*^9#dTgF?|DJU$FsSnzL8mm9-F8 zQJgt3!*15#(wo|}1b?D7kHo3_wlXc*)Mf64J3zU)gaZ(`MDd_31iaC?-V#Jz1*JVr zv|n~7>)q$D(1l(7{lS#%to|*hetIiKmy2TkSQU~!HZoo^UQ!PFTwRIjDjhGHZ>ro% z{`toJ77v-$$uSp0hQ9nYI@3TzJNY`*xuAjN=op|uOfFGP|7?P2?UE#DeE!aC&p!Hv zNbMm~X>9Hb6eTtY(mXNra#f(xXg&%ya?alPY<@|@o;{7!u&RuyZ@ztQCDKS|MldLd zMbQ_W&`$S^A-Gqev1Su^MZD^KqM&8_>E7d#XKm`O-fvV6N{-a)p0+*XbK(hRYfk%K z@;+RnF82vt1E0&Mj{@l|K)>3mcPL2@tSKAfX}<|Pp*tvuEEQeB{sC##Yt!$V9G@O~ z37b3(QlTYpfE>fLT}DZEYXvX*U!7}g&}27EM7g#c#FxbqP&ren*BQQbkFfm8_d#(H zWbE^kD{0}qby692XwMfo?&r9pYyk-M=(5(`KI5OQl#MC9I?rkvtDgz1UQXYghXoQ1 zdY`>Ms9YgNRfS&$nqDnC5}Xskk}h!tpKrjh3CtvR+0R9Dx0+MN7v~Q49{9gUV~$g^ z_4<@d;cs6meR`Dpikp^jhsM2ny3%f3N#RTD^J8T*i}jko;c0tFjLv<04t#5}Zb!7XG4lnBp z?R=j|WUkG=!P`&Yg&eZ>Y#DnSxt_(zzP6uvS{wD^5{vk|-Ua!$#~HKd#7A8ncx;P^ z_jZj@Bq|}MP`~C&GyBcN4K&d2;jWD&v%?bAd-m~X zczx{2V%|b+qTl8o%~LIqv$40aNC_yn1&PnqvcUw0%ZR^pi&giwa0nEi6Nx$0bZz3j zdu!^`VM!^>9(?7hrlP2#_{jZdXw3AIK}R4;ek`9flYDKE7*0ZAnKSrQ-Dkh0 ztE1^?*W}vtrn-t#El(awR(O=9x9Y}oMBcN*a+kvggw)@&2{_}1KaWU<2FwtxX}6%Q z7PW|?&ArCn#cpo0J9304neTI0K9Ty&Wdv|(OE}mdciQQ_`boI!z3CY4s~iaFkdxNa z9U8<(^?xz$(@I*?MESdwRi@lDD!~0P12a9vb=sa*^y0Bx{mc~u;jz3e`lIk`CAP)XjvA}ArE?EB7vLz{(&94u5SG4W8Ee}25yw~_dgvUMbD*lXJ#nE#F z0jl(_D^X2ItBbwqV-p50_3E1Vkn%*`L!6h0rrL*n$-Uz-SQI`AFPoxYt!X>qPO^Q# zR$a#7D|puBild02U;4hev>*^%9s)ozDBwu|LbczZDvBS2@U3oq>1)sm2k#74-mN9m zS2y3`0iNy1zdZQL8UTIr_CZp~G@h0L(H6>ZC2j|w(eDvsh7L1LOzMq$9_A-23UZs~ z{3LJr!I5q>ZcvHThYDY`f2i2?*O=IYY0}$S*VQ_1CzCny!X51aV{hOWDuB84(VwH?DE z?xYVQlY6G7Y12xtfR*=dd6S(BM_-d{rQA^@334r7N!{!|S*2C`cyIEd17UQrsTPWU zyX-9ft~*{Bz&t*y-yLqadEYq_sqzR_=Q;hA1vTYEDICsRll$qka=)kvRrX_(f`ei5 z@$-%l!x$Y;YPySLg4sAthTF6E4O49R4ltjg_ zh7M}+9^H6|-$7JB>Ja0T`30+fH#V!xBekIDP8s%c2woL)@A0rl=IIUaK%}%4il0{U zrCt%oxulhh65j$Uy!!Vmy;opKnhNA`*gH9;7 zZw@(ez}dAZ?rR5lNaIEQ!v9IRhT~A z-t@u_?EY>0WHA-9CxyGWdkFF?zILMVDe~ZtcUGHW2FEi+3kxT&|Ku*iL17_qhz!;G zzANj9+FsLoIo(VW9CUN^BmS5(bvnc<7$fcVV~VUy5N zlPb3XtBsjXxMxlOGTl+zJS)K_Af~N?bHcA$PSz{QQ?8V zMX*YtJcE%n#UcrqF^$bmC^5a@WF_|TyyIDT8qq_FH?{;D?6vfZ^e8;WmVhq~vy$*u zXD}c4Gt908ri60o=8y;v;4MN!Lvih$us<+1`jsn_G&$`d!VSCmIWnO20(%j|v!qnh z*Xr+300+r$UEDINCDkZ&U88o`|H?~Z*#^y}W^e}&+G z<`puMS0E!LrMWrP!9E#1<4)chI`C5BrP`G=IOBP+R-#ImLZ<{u&xUM(iCDE1#_fX* zAA~0=URLWHQwy>d=BueA{J=4ZDig9jyN0 zwk^{xaN2@GaUntZnAN>tg7|S6Yz~{>%)avGLzg$y+r(tjZzE^FGdy>rlxMqQj@Ceb z^RZyeE#6_)>D3)V-9PRK3-V@nWhYLq4sYK3QuqyKdw~D4)xl-a%v#0GHJw*JhFe;O zu?~l^b24J#03k3rCb>=K7S(9hTs>!&E&qQCIq#?@wyup!=ry5tX#uGrBoG2hy;LdE zM5;7t(u+tB(ghR*AylafO7B%lDAJTJA|N1Q5T&YgHcVcjCpT!V)eSXe|h5mzLd)_L)#m&eD>AyS`H*dA@MS`U|@Iy zlW9znK7@gI^xQQHm(sxtF_o9Zd%x(b87E zW~7Pr3-YP=W7Ji5{fHgjJ*|=?1q>%A=rIw)} z#FM*IZwTutT&9_?9v^+#aEvv&>ij^F!_V*IaP0Qt#L(;9k1LOkS0B}0yE$0#sy|nu zI8~ZObg=AViu<&R%*yn{$Ggub?HeAaQR4@A$5ewH!|8K!TCywjG2qOx>d{S}1{!wE z#%k-y!O+euY#?VczmQ55i>cK%;!ORyrMbTio_Iy3C{a+YaeBR??i9YKjL|_1?h<86 zsjX~72YN2@4K~`Uym_`U)@rQ0N%#E4swGY~x<2w26>uA5c}@3PqKv=H&P2=PMOE5X z!V+t@eRN02iFgC0we(SgpXs$}2mIV&>%E#*Bc~QQo@X^RjdQbdW)&myT3M-kB#QSu z?_=LZfzAJRC9+{nqrvp_@P7;8SON_&4>qA1P`3db@k^-x+kx!Zjc72K8vd6cF7sOu zcfAazXmpDKD-nR;s9$trIu=|a2!czU{tvns;Qs^*2*Av7MnOVqz_&p309XNJwf{pq zui3&x5X@LZe@5x)kI)8D870hejquLna3u<9lN`8d4}y)I^Wcv_Q2!$j!A}1U zDsT}13fxOz2=y79`;``n|E5JkrZHm)ya~XiGlHVnIJgW>W&9YUBQYn02vH;cN{7T*Is$>O)=?2bSebwOf`Agz)7bzZeb!oH zW1<0^VMhe90|D&2p9rFP+R_t4j89va-);BrmYw8Rzt8V>hy-GAI*yAJ!U;l1VdY67 zm$A=CAqW7V&yzyT2qds*GKc{XLp&Me54eM~<1KG{n4OKMud|0cOhO2Nc)wOJ_>Wx( zU>9&cE^-L`HwqK<1S)x9lBfN@mh!7Z3g|e)?FQU}K0wtka9&{;^0ZS1n@$dC5k#CW z^iODM@8jzOSP_0*!GF6L5z^Q|3J8Ms_qoM&ZQKBJ7vNdIKBs_`oYtEv-HyvGKL%7E}SMPZ22ipF6Q8G!PrC2r&%@H2|-~)ttS3 ze7U88x%34I&J!U3f6=e{&lNK_kGVeDpP?PL3U38fHeR(etDQ4#pd<9UY2v>8nu+4@ ze)4Lz#+@R%qlG=zjl0Du!GUtB?(FIoFd6BuU%qS%BcbIZEbSX=A2OAH`Ww#oQaZ|j15O8F(Evh1>~-uxj>CHW0asaW+wz>Nw=$E|dCDy%0YJ#D|*euAYr#fei5M1;_Ns7%R9IYiP;rsz<} z$WwLk;Z_3FWGP?3)*o{|)1KS!8NWZSJDkKKp~tLTfMCxuiFYxy)6?rFU*w`q@uGT3 zIP)T39b32IR@MC`T^3r4XK*9D|9iOUR!Fu=f&|s=bMhp9K9K8OxwH=Cq){_l(W%?6 z$iU;x8S->7t^(!O}j^HElUYg=2WYj($*#!ioOn<`*%B z@*IgV;b1YYcKIAk*;O@{hu;h1mdZG$XdJ)$FBeKC&qjku@s*Q>-#pVT36vK=J%JGs z!(;qiN>A9tM=A;JUWb{oaZvrJUR!2NsolOT&8EN+>>XtpG_x^#Rc__kz$f%OWa|`K z>q6n`w{fFe(m4wvoh?Jwjp{BP_k=UY48N`@vplc%)4z4i1heZR+99jZGvSrp`>Z~d z^_kk!Ph{1sl{ql{uAgUGtsneReQ9Af?80u%1~~hAk-x~!{4?gGvRwXmA#y^hKcAp` zyo#S5)ofg}2}uq%_j~DlQS9RV_3qT^l``Q+(pGu86!(Wvd_7;SjhJ}rlnp5v$yzl- zTPxTFeFr36Ll-e+eR#DeYLBB03JvWH>%{R|TA_HRH7uJU@}5KGkU!>&Uvb zG-BK*3DFM}Jl>klgj(x(E1U>nQaAG={F_J4r|+_eunx9;2wp}KdDN`pa$XIJ()r>h zigN=6Ut%iO7IrLccGew7NRS7+@W%(&f9^F_?r#q0FXUo;F)zP9$uX`ia|Zhir_bTI z+=%UeqfrnuYPo@mokO(nXQki*35IiZIa@UteG^f+x`(Un#qEMmLltn3*zt2 zwRwm+cnz81HzU&Oh81C*imoJ&he5~~i}kvp&T>@DHIE0>O#amKBdr;BxD|^hqJ|6^ ztXrg?ms~Mb(-(0$Qd8%c-nLVntkgE3S2qb6%;O-p%c8B6VhNWeacsV{G#jyBtJV^A z73?QIqCIp650M3?mcy+&i%{U3SrR#SfGgEs9@*Y*Gq6YMms`9zKnSNw3Gw}PyAptHq_lVpy_P^ZHa+31_C<7bsEZfH+D6R&N!^hUTPxri-x zf1r;!&1WHJ0-{FH5U1U}7W3`f;3*ZrX2*ius~c@SiXuwPluYmxwPu=i|KanNMttK#beRWWyqw(>TYl~(^U*T z_s{hS$49cS_tsbf<`Np;)xGqs{9yWWq3untnPA5^(a2{GqLu+6YB+iuxjqV*OIr1<^BD=nn`2nOep6VWt!8wF>m{Xi5Yg{aEp z3ynUKu6Zpfbb^&g_O>97CU%-Gl60`oE||EazsUXzpH|FPe;!q>3|j`n9bSzCj8hwN zDyhN`(i}cQ+ba%>#(BJ~nm;FWg4d$Q*#^Y8RZpynf>9VAJhJTw z8WQ>v81I* zJ)cP{e$y8@xv#iULnM7?5_VK_zblxFJ?*`E>2Zeg@V?Kv_#V2#P*D?EVNv0qH$W;h zrSI}+P(6tfRQtsz6J$ATvLcmyyqFu+VZAPwAGEW*R=VZ?Hf^R3ZQ8E20L2vgkeQSi zYswH2C4Pq?e`bb}Im?XOd09n7_1rnHPcX%J35tA0yeMR5`C^qOuA!Xe;dp2`N=B*#*6c(&BemxlJEubg-oc>*Xe#?FwCsoQl@|69Md_fpfFdyLcCcms} zc(=DL(|tPTE8W`8y8GZnCl#Ng>f8)3Brmuaeekr5CXq?H5u?BnmbXUIIo*q$G`w#= z^QBVrO_=7#P38`0{HyC$=T9O!99ViD6{}X;^GLH8P>`akbBX*!9vUc%wIf)VK-5%w zBgc-kY+tqpo&ONhRfqXtIGwUzs40@X4)h;irX9J9yJx`RmOXaoGcVSOkjh>2lVte3 znYPH|{SI*EM4&Oeg_Kv%?1lDM>AXC$64o$A#^z~u-0t%WvIXY+EmyZcC|>Bwuj_b< z){K{ept4qtQ(jxvyhGGmUJmk%;#Jrc#dHKkK>I5@JW_~H2xH-zl!t!_7g$I=m;=i| z523{%;is=*hR%WZ+z4PBc;5nyq}+h{^t9W+!^8K~^n3OcFRtnC-~s&l_0A%$d1^Vm zgjP{gMyW_6(P))FP^t*PJ_T6pQA$#hYH+j)N(m+Z{|q5H8)J#RqXp&yp^#W54KO2y zbJA4fp*{oYR*<&34@ZZWmZ6I_@-`E&t!3V@m)jXoTfVLIv$k-+?#p()*ikhrlt4dov!rhf<{X6GC%cU*A)B& zv0JHrm7$r=!Y^K(p9;&VP4bMamX(Z6RY?-0)3NOr0MxgB2f8bC8oV>g4asiMjH1}f z+p9aXn$abrs$Uj!P4NZ;G!pWc6gNy2n&`_exI#4(Y6<(rj&8vWpzD;lah^p?8`0pa z6k)k69d8ZT0H0s-6$(vSA^r^BnCR~{wg!qEP|6j;nPr)vHtU1p!)A(__$G~e{iAg_ zN9&yYtVufYr+hj^+C3e4l#?krMaf1FZjN)S?NPDS3w*{8jv087!uX#&klEF1c@a9U zMw9FyVR`Y{wYn6-qN@pgi6);9F@+cWQsiz_ruWI}dkj=U-+ZUA7y#8uTvm2nP`X`2 z=VO50;v{I)?>+unl?*@9}tbwp(;5tX&)5)*%rV+2c}LwQN~ zypl8*kMKAPy{RsJaEV`7W@+a-!FTujufj;!on#3zb2=x6XoLJ+_T>c>XW?7h0*ygR zc=W0aNnm`}Q;)pFCFfw`Webi20=0l>Sk>n0CGdA@IEqXFEfE^ecYbAndxvS*DD0`x zmOR05^=vS&$uSX=sI_x48?$31^Lr_ipvV(CR%nSMbrZO3YE9w2dZr@wF%yJ~5+RM0 Mrr_sSy`o0(ADCk)4*&oF