...
 
This diff is collapsed.
# 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.