# import modules
import numpy as np
def normalize(matrix_input,axis_input):
"""
function to normalize a vector
last modified: 2019-02-26
input
matrix_input: matrix to be normalized. (m x n) numpy array
axis_input: axis along which the matrix M will be normalized.
1 normalize along horizontal axis
2 normalize along vertical axis
output
n: matrix M normalized along axisdef normalize(M, axis):
"""
matrix_squared = np.square(matrix_input)
vector_sum = np.sum(matrix_squared,axis=axis_input)
vector_length = vector_sum**(1/2)
matrix_normalized = (matrix_input / vector_length[:,None])
return matrix_normalized
def func_force(X, L, D, N):
"""
Function to determine the total forces on all particles
Function to determine the total forces on all particles.
last modified: 20-03-2019
input
X: positions of all particles ([N,D] matrix)
L: size of the box (scalar)
D: number of spatial dimensions (scalar)
N: number of particles (scalar)
Parameters
----------
X : [N,D] matrix
positions of all particles
L : float
size of the box
D : float
number of spatial dimensions
N : float
number of particles
output
F: total forces on all particles ([N,D] matrix)
U: Total potential energy of the system (scalar)
dr_mag-np.eye(N): all mutual distances between all particles([N,N] matrix)
Returns
-------
F : [N,D] matrix
total forces on all particles
U : float
Total potential energy of the system
dr_mag-np.eye(N) : [N,N] matrix
all mutual distances between all particles
"""
x = X*np.ones([N,N,D]) # Positions of all particles ([N,N,D] array)
......@@ -54,78 +42,112 @@ def func_force(X, L, D, N):
F = np.sum(F, axis = 2) # All forces on all particles ([D,N] matrix)
F = F.transpose() # Transpose matrix to obtain desired [N,D] matrix
# return
return F, U, dr_mag-np.eye(N)
def func_position(x0, v0, f, L, h):
"""
Function to calculate the new positions of N particles based on the old positions and velocities
Function to calculate the new positions of N particles based on the old positions and velocities.
last modified: 22-03-2019
input
x0: positions of all particles at time 0 ([N,D] matrix)
v0: velocities of all particles at time 0 ([N,D] matrix)
f: forces on all particles at time 0 ([N,D] matrix)
L: size of the box (scalar)
h: timestep (scalar)
diff: differences between particle positions at times 0 and 1 ([N,D] matrix)
Parameters
----------
x0 : [N,D] matrix
positions of all particles at time 0
v0 : [N,D] matrix
velocities of all particles at time 0
f : [N,D] matrix
forces on all particles at time 0
L : float
size of the box
h : float
timestep
diff : [N,D] matrix
differences between particle positions at times 0 and 1
output
x1: positions of all particles at time 1 ([N,D] matrix)
Returns
-------
x1 : [N,D] matrix
positions of all particles at time 1
"""
x1 = x0 + h*v0 + ((h**2)/2)*f # New positions based on current positions, velocities and forces
diff = x1 - x0
x1 = np.remainder(x1,L) # Account for the periodic boundary conditions
# return
return x1, diff
def func_velocity(v0, f0, f1, h):
"""
Function to calculate the new velocities of N particles based on the old velocities and forces and the new forces
Function to calculate the new velocities of N particles based on the old velocities and forces and the new forces.
last modified: 05-03-2019
input
v0: velocities of all particles at time 0 ([N,D] matrix)
f0: forces on all particles at time 0 ([N,D] matrix)
f1: forces on all particles at time 1 ([N,D] martix)
h: timestep (scalar)
Parameters
----------
v0 : [N,D] matrix
velocities of all particles at time 0
f0 : [N,D] matrix
forces on all particles at time 0
f1 : [N,D] matrix
forces on all particles at time 1
h : float
timestep
output
v1: velocities of all particles at time 1 ([N,D] matrix)
E_kin: total kinetic energy of all particles at time 1 (scalar)
Returns
-------
v1 : [N,D] matrix
velocities of all particles at time 1
E_kin : float
total kinetic energy of all particles at time 1
"""
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 = np.sum(E_kin_array) # Total kinetic energy
# return
return v1, E_kin
def steps(x0, v0, f0, L, h):
"""
Function to calculate the positions and velocities of all particles at time 1 based on the positions and velocities of all particles at time 0
Function to calculate the positions and velocities of all particles at time 1 based on the positions and velocities of all particles at time 0.
last modified: 22-03-2019
input
x0: positions of all particles at time 0 ([N,D] matrix)
v0: velocities of all particles at time 0 ([N,D] matrix)
f: forces on all particles at time 0 ([N,D] matrix)
L: size of the box (scalar)
h: timestep (scalar)
output
x1: positions of all particles at time 1 ([N,D] matrix)
v1: velocities of all particles at time 1 ([N,D] matrix)
E: total kinetic energy of all particles at time 1 (float)
U1: total potential energy at time 1 (float)
f1: total forces on all particles at time 1 ([N,D] matrix)
dr1: all mutial distances between particles ([N,N] matrix)
diff: differences between particle positions at times 0 and 1 ([N,D] matrix)
Parameters
----------
x0 : [N,D] matrix
positions of all particles at time 0
v0 : [N,D] matrix
velocities of all particles at time 0
f : [N,D] matrix
forces on all particles at time 0
L : float
size of the box
h : float
timestep
Returns
-------
x1 : [N,D] matrix
positions of all particles at time 1
v1 : [N,D] matrix
velocities of all particles at time 1
E : float
total kinetic energy of all particles at time 1
U1 : float
total potential energy at time 1
f1 : [N,D] matrix
total forces on all particles at time 1
dr1 : [N,N] matrix
all mutial distances between particles
diff : [N,D] matrix
differences between particle positions at times 0 and 1
"""
D = x0.shape[1] # Number of spatial dimensions
......@@ -134,7 +156,8 @@ def steps(x0, v0, f0, L, h):
x1, diff = func_position(x0, v0, f0, L, h) # Positions of all particles at time 1
f1, U1, dr1 = func_force(x1, L, D, N) # Forces on all particles at time 1
v1, E = func_velocity(v0, f0, f1, h) # Velocities (v1) and total kinetic energy (E) of all particles at time 1
# return
return x1, v1, E, U1, f1, dr1, diff
def scalefactor(N, kbT, E_kin_nat, E_kin_selection):
......@@ -142,14 +165,19 @@ def scalefactor(N, kbT, E_kin_nat, E_kin_selection):
function to generate a scalefactor when imposing a different temperature. Works for both physical and natural units
last modified: 2019-03-05
input
N: [-] number of particles (scalar)
kbT: [-] or [J = m^2 kg s^-2] boltzmann constant times imposed temperature (scalar)
E_kin: [-] or [J] (scalar) with computed kinetic energy
Parameters
----------
N : float
number of particles
kbT : float
boltzmann constant times imposed temperature, [-] or [J = m^2 kg s^-2]
E_kin : float
with computed kinetic energy, [-] or [J]
output
scalefactor: [-] (scalar) lambda: velocity scalefactor
Returns
-------
scalefactor : float
lambda: velocity scalefactor
"""
# average energy
......@@ -161,32 +189,54 @@ def scalefactor(N, kbT, E_kin_nat, E_kin_selection):
# return
return(scalefactor)
def assemble(x0, v0, N_steps, N_argon, L, h, h_settemp, temp, temp_final, temp_delta, E_kin_selection):
def assemble(x0, v0,
N_steps, N_argon, L,
h, h_settemp,
temp, temp_final, temp_delta,
E_kin_selection):
"""
Assembles two (N_steps, N_argon, topology) arrays with the time evolution of the particles' positions and velocities.
Here N_steps corresponds to the number of time steps to be considered.
Both input and output in natural units.
input
x0: positions of all N particles at the initial time ([N,D] matrix)
v0: velocities of all N particles at the initial time ([N,D] matrix)
N_steps: number of timesteps
L: size of the box (scalar)
h: timestep (scalar)
h_settemp: number of timesteps after which the temperature is scaled (scalar)
temp_current: current temperature (scalar)
temp_final: goal temperature (scalar)
temp_delta: stepsize temperature at h_settemp (scalar)
Parameters
----------
x0 : [N,D] matrix
positions of all N particles at the initial time
v0 : [N,D] matrix
velocities of all N particles at the initial time
N_steps : float
number of timesteps
L : float
size of the box
h : float
timestep
h_settemp : float
number of timesteps after which the temperature is scaled
temp_current : float
current temperature
temp_final : float
goal temperature
temp_delta : float
stepsize temperature at h_settemp
output
X: (N_steps, N_argon, topology) array with the time evolution of the particle positions
V: (N_steps, N_argon, topology) array with the time evolution of the particle velocities
E_kin: total kinetic energy of the system for every step ([1,N_steps - 1] array)
U: total potential energy of the system for every step ([1,N_steps - 1] array)
D: (N_steps, N_argon, topology) array with the time evolution of the mutual distances between particles
diffusion: averaged diffusion of all particles ([1,N_steps] array)
Returns
-------
X : (N_steps, N_argon, topology) array
with the time evolution of the particle positions
V : (N_steps, N_argon, topology) array
with the time evolution of the particle velocities
E_kin : ([1,N_steps - 1] array)
total kinetic energy of the system for every step
U : ([1,N_steps - 1] array)
total potential energy of the system for every step
D : (N_steps, N_argon, topology) array
with the time evolution of the mutual distances between particles
diffusion : ([1,N_steps] array)
averaged diffusion of all particles
"""
# allocate
X = np.zeros((N_steps,x0.shape[0],v0.shape[1]))
V = np.zeros((N_steps,x0.shape[0],v0.shape[1]))
D = np.zeros((N_steps,x0.shape[0],x0.shape[0]))
......
# import modules
import numpy as np
import math
......@@ -6,22 +7,29 @@ def pressure(X, F, L, kbT):
"""
Calculates the pressure inside the particle box at each time step. Formula from Exercize 8.3.c from the book "Molecular Dynamics" (see background reading at GitLab).
Last modified: 6-3-2019
Status: does not work, returns also negative and very fluctuating values because of periodic boundary conditions.
input:
X Matrix with the position of each particle at each time step ((N_time_steps, N_particles, N_dimensions) np.array, natural units)
F Matrix with the force at each particle at each time step ((N_time_steps, N_particles, N_dimensions) np.array, natural units)
L Length of the particle box (natural units)
kbT Energy, kb*T/epsilon (natural units)
output:
P Array with the pressure inside the particle box at each time step ((N_time_steps) np.array)
Parameters
----------
X : ((N_time_steps, N_particles, N_dimensions) np.array, natural units)
Matrix with the position of each particle at each time step
F : ((N_time_steps, N_particles, N_dimensions) np.array, natural units)
Matrix with the force at each particle at each time step
L : scalar
Length of the particle box (natural units)
kbT : scalar
Energy, kb*T/epsilon (natural units)
Returns
-------
P : ((N_time_steps) np.array)
Array with the pressure inside the particle box at each time step
"""
P = np.zeros(X.shape[0])
for t in range(X.shape[0]):
P[t] = X.shape[1]*kbT/(L**3) + 1/(3*L**3)*np.sum(np.sum(X[t]*F[t], axis=1))
# return
return P
......@@ -29,14 +37,21 @@ def pressure_coll(D, V, L, kbT, h):
"""
Calculates the pressure inside the particle box at each time step, based on collisions. Formula from Exercize 8.3.c from the book "Molecular Dynamics" (see background reading at GitLab).
input:
D Matrix with the mutual distances between particles at each time step ((N_time_steps, N_particles, N_particles) np.array, natural units)
V: Matrix with the time evolution of the particle velocities ((N_steps, N_argon, topology) array, natural units)
L Length of the particle box (natural units)
kbT Energy, kb*T/epsilon (natural units)
output:
P Array with the pressure inside the particle box at each time step ((N_time_steps) np.array)
Parameters
----------
D : ((N_time_steps, N_particles, N_particles) np.array, natural units)
Matrix with the mutual distances between particles at each time step
V : ((N_steps, N_argon, topology) array, natural units)
Matrix with the time evolution of the particle velocities
L : float
Length of the particle box (natural units)
kbT : float
Energy, kb*T/epsilon (natural units)
Returns
-------
P : ((N_time_steps) np.array)
Array with the pressure inside the particle box at each time step
"""
r_coll = 1 # radius sigma from particle position defines edge hard core repulsion
......@@ -57,6 +72,7 @@ def pressure_coll(D, V, L, kbT, h):
for t in range(D_coll.shape[0]):
P[t] = 1 + (V.shape[1]*v_av[t]*h)**(-1)*np.sum(D_coll[t]*V_diff[t])
# return
return P
......@@ -64,12 +80,17 @@ def pressure_scalar(D, L, kbT):
"""
Calculates time-averaged pressure inside the particle box.
Input
D Matrix with the mutual distances between particles at each time step ((N_time_steps, N_particles, N_particles) np.array, natural units)
L Length of the particle box (natural units)
kbT Energy, kb*T/epsilon (natural units)
Output
Parameters
----------
D : ((N_time_steps, N_particles, N_particles) np.array, natural units)
Matrix with the mutual distances between particles at each time step
L : float
Length of the particle box (natural units)
kbT : float
Energy, kb*T/epsilon (natural units)
Returns
-------
P Time-averaged pressure (natural units)
"""
D += np.eye(D.shape[1]) # Set diagonal elements to 1
......@@ -78,6 +99,7 @@ def pressure_scalar(D, L, kbT):
S = np.sum(24*D*dUdr) # also over time
P = 1 - (3*D.shape[1]*kbT)**(-1)*0.5*S/D.shape[0]
# return
return P
......@@ -86,29 +108,40 @@ def specific_heat(E_kin, N):
function to calculate the specific heat
last modified: 22-03-2019
input
E_kin: array of total kinetic energy for every timestep
N: number of argon particles
Parameters
----------
E_kin : array
array of total kinetic energy for every timestep
N : float
number of argon particles
output
Cv: specific heat
n: number of timesteps over which the energies are averaged
Returns
-------
Cv : array
specific heat
n : float
number of timesteps over which the energies are averaged
"""
# set number of timesteps over which the energies are averaged
if E_kin.shape[0]>6:
n = 5
else:
n = 1
# allocate
T = E_kin.shape[0]
dE_mean = ([])
E_mean = ([])
# calculate average energy
for t in range(n, T):
E_mean = np.append(E_mean, np.mean(E_kin[t-n:t])**2)
dE_mean = np.append(dE_mean, np.mean(E_kin[t-n:t]**2) - np.mean(E_kin[t-n:t])**2)
# calculate specific heat
Cv = (2/(3*N) - dE_mean/E_mean)**(-1)
# return
return Cv, n
......@@ -117,21 +150,34 @@ def pair_correlation(L, N, R):
function to calculate the pair correlation function
last modified: 19-03-2019
input
L: Size of the box
N: Number of argon particles
R: matrix of all distances between particles at all times
Parameters
----------
L : float
Size of the box
N : float
Number of argon particles
R : [N,N,timesteps]
matrix of all distances between particles at all times
output
G: Pair correlation function
r: Distances for which the pair correlation function is defined
dA: All distances between all peaks in the pair correlation function
Returns
-------
G : array
Pair correlation function
r : array
Distances for which the pair correlation function is defined
dA : array
All distances between all peaks in the pair correlation function
Raises
------
print("Not enough particles to determine a pair correlation function")
When not enough particles to determine a pair correlation function (obviously)
"""
if N>4:
M = N # Number of points in linspace
M = N # Number of points in linspace
R = np.mean(R, axis=0) # Mean of all mutual distances
R = np.mean(R, axis=0) # Mean of all mutual distances
dr = L/M
r = np.linspace(0,L,L/dr)
......@@ -148,17 +194,18 @@ def pair_correlation(L, N, R):
r_peak = r[peaks] # Approximate value of mutual distance r at peak
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 = 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[np.where(dA != 0)[0]] # Remove all zeros
dA = np.sort(dA) # Sort all values
dA = A - A.transpose() # Difference between all positions of all peaks
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[np.where(dA != 0)[0]] # Remove all zeros
dA = np.sort(dA) # Sort all values
else:
dA = 0
G = 0
r = 0
print("Not enough particles to determine a pair correlation function")
# return
return G,r, dA
......@@ -167,16 +214,22 @@ def specific_heat_total(E_kin, N):
function to calculate the specific heat averaged over total time
last modified: 22-03-2019
input
E_kin: Kinetic energy of the system at all times
N: Number of argon particles
Parameters
----------
E_kin : array
Kinetic energy of the system at all times
N : float
Number of argon particles
output
Cv: Heat capacity, calculating using averages over total simulated time
Returns
-------
Cv : float
Heat capacity, calculated using averages over total simulated time
"""
E_mean = np.mean(E_kin)**2
dE_mean = np.mean(E_kin**2) - E_mean
Cv = (2/(3*N) - dE_mean/E_mean)**(-1)
# return
return Cv
\ No newline at end of file
"""Non-tunable parameters for conversion between physical and natural units"""
sigma = 3.405E-10 # Natural scalefactor for distance (meters)
tau = 2.15E-12 # Natural scalefactor for time (seconds)
kb = 1.381E-23 # Boltzmann constant (Joule/Kelvin)
T_ref = 119.8 # Reference temperature (Kelvin)
mass = 6.6335209*10**-26 # Mass of the argon atom (kg)
epsilon = T_ref*kb # Natural units for energy (Joule)
topology = 3 # Number of dimensions
\ No newline at end of file
"""Tunable parameters for the looks of the graphs"""
error_alpha = 0.2 # transparancy of errorbar shade
\ No newline at end of file
......@@ -3,9 +3,6 @@ import numpy as np
import math
import random
# import functions
from comp.movements import normalize
def build(
N, L, T_init,
topology, org='fcc'):
......
This source diff could not be displayed because it is too large. You can view the blob instead.