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