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

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

5logger = logger.getChild(__name__)

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.

13 Configurations are generated by superimposing harmonic phonon

14 eigenmodes with random amplitudes and phase factors consistent

15 with a certain temperature.

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

22 .. math::

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)

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

32 or

34 .. math::

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)

43 where

45 .. math::

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

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)

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

80def _n_BE(T, w_s):

81 """

82 Bose-Einstein distribution function.

84 Parameters

85 ---------

86 T : float

87 Temperature in Kelvin

88 w_s: numpy.ndarray

89 frequencies in eV (3*N,)

91 Returns

92 ------

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

94 """

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

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).

108 _s is a mode index

109 _i is a Carteesian index

110 _a is an atom index

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

126 Returns

127 -------

128 displacements : numpy.ndarray

129 shape (N, 3)

130 """

131 n_modes = 3 * len(m_a)

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:]

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

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)

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)))

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

156class _PhononRattler:

157 """

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

159 for phonon rattle.

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)

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}')

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')

198 # treat imaginary modes as real

199 imag_mask = w2_s < -frequency_tol

202 self.w2_s = w2_s

203 self.e_sai = e_sai

204 self.masses = masses

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

207 """ rattle atoms by adding displacements

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