Coverage for /builds/kinetik161/ase/ase/vibrations/raman.py: 89.95%

1import numpy as np 

2 

3import ase.units as u 

4from ase.dft import monkhorst_pack 

5from ase.parallel import world 

6from ase.phonons import Phonons 

7from ase.utils import IOContext 

8from ase.vibrations.vibrations import AtomicDisplacements, Vibrations 

9 

10 

11class RamanCalculatorBase(IOContext): 

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

13 *args, 

14 name='raman', 

15 exext='.alpha', 

16 txt='-', 

17 verbose=False, 

18 comm=world, 

19 **kwargs): 

20 """ 

21 Parameters 

22 ---------- 

23 atoms: ase Atoms object 

24 exext: string 

25 Extension for excitation filenames 

26 txt: 

27 Output stream 

28 verbose: 

29 Verbosity level of output 

30 comm: 

31 Communicator, default world 

32 """ 

33 kwargs['name'] = name 

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

35 

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

37 

38 self.exext = exext 

39 

40 self.txt = self.openfile(txt, comm) 

41 self.verbose = verbose 

42 

43 self.comm = comm 

44 

45 

46class StaticRamanCalculatorBase(RamanCalculatorBase): 

47 """Base class for Raman intensities derived from 

48 static polarizabilities""" 

49 

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

51 self.exobj = exobj 

52 if exkwargs is None: 

53 exkwargs = {} 

54 self.exkwargs = exkwargs 

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

56 

57 def _new_exobj(self): 

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

59 

60 def calculate(self, atoms, disp): 

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

62 disp.calculate_and_save_static_polarizability(atoms) 

63 return returnvalue 

64 

65 

66class StaticRamanCalculator(StaticRamanCalculatorBase, Vibrations): 

67 pass 

68 

69 

70class StaticRamanPhononsCalculator(StaticRamanCalculatorBase, Phonons): 

71 pass 

72 

73 

74class RamanBase(AtomicDisplacements, IOContext): 

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

76 *args, 

77 name='raman', 

78 exname=None, 

79 exext='.alpha', 

80 txt='-', 

81 verbose=False, 

82 comm=world, 

83 **kwargs): 

84 """ 

85 Parameters 

86 ---------- 

87 atoms: ase Atoms object 

88 exext: string 

89 Extension for excitation filenames 

90 txt: 

91 Output stream 

92 verbose: 

93 Verbosity level of output 

94 comm: 

95 Communicator, default world 

96 """ 

97 self.atoms = atoms 

98 

99 self.name = name 

100 if exname is None: 

101 self.exname = name 

102 else: 

103 self.exname = exname 

104 self.exext = exext 

105 

106 self.txt = self.openfile(txt, comm) 

107 self.verbose = verbose 

108 

109 self.comm = comm 

110 

111 

112class RamanData(RamanBase): 

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

114 

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

116 *args, 

117 exname=None, # name for excited state calculations 

118 **kwargs): 

119 """ 

120 Parameters 

121 ---------- 

122 atoms: ase Atoms object 

123 exname: string 

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

125 used for reading excitations 

126 """ 

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

128 

129 if exname is None: 

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

131 self.exname = exname 

132 self._already_read = False 

133 

134 def get_energies(self): 

135 self.calculate_energies_and_modes() 

136 return self.om_Q 

137 

138 def init_parallel_read(self): 

139 """Initialize variables for parallel read""" 

140 rank = self.comm.rank 

141 indices = self.indices 

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

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

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

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

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

147 self.mynd = len(self.myr) 

148 

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

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

151 if self._already_read: 

152 return 

153 

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

155 self.init_parallel_read() 

156 self.read_excitations() 

157 

158 self._already_read = True 

159 

160 @staticmethod 

161 def m2(z): 

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

163 

164 def map_to_modes(self, V_rcc): 

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

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

167 return V_Qcc 

168 

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

170 """Full matrix element 

171 

172 Returns 

173 ------- 

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

175 """ 

176 # Angstrom^2 / sqrt(amu) 

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

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

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

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

181 

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

183 """Absolute Raman intensity or Raman scattering factor 

184 

185 Parameter 

186 --------- 

187 delta: float 

188 pre-factor for asymmetric anisotropy, default 0 

189 

190 References 

191 ---------- 

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

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

194 

195 Returns 

196 ------- 

197 raman intensity, unit Ang**4/amu 

198 """ 

