Commit 4dac1845 authored by dmar's avatar dmar
Browse files

add 5th exercise

parent d4cd2666
%% Cell type:code id: tags:
``` python
import numpy as np
from ase.io import read
from ase.visualize import view
import matplotlib.pyplot as plt
```
%% Cell type:code id: tags:
``` python
sigma = 3.405
epsilon = 119.8*8.616733e-5 # eV
```
%% Cell type:code id: tags:
``` python
def g_step(box, n_bin, d):
box += 1
box *= np.sqrt(2)/sigma
del_bin = box/n_bin #width of a bin
g = np.zeros([n_bin,n_bin])
g[0] = np.linspace(0,box,n_bin)
g[0]+= del_bin*0.5
for i in range(len(d)):
for j in range(i+1,len(d)):
g_index = d[i][j]/(del_bin)
g[1,g_index.astype(int)] += 2
# we count both for i and j.
#nomalized by number of particles
g[1] /= len(d[0])
#and by bin volume
for k in range(len(g[1])):
g[1][k] /= ((k+1)**3 -k**3)*del_bin**3
#side of the optimized cell: 5,269
vol = (5.269/sigma)**3
rho = 4./vol
g[1] /= np.pi*(4/3)*rho
return g;
```
%% Cell type:markdown id: tags:
### 38 atom
%% Cell type:code id: tags:
``` python
trajectory = read('production_38-pos-1.xyz', index='::2')
N = len(trajectory[0]) #number of atoms
n_step = len(trajectory) #number of step
#create a distance array
dist = np.empty([0,N,N])
for frame in trajectory:
dist = np.append(dist, [frame.get_all_distances()],axis=0)
#Lennard Jones units
dist /= sigma
```
%% Cell type:code id: tags:
``` python
#understand what the distance array is:
#dist0 = trajectory[0].get_all_distances()
# print(dist0.shape)
# print(dist0)
# print(dist0[0])
```
%% Cell type:code id: tags:
``` python
n_bin = 200 #number of bin in the graph
box = 25 #side of the bon in the simulation
g_38 = np.zeros([n_bin,n_bin])
#average during the thermalized part of the simulation
for i in range(n_step):
g_38 += g_step(box,n_bin,dist[i])
g_38 /= n_step
plt.plot(g_38[0],g_38[1])
plt.xlim(0,box/(2*sigma))
plt.xlabel('distance [ units of $\sigma$]')
plt.ylabel('radial distribution function')
```
%% Cell type:markdown id: tags:
### 150 atoms
%% Cell type:code id: tags:
``` python
trajectory = read('production_150-pos-1.xyz', index='::2')
N = len(trajectory[0]) #number of atoms
n_step = len(trajectory) #number of step
dist = np.empty([0,N,N])
for frame in trajectory:
dist = np.append(dist, [frame.get_all_distances()],axis=0)
dist /= sigma
```
%% Cell type:code id: tags:
``` python
n_bin = 400
box = 30
g_150 = np.zeros([n_bin,n_bin])
for i in range(n_step):
g_150 += g_step(box,n_bin,dist[i])
g_150 /= n_step
plt.plot(g_150[0],g_150[1])
plt.xlim(0,box/(2*sigma))
plt.xlabel('distance [ units of $\sigma$]')
plt.ylabel('radial distribution function')
```
%% Cell type:markdown id: tags:
### 450 atoms
%% Cell type:code id: tags:
``` python
trajectory = read('production_450-pos-1.xyz', index='::2')
N = len(trajectory[0]) #number of atoms
n_step = len(trajectory) #number of step
dist = np.empty([0,N,N])
for frame in trajectory:
dist = np.append(dist, [frame.get_all_distances()],axis=0)
dist /= sigma
```
%% Cell type:code id: tags:
``` python
n_bin = 400
box = 40
g_450 = np.zeros([n_bin,n_bin])
for i in range(n_step):
g_450 += g_step(box,n_bin,dist[i])
g_450 /= n_step
plt.plot(g_450[0],g_450[1])
plt.xlim(0,box/(2*sigma))
plt.xlabel('distance [ units of $\sigma$]')
plt.ylabel('radial distribution function')
```
%% Cell type:markdown id: tags:
### 817 atoms
%% Cell type:code id: tags:
``` python
trajectory = read('production_817-pos-1.xyz', index='::2')
N = len(trajectory[0]) #number of atoms
n_step = len(trajectory) #number of step
dist = np.empty([0,N,N])
for frame in trajectory:
dist = np.append(dist, [frame.get_all_distances()],axis=0)
dist /= sigma
```
%% Cell type:code id: tags:
``` python
n_bin = 400
box = 45
g_817 = np.zeros([n_bin,n_bin])
for i in range(n_step):
g_817 += g_step(box,n_bin,dist[i])
g_817 /= n_step
plt.plot(g_817[0],g_817[1])
plt.xlim(0,box/(2*sigma))
plt.xlabel('distance [ units of $\sigma$]')
plt.ylabel('radial distribution function')
```
%% Cell type:markdown id: tags:
### temperature
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
### Thermalization plots (NVT)
in the .energ file you can find
- column 0: step
- column 1: time
- column 2: kinetic energy
- column 3: temperature
- column 4: potential energy
Try to plot these values and describe what you see and what is conserved
%% Cell type:code id: tags:
``` python
data38 = np.loadtxt('thermalization_38-1.ener' )
data150 = np.loadtxt('thermalization_150-1.ener')
data450 = np.loadtxt('thermalization_450-1.ener')
data817 = np.loadtxt('thermalization_817-1.ener')
data38 = np.transpose(data38)
data150 = np.transpose(data150)
data450 = np.transpose(data450)
data817 = np.transpose(data817)
plt.plot(data38[1],data38[3], label='38')
plt.plot(data150[1],data150[3], label='150')
plt.plot(data450[1],data450[3], label='450')
plt.plot(data817[1],data817[3], label='817')
plt.legend( loc='upper right')
plt.xlabel(' time [fs]')
plt.ylabel('T [K]')
```
%% Cell type:markdown id: tags:
### Production plots (NVE)
in the .energ file you can find
- column 0: step
- column 1: time
- column 2: kinetic energy
- column 3: temperature
- column 4: potential energy
Try to plot these values and describe what you see and what is conserved
%% Cell type:code id: tags:
``` python
data38 = np.loadtxt('production_38-1.ener' )
data150 = np.loadtxt('production_150-1.ener')
data450 = np.loadtxt('production_450-1.ener')
data817 = np.loadtxt('production_817-1.ener')
data38 = np.transpose(data38)
data150 = np.transpose(data150)
data450 = np.transpose(data450)
data817 = np.transpose(data817)
plt.plot(data38[1],data38[3], label='38')
plt.plot(data150[1],data150[3], label='150')
plt.plot(data450[1],data450[3], label='450')
plt.plot(data817[1],data817[3], label='817')
plt.legend( loc='upper right')
plt.xlabel(' time [fs]')
plt.ylabel('T [K]')
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment