run

Short description

@author: Moritz F P Becker

class pyrid.run.Simulation(box_lengths=array([50., 50., 50.]), dt=0.1, Temp=293.15, eta=1e-21, stride=100, write_trajectory=True, file_path='Files/', file_name='PyRID', fig_path='Figures/', boundary_condition='periodic', nsteps=10000, sim_time=None, wall_force=100.0, seed=None, length_unit='nanometer', time_unit='ns', cells_per_dim=None, max_cells_per_dim=50)

The Simulation class takes care of most of the user communication with PyRIDs numba compiled functions that are the core of the PyRID simulator. The Simulation class thereby represents the biggest part of the user API. By creating an instance of the Simulation class and calling its methods, the user can define molecule types, setup the simulation volume, define reactions etc. and, last but not least, start the simulation.

Attributes

Systemobject

Instance of the System class

Particlesobject

Instance of the Particles class

RBsobject

Instance of the RBs class

file_pathobject

Directory for simulation results. Uses the pathlib library.

file_namestring

Project name for the simulation files.

fig_pathobject

Directory for any figures created with PyRID. Uses the pathlib library.

write_trajectoryboolean

If false, PyRID will not save any molecule trajectory data. Default = True

strideint64

Stride used for saving molecule trajectories to the hard drive.

current_stepint64

Current simulation time step.

release_events = `list`

List of molecule release events

nstepsint64

Number of simulation steps.

sim_timefloat64

Instead of a number of simulation steps, the user can also define the simulation time/duration. This option is useful, e.g., when running PyRID on a cluster with maximum runtime for users.

Observerobject

Instance of the Observables class.

Measureslist

List of all measures / observables supported by PyRID.

length_unitsdict

Dictionary assigning written out unit name and unit short form.

unitsdict

Dictionary containing for each system property / observable the corresponding unit.

time_unit‘ns’, ‘mus’, ‘ms’, ‘s’

Unit of time

checkpointarray like

Structured array that contains fields for different checkpoint parameters (stride, directory, maximum number of saves).

Timerdict

Dictionary that saves the runtime for different parts of the simulation loop such as integration, reaction handling and writing data to files.

progress_propertieslist

List of simulation properties to print/log while the simulation is running. Supported are ‘it/s’ (iterations per second), ‘pu/s’ (particle updates per second), ‘Time passed’, ‘step’, ‘N’ (Number of molecules), ‘Pressure’, ‘Volume’ (Simulation box volume), ‘Vol_Frac’ (Volume fraction per molecule type), ‘Bonds’ (Number of particle-particle bonds).

hdfobject

hdf5 file that PyRID writes all observable data to.

HGridarray like

Hierarchical grid represented by a numpy structures array.

particles_left_boxlist

List of particles that have left the simulation volume. If this list is not empty this is always a sign that some system parameters need to changed because the simulation is unstable.

Methods

print_timer(self)	Prints the runtime of different parts of the simulation loop such as integration, reaction handling and writing data to files.
k_micro(self, mol_type_1, mol_type_2, k_macro, R_react, loc_type = ‘Volume’)	Calculates the microscopic reaction rate from a given macroscopic rate, reaction radius and the educts´ diffusion coefficients. Calls pyrid.reactions.reactions_util.k_micro.
k_macro(self, D1, D2, k_micro, R_react)	Calculates the macroscopic reaction rate from a given microscopic rate, reaction radius and the educts´ diffusion coefficients. Calls pyrid.reactions.reactions_util.k_macro.
add_checkpoints(self, stride, directory, max_saves)	Registers the parameters for saving checkpoints.
add_barostat_berendsen(self, Tau_P, P0, start)	Adds a Berendsen barostat. Calls pyrid.system.system_util.System.add_barostat_berendsen.
set_compartments(self, Compartments, triangles, vertices, mesh_scale = 1.0, adapt_box = False)	Sets up 3d mesh compartments.
register_particle_type(self, Type, radius)	Registers a particle type. Calls pyrid.system.system_util.System.register_particle_type.
add_interaction(self, interaction_type, type1, type2, parameters, bond = False)	Assigns a particle-particle interaction potential to a particle type pair. Calls pyrid.system.system_util.System.add_interaction.
add_up_reaction(self, reaction_type, educt_type, rate, product_types = None, product_loc = None, product_direction = None, radius = None)	Registeres a Uni-Particle (UP) reaction for an educt particle. Calls pyrid.system.system_util.System.add_up_reaction.
add_um_reaction(self, reaction_type, educt_type, rate, product_types = None, product_loc = None, product_direction = None, radius = None)	Registeres a Uni-Molecular (UM) reaction for an educt Molecule. Calls pyrid.system.system_util.System.add_um_reaction.
add_bp_reaction(self, reaction_type, educt_types, product_types, rate, radius, interaction_type = None, interaction_parameters = None)	Registeres a Bi-Particle (BP) reaction for an educt particle. Calls pyrid.system.system_util.System.add_bp_reaction.
add_bm_reaction(self, reaction_type, educt_types, product_types, particle_pairs, pair_rates, pair_radii, placement_factor = 0.5)	Registeres a Bi-Molecular (BM) reaction for an educt molecule. Calls pyrid.system.system_util.System.add_bm_reaction.
register_molecule_type(self, molecule_name, particle_pos, particle_types, collision_type = 0, h_membrane = None)	Regsiters a molecule type. Calls pyrid.system.system_util.System.register_molecule_type.
set_diffusion_tensor(self, molecule_name, D_tt, D_rr)	Checks whether a given diffusion tensor is valid (physically meaningful) and sets the diffusion tensor for a given molecule type. Calls pyrid.system.system_util.System.set_diffusion_tensor.
fixed_concentration_at_boundary(self, molecule_name, C, comp_name, location)	Sets the fixed concentration boundary given a concentartion for each defined molecule type and for each compartents volume and surface. Calls pyrid.system.system_util.System.fixed_concentration_at_boundary.
add_molecules(self, Location, Compartment_Number, points, quaternion, points_type, face_ids = None)	Places molecules in the simulation box given their location, orientation, type, and compartment number.
distribute(self, Method, Location, Compartment_Number, Types, Number, clustering_factor=1, max_trials=100, jitter = None, triangle_id = None, multiplier = 100, facegroup = None)	Calculates the distribution of molecules in the simulation box / the surface or volume of a compartment using one of several different methods.
observe(self, Measure, molecules = ‘All’, reactions = ‘All’, obs_stride = 1, bin_width = None, n_bins = None, stepwise = True, binned= False)	Sets up an observer for a certain system property / measure.
observe_rdf(self, rdf_pairs = None, rdf_bins = 100, rdf_cutoff = None, stride = None)	Sets up an observer for the radial distribution function.
add_release_site(self, Location, timepoint, compartment_id, Number = None, Types = None, origin = None, jitter = None, Molecules_id = None, Positions = None, Quaternion = None, triangle_id = None)	Adds a molecule release site that releases a number of molecules at a defined time point into the system.
progress(self, sim_time, nsteps, j, time_passed, it_per_s, Np, Pressure_total, N_bonds, Volume, out_linebreak, N, width=30)	prints / logs several different simulation parameters such as iterations per second, molecule number and pressure.
load_checkpoint(self, file_name, index, directory = ‘checkpoints/’)	Loads a system state from a checkpoint file.
run(self, progress_stride = 100, keep_time = False, out_linebreak = False, progress_bar_properties = None, start_reactions = 0)	Simulation loop. Runs the simulation.

add_barostat_berendsen(Tau_P, P0, start)

Adds a Berendsen barostat. Calls pyrid.system.system_util.System.add_barostat_berendsen.

Parameters

Tau_Pfloat

Time constant of berendsen barostat

P0float

Target pressure of berendsen barostat

startint

Simulation step from which to start simulating in NPT ensemble.

Notes

The Berendsen barostat is frequently used in molecular dynamics simulations due to its simplicity and fast equilibration. While the Berendsen barostat results in the correct particle density, it does not correctly sample from the (NPT) thermodynamical ensemble, i.e. the fluctuation in the pressure are not 100% accurate. Therefore, one needs to be carefull when to use this Barostat. However, if the barostat is used in the preparation step but not during sampling, e.g. when doing LLPs simualtions where the barostat is used to relax the system’s shape but afterwards one simulates in the NVT ensemble Muller et al. [32], this should not be a problem. Note: Baruffi et al. [33] introduced a barostat for overdamped Langevin dynamics that settles to the correct thermodynamic equilibrium distribution.

add_bm_reaction(reaction_type, educt_types, product_types, particle_pairs, pair_rates, pair_radii, placement_factor=0.5)

Registeres a Bi-Molecular (BM) reaction for an educt molecule. Calls pyrid.system.system_util.System.add_bm_reaction.

Parameters

reaction_typestr {‘fusion’, ‘enzymatic_mol’}

Name of the reaction type

educt_typearray_like

Names of the educt molecule types

particle_pairsarray_like

List of educt particle pairs belonging to each of the two educt molecules

pair_ratesarray_like

List of reaction rates for each particle pair

pair_Radiiarray_like

List of reaction radii for each particle pair

placement_factorfloat64 [0,1]

Only for fusion reactions: Affects where between the educt molecules the product molecule is placed. A factor of 0.5 means in the middle. A smaller factor will shift the product towards the first educt, a larger value towards the second educt molecule. A value different from 0.5 may increase accuracy in the case of a strong size difference between the two educts (a strong difference in the diffusion tensors). Default = 0.5

Notes

Bi-molecular reactions are evaluated using the Doi scheme. Thereby, two educt particles react with a reaction rate k if the inter-particle distance is below the reaction radius. PyRID only evaluates the eucludean distances between particles but not molecule centers. Thereby, also for bi-molecular reactions educt particles need to be set. Since molecules can consist of several particles, a list of particle pairs can be passed and for each particle pair a reaction rate and radius is defined. If one of the pair particle reactions is successful, the corresponding molecule that the particle belongs to is converted as defined by the reaction type. If the molecule cosnsists of a single particle, the reaction will be handled just as the bi-particle reactions. The user needs to take care of calculating accurate reaction radii and rates such that the desired reaction dynamics are simulated, which can become difficult for large molecules. Therefore it is sometimes more convenient to place an interaction free particle at the origin of the molecule and execute the bimolecular reaction via this particle. This method is equivalent to the case where one would evaluate reactions based on the distance between molecule centers. As for the uni-molecular reactions, several reaction paths can be defined. Each reaction path may involve different particle pairs of the educt molecules. On a side note: Thereby, orientation dependent bimolecular reactions are in principle possible. However, since the reaction probability does not scale with the distance, directionality is only really accounted for if the reaction radius is small compared to the molecule size.

add_bp_reaction(reaction_type, educt_types, product_types, rate, radius, interaction_type=None, interaction_parameters=None)

Registeres a Bi-Particle (BP) reaction for an educt particle. Calls pyrid.system.system_util.System.add_bp_reaction.

Parameters

reaction_typestr {‘bind’, ‘absorption’, ‘enzymatic’}

Name of the reaction type

educt_typearray_like

Names of the educt particle types

product_typesarray_like

List of products (only used for binding reactions)

ratefloat64

Reaction rate

radiusfloat64

Reaction radius.

interaction_typestr

Name of the energy potential function in case ‘reaction_type’ is a binding reaction.

interaction_parametersdict

Parameters for the interaction function in case ‘reaction_type’ is a binding reaction.

Notes

Bi-particle reactions are evaluated using the Doi scheme. Thereby, two educt particles react with a reaction rate k if the inter-particle distance is below the reaction radius.

Binding Reactions

A particle pair will only enter a bound state if the binding reaction was succesfull. Otherwise two particles won’t experience any pairwise interaction force! A reaction occurs if two particles for whose type a binding reaction has been defined are within the reaction radius. The reaction probability is evaluated as:

\[1-exp(\lambda \cdot \Delta t),\]

where \(\lambda\) is the reaction rate and \(\Delta t\) is the integration time step.

add_checkpoints(stride, directory, max_saves)

Registers the parameters for saving checkpoints such as stride, maximum number of saves and target directory.

Parameters

strideint64

Stride for making checkpoints.

directorystring

Target directory for the checkpoint files.

max_savesint64

Maximum number of files that are saved. Old files will be overwritten to keep the number of checkpoint files at the value of max_saves.

add_interaction(interaction_type, type1, type2, parameters, bond=False)

Assigns a particle-particle interaction potential to a particle type pair. Calls pyrid.system.system_util.System.add_interaction.

Parameters

interaction_typestr {‘harmonic_repulsion’, ‘PHS’, ‘harmonic_attraction’, ‘CSW’}

Type of interaction potential

type1str

Names of the educt particle type 1

type2str

Names of the educt particle type 2

parametersdict

Parameters for the interaction potential function.

bondbool

True if interaction is initialized as part of a bindign reaction. (Default: False)

breakablebool

If set to True, bonds that have been formed by a binding reaction are deleted if the interparticle distance is above the cutoff radius of the interaction potenial.

Notes

Bi-particle reactions are evaluated using the Doi scheme. Thereby, two educt particles react with a reaction rate k if the inter-particle distance is below the reaction radius.

add_molecules(Location, Compartment_Number, points, quaternion, points_type, face_ids=None)

Places molecules in the simulation box given their location, orientation, type, and compartment number.

Parameters

Location‘Volume’ or ‘Surface’’

Determines whether the molecules are surface or volume molecules.

Compartment_Numberint64

The number of the compartment in which to place the molecules

pointsfloat64[:,3]

Molecule coordinates

quaternionfloat64[:,4]

Orientation of the molecules represented by rotation /orientation quaternions.

points_typeint64[:]

Molecule type ids.

face_idsint64[:]

Indices of the mesh triangle/faces on which the molecules are placed. Default = None

Raises

raise IndexError(‘Compartment number must be a positive integer!’)

Compartment number must be an integer value.

raise IndexError(‘Compartment number must be greater zero!’)

Compartment numbers must be greater than zero.

raise TypeError(‘Missing required argument face_ids’)

In the case of Location = ‘Surface’ face_ids must be passed.

NameError(‘Location must be either Volume or Surface!’)

There are only two locations a molecule can assigned to: the surface or the volume of a compartment.

add_release_site(Location, timepoint, compartment_number, Number=None, Types=None, origin=None, jitter=None, Molecules_id=None, Positions=None, Quaternion=None, triangle_id=None)

Adds a molecule release site that releases a number of molecules at a defined time point into the system. The molecules are distributed normally around a point in the volume or a mesh surface.

Parameters

Location‘Volume’ or ‘Surface’

Determines of whether the molecules are released on the compartment surface or in its volume.

timepointint64

Time step at which to release the molecules.

compartment_numberint64

Number of the mesh compartment. 0 means simulation box, i.e. outside any mesh compartment.

Numberint64[:]

Number of molecules per molecule type.

Typeslist of strings

List of molecule types to release.

originfloat64[3]

Only if Location = ‘Volume’. Position around which to distribute the molecules upon release.

jitterfloat64

Initial width of the molecule distribution. Molecules are distributed normally around the origin.

triangle_idint64

Only if Location=’Surface’. Index of the triangle around whose centroid the molecules are distributed.

Molecules_idint64[:]

Only for Location=’Volume’. In case you want a custom distribution of molecules this list contains the type id of each molecule.

Positionsfloat64[:,3]

Only for Location=’Volume’. In case you want a custom distribution of molecules this list contains the coordinates of each molecule.

Quaternionfloat64[:,4]

Only for Location=’Volume’. In case you want a custom distribution of molecules this list contains the orientation of each molecule.

add_um_reaction(reaction_type, educt_type, rate, product_types=None, product_loc=None, product_direction=None, radius=None)

Registeres a Uni-Molecular (UM) reaction for an educt Molecule. Calls pyrid.system.system_util.System.add_um_reaction.

Parameters

reaction_typestr {‘production’, ‘fission’, ‘conversion_mol’, ‘decay’}

Name of the reaction type

educt_typestr

Name of the educt molecule type

ratefloat64

Reaction rate

product_typesarray_like

List of products

product_locarray_like

List specifying whether the individual products are volume (0) or surface (1) molecules

product_directionarray_like

Specifies whether, in case of a surface molecule releasing products into the volume, the realease occurs inside (-1) or outside (1) the mesh compartment.

radiusfloat64

The radius within which molecules are distributed upon a successfull release reaction.

Notes

Uni-molecular reactions are evaluated using the Gillespie stochastic simulation algorithm. Thereby, for a molecule type, the time point of the next uni-molecular reaction occuring is calculated in advance and executed when the simualtion reaches the respective time point. A uni-molecular reaction can consist of several reaction paths. Each molecule type has exactly one uni-molecular reaction id registered for which the time point is drawn randomly based on the total reaction rate which is the sum of all path reaction rates. At the time point of the reaction occuring, a path is randomly drawn from the set of possible paths weighted by the individual path reaction rates. To illustrate the above by an example: For Molecule ‘A’ a we may define two convcersion reactions where A -k1-> B, A -k2-> C. In addition, we define a fission reaction A -k3-> D+E. The next reaction will occur at a time point T that is calculated based on the total rate k1+k2+k3. At time T, one of the three reactions is executed. The probability for each reaction to be selected being k1/(k1+k2+k3), k2/(k1+k2+k3), k3/(k1+k2+k3).

add_up_reaction(reaction_type, educt_type, rate, product_types=None, product_loc=None, product_direction=None, radius=None)

Registeres a Uni-Particle (UP) reaction for an educt particle. Calls pyrid.system.system_util.System.add_up_reaction.

Parameters

reaction_typestr {‘relase’, ‘conversion’}

Name of the reaction type

educt_typestr

Name of the educt particle type

ratefloat64

Reaction rate

product_typesarray_like

List of products

product_locarray_like

List specifying whether the individual products are volume (0) or surface (1) molecules

product_directionarray_like

Specifies whether, in case of a surface particle releasing products into the volume, the realease occurs inside (-1) or outside (1) the mesh compartment.

radiusfloat64

The radius within which molecules are distributed upon a successfull release reaction.

Notes

