Coverage for hiphive/structure_generation/phonon.py: 95%

Shortcuts on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

61 statements  

1import numpy as np 

2import ase.units as aunits 

3from hiphive.input_output.logging_tools import logger 

4 

5logger = logger.getChild(__name__) 

6 

7 

8def generate_phonon_rattled_structures(atoms, fc2, n_structures, temperature, 

9 QM_statistics=False, imag_freq_factor=1.0): 

10 """ 

11 Returns list of phonon-rattled configurations. 

12 

13 Configurations are generated by superimposing harmonic phonon 

14 eigenmodes with random amplitudes and phase factors consistent 

15 with a certain temperature. 

16 

17 Let :math:`\\boldsymbol{X}_{ai}` be the phonon modes indexed by atom 

18 :math:`a` and mode :math:`i`, :math:`\\omega_i` the phonon 

19 frequencies, and let :math:`0 < Q_i \\leq 1` and :math:`0 \\leq 

20 U_i < 1` be uniformly random numbers. Then 

21 

22 .. math:: 

23 

24 \\boldsymbol{R}_a 

25 = \\boldsymbol{R}^0_a 

26 + \\left<\\frac{k_B T}{m_a} \\right>^{1/2} 

27 \\sum_i \\frac{1}{\\omega_i} \\boldsymbol{X}_{ai} 

28 \\sqrt{-2 \\ln Q_i} \\cos(\\pi \\omega_i U_i) 

29 

30 See: West and Estreicher, Physical Review Letters **96**, 115504 (2006) 

31 

32 or 

33 

34 .. math:: 

35 

36 \\boldsymbol{R}_a 

37 = \\boldsymbol{R}^0_a 

38 + \\sum_i 

39 \\left(\\frac{\\hbar(0.5 + n(\\omega_i, T))}{m_a\\omega_i} \\right)^{1/2} 

40 \\boldsymbol{X}_{ai} 

41 \\sqrt{-2 \\ln Q_i} \\cos(\\pi \\omega_i U_i) 

42 

43 where 

44 

45 .. math:: 

46 n = \\frac{1}{e^{\\hbar\\omega_i/k_BT}-1} 

47 

48 Parameters 

49 ---------- 

50 atoms : ase.Atoms 

51 prototype structure 

52 fc2 : numpy.ndarray 

53 second order force constant matrix, with shape `(3N, 3N)` or 

54 `(N, N, 3, 3)`. The conversion will be done internally if. 

55 n_structures : int 

56 number of structures to generate 

57 temperature : float 

58 temperature in Kelvin 

59 QM_statistics : bool 

60 if the amplitude of the quantum harmonic oscillator shoule be used 

61 instead of the classical amplitude 

62 imag_freq_factor: float 

63 If a squared frequency, w2, is negative then it is set to 

64 w2 = imag_freq_factor * np.abs(w2) 

65 

66 Returns 

67 ------- 

68 list(ase.Atoms) 

69 generated structures 

70 """ 

71 structures = [] 

72 pr = _PhononRattler(atoms.get_masses(), fc2, imag_freq_factor) 

73 for _ in range(n_structures): 

74 atoms_tmp = atoms.copy() 

75 pr(atoms_tmp, temperature, QM_statistics) 

76 structures.append(atoms_tmp) 

77 return structures 

78 

79 

80def _n_BE(T, w_s): 

81 """ 

82 Bose-Einstein distribution function. 

83 

84 Parameters 

85 --------- 

86 T : float 

87 Temperature in Kelvin 

88 w_s: numpy.ndarray 

89 frequencies in eV (3*N,) 

90 

91 Returns 

92 ------ 

93 Bose-Einstein distribution for each energy at a given temperature 

94 """ 

95 

96 with np.errstate(divide="raise", over="raise"): 

97 try: 

98 n = 1 / (np.exp(w_s / (aunits.kB * T)) - 1) 

99 except Exception: 

100 n = np.zeros_like(w_s) 

101 return n 

102 

103 

104def _phonon_rattle(m_a, T, w2_s, e_sai, QM_statistics): 

105 """ Thermal excitation of phonon modes as described by West and 

106 Estreicher, Physical Review Letters **96**, 115504 (2006). 

107 

108 _s is a mode index 

109 _i is a Carteesian index 

110 _a is an atom index 

111 

112 Parameters 

113 ---------- 

114 m_a : numpy.ndarray 

115 masses (N,) 

116 T : float 

117 temperature in Kelvin 

118 w2_s : numpy.ndarray 

119 the squared frequencies from the eigenvalue problem (3*N,) 

120 e_sai : numpy.ndarray 

121 polarizations (3*N, N, 3) 

122 QM_statistics : bool 

123 if the amplitude of the quantum harmonic oscillator shoule be used 

124 instead of the classical amplitude 

125 

126 Returns 

127 ------- 

128 displacements : numpy.ndarray 

129 shape (N, 3) 

130 """ 

