 # 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 # 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,v0.shape)) V = np.zeros((N_steps,x0.shape,v0.shape)) D = np.zeros((N_steps,x0.shape,x0.shape)) ... ...
 # 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) for t in range(X.shape): P[t] = X.shape*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): P[t] = 1 + (V.shape*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) # 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*kbT)**(-1)*0.5*S/D.shape # 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>6: n = 5 else: n = 1 # allocate T = E_kin.shape 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, r_peak.shape]) 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*dA.shape,) # Change schape into an array dA = dA[np.where(dA != 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*dA.shape,) # Change schape into an array dA = dA[np.where(dA != 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.