1from typing import Tuple 2import numpy as np 3 4from ase.units import Bohr, Ha 5from ase.data import covalent_radii 6from ase.neighborlist import NeighborList 7 8 9class LippincottStuttman: 10 # atomic polarizability values from: 11 # Lippincott and Stutman J. Phys. Chem. 68 (1964) 2926-2940 12 # DOI: 10.1021/j100792a033 13 # see also: 14 # Marinov and Zotov Phys. Rev. B 55 (1997) 2938-2944 15 # DOI: 10.1103/PhysRevB.55.2938 16 # unit: Angstrom^3 17 atomic_polarizability = { 18 'B': 1.358, 19 'C': 0.978, 20 'N': 0.743, 21 'O': 0.592, 22 'Al': 3.918, 23 'Si': 2.988, 24 } 25 # reduced electronegativity Table I 26 reduced_eletronegativity = { 27 'B': 0.538, 28 'C': 0.846, 29 'N': 0.927, 30 'O': 1.0, 31 'Al': 0.533, 32 'Si': 0.583, 33 } 34 35 def __call__(self, el1: str, el2: str, 36 length: float) -> Tuple[float, float]: 37 """Bond polarizability 38 39 Parameters 40 ---------- 41 el1: element string 42 el2: element string 43 length: float 44 45 Returns 46 ------- 47 alphal: float 48 Parallel component 49 alphap: float 50 Perpendicular component 51 """ 52 alpha1 = self.atomic_polarizability[el1] 53 alpha2 = self.atomic_polarizability[el2] 54 ren1 = self.reduced_eletronegativity[el1] 55 ren2 = self.reduced_eletronegativity[el2] 56 57 sigma = 1. 58 if el1 != el2: 59 sigma = np.exp(- (ren1 - ren2)**2 / 4) 60 61 # parallel component 62 alphal = sigma * length**4 / (4**4 * alpha1 * alpha2)**(1. / 6) 63 # XXX consider fractional covalency ? 64 65 # prependicular component 66 alphap = ((ren1**2 * alpha1 + ren2**2 * alpha2) 67 / (ren1**2 + ren2**2)) 68 # XXX consider fractional covalency ? 69 70 return alphal, alphap 71 72 73class Linearized: 74 def __init__(self): 75 self._data = { 76 # L. Wirtz, M. Lazzeri, F. Mauri, A. Rubio, 77 # Phys. Rev. B 2005, 71, 241402. 78 # R0 al al' ap ap' 79 'CC': (1.53, 1.69, 7.43, 0.71, 0.37), 80 'BN': (1.56, 1.58, 4.22, 0.42, 0.90), 81 } 82 83 def __call__(self, el1: str, el2: str, 84 length: float) -> Tuple[float, float]: 85 """Bond polarizability 86 87 Parameters 88 ---------- 89 el1: element string 90 el2: element string 91 length: float 92 93 Returns 94 ------- 95 alphal: float 96 Parallel component 97 alphap: float 98 Perpendicular component 99 """ 100 if el1 > el2: 101 bond = el2 + el1 102 else: 103 bond = el1 + el2 104 assert bond in self._data 105 length0, al, ald, ap, apd = self._data[bond] 106 107 return al + ald * (length - length0), ap + apd * (length - length0) 108 109 110class BondPolarizability: 111 def __init__(self, model=LippincottStuttman()): 112 self.model = model 113 114 def __call__(self, *args, **kwargs): 115 """Shorthand for calculate""" 116 return self.calculate(*args, **kwargs) 117 118 def calculate(self, atoms, radiicut=1.5): 119 """Sum up the bond polarizability from all bonds 120 121 Parameters 122 ---------- 123 atoms: Atoms object 124 radiicut: float 125 Bonds are counted up to 126 radiicut * (sum of covalent radii of the pairs) 127 Default: 1.5 128 129 Returns 130 ------- 131 polarizability tensor with unit (e^2 Angstrom^2 / eV). 132 Multiply with Bohr * Ha to get (Angstrom^3) 133 """ 134 radii = np.array([covalent_radii[z] 135 for z in atoms.numbers]) 136 nl = NeighborList(radii * 1.5, skin=0, 137 self_interaction=False) 138 nl.update(atoms) 139 pos_ac = atoms.get_positions() 140 141 alpha = 0 142 for ia, atom in enumerate(atoms): 143 indices, offsets = nl.get_neighbors(ia) 144 pos_ac = atoms.get_positions() - atoms.get_positions()[ia] 145 146 for ib, offset in zip(indices, offsets): 147 weight = 1 148 if offset.any(): # this comes from a periodic image 149 weight = 0.5 # count half the bond only 150 151 dist_c = pos_ac[ib] + np.dot(offset, atoms.get_cell()) 152 dist = np.linalg.norm(dist_c) 153 al, ap = self.model(atom.symbol, atoms[ib].symbol, dist) 154 155 eye3 = np.eye(3) / 3 156 alpha += weight * (al + 2 * ap) * eye3 157 alpha += weight * (al - ap) * ( 158 np.outer(dist_c, dist_c) / dist**2 - eye3) 159 return alpha / Bohr / Ha 160