The goal here is to construct the one-particle distribution function of a system evolving under Brownian dynamics; one has to produce a random number drawn from a Gaussian distribution. To construct this quantity, I am thinking of running several simulations, and for specific times in each of the simulations, save the distances of each particle from the center of the 2D square and only in the end, create a histogram of all the values.
My problem is that during each of the simulations, time begins from zero and goes on with a certain time-step, for each of which the particles move randomly. So, the distances to be saved have to be labeled correctly for their corresponding times.
So, my thought was to create an array that will have in each row, 5 sub-arrays, one for each time I want to save the distances of the particles from the center of the square. I am trying to work this with numpy, but with no success. For each simulation, and for specific times, I create an array with all the distances, and I try to append it with numpy.append to the specific subarray, but this doesn't work correctly; as I understand the problem lies in the fact that I don't know how to index properly( and for all the simulations) the sub-arrays.
Beyond that, I think that the approach is not the best: either I will have to abandon the idea of using numpy and figure out how I can index with two indices the array properly, or figure out a way to use numpy more effectively.
So, to the point, the general question here is how I could add/append values to specific sub-arrays of an array ( either pre-constructed with numpy or not and treated as a list).
The alternative would be for someone to mention a more efficient way of creating the one-particle distribution function of a Brownian motion problem, which would be really helpful.
I am adding the relevant code below. Thank you all in advance.
Code:
import random
import math
import matplotlib.pyplot as plt
import numpy as np
# def dump(particles,step,n):
# fileoutput = open('coord.txt', 'a')
# fileoutput.write("ITEM: TIMESTEP \n")
# fileoutput.write("%i \n" % step)
# fileoutput.write("ITEM: NUMBER OF ATOMS \n")
# fileoutput.write("%i \n" % n)
# fileoutput.write("ITEM: BOX BOUNDS \n")
# fileoutput.write("%e %e xlo xhi \n" % (0.0, 100))
# fileoutput.write("%e %e xlo xhi \n" % (0.0, 100))
# fileoutput.write("%e %e xlo xhi \n" % (-0.25, 0.25))
# fileoutput.write("ITEM: ATOMS id type x y z \n")
# i = 0
# while i < n:
# x = particles[i][0]
# y = particles[i][1]
# #fileoutput.write("%i %i %f %f %f \n" % (i, 1, x*1e10, y*1e10, z*1e10))
# fileoutput.write("%i %i %f %f %f \n" % (i, 1, x, y, 0))
# i += 1
# fileoutput.close()
num_sims = 2
N = 49
L = 10
meanz = 0
varz = 1
sigma = 1
# tau = sigma**2*ksi/(kT)
# Starting time
t_0 = 0
# Time increments
dt = 10**(-4) # dt/tau
# Ending time
T = 10**2 # T/tau
# Produce random particles and avoid overlap:
particles = np.full((N, 2), L/2)
times = np.arange(t_0, T, dt)
check = 0
distances = np.empty([50*num_sims, 5])
for sim in range(0, num_sims):
step = 0
t_index = 0
for t in times:
r=[]
for i in range(0,N):
z = np.random.normal(meanz, varz)
particles[i][0] = particles[i][0] + ((2*dt*sigma**2)**(1/2))*z
z = random.gauss(meanz, varz)
particles[i][1] = particles[i][1] + ((2*dt*sigma**2)**(1/2))*z
if (t%(2*(10**5)*dt) == 0):
for j in range (0,N):
rj = ((particles[j][0]-L/2)**2 + (particles[j][1]-L/2)**2)**(1/2)
r.append(rj)
distances[t_index] = np.append(distances[t_index],r)
t_index += 1
