Coverage for ase / calculators / bond_polarizability.py: 100.00%

1# fmt: off 

2 

3 

4import numpy as np 

5 

6from ase.data import covalent_radii 

7from ase.neighborlist import NeighborList 

8from ase.units import Bohr, Ha 

9 

10from .polarizability import StaticPolarizabilityCalculator 

11 

12 

13class LippincottStuttman: 

14 # atomic polarizability values from: 

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

16 # DOI: 10.1021/j100792a033 

17 # see also: 

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

19 # DOI: 10.1103/PhysRevB.55.2938 

20 # unit: Angstrom^3 

21 atomic_polarizability = { 

22 'H': 0.592, 

23 'Be': 3.802, 

24 'B': 1.358, 

25 'C': 0.978, 

26 'N': 0.743, 

27 'O': 0.592, 

28 'Al': 3.918, 

29 'Si': 2.988, 

30 'P': 2.367, 

31 'S': 1.820, 

32 } 

33 # reduced electronegativity Table I 

34 reduced_eletronegativity = { 

35 'H': 1.0, 

36 'Be': 0.538, 

37 'B': 0.758, 

38 'C': 0.846, 

39 'N': 0.927, 

40 'O': 1.0, 

41 'Al': 0.533, 

42 'Si': 0.583, 

43 'P': 0.630, 

44 'S': 0.688, 

45 } 

46 

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

48 length: float) -> tuple[float, float]: 

49 """Bond polarizability 

50 

51 Parameters 

52 ---------- 

53 el1: element string 

54 el2: element string 

55 length: float 

56 

57 Returns 

58 ------- 

59 alphal: float 

60 Parallel component 

61 alphap: float 

62 Perpendicular component 

63 """ 

64 alpha1 = self.atomic_polarizability[el1] 

65 alpha2 = self.atomic_polarizability[el2] 

66 ren1 = self.reduced_eletronegativity[el1] 

67 ren2 = self.reduced_eletronegativity[el2] 

68 

69 sigma = 1. 

70 if el1 != el2: 

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

72 

73 # parallel component 

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

75 # XXX consider fractional covalency ? 

76 

77 # prependicular component 

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

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

80 # XXX consider fractional covalency ? 

81 

82 return alphal, alphap 

83 

84 

85class Linearized: 

86 def __init__(self): 

87 self._data = { 

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

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

90 # R0 al al' ap ap' 

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

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

93 } 

94 

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

96 length: float) -> tuple[float, float]: 

97 """Bond polarizability 

98 

99 Parameters 

100 ---------- 

101 el1: element string 

102 el2: element string 

103 length: float 

104 

105 Returns 

106 ------- 

107 alphal: float 

108 Parallel component 

109 alphap: float 

110 Perpendicular component 

111 """ 

112 if el1 > el2: 

113 bond = el2 + el1 

114 else: 

115 bond = el1 + el2 

116 assert bond in self._data 

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

118 

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

120 

121 

122class BondPolarizability(StaticPolarizabilityCalculator): 

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

124 self.model = model 

125 

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

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

128 

129 Parameters 

130 ---------- 

131 atoms: Atoms object 

132 radiicut: float 

133 Bonds are counted up to 

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

135 Default: 1.5 

136 

137 Returns 

138 ------- 

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

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

141 """ 

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

143 for z in atoms.numbers]) 

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

145 self_interaction=False) 

146 nl.update(atoms) 

147 pos_ac = atoms.get_positions() 

148 

149 alpha = 0 

150 for ia, atom in enumerate(atoms): 

151 indices, offsets = nl.get_neighbors(ia) 

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

153 

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

155 weight = 1 

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

157 weight = 0.5 # count half the bond only 

158 

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

160 dist = np.linalg.norm(dist_c) 

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

162 

163 eye3 = np.eye(3) / 3 

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

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

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

167 return alpha / Bohr / Ha