1#!/usr/bin/evn python
2
3import sys
4import spglib
5import numpy as np
6
7#######################################
8# Uncomment to use Atoms class in ASE #
9#######################################
10# from ase import Atoms
11
12#######################################################################
13# Using the local Atoms-like class (BSD license) where a small set of #
14# ASE Atoms features is comaptible but enough for this example.       #
15#######################################################################
16from atoms import Atoms
17
18def show_symmetry(symmetry):
19    for i in range(symmetry['rotations'].shape[0]):
20        print("  --------------- %4d ---------------" % (i + 1))
21        rot = symmetry['rotations'][i]
22        trans = symmetry['translations'][i]
23        print("  rotation:")
24        for x in rot:
25            print("     [%2d %2d %2d]" % (x[0], x[1], x[2]))
26        print("  translation:")
27        print("     (%8.5f %8.5f %8.5f)" % (trans[0], trans[1], trans[2]))
28
29def show_lattice(lattice):
30    print("Basis vectors:")
31    for vec, axis in zip(lattice, ("a", "b", "c")):
32        print("%s %10.5f %10.5f %10.5f" % (tuple(axis,) + tuple(vec)))
33
34def show_cell(lattice, positions, numbers):
35    show_lattice(lattice)
36    print("Atomic points:")
37    for p, s in zip(positions, numbers):
38        print("%2d %10.5f %10.5f %10.5f" % ((s,) + tuple(p)))
39
40silicon_ase = Atoms(symbols=['Si'] * 8,
41                    cell=[(4, 0, 0),
42                          (0, 4, 0),
43                          (0, 0, 4)],
44                    scaled_positions=[(0, 0, 0),
45                                      (0, 0.5, 0.5),
46                                      (0.5, 0, 0.5),
47                                      (0.5, 0.5, 0),
48                                      (0.25, 0.25, 0.25),
49                                      (0.25, 0.75, 0.75),
50                                      (0.75, 0.25, 0.75),
51                                      (0.75, 0.75, 0.25)],
52                    pbc=True)
53
54silicon = ([(4, 0, 0),
55            (0, 4, 0),
56            (0, 0, 4)],
57           [(0, 0, 0),
58            (0, 0.5, 0.5),
59            (0.5, 0, 0.5),
60            (0.5, 0.5, 0),
61            (0.25, 0.25, 0.25),
62            (0.25, 0.75, 0.75),
63            (0.75, 0.25, 0.75),
64            (0.75, 0.75, 0.25)],
65           [14,] * 8)
66
67silicon_dist = ([(4.01, 0, 0),
68                 (0, 4, 0),
69                 (0, 0, 3.99)],
70                [(0.001, 0, 0),
71                 (0, 0.5, 0.5),
72                 (0.5, 0, 0.5),
73                 (0.5, 0.5, 0),
74                 (0.25, 0.25, 0.251),
75                 (0.25, 0.75, 0.75),
76                 (0.75, 0.25, 0.75),
77                 (0.75, 0.75, 0.25)],
78                [14,] * 8)
79
80silicon_prim = ([(0, 2, 2),
81                 (2, 0, 2),
82                 (2, 2, 0)],
83                [(0, 0, 0),
84                 (0.25, 0.25, 0.25)],
85                [14, 14])
86
87rutile = ([(4, 0, 0),
88           (0, 4, 0),
89           (0, 0, 3)],
90          [(0, 0, 0),
91           (0.5, 0.5, 0.5),
92           (0.3, 0.3, 0.0),
93           (0.7, 0.7, 0.0),
94           (0.2, 0.8, 0.5),
95           (0.8, 0.2, 0.5)],
96          [14, 14, 8, 8, 8, 8])
97
98rutile_dist = ([(3.97, 0, 0),
99                (0, 4.03, 0),
100                (0, 0, 3.0)],
101               [(0, 0, 0),
102                (0.5001, 0.5, 0.5),
103                (0.3, 0.3, 0.0),
104                (0.7, 0.7, 0.002),
105                (0.2, 0.8, 0.5),
106                (0.8, 0.2, 0.5)],
107               [14, 14, 8, 8, 8, 8])
108
109a = 3.07
110c = 3.52
111MgB2 = ([(a, 0, 0),
112         (-a/2, a/2*np.sqrt(3), 0),
113         (0, 0, c)],
114        [(0, 0, 0),
115         (1.0/3, 2.0/3, 0.5),
116         (2.0/3, 1.0/3, 0.5)],
117        [12, 5, 5])
118
119a = [3., 0., 0.]
120b = [-3.66666667, 3.68178701, 0.]
121c = [-0.66666667, -1.3429469, 1.32364995]
122niggli_lattice = np.array([a, b, c])
123
124# For VASP case
125# import vasp
126# bulk = vasp.read_vasp(sys.argv[1])
127
128print("[get_spacegroup]")
129print("  Spacegroup of Silicon is %s." % spglib.get_spacegroup(silicon))
130print('')
131print("[get_spacegroup]")
132print("  Spacegroup of Silicon (ASE Atoms-like format) is %s." %
133      spglib.get_spacegroup(silicon_ase))
134print('')
135print("[get_spacegroup]")
136print("  Spacegroup of Rutile is %s." % spglib.get_spacegroup(rutile))
137print('')
138print("[get_spacegroup]")
139print("  Spacegroup of MgB2 is %s." % spglib.get_spacegroup(MgB2))
140print('')
141print("[get_symmetry]")
142print("  Symmetry operations of Rutile unitcell are:")
143print('')
144symmetry = spglib.get_symmetry(rutile)
145show_symmetry(symmetry)
146print('')
147print("[get_symmetry]")
148print("  Symmetry operations of MgB2 are:")
149print('')
150symmetry = spglib.get_symmetry(MgB2)
151show_symmetry(symmetry)
152print('')
153print("[get_pointgroup]")
154print("  Pointgroup of Rutile is %s." %
155      spglib.get_pointgroup(symmetry['rotations'])[0])
156print('')
157
158dataset = spglib.get_symmetry_dataset( rutile )
159print("[get_symmetry_dataset] ['international']")
160print("  Spacegroup of Rutile is %s (%d)." % (dataset['international'],
161                                              dataset['number']))
162print('')
163print("[get_symmetry_dataset] ['pointgroup']")
164print("  Pointgroup of Rutile is %s." % (dataset['pointgroup']))
165print('')
166print("[get_symmetry_dataset] ['hall']")
167print("  Hall symbol of Rutile is %s (%d)." % (dataset['hall'],
168                                               dataset['hall_number']))
169print('')
170print("[get_symmetry_dataset] ['wyckoffs']")
171alphabet = "abcdefghijklmnopqrstuvwxyz"
172print("  Wyckoff letters of Rutile are: ", dataset['wyckoffs'])
173print('')
174print("[get_symmetry_dataset] ['equivalent_atoms']")
175print("  Mapping to equivalent atoms of Rutile are: ")
176for i, x in enumerate(dataset['equivalent_atoms']):
177  print("  %d -> %d" % (i + 1, x + 1))
178print('')
179print("[get_symmetry_dataset] ['rotations'], ['translations']")
180print("  Symmetry operations of Rutile unitcell are:")
181for i, (rot,trans) in enumerate(zip(dataset['rotations'],
182                                    dataset['translations'])):
183    print("  --------------- %4d ---------------" % (i + 1))
184    print("  rotation:")
185    for x in rot:
186        print("     [%2d %2d %2d]" % (x[0], x[1], x[2]))
187    print("  translation:")
188    print("     (%8.5f %8.5f %8.5f)" % (trans[0], trans[1], trans[2]))
189print('')
190
191print("[refine_cell]")
192print(" Refine distorted rutile structure")
193lattice, positions, numbers = spglib.refine_cell(rutile_dist, symprec=1e-1)
194show_cell(lattice, positions, numbers)
195print('')
196
197print("[find_primitive]")
198print(" Fine primitive distorted silicon structure")
199lattice, positions, numbers = spglib.find_primitive(silicon_dist, symprec=1e-1)
200show_cell(lattice, positions, numbers)
201print('')
202
203print("[standardize_cell]")
204print(" Standardize distorted rutile structure:")
205print(" (to_primitive=0 and no_idealize=0)")
206lattice, positions, numbers = spglib.standardize_cell(rutile_dist,
207                                                      to_primitive=0,
208                                                      no_idealize=0,
209                                                      symprec=1e-1)
210show_cell(lattice, positions, numbers)
211print('')
212
213print("[standardize_cell]")
214print(" Standardize distorted rutile structure:")
215print(" (to_primitive=0 and no_idealize=1)")
216lattice, positions, numbers = spglib.standardize_cell(rutile_dist,
217                                                      to_primitive=0,
218                                                      no_idealize=1,
219                                                      symprec=1e-1)
220show_cell(lattice, positions, numbers)
221print('')
222
223print("[standardize_cell]")
224print(" Standardize distorted silicon structure:")
225print(" (to_primitive=1 and no_idealize=0)")
226lattice, positions, numbers = spglib.standardize_cell(silicon_dist,
227                                                      to_primitive=1,
228                                                      no_idealize=0,
229                                                      symprec=1e-1)
230show_cell(lattice, positions, numbers)
231print('')
232
233print("[standardize_cell]")
234print(" Standardize distorted silicon structure:")
235print(" (to_primitive=1 and no_idealize=1)")
236lattice, positions, numbers = spglib.standardize_cell(silicon_dist,
237                                                      to_primitive=1,
238                                                      no_idealize=1,
239                                                      symprec=1e-1)
240show_cell(lattice, positions, numbers)
241print('')
242
243symmetry = spglib.get_symmetry(silicon)
244print("[get_symmetry]")
245print("  Number of symmetry operations of silicon conventional")
246print("  unit cell is %d (192)." % len(symmetry['rotations']))
247show_symmetry(symmetry)
248print('')
249
250symmetry = spglib.get_symmetry_from_database(525)
251print("[get_symmetry_from_database]")
252print("  Number of symmetry operations of silicon conventional")
253print("  unit cell is %d (192)." % len(symmetry['rotations']))
254show_symmetry(symmetry)
255print('')
256
257reduced_lattice = spglib.niggli_reduce(niggli_lattice)
258print("[niggli_reduce]")
259print("  Original lattice")
260show_lattice(niggli_lattice)
261print("  Reduced lattice")
262show_lattice(reduced_lattice)
263print('')
264
265
266mapping, grid = spglib.get_ir_reciprocal_mesh([11, 11, 11],
267                                              silicon_prim,
268                                              is_shift=[0, 0, 0])
269num_ir_kpt = len(np.unique(mapping))
270print("[get_ir_reciprocal_mesh]")
271print("  Number of irreducible k-points of primitive silicon with")
272print("  11x11x11 Monkhorst-Pack mesh is %d (56)." % num_ir_kpt)
273print('')
274
275mapping, grid = spglib.get_ir_reciprocal_mesh([8, 8, 8],
276                                              rutile,
277                                              is_shift=[1, 1, 1])
278num_ir_kpt = len(np.unique(mapping))
279print("[get_ir_reciprocal_mesh]")
280print("  Number of irreducible k-points of Rutile with")
281print("  8x8x8 Monkhorst-Pack mesh is %d (40)." % num_ir_kpt)
282print('')
283
284mapping, grid = spglib.get_ir_reciprocal_mesh([9, 9, 8],
285                                              MgB2,
286                                              is_shift=[0, 0, 1])
287num_ir_kpt = len(np.unique(mapping))
288print("[get_ir_reciprocal_mesh]")
289print("  Number of irreducible k-points of MgB2 with")
290print("  9x9x8 Monkhorst-Pack mesh is %s (48)." % num_ir_kpt)
291print('')
292