199 alpha2_r, gamma2_r, delta2_r = self._invariants( 

200 self.electronic_me_Qcc(**kwargs)) 

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

202 

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

204 """Raman intensity 

205 

206 Returns 

207 ------- 

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

209 """ 

210 self.calculate_energies_and_modes() 

211 

212 m2 = Raman.m2 

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

214 if not self.observation: # XXXX remove 

215 """Simple sum, maybe too simple""" 

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

217 # XXX enable when appropriate 

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

219 # raise NotImplementedError('not yet') 

220 

221 # random orientation of the molecular frame 

222 # Woodward & Long, 

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

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

225 

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

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

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

229 return gamma2_r / 15. 

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

231 # scattered light in direction of incoming light 

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

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

234 # scattered light perendicular and 

235 # polarization in plane 

236 return 6 * gamma2_r / 45. 

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

238 # scattered light perendicular and 

239 # polarization out of plane 

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

241 else: 

242 raise NotImplementedError 

243 

244 def _invariants(self, alpha_Qcc): 

245 """Raman invariants 

246 

247 Parameter 

248 --------- 

249 alpha_Qcc: array 

250 Matrix element or polarizability tensor 

251 

252 Reference 

253 --------- 

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

255 

256 Returns 

257 ------- 

258 mean polarizability, anisotropy, asymmetric anisotropy 

259 """ 

260 m2 = Raman.m2 

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

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

263 delta2_r = 3 / 4. * ( 

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

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

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

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

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

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

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

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

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

273 return alpha2_r, gamma2_r, delta2_r 

274 

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

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

277 with IOContext() as io: 

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

279 return self._summary(log) 

280 

281 def _summary(self, log): 

282 hnu = self.get_energies() 

283 intensities = self.get_absolute_intensities() 

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

285 scale = 10**(-te) 

286 

287 if not te: 

288 ts = '' 

289 elif te > -2 and te < 3: 

290 ts = str(10**te) 

291 else: 

292 ts = f'10^{te}' 

293 

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

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

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

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

298 for n, e in enumerate(hnu): 

299 if e.imag != 0: 

300 c = 'i' 

301 e = e.imag 

302 else: 

303 c = ' ' 

304 e = e.real 

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

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

307 file=log) 

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

309 # XXX enable this in phonons 

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

311 # self.vibrations.get_zero_point_energy(), 

312 # file=log) 

313 

314 

315class Raman(RamanData): 

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

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

318 

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

320 kwargs.pop(key, None) 

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

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

323 

324 self.delta = self.vibrations.delta 

325 self.indices = self.vibrations.indices 

326 

327 def calculate_energies_and_modes(self): 

328 if hasattr(self, 'im_r'): 

329 return 

330 

331 self.read() 

332 

333 self.im_r = self.vibrations.im 

334 self.modes_Qq = self.vibrations.modes 

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

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

337 # pre-factors for one vibrational excitation 

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

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

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

341 # -> sqrt(amu) * Angstrom 

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

343 

344 

345class RamanPhonons(RamanData): 

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

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

348 

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

350 kwargs.pop(key, None) 

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

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

353 

354 self.delta = self.vibrations.delta 

355 self.indices = self.vibrations.indices 

356 

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

358 

359 @property 

360 def kpts(self): 

361 return self._kpts 

362 

363 @kpts.setter 

364 def kpts(self, kpts): 

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

366 self._kpts = kpts 

367 self.kpts_kc = monkhorst_pack(self.kpts) 

368 if hasattr(self, 'im_r'): 

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

370 

371 def calculate_energies_and_modes(self): 

372 if not self._already_read: 

373 if hasattr(self, 'im_r'): 

374 del self.im_r 

375 self.read() 

376 

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

378 omega_kl, u_kl = self.vibrations.band_structure( 

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

380 

381 self.im_r = self.vibrations.m_inv_x 

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

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

384 3 * len(self.atoms)) 

385 self.modes_Qq /= self.im_r 

386 self.om_v = self.om_Q 

387 

388 # pre-factors for one vibrational excitation 

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

390 self.vib01_Q = np.where( 

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

392 # -> sqrt(amu) * Angstrom 

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