1"""Force constants calculation.""" 2# Copyright (C) 2011 Atsushi Togo 3# All rights reserved. 4# 5# This file is part of phonopy. 6# 7# Redistribution and use in source and binary forms, with or without 8# modification, are permitted provided that the following conditions 9# are met: 10# 11# * Redistributions of source code must retain the above copyright 12# notice, this list of conditions and the following disclaimer. 13# 14# * Redistributions in binary form must reproduce the above copyright 15# notice, this list of conditions and the following disclaimer in 16# the documentation and/or other materials provided with the 17# distribution. 18# 19# * Neither the name of the phonopy project nor the names of its 20# contributors may be used to endorse or promote products derived 21# from this software without specific prior written permission. 22# 23# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 24# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 25# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS 26# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE 27# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 28# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 29# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; 30# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER 31# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 32# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 33# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 34# POSSIBILITY OF SUCH DAMAGE. 35 36import textwrap 37import numpy as np 38from phonopy.structure.cells import (get_smallest_vectors, 39 compute_permutation_for_rotation) 40 41 42def get_force_constants(set_of_forces, 43 symmetry, 44 supercell, 45 atom_list=None, 46 decimals=None): 47 """Calculate force constants from disp-force dataset.""" 48 first_atoms = [{'number': x.get_atom_number(), 49 'displacement': x.get_displacement(), 50 'forces': x.get_forces()} for x in set_of_forces] 51 dataset = {'natom': supercell.get_number_of_atoms(), 52 'first_atoms': first_atoms} 53 force_constants = get_fc2(supercell, 54 symmetry, 55 dataset, 56 atom_list=atom_list) 57 return force_constants 58 59 60def get_fc2(supercell, 61 symmetry, 62 dataset, 63 atom_list=None, 64 decimals=None): 65 """Force constants are computed. 66 67 Force constants, Phi, are calculated from sets for forces, F, and 68 atomic displacement, d: 69 Phi = -F / d 70 This is solved by matrix pseudo-inversion. 71 Crystal symmetry is included when creating F and d matrices. 72 73 Returns 74 ------- 75 ndarray 76 Force constants[ i, j, a, b ] 77 i: Atom index of finitely displaced atom. 78 j: Atom index at which force on the atom is measured. 79 a, b: Cartesian direction indices = (0, 1, 2) for i and j, respectively 80 dtype=double 81 shape=(len(atom_list),n_satom,3,3), 82 83 """ 84 if atom_list is None: 85 fc_dim0 = len(supercell) 86 else: 87 fc_dim0 = len(atom_list) 88 89 force_constants = np.zeros((fc_dim0, len(supercell), 3, 3), 90 dtype='double', order='C') 91 92 # Fill force_constants[ displaced_atoms, all_atoms_in_supercell ] 93 atom_list_done = _get_force_constants_disps( 94 force_constants, 95 supercell, 96 dataset, 97 symmetry, 98 atom_list=atom_list) 99 100 rotations = symmetry.get_symmetry_operations()['rotations'] 101 lattice = np.array(supercell.cell.T, dtype='double', order='C') 102 permutations = symmetry.atomic_permutations 103 distribute_force_constants(force_constants, 104 atom_list_done, 105 lattice, 106 rotations, 107 permutations, 108 atom_list=atom_list) 109 110 if decimals: 111 force_constants = force_constants.round(decimals=decimals) 112 113 return force_constants 114 115 116def compact_fc_to_full_fc(phonon, compact_fc, log_level=0): 117 """Transform compact fc to full fc.""" 118 fc = np.zeros((compact_fc.shape[1], compact_fc.shape[1], 3, 3), 119 dtype='double', order='C') 120 fc[phonon.primitive.p2s_map] = compact_fc 121 distribute_force_constants_by_translations( 122 fc, phonon.primitive, phonon.supercell) 123 if log_level: 124 print("Force constants were expanded to full format.") 125 126 return fc 127 128 129def cutoff_force_constants(force_constants, 130 supercell, 131 primitive, 132 cutoff_radius, 133 symprec=1e-5): 134 """Set zero to force constants outside of cutoff distance. 135 136 Note 137 ---- 138 `force_constants` is overwritten. 139 140 """ 141 fc_shape = force_constants.shape 142 143 if fc_shape[0] == fc_shape[1]: 144 svecs, multi = get_smallest_vectors( 145 supercell.cell, 146 supercell.scaled_positions, 147 supercell.scaled_positions, 148 symprec=symprec, 149 store_dense_svecs=primitive.store_dense_svecs) 150 lattice = supercell.cell 151 else: 152 svecs, multi = primitive.get_smallest_vectors() 153 lattice = primitive.cell 154 155 if primitive.store_dense_svecs: 156 _svecs = svecs[multi[:, :, 1]] 157 else: 158 _svecs = svecs[:, :, 0, :] 159 160 min_distances = np.sqrt(np.sum(np.dot(_svecs, lattice) ** 2, axis=-1)) 161 162 for i in range(fc_shape[0]): 163 for j in range(fc_shape[1]): 164 if min_distances[j, i] > cutoff_radius: 165 force_constants[i, j] = 0.0 166 167 168def symmetrize_force_constants(force_constants, level=1): 169 """Symmetry force constants by translational and permutation symmetries. 170 171 Note 172 ---- 173 Schemes of symmetrization are slightly different between C and 174 python implementations. If these give very different results, the 175 original force constants are not reliable anyway. 176 177 Parameters 178 ---------- 179 force_constants: ndarray 180 Force constants. Symmetrized force constants are overwritten. 181 dtype=double 182 shape=(n_satom,n_satom,3,3) 183 primitive: Primitive 184 Primitive cell 185 level: int 186 Controls the number of times the following steps repeated: 187 1) Subtract drift force constants along row and column 188 2) Average fc and fc.T 189 190 """ 191 try: 192 import phonopy._phonopy as phonoc 193 phonoc.perm_trans_symmetrize_fc(force_constants, level) 194 except ImportError: 195 for i in range(level): 196 set_translational_invariance(force_constants) 197 set_permutation_symmetry(force_constants) 198 set_translational_invariance(force_constants) 199 200 201def symmetrize_compact_force_constants(force_constants, 202 primitive, 203 level=1): 204 """Symmetry force constants by translational and permutation symmetries. 205 206 Parameters 207 ---------- 208 force_constants: ndarray 209 Compact force constants. Symmetrized force constants are overwritten. 210 dtype=double 211 shape=(n_patom,n_satom,3,3) 212 primitive: Primitive 213 Primitive cell 214 level: int 215 Controls the number of times the following steps repeated: 216 1) Subtract drift force constants along row and column 217 2) Average fc and fc.T 218 219 """ 220 s2p_map = primitive.s2p_map 221 p2s_map = primitive.p2s_map 222 p2p_map = primitive.p2p_map 223 permutations = primitive.atomic_permutations 224 s2pp_map, nsym_list = get_nsym_list_and_s2pp(s2p_map, 225 p2p_map, 226 permutations) 227 try: 228 import phonopy._phonopy as phonoc 229 phonoc.perm_trans_symmetrize_compact_fc(force_constants, 230 permutations, 231 s2pp_map, 232 p2s_map, 233 nsym_list, 234 level) 235 except ImportError: 236 text = ("Import error at phonoc.perm_trans_symmetrize_compact_fc. " 237 "Corresponding pytono code is not implemented.") 238 raise RuntimeError(text) 239 240 241def distribute_force_constants(force_constants, 242 atom_list_done, 243 lattice, # column vectors 244 rotations, # scaled (fractional) 245 permutations, 246 atom_list=None): 247 """Fill force constants elements by symmetry.""" 248 map_atoms, map_syms = _get_sym_mappings_from_permutations( 249 permutations, atom_list_done) 250 rots_cartesian = np.array([similarity_transformation(lattice, r) 251 for r in rotations], 252 dtype='double', order='C') 253 if atom_list is None: 254 targets = np.arange(force_constants.shape[1], dtype='intc') 255 else: 256 targets = np.array(atom_list, dtype='intc') 257 import phonopy._phonopy as phonoc 258 259 phonoc.distribute_fc2(force_constants, 260 targets, 261 rots_cartesian, 262 permutations, 263 np.array(map_atoms, dtype='intc'), 264 np.array(map_syms, dtype='intc')) 265 266 267def distribute_force_constants_by_translations(fc, primitive, supercell): 268 """Distribute compact fc data to full fc by pure translations. 269 270 For example, the input fc has to be prepared in the following way 271 in advance: 272 273 fc = np.zeros((compact_fc.shape[1], compact_fc.shape[1], 3, 3), 274 dtype='double', order='C') 275 fc[primitive.p2s_map] = compact_fc 276 277 """ 278 s2p = primitive.s2p_map 279 p2s = primitive.p2s_map 280 positions = supercell.scaled_positions 281 lattice = supercell.cell.T 282 diff = positions - positions[p2s[0]] 283 trans = np.array(diff[np.where(s2p == p2s[0])[0]], 284 dtype='double', order='C') 285 rotations = np.array([np.eye(3, dtype='intc')] * len(trans), 286 dtype='intc', order='C') 287 permutations = primitive.atomic_permutations 288 distribute_force_constants(fc, p2s, lattice, rotations, permutations) 289 290 291def solve_force_constants(force_constants, 292 disp_atom_number, 293 displacements, 294 sets_of_forces, 295 supercell, 296 site_symmetry, 297 symprec, 298 atom_list=None): 299 """Calculate force constants elements of pairs from an atom.""" 300 if atom_list is None: 301 fc_index = disp_atom_number 302 else: 303 fc_index = np.where(disp_atom_number == atom_list)[0] 304 if len(fc_index) == 0: 305 raise RuntimeError 306 else: 307 fc_index = fc_index[0] 308 force_constants[fc_index] = _solve_force_constants_svd( 309 disp_atom_number, 310 displacements, 311 sets_of_forces, 312 supercell, 313 site_symmetry, 314 symprec) 315 return None 316 317 318def get_positions_sent_by_rot_inv(lattice, # column vectors 319 positions, 320 site_symmetry, 321 symprec): 322 """Return atom indices of positions sent by inverse site symmetries. 323 324 Rotated_positions[rot_map] == positions. 325 326 Note 327 ---- 328 This method is public because of being used by phono3py. 329 330 """ 331 rot_map_syms = [] 332 for sym in site_symmetry: 333 rot_map = compute_permutation_for_rotation(np.dot(positions, sym.T), 334 positions, 335 lattice, 336 symprec) 337 rot_map_syms.append(rot_map) 338 339 return np.array(rot_map_syms, dtype='intc', order='C') 340 341 342def get_rotated_displacement(displacements, site_sym_cart): 343 """Rotate displacements by site symmetry. 344 345 Note 346 ---- 347 This method is public because of being used by phono3py. 348 349 """ 350 rot_disps = [] 351 for u in displacements: 352 rot_disps.append([np.dot(sym, u) for sym in site_sym_cart]) 353 return np.array(np.reshape(rot_disps, (-1, 3)), dtype='double', order='C') 354 355 356def set_tensor_symmetry(force_constants, 357 lattice, # column vectors 358 positions, 359 symmetry): 360 """Full force constants are symmetrized using crystal symmetry. 361 362 This method extracts symmetrically equivalent sets of atomic pairs and 363 take sum of their force constants and average the sum. 364 365 """ 366 rotations = symmetry.get_symmetry_operations()['rotations'] 367 translations = symmetry.get_symmetry_operations()['translations'] 368 map_atoms = symmetry.get_map_atoms() 369 symprec = symmetry.tolerance 370 cart_rot = np.array([similarity_transformation(lattice, rot) 371 for rot in rotations]) 372 373 mapa = _get_atom_indices_by_symmetry(lattice, 374 positions, 375 rotations, 376 translations, 377 symprec) 378 fc_new = np.zeros_like(force_constants) 379 indep_atoms = symmetry.get_independent_atoms() 380 381 for i in indep_atoms: 382 fc_combined = np.zeros(force_constants.shape[1:], dtype='double') 383 num_equiv_atoms = _combine_force_constants_equivalent_atoms( 384 fc_combined, 385 force_constants, 386 i, 387 cart_rot, 388 map_atoms, 389 mapa) 390 num_sitesym = _average_force_constants_by_sitesym(fc_new, 391 fc_combined, 392 i, 393 cart_rot, 394 mapa) 395 396 assert num_equiv_atoms * num_sitesym == len(rotations) 397 398 permutations = symmetry.atomic_permutations 399 distribute_force_constants(fc_new, 400 indep_atoms, 401 lattice, 402 rotations, 403 permutations) 404 405 force_constants[:] = fc_new 406 407 408def set_tensor_symmetry_PJ(force_constants, 409 lattice, 410 positions, 411 symmetry): 412 """Full force constants are symmetrized using crystal symmetry. 413 414 This method extracts symmetrically equivalent sets of atomic pairs and 415 take sum of their force constants and average the sum. 416 417 """ 418 rotations = symmetry.get_symmetry_operations()['rotations'] 419 translations = symmetry.get_symmetry_operations()['translations'] 420 symprec = symmetry.tolerance 421 422 N = len(rotations) 423 424 mapa = _get_atom_indices_by_symmetry(lattice, 425 positions, 426 rotations, 427 translations, 428 symprec) 429 cart_rot = np.array([similarity_transformation(lattice, rot).T 430 for rot in rotations]) 431 cart_rot_inv = np.array([np.linalg.inv(rot) for rot in cart_rot]) 432 fcm = np.array([force_constants[mapa[n], :, :, :][:, mapa[n], :, :] 433 for n in range(N)]) 434 s = np.transpose(np.array([np.dot(cart_rot[n], 435 np.dot(fcm[n], cart_rot_inv[n])) 436 for n in range(N)]), (0, 2, 3, 1, 4)) 437 force_constants[:] = np.array(np.average(s, axis=0), 438 dtype='double', 439 order='C') 440 441 442def set_translational_invariance(force_constants): 443 """Impose translational invariance (Python version). 444 445 The type1 is quite simple implementation, which is just taking sum of the 446 force constants in an axis and an atom index. The sum has to be zero due to 447 the translational invariance. If the sum is not zero, this error is 448 uniformly subtracted from force constants. 449 450 """ 451 for i in range(2): 452 set_translational_invariance_per_index(force_constants, index=i) 453 454 455def set_translational_invariance_per_index(fc2, index=0): 456 """Impose translational invariance per index (Python version).""" 457 for i in range(fc2.shape[1 - index]): 458 for j, k in list(np.ndindex(3, 3)): 459 if index == 0: 460 fc2[:, i, j, k] -= np.sum( 461 fc2[:, i, j, k]) / fc2.shape[0] 462 else: 463 fc2[i, :, j, k] -= np.sum( 464 fc2[i, :, j, k]) / fc2.shape[1] 465 466 467def set_permutation_symmetry(force_constants): 468 """Enforce permutation symmetry to force constants (Python version). 469 470 This is done by 471 472 Phi_ij_ab = Phi_ji_ba 473 474 i, j: atom index 475 a, b: Cartesian axis index 476 477 This is not necessary for harmonic phonon calculation because this 478 condition is imposed when making dynamical matrix Hermite in 479 dynamical_matrix.py. 480 481 """ 482 fc_copy = force_constants.copy() 483 for i in range(force_constants.shape[0]): 484 for j in range(force_constants.shape[1]): 485 force_constants[i, j] = (force_constants[i, j] + 486 fc_copy[j, i].T) / 2 487 488 489def similarity_transformation(rot, mat): 490 """Similarity transformation by R x M x R^-1.""" 491 return np.dot(rot, np.dot(mat, np.linalg.inv(rot))) 492 493 494def show_drift_force_constants(force_constants, 495 primitive=None, 496 name="force constants", 497 values_only=False): 498 """Show force constants drift.""" 499 if force_constants.shape[0] == force_constants.shape[1]: 500 num_atom = force_constants.shape[0] 501 maxval1 = 0 502 maxval2 = 0 503 jk1 = [0, 0] 504 jk2 = [0, 0] 505 for i, j, k in list(np.ndindex((num_atom, 3, 3))): 506 val1 = force_constants[:, i, j, k].sum() 507 val2 = force_constants[i, :, j, k].sum() 508 if abs(val1) > abs(maxval1): 509 maxval1 = val1 510 jk1 = [j, k] 511 if abs(val2) > abs(maxval2): 512 maxval2 = val2 513 jk2 = [j, k] 514 else: 515 s2p_map = primitive.s2p_map 516 p2s_map = primitive.p2s_map 517 p2p_map = primitive.p2p_map 518 permutations = primitive.atomic_permutations 519 s2pp_map, nsym_list = get_nsym_list_and_s2pp(s2p_map, 520 p2p_map, 521 permutations) 522 523 try: 524 import phonopy._phonopy as phonoc 525 phonoc.transpose_compact_fc(force_constants, 526 permutations, 527 s2pp_map, 528 p2s_map, 529 nsym_list) 530 maxval1, jk1 = _get_drift_per_index(force_constants) 531 phonoc.transpose_compact_fc(force_constants, 532 permutations, 533 s2pp_map, 534 p2s_map, 535 nsym_list) 536 maxval2, jk2 = _get_drift_per_index(force_constants) 537 538 except ImportError: 539 text = ("Import error at phonoc.tranpose_compact_fc. " 540 "Corresponding python code is not implemented.") 541 raise RuntimeError(text) 542 543 if values_only: 544 text = "" 545 else: 546 text = "Max drift of %s: " % name 547 text += "%f (%s%s) %f (%s%s)" % (maxval1, "xyz"[jk1[0]], "xyz"[jk1[1]], 548 maxval2, "xyz"[jk2[0]], "xyz"[jk2[1]]) 549 print(text) 550 551 552def get_nsym_list_and_s2pp(s2p_map, 553 p2p_map, 554 permutations): 555 """Find lattice points corresponding to atoms in s2p_map. 556 557 Parameters 558 ---------- 559 s2p_map : array_like 560 See Primitive class. 561 p2p_map : dict 562 See Primitive class. 563 permutations : ndarray 564 See Primitive.atomic_permutations. 565 566 Returns 567 ------- 568 s2pp : ndarray 569 Atom indices in primitive cell that correspond to supercell atoms. 570 shape=(num_atoms_in_supercell, ), dtype='intc' 571 nsym_list : ndarray 572 Pure translation indices that map atoms in supercell to those in 573 primitive cell. 574 shape=(num_pure_translation, ), dtype='intc' 575 576 Note 577 ---- 578 This method is public because of being used by phono3py. 579 580 """ 581 s2pp = np.array([p2p_map[i] for i in s2p_map], dtype='intc') 582 nsym_list = np.array([np.where(permutations[:, i] == target)[0][0] 583 for i, target in enumerate(s2p_map)], dtype='intc') 584 return s2pp, nsym_list 585 586 587def get_harmonic_potential_energy(force_constants, displacements): 588 """Calculate harmonic potential energy of displacements. 589 590 Parameters 591 ---------- 592 force_constants : ndarray 593 Full shape force constants. 594 shape=(num_supercell_atoms, num_supercell_atoms, 3, 3), dtype='double' 595 displacements : ndarray 596 Displacements of atoms. 597 shape=(num_supercell_atoms, 3) or shape=(N, num_supercell_atoms, 3), 598 dtype='double' 599 N in shape means number of snapshots. 600 601 Returns 602 ------- 603 float or list of float 604 Increase of harmonic potential energy by displacements. 605 606 Note 607 ---- 608 This is not directly used in phonopy, but is kept useful. 609 610 """ 611 if force_constants.shape[0] != force_constants.shape[1]: 612 raise RuntimeError("Full shape force constants are necessary.") 613 614 def _get_harm_pot(fc, d): 615 return np.dot(d, np.dot(fc, d)) / 2 616 617 n = force_constants.shape[0] 618 fc = np.swapaxes(force_constants, 1, 2).reshape(n * 3, n * 3) 619 if displacements.ndim == 3: 620 return [_get_harm_pot(fc, d.ravel()) for d in displacements] 621 elif displacements.ndim == 2: 622 d = displacements.ravel() 623 return _get_harm_pot(fc, d) 624 else: 625 raise RuntimeError("Array shape of displacements is wrong.") 626 627 628def _get_rotated_forces(forces_syms, site_sym_cart): 629 rot_forces = [] 630 for forces, sym_cart in zip(forces_syms, site_sym_cart): 631 rot_forces.append(np.dot(forces, sym_cart.T)) 632 633 return rot_forces 634 635 636def _get_drift_per_index(force_constants): 637 num_atom = force_constants.shape[0] 638 maxval = 0 639 jk = [0, 0] 640 for i, j, k in list(np.ndindex((num_atom, 3, 3))): 641 val = force_constants[i, :, j, k].sum() 642 if abs(val) > abs(maxval): 643 maxval = val 644 jk = [j, k] 645 return maxval, jk 646 647 648def _solve_force_constants_svd(disp_atom_number, 649 displacements, 650 sets_of_forces, 651 supercell, 652 site_symmetry, 653 symprec): 654 lattice = supercell.get_cell().T 655 positions = supercell.get_scaled_positions() 656 pos_center = positions[disp_atom_number] 657 positions -= pos_center 658 rot_map_syms = get_positions_sent_by_rot_inv(lattice, 659 positions, 660 site_symmetry, 661 symprec) 662 site_sym_cart = [similarity_transformation(lattice, sym) 663 for sym in site_symmetry] 664 rot_disps = get_rotated_displacement(displacements, site_sym_cart) 665 inv_displacements = np.linalg.pinv(rot_disps) 666 667 fc = np.zeros((supercell.get_number_of_atoms(), 3, 3), 668 dtype='double', order='C') 669 for i in range(supercell.get_number_of_atoms()): 670 combined_forces = [] 671 for forces in sets_of_forces: 672 combined_forces.append( 673 _get_rotated_forces(forces[rot_map_syms[:, i]], 674 site_sym_cart)) 675 676 combined_forces = np.reshape(combined_forces, (-1, 3)) 677 fc[i] = -np.dot(inv_displacements, combined_forces) 678 return fc 679 680 681def _get_force_constants_disps(force_constants, 682 supercell, 683 dataset, 684 symmetry, 685 atom_list=None): 686 """Calculate force constants Phi = -F / d. 687 688 Force constants are obtained by one of the following algorithm. 689 690 Parameters 691 ---------- 692 force_constants: ndarray 693 Force constants 694 shape=(len(atom_list),n_satom,3,3) 695 dtype=double 696 supercell: Supercell 697 Supercell 698 dataset: dict 699 Distplacement dataset. Forces are also stored. 700 symmetry: Symmetry 701 Symmetry information of supercell 702 atom_list: list 703 List of atom indices corresponding to the first index of force 704 constants. None assigns all atoms in supercell. 705 706 """ 707 symprec = symmetry.tolerance 708 disp_atom_list = np.unique([x['number'] for x in dataset['first_atoms']]) 709 for disp_atom_number in disp_atom_list: 710 disps = [] 711 sets_of_forces = [] 712 713 for x in dataset['first_atoms']: 714 if x['number'] != disp_atom_number: 715 continue 716 disps.append(x['displacement']) 717 sets_of_forces.append(x['forces']) 718 719 site_symmetry = symmetry.get_site_symmetry(disp_atom_number) 720 721 solve_force_constants(force_constants, 722 disp_atom_number, 723 disps, 724 sets_of_forces, 725 supercell, 726 site_symmetry, 727 symprec, 728 atom_list=atom_list) 729 730 return disp_atom_list 731 732 733def _combine_force_constants_equivalent_atoms(fc_combined, 734 force_constants, 735 i, 736 cart_rot, 737 map_atoms, 738 mapa): 739 num_equiv_atoms = 0 740 for j, k in enumerate(map_atoms): 741 if k != i: 742 continue 743 744 num_equiv_atoms += 1 745 r_i = (mapa[:, j] == i).nonzero()[0][0] 746 for k, l in enumerate(mapa[r_i]): 747 fc_combined[l] += similarity_transformation( 748 cart_rot[r_i], force_constants[j, k]) 749 750 fc_combined /= num_equiv_atoms 751 752 return num_equiv_atoms 753 754 755def _average_force_constants_by_sitesym(fc_new, 756 fc_i, 757 i, 758 cart_rot, 759 mapa): 760 num_sitesym = 0 761 for r_i, mapa_r in enumerate(mapa): 762 if mapa_r[i] != i: 763 continue 764 num_sitesym += 1 765 for j in range(fc_i.shape[0]): 766 fc_new[i, j] += similarity_transformation( 767 cart_rot[r_i].T, fc_i[mapa[r_i, j]]) 768 769 fc_new[i] /= num_sitesym 770 771 return num_sitesym 772 773 774def _get_atom_indices_by_symmetry(lattice, 775 positions, 776 rotations, 777 translations, 778 symprec): 779 # To understand this method, understanding numpy broadcasting is mandatory. 780 781 K = len(positions) 782 # positions[K, 3] 783 # dot()[K, N, 3] where N is number of sym opts. 784 # translation[N, 3] is added to the last two dimenstions after dot(). 785 rpos = np.dot(positions, np.transpose(rotations, (0, 2, 1))) + translations 786 787 # np.tile(rpos, (K, 1, 1, 1))[K(2), K(1), N, 3] 788 # by adding one dimension in front of [K(1), N, 3]. 789 # np.transpose(np.tile(rpos, (K, 1, 1, 1)), (2, 1, 0, 3))[N, K(1), K(2), 3] 790 diff = positions - np.transpose(np.tile(rpos, (K, 1, 1, 1)), (2, 1, 0, 3)) 791 diff -= np.rint(diff) 792 diff = np.dot(diff, lattice.T) 793 # m[N, K(1), K(2)] 794 m = (np.sqrt(np.sum(diff ** 2, axis=3)) < symprec) 795 # index_array[K(1), K(2)] 796 index_array = np.tile(np.arange(K, dtype='intc'), (K, 1)) 797 # Understanding numpy boolean array indexing (extract True elements) 798 # mapa[N, K(1)] 799 mapa = np.array([index_array[mr] for mr in m]) 800 return mapa 801 802 803def _get_shortest_distance_in_PBC(pos_i, pos_j, reduced_bases): 804 distances = [] 805 for k in (-1, 0, 1): 806 for l in (-1, 0, 1): 807 for m in (-1, 0, 1): 808 diff = pos_j + np.array([k, l, m]) - pos_i 809 distances.append(np.linalg.norm(np.dot(diff, reduced_bases))) 810 return np.min(distances) 811 812 813def _get_sym_mappings_from_permutations(permutations, atom_list_done): 814 """Find atomic indices where force constants have not yet calculated. 815 816 This can be thought of as computing 'map_atom_disp' and 'map_sym' 817 for all atoms, except done using permutations instead of by 818 computing overlaps. 819 820 Parameters 821 ---------- 822 permutations : ndarray 823 Atomic index permutation table by space group operations. 824 shape=(operations, positions) 825 atom_list_done : array_like 826 Atomic indices where force constants (first index) were already 827 calculated. 828 829 Returns 830 ------- 831 map_atoms : ndarray 832 Maps each atom in the full structure to its equivalent atom in 833 atom_list_done. (each entry will be an integer found in 834 atom_list_done) 835 shape=(positions, ), dtype='intc' 836 map_syms : ndarray 837 For each atom, provides the index of a rotation that maps it 838 into atom_list_done. (there might be more than one such 839 rotation, but only one will be returned) (each entry will be 840 an integer 0 <= i < num_rot) 841 shape=(positions, ), dtype='intc' 842 843 """ 844 assert permutations.ndim == 2 845 num_pos = permutations.shape[1] 846 847 # filled with -1 848 map_atoms = np.zeros((num_pos,), dtype='intc') - 1 849 map_syms = np.zeros((num_pos,), dtype='intc') - 1 850 851 atom_list_done = set(atom_list_done) 852 for atom_todo in range(num_pos): 853 for (sym_index, permutation) in enumerate(permutations): 854 if permutation[atom_todo] in atom_list_done: 855 map_atoms[atom_todo] = permutation[atom_todo] 856 map_syms[atom_todo] = sym_index 857 break 858 else: 859 text = ("Input forces are not enough to calculate force constants," 860 "or something wrong (e.g. crystal structure does not " 861 "match).") 862 print(textwrap.fill(text)) 863 raise ValueError 864 865 assert set(map_atoms) & set(atom_list_done) == set(map_atoms) 866 assert -1 not in map_atoms 867 assert -1 not in map_syms 868 return map_atoms, map_syms 869