1import numpy as np 

2 

3import ase.units as u 

4from ase.vibrations.raman import Raman, RamanPhonons 

5from ase.vibrations.resonant_raman import ResonantRaman 

6from ase.calculators.excitation_list import polarizability 

7 

8 

9class Placzek(ResonantRaman): 

10 """Raman spectra within the Placzek approximation.""" 

11 def __init__(self, *args, **kwargs): 

12 self._approx = 'PlaczekAlpha' 

13 ResonantRaman.__init__(self, *args, **kwargs) 

14 

15 def set_approximation(self, value): 

16 raise ValueError('Approximation can not be set.') 

17 

18 def _signed_disps(self, sign): 

19 for a, i in zip(self.myindices, self.myxyz): 

20 yield self._disp(a, i, sign) 

21 

22 def _read_exobjs(self, sign): 

23 return [disp.read_exobj() for disp in self._signed_disps(sign)] 

24 

25 def read_excitations(self): 

26 """Read excitations from files written""" 

27 self.ex0E_p = None # mark as read 

28 self.exm_r = self._read_exobjs(sign=-1) 

29 self.exp_r = self._read_exobjs(sign=1) 

30 

31 def electronic_me_Qcc(self, omega, gamma=0): 

32 self.calculate_energies_and_modes() 

33 

34 V_rcc = np.zeros((self.ndof, 3, 3), dtype=complex) 

35 pre = 1. / (2 * self.delta) 

36 pre *= u.Hartree * u.Bohr # e^2Angstrom^2 / eV -> Angstrom^3 

37 

38 om = omega 

39 if gamma: 

40 om += 1j * gamma 

41 

42 for i, r in enumerate(self.myr): 

43 V_rcc[r] = pre * ( 

44 polarizability(self.exp_r[i], om, 

45 form=self.dipole_form, tensor=True) - 

46 polarizability(self.exm_r[i], om, 

47 form=self.dipole_form, tensor=True)) 

48 self.comm.sum(V_rcc) 

49 

50 return self.map_to_modes(V_rcc) 

51 

52 

53class PlaczekStatic(Raman): 

54 def read_excitations(self): 

55 """Read excitations from files written""" 

56 self.al0_rr = None # mark as read 

57 self.alm_rr = [] 

58 self.alp_rr = [] 

59 for a, i in zip(self.myindices, self.myxyz): 

60 for sign, al_rr in zip([-1, 1], [self.alm_rr, self.alp_rr]): 

61 disp = self._disp(a, i, sign) 

62 al_rr.append(disp.load_static_polarizability()) 

63 

64 def electronic_me_Qcc(self): 

65 self.calculate_energies_and_modes() 

66 

67 V_rcc = np.zeros((self.ndof, 3, 3), dtype=complex) 

68 pre = 1. / (2 * self.delta) 

69 pre *= u.Hartree * u.Bohr # e^2Angstrom^2 / eV -> Angstrom^3 

70 

71 for i, r in enumerate(self.myr): 

72 V_rcc[r] = pre * (self.alp_rr[i] - self.alm_rr[i]) 

73 self.comm.sum(V_rcc) 

74 

75 return self.map_to_modes(V_rcc) 

76 

77 

78class PlaczekStaticPhonons(RamanPhonons, PlaczekStatic): 

79 pass 

80 

81 

82class Profeta(ResonantRaman): 

83 """Profeta type approximations. 

84 

85 Reference 

86 --------- 

87 Mickael Profeta and Francesco Mauri 

88 Phys. Rev. B 63 (2000) 245415 

89 """ 

90 def __init__(self, *args, **kwargs): 

91 self.set_approximation(kwargs.pop('approximation', 'Profeta')) 

92 self.nonresonant = kwargs.pop('nonresonant', True) 

93 ResonantRaman.__init__(self, *args, **kwargs) 

94 

95 def set_approximation(self, value): 

96 approx = value.lower() 

97 if approx in ['profeta', 'placzek', 'p-p']: 

98 self._approx = value 

99 else: 

100 raise ValueError('Please use "Profeta", "Placzek" or "P-P".') 

101 

102 def electronic_me_profeta_rcc(self, omega, gamma=0.1, 

103 energy_derivative=False): 

104 """Raman spectra in Profeta and Mauri approximation 

105 

106 Returns 

107 ------- 

108 Electronic matrix element, unit Angstrom^2 

109 """ 

110 self.calculate_energies_and_modes() 

111 

112 V_rcc = np.zeros((self.ndof, 3, 3), dtype=complex) 

113 pre = 1. / (2 * self.delta) 

114 pre *= u.Hartree * u.Bohr # e^2Angstrom^2 / eV -> Angstrom^3 

115 

116 def kappa_cc(me_pc, e_p, omega, gamma, form='v'): 

117 """Kappa tensor after Profeta and Mauri 

118 PRB 63 (2001) 245415""" 

119 k_cc = np.zeros((3, 3), dtype=complex) 

120 for p, me_c in enumerate(me_pc): 

121 me_cc = np.outer(me_c, me_c.conj()) 

122 k_cc += me_cc / (e_p[p] - omega - 1j * gamma) 

123 if self.nonresonant: 

124 k_cc += me_cc.conj() / (e_p[p] + omega + 1j * gamma) 

125 return k_cc 

126 

127 mr = 0 

128 for a, i, r in zip(self.myindices, self.myxyz, self.myr): 

129 if not energy_derivative < 0: 

130 V_rcc[r] += pre * ( 

131 kappa_cc(self.expm_rpc[mr], self.ex0E_p, 

132 omega, gamma, self.dipole_form) - 

133 kappa_cc(self.exmm_rpc[mr], self.ex0E_p, 

134 omega, gamma, self.dipole_form)) 

135 if energy_derivative: 

136 V_rcc[r] += pre * ( 

137 kappa_cc(self.ex0m_pc, self.expE_rp[mr], 

138 omega, gamma, self.dipole_form) - 

139 kappa_cc(self.ex0m_pc, self.exmE_rp[mr], 

140 omega, gamma, self.dipole_form)) 

141 mr += 1 

142 self.comm.sum(V_rcc) 

143 

144 return V_rcc 

145 

146 def electronic_me_Qcc(self, omega, gamma): 

147 self.read() 

148 Vel_rcc = np.zeros((self.ndof, 3, 3), dtype=complex) 

149 approximation = self.approximation.lower() 

150 if approximation == 'profeta': 

151 Vel_rcc += self.electronic_me_profeta_rcc(omega, gamma) 

152 elif approximation == 'placzek': 

153 Vel_rcc += self.electronic_me_profeta_rcc(omega, gamma, True) 

154 elif approximation == 'p-p': 

155 Vel_rcc += self.electronic_me_profeta_rcc(omega, gamma, -1) 

156 else: 

157 raise RuntimeError( 

158 'Bug: call with {0} should not happen!'.format( 

159 self.approximation)) 

160 

161 return self.map_to_modes(Vel_rcc)