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 

7from .polarizability import StaticPolarizabilityCalculator 

8 

9 

10class LippincottStuttman: 

11 # atomic polarizability values from: 

12 # Lippincott and Stutman J. Phys. Chem. 68 (1964) 2926-2940 

13 # DOI: 10.1021/j100792a033 

14 # see also: 

15 # Marinov and Zotov Phys. Rev. B 55 (1997) 2938-2944 

16 # DOI: 10.1103/PhysRevB.55.2938 

17 # unit: Angstrom^3 

18 atomic_polarizability = { 

19 'H': 0.592, 

20 'Be': 3.802, 

21 'B': 1.358, 

22 'C': 0.978, 

23 'N': 0.743, 

24 'O': 0.592, 

25 'Al': 3.918, 

26 'Si': 2.988, 

27 'P': 2.367, 

28 'S': 1.820, 

29 } 

30 # reduced electronegativity Table I 

31 reduced_eletronegativity = { 

32 'H': 1.0, 

33 'Be': 0.538, 

34 'B': 0.758, 

35 'C': 0.846, 

36 'N': 0.927, 

37 'O': 1.0, 

38 'Al': 0.533, 

39 'Si': 0.583, 

40 'P': 0.630, 

41 'S': 0.688, 

42 } 

43 

44 def __call__(self, el1: str, el2: str, 

45 length: float) -> Tuple[float, float]: 

46 """Bond polarizability 

47 

48 Parameters 

49 ---------- 

50 el1: element string 

51 el2: element string 

52 length: float 

53 

54 Returns 

55 ------- 

56 alphal: float 

57 Parallel component 

58 alphap: float 

59 Perpendicular component 

60 """ 

61 alpha1 = self.atomic_polarizability[el1] 

62 alpha2 = self.atomic_polarizability[el2] 

63 ren1 = self.reduced_eletronegativity[el1] 

64 ren2 = self.reduced_eletronegativity[el2] 

65 

66 sigma = 1. 

67 if el1 != el2: 

68 sigma = np.exp(- (ren1 - ren2)**2 / 4) 

69 

70 # parallel component 

71 alphal = sigma * length**4 / (4**4 * alpha1 * alpha2)**(1. / 6) 

72 # XXX consider fractional covalency ? 

73 

74 # prependicular component 

75 alphap = ((ren1**2 * alpha1 + ren2**2 * alpha2) 

76 / (ren1**2 + ren2**2)) 

77 # XXX consider fractional covalency ? 

78 

79 return alphal, alphap 

80 

81 

82class Linearized: 

83 def __init__(self): 

84 self._data = { 

85 # L. Wirtz, M. Lazzeri, F. Mauri, A. Rubio, 

86 # Phys. Rev. B 2005, 71, 241402. 

87 # R0 al al' ap ap' 

88 'CC': (1.53, 1.69, 7.43, 0.71, 0.37), 

89 'BN': (1.56, 1.58, 4.22, 0.42, 0.90), 

90 } 

91 

92 def __call__(self, el1: str, el2: str, 

93 length: float) -> Tuple[float, float]: 

94 """Bond polarizability 

95 

96 Parameters 

97 ---------- 

98 el1: element string 

99 el2: element string 

100 length: float 

101 

102 Returns 

103 ------- 

104 alphal: float 

105 Parallel component 

106 alphap: float 

107 Perpendicular component 

108 """ 

109 if el1 > el2: 

110 bond = el2 + el1 

111 else: 

112 bond = el1 + el2 

113 assert bond in self._data 

114 length0, al, ald, ap, apd = self._data[bond] 

115 

116 return al + ald * (length - length0), ap + apd * (length - length0) 

117 

118 

119class BondPolarizability(StaticPolarizabilityCalculator): 

120 def __init__(self, model=LippincottStuttman()): 

121 self.model = model 

122 

123 def __call__(self, atoms, radiicut=1.5): 

124 """Sum up the bond polarizability from all bonds 

125 

126 Parameters 

127 ---------- 

128 atoms: Atoms object 

129 radiicut: float 

130 Bonds are counted up to 

131 radiicut * (sum of covalent radii of the pairs) 

132 Default: 1.5 

133 

134 Returns 

135 ------- 

136 polarizability tensor with unit (e^2 Angstrom^2 / eV). 

137 Multiply with Bohr * Ha to get (Angstrom^3) 

138 """ 

139 radii = np.array([covalent_radii[z] 

140 for z in atoms.numbers]) 

141 nl = NeighborList(radii * 1.5, skin=0, 

142 self_interaction=False) 

143 nl.update(atoms) 

144 pos_ac = atoms.get_positions() 

145 

146 alpha = 0 

147 for ia, atom in enumerate(atoms): 

148 indices, offsets = nl.get_neighbors(ia) 

149 pos_ac = atoms.get_positions() - atoms.get_positions()[ia] 

150 

151 for ib, offset in zip(indices, offsets): 

152 weight = 1 

153 if offset.any(): # this comes from a periodic image 

154 weight = 0.5 # count half the bond only 

155 

156 dist_c = pos_ac[ib] + np.dot(offset, atoms.get_cell()) 

157 dist = np.linalg.norm(dist_c) 

158 al, ap = self.model(atom.symbol, atoms[ib].symbol, dist) 

159 

160 eye3 = np.eye(3) / 3 

161 alpha += weight * (al + 2 * ap) * eye3 

162 alpha += weight * (al - ap) * ( 

163 np.outer(dist_c, dist_c) / dist**2 - eye3) 

164 return alpha / Bohr / Ha