The Weighted Ensemble Dynamics algorithm (Huber and Kim, Biophys. J, 1996) is useful for finding rates of rare events. It maintains several copies of the system of interest, and divides and merges copies to enhance sampling over free energy barriers along the reaction coordinate. It makes use of divisions of configuration space usually called bins, and approximately enforces the condition of having equal number of system copies in each bin.
This code can run several threads
on any machine that has pthreads,
which is standard on most architectures.
Upon the first call to step_objects_forward,
the threads are created and used to independently
call step_forward on the objects. The threads
are kept waiting between calls to step_objects_forward,
and are sychronized so that the main thread does not
return until all the threads are finished.
In addition to the simulation, this code can also generate sequential bins along a user-defined reaction coordinate. The user starts the system copies at one end of the coordinate, and the bins are generated by stepping the copies forward and forcing the ensemble along the coordinate.
The algorithm is wrapped in a generic
class Sampler< Interface>, which sits in the namespace
Weighted_Ensemble_Dynamics::Multithreaded. (I still
have a single-threaded version which does not use pthreads at all;
the multithreaded version sits on top of it.) The header
file we_multithread.hh is included (which then pulls
in other header files in the same directory),
and the code
must also be linked to the compiled result of barrier.cc.
Both of these files are in the src directory.
The Interface class, provided by the
user, must have the following members:
typename Info - reference to an object that
contains all the information about the system to be simulated, including
the system copies and their weights. In most cases this will be pointer.
typename Object - reference to an individual object;
used in conjunction with the Info reference. Often this
will be a pointer, but one could also use an integer index, for example.
typename Exception - type of exception that is thrown
from inside
void initialize_object( Info info, Object obj) -
called on obj at the beginning of the simulation and
whenever it is restarted after crossing a reaction boundary.
void step_forward( Info info, unsigned int ithread, Object obj) - the configuration of obj is stepped forward by one time step. Thread
number ithread is the one doing the work. Threads are numbered starting with 0.
unsigned int bin( Info info, Object obj) -
returns the integer designating the bin containing obj
unsigned int bin( Info info, Object obj, bool& exists) -
returns the integer designating the bin containing obj.
Also outputs into exists whether a bin actually exists.
If exists is false after the function call,
the integer return value is undefined. This function is necessary because
no bins are defined yet during the bin-building process.
Object_Ref new_object( Info info) -
creates a new object and returns
its reference. It is the user's responsibility to assign an appropriate weight
to the object.
void get_duplicated_objects( Info info, Object obj, unsigned int ncopies,
std::list< Object>& objs)-
places ncopies copies of obj into the list objs,
which is initially empty.
Also, the weight originally in obj is split evenly among obj
and its copies. So, if obj initially has a weight of 1 and
ncopies is 3, then after the call obj and the three new
copies each have a weight of 0.25.
void remove( Objects objs, Object obj) - removes obj from the container.
void get_crossing( Info info, Object obj, bool& crossed, unsigned int& crossing) - We allow for more than one reaction condition, and each reaction
condition is designated by a non-negative integer, starting at zero.
This function tests the object obj to see if it has crossed.
(crossed is set true if so, false otherwise).
If there is a crossing, then crossing is set to the integer
for that crossing.
unsigned int number_of_bins( Info info)
unsigned int number_of_bin_neighbors( Info info, unsigned int bin) -
Each bin is adjacent to one or more other bins. For example, in
a conventional reaction coordinate approach, each bin would have
two neighbors, except for the ones at either end. We allow for more
general bin connectivity, so this function returns the number of
neighbors of bin.
template< class Function> void apply_to_bin_neighbors( Info info,
unsigned int bin, Function&) - applies an arbitary
function object Function to each reference of
of neighbors of bin. Function will have
the () operator overloaded to take an unsigned integer.
bool have_population_reversal( Info info, unsigned int bin0, double wt0,
unsigned int bin1, double wt1) - At steady state, one expects some
bins to have more system statistical weight than others. Due to
fluctuations in the simulation, this isn't always so, so the algorithm
will merge bins in order to get better convergence. This function
tests the weights (wt0 and wt1) of the
bins (bin0 and bin1) and returns true if
the relative weights aren't what one would expect in steady state.
This is optional, and one can always just return false
and the algorithm will still work, but maybe not as efficiently.
void set_null( Info& info) - sets reference
to a "null" value, which is used to denote an uninitialized
reference.
bool is_null( Info info) -
returns true if reference is
to a "null" value.
double uniform( Info info) - returns a random number uniformly
distributed between 0 and 1.
double weight( Info info, Obj obj) - weight of obj
void set_weight( Info info, double weight, Obj obj) - sets
weight of obj to weight.
unsigned int number_of_threads( Info info) - number of threads
used by system.
void output( Info info, unsigned int crossing, double prob) - informs
the user's software that an amount prob of probability has crossed
crossing since the last time step.
template< class Function> void apply_to_object_refs( Function& f, Info info) -
templated function that applies an arbitrary function object f to
each object.
Function will have
the () operator overloaded to take an object reference.
double reaction_coordinate( Info info, Obj obj) - returns value
of reaction coordinate of obj. This is used only by the bin-building
process.
Sampler class itself has the following member functions:
void set_number_of_objects( unsigned int n) -
this is the target
number about which the total number of objects
fluctuates at steady-state.
void set_number_of_threads( unsigned int) - sets number of threads.
It is probably best not to use any more threads than processors.
If the number of threads is not set, then only
one thread is created.
void set_number_of_fluxes( unsigned int n) - number of different
crossing, or reaction criteria. If not called, the number of
crossings is one.
void set_random_number_generator( RNG_Ref)
void set_bins( Bins_Ref)
void set_object_container( Objects_Ref objs) - initially
objs should be empty
void initialize_objects() -
populates the object list with objects at the start of the
process.
void step_objects_forward() - once everything is
initialized, each object is stepped forward.
void set_number_of_moves_per_output( unsigned int n) -
If set, flux information is output every n calls
to step_objects_forward; otherwise,
it is output after every call.
void renormalize_objects() - splits and deletes
objects in order to maintain bin populations.
void generate_bins( Info info, std::list< double>& bins) -
places the bin partitions
into the initially empty list bins.
At the end, the reaction coordinates are ordered from highest to zero.
Before calling the function, the number of objects is set
by set_number_of_objects, but the objects themselves
are not yet created. Then, the function is called once.
Normally, after initializing everything, most of the work will
take place in a loop such as this:
for( int i = 0; i < n_moves; i++){
sampler.step_objects_forward();
sampler.renormalize_objects();
}