Commit e4b81861 authored by ggandus's avatar ggandus
Browse files

Merge branch 'master' of gitlab.ethz.ch:danielep/mmm_2019

parents 1ceabedd 67cb68a7
%% Cell type:code id: tags:
``` python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as tck
import scipy
from scipy.special import sph_harm
from ase.io import read
from ase import neighborlist
from ase.visualize import view
from pylab import *
```
%% Cell type:markdown id: tags:
##### plot preferences
%% Cell type:code id: tags:
``` python
SMALL_SIZE = 8 #8
MEDIUM_SIZE = 12 #20
BIGGER_SIZE = 12 #12
plt.rc('font', size=SMALL_SIZE) # controls default text sizes
plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
```
%% Cell type:markdown id: tags:
##### Read files
%% Cell type:code id: tags:
``` python
trajectory = read('trajectory.xyz', index=':')
energy = np.loadtxt('energy.dat')
```
%% Cell type:code id: tags:
``` python
view(trajectory)
```
%% Cell type:markdown id: tags:
##### Compute distances between center of mass and the atom that moves
%% Cell type:code id: tags:
``` python
distance = np.empty(0)
for frame in trajectory:
distance=np.append(distance,np.linalg.norm(frame.get_positions()[3]))
distance/=3.405
```
%% Cell type:markdown id: tags:
##### Compute $Q_l$ with $l=6$
$
q_{lm} = \frac{1}{N_{bonds}} \sum_{N_n} y_{lm} (\theta_{ij}, \phi_{ij})
$
$
Q_l^2 = \frac{4 \pi}{2l+1} \sum_{m=-l}^l |q_{lm}|^2
$
%% Cell type:code id: tags:
``` python
#the desired Ql
l = 6
rcut = 3.405*1.391
#update the cutoff for each frame
for frame in trajectory:
array_rcut = [ rcut for number_of_atoms in frame ]
new_neighbour_list = neighborlist.NeighborList(array_rcut, self_interaction=False, bothways=True)
new_neighbour_list.update(frame)
#compute Ql for each frame
Ql = np.empty(len(trajectory))
i = 0
for frame in trajectory:
nbonds = 0
qlm = np.zeros(2*l+1)
for atom in frame:
nlist = new_neighbour_list.get_neighbors(atom.index)[0]
for theneig in nlist: #cycle over the neighbours
#get angles and distances
nbonds = nbonds+1
rij = frame[theneig].position - atom.position
dist = np.linalg.norm(rij)
phi_ij = np.arccos(rij[2]/dist)
theta_ij = np.arctan2(rij[1],rij[0])
if theta_ij < 0:
theta_ij += 2*np.pi
#move in spherical coordinates space
# In a like-oriented coordinate system at j,
#the spherical coordinates of atom i are:
if theta_ij <= np.pi:
theta_ji = theta_ij + np.pi
elif theta_ij > np.pi:
theta_ji = theta_ij - np.pi
if np.absolute(theta_ji-2*np.pi)<0.0001:
theta_ji=0.0
phi_ji = np.pi-phi_ij
#compute spherical harmonics and perform qml summation
qlm = qlm + np.array([ sph_harm(m,l,theta_ij,phi_ij) for m in range(-l,l+1) ])
qlm = np.real(qlm*np.conj(qlm)/(nbonds*nbonds))
#prefactor and second summation
Ql[i] = np.sqrt(np.pi *4 /(2*l+1)*np.sum(qlm))
i += 1
```
%% Cell type:markdown id: tags:
##### Compute asphericity
%% Cell type:code id: tags:
``` python
I = np.empty([0,3])
for frame in trajectory:
I = np.append(I,[frame.get_moments_of_inertia()],axis=0)
I = I.transpose()
asph = ((I[0]- I[1])**2 + (I[0]- I[2])**2 + (I[2]- I[1])**2 )/ (I[0]**2 + I[1]**2 + I[2]**2)
```
%% Cell type:markdown id: tags:
##### Energy representation
%% Cell type:code id: tags:
``` python
energy = np.loadtxt('energy.dat')
energy = np.transpose(energy)
epsilon = 119*8.616733e-5 #eV
energy[1] *= 27.1442/epsilon
min_e = np.min(energy[1])
energy[1] -= min_e
fig = plt.plot(energy[0],energy[1],color='crimson', lw=2)
#get the axis just defined
ax = plt.gca()
ax.set_ylabel("V/$\epsilon$")
ax.set_xlabel("replicas")
# ---- create a second axis (ax2) ----
# axis ticks
vec_tick =[]
vec_tick_values=[]
for i in range (0,len(energy[0])):
vec_tick.append(i)
vec_tick_values.append(round(Ql[i]*100.,2))
ax2 = ax.twiny()
ax2.set_xticks(vec_tick)
ax2.set_xticklabels(vec_tick_values)
ax2.set_xlabel("$Q_l$ (x 100)")
# ---- create a third axis (ax3) ----
vec_tick =[]
vec_tick_values=[]
for i in range (0,len(energy[0])):
vec_tick.append(i)
vec_tick_values.append(round(asph[i]*100.,2))
ax3 = ax.twiny()
ax3.set_xticks(vec_tick)
ax3.set_xticklabels(vec_tick_values)
ax3.set_xlabel("asphericity (x 100)")
# customising position of this third axes
ax3.xaxis.set_ticks_position("top")
ax3.xaxis.set_label_position("top")