Uni-particle reactions are evaluated using a variant of the Gillespie stochastic simulation algorithm. Thereby, for a particle type, the time point of the next uni-particle reaction occuring is calculated in advance and executed when the simualtion reaches the respective time point. For unimolecular reactions we can draw the time point of the next reaction from a distribution using a variant of the Gillespie Stochastic Simulation Algorithm (SSA) Stiles and Bartol [11], Erban et al. [12]. For a single molecule having \(n\) possible transition reactions/reaction paths each having a reaction rate \(k_i\), let \(k_t = \sum_i^n k_i\) be the total reaction rate.

Now, let \(\rho(\tau) d\tau\) be the probability that the next reaction occurs within \([t+\tau,t+\tau+d\tau)\) and let \(g(\tau)\) be the probability that no reaction occurs within \([t,t+\tau)\). The probability that a reaction occurs within the time interval \(d\tau\) is simply \(k_t d\tau\). Thereby

\[\rho(\tau) d\tau = g(\tau) k_t d\tau\]

where \(g(\tau) = e^{-k_t \tau}\) (see also Poisson process). From the above equation we easily find \(P(\tau) = 1-e^{-k_t \tau}\) by integration. To sample from this distribution, we can use the inverse distribution function.

\[\tau = P^{-1}(U)\]

where \(U\) i uniformly distributed in \((0,1)\). Given that \(U = P(\tau) = 1-e^{-k_t \tau}\), we find \(P^{-1}(U) = \frac{-log(1-U)}{k_t}\). Since \(U\) is uniformly distributed in 0,1, so is \(1-U\). Thereby, we can draw the time point of the next reaction from:

\[\tau = \frac{1}{k_t} \ln\Big[\frac{1}{U}\Big],\]

With the above method, we accurately sample from the distribution of expected molecule lifetimes \(\rho(\tau) = k_t e^{-k_t \tau}\).

A uni-particle reaction can consist of several reaction paths. Each particle type has exactly one uni-particle reaction id registered for which the time point is drawn randomly based on the total reaction rate which is the sum of all path reaction rates. At the time point of the reaction occurring, a path is randomly drawn from the set of possible paths weighted by the individual path reaction rates. To get the reaction path we compare a second random number, uninformly distributed in (0,k_t), with the cumulative set of reaction rates \((k_1, k_1+k_2, ..., k_t)\). The comparison is done via a bisection algorithm (see rand_util.random_util.bisect_right()).

To illustrate the above by an example: For particle ‘A’ a we may define two convcersion reactions where A -k1-> B, A -k2-> C. In addition, we define a release reaction A -k3-> 2*D+5*E. The next reaction will occur at a time point T that is calculated based on the total rate k1+k2+k3. At time T, one of the three reactions is executed. The probability for each reaction to be selected being k1/(k1+k2+k3), k2/(k1+k2+k3), k3/(k1+k2+k3).

distribute(Method, Location, Compartment_Number, Types, Number, clustering_factor=1, max_trials=100, jitter=None, triangle_id=None, multiplier=100, facegroup=None)

Calculates the distribution of molecules in the simulation box / the surface or volume of a compartment using one of several different methods. The available methods are:

• ‘PDS’: The ‘PDS’ (Poisson-Disc-Sampling) method does resolve collisions based on the molecule radius but can only reach volume fractions of <30% for mono-disperse systems. Note that computation times can become very long if the target number of molecules is much larger than would fit into the compartment. The molecule density can also be increased by increasing max_trials at the cost of longer runtimes.In additions one can play around with the clustering factor that affects the distribution of molecules. A larger clustering factor will result in a more uniform distribution of molecules. See pyrid.system.distribute_vol_util.pds

• ‘PDS_unfiform’: The same holds true for ‘PDS_unfiform’ method, however, with the addition that the molecules are uniformly distributed. The method first creates a large number of uniformly distributed points (using the ‘MC’ method). From this sampling pool only those points that are a minimum distance apart from each other are kept. The molecule density can be increased by the multiplier parameter at the cost of longer runtimes.See pyrid.system.distribute_vol_util.pds_uniform and pyrid.system.distribute_surface_util.pds

• ‘Gauss’: The ‘Gauss’ method can be used to normally distribute molecules around the centroid of a triangle (given by triangle_id) on a mesh surface. The distribution width can be set by the jitter parameter. See pyrid.system.distribute_surface_util.normal

• ‘MC’: The ‘MC’ (Monte-Carlo) method is fast and can reach arbitrary molecule densities, however, it does not resolve any collision. See pyrid.system.distribute_surface_util.mc and pyrid.system.distribute_vol_util.mc

Parameters

Method‘PDS’, ‘PDS uniform’, ‘Gauss’, ‘MC’

Method for calculating the molecule distribution.

Location‘Volume’ or ‘Surface’

