Coverage for hiphive/structure_generation/phonon.py: 94%
62 statements
« prev ^ index » next coverage.py v7.10.1, created at 2025-08-01 17:04 +0000
1import numpy as np
2from ase import Atoms
3import ase.units as aunits
4from hiphive.input_output.logging_tools import logger
6logger = logger.getChild(__name__)
9def generate_phonon_rattled_structures(
10 atoms: Atoms,
11 fc2: np.ndarray,
12 n_structures: int,
13 temperature: float,
14 QM_statistics: bool = False,
15 imag_freq_factor: float = 1.0,
16) -> list[Atoms]:
17 r"""
18 Returns list of phonon-rattled configurations.
20 Configurations are generated by superimposing harmonic phonon
21 eigenmodes with random amplitudes and phase factors consistent
22 with a certain temperature.
24 Let :math:`\boldsymbol{X}_{ai}` be the phonon modes indexed by atom
25 :math:`a` and mode :math:`i`, :math:`\omega_i` the phonon
26 frequencies, and let :math:`0 < Q_i \leq 1` and :math:`0 \leq
27 U_i < 1` be uniformly random numbers. Then
29 .. math::
31 \boldsymbol{R}_a
32 = \boldsymbol{R}^0_a
33 + \left<\frac{k_B T}{m_a} \right>^{1/2}
34 \sum_i \frac{1}{\omega_i} \boldsymbol{X}_{ai}
35 \sqrt{-2 \ln Q_i} \cos(\pi \omega_i U_i)
37 See: West and Estreicher, Physical Review Letters **96**, 115504 (2006)
39 or
41 .. math::
43 \boldsymbol{R}_a
44 = \boldsymbol{R}^0_a
45 + \sum_i
46 \left(\frac{\hbar(0.5 + n(\omega_i, T))}{m_a\omega_i} \right)^{1/2}
47 \boldsymbol{X}_{ai}
48 \sqrt{-2 \ln Q_i} \cos(\pi \omega_i U_i)
50 where
52 .. math::
53 n = \frac{1}{e^{\hbar\omega_i/k_BT}-1}
55 Parameters
56 ----------
57 atoms
58 Prototype structure.
59 fc2
60 Second order force constant matrix, with shape `(3N, 3N)` or
61 `(N, N, 3, 3)`. The conversion will be done internally if.
62 n_structures
63 Number of structures to generate.
64 temperature
65 Temperature in Kelvin.
66 QM_statistics
67 If the amplitude of the quantum harmonic oscillator shoule be used
68 instead of the classical amplitude.
69 imag_freq_factor
70 If a squared frequency, `w2`, is negative it is set to
71 `w2 = imag_freq_factor * np.abs(w2)`.
73 Returns
74 -------
75 Generated structures.
76 """
77 structures = []
78 pr = _PhononRattler(atoms.get_masses(), fc2, imag_freq_factor)
79 for _ in range(n_structures):
80 atoms_tmp = atoms.copy()
81 pr(atoms_tmp, temperature, QM_statistics)
82 structures.append(atoms_tmp)
83 return structures
86def _n_BE(T, w_s):
87 """
88 Bose-Einstein distribution function.
90 Parameters
91 ---------
92 T : float
93 Temperature in Kelvin
94 w_s: numpy.ndarray
95 frequencies in eV (3*N,)
97 Returns
98 ------
99 Bose-Einstein distribution for each energy at a given temperature
100 """
102 with np.errstate(divide='raise', over='raise'):
103 try:
104 n = 1 / (np.exp(w_s / (aunits.kB * T)) - 1)
105 except Exception:
106 n = np.zeros_like(w_s)
107 return n
110def _phonon_rattle(m_a, T, w2_s, e_sai, QM_statistics):
111 """ Thermal excitation of phonon modes as described by West and
112 Estreicher, Physical Review Letters **96**, 115504 (2006).
114 _s is a mode index
115 _i is a Carteesian index
116 _a is an atom index
118 Parameters
119 ----------
120 m_a : numpy.ndarray
121 masses (N,)
122 T : float
123 temperature in Kelvin
124 w2_s : numpy.ndarray
125 the squared frequencies from the eigenvalue problem (3*N,)
126 e_sai : numpy.ndarray
127 polarizations (3*N, N, 3)
128 QM_statistics : bool
129 if the amplitude of the quantum harmonic oscillator shoule be used
130 instead of the classical amplitude
132 Returns
133 -------
134 displacements : numpy.ndarray
135 shape (N, 3)
136 """
137 n_modes = 3 * len(m_a)
139 # skip 3 translational modes
140 argsort = np.argsort(np.abs(w2_s))
141 e_sai = e_sai[argsort][3:]
142 w2_s = w2_s[argsort][3:]
144 w_s = np.sqrt(np.abs(w2_s))
146 prefactor_a = np.sqrt(1 / m_a).reshape(-1, 1)
147 if QM_statistics:
148 hbar = aunits._hbar * aunits.J * aunits.s
149 frequencyfactor_s = np.sqrt(hbar * (0.5 + _n_BE(T, hbar * w_s)) / w_s)
150 else:
151 frequencyfactor_s = 1 / w_s
152 prefactor_a *= np.sqrt(aunits.kB * T)
154 phases_s = np.random.uniform(0, 2 * np.pi, size=n_modes - 3)
155 amplitudes_s = np.sqrt(-2 * np.log(1 - np.random.random(n_modes - 3)))
157 u_ai = prefactor_a * np.tensordot(
158 amplitudes_s * np.cos(phases_s) * frequencyfactor_s, e_sai, (0, 0))
159 return u_ai # displacements
162class _PhononRattler:
163 """
164 Class to be able to conveniently save modes and frequencies needed
165 for phonon rattle.
167 Parameters
168 ----------
169 masses : numpy.ndarray
170 masses (N,)
171 force_constants : numpy.ndarray
172 second order force constant matrix, with shape `(3N, 3N)` or
173 `(N, N, 3, 3)`. The conversion will be done internally if.
174 imag_freq_factor: float
175 If a squared frequency, w2, is negative then it is set to
176 w2 = imag_freq_factor * np.abs(w2)
177 """
178 def __init__(self, masses, force_constants, imag_freq_factor=1.0):
179 n_atoms = len(masses)
180 if len(force_constants.shape) == 4: # assume shape = (n_atoms, n_atoms, 3, 3)
181 force_constants = force_constants.transpose(0, 2, 1, 3)
182 force_constants = force_constants.reshape(3 * n_atoms, 3 * n_atoms)
183 # Now the fc should have shape = (n_atoms * 3, n_atoms * 3)
184 # Construct the dynamical matrix
185 inv_root_masses = (1 / np.sqrt(masses)).repeat(3)
186 D = np.outer(inv_root_masses, inv_root_masses)
187 D *= force_constants
188 # find frequnecies and eigenvectors
189 w2_s, e_sai = np.linalg.eigh(D)
190 # reshape to get atom index and Cartesian index separate
191 e_sai = e_sai.T.reshape(-1, n_atoms, 3)
193 # The three modes closest to zero are assumed to be zero, ie acoustic sum rules are assumed
194 frequency_tol = 1e-6
195 argsort = np.argsort(np.abs(w2_s))
196 w2_gamma = w2_s[argsort][:3]
197 if np.any(np.abs(w2_gamma) > frequency_tol): 197 ↛ 198line 197 didn't jump to line 198 because the condition on line 197 was never true
198 logger.warning(f'Acoustic sum rules not enforced, squared frequencies: {w2_gamma}')
200 # warning if any imaginary modes
201 if np.any(w2_s < -frequency_tol): 201 ↛ 202line 201 didn't jump to line 202 because the condition on line 201 was never true
202 logger.warning('Imaginary modes present')
204 # treat imaginary modes as real
205 imag_mask = w2_s < -frequency_tol
206 w2_s[imag_mask] = imag_freq_factor * np.abs(w2_s[imag_mask])
208 self.w2_s = w2_s
209 self.e_sai = e_sai
210 self.masses = masses
212 def __call__(self, atoms, T, QM_statistics):
213 """ rattle atoms by adding displacements
215 Parameters
216 ----------
217 atoms : ase.Atoms
218 Ideal structure to add displacements to.
219 T : float
220 temperature in Kelvin
221 """
222 u_ai = _phonon_rattle(self.masses, T, self.w2_s, self.e_sai,
223 QM_statistics)
224 atoms.positions += u_ai