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

Exercises/Exercise-10/.ipynb_checkpoints/Exercise-10-2-checkpoint.ipynb

+{
+ "cells": [],
+ "metadata": {},
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

Exercises/Exercise-10/.ipynb_checkpoints/Exercise-10-checkpoint.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import itertools
+import numpy as np
+from numpy.linalg import norm
+import spglib
+from glob import glob
+
+from ase import Atoms
+from ase.io import read
+from ase.spacegroup import crystal
+from ase.spacegroup import Spacegroup
+from ase.data import atomic_numbers, atomic_names
+import nglview
+```
+
+%% Cell type:markdown id: tags:
+
+## Crystallographic point groups
+
+%% Cell type:code id: tags:
+
+``` python
+a=5.117
+b=5.175
+c=5.291
+alpha=90
+beta=99.22
+gamma=90
+thespacegroup=14
+print(Spacegroup(thespacegroup))
+
+hfo2 = crystal(symbols=['Hf','O','O'],
+               basis=[(0.276,0.04,0.208),(0.074,0.332,0.347),(0.449,0.758,0.480)],
+               spacegroup=thespacegroup, 
+               cellpar=[a, b, c, alpha, beta, gamma])
+#hfo2.write('hfo2.xyz')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bondatoms = []
+symbols = hfo2.get_chemical_symbols()
+for i in range(len(hfo2)):
+    for j in range(i):
+        if (symbols[i] == 'Hf' and symbols[j] == 'O' and hfo2.get_distance(i, j) < 2.6):
+            bondatoms.append((i, j))
+        elif (symbols[i] == 'O' and symbols[j] == 'Hf' and hfo2.get_distance(i, j) < 2.6):
+            bondatoms.append((i, j))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+nglview.show_ase(hfo2)
+```
+
+%%%% Output: display_data
+
+
+%% Cell type:markdown id: tags:
+
+Questions:
+
+    -how many atoms are contained in the unit cell?
+    -compute the volume of the unit cell 
+
+%% Cell type:markdown id: tags:
+
+## 
+
+%% Cell type:code id: tags:
+
+``` python
+def ase_to_spgcell(ase_atoms=None, cell=None, inverse=False):
+    if not inverse:
+        assert ase_atoms is not None
+        return (ase_atoms.get_cell(),
+                ase_atoms.get_scaled_positions(),
+                ase_atoms.get_atomic_numbers())
+    else:
+        assert cell is not None
+        return Atoms(cell=cell[0],
+                     scaled_positions=cell[1],
+                     numbers=cell[2])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def a_equiv_b(a,b):
+    """Function that identifies whether two crystals are equivalent"""
+
+    # getting symmetry datasets for both crystals
+    cryst_a = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=a), symprec=1e-5, angle_tolerance=-1.0, hall_number=0)
+    cryst_b = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=b), symprec=1e-5, angle_tolerance=-1.0, hall_number=0)
+
+    samecell = np.allclose(cryst_a['std_lattice'], cryst_b['std_lattice'], atol=1e-5)
+    samenatoms = len(cryst_a['std_positions']) == len(cryst_b['std_positions'])
+    samespg = cryst_a['number'] == cryst_b['number']
+    
+    def test_rotations_translations(cryst_a, cryst_b, repeat):
+        cell = cryst_a['std_lattice']
+        pristine = crystal('Mg', [(0, 0., 0.)], 
+                           spacegroup=int(cryst_a['number']),
+                           cellpar=[cell[0]/repeat[0], cell[1]/repeat[1], cell[2]/repeat[2]]).repeat(repeat)
+
+        sym_set_p = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=pristine), symprec=1e-5,
+                                               angle_tolerance=-1.0, hall_number=0)
+
+        for _,trans in enumerate(zip(sym_set_p['rotations'], sym_set_p['translations'])):
+            pnew=(np.matmul(trans[0],cryst_a['std_positions'].T).T + trans[1]) % 1.0
+            fulln = np.concatenate([cryst_a['std_types'][:, None], pnew], axis=1)
+            fullb = np.concatenate([cryst_b['std_types'][:, None], cryst_b['std_positions']], axis=1)
+            sorted_n = np.array(sorted([ list(row) for row in list(fulln) ]))
+            sorted_b = np.array(sorted([ list(row) for row in list(fullb) ]))
+            if np.allclose(sorted_n, sorted_b, atol=1e-5):
+                return True
+        return False
+
+    if samecell and samenatoms and samespg:
+        cell = cryst_a['std_lattice']
+        rng1 = range(1, int(norm(cell[0])/2.))
+        rng2 = range(1, int(norm(cell[1])/2.))
+        rng3 = range(1, int(norm(cell[2])/2.))
+
+        for repeat in itertools.product(rng1, rng2, rng3):
+            if test_rotations_translations(cryst_a, cryst_b, repeat):
+                return True
+
+    return False
+```
+
+%% Cell type:markdown id: tags:
+
+## Task 1
+
+%% Cell type:code id: tags:
+
+``` python
+atoms = read('./quartz_alpha.xyz')
+nglview.show_ase(atoms)
+```
+
+%%%% Output: display_data
+
+
+%% Cell type:code id: tags:
+
+``` python
+cell = ase_to_spgcell(ase_atoms=atoms)
+lattice, scaled_positions, numbers = spglib.find_primitive(cell, symprec=1e-5)
+reduced = ase_to_spgcell(cell=(lattice, scaled_positions, numbers),inverse=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+nglview.show_ase(reduced)
+```
+
+%%%% Output: display_data
+
+
+%% Cell type:markdown id: tags:
+
+## Task 2
+
+%% Cell type:code id: tags:
+
+``` python
+versions = [read(atoms) for atoms in glob('./quartz_alpha*')]
+
+for outer in range(len(versions)-1):
+    for inner in range(outer+1,len(versions)):
+        print('Comparing: ',outer,inner)
+        print('Equivalent: ',a_equiv_b(versions[outer],versions[inner]))
+```

Exercises/Exercise-10/.ipynb_checkpoints/Solution-10-checkpoint.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import itertools
+import numpy as np
+from numpy.linalg import norm
+import spglib
+from glob import glob
+
+from ase import Atoms
+from ase.io import read
+from ase.spacegroup import crystal
+from ase.spacegroup import Spacegroup
+from ase.data import atomic_numbers, atomic_names
+import nglview
+```
+
+%% Cell type:markdown id: tags:
+
+## Crystallographic point groups
+
+%% Cell type:code id: tags:
+
+``` python
+a=...
+b=...
+c=...
+alpha=...
+beta=...
+gamma=...
+thespacegroup=...
+print Spacegroup(thespacegroup)
+
+hfo2 = crystal(['...','...','...'],basis=[(...,...,...),(...,...,...),(...,...,...)], 
+               spacegroup=thespacegroup, cellpar=[a, b, c, alpha, beta, gamma])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bondatoms = []
+symbols = hfo2.get_chemical_symbols()
+for i in range(len(hfo2)):
+    for j in range(i):
+        if (symbols[i] == 'Hf' and symbols[j] == 'O' and hfo2.get_distance(i, j) < 2.6):
+            bondatoms.append((i, j))
+        elif (symbols[i] == 'O' and symbols[j] == 'Hf' and hfo2.get_distance(i, j) < 2.6):
+            bondatoms.append((i, j))
+```
+
+%% Cell type:markdown id: tags:
+
+Questions:
+
+    -how many atoms are contained in the unit cell?
+    -compute the volume of the unit cell 
+
+%% Cell type:markdown id: tags:
+
+## 
+
+%% Cell type:code id: tags:
+
+``` python
+def ase_to_spgcell(ase_atoms=None, cell=None, inverse=False):
+    if not inverse:
+        assert ase_atoms is not None
+        return (ase_atoms.get_cell(),
+                ase_atoms.get_scaled_positions(),
+                ase_atoms.get_atomic_numbers())
+    else:
+        assert cell is not None
+        return Atoms(cell=cell[0],
+                     scaled_positions=cell[1],
+                     numbers=cell[2])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def a_equiv_b(a,b):
+    """Function that identifies whether two crystals are equivalent"""
+
+    # getting symmetry datasets for both crystals
+    cryst_a = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=a), symprec=1e-5, angle_tolerance=-1.0, hall_number=0)
+    cryst_b = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=b), symprec=1e-5, angle_tolerance=-1.0, hall_number=0)
+
+    samecell = np.allclose(cryst_a['std_lattice'], cryst_b['std_lattice'], atol=1e-5)
+    samenatoms = len(cryst_a['std_positions']) == len(cryst_b['std_positions'])
+    samespg = cryst_a['number'] == cryst_b['number']
+    
+    def test_rotations_translations(cryst_a, cryst_b, repeat):
+        cell = cryst_a['std_lattice']
+        pristine = crystal('Mg', [(0, 0., 0.)], 
+                           spacegroup=int(cryst_a['number']),
+                           cellpar=[cell[0]/repeat[0], cell[1]/repeat[1], cell[2]/repeat[2]]).repeat(repeat)
+
+        sym_set_p = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=pristine), symprec=1e-5,
+                                               angle_tolerance=-1.0, hall_number=0)
+
+        for _,trans in enumerate(zip(sym_set_p['rotations'], sym_set_p['translations'])):
+            pnew=(np.matmul(trans[0],cryst_a['std_positions'].T).T + trans[1]) % 1.0
+            fulln = np.concatenate([cryst_a['std_types'][:, None], pnew], axis=1)
+            fullb = np.concatenate([cryst_b['std_types'][:, None], cryst_b['std_positions']], axis=1)
+            sorted_n = np.array(sorted([ list(row) for row in list(fulln) ]))
+            sorted_b = np.array(sorted([ list(row) for row in list(fullb) ]))
+            if np.allclose(sorted_n, sorted_b, atol=1e-5):
+                return True
+        return False
+
+    if samecell and samenatoms and samespg:
+        cell = cryst_a['std_lattice']
+        rng1 = range(1, int(norm(cell[0])/2.))
+        rng2 = range(1, int(norm(cell[1])/2.))
+        rng3 = range(1, int(norm(cell[2])/2.))
+
+        for repeat in itertools.product(rng1, rng2, rng3):
+            if test_rotations_translations(cryst_a, cryst_b, repeat):
+                return True
+
+    return False
+```
+
+%% Cell type:code id: tags:
+
+``` python
+atoms = read('./quartz_alpha.xyz')
+```
+
+%% Cell type:code id: tags:
+
+``` python
+view(atoms)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+cell = ase_to_spgcell(ase_atoms=atoms)
+lattice, scaled_positions, numbers = spglib.find_primitive(cell, symprec=1e-5)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+reduced = ase_to_spgcell(cell=(lattice, scaled_positions, numbers),inverse=True)
+```
+
+%% Cell type:code id: tags:
+
+``` python
+versions = [read(atoms) for atoms in glob('./quartz_alpha*')]
+```
+
+%% Cell type:code id: tags:
+
+``` python
+for outer in range(len(versions)-1):
+    for inner in range(outer+1,len(versions)):
+        print('Comparing: ',outer,inner)
+        print('Equivalent: ',a_equiv_b(versions[outer],versions[inner]))
+```

Exercises/Exercise-10/Exercise-10.ipynb

+%% Cell type:code id: tags:
+
+``` python
+import itertools
+import numpy as np
+from numpy.linalg import norm
+import spglib
+from glob import glob
+
+from ase import Atoms
+from ase.io import read
+from ase.spacegroup import crystal
+from ase.spacegroup import Spacegroup
+from ase.data import atomic_numbers, atomic_names
+import nglview
+```
+
+%% Cell type:markdown id: tags:
+
+## Crystallographic point groups
+
+%% Cell type:code id: tags:
+
+``` python
+a=...
+b=...
+c=...
+alpha=...
+beta=...
+gamma=...
+thespacegroup=...
+print(Spacegroup(thespacegroup))
+
+hfo2 = crystal(['...','...','...'],basis=[(...,...,...),(...,...,...),(...,...,...)], 
+               spacegroup=thespacegroup, cellpar=[a, b, c, alpha, beta, gamma])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+bondatoms = []
+symbols = hfo2.get_chemical_symbols()
+for i in range(len(hfo2)):
+    for j in range(i):
+        if (symbols[i] == 'Hf' and symbols[j] == 'O' and hfo2.get_distance(i, j) < 2.6):
+            bondatoms.append((i, j))
+        elif (symbols[i] == 'O' and symbols[j] == 'Hf' and hfo2.get_distance(i, j) < 2.6):
+            bondatoms.append((i, j))
+```
+
+%% Cell type:code id: tags:
+
+``` python
+nglview.show_ase(hfo2)
+```
+
+%%%% Output: display_data
+
+
+%% Cell type:markdown id: tags:
+
+Questions:
+
+    -how many atoms are contained in the unit cell?
+    -compute the volume of the unit cell 
+
+%% Cell type:markdown id: tags:
+
+## 
+
+%% Cell type:code id: tags:
+
+``` python
+def ase_to_spgcell(ase_atoms=None, cell=None, inverse=False):
+    if not inverse:
+        assert ase_atoms is not None
+        return (ase_atoms.get_cell(),
+                ase_atoms.get_scaled_positions(),
+                ase_atoms.get_atomic_numbers())
+    else:
+        assert cell is not None
+        return Atoms(cell=cell[0],
+                     scaled_positions=cell[1],
+                     numbers=cell[2])
+```
+
+%% Cell type:code id: tags:
+
+``` python
+def a_equiv_b(a,b):
+    """Function that identifies whether two crystals are equivalent"""
+
+    # getting symmetry datasets for both crystals
+    cryst_a = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=a), symprec=1e-5, angle_tolerance=-1.0, hall_number=0)
+    cryst_b = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=b), symprec=1e-5, angle_tolerance=-1.0, hall_number=0)
+
+    samecell = np.allclose(cryst_a['std_lattice'], cryst_b['std_lattice'], atol=1e-5)
+    samenatoms = len(cryst_a['std_positions']) == len(cryst_b['std_positions'])
+    samespg = cryst_a['number'] == cryst_b['number']
+    
+    def test_rotations_translations(cryst_a, cryst_b, repeat):
+        cell = cryst_a['std_lattice']
+        pristine = crystal('Mg', [(0, 0., 0.)], 
+                           spacegroup=int(cryst_a['number']),
+                           cellpar=[cell[0]/repeat[0], cell[1]/repeat[1], cell[2]/repeat[2]]).repeat(repeat)
+
+        sym_set_p = spglib.get_symmetry_dataset(ase_to_spgcell(ase_atoms=pristine), symprec=1e-5,
+                                               angle_tolerance=-1.0, hall_number=0)
+
+        for _,trans in enumerate(zip(sym_set_p['rotations'], sym_set_p['translations'])):
+            pnew=(np.matmul(trans[0],cryst_a['std_positions'].T).T + trans[1]) % 1.0
+            fulln = np.concatenate([cryst_a['std_types'][:, None], pnew], axis=1)
+            fullb = np.concatenate([cryst_b['std_types'][:, None], cryst_b['std_positions']], axis=1)
+            sorted_n = np.array(sorted([ list(row) for row in list(fulln) ]))
+            sorted_b = np.array(sorted([ list(row) for row in list(fullb) ]))
+            if np.allclose(sorted_n, sorted_b, atol=1e-5):
+                return True
+        return False
+
+    if samecell and samenatoms and samespg:
+        cell = cryst_a['std_lattice']
+        rng1 = range(1, int(norm(cell[0])/2.))
+        rng2 = range(1, int(norm(cell[1])/2.))
+        rng3 = range(1, int(norm(cell[2])/2.))
+
+        for repeat in itertools.product(rng1, rng2, rng3):
+            if test_rotations_translations(cryst_a, cryst_b, repeat):
+                return True
+
+    return False
+```
+
+%% Cell type:markdown id: tags:
+
+## Task 1
+
+%% Cell type:code id: tags:
+
+``` python
+atoms = read('./quartz_alpha.xyz')
+nglview.show_ase(atoms)
+```
+
+%%%% Output: display_data
+
+
+%% Cell type:code id: tags:
+
+``` python
+''' 
+
+    ..TO DO..
+
+'''
+```
+
+%% Cell type:code id: tags:
+
+``` python
+nglview.show_ase(reduced)
+```
+
+%%%% Output: display_data
+
+
+%% Cell type:markdown id: tags:
+
+## Task 2
+
+%% Cell type:code id: tags:
+
+``` python
+''' 
+
+    ..TO DO..
+
+'''
+```

Exercises/Exercise-10/hfo2.xyz

+12
+Lattice="5.117 0.0 0.0 0.0 5.175 0.0 -0.8477544617904331 0.0 5.222642374556621" Properties=species:S:1:pos:R:3:tags:I:1 spacegroup="P 21/c" unit_cell=conventional pbc="T T T"
+Hf       1.23595907       0.20700000       1.08630961        0 
+Hf       3.45716370       2.79450000       1.52501157        0 
+Hf       3.03328647       4.96800000       4.13633276        0 
+Hf       0.81208184       2.38050000       3.69763080        0 
+O       0.08448720       1.71810000       1.81225690        1 
+O       4.60863557       4.30560000       0.79906428        1 
+O       4.18475834       3.45690000       3.41038547        1 
+O      -0.33939003       0.86940000       4.42357809        1 
+O       1.89061086       3.92265000       2.50686834        2 
+O       2.80251191       1.33515000       0.10445285        2 
+O       2.37863468       1.25235000       2.71577403        2 
+O       1.46673363       3.83985000       5.11818953        2 

Exercises/Exercise-10/quartz_alpha.xyz

+144
+Lattice="9.81999180292937 -1.57917512379296e-15 0.0 0.0 17.0087247325836 0.0 0.0 0.0 10.804" Properties=species:S:1:pos:R:3 pbc="T T T"
+O       1.38650516       1.12431739       0.62933000 
+O       3.84150311       5.37649857       0.62933000 
+O       0.78805795       4.89077118       4.23066333 
+O       3.24305590       0.63859000       4.23066333 
+O       0.28043484       6.74145498       2.42999667 
+O       2.73543279       2.48927380       2.42999667 
+O       1.38650516       7.38004498       4.77267000 
+O       3.84150311       3.12786380       4.77267000 
+O       0.78805795       3.61359118       1.17133667 
+O       3.24305590       7.86577237       1.17133667 
+O       0.28043484       1.76290739       2.97200333 
+O       2.73543279       6.01508857       2.97200333 
+Si       2.36063135      -0.00000000       0.00000000 
+Si       4.81562931       4.25218118       0.00000000 
+Si       1.27468227       6.29654791       3.60133333 
+Si       3.72968022       2.04436672       3.60133333 
+Si       1.27468227       2.20781446       1.80066667 
+Si       3.72968022       6.45999564       1.80066667 
+O       6.29650106       1.12431739       0.62933000 
+O       8.75149901       5.37649857       0.62933000 
+O       5.69805385       4.89077118       4.23066333 
+O       8.15305180       0.63859000       4.23066333 
+O       5.19043074       6.74145498       2.42999667 
+O       7.64542869       2.48927380       2.42999667 
+O       6.29650106       7.38004498       4.77267000 
+O       8.75149901       3.12786380       4.77267000 
+O       5.69805385       3.61359118       1.17133667 
+O       8.15305180       7.86577237       1.17133667 
+O       5.19043074       1.76290739       2.97200333 
+O       7.64542869       6.01508857       2.97200333 
+Si       7.27062726      -0.00000000       0.00000000 
+Si       9.72562521       4.25218118       0.00000000 
+Si       6.18467817       6.29654791       3.60133333 
+Si       8.63967613       2.04436672       3.60133333