# Harmonic repulsion - Boltzmann distribution

In [2]:
import os
import numpy as np
import matplotlib.pyplot as plt


In [3]:
readdy.__version__

Out[3]:
'v1.0.1-0'

Prepare the system

In [4]:
system = readdy.ReactionDiffusionSystem(
box_size=(4, 4, 4), periodic_boundary_conditions=[True, True, True], unit_system=None)


Prepare the simulation

In [5]:
simulation = system.simulation(kernel="SingleCPU")

simulation.output_file = "out_rdf.h5"
simulation.observe.rdf(1000, np.linspace(0., 2., 10), ["A"], ["B"], 1. / system.box_volume)


Run the simulation

In [7]:
if os.path.exists(simulation.output_file):
os.remove(simulation.output_file)

simulation.progress_output_stride = 1000
simulation.run(n_steps=100000000, timestep=1e-4, show_system=False)

In [8]:
traj = readdy.Trajectory(simulation.output_file)

In [9]:
def average_across_first_axis(values):
n_values = len(values)
mean = np.sum(values, axis=0) / n_values  # shape = n_bins
difference = values - mean  # broadcasting starts with last axis
std_dev = np.sqrt(np.sum(difference * difference, axis=0) / n_values)
std_err = np.sqrt(np.sum(difference * difference, axis=0) / n_values ** 2)
return mean, std_dev, std_err

return 0.5 * force_const * np.power(r - interaction_radius, 2.)
else:
return 0.

boltz = lambda r: np.exp(-1. * potential(r, force_const, interaction_radius))
r_range = np.linspace(0.1, 2., 100)
b_range = np.fromiter(map(boltz, r_range), dtype=float)
plt.plot(r_range, b_range, label=r"Boltzmann correlation $e^{-\beta U(r)}$")

mean, std_dev, std_err = average_across_first_axis(rdf_values)
plt.xlabel(r"Distance $r$ of A and B")
plt.ylabel(r"Radial distribution function $g(r)$")