...
 
from . import initial_state
from . import convert
from . import graphics
\ No newline at end of file
from . import parameters
\ No newline at end of file
# import modules
import numpy as np
import math
import random
# import functions
from comp.movements import normalize
def build(N, L, T_init, topology, org='fcc'):
def build(
N, L, T_init,
topology, org='fcc'):
"""
Sets initial particle positions and velocities.
Input
N Number of particles inside the particle box
L Size of the particle box
T_init kbT/epsilon, initial temperature
topology Number of dimensions
org Way to spatially organize the particles, either 'fcc' or 'cubic'
Output
x0 Particle positions ((N, topology) np.array)
v0 Particle velocities ((N, topology) np.array)
density Particle density after initialization
Parameters
----------
N : float
Number of particles inside the particle box
L : float
Size of the particle box
T_init : float
kbT/epsilon, initial temperature
topology : float
Number of dimensions
org : string
Way to spatially organize the particles, either 'fcc' or 'cubic'
Returns
-------
x0 : (N, topology) np.array
Particle positions
v0 : (N, topology) np.array
Particle velocities
density : float
Particle density after initialization, as fcc(N,L) can change the number of particles
Raise
-----
IOError
org has not been given a correct string input
"""
if org == "fcc": # Initialize particles in an fcc lattice
# initial positions
if org == "fcc":
# Initialize particles in an fcc lattice and calculate new density
N, x0 = fcc(N, L)
density = N/(L**3) # calculate new density, as fcc(N,L) can change number of particles
elif org == "cubic": # Initialize particles in a cubic lattice, randomly distributed
density = N/(L**3)
elif org == "cubic":
# Initialize particles in a cubic lattice, randomly distributed
x0 = cubic(N, topology, L, T_init)
density = N/(L**3)
else:
raise IOError("org has not been given a correct string input, should be either 'fcc' or 'cubic'")
# initial velocities
v0 = velocity(T_init, N, topology)
# return
return x0, v0, density
......@@ -39,60 +65,66 @@ def cubic(
temperature):
"""
Initializes particles in a cubic lattice, randomly distributed.
input
N_particles # [-] number of Argonatoms
topology # [-] Number of dimensions
L_box # [-] length of box; size to be determined
N_grid # [-] number of gridpoints in one direction
temperature # [-] Initial temperature
A box of size (L_box,L_box,L_box (in 3D)) is divided in N_grid**3 gridpoints.
All this gridpoints have an internal distance d = L_box/N_grid, measured in orthogonal directions.
These gridpoints have positions n_x = 0, 1, 2, ... . Similar for n_y and n_z.
Using random.sample() every particle is assigned a unique gridpoint.
Subsequently these gridpoints are translated to their corresponding x, y and z coordinates.
Parameters
----------
N_particles : float
number of Argonatoms
topology: float
Number of dimensions
L_box: float
length of box; size to be determined
N_grid: float
number of gridpoints in one direction
temperature: float
Initial temperature
Returns
-------
x: (N_particles, topology) array
x, y, and z positions for every particle
output
x # [-] (N_particles, topology) array with for every particle the x, y, and z positions
description
A box of size (L_box,L_box,L_box (in 3D)) is divided in N_grid**3 gridpoints.
All this gridpoints have an internal distance d = L_box/N_grid, measured in orthogonal directions.
These gridpoints have positions n_x = 0, 1, 2, ... . Similar for n_y and n_z.
Using random.sample() every particle is assigned a unique gridpoint.
Subsequently these gridpoints are translated to their corresponding x, y and z coordinates.
Raises
------
ValueError
When the number of particles is larger than the number of gridpoints.
ValueError
When the topology is not equal to 1, 2 or 3.
"""
# import modules
import numpy as np
import math
import random
# import functions
from comp.movements import normalize
# Position
N_grid = math.ceil(N_particles**(1/3))
d = L_box/N_grid # [-] distance between two gridpoints
d = L_box/N_grid
# # create random sample
# Create random sample
if N_particles > (N_grid**topology):
raise ValueError(f"N_particles ({N_particles}) is larger than number of gridpoints ({N_grid}**{topology}={N_grid**topology})")
else:
# (N_particles, 1) array with random sample of unique position numbers
n_positions = np.array(random.sample((range(0,N_grid**topology)),N_particles))
# # select position number
if topology == 3: # 3D
n_x = (n_positions-(n_positions % N_grid**2))/N_grid**2 # x position number
n_y = n_positions % N_grid # y position number
n_z = np.floor(n_positions / N_grid) % N_grid # z position number
n = np.vstack([n_x, n_y, n_z]) # (topology, N_particles) array with position numbers
elif topology == 2: #2D
n_x = (n_positions-(n_positions % N_grid))/N_grid # x position number
n_y = n_positions % N_grid # y position number
n = np.vstack([n_x, n_y]) # (topology, N_particles) array with position numbers
elif topology == 1: #1D
n_x = n_positions # x position number
n = np.vstack([n_x]) # (topology, N_particles) array with position numbers
# Select position number
if topology == 3:
n_x = (n_positions-(n_positions % N_grid**2))/N_grid**2
n_y = n_positions % N_grid
n_z = np.floor(n_positions / N_grid) % N_grid
n = np.vstack([n_x, n_y, n_z])
elif topology == 2:
n_x = (n_positions-(n_positions % N_grid))/N_grid
n_y = n_positions % N_grid
n = np.vstack([n_x, n_y])
elif topology == 1:
n_x = n_positions
n = np.vstack([n_x])
else:
raise ValueError("Topology is not equal to 1, 2, or 3")
# # generate position coordinates
x = np.transpose(n*d) # (N_particles, topology) array with for every particle a unique x, y and z position
# Generate position coordinates
x = np.transpose(n*d)
# return
return x
......@@ -101,15 +133,21 @@ def cubic(
def fcc(N,L):
"""
Function to create initial position of N particles in a box of (L,L,L), ordered in a FCC-lattice.
last modified: 2019
last modified: 2019/03/26
input
N: number of particles (scalar)
L: size of the box (scalar)
Parameters
----------
N : float
number of particles
L : float
size of the box
output
N: number of particles that is in accordance with a FCC-latice (scalar)
pos_initial: (x,y,z) positions of particles (N,[x,y,z])-array
Returns
-------
N: float
number of particles that is in accordance with a FCC-latice
pos_initial : (N,[x,y,z])-array
(x,y,z) positions of particles
"""
......@@ -137,33 +175,40 @@ def fcc(N,L):
# positions for N particles: (N,[x,y,z])-array
pos_initial = np.concatenate((pos_basis, pos_1, pos_2, pos_3),axis=0)
# returns
return N, pos_initial
def velocity(T, N, topology):
"""
function to compute initial velocity for N particles in (x,y,z) direction
Function to compute initial velocity for N particles in (x,y,z) direction.
last modified: 2019-03-15
input
T: temperature (unitless) (scalar)
N: number of particles (scalar)
topology: number of dimensions
Parameters
----------
T : float
temperature
N : float
number of particles
topology : float
number of dimensions
output
v: initial velocity (N,[x,y,z])-array
Returns
-------
v : (N,[x,y,z])-array
initial velocity
"""
# MAGNITUDE
# Magnitude
v_magnitude = np.random.normal(0, math.sqrt(2*T), N)
# DIRECTION
# Direction
v_direction = np.random.rand(N,topology)
# # normalize direction
v_length = np.linalg.norm(v_direction,axis=1)
v_normalized = v_direction / v_length[:,None]
# VELOCITY (N,[x,y,z])-array
# Velocity
v = np.multiply(v_normalized, v_magnitude[:, np.newaxis])
# return
......
"""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)
"""Semi-tunable model parameterss"""
topology = 3 # Number of dimensions
error_alpha = 0.2 # Transparancy of errorbar shade
This source diff could not be displayed because it is too large. You can view the blob instead.