Dear developers, thank you for the nice simulation framework. I'm interested in simulating cross-linked biopolymer networks with readdy. I added one feature for my simulations, which I happily share, if you think others might find it useful:

• The descriptor language was extended for topology-topology fusion reactions. One can now specify for example [network>5], which is similar to self fusion with [self=true], but checks whether two particles within the same topology are separated by at least 6 edges from each other.
• Only if the number of edges is larger than the specified value, the fusion can happen.
• This can be used to e.g. have all initial separate polymers fuse to one single network topology.
• With [self=true], I wouldn't be able to do that, since I would not be able to stop direct neighbors within one topology from fusing.
• An example can be seen in the added tests in readdy/test/TestTopologyReactions.cpp, where I created a structure of 4 rectangular polymers, that touch each other on the 4 corners and can fuse to one connected topology.

Units

Introduced a runtime dependency to pint. Units can now be used in all relevant api functions:

Construct a system with a simulation box size of (10, 10, 10) and units of: - length in nanometers - time in nanoseconds - energy in kJ/mol The temperature is defaulted to 2.437 kJ/mol, all boundaries are set to be periodic.

import readdy
system = readdy.ReactionDiffusionSystem([10., 10., 10.])

The simulation box size is the only required constructor argument of a reaction diffusion system.

Construct a system with a simulation box size of (10, 10, 10) and units of: - length in kilometers - time in hours - energy in kcal/mol

import readdy
system = readdy.ReactionDiffusionSystem([10, 10, 10] * readdy.units.meters,
            length_unit='kilometer', time_unit='hour', energy_unit='kilocal/mol')
print(system.box_size)
>>> [ 0.01  0.01  0.01] kilometer
print(system.kbt)
>>> 0.5824569789674953 kilocalorie / mole

The room temperature as well as any other units of time/length/energy will be automatically converted to this particular set of system units.

Construct a unitless system with a simulation box size of (10, 10, 10):

import readdy
system = readdy.ReactionDiffusionSystem([10, 10, 10], length_unit=None, time_unit=None,
            energy_unit=None)
print(system.box_size)
>>> [ 10.  10.  10.]
print(system.kbt)
>>> 2.437

If time, length, and energy units are set to None or empty strings, the system will be treated as completely unitless system.

If properties are set without providing a particular unit, e.g.,

system.potentials.add_harmonic_repulsion("A", "B", force_constant=10., 
    interaction_distance=5.)

the system's units are assumed without further conversion, i.e., energy/length**2 for the force constant, and length for the interaction distance. Providing a unit will trigger a conversion:

from readdy import units as units
sys = readdy.ReactionDiffusionSystem(box_size=[1., 1., 1.])
sys.add_species("A", 1.)
sys.add_species("B", 1.)
sys.potentials.add_harmonic_repulsion("A", "B", 
    force_constant=10. * units.kilocal / (units.mol * units.mile**2), 
    interaction_distance=5.)

or also trigger an error if the dimensionality does not match:

sys.potentials.add_harmonic_repulsion("A", "B", 
    force_constant=10. * units.kilocal / (units.mol * units.mile**3), 
    interaction_distance=5.)

raises an

DimensionalityError: Cannot convert from 'kilocalorie / mile ** 3 / mole' ([mass] / [length] / [substance] / [time] ** 2) to 'kilojoule / mole / nanometer ** 2' ([mass] / [substance] / [time] ** 2)

PBC and box potentials

Upon simulation start there will be a check on the particle types together with periodicity of the simulation box and box potentials, in particular:

If the simulation box is not completely periodic and if the particle types in that box have no box potential to keep them contained in the non-periodic directions, i.e.,

box_lower_left[dim] >= potential_lower_left[dim] 
    or box_upper_right[dim] <= potential_upper_right[dim] 

is true for any non-periodic dim, an exception is thrown.

This is what I got so far after our discussion in issue #191 . It's not finished yet, so no need to merge it right away, but maybe you can already take a look at it. Rate functions for spatial reactions works only for TT reactions yet, because I ran into some issues with std::variant and pushed resolving them aside for now. But I can try to get it working next.

For structural reactions, there a _cumulativeRate gets updated with calls to GraphTopology::updateReactions. I implemented an analogous GraphTopolgy::updateSpatialReactions, which do not update the _cumulativeRate. Instead, the current rates get added to _cumulativeRate in the kernels' EvaluateTopologyReactions::gatherEvents methods, just as it was with constant rates as well. I hope that is o.k. for dynamic rates, too. I'm not completely sure, because the tests that I implemented in TestTopologyReactions.cpp only produce constant rates (the functions simply return a constant). I'll probably have to come up with some distinct tests.

SpatialTopologyReaction::rate() throws an error, if the evaluated rate functions leads to a negative rate. First I had a variant where I instead cut it at 0.0, but I think it's more user friendly if the user gets an error, so he/she gets some feedback that there is something wrong with the rate function provided.

It is currently not possible to define a rate function for spatial topology reactions. I would like to use a function there, to account for an (implicit) background concentration of some protein. Would it be possible to implement it such that a function as well as a float can be used as the rate argument in system.topologies.add_spatial_reaction ?

The KernelContext does not depend on the ReactionFactory anymore and the Kernel itself dispatches the reactions into the context. Further, the argument sequence for reactions has been unified: First the name, then educts, then products, then the educt distance, then the product distance and then weights (some of them omitted depending on which reaction type is considered). Also, with respect to issue #17, added an action parameter (create=default/clear) to the UpdateNeighborList program such that the neighbor list is cleared after a simulation run().

Hi, I'm trying to recover data from a run that was killed (due to timeout). When I try to open the file I get:

t = readdy.Trajectory(file) OSError: Unable to open file (file is already open for write (may use to clear file consistency flags))

I tried: $ h5clear -s file

But now I get: t = readdy.Trajectory(file) RuntimeError: Unable to get link info (bad object header version number)

Is this expected and/or is there some other way I could recover the trajectory?

Thanks Jan

I encountered some issue in long simulations of linear polymers. I wanted to validate that my polymers show expected diffusion of polymers. For that I revert the minimum image projection similar to the way in your MSD validation example. I look at the trajectories of a beads that are rather central in the polymer for my analysis. But the trajectories show a repeating pattern:

I was certain I must have done something stupid in the minimum image reversion, or in interactions between the particles and the polymer topolgy or something like that. But the pattern shows also with the minimum image projection (in some examples where the long-term direction is along one of the cardinal axes, this is nicely visible; sorry the plots are not that nice, color coded is the time, dark blue at the beginning of the simulation, light yellow at the end):

I tried to reproduce this with rather simple code as well, e.g. with:

import numpy as np
import readdy
import os

stride = int(1e5)
n_steps = int(1e8)
progress_output_stride = int(1e5)
dt = 1e-3


system = readdy.ReactionDiffusionSystem(box_size=[60., 60., 60.],
                                        unit_system=None)

system.topologies.add_type("Polymer")
system.add_topology_species("Head", 1.)
system.add_topology_species("Core", 1.)

system.topologies.configure_harmonic_bond("Head", "Core", force_constant=70, length=1.)
system.topologies.configure_harmonic_bond("Core", "Core", force_constant=70, length=1.)
system.topologies.configure_harmonic_angle("Core", "Core", "Core", force_constant=5., equilibrium_angle=np.pi)
system.topologies.configure_harmonic_angle("Core", "Core", "Head", force_constant=5., equilibrium_angle=np.pi)
system.topologies.configure_harmonic_angle("Head", "Core", "Head", force_constant=5., equilibrium_angle=np.pi)


simulation = system.simulation(kernel="SingleCPU")

init_top_pos = np.zeros((10, 3))
init_top_pos[:, 0] = np.arange(-5, 5)
top = simulation.add_topology("Polymer", ["Head"] + ["Core"]*8 + ["Head"], init_top_pos)
for i in range(9):
    top.get_graph().add_edge(i, i+1)

simulation.output_file = "data.h5"
simulation.observe.particle_positions(stride=stride)
simulation.record_trajectory(stride=stride)
simulation.progress_output_stride = progress_output_stride

if os.path.exists(simulation.output_file):
    os.remove(simulation.output_file)

simulation.run(n_steps=n_steps, timestep=dt)

This code leads to similar patterns. And even for diffusion of single particles without any interactions (I had 9 particles in a box simulated for 1e10 time steps, recorded every 1e4 steps, D=1, dt=1e-3, box=[50, 50, 50]), here shown for the first 60000 recorded frames of one of the particles:

For other particles in the box, the same pattern shows, but they start at a different point of the cycle:

I am not sure what to make of this. Did I miss something obvious? Do I have to do something like choose a proper random number generator for rather long simulations? Or is this an expected behaviour in simulations of free diffusion for many time steps due to limitations of pseudo random number generators? Can you reproduce this or is this an artifact specific to my set up (I already tried it on 3 devices, but the set up might be very similar on all)?

Mainly: Decoupled CPU and SingleCPU kernel

Other than that, added a small performance test suite, made core library tests parameterizable with respect to the kernel against which is to be tested.

Hi! I have faced with a non trivial problem to me while trying to run simulation with different initial conditions. I was trying to simulate simple bimolecular reaction A+B->C in two cases:

1. Equal initial concentrations of A and B species D_A = 1, D_B = 1, D_C = 1 L_x, L_y, L_z = 300, 300, 300 R_f = 500 N_A_0 = 50000, N_B_0 = 50000 T = 300 R_AB = 1
2. Species B in excess as compared to A species D_A = 1, D_B = 1, D_C = 1 L_x, L_y, L_z = 1300, 1300, 1300 R_f = 500 N_A_0 = 50000, N_B_0 = 5000000 T = 300 R_AB = 1 In both cases other parameters were identical: simulation = system.simulation(kernel="CPU") simulation.output_file = "A+B-C.h5" simulation.reaction_handler = "Gillespie"

