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