Coverage for /builds/ase/ase/ase/vibrations/raman.py: 89.47%

1# fmt: off 

2 

3import numpy as np 

4 

5import ase.units as u 

6from ase.dft import monkhorst_pack 

7from ase.parallel import world 

8from ase.phonons import Phonons 

9from ase.utils import IOContext 

10from ase.vibrations.vibrations import AtomicDisplacements, Vibrations 

11 

12 

13class RamanCalculatorBase(IOContext): 

14 def __init__(self, atoms, # XXX do we need atoms at this stage ? 

15 *args, 

16 name='raman', 

17 exext='.alpha', 

18 txt='-', 

19 verbose=False, 

20 comm=world, 

21 **kwargs): 

22 """ 

23 Parameters 

24 ---------- 

25 atoms: ase Atoms object 

26 exext: string 

27 Extension for excitation filenames 

28 txt: 

29 Output stream 

30 verbose: 

31 Verbosity level of output 

32 comm: 

33 Communicator, default world 

34 """ 

35 kwargs['name'] = name 

36 self.exname = kwargs.pop('exname', name) 

37 

38 super().__init__(atoms, *args, **kwargs) 

39 

40 self.exext = exext 

41 

42 self.txt = self.openfile(file=txt, comm=comm) 

43 self.verbose = verbose 

44 

45 self.comm = comm 

46 

47 

48class StaticRamanCalculatorBase(RamanCalculatorBase): 

49 """Base class for Raman intensities derived from 

50 static polarizabilities""" 

51 

52 def __init__(self, atoms, exobj, exkwargs=None, *args, **kwargs): 

53 self.exobj = exobj 

54 if exkwargs is None: 

55 exkwargs = {} 

56 self.exkwargs = exkwargs 

57 super().__init__(atoms, *args, **kwargs) 

58 

59 def _new_exobj(self): 

60 return self.exobj(**self.exkwargs) 

61 

62 def calculate(self, atoms, disp): 

63 returnvalue = super().calculate(atoms, disp) 

64 disp.calculate_and_save_static_polarizability(atoms) 

65 return returnvalue 

66 

67 

68class StaticRamanCalculator(StaticRamanCalculatorBase, Vibrations): 

69 pass 

70 

71 

72class StaticRamanPhononsCalculator(StaticRamanCalculatorBase, Phonons): 

73 pass 

74 

75 

76class RamanBase(AtomicDisplacements, IOContext): 

77 def __init__(self, atoms, # XXX do we need atoms at this stage ? 

78 *args, 

79 name='raman', 

80 exname=None, 

81 exext='.alpha', 

82 txt='-', 

83 verbose=False, 

84 comm=world, 

85 **kwargs): 

86 """ 

87 Parameters 

88 ---------- 

89 atoms: ase Atoms object 

90 exext: string 

91 Extension for excitation filenames 

92 txt: 

93 Output stream 

94 verbose: 

95 Verbosity level of output 

96 comm: 

97 Communicator, default world 

98 """ 

99 self.atoms = atoms 

100 

101 self.name = name 

102 if exname is None: 

103 self.exname = name 

104 else: 

105 self.exname = exname 

106 self.exext = exext 

107 

108 self.txt = self.openfile(file=txt, comm=comm) 

109 self.verbose = verbose 

110 

111 self.comm = comm 

112 

113 

114class RamanData(RamanBase): 

115 """Base class to evaluate Raman spectra from pre-computed data""" 

116 

117 def __init__(self, atoms, # XXX do we need atoms at this stage ? 

118 *args, 

119 exname=None, # name for excited state calculations 

120 **kwargs): 

121 """ 

122 Parameters 

123 ---------- 

124 atoms: ase Atoms object 

125 exname: string 

126 name for excited state calculations (defaults to name), 

127 used for reading excitations 

128 """ 

129 super().__init__(atoms, *args, **kwargs) 

130 