Determines of whether the molecules are distributed on the compartment surface or in its volume. Note that the method ‘pds_uniform’ is only available for Location=’Volume’ and the method ‘Gauss’ only for Location=’Surface’.

Compartment_Numberint64

Number of the mesh compartment. 0 means simulation box, i.e. outside any mesh compartment.

Typesint64[:]

List of molecule type ids.

Numberint64[:]

Number of molecules to distribute per type.

clustering_factorfloat64

Only for Location=’Volume’ and Method=’PDS’. The clustering factor affects the distribution of molecules. A larger clustering factor will result in a more uniform distribution of molecules.

max_trialsint64

Only for Location=’Volume’ and Method=’PDS’. Increases the number attempts to randomly place a molecule in the vicinity of another molecule.

jitterfloat64

Only for Location=’Surface’ and Method=’Gauss’. Width of the normal distribution.

triangle_idint64

Only for Location=’Surface’ and Method=’Gauss’. Index of the triangle around whose centroid molecules are normally distributed.

multiplierint64

Only for (Location=’Volume’ and Method=’PDS_uniform’) or (Location=’Surface’ and Method=’PDS’). Increases the sampling pool size.

facegroupstring

Only for Location=’Surface’ and Method=’PDS’ or ‘MC’. Name of the face group on whose triangles to distribute the molecules.

Returns

dtype

Some information

Raises

raise NameError(‘Method {} does not exist!’.format(Method))

Available methods are ‘PDS’, ‘PDS uniform’, ‘Gauss’, ‘MC’.

raise IndexError(‘Compartment number must be a positive integer!’)

Compartment number must be an integer value.

raise IndexError(‘Compartment number must be greater zero!’)

Compartment numbers must be greater than zero.

raise TypeError(‘Missing required argument face_ids’)

In the case of Location = ‘Surface’ face_ids must be passed.

NameError(‘Location must be either Volume or Surface!’)

There are only two locations a molecule can assigned to: the surface or the volume of a compartment.

fixed_concentration_at_boundary(molecule_name, C, comp_name, location)

Sets the fixed concentration boundary given a concentartion for each defined molecule type and for each compartents volume and surface. Calls pyrid.system.system_util.System.fixed_concentration_at_boundary.

Parameters

molecule_namestr

Name of the molecule

Cfloat

concentration of molecule type outside simulation box (‘virtual molecules’)

comp_namestr

Name of the compartment

locationint {0,1}

Location of the molecule. 0: Volume, 1: Surface

Notes

fixed_concentration_at_boundary() calculates some properties that are necessary to properly distribute molecules inside the simulation box that hit the simulation box boundary from the outside (Thereby, ‘virtual molecules’ become real molecules in our simualtion). The number of hits per time step a boundary of area A experiences is \(N = l_{perp}*A*C\). Where \(C\) is the concentration in molecules per volume and \(l_{perp}\) is the average net displacement in one tiem step towards or away from any plane, where \(l_{perp} = \sqrt{(4*D*\Delta t/\pi)}\) Stiles and Bartol [11].

k_macro(D1, D2, k_micro, R_react)

Calculates the macroscopic reaction rate from a given microscopic rate, reaction radius and the educts´ diffusion coefficients. Calls pyrid.reactions.reactions_util.k_macro.

Parameters

D1float

Translational diffusion coefficient of educt 1.

D2float

Translational diffusion coefficient of educt 2.

k_microfloat

Microscopic reaction rate.

R_reactfloat

Reaction radius.

Returns

float64

Macroscopic reaction rate

k_micro(mol_type_1, mol_type_2, k_macro, R_react, loc_type='Volume')

Calculates the microscopic reaction rate from a given macroscopic rate, reaction radius and the educts´ diffusion coefficients. Calls pyrid.reactions.reactions_util.k_micro.

Parameters

mol_type_1string

Name of educt moelcule 1.

mol_type_2string

Name of educt moelcule 2.

k_macrofloat

Macroscopic reaction rate.

R_reactfloat

Reaction radius.

loc_typestring

Location type of the educts (‘Volume’ or ‘Surface’).

Returns

float64

Microscopic reaction rate

Notes

The method used here is only valid for volume molecules, i.e. diffusion in 3D. For surface molecules a value is also returned, however with a warning.

load_checkpoint(file_name, index, directory='checkpoints/')

Loads a system state from a checkpoint file.

Parameters

file_namestring

Checkpoint file name.

indexint64

Index of the checkpoint file (index goes from 0 to maximum number of saves).

directorystring

Directory of the checkpoint files. Default = ‘checkpoints/’

observe(Measure, molecules='All', reactions='All', obs_stride=1, stepwise=True, binned=False)

Sets up an observer for a certain system property / measure.

Parameters

