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

1from typing import Tuple 

2 

3import numpy as np 

4 

5from ase.data import covalent_radii 

6from ase.neighborlist import NeighborList 

7from ase.units import Bohr, Ha 

8 

9from .polarizability import StaticPolarizabilityCalculator 

10 

11 

12class LippincottStuttman: 

13 # atomic polarizability values from: 

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

15 # DOI: 10.1021/j100792a033 

16 # see also: 

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

18 # DOI: 10.1103/PhysRevB.55.2938 

19 # unit: Angstrom^3 

20 atomic_polarizability = { 

21 'H': 0.592, 

22 'Be': 3.802, 

23 'B': 1.358, 

24 'C': 0.978, 

25 'N': 0.743, 

26 'O': 0.592, 

27 'Al': 3.918, 

28 'Si': 2.988, 

29 'P': 2.367, 

30 'S': 1.820, 

31 } 

32 # reduced electronegativity Table I 

33 reduced_eletronegativity = { 

34 'H': 1.0, 

35 'Be': 0.538, 

36 'B': 0.758, 

37 'C': 0.846, 

38 'N': 0.927, 

39 'O': 1.0, 

40 'Al': 0.533, 

41 'Si': 0.583, 

42 'P': 0.630, 

43 'S': 0.688, 

44 } 

45 

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

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

48 """Bond polarizability 

49 

50 Parameters 

51 ---------- 

52 el1: element string 

53 el2: element string 

54 length: float 

55 

56 Returns 

57 ------- 

58 alphal: float 

59 Parallel component 

60 alphap: float 

61 Perpendicular component 

62 """ 

63 alpha1 = self.atomic_polarizability[el1] 

64 alpha2 = self.atomic_polarizability[el2] 

65 ren1 = self.reduced_eletronegativity[el1] 

66 ren2 = self.reduced_eletronegativity[el2] 

67 

68 sigma = 1. 

69 if el1 != el2: 

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

71 

72 # parallel component 

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

74 # XXX consider fractional covalency ? 

75 

76 # prependicular component 

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

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

79 # XXX consider fractional covalency ? 

80 

81 return alphal, alphap 

82 

83 

84class Linearized: 

85 def __init__(self): 

86 self._data = { 

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

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

89 # R0 al al' ap ap' 

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

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

92 } 

93 

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

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

96 """Bond polarizability 

97 

98 Parameters 

99 ---------- 

100 el1: element string 

101 el2: element string 

102 length: float 

103 

104 Returns 

105 ------- 

106 alphal: float 

107 Parallel component 

108 alphap: float 

109 Perpendicular component 

110 """ 

111 if el1 > el2: 

112 bond = el2 + el1 

113 else: 

114 bond = el1 + el2 

115 assert bond in self._data 

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

117 

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

119 

120 

121class BondPolarizability(StaticPolarizabilityCalculator): 

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

123 self.model = model 

124 

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

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

127 

128 Parameters 

129 ---------- 

130 atoms: Atoms object 

131 radiicut: float 

132 Bonds are counted up to 

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

134 Default: 1.5 

135 

136 Returns 

137 ------- 

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

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

140 """ 

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

142 for z in atoms.numbers]) 

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

144 self_interaction=False) 

145 nl.update(atoms) 

146 pos_ac = atoms.get_positions() 

147 

148 alpha = 0 

149 for ia, atom in enumerate(atoms): 

150 indices, offsets = nl.get_neighbors(ia) 

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

152 

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

154 weight = 1 

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

156 weight = 0.5 # count half the bond only 

157 

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

159 dist = np.linalg.norm(dist_c) 

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

161 

162 eye3 = np.eye(3) / 3 

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

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

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

166 return alpha / Bohr / Ha