1# Copyright (C) 2011 Atsushi Togo 2# All rights reserved. 3# 4# This file is part of phonopy. 5# 6# Redistribution and use in source and binary forms, with or without 7# modification, are permitted provided that the following conditions 8# are met: 9# 10# * Redistributions of source code must retain the above copyright 11# notice, this list of conditions and the following disclaimer. 12# 13# * Redistributions in binary form must reproduce the above copyright 14# notice, this list of conditions and the following disclaimer in 15# the documentation and/or other materials provided with the 16# distribution. 17# 18# * Neither the name of the phonopy project nor the names of its 19# contributors may be used to endorse or promote products derived 20# from this software without specific prior written permission. 21# 22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 26# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 29# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 30# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 31# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 32# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 33# POSSIBILITY OF SUCH DAMAGE. 34 35import numpy as np 36 37directions_axis = np.array([[1, 0, 0], 38 [0, 1, 0], 39 [0, 0, 1]]) 40 41directions_diag = np.array([[1, 0, 0], 42 [0, 1, 0], 43 [0, 0, 1], 44 [1, 1, 0], 45 [1, 0, 1], 46 [0, 1, 1], 47 [1, -1, 0], 48 [1, 0, -1], 49 [0, 1, -1], 50 [1, 1, 1], 51 [1, 1, -1], 52 [1, -1, 1], 53 [-1, 1, 1]]) 54 55 56def directions_to_displacement_dataset(displacement_directions, 57 distance, 58 supercell): 59 lattice = supercell.get_cell() 60 first_atoms = [] 61 for disp in displacement_directions: 62 direction = disp[1:] 63 disp_cartesian = np.dot(direction, lattice) 64 disp_cartesian *= distance / np.linalg.norm(disp_cartesian) 65 first_atoms.append({'number': int(disp[0]), 66 'displacement': disp_cartesian.tolist()}) 67 displacement_dataset = { 68 'natom': supercell.get_number_of_atoms(), 69 'first_atoms': first_atoms} 70 71 return displacement_dataset 72 73 74def get_least_displacements(symmetry, 75 is_plusminus='auto', 76 is_diagonal=True, 77 is_trigonal=False, 78 log_level=0): 79 """Return a set of displacements 80 81 Returns 82 ------- 83 array_like 84 List of directions with respect to axes. This gives only the 85 symmetrically non equivalent directions. The format is like: 86 [[0, 1, 0, 0], 87 [7, 1, 0, 1], ...] 88 where each list is defined by: 89 First value: Atom index in supercell starting with 0 90 Second to fourth: If the direction is displaced or not (1, 0, or -1) 91 with respect to the axes. 92 93 """ 94 displacements = [] 95 if is_diagonal: 96 directions = directions_diag 97 else: 98 directions = directions_axis 99 100 if log_level > 2: 101 print("Site point symmetry:") 102 103 for atom_num in symmetry.get_independent_atoms(): 104 site_symmetry = symmetry.get_site_symmetry(atom_num) 105 106 if log_level > 2: 107 print("Atom %d" % (atom_num + 1)) 108 for i, rot in enumerate(site_symmetry): 109 print("----%d----" % (i + 1)) 110 for v in rot: 111 print("%2d %2d %2d" % tuple(v)) 112 113 for disp in get_displacement(site_symmetry, 114 directions, 115 is_trigonal, 116 log_level): 117 displacements.append([atom_num, 118 disp[0], disp[1], disp[2]]) 119 if is_plusminus == 'auto': 120 if is_minus_displacement(disp, site_symmetry): 121 displacements.append([atom_num, 122 -disp[0], -disp[1], -disp[2]]) 123 elif is_plusminus is True: 124 displacements.append([atom_num, 125 -disp[0], -disp[1], -disp[2]]) 126 127 return displacements 128 129 130def get_displacement(site_symmetry, 131 directions=directions_diag, 132 is_trigonal=False, 133 log_level=0): 134 # One 135 sitesym_num, disp = get_displacement_one(site_symmetry, 136 directions) 137 if disp is not None: 138 if log_level > 2: 139 print("Site symmetry used to expand a direction %s" % disp[0]) 140 print(site_symmetry[sitesym_num]) 141 return disp 142 # Two 143 sitesym_num, disps = get_displacement_two(site_symmetry, 144 directions) 145 if disps is not None: 146 if log_level > 2: 147 print("Site symmetry used to expand directions %s %s" % 148 (disps[0], disps[1])) 149 print(site_symmetry[sitesym_num]) 150 151 if is_trigonal: 152 disps_new = [disps[0]] 153 if is_trigonal_axis(site_symmetry[sitesym_num]): 154 if log_level > 2: 155 print("Trigonal axis is found.") 156 disps_new.append(np.dot(disps[0], 157 site_symmetry[sitesym_num].T)) 158 disps_new.append(np.dot(disps_new[1], 159 site_symmetry[sitesym_num].T)) 160 disps_new.append(disps[1]) 161 return disps_new 162 else: 163 return disps 164 # Three 165 return [directions[0], directions[1], directions[2]] 166 167 168def get_displacement_one(site_symmetry, 169 directions=directions_diag): 170 for direction in directions: 171 rot_directions = [] 172 for r in site_symmetry: 173 rot_directions.append(np.dot(direction, r.T)) 174 num_sitesym = len(site_symmetry) 175 for i in range(num_sitesym): 176 for j in range(i+1, num_sitesym): 177 det = determinant(direction, 178 rot_directions[i], 179 rot_directions[j]) 180 if det != 0: 181 return i, [direction] 182 return None, None 183 184 185def get_displacement_two(site_symmetry, 186 directions=directions_diag): 187 for direction in directions: 188 rot_directions = [] 189 for r in site_symmetry: 190 rot_directions.append(np.dot(direction, r.T)) 191 num_sitesym = len(site_symmetry) 192 for i in range(num_sitesym): 193 for second_direction in directions: 194 det = determinant(direction, 195 rot_directions[i], 196 second_direction) 197 if det != 0: 198 return i, [direction, second_direction] 199 return None, None 200 201 202def is_minus_displacement(direction, site_symmetry): 203 is_minus = True 204 for r in site_symmetry: 205 rot_direction = np.dot(direction, r.T) 206 if (rot_direction + direction).any(): 207 continue 208 else: 209 is_minus = False 210 break 211 return is_minus 212 213 214def is_trigonal_axis(r): 215 r3 = np.dot(np.dot(r, r), r) 216 if (r3 == np.eye(3, dtype=int)).all(): 217 return True 218 else: 219 return False 220 221 222def determinant(a, b, c): 223 det = (a[0] * b[1] * c[2] - a[0] * b[2] * c[1] 224 + a[1] * b[2] * c[0] - a[1] * b[0] * c[2] 225 + a[2] * b[0] * c[1] - a[2] * b[1] * c[0]) 226 return det 227 228 229def get_random_displacements_dataset(num_supercells, 230 distance, 231 num_atoms, 232 random_seed=None): 233 if (np.issubdtype(type(random_seed), np.integer) and 234 random_seed >= 0 and random_seed < 2 ** 32): 235 seed = random_seed 236 else: 237 seed = None 238 239 disps = get_random_directions(num_atoms * num_supercells, 240 random_seed=random_seed) * distance 241 supercell_disps = np.array(disps.reshape(num_supercells, num_atoms, 3), 242 dtype='double', order='C') 243 dataset = {'displacements': supercell_disps} 244 245 if seed is not None: 246 dataset['random_seed'] = seed 247 return dataset 248 249 250def get_random_directions(num_atoms, random_seed=None): 251 """Returns random directions in sphere with radius 1""" 252 253 if (np.issubdtype(type(random_seed), np.integer) and 254 random_seed >= 0 and random_seed < 2 ** 32): 255 np.random.seed(random_seed) 256 257 xy = np.random.randn(3, num_atoms) 258 r = np.sqrt((xy ** 2).sum(axis=0)) 259 return (xy / r).T 260 261 262def print_displacements(symmetry, 263 directions=directions_diag): 264 displacements = get_least_displacements(symmetry, directions) 265 print("Least displacements:") 266 print("Atom Directions") 267 print("----------------------------") 268 for key in displacements: 269 print("%4d %s" % (key + 1, displacements[key])) 270