import numpy as np
from ase import Atoms
import ase.units as aunits
from hiphive.input_output.logging_tools import logger
logger = logger.getChild(__name__)
[docs]
def generate_phonon_rattled_structures(
atoms: Atoms,
fc2: np.ndarray,
n_structures: int,
temperature: float,
QM_statistics: bool = False,
imag_freq_factor: float = 1.0,
) -> list[Atoms]:
r"""
Returns list of phonon-rattled configurations.
Configurations are generated by superimposing harmonic phonon
eigenmodes with random amplitudes and phase factors consistent
with a certain temperature.
Let :math:`\boldsymbol{X}_{ai}` be the phonon modes indexed by atom
:math:`a` and mode :math:`i`, :math:`\omega_i` the phonon
frequencies, and let :math:`0 < Q_i \leq 1` and :math:`0 \leq
U_i < 1` be uniformly random numbers. Then
.. math::
\boldsymbol{R}_a
= \boldsymbol{R}^0_a
+ \left<\frac{k_B T}{m_a} \right>^{1/2}
\sum_i \frac{1}{\omega_i} \boldsymbol{X}_{ai}
\sqrt{-2 \ln Q_i} \cos(\pi \omega_i U_i)
See: West and Estreicher, Physical Review Letters **96**, 115504 (2006)
or
.. math::
\boldsymbol{R}_a
= \boldsymbol{R}^0_a
+ \sum_i
\left(\frac{\hbar(0.5 + n(\omega_i, T))}{m_a\omega_i} \right)^{1/2}
\boldsymbol{X}_{ai}
\sqrt{-2 \ln Q_i} \cos(\pi \omega_i U_i)
where
.. math::
n = \frac{1}{e^{\hbar\omega_i/k_BT}-1}
Parameters
----------
atoms
Prototype structure.
fc2
Second order force constant matrix, with shape `(3N, 3N)` or
`(N, N, 3, 3)`. The conversion will be done internally if.
n_structures
Number of structures to generate.
temperature
Temperature in Kelvin.
QM_statistics
If the amplitude of the quantum harmonic oscillator shoule be used
instead of the classical amplitude.
imag_freq_factor
If a squared frequency, `w2`, is negative it is set to
`w2 = imag_freq_factor * np.abs(w2)`.
Returns
-------
Generated structures.
"""
structures = []
pr = _PhononRattler(atoms.get_masses(), fc2, imag_freq_factor)
for _ in range(n_structures):
atoms_tmp = atoms.copy()
pr(atoms_tmp, temperature, QM_statistics)
structures.append(atoms_tmp)
return structures
def _n_BE(T, w_s):
"""
Bose-Einstein distribution function.
Parameters
---------
T : float
Temperature in Kelvin
w_s: numpy.ndarray
frequencies in eV (3*N,)
Returns
------
Bose-Einstein distribution for each energy at a given temperature
"""
with np.errstate(divide='raise', over='raise'):
try:
n = 1 / (np.exp(w_s / (aunits.kB * T)) - 1)
except Exception:
n = np.zeros_like(w_s)
return n
def _phonon_rattle(m_a, T, w2_s, e_sai, QM_statistics):
""" Thermal excitation of phonon modes as described by West and
Estreicher, Physical Review Letters **96**, 115504 (2006).
_s is a mode index
_i is a Carteesian index
_a is an atom index
Parameters
----------
m_a : numpy.ndarray
masses (N,)
T : float
temperature in Kelvin
w2_s : numpy.ndarray
the squared frequencies from the eigenvalue problem (3*N,)
e_sai : numpy.ndarray
polarizations (3*N, N, 3)
QM_statistics : bool
if the amplitude of the quantum harmonic oscillator shoule be used
instead of the classical amplitude
Returns
-------
displacements : numpy.ndarray
shape (N, 3)
"""
n_modes = 3 * len(m_a)
# skip 3 translational modes
argsort = np.argsort(np.abs(w2_s))
e_sai = e_sai[argsort][3:]
w2_s = w2_s[argsort][3:]
w_s = np.sqrt(np.abs(w2_s))
prefactor_a = np.sqrt(1 / m_a).reshape(-1, 1)
if QM_statistics:
hbar = aunits._hbar * aunits.J * aunits.s
frequencyfactor_s = np.sqrt(hbar * (0.5 + _n_BE(T, hbar * w_s)) / w_s)
else:
frequencyfactor_s = 1 / w_s
prefactor_a *= np.sqrt(aunits.kB * T)
phases_s = np.random.uniform(0, 2 * np.pi, size=n_modes - 3)
amplitudes_s = np.sqrt(-2 * np.log(1 - np.random.random(n_modes - 3)))
u_ai = prefactor_a * np.tensordot(
amplitudes_s * np.cos(phases_s) * frequencyfactor_s, e_sai, (0, 0))
return u_ai # displacements
class _PhononRattler:
"""
Class to be able to conveniently save modes and frequencies needed
for phonon rattle.
Parameters
----------
masses : numpy.ndarray
masses (N,)
force_constants : numpy.ndarray
second order force constant matrix, with shape `(3N, 3N)` or
`(N, N, 3, 3)`. The conversion will be done internally if.
imag_freq_factor: float
If a squared frequency, w2, is negative then it is set to
w2 = imag_freq_factor * np.abs(w2)
"""
def __init__(self, masses, force_constants, imag_freq_factor=1.0):
n_atoms = len(masses)
if len(force_constants.shape) == 4: # assume shape = (n_atoms, n_atoms, 3, 3)
force_constants = force_constants.transpose(0, 2, 1, 3)
force_constants = force_constants.reshape(3 * n_atoms, 3 * n_atoms)
# Now the fc should have shape = (n_atoms * 3, n_atoms * 3)
# Construct the dynamical matrix
inv_root_masses = (1 / np.sqrt(masses)).repeat(3)
D = np.outer(inv_root_masses, inv_root_masses)
D *= force_constants
# find frequnecies and eigenvectors
w2_s, e_sai = np.linalg.eigh(D)
# reshape to get atom index and Cartesian index separate
e_sai = e_sai.T.reshape(-1, n_atoms, 3)
# The three modes closest to zero are assumed to be zero, ie acoustic sum rules are assumed
frequency_tol = 1e-6
argsort = np.argsort(np.abs(w2_s))
w2_gamma = w2_s[argsort][:3]
if np.any(np.abs(w2_gamma) > frequency_tol):
logger.warning(f'Acoustic sum rules not enforced, squared frequencies: {w2_gamma}')
# warning if any imaginary modes
if np.any(w2_s < -frequency_tol):
logger.warning('Imaginary modes present')
# treat imaginary modes as real
imag_mask = w2_s < -frequency_tol
w2_s[imag_mask] = imag_freq_factor * np.abs(w2_s[imag_mask])
self.w2_s = w2_s
self.e_sai = e_sai
self.masses = masses
def __call__(self, atoms, T, QM_statistics):
""" rattle atoms by adding displacements
Parameters
----------
atoms : ase.Atoms
Ideal structure to add displacements to.
T : float
temperature in Kelvin
"""
u_ai = _phonon_rattle(self.masses, T, self.w2_s, self.e_sai,
QM_statistics)
atoms.positions += u_ai
