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

1# fmt: off 

2 

3from typing import Tuple 

4 

5import numpy as np 

6 

7from ase.data import covalent_radii 

8from ase.neighborlist import NeighborList 

9from ase.units import Bohr, Ha 

10 

11from .polarizability import StaticPolarizabilityCalculator 

12 

13 

14class LippincottStuttman: 

15 # atomic polarizability values from: 

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

17 # DOI: 10.1021/j100792a033 

18 # see also: 

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

20 # DOI: 10.1103/PhysRevB.55.2938 

21 # unit: Angstrom^3 

22 atomic_polarizability = { 

23 'H': 0.592, 

24 'Be': 3.802, 

25 'B': 1.358, 

26 'C': 0.978, 

27 'N': 0.743, 

28 'O': 0.592, 

29 'Al': 3.918, 

30 'Si': 2.988, 

31 'P': 2.367, 

32 'S': 1.820, 

33 } 

34 # reduced electronegativity Table I 

35 reduced_eletronegativity = { 

36 'H': 1.0, 

37 'Be': 0.538, 

38 'B': 0.758, 

39 'C': 0.846, 

40 'N': 0.927, 

41 'O': 1.0, 

42 'Al': 0.533, 

43 'Si': 0.583, 

44 'P': 0.630, 

45 'S': 0.688, 

46 } 

47 

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

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

50 """Bond polarizability 

51 

52 Parameters 

53 ---------- 

54 el1: element string 

55 el2: element string 

56 length: float 

57 

58 Returns 

59 ------- 

60 alphal: float 

61 Parallel component 

62 alphap: float 

63 Perpendicular component 

64 """ 

65 alpha1 = self.atomic_polarizability[el1] 

66 alpha2 = self.atomic_polarizability[el2] 

67 ren1 = self.reduced_eletronegativity[el1] 

68 ren2 = self.reduced_eletronegativity[el2] 

69 

70 sigma = 1. 

71 if el1 != el2: 

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

73 

74 # parallel component 

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

76 # XXX consider fractional covalency ? 

77 

78 # prependicular component 

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

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

81 # XXX consider fractional covalency ? 

82 

83 return alphal, alphap 

84 

85 

86class Linearized: 

87 def __init__(self): 

88 self._data = { 

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

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

91 # R0 al al' ap ap' 

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

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

94 } 

95 

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

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

98 """Bond polarizability 

99 

100 Parameters 

101 ---------- 

102 el1: element string 

103 el2: element string 

104 length: float 

105 

106 Returns 

107 ------- 

108 alphal: float 

109 Parallel component 

110 alphap: float 

111 Perpendicular component 

112 """ 

113 if el1 > el2: 

114 bond = el2 + el1 

115 else: 

116 bond = el1 + el2 

117 assert bond in self._data 

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

119 

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

121 

122 

123class BondPolarizability(StaticPolarizabilityCalculator): 

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

125 self.model = model 

126 

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

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

129 

130 Parameters 

131 ---------- 

132 atoms: Atoms object 

133 radiicut: float 

134 Bonds are counted up to 

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

136 Default: 1.5 

137 

138 Returns 

139 ------- 

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

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

142 """ 

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

144 for z in atoms.numbers]) 

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

146 self_interaction=False) 

147 nl.update(atoms) 

148 pos_ac = atoms.get_positions() 

149 

150 alpha = 0 

151 for ia, atom in enumerate(atoms): 

152 indices, offsets = nl.get_neighbors(ia) 

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

154 

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

156 weight = 1 

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

158 weight = 0.5 # count half the bond only 

159 

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

161 dist = np.linalg.norm(dist_c) 

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

163 

164 eye3 = np.eye(3) / 3 

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

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

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

168 return alpha / Bohr / Ha