Source code for hiphive.structure_generation.phonon

import numpy as np
import ase.units as aunits
from hiphive.input_output.logging_tools import logger

logger = logger.getChild(__name__)


[docs]def generate_phonon_rattled_structures(atoms, fc2, n_structures, temperature, QM_statistics=False, imag_freq_factor=1.0): """ 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 : ase.Atoms prototype structure fc2 : numpy.ndarray second order force constant matrix, with shape `(3N, 3N)` or `(N, N, 3, 3)`. The conversion will be done internally if. n_structures : int number of structures to generate temperature : float temperature in Kelvin QM_statistics : bool if the amplitude of the quantum harmonic oscillator shoule be used instead of the classical amplitude imag_freq_factor: float If a squared frequency, w2, is negative then it is set to w2 = imag_freq_factor * np.abs(w2) Returns ------- list(ase.Atoms) 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