import os import numpy as np import matplotlib.pyplot as plt import readdy # set up reaction diffusion system system = readdy.ReactionDiffusionSystem( box_size=[200.,200.,200.], unit_system={"length_unit": "nanometer", "time_unit": "microsecond"}) system.add_species("A", diffusion_constant=1.0) system.add_species("B", diffusion_constant=1.0) system.add_species("C", diffusion_constant=1.0) # Reaction rate value is lambda for Erban-Chapman equation for macroscopic # For k = 1 nm^3/us, use radius = 1 nm, lambda = 0.250695 # For k = 10 nm^3/us, use radius = 1 nm, lambda = 4.52817 # For k = 100 nm^3/us, use radius = 5 nm, lambda = 1.91767 system.reactions.add("myfusion: A +(1) B -> C", 0.250695) #system.reactions.add("myfusion: A +(1) B -> C", 4.52817) #system.reactions.add("myfusion: A +(5) B -> C", 1.91767) # set up corresponding simulation simulation = system.simulation(kernel="CPU") simulation.output_file = "out.h5" simulation.reaction_handler = "Gillespie" n_particles = 8000 initial_positions_a = np.random.random(size=(n_particles, 3)) * 200. - 100. initial_positions_b = np.random.random(size=(n_particles, 3)) * 200. - 100. simulation.add_particles("A", initial_positions_a) simulation.add_particles("B", initial_positions_b) simulation.observe.number_of_particles(stride=1, types=["A"]) # simulate if os.path.exists(simulation.output_file): os.remove(simulation.output_file) dt = 0.1 #dt = 0.01 simulation.run(n_steps= int(200./dt), timestep=dt) traj = readdy.Trajectory(simulation.output_file) time, counts = traj.read_observable_number_of_particles() np.savetxt('ReaDDyBireact_out.txt',counts)