131 n_modes = 3 * len(m_a) 

132 

133 # skip 3 translational modes 

134 argsort = np.argsort(np.abs(w2_s)) 

135 e_sai = e_sai[argsort][3:] 

136 w2_s = w2_s[argsort][3:] 

137 

138 w_s = np.sqrt(np.abs(w2_s)) 

139 

140 prefactor_a = np.sqrt(1 / m_a).reshape(-1, 1) 

141 if QM_statistics: 

142 hbar = aunits._hbar * aunits.J * aunits.s 

143 frequencyfactor_s = np.sqrt(hbar * (0.5 + _n_BE(T, hbar * w_s)) / w_s) 

144 else: 

145 frequencyfactor_s = 1 / w_s 

146 prefactor_a *= np.sqrt(aunits.kB * T) 

147 

148 phases_s = np.random.uniform(0, 2 * np.pi, size=n_modes - 3) 

149 amplitudes_s = np.sqrt(-2 * np.log(1 - np.random.random(n_modes - 3))) 

150 

151 u_ai = prefactor_a * np.tensordot( 

152 amplitudes_s * np.cos(phases_s) * frequencyfactor_s, e_sai, (0, 0)) 

153 return u_ai # displacements 

154 

155 

156class _PhononRattler: 

157 """ 

158 Class to be able to conveniently save modes and frequencies needed 

159 for phonon rattle. 

160 

161 Parameters 

162 ---------- 

163 masses : numpy.ndarray 

164 masses (N,) 

165 force_constants : numpy.ndarray 

166 second order force constant matrix, with shape `(3N, 3N)` or 

167 `(N, N, 3, 3)`. The conversion will be done internally if. 

168 imag_freq_factor: float 

169 If a squared frequency, w2, is negative then it is set to 

170 w2 = imag_freq_factor * np.abs(w2) 

171 """ 

172 def __init__(self, masses, force_constants, imag_freq_factor=1.0): 

173 n_atoms = len(masses) 

174 if len(force_constants.shape) == 4: # assume shape = (n_atoms, n_atoms, 3, 3) 

175 force_constants = force_constants.transpose(0, 2, 1, 3) 

176 force_constants = force_constants.reshape(3 * n_atoms, 3 * n_atoms) 

177 # Now the fc should have shape = (n_atoms * 3, n_atoms * 3) 

178 # Construct the dynamical matrix 

179 inv_root_masses = (1 / np.sqrt(masses)).repeat(3) 

180 D = np.outer(inv_root_masses, inv_root_masses) 

181 D *= force_constants 

182 # find frequnecies and eigenvectors 

183 w2_s, e_sai = np.linalg.eigh(D) 

184 # reshape to get atom index and Cartesian index separate 

185 e_sai = e_sai.T.reshape(-1, n_atoms, 3) 

186 

187 # The three modes closest to zero are assumed to be zero, ie acoustic sum rules are assumed 

188 frequency_tol = 1e-6 

189 argsort = np.argsort(np.abs(w2_s)) 

190 w2_gamma = w2_s[argsort][:3] 

191 if np.any(np.abs(w2_gamma) > frequency_tol): 191 ↛ 192line 191 didn't jump to line 192, because the condition on line 191 was never true

192 logger.warning(f'Acoustic sum rules not enforced, squared frequencies: {w2_gamma}') 

193 

194 # warning if any imaginary modes 

195 if np.any(w2_s < -frequency_tol): 195 ↛ 196line 195 didn't jump to line 196, because the condition on line 195 was never true

196 logger.warning('Imaginary modes present') 

197 

198 # treat imaginary modes as real 

199 imag_mask = w2_s < -frequency_tol 

200 w2_s[imag_mask] = imag_freq_factor * np.abs(w2_s[imag_mask]) 

201 

202 self.w2_s = w2_s 

203 self.e_sai = e_sai 

204 self.masses = masses 

205 

206 def __call__(self, atoms, T, QM_statistics): 

207 """ rattle atoms by adding displacements 

208 

209 Parameters 

210 ---------- 

211 atoms : ase.Atoms 

212 Ideal structure to add displacements to. 

213 T : float 

214 temperature in Kelvin 

215 """ 

216 u_ai = _phonon_rattle(self.masses, T, self.w2_s, self.e_sai, 

217 QM_statistics) 

218 atoms.positions += u_ai