Measure‘Position’, ‘Energy’, ‘Pressure’, ‘Volume’, ‘Virial’, ‘Virial Tensor’, ‘Number’, ‘Orientation’, ‘Bonds’, ‘Reactions’, ‘Force’, ‘Torque’

System property which to observe.

moleculeslist or ‘All’

List of molecule types. Default = ‘All’

reactionslist or ‘All’

List of reaction types to observe. Default = ‘All’

obs_strideint64

Stride by which observables are saved to the hdf5 file.

stepwiseboolean

If True, observables are saved in a stepwise manner, i.e. the current value at the respective timestep is saved. Default = True

binnedboolean

If True, observables are binned where the bin size is equal to stride.

Raises

warnings.warn(‘Warning: Binning is not supported for property {0}. PyRID will instead sample {0} stepwise!’.format(Measure))

Binning is not supported for ‘Force’, ‘Torque’, ‘Orientation’, ‘Position’.

observe_rdf(rdf_pairs=None, rdf_bins=100, rdf_cutoff=None, stride=None)

Sets up an observer for the radial distribution function.

Parameters

rdf_pairsnested list

List of molecule type pairs for which to calculate the radial distribution function. Example: [[‘A’,’A’],[‘A’,’B’]].

rdf_binsint64[:]

Number of bins used for each respective molecule pair (resolution of the rdf histogram).

rdf_cutofffloat64[:]

Cutoff distance for each respective molecule pair (choose wisely).

strideint64

Stride by which radial distribution function is sampled and saved to the hdf5 file. Note: calculation of the rdf is computationally expensive!

print_timer()

Prints the runtime of different parts of the simulation loop such as integration, reaction handling and writing data to files.

progress(sim_time, nsteps, j, time_passed, it_per_s, Np, Pressure_total, N_bonds, Volume, out_linebreak, N, width=30)

prints / logs several different simulation parameters such as iterations per second, molecule number and pressure. Also prints a progress bar.

Parameters

sim_timefloat64

Target simulation runtime.

nstepsint64

Target number of simulation steps.

jint64

Current time step.

time_passedfloat

Current simulation runtime.

it_per_sfloat64

Iterations per second.

Npint64

Number of particles / beads in the system.

Pressure_totalfloat64

Current pressure.

N_bondsint64

Current number of particle-particle bonds.

Volumefloat64

Current simulation box volume.

out_linebreakboolean

If True, after each progress print there is a line break.

Nint64

Current number of molecules in the system.

widthint64

Width of the progress bar.

register_molecule_type(molecule_name, particle_pos, particle_types, collision_type=0, h_membrane=None)

Regsiters a molecule type. Calls pyrid.system.system_util.System.register_molecule_type.

Parameters

molecule_namestr

Name of the molecule.

particle_posint[:,3]

Positions of each of the molecule’s particles.

particle_typesarray_like (dtype = ‘U20’)

Array of names for each of the molecule’s particles.

collision_typeint {0,1}

Collision type of the molecule (Default = 0). 0: Collisions with compartments are handled via interaction force. 1: Collisions with compartments are handled by raytracing algorithm (recommended for molcules whose diffusion speed and/or interaction radius is small compared to the chosen integration time step).

h_membranefloat64

Distance of the molecules´ center of diffusion from the membrane surface. Default = None

register_particle_type(Type, radius)

Registers a particle type. Calls pyrid.system.system_util.System.register_particle_type.

Parameters

Typestr

Name of the particle type.

radiusfloat

radius of the particle.

release_molecules(release_event)

A brief description of what the function (method in case of classes) is and what it’s used for

Parameters

parameter_1dtype

Some Information

parameter_2dtype

Some Information

Returns

dtype

Some information

Raises

NotImplementedError (just an example)

Brief explanation of why/when this exception is raised

run(progress_stride=100, keep_time=False, out_linebreak=False, progress_bar_properties=None, start_reactions=0)

Simulation loop. Runs the simulation.

Parameters

progress_strideint64

Stride by which the progress bar is updated. Default = 100

keep_timeboolean

If true, different parts of the simulation loop are timed such as integration of the equations of motion or writing data to files. Default = False

out_linebreak

If True, after each progress print there is a line break. Default = False

progress_bar_propertieslist of strings

List of properties to print/log with the progress bar during simulation. Available are: ‘it/s’, ‘pu/s’, ‘Time passed’, ‘step’, ‘N’, ‘Pressure’, ‘Volume’, ‘Vol_Frac’, ‘Bonds’ or ‘All’.

start_reactions = `int64`

Set the time step from which reactions are allowed. Default = 0

Raises

raise ValueError(‘progress_bar_properties must be a list of strings!’)

If progress_bar_properties is not None or ‘All’ a list of strings must be passed.

