import sys import numpy as np import ase.units as u from ase.parallel import world, parprint, paropen from ase.phonons import Phonons from ase.vibrations import Vibrations from ase.utils.timing import Timer from ase.utils import convert_string_to_fd from ase.dft import monkhorst_pack class RamanCalculatorBase: def __init__(self, atoms, # XXX do we need atoms at this stage ? *args, name='raman', exext='.alpha', txt='-', verbose=False, comm=world, **kwargs): """ Parameters ---------- atoms: ase Atoms object exext: string Extension for excitation filenames txt: Output stream verbose: Verbosity level of output comm: Communicator, default world """ kwargs['name'] = name self.exname = kwargs.pop('exname', name) super().__init__(atoms, *args, **kwargs) self.exext = exext self.timer = Timer() self.txt = convert_string_to_fd(txt) self.verbose = verbose self.comm = comm def log(self, message, pre='# ', end='\n'): if self.verbose: self.txt.write(pre + message + end) self.txt.flush() class StaticRamanCalculatorBase(RamanCalculatorBase): """Base class for Raman intensities derived from static polarizabilities""" def __init__(self, atoms, exobj, exkwargs=None, *args, **kwargs): self.exobj = exobj if exkwargs is None: exkwargs = {} self.exkwargs = exkwargs super().__init__(atoms, *args, **kwargs) def calculate(self, atoms, filename, fd): # write forces super().calculate(atoms, filename, fd) # write static polarizability fname = filename.replace('.pckl', self.exext) np.savetxt(fname, self.exobj(**self.exkwargs).calculate(atoms)) class StaticRamanCalculator(StaticRamanCalculatorBase, Vibrations): pass class StaticRamanPhononsCalculator(StaticRamanCalculatorBase, Phonons): pass class RamanBase: def __init__(self, atoms, # XXX do we need atoms at this stage ? *args, name='raman', exname=None, exext='.alpha', txt='-', verbose=False, comm=world, **kwargs): """ Parameters ---------- atoms: ase Atoms object exext: string Extension for excitation filenames txt: Output stream verbose: Verbosity level of output comm: Communicator, default world """ self.atoms = atoms self.name = name if exname is None: self.exname = name else: self.exname = exname self.exext = exext self.timer = Timer() self.txt = convert_string_to_fd(txt) self.verbose = verbose self.comm = comm def log(self, message, pre='# ', end='\n'): if self.verbose: self.txt.write(pre + message + end) self.txt.flush() class RamanData(RamanBase): """Base class to evaluate Raman spectra from pre-computed data""" def __init__(self, atoms, # XXX do we need atoms at this stage ? *args, exname=None, # name for excited state calculations **kwargs): """ Parameters ---------- atoms: ase Atoms object exname: string name for excited state calculations (defaults to name), used for reading excitations """ super().__init__(atoms, *args, **kwargs) if exname is None: exname = kwargs.get('name', self.name) self.exname = exname self._already_read = False def get_energies(self): self.calculate_energies_and_modes() return self.om_Q def init_parallel_read(self): """Initialize variables for parallel read""" rank = self.comm.rank indices = self.indices self.ndof = 3 * len(indices) myn = -(-self.ndof // self.comm.size) # ceil divide self.slize = s = slice(myn * rank, myn * (rank + 1)) self.myindices = np.repeat(indices, 3)[s] self.myxyz = ('xyz' * len(indices))[s] self.myr = range(self.ndof)[s] self.mynd = len(self.myr) def read(self, *args, **kwargs): """Read data from a pre-performed calculation.""" if self._already_read: return self.timer.start('read') self.timer.start('vibrations') self.vibrations.read(*args, **kwargs) self.timer.stop('vibrations') self.timer.start('excitations') self.init_parallel_read() self.read_excitations() self.timer.stop('excitations') self._already_read = True self.timer.stop('read') @staticmethod def m2(z): return (z * z.conj()).real def map_to_modes(self, V_rcc): self.timer.start('map R2Q') V_qcc = (V_rcc.T * self.im_r).T # units Angstrom^2 / sqrt(amu) V_Qcc = np.dot(V_qcc.T, self.modes_Qq.T).T self.timer.stop('map R2Q') return V_Qcc def me_Qcc(self, *args, **kwargs): """Full matrix element Returns ------- Matrix element in e^2 Angstrom^2 / eV """ # Angstrom^2 / sqrt(amu) elme_Qcc = self.electronic_me_Qcc(*args, **kwargs) # Angstrom^3 -> e^2 Angstrom^2 / eV elme_Qcc /= u.Hartree * u.Bohr # e^2 Angstrom / eV / sqrt(amu) return elme_Qcc * self.vib01_Q[:, None, None] def get_absolute_intensities(self, delta=0, **kwargs): """Absolute Raman intensity or Raman scattering factor Parameter --------- delta: float pre-factor for asymmetric anisotropy, default 0 References ---------- Porezag and Pederson, PRB 54 (1996) 7830-7836 (delta=0) Baiardi and Barone, JCTC 11 (2015) 3267-3280 (delta=5) Returns ------- raman intensity, unit Ang**4/amu """ alpha2_r, gamma2_r, delta2_r = self._invariants( self.electronic_me_Qcc(**kwargs)) return 45 * alpha2_r + delta * delta2_r + 7 * gamma2_r def intensity(self, *args, **kwargs): """Raman intensity Returns ------- unit e^4 Angstrom^4 / eV^2 """ self.calculate_energies_and_modes() m2 = Raman.m2 alpha_Qcc = self.me_Qcc(*args, **kwargs) if not self.observation: # XXXX remove """Simple sum, maybe too simple""" return m2(alpha_Qcc).sum(axis=1).sum(axis=1) # XXX enable when appropriate # if self.observation['orientation'].lower() != 'random': # raise NotImplementedError('not yet') # random orientation of the molecular frame # Woodward & Long, # Guthmuller, J. J. Chem. Phys. 2016, 144 (6), 64106 alpha2_r, gamma2_r, delta2_r = self._invariants(alpha_Qcc) if self.observation['geometry'] == '-Z(XX)Z': # Porto's notation return (45 * alpha2_r + 5 * delta2_r + 4 * gamma2_r) / 45. elif self.observation['geometry'] == '-Z(XY)Z': # Porto's notation return gamma2_r / 15. elif self.observation['scattered'] == 'Z': # scattered light in direction of incoming light return (45 * alpha2_r + 5 * delta2_r + 7 * gamma2_r) / 45. elif self.observation['scattered'] == 'parallel': # scattered light perendicular and # polarization in plane return 6 * gamma2_r / 45. elif self.observation['scattered'] == 'perpendicular': # scattered light perendicular and # polarization out of plane return (45 * alpha2_r + 5 * delta2_r + 7 * gamma2_r) / 45. else: raise NotImplementedError def _invariants(self, alpha_Qcc): """Raman invariants Parameter --------- alpha_Qcc: array Matrix element or polarizability tensor Reference --------- Derek A. Long, The Raman Effect, ISBN 0-471-49028-8 Returns ------- mean polarizability, anisotropy, asymmetric anisotropy """ m2 = Raman.m2 alpha2_r = m2(alpha_Qcc[:, 0, 0] + alpha_Qcc[:, 1, 1] + alpha_Qcc[:, 2, 2]) / 9. delta2_r = 3 / 4. * ( m2(alpha_Qcc[:, 0, 1] - alpha_Qcc[:, 1, 0]) + m2(alpha_Qcc[:, 0, 2] - alpha_Qcc[:, 2, 0]) + m2(alpha_Qcc[:, 1, 2] - alpha_Qcc[:, 2, 1])) gamma2_r = (3 / 4. * (m2(alpha_Qcc[:, 0, 1] + alpha_Qcc[:, 1, 0]) + m2(alpha_Qcc[:, 0, 2] + alpha_Qcc[:, 2, 0]) + m2(alpha_Qcc[:, 1, 2] + alpha_Qcc[:, 2, 1])) + (m2(alpha_Qcc[:, 0, 0] - alpha_Qcc[:, 1, 1]) + m2(alpha_Qcc[:, 0, 0] - alpha_Qcc[:, 2, 2]) + m2(alpha_Qcc[:, 1, 1] - alpha_Qcc[:, 2, 2])) / 2) return alpha2_r, gamma2_r, delta2_r def summary(self, log=sys.stdout): """Print summary for given omega [eV]""" hnu = self.get_energies() intensities = self.get_absolute_intensities() te = int(np.log10(intensities.max())) - 2 scale = 10**(-te) if not te: ts = '' elif te > -2 and te < 3: ts = str(10**te) else: ts = '10^{0}'.format(te) if isinstance(log, str): log = paropen(log, 'a') parprint('-------------------------------------', file=log) parprint(' Mode Frequency Intensity', file=log) parprint(' # meV cm^-1 [{0}A^4/amu]'.format(ts), file=log) parprint('-------------------------------------', file=log) for n, e in enumerate(hnu): if e.imag != 0: c = 'i' e = e.imag else: c = ' ' e = e.real parprint('%3d %6.1f%s %7.1f%s %9.2f' % (n, 1000 * e, c, e / u.invcm, c, intensities[n] * scale), file=log) parprint('-------------------------------------', file=log) # XXX enable this in phonons # parprint('Zero-point energy: %.3f eV' % # self.vibrations.get_zero_point_energy(), # file=log) class Raman(RamanData): def __init__(self, atoms, *args, **kwargs): super().__init__(atoms, *args, **kwargs) for key in ['txt', 'exext', 'exname']: kwargs.pop(key, None) kwargs['name'] = kwargs.get('name', self.name) self.vibrations = Vibrations(atoms, *args, **kwargs) self.delta = self.vibrations.delta self.indices = self.vibrations.indices def calculate_energies_and_modes(self): if hasattr(self, 'im_r'): return self.read() self.timer.start('energies_and_modes') self.im_r = self.vibrations.im self.modes_Qq = self.vibrations.modes self.om_Q = self.vibrations.hnu.real # energies in eV self.H = self.vibrations.H # XXX used in albrecht.py # pre-factors for one vibrational excitation with np.errstate(divide='ignore'): self.vib01_Q = np.where(self.om_Q > 0, 1. / np.sqrt(2 * self.om_Q), 0) # -> sqrt(amu) * Angstrom self.vib01_Q *= np.sqrt(u.Ha * u._me / u._amu) * u.Bohr self.timer.stop('energies_and_modes') class RamanPhonons(RamanData): def __init__(self, atoms, *args, **kwargs): RamanData.__init__(self, atoms, *args, **kwargs) for key in ['txt', 'exext', 'exname']: kwargs.pop(key, None) kwargs['name'] = kwargs.get('name', self.name) self.vibrations = Phonons(atoms, *args, **kwargs) self.delta = self.vibrations.delta self.indices = self.vibrations.indices self.kpts = (1, 1, 1) @property def kpts(self): return self._kpts @kpts.setter def kpts(self, kpts): if not hasattr(self, '_kpts') or kpts != self._kpts: self._kpts = kpts self.kpts_kc = monkhorst_pack(self.kpts) if hasattr(self, 'im_r'): del self.im_r # we'll have to recalculate def calculate_energies_and_modes(self): if not self._already_read: if hasattr(self, 'im_r'): del self.im_r self.read() if not hasattr(self, 'im_r'): self.timer.start('band_structure') omega_kl, u_kl = self.vibrations.band_structure( self.kpts_kc, modes=True, verbose=self.verbose) self.im_r = self.vibrations.m_inv_x self.om_Q = omega_kl.ravel().real # energies in eV self.modes_Qq = u_kl.reshape(len(self.om_Q), 3 * len(self.atoms)) self.modes_Qq /= self.im_r self.om_v = self.om_Q # pre-factors for one vibrational excitation with np.errstate(divide='ignore', invalid='ignore'): self.vib01_Q = np.where( self.om_Q > 0, 1. / np.sqrt(2 * self.om_Q), 0) # -> sqrt(amu) * Angstrom self.vib01_Q *= np.sqrt(u.Ha * u._me / u._amu) * u.Bohr self.timer.stop('band_structure')