initial_coords_A = L_xnp.random.random((N_A_0, 3)) - L_x/2 initial_coords_B = L_xnp.random.random((N_B_0, 3)) - L_x/2

simulation.add_particles(type = "A", positions = initial_coords_A) simulation.add_particles(type = "B", positions = initial_coords_B)

simulation.observe.number_of_particles(stride=1, types = ["A"])

simulation.run(n_steps=8600, timestep=1e-2) But final step (simulation.run(n_steps=8600, timestep=1e-2)) was failed in the second case (in the excess of B particles). Timestep=1e-2 was the same in the first case and simulation showed nice results.

Hello there,

I've found an error while trying to get the index of the neighbors of a particular vertex. Here is an example of the code I wrote to perform this action:

for v in vertices: a = vertices[v].neighbors()[0] b = vertices.index(a)

The error says: ValueError: Vertex[particleIndex= n, data=, neighbors=[n, n]] is not in list

but if I print "vertices" I see that this vertex ("a") is indeed in the list.

Is it possible to check the index of a particular vertex's neighbor in "vertices"?

Hope it's clear!

Thanks!

Best

Matias

Hi,

from the online documentation I understand, that the function should work as follows:

simulation.make_checkpoints(1000, output_directory="checkpoints/", max_n_saves=10) simulation.make_checkpoints(stride=1000, output_directory="checkpoints/", max_n_saves=10)

However, I get the error, that "output_directory" is an unexpected keyword and if I leave it out (by writing simulation.make_checkpoints(1000)), the simulation runs but no checkpoints are saved.

Is there anything obvious that I am doing wrong?

Thanks in advance, Julia

Hi,

I am using Readdy for a reversible A + B <-> C simulation. I want to compute some statistics to high accuracy and so I am running a lot of trials, which would best be done as a batch job. Is there a way to set the seed for the random number generator for each job?

Sincerely,

Max Heldman

Is it possible to recentre the whole system so that all of the particles a translated by a distance vector [x, y, z] so that a single topology's centre of mass is at [0, 0, 0] ?

Hi there, we apply the simulation for different types of bimolecular reactions, and discovered some problems in the reaction of identical particles like A+A->Product. We use the following setup of ReaDDy sistem: `system = readdy.ReactionDiffusionSystem([L_x, L_y, L_z], temperature=T*readdy.units.kelvin) system.add_species("A", diffusion_constant = D_A) system.add_species("C", diffusion_constant = D_C)

system.reactions.add_fusion(name="fus", type_from1="A", type_from2="A", type_to="C", rate=R_f, educt_distance = R_AA)

initial_coords_A = L_x*np.random.random((N_A_0, 3)) - L_x/2 simulation.add_particles(type = "A", positions = initial_coords_A)`

But analysing the deviation of observable kinetics from the kinetics of the mass action law we found fundamental contradiction. The results shows that the simulation system behaves like there two kinds of particles A (type_from1 and type_from2, and they interacts just between type_from1 and type_from2, but doesn't interact between type_from1 and type_from1, and type_from2 and type_from2 ), while the total number of A particles is equal to 2*N_A_0. It's worth noting that the kinetics of the mass action law in this case is correct. Could you clarify the issue? Thanks!

Hi @clonker and @chrisfroe, have you considered implementing hierarchical grids for polydisperse systems, i.e. linked cells of different cuboid sizes for particles/reactions with different interaction radii? Ogarko2012 Iwai1999 I guess that would be a not so small incision into the core algorithm? My simulations are largely monodisperse, but with a few large particles. Run times with size ratio of 1:7 are still feasible, but this already has a huge impact on performance. I guess the cuboid size is determined by the largest interaction radius? I think this woul lead to a exponential increase in run time with the size ratio, which is what I see in my simulations:

Is it possible to run simulations without any reactions occurring for equilibration purposes.

e.g. something like:

# get Simulation object
simulation = system.simulation()

# equilibrate
simulation.evaluate_topology_reactions = False
simulation.reaction_handler = None
simulation.run(1000)

# run with reactions
simulation.evaluate_topology_reactions = True
simulation.reaction_handler = "Gillespie"
simulation.run(1000000)

What @chrisfroe and I came up with during a discussion:

• the context should be restructured so that there is a separation between system and simulation (similar to python api)
• the kernel is created with a context and instantiates everything as necessary (eg observables etc)
• since the context is then also carrying information about the observables that should be evaluated, it can be logically consistent (things like evaluateVirial() etc can be set automatically in the context)
• simulation is performed by creating actions and performing them or running a default loop