# -*- coding: utf-8 -*-
# Copyright (c) 2023, PyRETIS Development Team.
# Distributed under the LGPLv2.1+ License. See LICENSE for more info.
"""This module contains functions for RETIS.
This module defines functions that are needed to perform Replica
Exchange Transition Interface Sampling (RETIS). The RETIS algorithm
was first described by van Erp [RETIS]_.
Methods defined here
~~~~~~~~~~~~~~~~~~~~
make_retis_step (:py:func:`.make_retis_step`)
Function to select and execute the RETIS move.
retis_tis_moves (:py:func:`.retis_tis_moves`)
Function to execute the TIS steps in the RETIS algorithm.
retis_moves (:py:func:`.retis_moves`)
Function to perform RETIS swapping moves - it selects what scheme
to use, i.e. ``[0^-] <-> [0^+], [1^+] <-> [2^+], ...`` or
``[0^+] <-> [1^+], [2^+] <-> [3^+], ...``.
retis_swap_wrapper (:py:func:`.retis_swap_wrapper`)
The function that actually swaps two path ensembles.
This function decides which swap will be done:
it can be a typical RETIS swap (retis_swap), or a swap between two PPTIS
ensembles with limited memory (prretis_swap).
retis_swap (:py:func:`.retis_swap`)
The function that actually swaps two path ensembles.
A swap between two typical RETIS ensembles, not PPTIS ensembles.
repptis_swap (:py:func:`.repptis_swap`)
The function that actually swaps two path ensembles.
A swap between two PPTIS ensembles with truncated paths,
so this is not the typical RETIS swap.
retis_swap_zero (:py:func:`.retis_swap_zero`)
The function that performs the swapping for the
``[0^-] <-> [0^+]`` swap.
high_acc_swap (:py:func:`.high_acc_wap`)
The function coputes if a path generated via SS can be accepted
for swapping in accordance to super detail balance.
References
~~~~~~~~~~
.. [RETIS] T. S. van Erp,
Phys. Rev. Lett. 98, 26830 (2007),
http://dx.doi.org/10.1103/PhysRevLett.98.268301
"""
import copy
import logging
from pyretis.core.common import (null_move,
compute_weight,
priority_checker)
from pyretis.core.tis import make_tis, paste_paths
from pyretis.inout.restart import write_ensemble_restart
logger = logging.getLogger(__name__) # pylint: disable=invalid-name
logger.addHandler(logging.NullHandler())
# pylint: disable=W0106
__all__ = ['high_acc_swap',
'make_retis_step',
'repptis_swap',
'retis_moves',
'retis_swap_wrapper',
'retis_swap',
'retis_swap_zero']
[docs]def make_retis_step(ensembles, rgen, settings, cycle):
"""Determine and execute the appropriate RETIS move.
Here we will determine what kind of RETIS moves we should do.
We have two options:
1) Do the RETIS swapping moves. This is done by calling
:py:func:`.retis_moves`.
2) Do TIS moves, either for all ensembles or for just one, based on
values of relative shoot frequencies. This is done by calling
:py:func:`.retis_tis_moves`.
This function will just determine and execute the appropriate move
(1 or 2) based on the given swapping frequencies in the `settings`
and drawing a random number from the random number generator `rgen`.
Parameters
----------
ensembles : list of dicts of objects
This is a list of the ensembles we are using in the RETIS method.
It contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list.
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
rgen : object like :py:class:`.RandomGenerator`
This is a random generator. Here we assume that we can call
`rgen.rand()` to draw random uniform numbers.
settings : dict
Contains the settings for the simulation.
cycle : integer
The current cycle number.
Returns
-------
out : generator
This generator yields the results after performing the
RETIS moves.
"""
prio_skip = priority_checker(ensembles, settings)
swap_freq = settings['retis']['swapfreq']
swap = False if True in prio_skip else rgen.rand() < swap_freq
if swap:
# Do RETIS moves
logger.info('Performing RETIS swapping move(s).')
results = retis_moves(ensembles, rgen, settings, cycle)
else:
logger.info('Performing RETIS TIS move(s)')
results = make_tis(ensembles, rgen, settings, cycle)
return results
[docs]def retis_moves(ensembles, rgen, settings, cycle):
"""Perform RETIS moves on the given ensembles.
This function will perform RETIS moves on the given ensembles.
First we have two strategies based on
`settings['retis']['swapsimul']`:
1) If `settings['retis']['swapsimul']` is True we will perform
several swaps, either ``[0^-] <-> [0^+], [1^+] <-> [2^+], ...``
or ``[0^+] <-> [1^+], [2^+] <-> [3^+], ...``. Which one of these
two swap options we use is determined randomly and they have
equal probability.
2) If `settings['retis']['swapsimul']` is False we will just
perform one swap for randomly chosen ensembles, i.e. we pick a
random ensemble and try to swap with the ensemble to the right.
Here we may also perform null moves if the
`settings['retis']['nullmove']` specifies so.
Parameters
----------
ensembles : list of dicts of objects
Each contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
rgen : object like :py:class:`.RandomGenerator`
This is a random generator. Here we assume that we can call
`rgen.rand()` to draw random uniform numbers.
settings : dict
This dict contains the settings for the RETIS method.
cycle : integer
The current cycle number.
Yields
------
out : dict
This dictionary contains the result of the RETIS moves.
"""
if settings['retis']['swapsimul']:
# Here we have two schemes:
# 1) scheme == 0: [0^-] <-> [0^+], [1^+] <-> [2^+], ...
# 2) scheme == 1: [0^+] <-> [1^+], [2^+] <-> [3^+], ...
if len(ensembles) < 3:
# Small number of ensembles, can only do the [0^-] <-> [0^+] swap:
scheme = 0
else:
scheme = 0 if rgen.rand() < 0.5 else 1
for i in range(scheme, len(ensembles) - 1, 2):
idx = ensembles[i]['path_ensemble'].ensemble_number
accept, trial, status = retis_swap_wrapper(
ensembles, idx, settings, cycle
)
result = {
'ensemble_number': idx,
'mc-move': 'swap',
'status': status,
'trial': trial[0],
'accept': accept,
'swap-with': idx + 1,
}
yield result
result = {
'ensemble_number': idx + 1,
'mc-move': 'swap',
'status': status,
'trial': trial[1],
'accept': accept,
'swap-with': idx,
}
yield result
# We might have missed some ensembles in the two schemes.
# Here, we do null moves in these, if requested:
if settings['retis']['nullmoves']:
if len(ensembles) % 2 != scheme: # Missed last ensemble:
# This is perhaps strange but it is equivalent to:
# (scheme == 0 and len(ensembles) % 2 != 0) or
# (scheme == 1 and len(ensembles) % 2 == 0)
accept, trial, status = null_move(ensembles[-1], cycle)
result = {
'ensemble_number':
ensembles[-1]['path_ensemble'].ensemble_number,
'mc-move': 'nullmove',
'status': status,
'trial': trial,
'accept': accept,
}
yield result
# We always miss the first ensemble in scheme 1:
if scheme == 1:
accept, trial, status = null_move(ensembles[0], cycle)
result = {
'ensemble_number':
ensembles[0]['path_ensemble'].ensemble_number,
'mc-move': 'nullmove',
'status': status,
'trial': trial,
'accept': accept,
}
yield result
else: # Just swap two ensembles:
idx = rgen.random_integers(0, len(ensembles) - 2)
accept, trial, status = retis_swap_wrapper(
ensembles, idx, settings, cycle
)
result = {
'ensemble_number': idx,
'mc-move': 'swap',
'status': status,
'trial': trial[0],
'accept': accept,
'swap-with': idx + 1,
}
yield result
result = {
'ensemble_number': idx + 1,
'mc-move': 'swap',
'status': status,
'trial': trial[1],
'accept': accept,
'swap-with': idx,
}
yield result
# Do null moves in the other ensembles, if requested:
if settings['retis']['nullmoves']:
for ensemble in ensembles:
idx2 = ensemble['path_ensemble'].ensemble_number
if idx2 not in (idx, idx + 1):
accept, trial, status = null_move(ensemble, cycle)
result = {
'ensemble_number': idx2,
'mc-move': 'nullmove',
'status': status,
'trial': trial,
'accept': accept,
}
yield result
def retis_swap_wrapper(ensembles, idx, settings, cycle):
"""Swap the last accepted paths of two TIS or two PPTIS ensembles.
This function selects the correct swap function to use, based on the
simulation task (retis or repptis).
Parameters
----------
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the RE(PP)TIS method.
Each one contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
idx : integer
Definition of what path ensembles to swap. We will swap
`ensembles[idx]` with `ensembles[idx+1]`.
settings : dict
This dict contains the settings for the RE(PP)TIS method.
cycle : integer
Current cycle number.
Returns
-------
out[0] : boolean
Should the path be accepted or not?
out[1] : list of object like :py:class:`.PathBase`
The trial paths.
out[2] : string
The status for the trial paths.
"""
if settings['simulation']['task'] == "repptis":
return repptis_swap(ensembles, idx, settings, cycle)
return retis_swap(ensembles, idx, settings, cycle)
[docs]def retis_swap(ensembles, idx, settings, cycle):
"""Perform a RETIS swapping move for two ensembles.
These ensembles are not PPTIS ensembles.
The RETIS swapping move will attempt to swap accepted paths between
two ensembles in the hope that path from [i^+] is an acceptable path
for [(i+1)^+] as well. We have two cases:
1) If we try to swap between [0^-] and [0^+] we need to integrate
the equations of motion.
2) Otherwise, we can just swap and accept if the path from [i^+] is
an acceptable path for [(i+1)^+]. The path from [(i+1)^+] is
always acceptable for [i^+] (by construction).
Parameters
----------
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the RETIS method.
Each one contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
idx : integer
Definition of what path ensembles to swap. We will swap
`ensembles[idx]` with `ensembles[idx+1]`. If `idx == 0` we have
case 1) defined above.
settings : dict
This dict contains the settings for the RETIS method.
cycle : integer
Current cycle number.
Returns
-------
out[0] : boolean
Should the path be accepted or not?
out[1] : list of object like :py:class:`.PathBase`
The trial paths.
out[2] : string
The status for the trial paths.
Note
----
Note that path.generated is **NOT** updated here. This is because
we are just swapping references and not the paths. In case the
swap is rejected updating this would invalidate the last accepted
path.
"""
logger.info("Swapping: %s <-> %s",
ensembles[idx]['path_ensemble'].ensemble_name,
ensembles[idx+1]['path_ensemble'].ensemble_name)
if idx == 0:
return retis_swap_zero(ensembles, settings, cycle)
path_ensemble1 = ensembles[idx]['path_ensemble']
path_ensemble2 = ensembles[idx+1]['path_ensemble']
path1 = path_ensemble1.last_path
path2 = path_ensemble2.last_path
ens_moves = [settings['ensemble'][i]['tis'].get('shooting_move', 'sh')
for i in [idx, idx+1]]
intf_w = [list(i) for i in (path_ensemble1.interfaces,
path_ensemble2.interfaces)]
for i, j in enumerate([settings['ensemble'][k] for k in (idx, idx+1)]):
if ens_moves[i] == 'wf':
intf_w[i][2] = j['tis'].get('interface_cap', intf_w[i][2])
# Check if path1 can be accepted in ensemble 2:
cross = path1.check_interfaces(path_ensemble2.interfaces)[-1]
accept = cross[1]
status = 'ACC' if accept else 'NCR'
if accept and settings['tis'].get('high_accept', False):
if 'ss' in ens_moves or 'wf' in ens_moves:
accept, status = high_acc_swap([path1, path2],
ensembles[idx]['rgen'],
intf_w[0], intf_w[1],
ens_moves)
if accept: # Accept the swap:
logger.info('Swap was accepted.')
# To avoid overwriting files, we move the paths to the
# generate directory here.
path_ensemble2.move_path_to_generate(path1)
path_ensemble1.move_path_to_generate(path2)
for i, path, path_ensemble, flag in ((0, path2, path_ensemble1, 's+'),
(1, path1, path_ensemble2, 's-')):
path.status = status
ens_set = settings['ensemble'][idx + i]
move = ens_moves[i]
path.weight = compute_weight(path, intf_w[i], move)\
if (ens_set['tis'].get('high_accept', False) and
move in ('wf', 'ss')) else 1
path.set_move(flag) if path.get_move() != 'ld' else \
path.set_move('ld')
# Then moved into accepted directory by the `add_path_data` below.
path_ensemble.add_path_data(path, status, cycle=cycle)
if cycle % ens_set.get('output', {}).get('restart-file', 1) == 0:
write_ensemble_restart(ensembles[idx+i], ens_set)
return accept, (path2, path1), status
logger.info('Swap was rejected. (%s)', status)
# Make shallow copies:
trial1 = copy.copy(path2)
trial2 = copy.copy(path1)
trial1.set_move('s+') # Came from right.
trial2.set_move('s-') # Came from left.
# Calculate weights:
for i, trial in ((0, trial1), (1, trial2)):
trial.weight = compute_weight(trial, intf_w[i], ens_moves[i])\
if (settings['ensemble'][idx+i]['tis'].get('high_accept', False)
and ens_moves[i] in ('wf', 'ss')) else 1.
path_ensemble1.add_path_data(trial1, status, cycle=cycle)
path_ensemble2.add_path_data(trial2, status, cycle=cycle)
return accept, (trial1, trial2), status
[docs]def retis_swap_zero(ensembles, settings, cycle):
"""Perform the RETIS swapping for ``[0^-] <-> [0^+]`` swaps.
The RETIS swapping move for ensembles [0^-] and [0^+] requires some
extra integration. Here we are generating new paths for [0^-] and
[0^+] in the following way:
1) For [0^-] we take the initial point in [0^+] and integrate
backward in time. This is merged with the second point in [0^+]
to give the final path. The initial point in [0^+] starts to the
left of the interface and the second point is on the right
side - i.e. the path will cross the interface at the end points.
If we let the last point in [0^+] be called ``A_0`` and the
second last point ``B``, and we let ``A_1, A_2, ...`` be the
points on the backward trajectory generated from ``A_0`` then
the final path will be made up of the points
``[..., A_2, A_1, A_0, B]``. Here, ``B`` will be on the right
side of the interface and the first point of the path will also
be on the right side.
2) For [0^+] we take the last point of [0^-] and use that as an
initial point to generate a new trajectory for [0^+] by
integration forward in time. We also include the second last
point of the [0^-] trajectory which is on the left side of the
interface. We let the second last point be ``B`` (this is on the
left side of the interface), the last point ``A_0`` and the
points generated from ``A_0`` we denote by ``A_1, A_2, ...``.
Then the resulting path will be ``[B, A_0, A_1, A_2, ...]``.
Here, ``B`` will be on the left side of the interface and the
last point of the path will also be on the left side of the
interface.
Parameters
----------
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the RETIS method.
It contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list.
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
settings : dict
This dict contains the settings for the RETIS method.
cycle : integer
The current cycle number.
Returns
-------
out : string
The result of the swapping move.
"""
path_ensemble0 = ensembles[0]['path_ensemble']
path_ensemble1 = ensembles[1]['path_ensemble']
engine0, engine1 = ensembles[0]['engine'], ensembles[1]['engine']
maxlen0 = settings['ensemble'][0]['tis']['maxlength']
maxlen1 = settings['ensemble'][1]['tis']['maxlength']
ens_moves = [settings['ensemble'][i]['tis'].get('shooting_move', 'sh')
for i in [0, 1]]
intf_w = [list(i) for i in (path_ensemble0.interfaces,
path_ensemble1.interfaces)]
for i, j in enumerate([settings['ensemble'][k] for k in (0, 1)]):
if ens_moves[i] == 'wf':
intf_w[i][2] = j['tis'].get('interface_cap', intf_w[i][2])
# 0. check if MD is allowed
allowed = (path_ensemble0.last_path.get_end_point(
path_ensemble0.interfaces[0],
path_ensemble0.interfaces[-1]) == 'R')
# 1. Generate path for [0^-] from [0^+]:
# We generate from the first point of the path in [0^+]:
logger.debug('Swapping [0^-] <-> [0^+]')
logger.debug('Creating path for [0^-]')
system = path_ensemble1.last_path.phasepoints[0].copy()
logger.debug('Initial point is: %s', system)
ensembles[0]['system'] = system
# Propagate it backward in time:
path_tmp = path_ensemble1.last_path.empty_path(maxlen=maxlen1-1)
if allowed:
logger.debug('Propagating for [0^-]')
engine0.propagate(path_tmp, ensembles[0], reverse=True)
else:
logger.debug('Not propagating for [0^-]')
path_tmp.append(system)
path0 = path_tmp.empty_path(maxlen=maxlen0)
for phasepoint in reversed(path_tmp.phasepoints):
path0.append(phasepoint)
# Add second point from [0^+] at the end:
logger.debug('Adding second point from [0^+]:')
# Here we make a copy of the phase point, as we will update
# the configuration and append it to the new path:
phase_point = path_ensemble1.last_path.phasepoints[1].copy()
logger.debug('Point is %s', phase_point)
engine1.dump_phasepoint(phase_point, 'second')
path0.append(phase_point)
if path0.length == maxlen0:
path0.status = 'BTX'
elif path0.length < 3:
path0.status = 'BTS'
elif ('L' not in set(path_ensemble0.start_condition) and
'L' in path0.check_interfaces(path_ensemble0.interfaces)[:2]):
path0.status = '0-L'
else:
path0.status = 'ACC'
# 2. Generate path for [0^+] from [0^-]:
logger.debug('Creating path for [0^+] from [0^-]')
# This path will be generated starting from the LAST point of [0^-] which
# should be on the right side of the interface. We will also add the
# SECOND LAST point from [0^-] which should be on the left side of the
# interface, this is added after we have generated the path and we
# save space for this point by letting maxlen = maxlen1-1 here:
path_tmp = path0.empty_path(maxlen=maxlen1-1)
# We start the generation from the LAST point:
# Again, the copy below is not needed as the propagate
# method will not alter the initial state.
system = path_ensemble0.last_path.phasepoints[-1].copy()
if allowed:
logger.debug('Initial point is %s', system)
ensembles[1]['system'] = system
logger.debug('Propagating for [0^+]')
engine1.propagate(path_tmp, ensembles[1], reverse=False)
# Ok, now we need to just add the SECOND LAST point from [0^-] as
# the first point for the path:
path1 = path_tmp.empty_path(maxlen=maxlen1)
phase_point = path_ensemble0.last_path.phasepoints[-2].copy()
logger.debug('Add second last point: %s', phase_point)
engine0.dump_phasepoint(phase_point, 'second_last')
path1.append(phase_point)
path1 += path_tmp # Add rest of the path.
else:
path1 = path_tmp
path1.append(system)
logger.debug('Skipping propagating for [0^+] from L')
if path_ensemble1.last_path.get_move() != 'ld':
path0.set_move('s+')
else:
path0.set_move('ld')
if path_ensemble0.last_path.get_move() != 'ld':
path1.set_move('s-')
else:
path1.set_move('ld')
if path1.length >= maxlen1:
path1.status = 'FTX'
elif path1.length < 3:
path1.status = 'FTS'
else:
path1.status = 'ACC'
logger.debug('Done with swap zero!')
# Final checks:
accept = path0.status == 'ACC' and path1.status == 'ACC'
status = 'ACC' if accept else (path0.status if path0.status != 'ACC' else
path1.status)
# High Acceptance swap is required when Wire Fencing are used
if accept and settings['tis'].get('high_accept', False):
if 'wf' in ens_moves:
accept, status = high_acc_swap([path1, path_ensemble1.last_path],
ensembles[0]['rgen'],
intf_w[0],
intf_w[1],
ens_moves)
for i, path, path_ensemble, flag in ((0, path0, path_ensemble0, 's+'),
(1, path1, path_ensemble1, 's-')):
if not accept and path.status == 'ACC':
path.status = status
# These should be 1 unless length of paths equals 3.
# This technicality is not yet fixed. (An issue is open as a reminder)
ens_set = settings['ensemble'][i]
move = ens_moves[i]
path.weight = compute_weight(path, intf_w[i], move)\
if (ens_set['tis'].get('high_accept', False) and
move in ('wf', 'ss')) else 1
if accept:
path_ensemble.move_path_to_generate(path)
else:
logger.debug("Rejected swap path in [0^%s], %s", flag[:-1], status)
path_ensemble.add_path_data(path, path.status, cycle=cycle)
if cycle % ens_set.get('output', {}).get('restart-file', 1) == 0:
write_ensemble_restart(ensembles[i], settings['ensemble'][i])
return accept, (path0, path1), status
def repptis_swap(ensembles, idx, settings, cycle, inc_sh=True):
"""Perform a replica exchange move between adjacent PPTIS ensembles.
First, propagation directions are randomly chosen for both paths, and it is
checked whether the pathtypes (i.e. LMR, RMR, ...) are compatible. If not,
the move is rejected. If the proposed directions are compatible with the
pathtypes, the path segments in the overlapping region are cut, and they
are propagated in the proposed directions.
When compatible, the paths are created as follows:
1) new [i^+-] path from old [(i+1)^+-] path. The path segment in the 'LM'
region is cut out of the old [(i+1)^+-] path (with one point left of the
L interface, and one point right of the M interface). This path segment is
then extended by propagating the phasepoint left of the L interface in the
proposed direction till a crossing condition of the [i^+-] interface occurs
2) new [(i+1)^+-] path from old [i^+-] path. The path segment in the 'MR'
region is cut out of the old [i^+-] path (with one point left of the M
interface, and one point right of the R interface). This path segment is
then extended by propagating the phasepoint right of the R interface in the
proposed direction until a crossing condition of the [(i+1)^+-] interface
occurs.
Floating point precision problem (for GRO):
The phase point to be shot from, is first dumped to a .gro file, where
precision goes from ~9 decimals to 3 decimals. This impacts the order
parameter of the dumped phasepoint, and if it was very closely located to
an interface, it might even switch sides w.r.t. that interface. The problem
has an easy fix. We include the shooting point in the overlap_path, such
that the order parameter does not change. When we shoot from this point,
the order parameter will still change due to precision decrease, and the
first point of the propagated_path will be slightly different. We discard
this point in the propagated_path, and keep the one from the overlap_path.
Whichever shootpoint you keep, there will be a mismatch. The choice we make
does not cause trouble with the swapping move. The above described behavior
(i.e. the fix) is enabled with the paramter inc_sh set to True, being
the standard behavior.
Parameters
----------
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the RETIS method.
It contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list.
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
idx : int
The index of the ensemble that will be swapped.
settings : dict
This dict contains the settings for the RETIS method.
cycle : integer
The current cycle number.
inc_sh : bool
If True, the phasepoint that is used for shooting is included in the
overlap path. If False, it will be included in the newly propagated
path (with lower .gro precision).
Returns
-------
out[0] : boolean
Should the paths be accepted or not?
out[1] : list of object like :py:class:`.PathBase`
The trial paths.
out[2] : string
The status for the trial paths.
"""
msg = f"Swapping: {ensembles[idx]['path_ensemble'].ensemble_name} <-> "
msg += f"{ensembles[idx+1]['path_ensemble'].ensemble_name}"
logger.info(msg)
# if idx == 0, we perform a swap_zero move
if idx == 0:
# Since we allow RML paths in [0+-'], we need to do an extra check,
# as this pathtype will result in a bogus swap...
# NB: if we don't check for this, the move will be rejected anyways,
# as the first propagated ph will be OOB for [0-], which results in
# an "FTS" or "BTS" status. However, it's cleaner to reject the move
# in advance...
start, end, _, _ =\
ensembles[1]['path_ensemble'].last_path.check_interfaces(
ensembles[1]['path_ensemble'].interfaces)
ptype = ""+start+"M"+end
if ptype == "RML":
msg = "Rejecting swap: RML in [0+-'] cannot be swapped with [0-]!"
logger.info(msg)
return reject_repptis_swap_wrong_prop(
ensembles[0]['path_ensemble'], ensembles[1]['path_ensemble'],
cycle, settings, idx, ensembles)
return retis_swap_zero(ensembles, settings, cycle)
# else we perform a repptis swap move
path_ensemble0 = ensembles[idx]['path_ensemble']
path_ensemble1 = ensembles[idx+1]['path_ensemble']
# Choose propagation directions for the two paths, and check if the
# proposed move is allowed or not
allowed, propdir0, propdir1, path0_type, path1_type = \
re_propagation_directions(ensembles, idx)
# If not allowed, we immediately reject the swap
if not allowed:
return reject_repptis_swap_wrong_prop(path_ensemble0, path_ensemble1,
cycle, settings, idx, ensembles)
# If allowed, the swap move is attempted
logger.info("Attempting swap.")
logger.debug("Attempting swap of %s and %s paths with propdirs %s and %s.",
path0_type, path1_type, propdir0, propdir1)
# First create path0_new, then path1_new.
# We will immediately reject the move if path0_new is not valid.
# -----------------------------------------------------------------------
# Create new path for [i+-]
maxlen0 = settings['ensemble'][idx]['tis']['maxlength']
engine0 = ensembles[idx]['engine']
# We distinguish between the part of the path that is copied (overlap,cut),
# and the part that is propagated (prop).
# First we cut:
maxlen_cut0 = maxlen0 if inc_sh else maxlen0 - 1
path_overlap0_new = path_ensemble0.last_path.empty_path(maxlen=maxlen_cut0)
path0_new_shootpoint = \
cut_overlap_phasepoints(path_ensemble1, propdir1, path_overlap0_new,
side='right', include_shootpoint=inc_sh)
# Then create empty path to be filled with propagated points:
# + 1 because the shootpoint will not be used in the propagated path.
# (Thus, one point of prop_path will be discarded)
maxlen_prop0_new = maxlen0 - path_overlap0_new.length + 1
path_prop0_new = path_ensemble0.last_path.empty_path(
maxlen=maxlen_prop0_new)
# Set the system to the shootpoint:
ensembles[idx]['system'] = path0_new_shootpoint
logger.debug("Initial point is %s", ensembles[idx]['system'])
# Propagate the path:
engine0.propagate(path_prop0_new, ensembles[idx],
reverse=propdir1 == 'backwards')
if inc_sh: # We need to remove the shootpoint from the propagated path
path_prop0_new.phasepoints.pop(0)
# Append the propagated path to the overlap path:
if propdir1 == 'backwards':
path0_new = paste_paths(path_prop0_new, path_overlap0_new,
overlap=False, maxlen=maxlen0)
elif propdir1 == 'forwards':
path0_new = paste_paths(path_overlap0_new, path_prop0_new,
overlap=False, maxlen=maxlen0)
# Check if the new path is valid:
path0_new.status = set_swap_status(path0_new, propdir1, maxlen0)
# If the new path is not valid, we stop and reject the entire swap move
if path0_new.status != 'ACC':
return reject_repptis_swap_half(path_ensemble0, path_ensemble1,
path0_new, cycle, settings, idx,
ensembles)
# -----------------------------------------------------------------------
# Create new path for [(i+1)+-]
maxlen1 = settings['ensemble'][idx+1]['tis']['maxlength']
engine1 = ensembles[idx+1]['engine']
# First we cut:
maxlen_cut1 = maxlen1 if inc_sh else maxlen1 - 1
path_overlap1_new = path_ensemble1.last_path.empty_path(maxlen=maxlen_cut1)
path1_new_shootpoint = \
cut_overlap_phasepoints(path_ensemble0, propdir0, path_overlap1_new,
side='left', include_shootpoint=inc_sh)
# Then create empty path to be filled with propagated points:
# + 1 because the shootpoint will not be used in the propagated path.
# (Thus, one point of prop_path will be discarded)
maxlen_prop1_new = maxlen1 - path_overlap1_new.length + 1
path_prop1_new = path_ensemble1.last_path.empty_path(
maxlen=maxlen_prop1_new)
# Set the system to the shootpoint:
ensembles[idx+1]['system'] = path1_new_shootpoint
logger.debug("Initial point is %s", ensembles[idx+1]['system'])
# Propagate the path:
engine1.propagate(path_prop1_new, ensembles[idx+1],
reverse=propdir0 == 'backwards')
if inc_sh: # We need to remove the shootpoint from the propagated path
path_prop1_new.phasepoints.pop(0)
# Append the propagated path to the overlap path:
if propdir0 == 'backwards':
path1_new = paste_paths(path_prop1_new, path_overlap1_new,
overlap=False, maxlen=maxlen1)
elif propdir0 == 'forwards':
path1_new = paste_paths(path_overlap1_new, path_prop1_new,
overlap=False, maxlen=maxlen1)
# Check if the new path is valid:
path1_new.status = set_swap_status(path1_new, propdir0, maxlen1)
# -----------------------------------------------------------------------
# Bookkeeping
input_tuple = (path_ensemble0, path_ensemble1, path0_new, path1_new,
settings, idx, cycle, ensembles)
return repptis_swap_bookkeeping(input_tuple)
def repptis_swap_bookkeeping(input_tuple):
"""Bookkeeping for the repptis swap move.
If you reached here, then the first half of the swapping move was already
successful. Here it is checked whether the second half of the move is
successful or not. Afterwards the final bookkeeping is performed.
The input_tuple contains the parameters below, which are assigned at the
beginning of the function.
Parameters
----------
path_ensemble0 : object like :py:class:`.PathEnsemble`
The path ensemble of the first trial path ([i+-]).
path_ensemble1 : object like :py:class:`.PathEnsemble`
The path ensemble of the second trial path ([(i+1)+-]).
path0_new : object like :py:class:`.PathBase`
The first trial path ([i+-]).
path1_new : object like :py:class:`.PathBase`
The second trial path ([(i+1)+-]).
settings : dict
This dict contains the settings for the REPPTIS method.
idx : int
The index of the ensemble that will be swapped with the next one.
cycle : integer
The current cycle number.
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the REPPTIS method.
Returns
-------
out[0] : boolean
Should the paths be accepted or not?
out[1] : list of object like :py:class:`.PathBase`
The trial paths.
out[2] : string
The status for the trial paths.
"""
# Get the input:
path_ensemble0, path_ensemble1, path0_new, path1_new, settings, idx, \
cycle, ensembles = input_tuple
status = 'ACC'
assert path0_new.status == 'ACC', "First half of swap should've been ACC."
if path1_new.status != 'ACC': # This can happen (FTX or BTX)
status = path1_new.status
for path_new, path_ensemble, flag in ((path0_new, path_ensemble0, 's+'),
(path1_new, path_ensemble1, 's-')):
path_new.status = status
path_new.weight = 1 # NB: no high acceptance allowed for repptis
path_new.set_move('ld' if path_ensemble.last_path.get_move() == 'ld'
else flag)
accept = path0_new.status == 'ACC' and path1_new.status == 'ACC'
# Here is where you would want to put 'plot_repptis_swap' for debugging
# plot_repptis_swap(path_ensemble0, path_ensemble1, path0_new, path1_new,
# idx, cycle, accept)
if accept:
path0_new = path_ensemble0.copy_path_to_generate(path0_new)
path1_new = path_ensemble1.copy_path_to_generate(path1_new)
for i, path_new, path_ensemble in ((0, path0_new, path_ensemble0),
(1, path1_new, path_ensemble1)):
path_ensemble.add_path_data(path_new, status, cycle=cycle)
ens_set = settings['ensemble'][idx+i]
if cycle % ens_set.get('output', {}).get('restart-file', 1) == 0:
write_ensemble_restart(ensembles[idx+i], ens_set)
return accept, (path0_new, path1_new), status
def set_swap_status(path, propdir, maxlen):
"""Set the status of a path after a swap-propagation in REPPTIS.
Normally, the path extension should automatically result in an acceptable
path of the adjacent ensemble, but there are two exceptions:
1) The path is too long (exceeds maxlen): breaks detailed balance.
2) The path is too short (length < 3): Not a path breaks detailed balance?
Parameters
----------
path : object like :py:class:`.PathBase`
The path that was extended by a REPPTIS swap move.
propdir : string
The propagation direction of the path. Either 'forwards' or 'backwards'
maxlen : int
The maximum length of the path.
Returns
-------
status : string
The status of the path after the REPPTIS swap move.
"""
# First we check if it's too long
status_map = {'backwards': {True: 'BTX', False: 'ACC'},
'forwards': {True: 'FTX', False: 'ACC'}}
status = status_map[propdir][path.length >= maxlen]
# Then we check if it's too short. This is theoretically not possible,
# as we already start from an overlap path that is by definition 2
# phasepoints long (and we add another one with shooting). However, we
# keep the check, just in case something terrible happened.
assert path.length >= 3, "Pathlength < 3 is NOT possible in REPPTIS swap."
return status
def reject_repptis_swap_wrong_prop(path_ensemble0, path_ensemble1, cycle,
settings, idx, ensembles):
"""Reject REPPTIS swap for incompatible propagation directions.
E.g. Trying to extend an LMR path of [i+-] into a [(i+1)+-] path by
backwards propagation.
Parameters
----------
path_ensemble0 : object like :py:class:`.PathEnsemble`
The path ensemble of the first trial path ([i+-]).
path_ensemble1 : object like :py:class:`.PathEnsemble`
The path ensemble of the second trial path ([(i+1)+-]).
cycle : integer
The current cycle number.
settings : dict
This dict contains the settings for the REPPTIS method.
idx : int
The index of the ensemble that will be swapped with the next one.
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the REPPTIS method.
Returns
-------
out[0] : boolean
Should the paths be accepted or not?
out[1] : list of object like :py:class:`.PathBase`
The trial paths.
out[2] : string
The status for the trial paths.
"""
logger.info("Swap rejected: incompatible propagation directions")
status = 'SWD' # swap wrong direction
accept = False
# Make shallow copies:
trial0 = copy.copy(path_ensemble1.last_path)
trial1 = copy.copy(path_ensemble0.last_path)
for trial, move in ((trial0, 's+'), (trial1, 's-')):
trial.set_move('ld' if trial.get_move() == 'ld' else move)
trial.weight = 1.0
path_ensemble0.add_path_data(trial0, status, cycle=cycle)
path_ensemble1.add_path_data(trial1, status, cycle=cycle)
for i in (0, 1):
ens_set = settings['ensemble'][idx+i]
if cycle % ens_set.get('output', {}).get('restart-file', 1) == 0:
write_ensemble_restart(ensembles[idx+i], ens_set)
return accept, (trial0, trial1), status
def reject_repptis_swap_half(path_ensemble0, path_ensemble1, path0_new, cycle,
settings, idx, ensembles):
"""Bookkeeping when the first REPPTIS extension is not acceptable.
If the first extension of the REPPTIS swap move does not result in an
acceptable path for the [(i+1)+-] ensemble, we do not waste computational
time on the second extension. We immediately reject the move.
Parameters
----------
path_ensemble0 : object like :py:class:`.PathEnsemble`
The path ensemble of the first trial path ([i+-]).
path_ensemble1 : object like :py:class:`.PathEnsemble`
The path ensemble of the second trial path ([(i+1)+-]).
path0_new : object like :py:class:`.PathBase`
The first trial path ([i+-]).
cycle : integer
The current cycle number.
settings : dict
This dict contains the settings for the REPPTIS method.
idx : int
The index of the ensemble that will be swapped with the next one.
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the REPPTIS method.
Returns
-------
out[0] : boolean
Should the paths be accepted or not?
out[1] : list of object like :py:class:`.PathBase`
The trial paths.
out[2] : string
The status for the trial paths.
"""
logger.info("Swap rejected: first path not accepted")
status = 'SWH' # swap rejected
accept = False
# Make shallow copies:
trial0 = copy.copy(path0_new)
trial1 = copy.copy(path_ensemble0.last_path)
for trial, move in ((trial0, 's+'), (trial1, 's-')):
trial.set_move('ld' if trial.get_move() == 'ld' else move)
trial.weight = 1.0
path_ensemble0.add_path_data(trial0, status, cycle=cycle)
path_ensemble1.add_path_data(trial1, status, cycle=cycle)
for i in (0, 1):
ens_set = settings['ensemble'][idx+i]
if cycle % ens_set.get('output', {}).get('restart-file', 1) == 0:
write_ensemble_restart(ensembles[idx+i], ens_set)
return accept, (trial0, trial1), status
def cut_overlap_phasepoints(path_ensemble, propdir, overlap_path, side,
include_shootpoint=True):
"""Cut path in overlapping region of adjacent ensembles.
Cuts the phasepoints in the overlapping region of two adjacent PPTIS
ensembles [i^+] and [(i+1)^-]. The overlapping region is [L,M] if the
we are cutting from an [(i+1)^+-] path (denoted by side = 'right'), and
[M,R] if we are cutting from an [i^+-] path (denoted by side = 'left).
The propagation direction is the direction in which the cut path is to be
propagated.
Parameters
----------
path_ensemble: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
propdir : string
This is either 'forwards' or 'backwards', and denotes the direction
in which the cut path is to be propagated.
overlap_path : object like :py:class:`.Path`
The path that is to be cut.
side : string
This is either 'left' or 'right', and denotes the side from which the
path originates. left: [i^+-] path to be propagated into [(i+1)^+-]
right: [(i+1)^+-] path to be propagated into [i^+-]
include_shootpoint : bool
If True, the shootpoint of the path is included in the cut path.
Returns
-------
shoot_point : object like :py:class:`.PhasePoint`
The phasepoint from which the cut path is to be propagated.
"""
# intf_M defines when to stop cutting phasepoints
# For now, this is equivalent for all PPTIS ensembles, but might change
# in the future.
intf_m = path_ensemble.interfaces[0] if path_ensemble.ensemble_number == 1\
else path_ensemble.interfaces[1]
stop_idx = 1 # to not include shootpoint in the cut path loop
# Different behavior depending on the side from wich the path originates,
# and the propagation direction. Note that, for the paths. Note that, if
# the propagation direction is forwards, we reverse the order in which we
# cut the phasepoints. We do not re-reverse this in the end, because the
# function 'paste_paths' will be used. As this cut-out path is the backward
# part of the merged paths, it will be reversed in that function.
# For the backwards propagation, the cut-out path is the forwards part of
# the merged paths, and will not be reversed in 'paste_paths' (or here).
if side == 'right':
if propdir == 'backwards':
shoot_point = path_ensemble.last_path.phasepoints[0].copy()
if include_shootpoint:
overlap_path.append(shoot_point.copy())
phase_points = path_ensemble.last_path.phasepoints[stop_idx:]
for phase_point in phase_points:
overlap_path.append(phase_point.copy())
if phase_point.order[0] >= intf_m: # if>=M
break
elif propdir == 'forwards':
shoot_point = path_ensemble.last_path.phasepoints[-1].copy()
if include_shootpoint:
overlap_path.append(shoot_point.copy())
phase_points = path_ensemble.last_path.phasepoints[:-stop_idx]
for phase_point in reversed(phase_points):
overlap_path.append(phase_point.copy())
if phase_point.order[0] >= intf_m: # if>=M
break
elif side == 'left':
if propdir == 'backwards':
shoot_point = path_ensemble.last_path.phasepoints[0].copy()
if include_shootpoint:
overlap_path.append(shoot_point.copy())
phase_points = path_ensemble.last_path.phasepoints[stop_idx:]
for phase_point in phase_points:
overlap_path.append(phase_point.copy())
if phase_point.order[0] <= intf_m: # if<=M
break
elif propdir == 'forwards':
shoot_point = path_ensemble.last_path.phasepoints[-1].copy()
if include_shootpoint:
overlap_path.append(shoot_point.copy())
phase_points = path_ensemble.last_path.phasepoints[:-stop_idx]
for phase_point in reversed(phase_points):
overlap_path.append(phase_point.copy())
if phase_point.order[0] <= intf_m: # if<=M
break
return shoot_point
def fix_directions(ensembles, idx): # pragma: no cover
"""Propose acceptable propagation direcs for [i^+-] and [(i+1)+-] paths.
If no acceptable directions are found, a rejection is returned.
Note that this will **break detailed balance**.
Parameters
----------
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the RETIS method.
It contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list.
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
idx : integer
Index of the ensemble, which will be swapped with the (idx+1) ensemble.
Returns
-------
accepted : bool
Whether the proposed directions can result in an acceptable swap.
propdir0 : string
The proposed propagation direction for the [i^+] path. This is either
'forwards' or 'backwards', or None if no acceptable direction was found
propdir1 : string
The proposed propagation direction for the [(i+1)^-] path. This is
either 'forwards' or 'backwards', or None if no acceptable direction
was found.
path0_type : string
The type of the [i^+] path. This is either 'LMR', 'RMR', 'LML' or 'RML'
path1_type : string
The type of the [(i+1)^-] path. This is either 'LMR', 'RMR', 'LML' or
'RML'.
Notes
-----
This function does **not** obey detailed balance (DB). Used for debugging.
It could allow detailed balance, but then the path ensemble definitions
should be changed (weight dependent on path type!)
"""
path_ensemble0 = ensembles[idx]['path_ensemble']
path_ensemble1 = ensembles[idx+1]['path_ensemble']
# Get the start and end points of the paths (L or R)
path0_old_start, path0_old_end, _, _ = \
path_ensemble0.last_path.check_interfaces(path_ensemble0.interfaces)
path1_old_start, path1_old_end, _, _ = \
path_ensemble1.last_path.check_interfaces(path_ensemble1.interfaces)
# Define the path types, for nice output
path0_type = ""+path0_old_start+"M"+path0_old_end
path1_type = ""+path1_old_start+"M"+path1_old_end
# Define the allowed propagation directions
propdirs = {'leftright': {'LMR': 'forwards',
'RMR': 'forwards', # backwards could work too
'LML': None,
'RML': 'backwards'},
'rightleft': {'LMR': 'backwards',
'RMR': None,
'LML': 'backwards', # forwards could work too
'RML': 'forwards'}}
# allocate directions
propdir0 = propdirs['leftright'][path0_type]
propdir1 = propdirs['rightleft'][path1_type]
# Check whether it is acceptable
accepted = propdir0 in ('forwards', 'backwards') and \
propdir1 in ('forwards', 'backwards')
# Log the move
msg = "Proposed swap move:\n"
msg += f"{path_ensemble0.ensemble_name} <-->" + \
f"{path_ensemble1.ensemble_name}\n" + \
f"{path0_type} <-----> {path1_type}\n" + \
f"{propdir0}<>{propdir1}"
logger.info(msg)
return accepted, propdir0, propdir1, path0_type, path1_type
def re_propagation_directions(ensembles, idx):
"""Proposes propagation directions for the [i^+-] and [(i+1)^+-] paths.
And checks whether this results in an acceptable swapping move.
Parameters
----------
ensembles : list of dictionaries of objects
This is a list of the ensembles we are using in the RETIS method.
It contains:
* `path_ensemble`: object like :py:class:`.PathEnsemble`
This is used for storing results for the simulation. It
is also used for defining the interfaces for this simulation.
* `system`: object like :py:class:`.System`
System is used here since we need access to the temperature
and to the particle list.
* `order_function`: object like :py:class:`.OrderParameter`
The class used for calculating the order parameter(s).
* `engine`: object like :py:class:`.EngineBase`
The engine to use for propagating a path.
idx : integer
Index of the ensemble, which will be swapped with the (idx+1) ensemble.
Returns
-------
allowed : bool
Whether the proposed directions can result in an acceptable swap.
propdir0 : string
The proposed propagation direction for the [i^+-] path. This is either
'forwards' or 'backwards', or None if no acceptable direction was found
propdir1 : string
The proposed propagation direction for the [(i+1)^+-] path. This is
either 'forwards' or 'backwards', or None if no acceptable direction
was found.
path0_type : string
The type of the [i^+-] path. This is 'LMR', 'RMR', 'LML' or 'RML'.
path1_type : string
The type of the [(i+1)^+-] path. This is either 'LMR', 'RMR', 'LML' or
'RML'.
"""
# Set path ensembles
path_ensemble0 = ensembles[idx]['path_ensemble']
path_ensemble1 = ensembles[idx+1]['path_ensemble']
# Check the start and end points of the paths (L or R)
path0_old_start, path0_old_end, _, _ = \
path_ensemble0.last_path.check_interfaces(path_ensemble0.interfaces)
path1_old_start, path1_old_end, _, _ = \
path_ensemble1.last_path.check_interfaces(path_ensemble1.interfaces)
# Define the path types, for nice output
path0_type = ""+path0_old_start+"M"+path0_old_end
path1_type = ""+path1_old_start+"M"+path1_old_end
# Pick two random numbers rand0 and rand1 between [0,1].
# If rand0 < 0.5: attempt backward propagation for path0. Otherwise forward
# If rand1 < 0.5: attempt backward propagation for path1. Otherwise forward
# For now, we will just pick a random number for each path from the
# random number generator of the 0 ensemble, because there are seed issues.
idx1 = int(round((ensembles[0]['rgen'].rand())[0]) + 0.1)
idx2 = int(round((ensembles[0]['rgen'].rand())[0]) + 0.1)
propdir0 = ['backwards', 'forwards'][idx1]
propdir1 = ['backwards', 'forwards'][idx2]
# Check if the proposed move is allowed
allow = {'leftright': {'L': False, 'R': True},
'rightleft': {'L': True, 'R': False}}
allowed0 = allow['leftright'][path0_old_start if propdir0 == 'backwards'
else path0_old_end]
allowed1 = allow['rightleft'][path1_old_start if propdir1 == 'backwards'
else path1_old_end]
# Both paths must be allowed to begin the swapping move
allowed = allowed0 and allowed1
# Log the move
msg = "Proposed swap move:\n"
msg += f"{path_ensemble0.ensemble_name} <-->" + \
f"{path_ensemble1.ensemble_name}\n" + \
f"{path0_type} <-----> {path1_type}\n" + \
f"{propdir0}<>{propdir1}\n" + \
f"Allowed: {allowed}"
logger.info(msg)
return allowed, propdir0, propdir1, path0_type, path1_type
def high_acc_swap(paths, rgen, intf0, intf1, ens_moves):
"""Accept or Reject a swap move using the High Acceptance weights.
Parameters
----------
paths: list of object like :py:class:`.PathBase`
The path in the LOWER and UPPER ensemble to exchange.
rgen : object like :py:class:`.RandomGenerator`
This is a random generator.
intf0: list of float
The interfaces of the LOWER ensemble.
intf1: list of float
The interfaces of the HIGHER ensemble.
ens_moves: list of string
The moves used in the two ensembles.
Returns
-------
out[0] : boolean
True if the move should be accepted.
Notes
-----
- This function is needed only when paths generated via Wire Fencing or
Stone Skipping are involved.
- In the case that a path bears a flag 'ld', the swap is accepted,
but the flag will be unchanged.
"""
# Crossing before the move
c1_old = compute_weight(paths[0], intf0, ens_moves[0])
c2_old = compute_weight(paths[1], intf1, ens_moves[1])
# Crossing if the move would be accepted
c1_new = compute_weight(paths[1], intf0, ens_moves[0])
c2_new = compute_weight(paths[0], intf1, ens_moves[1])
if c1_old == 0 or c2_old == 0:
logger.warning("div_by_zero. c1_old, c2_old, ens_moves: [%i,%i], %s",
c1_old, c2_old, str(ens_moves))
p_swap_acc = 1
else:
p_swap_acc = c1_new*c2_new/(c1_old*c2_old)
# Finally, randomly decide to accept or not:
if rgen.rand() < p_swap_acc:
return True, 'ACC' # Accepted
return False, 'HAS' # Rejected