Coverage for /builds/kinetik161/ase/ase/vibrations/placzek.py: 96.88%

1import numpy as np 

2 

3import ase.units as u 

4from ase.calculators.excitation_list import polarizability 

5from ase.vibrations.raman import Raman, RamanPhonons 

6from ase.vibrations.resonant_raman import ResonantRaman 

7 

8 

9class Placzek(ResonantRaman): 

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

11 

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

13 self._approx = 'PlaczekAlpha' 

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

15 

16 def set_approximation(self, value): 

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

18 

19 def _signed_disps(self, sign): 

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

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

22 

23 def _read_exobjs(self, sign): 

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

25 

26 def read_excitations(self): 

27 """Read excitations from files written""" 

28 self.ex0E_p = None # mark as read 

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

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

31 

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

33 self.calculate_energies_and_modes() 

34 

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

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

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

38 

39 om = omega 

40 if gamma: 

41 om += 1j * gamma 

42 

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

44 V_rcc[r] = pre * ( 

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

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

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

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

49 self.comm.sum(V_rcc) 

50 

51 return self.map_to_modes(V_rcc) 

52 

53 

54class PlaczekStatic(Raman): 

55 def read_excitations(self): 

56 """Read excitations from files written""" 

57 self.al0_rr = None # mark as read 

58 self.alm_rr = [] 

59 self.alp_rr = [] 

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

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

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

63 al_rr.append(disp.load_static_polarizability()) 

64 

65 def electronic_me_Qcc(self): 

66 self.calculate_energies_and_modes() 

67 

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

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

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

71 

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

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

74 self.comm.sum(V_rcc) 

75 

76 return self.map_to_modes(V_rcc) 

77 

78 

79class PlaczekStaticPhonons(RamanPhonons, PlaczekStatic): 

80 pass 

81 

82 

83class Profeta(ResonantRaman): 

84 """Profeta type approximations. 

85 

86 Reference 

87 --------- 

88 Mickael Profeta and Francesco Mauri 

89 Phys. Rev. B 63 (2000) 245415 

90 """ 

91 

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

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

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

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

96 

97 def set_approximation(self, value): 

98 approx = value.lower() 

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

100 self._approx = value 

101 else: 

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

103 

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

105 energy_derivative=False): 

106 """Raman spectra in Profeta and Mauri approximation 

107 

108 Returns 

109 ------- 

110 Electronic matrix element, unit Angstrom^2 

111 """ 

112 self.calculate_energies_and_modes() 

113 

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

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

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

117 

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

119 """Kappa tensor after Profeta and Mauri 

120 PRB 63 (2001) 245415""" 

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

122 for p, me_c in enumerate(me_pc): 

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

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

125 if self.nonresonant: 

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

127 return k_cc 

128 

129 mr = 0 

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

131 if not energy_derivative < 0: 

132 V_rcc[r] += pre * ( 

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

134 omega, gamma, self.dipole_form) - 

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

136 omega, gamma, self.dipole_form)) 

137 if energy_derivative: 

138 V_rcc[r] += pre * ( 

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

140 omega, gamma, self.dipole_form) - 

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

142 omega, gamma, self.dipole_form)) 

143 mr += 1 

144 self.comm.sum(V_rcc) 

145 

146 return V_rcc 

147 

148 def electronic_me_Qcc(self, omega, gamma): 

149 self.read() 

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

151 approximation = self.approximation.lower() 

152 if approximation == 'profeta': 

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

154 elif approximation == 'placzek': 

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

156 elif approximation == 'p-p': 

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

158 else: 

159 raise RuntimeError( 

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

161 self.approximation)) 

162 

163 return self.map_to_modes(Vel_rcc)