Commit aa0bd50a authored by Olaf's avatar Olaf

Remove trivial comments

parent 5d187866
...@@ -74,7 +74,7 @@ def func_position(x0, v0, f, L, h): ...@@ -74,7 +74,7 @@ def func_position(x0, v0, f, L, h):
""" """
x1 = x0 + h*v0 + ((h**2)/2)*f # New positions based on current positions, velocities and forces x1 = x0 + h*v0 + ((h**2)/2)*f # New positions based on current positions, velocities and forces
diff = x1 - x0 diff = x1 - x0 # Distance traveled. Used for diffusion
x1 = np.remainder(x1,L) # Account for the periodic boundary conditions x1 = np.remainder(x1,L) # Account for the periodic boundary conditions
# return # return
...@@ -107,7 +107,7 @@ def func_velocity(v0, f0, f1, h): ...@@ -107,7 +107,7 @@ def func_velocity(v0, f0, f1, h):
""" """
v1 = v0 + (h/2)*(f1 + f0) # New velocities based on current velocities and forces and new forces ([N,D] matrix) v1 = v0 + (h/2)*(f1 + f0) # New velocities based on current velocities and forces and new forces ([N,D] matrix)
E_kin_array = .5*(np.linalg.norm(v1, axis=1))**2 # New kinetic energy for every particle E_kin_array = .5*(np.linalg.norm(v1, axis=1))**2 # New kinetic energy for every particle
E_kin = np.sum(E_kin_array) # Total kinetic energy E_kin = np.sum(E_kin_array)
# return # return
return v1, E_kin return v1, E_kin
...@@ -153,9 +153,9 @@ def steps(x0, v0, f0, L, h): ...@@ -153,9 +153,9 @@ def steps(x0, v0, f0, L, h):
D = x0.shape[1] # Number of spatial dimensions D = x0.shape[1] # Number of spatial dimensions
N = x0.shape[0] # Number of particles N = x0.shape[0] # Number of particles
x1, diff = func_position(x0, v0, f0, L, h) # Positions of all particles at time 1 x1, diff = func_position(x0, v0, f0, L, h)
f1, U1, dr1 = func_force(x1, L, D, N) # Forces on all particles at time 1 f1, U1, dr1 = func_force(x1, L, D, N)
v1, E = func_velocity(v0, f0, f1, h) # Velocities (v1) and total kinetic energy (E) of all particles at time 1 v1, E = func_velocity(v0, f0, f1, h)
# return # return
return x1, v1, E, U1, f1, dr1, diff return x1, v1, E, U1, f1, dr1, diff
......
...@@ -61,12 +61,12 @@ def pressure_scalar(D, L, kbT): ...@@ -61,12 +61,12 @@ def pressure_scalar(D, L, kbT):
Returns Returns
------- -------
P Time-averaged pressure (natural units) P : Time-averaged pressure (natural units)
""" """
D += np.eye(D.shape[1]) # Set diagonal elements to 1 D += np.eye(D.shape[1]) # Set diagonal elements to 1
dUdr = D**(-7)-2*D**(-13) dUdr = D**(-7)-2*D**(-13)
dUdr += np.eye(D.shape[1]) # Set diagonal elements to 0 dUdr += np.eye(D.shape[1]) # Set diagonal elements to 0
S = np.sum(24*D*dUdr) # also over time S = np.sum(24*D*dUdr) # Note: also a summation over time
P = 1 - (3*D.shape[1]*kbT)**(-1)*0.5*S/D.shape[0] P = 1 - (3*D.shape[1]*kbT)**(-1)*0.5*S/D.shape[0]
# return # return
...@@ -151,24 +151,24 @@ def pair_correlation(L, N, R): ...@@ -151,24 +151,24 @@ def pair_correlation(L, N, R):
dr = L/M dr = L/M
r = np.linspace(0,L,L/dr) r = np.linspace(0,L,L/dr)
n = np.histogram(R,r)[0] # Make histogram n = np.histogram(R,r)[0] # Make histogram of mutual particle distances
n[0] = 0 # Set value at r=0 to zero to avoid "self-correlation" n[0] = 0 # Set value at r=0 to zero to avoid "self-correlation"
r = np.delete(r,0) # Delete r = 0 r = np.delete(r,0) # r[0] = 0
r -= (r[2]-r[1])/2 # Places points r on the middle of each interval dr r -= (r[2]-r[1])/2 # Places points r on the middle of each interval dr
G = 2*(L**3)/(N*(N-1))*(n/(4*math.pi*r**2*dr)) # Pair correlation function G = 2*(L**3)/(N*(N-1))*(n/(4*math.pi*r**2*dr)) # Pair correlation function
peaks = np.where((G[1:-1] > G[0:-2]) * (G[1:-1] > G[2:]))[0] + 1 # Finds locations of all peaks peaks = np.where((G[1:-1] > G[0:-2]) * (G[1:-1] > G[2:]))[0] + 1 # Finds locations of all peaks
G_peak = G[peaks] # Value of pair correlation at peak G_peak = G[peaks]
r_peak = r[peaks] # Approximate value of mutual distance r at peak r_peak = r[peaks]
A = r_peak*np.ones([r_peak.shape[0], r_peak.shape[0]]) A = r_peak*np.ones([r_peak.shape[0], r_peak.shape[0]])
dA = A - A.transpose() # Difference between all positions of all peaks dA = A - A.transpose() # Difference between all positions of all peaks
dA = np.tril(dA.transpose(), -1).transpose() # Only keep positive values dA = np.tril(dA.transpose(), -1).transpose() # Only keep positive values
dA = dA.reshape(dA.shape[0]*dA.shape[1],) # Change schape into an array dA = dA.reshape(dA.shape[0]*dA.shape[1],) # Change shape into an array
dA = dA[np.where(dA != 0)[0]] # Remove all zeros dA = dA[np.where(dA != 0)[0]]
dA = np.sort(dA) # Sort all values dA = np.sort(dA)
else: else:
dA = 0 dA = 0
G = 0 G = 0
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment