"""
Phonon unfolding: Reciprocal space method. The method is described in
P. B. Allen et al. Phys Rev B 87, 085322 (2013).
This method should be also applicable to other bloch waves on discrete grid, eg. electrons wave function in wannier basis set, magnons, etc. Now only phonon istested.
"""
from ase.build import make_supercell
from ase.atoms import Atoms
import numpy as np
[docs]class Unfolder:
""" phonon unfolding class"""
def __init__(
self,
cell,
basis,
positions,
supercell_matrix,
eigenvectors,
qpoints,
tol_r=0.1,
compare=None,
phase=True,
):
"""
Params:
===================
cell: cell matrix. [a,b,c]
basis: name of the basis. It's used to decide if two basis can
be identical by translation. eg. for phonons, the basis can be
['x','y','z']*natoms, for electrons, it can be
['Ni|dxy','Mn|dxy'] if the two dxy are seen as different, or
['dxy','dxy'] if they are seen as the same.
positions: positions(->basis).
supercell matrix: The matrix that convert the primitive cell
to supercell.
eigenvectors: The phonon eigenvectors. format np.array()
index=[ikpts, ifreq, 3*iatoms+j]. j=0..2
qpoints: list of q-points.
tol_r: tolerance. If abs(a-b) <r, they are seen as the same atom.
"""
self._cell = cell
self._basis = basis
self._positions = positions
self._scmat = supercell_matrix
self._evecs = eigenvectors
self._qpts = qpoints
self._tol_r = tol_r
self._trans_rs = None
self._trans_indices = None
self._make_translate_maps()
self._phase = phase
return
def _translate(self, evec, r):
"""
T(r) psi: r is integer numbers of primitive cell lattice matrix.
Params:
=================
evec: an eigen vector of supercell
r: The translate vector
Returns:
================
tevec: translated vector.
"""
pass
def _make_translate_maps(self):
"""
find the mapping between supercell and translated cell.
Returns:
===============
A N * nbasis array.
index[i] is the mapping from supercell to translated supercell so that
T(r_i) psi = psi[indices[i]].
TODO: vacancies/add_atoms not supported. How to do it? For
vacancies, a ghost atom can be added. For add_atom, maybe we
can just ignore them? Will it change the energy spectrum?
"""
a1 = Atoms(symbols="H", positions=[(0, 0, 0)], cell=[1, 1, 1])
sc = make_supercell(a1, self._scmat)
rs = sc.get_scaled_positions()
positions = self._positions
indices = np.zeros([len(rs), len(positions)], dtype="int32")
for i, ri in enumerate(rs):
inds = []
Tpositions = positions + np.array(ri)
close_to_int = lambda x: np.all(
np.abs(x - np.round(x)) < self._tol_r)
for i_basis, pos in enumerate(positions):
for j_basis, Tpos in enumerate(Tpositions):
dpos = Tpos - pos
if close_to_int(dpos) and (
self._basis[i_basis] == self._basis[j_basis]):
# indices[i, j_atom * self._ndim:j_atom * self._ndim + self._ndim] = np.arange(i_atom * self._ndim, i_atom * self._ndim + self._ndim)
indices[i, j_basis] = i_basis
self._trans_rs = rs
self._trans_indices = indices
# print(indices)
[docs] def get_weight(self, evec, qpt, G=None):
"""
get the weight of a mode which has the wave vector of qpt and
eigenvector of evec.
W= sum_1^N < evec| T(r_i)exp(-I (K+G) * r_i| evec>, here
G=0. T(r_i)exp(-I K r_i)| evec> = evec[indices[i]]
"""
if G is None:
G = np.zeros_like(qpt)
weight = 0j
N = len(self._trans_rs)
for r_i, ind in zip(self._trans_rs, self._trans_indices):
if self._phase:
weight += (np.vdot(evec, evec[ind]) *
np.exp(1j * 2 * np.pi * np.dot(qpt + G, r_i)) / N)
else:
weight += (np.vdot(evec, evec[ind]) / N *
np.exp(-1j * 2 * np.pi * np.dot(G, r_i)))
return weight.real
[docs] def get_weights(self):
"""
Get the weight for all the modes.
"""
nqpts, nfreqs = self._evecs.shape[0], self._evecs.shape[1]
weights = np.zeros([nqpts, nfreqs])
for iqpt in range(nqpts):
for ifreq in range(nfreqs):
weights[iqpt, ifreq] = self.get_weight(
self._evecs[iqpt, ifreq, :], self._qpts[iqpt])
self._weights = weights
return self._weights