Coverage for hiphive/structure_generation/phonon.py: 94%
61 statements
« prev ^ index » next coverage.py v7.6.8, created at 2024-11-28 11:20 +0000
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
200 w2_s[imag_mask] = imag_freq_factor * np.abs(w2_s[imag_mask])
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