In [38]:
import numpy as np
from boxvectors import directions as directions
import Initial_Parameters as ip
from md import System
from md import md
from distribution import maxwellboltzmann
import matplotlib.pyplot as plt
from scipy.special import erf
from scipy.special import erfc
from scipy.constants import epsilon_0
%matplotlib inline
In [2]:
Symbols = ip.Symbols
Coefficients = ip.Coefficients
Charges = ip.Charges
N = ip.N*np.sum(Coefficients)
L = ip.L
T = ip.T
dt = ip.dt
p_rea = ip.p_rea
p = ip.p
std = ip.std
k_cut = ip.k_cut
k_max = ip.k_max_long_range
n_boxes_LJ = ip.n_boxes_LJ
n_boxes_coulomb = ip.n_boxes_short_range
Sys= System(Symbols, Coefficients, Charges, N/2)
Labels = Sys.get_Labels()
Sigma, Epsilon = Sys.get_LJ_parameter()
m = Labels[:,0]
In [3]:
switch_parameter = np.array([1,-1,0,0])
r_switch = 1.5
In [4]:
def get_random_starting_Positions(N,L):
Positions = np.zeros((N,3))
Positions[:,0] = np.linspace(0.1,L[0],N, endpoint = False)
Positions[:,1] = np.linspace(0.1,L[1],N, endpoint = False)
Positions[:,2] = np.linspace(0.1,L[2],N, endpoint = False)
np.random.shuffle(Positions[:,0])
np.random.shuffle(Positions[:,1])
np.random.shuffle(Positions[:,2])
return Positions
Positions = get_random_starting_Positions(N,L)
Velocities = maxwellboltzmann().sample_distribution(N,m,T)
Forces = np.zeros((N,3))
R = np.linalg.norm(Positions,axis=1)
In [5]:
MD = md(
Positions,
R,
Labels,
Velocities,
Forces,
L,
T,
std,
Sigma,
Epsilon,
switch_parameter,
r_switch,
n_boxes_coulomb,
k_max,
dt,
p_rea,
k_cut)
In [6]:
MD.forces = MD.get_forces()
In [7]:
Temperature = np.zeros(10)
for i in np.arange(10):
Positions_New, Velocities_New, Forces_New = MD.propagte_system()
MD.positions = Positions_New
MD.velocities = Velocities_New
MD.forces = Forces_New
Temperature[i] = MD.get_Temperature()
plt.plot(Temperature)
Out[7]:
In [19]:
Temperature
Out[19]:
In [17]:
Velocities = maxwellboltzmann().sample_distribution(N,m,T)
Velocities
MD.velocities
Out[17]:
In [20]:
from particle_interaction import lennard_jones
In [35]:
def __short_range_forces(self,Positions,R, Labels,L):
''' Calculate the Force resulting from the short range coulomb interaction between the Particles
Parameters
---------------
Positions: Nx3 Array
Array with N rows and 3 Columns, each row i contains the x,y and z coordinates of particle i.
R: Nx1 Array
Array with N entries. Contains each Particles Distance to the coordinate origin.
Labels: Nx? Array
Array with N rows and ? Columns. The third Column should contain labels, that specify the chemical species of the Particles.
Particle A should have the label 1 and Particle B should have the label 0.
L:3x1 Array
Array containg the Dimensions of the Simulation box
Returns
--------------
Force_short_range: Nx3 Array
Array with N rows and 3 Columns. Each Row i contains the short-range-Force acting upon Particle i componentwise.
'''
K = self.LinComb(self.n_boxes_short_range)
K[:,0] *=L[0]
K[:,1] *=L[1]
K[:,2] *=L[2]
N = np.size(Positions[:,0])
Force_short_range = np.zeros((N,3))
k = np.size(K[:,0])
charges_pre_delete = np.zeros((N,3))
charges_pre_delete[:,0]=Labels[:,1]
charges_pre_delete[:,1]=Labels[:,1]
charges_pre_delete[:,2]=Labels[:,1]
for i in np.arange(N):
dists_single_cell_pre_delete = Positions[i,:]-Positions
dists_single_cell = np.delete(dists_single_cell_pre_delete,i,0)
charges = np.delete(charges_pre_delete,i,0)
#This Skript paralellizes the sum over j and executes the sum over vectors n within the loop
for j in np.arange(k):
dists = dists_single_cell+K[j]
norm_dists = np.linalg.norm(dists)
Force_short_range[i,:] += np.sum(charges*dists/norm_dists**2
*( erfc( dists/np.sqrt(2)/self.std )/norm_dists
+ np.sqrt(2.0/np.pi)/self.std*np.exp(-norm_dists**2/(2*self.std**2) )) ,0)
#Getting the Pre-factor right
Force_short_range = Force_short_range* charges_pre_delete /(8*np.pi*epsilon_0)
return Force_short_range
def __long_range_forces(self,Positions,Labels,L):
''' Calculate the Force resulting from the long range Coulomb interaction between the Particles
Parameters
---------------
Positions: Nx3 Array
Array with N rows and 3 Columns, each row i contains the x,y and z coordinates of particle i.
Labels: Nx? Array
Array with N rows and ? Columns. The third Column should contain labels, that specify the chemical species of the Particles.
Particle A should have the label 1 and Particle B should have the label 0.
L:3x1 Array
Array containg the Dimensions of the Simulation box
Returns
--------------
Force_long_range: Nx3 Array
Array with N rows and 3 Columns. Each Row i contains the long-range-Force acting upon Particle i componentwise.
'''
i = np.complex(0,1)
# setup k-vector matrix
k = self.LinComb(self.k_max_long_range)*(2*np.pi)/L[0]
# two k-vectors i,j have property k_i = -k_j respectively, delete one of them and delete k = (0,0,0)
k = np.delete(k, (np.arange((np.shape(k)[0]+1)/2)), axis=0)
# delete all k-vectors that are longer than cutoff
k = np.delete(k, (np.where(np.sqrt(sum(np.transpose(k**2))) > self.k_cut)[0]), axis=0)
# setup data needed for calculation
k_betqua = sum(np.transpose(k**2)) # |k|^2
num_k = np.shape(k_betqua)[0] # number of k-vectors
num_Par = np.shape(Positions)[0] # number of particles
f_1 = np.zeros((num_k,3)) # placeholder for Forces during summation over k outside of for second loop
f_2 = np.zeros((num_k)) # placeholder for Forces during summation over k in second for-loop
Force_long_range = np.zeros((num_Par,3)) # placeholder for long ranged forces
charges = np.zeros((num_Par,3)) # create num_Par x 3 matrix with charges
charges[:,0]=Labels[:,1] # in each column L[:,1] repeated
charges[:,1]=Labels[:,1]
charges[:,2]=Labels[:,1]
for h in np.arange(num_Par): # loop over all particles
for j in np.arange(num_k): # loop over all k-vectors
# "structure factor" sum (right sum in equation)
f_2[j] = sum(Labels[:,1]*np.sin(np.dot(k[j,:],(Positions[h,:]*np.ones((num_Par,3)) - Positions).transpose())))
# complete sum over all k, (left part of equation)
f_1 = (((np.exp(-self.std**2/2*k_betqua)/k_betqua)*f_2)*np.ones((3,num_k))).transpose()*k
# actually sum over all k
Force_long_range[h,:] = sum(f_1)
# get prefactor right
Force_long_range *= charges/(L[0]**3*epsilon_0)*2 # multiply by 2 because of symmetry properties of sine
return Force_long_range
In [47]:
print (__short_range_forces(MD.coulomb, Positions, R, Labels, L))
print (__long_range_forces(MD.coulomb, Positions, Labels, L))
print np.average(__long_range_forces(MD.coulomb, Positions, Labels, L)+__short_range_forces(MD.coulomb, Positions, R, Labels, L))
In [75]:
K = MD.coulomb.LinComb(MD.coulomb.n_boxes_short_range)
K[:,0] *=L[0]
K[:,1] *=L[1]
K[:,2] *=L[2]
N = np.size(Positions[:,0])
Force_short_range = np.zeros((N,3))
k = np.size(K[:,0])
charges_pre_delete = np.zeros((N,3))
charges_pre_delete[:,0]=Labels[:,1]
charges_pre_delete[:,1]=Labels[:,1]
charges_pre_delete[:,2]=Labels[:,1]
for i in np.arange(N):
dists_single_cell_pre_delete = Positions[i,:]-Positions
dists_single_cell = np.delete(dists_single_cell_pre_delete,i,0)
charges = np.delete(charges_pre_delete,i,0)
#This Skript paralellizes the sum over j and executes the sum over vectors n within the loop
for j in np.arange(k):
dists = dists_single_cell+K[j]
norm_dists = np.linalg.norm(dists)
Force_short_range[i,:] += np.sum(charges*dists/norm_dists**2
*( erfc( dists/np.sqrt(2)/std )/norm_dists
+ np.sqrt(2.0/np.pi)/std * np.exp(-norm_dists**2/(2*std**2) )) ,0)
#Getting the Pre-factor right
Force_short_range = Force_short_range* charges_pre_delete /(8*np.pi*epsilon_0)
In [76]:
charges_pre_delete /(8*np.pi*epsilon_0)
Out[76]: