1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179
|
import numpy as np
import ase.units as u
from ase.vibrations.raman import Raman, RamanPhonons
from ase.vibrations.resonant_raman import ResonantRaman
from ase.calculators.excitation_list import polarizability
class Placzek(ResonantRaman):
"""Raman spectra within the Placzek approximation."""
def __init__(self, *args, **kwargs):
self._approx = 'PlaczekAlpha'
ResonantRaman.__init__(self, *args, **kwargs)
def set_approximation(self, value):
raise ValueError('Approximation can not be set.')
def read_excitations(self):
"""Read excitations from files written"""
self.ex0E_p = None # mark as read
self.exm_r = []
self.exp_r = []
for a, i in zip(self.myindices, self.myxyz):
exname = '%s.%d%s-' % (self.exname, a, i) + self.exext
self.log('reading ' + exname)
self.exm_r.append(self.exobj.read(exname, **self.exkwargs))
exname = '%s.%d%s+' % (self.exname, a, i) + self.exext
self.log('reading ' + exname)
self.exp_r.append(self.exobj.read(exname, **self.exkwargs))
def electronic_me_Qcc(self, omega, gamma=0):
self.calculate_energies_and_modes()
self.timer.start('init')
V_rcc = np.zeros((self.ndof, 3, 3), dtype=complex)
pre = 1. / (2 * self.delta)
pre *= u.Hartree * u.Bohr # e^2Angstrom^2 / eV -> Angstrom^3
om = omega
if gamma:
om += 1j * gamma
self.timer.stop('init')
self.timer.start('alpha derivatives')
for i, r in enumerate(self.myr):
V_rcc[r] = pre * (
polarizability(self.exp_r[i], om,
form=self.dipole_form, tensor=True) -
polarizability(self.exm_r[i], om,
form=self.dipole_form, tensor=True))
self.comm.sum(V_rcc)
self.timer.stop('alpha derivatives')
return self.map_to_modes(V_rcc)
class PlaczekStatic(Raman):
def read_excitations(self):
"""Read excitations from files written"""
self.al0_rr = None # mark as read
self.alm_rr = []
self.alp_rr = []
for a, i in zip(self.myindices, self.myxyz):
exname = '%s.%d%s-' % (self.exname, a, i) + self.exext
self.log('reading ' + exname)
self.alm_rr.append(np.loadtxt(exname))
exname = '%s.%d%s+' % (self.exname, a, i) + self.exext
self.log('reading ' + exname)
self.alp_rr.append(np.loadtxt(exname))
def electronic_me_Qcc(self):
self.calculate_energies_and_modes()
self.timer.start('init')
V_rcc = np.zeros((self.ndof, 3, 3), dtype=complex)
pre = 1. / (2 * self.delta)
pre *= u.Hartree * u.Bohr # e^2Angstrom^2 / eV -> Angstrom^3
self.timer.stop('init')
self.timer.start('alpha derivatives')
for i, r in enumerate(self.myr):
V_rcc[r] = pre * (self.alp_rr[i] - self.alm_rr[i])
self.comm.sum(V_rcc)
self.timer.stop('alpha derivatives')
return self.map_to_modes(V_rcc)
class PlaczekStaticPhonons(RamanPhonons, PlaczekStatic):
pass
class Profeta(ResonantRaman):
"""Profeta type approximations.
Reference
---------
Mickael Profeta and Francesco Mauri
Phys. Rev. B 63 (2000) 245415
"""
def __init__(self, *args, **kwargs):
self.set_approximation(kwargs.pop('approximation', 'Profeta'))
self.nonresonant = kwargs.pop('nonresonant', True)
ResonantRaman.__init__(self, *args, **kwargs)
def set_approximation(self, value):
approx = value.lower()
if approx in ['profeta', 'placzek', 'p-p']:
self._approx = value
else:
raise ValueError('Please use "Profeta", "Placzek" or "P-P".')
def electronic_me_profeta_rcc(self, omega, gamma=0.1,
energy_derivative=False):
"""Raman spectra in Profeta and Mauri approximation
Returns
-------
Electronic matrix element, unit Angstrom^2
"""
self.calculate_energies_and_modes()
self.timer.start('amplitudes')
self.timer.start('init')
V_rcc = np.zeros((self.ndof, 3, 3), dtype=complex)
pre = 1. / (2 * self.delta)
pre *= u.Hartree * u.Bohr # e^2Angstrom^2 / eV -> Angstrom^3
self.timer.stop('init')
def kappa_cc(me_pc, e_p, omega, gamma, form='v'):
"""Kappa tensor after Profeta and Mauri
PRB 63 (2001) 245415"""
k_cc = np.zeros((3, 3), dtype=complex)
for p, me_c in enumerate(me_pc):
me_cc = np.outer(me_c, me_c.conj())
k_cc += me_cc / (e_p[p] - omega - 1j * gamma)
if self.nonresonant:
k_cc += me_cc.conj() / (e_p[p] + omega + 1j * gamma)
return k_cc
self.timer.start('kappa')
mr = 0
for a, i, r in zip(self.myindices, self.myxyz, self.myr):
if not energy_derivative < 0:
V_rcc[r] += pre * (
kappa_cc(self.expm_rpc[mr], self.ex0E_p,
omega, gamma, self.dipole_form) -
kappa_cc(self.exmm_rpc[mr], self.ex0E_p,
omega, gamma, self.dipole_form))
if energy_derivative:
V_rcc[r] += pre * (
kappa_cc(self.ex0m_pc, self.expE_rp[mr],
omega, gamma, self.dipole_form) -
kappa_cc(self.ex0m_pc, self.exmE_rp[mr],
omega, gamma, self.dipole_form))
mr += 1
self.comm.sum(V_rcc)
self.timer.stop('kappa')
self.timer.stop('amplitudes')
return V_rcc
def electronic_me_Qcc(self, omega, gamma):
self.read()
Vel_rcc = np.zeros((self.ndof, 3, 3), dtype=complex)
approximation = self.approximation.lower()
if approximation == 'profeta':
Vel_rcc += self.electronic_me_profeta_rcc(omega, gamma)
elif approximation == 'placzek':
Vel_rcc += self.electronic_me_profeta_rcc(omega, gamma, True)
elif approximation == 'p-p':
Vel_rcc += self.electronic_me_profeta_rcc(omega, gamma, -1)
else:
raise RuntimeError(
'Bug: call with {0} should not happen!'.format(
self.approximation))
return self.map_to_modes(Vel_rcc)
|