raise ValueError(str(prop)+’is not a supported property to be displayed in the progress bar!’)

A progress bar property has been passed that is not available.

set_compartments(Compartments, triangles, vertices, mesh_scale=1.0, adapt_box=False)

Sets up 3d mesh compartments.

Parameters

Compartmentsdict

Dictionary containing for each compartment number the triangle indices and face groups. Also see pyrid.geometry.load_wavefront.load_compartments.

trianglesint64[:,3]

List of vertex indices (3) per triangle.

verticesfloat64[:,3]

List of vertex coordinates.

meshscalefloat64

Factor by which to scale the 3d mesh. Default = 1.0

adapt_boxboolean

If true, the simulation box size is adapted such that the mesh fits exactly in the simulation box. Default = False

set_diffusion_tensor(molecule_name, D_tt, D_rr)

Checks whether a given diffusion tensor is valid (physically meaningful) and sets the diffusion tensor for a given molecule type. Calls pyrid.system.system_util.System.set_diffusion_tensor.

Parameters

molecule_namestr

Some Information

D_ttfloat[3,3]

Translational diffusion tensor ()

D_rrfloat[3,3]

Rotational diffusion tensor

Notes

PyRID offers a module for the calculation of the diffusion tensors of rigid bead models (see molecules_util.hydro_util). It uses the Kirkwood-Riseman calculation with modified Oseen tensor as introduced by Carrasco and de la Torre [2], Carrasco and de la Torre [16]. The same method has been implemented in the HYDRO++ tool, which is free to use but not open source! As such, you may also use HYDRO++ for the calculation of the diffusion tensors. Or even better, use both and tell me if you you get inconsistent results ;) !

Note

Diffusion tensor should be a real positive semidefinite matrix!

pyrid.run.convert_dict(untyped_dict)

Converts a python dictionary to a Numba typed dictionary.

Parameters

untyped_dictdict

Untyped python dictionary

Returns

nb.types.DictType

Numba typed dictionary

Raises

NotImplementedError (just an example)

Brief explanation of why/when this exception is raised

pyrid.run.random_quaternion()

Returns a random rotation quaternion.

Returns

float64[4]

Unit quaternion representing a uniformly distributed random rotation / orientation.

Notes

The method can be found in K. Shoemake, Graphics Gems III by D. Kirk, pages 124-132 [28] and on http://planning.cs.uiuc.edu/node198.html:

\[\boldsymbol{q} = [\sqrt{1-u_1} \, \sin(2 \pi u_2), \sqrt{1-u_1} \, \cos(2 \pi u_2), \sqrt{u_1} \, \sin(2 \pi u_3), \sqrt{u_1} \, \cos(2 \pi u_3)] ,\]

where \(u_1,u_2,u_3 \in [0,1]\) are uniformly distributed random numbers. \(\boldsymbol{q}\) is a uniformly distributed, random rotation quaternion.

pyrid.run.update_pressure(Simulation)

Calculates the pressure from the virial tensor.

Parameters

Simulationobject

Instance of the Simulation class.

Returns

tuple(float, float64[3,3])

Total pressure and pressure tensor.

Notes

When doing rigid body molecular dynamics or Brownian dynamics simulations, we need to be careful how to calculate the virial tensor. Taking the interparticle distances will result in the wrong pressure. Instead, one needs to use the pairwise distance between the centers of mass of the rigid bodies [35] :

(31)\[P_{mol} = P_{mol}^{kin} + \frac{1}{6 V} \sum_{i=1}^{N} \sum_{j \neq}^{N} \langle \boldsymbol{F}_{ij} \cdot (\boldsymbol{R}_i - \boldsymbol{R}_j) \rangle,\]

where \(V\) is the total volume of the simulation box, \(\boldsymbol{F}_{ij}\) is the force on particle i exerted by particle j and \(\boldsymbol{R}_{i}, \boldsymbol{R}_{j}\) are the center of mass of the rigid body molecules, not the center of particles i and j! In Brownian dynamics simulations, \(P_{mol}^{kin} = N_{mol} k_B T\), where \(N_{mol}\) is the number of molecules. Also, the origin of molecules is represented by the center of diffusion around which the molecule rotates about, which is not the center of mass [15]! The net frictional force and troque act thorugh the center of diffusion. This is because when doing Brownian dynamics (and equaly for Langevin dynamics), we do account for the surrounding fluid. Different parts of the molecule will therefore interact with each other via hydrodynamic interactions/coupling. As a result, the center of the molecule (around which the molecule rotates in response to external forces) is not the same as the center of mass, which assumes no such interactions (the molecule sites in empry space). However, for symmetric molecules, the center of mass and the center of diffusion are the same!