131 if exname is None: 

132 exname = kwargs.get('name', self.name) 

133 self.exname = exname 

134 self._already_read = False 

135 

136 def get_energies(self): 

137 self.calculate_energies_and_modes() 

138 return self.om_Q 

139 

140 def init_parallel_read(self): 

141 """Initialize variables for parallel read""" 

142 rank = self.comm.rank 

143 indices = self.indices 

144 myn = -(-self.ndof // self.comm.size) # ceil divide 

145 self.slize = s = slice(myn * rank, myn * (rank + 1)) 

146 self.myindices = np.repeat(indices, 3)[s] 

147 self.myxyz = ('xyz' * len(indices))[s] 

148 self.myr = range(self.ndof)[s] 

149 self.mynd = len(self.myr) 

150 

151 def read(self, *args, **kwargs): 

152 """Read data from a pre-performed calculation.""" 

153 if self._already_read: 

154 return 

155 

156 self.vibrations.read(*args, **kwargs) 

157 self.init_parallel_read() 

158 self.read_excitations() 

159 

160 self._already_read = True 

161 

162 @staticmethod 

163 def m2(z): 

164 return (z * z.conj()).real 

165 

166 def map_to_modes(self, V_rcc): 

167 V_qcc = (V_rcc.T * self.im_r).T # units Angstrom^2 / sqrt(amu) 

168 V_Qcc = np.dot(V_qcc.T, self.modes_Qq.T).T 

169 return V_Qcc 

170 

171 def me_Qcc(self, *args, **kwargs): 

172 """Full matrix element 

173 

174 Returns 

175 ------- 

176 Matrix element in e^2 Angstrom^2 / eV 

177 """ 

178 # Angstrom^2 / sqrt(amu) 

179 elme_Qcc = self.electronic_me_Qcc(*args, **kwargs) 

180 # Angstrom^3 -> e^2 Angstrom^2 / eV 

181 elme_Qcc /= u.Hartree * u.Bohr # e^2 Angstrom / eV / sqrt(amu) 

182 return elme_Qcc * self.vib01_Q[:, None, None] 

183 

184 def get_absolute_intensities(self, delta=0, **kwargs): 

185 """Absolute Raman intensity or Raman scattering factor 

186 

187 Parameter 

188 --------- 

189 delta: float 

190 pre-factor for asymmetric anisotropy, default 0 

191 

192 References 

193 ---------- 

194 Porezag and Pederson, PRB 54 (1996) 7830-7836 (delta=0) 

195 Baiardi and Barone, JCTC 11 (2015) 3267-3280 (delta=5) 

196 

197 Returns 

198 ------- 

199 raman intensity, unit Ang**4/amu 

200 """ 

201 alpha2_r, gamma2_r, delta2_r = self._invariants( 

202 self.electronic_me_Qcc(**kwargs)) 

203 return 45 * alpha2_r + delta * delta2_r + 7 * gamma2_r 

204 

205 def intensity(self, *args, **kwargs): 

206 """Raman intensity 

207 

208 Returns 

209 ------- 

210 unit e^4 Angstrom^4 / eV^2 

211 """ 

212 self.calculate_energies_and_modes() 

213 

214 m2 = Raman.m2 

215 alpha_Qcc = self.me_Qcc(*args, **kwargs) 

216 if not self.observation: # XXXX remove 

217 """Simple sum, maybe too simple""" 

218 return m2(alpha_Qcc).sum(axis=1).sum(axis=1) 

219 # XXX enable when appropriate 

220 # if self.observation['orientation'].lower() != 'random': 

221 # raise NotImplementedError('not yet') 

222 

223 # random orientation of the molecular frame 

224 # Woodward & Long, 

225 # Guthmuller, J. J. Chem. Phys. 2016, 144 (6), 64106 

226 alpha2_r, gamma2_r, delta2_r = self._invariants(alpha_Qcc) 

227 

228 if self.observation['geometry'] == '-Z(XX)Z': # Porto's notation 

229 return (45 * alpha2_r + 5 * delta2_r + 4 * gamma2_r) / 45. 

230 elif self.observation['geometry'] == '-Z(XY)Z': # Porto's notation 

231 return gamma2_r / 15. 

232 elif self.observation['scattered'] == 'Z': 

233 # scattered light in direction of incoming light 

234 return (45 * alpha2_r + 5 * delta2_r + 7 * gamma2_r) / 45. 

235 elif self.observation['scattered'] == 'parallel': 

236 # scattered light perendicular and 

237 # polarization in plane 

238 return 6 * gamma2_r / 45. 

239 elif self.observation['scattered'] == 'perpendicular': 

240 # scattered light perendicular and 

241 # polarization out of plane 

242 return (45 * alpha2_r + 5 * delta2_r + 7 * gamma2_r) / 45. 

243 else: 

244 raise NotImplementedError 

245 

246 def _invariants(self, alpha_Qcc): 

247 """Raman invariants 

248 

249 Parameter 

250 --------- 

251 alpha_Qcc: array 

252 Matrix element or polarizability tensor 

253 

254 Reference 

255 --------- 

256 Derek A. Long, The Raman Effect, ISBN 0-471-49028-8 

257 

258 Returns 

259 ------- 

260 mean polarizability, anisotropy, asymmetric anisotropy 

261 """ 

262 m2 = Raman.m2 

263 alpha2_r = m2(alpha_Qcc[:, 0, 0] + alpha_Qcc[:, 1, 1] + 

264 alpha_Qcc[:, 2, 2]) / 9. 

265 delta2_r = 3 / 4. * ( 

266 m2(alpha_Qcc[:, 0, 1] - alpha_Qcc[:, 1, 0]) + 

267 m2(alpha_Qcc[:, 0, 2] - alpha_Qcc[:, 2, 0]) + 

268 m2(alpha_Qcc[:, 1, 2] - alpha_Qcc[:, 2, 1])) 

269 gamma2_r = (3 / 4. * (m2(alpha_Qcc[:, 0, 1] + alpha_Qcc[:, 1, 0]) + 

270 m2(alpha_Qcc[:, 0, 2] + alpha_Qcc[:, 2, 0]) + 

271 m2(alpha_Qcc[:, 1, 2] + alpha_Qcc[:, 2, 1])) + 

272 (m2(alpha_Qcc[:, 0, 0] - alpha_Qcc[:, 1, 1]) + 

273 m2(alpha_Qcc[:, 0, 0] - alpha_Qcc[:, 2, 2]) + 

274 m2(alpha_Qcc[:, 1, 1] - alpha_Qcc[:, 2, 2])) / 2) 

275 return alpha2_r, gamma2_r, delta2_r 

276 

277 def summary(self, log='-'): 

278 """Print summary for given omega [eV]""" 

279 with IOContext() as io: 

280 log = io.openfile(file=log, mode='a', comm=self.comm) 

281 return self._summary(log) 

282 

283 def _summary(self, log): 

284 hnu = self.get_energies() 

285 intensities = self.get_absolute_intensities() 

286 te = int(np.log10(intensities.max())) - 2 

287 scale = 10**(-te) 

288 

289 if not te: 

290 ts = '' 

291 elif te > -2 and te < 3: 

292 ts = str(10**te) 

293 else: 

294 ts = f'10^{te}' 

295 

296 print('-------------------------------------', file=log) 

297 print(' Mode Frequency Intensity', file=log) 

298 print(f' # meV cm^-1 [{ts}A^4/amu]', file=log) 

299 print('-------------------------------------', file=log) 

300 for n, e in enumerate(hnu): 

301 if e.imag != 0: 

302 c = 'i' 

303 e = e.imag 

304 else: 

305 c = ' ' 

306 e = e.real 

307 print('%3d %6.1f%s %7.1f%s %9.2f' % 

308 (n, 1000 * e, c, e / u.invcm, c, intensities[n] * scale), 

309 file=log) 

310 print('-------------------------------------', file=log) 

311 # XXX enable this in phonons 

312 # parprint('Zero-point energy: %.3f eV' % 

313 # self.vibrations.get_zero_point_energy(), 

314 # file=log) 

315 

316 

317class Raman(RamanData): 

318 def __init__(self, atoms, *args, **kwargs): 

319 super().__init__(atoms, *args, **kwargs) 

320 

321 for key in ['txt', 'exext', 'exname']: 

322 kwargs.pop(key, None) 

323 kwargs['name'] = kwargs.get('name', self.name) 

324 self.vibrations = Vibrations(atoms, *args, **kwargs) 

325 

326 self.delta = self.vibrations.delta 

327 self.indices = self.vibrations.indices 

328 

329 def calculate_energies_and_modes(self): 

330 if hasattr(self, 'im_r'): 

331 return 

332 

333 self.read() 

334 

335 self.im_r = self.vibrations.im 

336 self.modes_Qq = self.vibrations.modes 

337 self.om_Q = self.vibrations.hnu.real # energies in eV 

338 self.H = self.vibrations.H # XXX used in albrecht.py 

339 # pre-factors for one vibrational excitation 

340 with np.errstate(divide='ignore'): 

341 self.vib01_Q = np.where(self.om_Q > 0, 

342 1. / np.sqrt(2 * self.om_Q), 0) 

343 # -> sqrt(amu) * Angstrom 

344 self.vib01_Q *= np.sqrt(u.Ha * u._me / u._amu) * u.Bohr 

345 

346 

347class RamanPhonons(RamanData): 

348 def __init__(self, atoms, *args, **kwargs): 

349 RamanData.__init__(self, atoms, *args, **kwargs) 

350 

351 for key in ['txt', 'exext', 'exname']: 

352 kwargs.pop(key, None) 

353 kwargs['name'] = kwargs.get('name', self.name) 

354 self.vibrations = Phonons(atoms, *args, **kwargs) 

355 

356 self.delta = self.vibrations.delta 

357 self.indices = self.vibrations.indices 

358 

359 self.kpts = (1, 1, 1) 

360 

361 @property 

362 def kpts(self): 

363 return self._kpts 

364 

365 @kpts.setter 

366 def kpts(self, kpts): 

367 if not hasattr(self, '_kpts') or kpts != self._kpts: 

368 self._kpts = kpts 

369 self.kpts_kc = monkhorst_pack(self.kpts) 

370 if hasattr(self, 'im_r'): 

371 del self.im_r # we'll have to recalculate 

372 

373 def calculate_energies_and_modes(self): 

374 if not self._already_read: 

375 if hasattr(self, 'im_r'): 

376 del self.im_r 

377 self.read() 

378 

379 if not hasattr(self, 'im_r'): 

380 omega_kl, u_kl = self.vibrations.band_structure( 

381 self.kpts_kc, modes=True, verbose=self.verbose) 

382 

383 self.im_r = self.vibrations.m_inv_x 

384 self.om_Q = omega_kl.ravel().real # energies in eV 

385 self.modes_Qq = u_kl.reshape(len(self.om_Q), 

386 3 * len(self.atoms)) 

387 self.modes_Qq /= self.im_r 

388 self.om_v = self.om_Q 

389 

390 # pre-factors for one vibrational excitation 

391 with np.errstate(divide='ignore', invalid='ignore'): 

392 self.vib01_Q = np.where( 

393 self.om_Q > 0, 1. / np.sqrt(2 * self.om_Q), 0) 

394 # -> sqrt(amu) * Angstrom 

395 self.vib01_Q *= np.sqrt(u.Ha * u._me / u._amu) * u.Bohr