Coverage for gpaw/lrtddft/omega_matrix.py: 84%
414 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1import numpy as np
2from scipy.linalg import eigh
4from ase.units import Hartree
5from ase.utils.timing import Timer
7import gpaw.mpi as mpi
8from .kssingle import KSSingles
9from gpaw.setup import CachedYukawaInteractions
10from gpaw.transformers import Transformer
11from gpaw.utilities import pack_density
12from gpaw.xc import XC
14"""This module defines a Omega Matrix class."""
17class OmegaMatrix:
19 """
20 Omega matrix in Casidas linear response formalism
22 Parameters
23 - calculator: the calculator object the ground state calculation
24 - kss: the Kohn-Sham singles object
25 - xc: the exchange correlation approx. to use
26 - derivativeLevel: which level i of d^i Exc/dn^i to use
27 - numscale: numeric epsilon for derivativeLevel=0,1
28 - filehandle: the oject can be read from a filehandle
29 - txt: output stream or file name
30 - finegrid: level of fine grid to use. 0: nothing, 1 for poisson only,
31 2 everything on the fine grid
32 """
34 def __init__(self,
35 calculator=None,
36 kss=None,
37 xc=None,
38 derivativeLevel=None,
39 numscale=0.001,
40 poisson=None,
41 filehandle=None,
42 log=None,
43 finegrid=2,
44 eh_comm=None):
45 self.log = log
47 if eh_comm is None:
48 eh_comm = mpi.serial_comm
50 self.eh_comm = eh_comm
51 self.fullkss = kss
53 if filehandle is not None:
54 self.read(fh=filehandle)
55 return None
57 self.finegrid = finegrid
59 if calculator is None:
60 return
62 self.paw = calculator
63 wfs = self.paw.wfs
65 # handle different grid possibilities
66 self.restrict = None
67 # self.poisson = PoissonSolver(nn=self.paw.hamiltonian.poisson.nn)
68 if poisson is None:
69 self.poisson = calculator.hamiltonian.poisson
70 else:
71 self.poisson = poisson
72 if finegrid:
73 self.poisson.set_grid_descriptor(self.paw.density.finegd)
75 self.gd = self.paw.density.finegd
76 if finegrid == 1:
77 self.gd = wfs.gd
78 else:
79 self.poisson.set_grid_descriptor(wfs.gd)
80 self.gd = wfs.gd
81 self.restrict = Transformer(self.paw.density.finegd, wfs.gd,
82 self.paw.density.stencil).apply
84 if xc == 'RPA':
85 xc = None # enable RPA as keyword
86 if xc is not None:
87 if isinstance(xc, str):
88 self.xc = XC(xc)
89 else:
90 self.xc = xc
91 self.xc.initialize(self.paw.density, self.paw.hamiltonian, wfs)
93 # check derivativeLevel
94 if derivativeLevel is None:
95 derivativeLevel = \
96 self.xc.get_functional().get_max_derivative_level()
97 self.derivativeLevel = derivativeLevel
98 else:
99 self.xc = None
101 if getattr(self.xc, 'omega', 0): # can be None or int
102 self.yukawa_interactions = CachedYukawaInteractions(self.xc.omega)
103 else:
104 self.yukawa_interactions = None
106 self.numscale = numscale
108 self.singletsinglet = False
109 if kss.nvspins < 2 and kss.npspins < 2:
110 # this will be a singlet to singlet calculation only
111 self.singletsinglet = True
113 nij = len(kss)
114 self.Om = np.zeros((nij, nij))
115 self.get_full()
117 def get_full(self):
118 """Evaluate full omega matrix"""
119 self.paw.timer.start('Omega RPA')
120 self.log()
121 self.log('Linear response TDDFT calculation')
122 self.log()
123 self.get_rpa()
124 self.paw.timer.stop()
126 if self.xc is not None:
127 self.paw.timer.start('Omega XC')
128 self.get_xc()
129 self.paw.timer.stop()
131 self.eh_comm.sum(self.Om)
132 self.full = self.Om
134 def get_xc(self):
135 """Add xc part of the coupling matrix"""
137 # shorthands
138 paw = self.paw
139 wfs = paw.wfs
140 gd = paw.density.finegd
141 eh_comm = self.eh_comm
143 fg = self.finegrid == 2
144 kss = self.fullkss
145 nij = len(kss)
147 Om_xc = self.Om
148 # initialize densities
149 # nt_sg is the smooth density on the fine grid with spin index
151 if kss.nvspins == 2:
152 # spin polarised ground state calc.
153 nt_sg = paw.density.nt_sg
154 else:
155 # spin unpolarised ground state calc.
156 if kss.npspins == 2:
157 # construct spin polarised densities
158 nt_sg = np.array([.5 * paw.density.nt_sg[0],
159 .5 * paw.density.nt_sg[0]])
160 else:
161 nt_sg = paw.density.nt_sg
162 # check if D_sp have been changed before
163 D_asp = {}
164 for a, D_sp in self.paw.density.D_asp.items():
165 if len(D_sp) != kss.npspins:
166 if len(D_sp) == 1:
167 D_asp[a] = np.array([0.5 * D_sp[0], 0.5 * D_sp[0]])
168 else:
169 D_asp[a] = np.array([D_sp[0] + D_sp[1]])
170 else:
171 D_asp[a] = D_sp.copy()
173 # restrict the density if needed
174 if fg:
175 nt_s = nt_sg
176 else:
177 nt_s = self.gd.zeros(nt_sg.shape[0])
178 for s in range(nt_sg.shape[0]):
179 self.restrict(nt_sg[s], nt_s[s])
180 gd = paw.density.gd
182 # initialize vxc or fxc
184 if self.derivativeLevel == 0:
185 raise NotImplementedError
186 if kss.npspins == 2:
187 v_g = nt_sg[0].copy()
188 else:
189 v_g = nt_sg.copy()
190 elif self.derivativeLevel == 1:
191 pass
192 elif self.derivativeLevel == 2:
193 fxc_sg = np.zeros(nt_sg.shape)
194 self.xc.calculate_fxc(gd, nt_sg, fxc_sg)
195 else:
196 raise ValueError('derivativeLevel can only be 0,1,2')
198 ns = self.numscale
199 xc = self.xc
200 self.log('XC', nij, 'transitions')
201 for ij in range(eh_comm.rank, nij, eh_comm.size):
202 self.log('XC kss[' + '%d' % ij + ']')
204 timer = Timer()
205 timer.start('init')
206 timer2 = Timer()
208 if self.derivativeLevel >= 1:
209 # vxc is available
210 # We use the numerical two point formula for calculating
211 # the integral over fxc*n_ij. The results are
212 # vvt_s smooth integral
213 # nucleus.I_sp atom based correction matrices (pack_hermitian)
214 # stored on each nucleus
215 timer2.start('init v grids')
216 vp_s = np.zeros(nt_s.shape, nt_s.dtype.char)
217 vm_s = np.zeros(nt_s.shape, nt_s.dtype.char)
218 if kss.npspins == 2: # spin polarised
219 nv_s = nt_s.copy()
220 nv_s[kss[ij].pspin] += ns * kss[ij].get(fg)
221 xc.calculate(gd, nv_s, vp_s)
222 nv_s = nt_s.copy()
223 nv_s[kss[ij].pspin] -= ns * kss[ij].get(fg)
224 xc.calculate(gd, nv_s, vm_s)
225 else: # spin unpolarised
226 nv = nt_s + ns * kss[ij].get(fg)
227 xc.calculate(gd, nv, vp_s)
228 nv = nt_s - ns * kss[ij].get(fg)
229 xc.calculate(gd, nv, vm_s)
230 vvt_s = (0.5 / ns) * (vp_s - vm_s)
231 timer2.stop()
233 # initialize the correction matrices
234 timer2.start('init v corrections')
235 I_asp = {}
236 for a, P_ni in wfs.kpt_u[kss[ij].spin].P_ani.items():
237 # create the modified density matrix
238 Pi_i = P_ni[kss[ij].i]
239 Pj_i = P_ni[kss[ij].j]
240 P_ii = np.outer(Pi_i, Pj_i)
241 # we need the symmetric form, hence we can pack
242 P_p = pack_density(P_ii)
243 D_sp = D_asp[a].copy()
244 D_sp[kss[ij].pspin] -= ns * P_p
245 setup = wfs.setups[a]
246 I_sp = np.zeros_like(D_sp)
247 self.xc.calculate_paw_correction(setup, D_sp, I_sp)
248 I_sp *= -1.0
249 D_sp = D_asp[a].copy()
250 D_sp[kss[ij].pspin] += ns * P_p
251 self.xc.calculate_paw_correction(setup, D_sp, I_sp)
252 I_sp /= 2.0 * ns
253 I_asp[a] = I_sp
254 timer2.stop()
256 timer.stop()
257 t0 = timer.get_time('init')
258 timer.start(ij)
260 for kq in range(ij, nij):
261 weight = self.weight_Kijkq(ij, kq)
263 if self.derivativeLevel == 0:
264 # only Exc is available
266 if kss.npspins == 2: # spin polarised
267 nv_g = nt_sg.copy()
268 nv_g[kss[ij].pspin] += kss[ij].get(fg)
269 nv_g[kss[kq].pspin] += kss[kq].get(fg)
270 Excpp = xc.get_energy_and_potential(
271 nv_g[0], v_g, nv_g[1], v_g)
272 nv_g = nt_sg.copy()
273 nv_g[kss[ij].pspin] += kss[ij].get(fg)
274 nv_g[kss[kq].pspin] -= kss[kq].get(fg)
275 Excpm = xc.get_energy_and_potential(
276 nv_g[0], v_g, nv_g[1], v_g)
277 nv_g = nt_sg.copy()
278 nv_g[kss[ij].pspin] -=\
279 kss[ij].get(fg)
280 nv_g[kss[kq].pspin] +=\
281 kss[kq].get(fg)
282 Excmp = xc.get_energy_and_potential(
283 nv_g[0], v_g, nv_g[1], v_g)
284 nv_g = nt_sg.copy()
285 nv_g[kss[ij].pspin] -= \
286 kss[ij].get(fg)
287 nv_g[kss[kq].pspin] -=\
288 kss[kq].get(fg)
289 Excpp = xc.get_energy_and_potential(
290 nv_g[0], v_g, nv_g[1], v_g)
291 else: # spin unpolarised
292 nv_g = nt_sg + ns * kss[ij].get(fg)\
293 + ns * kss[kq].get(fg)
294 Excpp = xc.get_energy_and_potential(nv_g, v_g)
295 nv_g = nt_sg + ns * kss[ij].get(fg)\
296 - ns * kss[kq].get(fg)
297 Excpm = xc.get_energy_and_potential(nv_g, v_g)
298 nv_g = nt_sg - ns * kss[ij].get(fg)\
299 + ns * kss[kq].get(fg)
300 Excmp = xc.get_energy_and_potential(nv_g, v_g)
301 nv_g = nt_sg - ns * kss[ij].get(fg)\
302 - ns * kss[kq].get(fg)
303 Excmm = xc.get_energy_and_potential(nv_g, v_g)
305 Om_xc[ij, kq] += weight *\
306 0.25 * \
307 (Excpp - Excmp - Excpm + Excmm) / (ns * ns)
309 elif self.derivativeLevel == 1:
310 # vxc is available
312 timer2.start('integrate')
313 Om_xc[ij, kq] += weight *\
314 self.gd.integrate(kss[kq].get(fg) *
315 vvt_s[kss[kq].pspin])
316 timer2.stop()
318 timer2.start('integrate corrections')
319 Exc = 0.
320 for a, P_ni in wfs.kpt_u[kss[kq].spin].P_ani.items():
321 # create the modified density matrix
322 Pk_i = P_ni[kss[kq].i]
323 Pq_i = P_ni[kss[kq].j]
324 P_ii = np.outer(Pk_i, Pq_i)
325 # we need the symmetric form, hence we can pack
326 # use pack_density as I_sp used pack_hermitian
327 P_p = pack_density(P_ii)
328 Exc += np.dot(I_asp[a][kss[kq].pspin], P_p)
329 Om_xc[ij, kq] += weight * self.gd.comm.sum_scalar(Exc)
330 timer2.stop()
332 elif self.derivativeLevel == 2:
333 # fxc is available
334 if kss.npspins == 2: # spin polarised
335 Om_xc[ij, kq] += weight *\
336 gd.integrate(kss[ij].get(fg) *
337 kss[kq].get(fg) *
338 fxc_sg[kss[ij].pspin, kss[kq].pspin])
339 else: # spin unpolarised
340 Om_xc[ij, kq] += weight *\
341 gd.integrate(kss[ij].get(fg) *
342 kss[kq].get(fg) *
343 fxc_sg)
345 # XXX still numeric derivatives for local terms
346 timer2.start('integrate corrections')
347 Exc = 0.
348 for a, P_ni in wfs.kpt_u[kss[kq].spin].P_ani.items():
349 # create the modified density matrix
350 Pk_i = P_ni[kss[kq].i]
351 Pq_i = P_ni[kss[kq].j]
352 P_ii = np.outer(Pk_i, Pq_i)
353 # we need the symmetric form, hence we can pack
354 # use pack_density as I_sp used pack_hermitian
355 P_p = pack_density(P_ii)
356 Exc += np.dot(I_asp[a][kss[kq].pspin], P_p)
357 Om_xc[ij, kq] += weight * self.gd.comm.sum(Exc)
358 timer2.stop()
360 if ij != kq:
361 Om_xc[kq, ij] = Om_xc[ij, kq]
363 timer.stop()
364# timer2.write()
365 if ij < (nij - 1):
366 self.log('XC estimated time left',
367 self.time_left(timer, t0, ij, nij))
369 def Coulomb_integral_kss(self, kss_ij, kss_kq, phit, rhot,
370 timer=None, yukawa=False):
371 # smooth part
372 if timer:
373 timer.start('integrate')
374 I = self.gd.integrate(rhot * phit)
375 if timer:
376 timer.stop()
377 timer.start('integrate corrections 2')
379 wfs = self.paw.wfs
380 Pij_ani = wfs.kpt_u[kss_ij.spin].P_ani
381 Pkq_ani = wfs.kpt_u[kss_kq.spin].P_ani
383 # Add atomic corrections
384 Ia = 0.0
385 for a, Pij_ni in Pij_ani.items():
386 Pi_i = Pij_ni[kss_ij.i]
387 Pj_i = Pij_ni[kss_ij.j]
388 Dij_ii = np.outer(Pi_i, Pj_i)
389 Dij_p = pack_density(Dij_ii)
390 Pk_i = Pkq_ani[a][kss_kq.i]
391 Pq_i = Pkq_ani[a][kss_kq.j]
392 Dkq_ii = np.outer(Pk_i, Pq_i)
393 Dkq_p = pack_density(Dkq_ii)
394 if yukawa:
395 assert abs(
396 self.yukawa_interactions.omega - self.xc.omega) < 1e-14
397 C_pp = self.yukawa_interactions.get_Mg_pp(wfs.setups[a])
398 else:
399 C_pp = wfs.setups[a].M_pp
400 # ----
401 # 2 > P P C P P
402 # ---- ip jr prst ks qt
403 # prst
404 Ia += 2.0 * np.dot(Dkq_p, np.dot(C_pp, Dij_p))
405 I += self.gd.comm.sum_scalar(Ia)
406 if timer:
407 timer.stop()
409 return I
411 def get_rpa(self):
412 """calculate RPA part of the omega matrix"""
414 # shorthands
415 kss = self.fullkss
416 finegrid = self.finegrid
417 eh_comm = self.eh_comm
419 # calculate omega matrix
420 nij = len(kss)
421 self.log('RPA', nij, 'transitions')
423 Om = self.Om
425 for ij in range(eh_comm.rank, nij, eh_comm.size):
426 self.log('RPA kss[' + '%d' % ij + ']=', kss[ij])
428 timer = Timer()
429 timer.start('init')
430 timer2 = Timer()
432 # smooth density including compensation charges
433 timer2.start('with_compensation_charges 0')
434 rhot_p = kss[ij].with_compensation_charges(
435 finegrid != 0)
436 timer2.stop()
438 # integrate with 1/|r_1-r_2|
439 timer2.start('poisson')
440 phit_p = np.zeros(rhot_p.shape, rhot_p.dtype.char)
441 self.poisson.solve(phit_p, rhot_p, charge=None)
442 timer2.stop()
444 timer.stop()
445 t0 = timer.get_time('init')
446 timer.start(ij)
448 if finegrid == 1:
449 rhot = kss[ij].with_compensation_charges()
450 phit = self.gd.zeros()
451 self.restrict(phit_p, phit)
452 else:
453 phit = phit_p
454 rhot = rhot_p
456 for kq in range(ij, nij):
457 if kq != ij:
458 # smooth density including compensation charges
459 timer2.start('kq with_compensation_charges')
460 rhot = kss[kq].with_compensation_charges(
461 finegrid == 2)
462 timer2.stop()
464 pre = 2 * np.sqrt(kss[ij].get_energy() *
465 kss[kq].get_energy() *
466 kss[ij].get_weight() *
467 kss[kq].get_weight())
468 I = self.Coulomb_integral_kss(kss[ij], kss[kq],
469 rhot, phit, timer2)
470 Om[ij, kq] = pre * I
472 if ij == kq:
473 Om[ij, kq] += kss[ij].get_energy() ** 2
474 else:
475 Om[kq, ij] = Om[ij, kq]
477 timer.stop()
478# timer2.write()
479 if ij < (nij - 1):
480 t = timer.get_time(ij) # time for nij-ij calculations
481 t = .5 * t * \
482 (nij - ij) # estimated time for n*(n+1)/2, n=nij-(ij+1)
483 self.log('RPA estimated time left',
484 self.timestring(t0 * (nij - ij - 1) + t))
486 def singlets_triplets(self):
487 """Split yourself into singlet and triplet transitions"""
489 assert self.fullkss.npspins == 2
490 assert self.fullkss.nvspins == 1
492 # strip kss from down spins
493 skss = KSSingles()
494 skss.dtype = self.fullkss.dtype
495 tkss = KSSingles()
496 tkss.dtype = self.fullkss.dtype
497 map = []
498 for ij, ks in enumerate(self.fullkss):
499 if ks.pspin == ks.spin:
500 skss.append((ks + ks) / np.sqrt(2))
501 tkss.append((ks - ks) / np.sqrt(2))
502 map.append(ij)
503 skss.istart = tkss.istart = self.fullkss.restrict['istart']
504 skss.jend = tkss.jend = self.fullkss.restrict['jend']
505 nkss = len(skss)
507 # define the singlet and the triplet omega-matrices
508 sOm = OmegaMatrix(kss=skss, log=self.log)
509 sOm.full = np.empty((nkss, nkss))
510 tOm = OmegaMatrix(kss=tkss, log=self.log)
511 tOm.full = np.empty((nkss, nkss))
512 for ij in range(nkss):
513 for kl in range(nkss):
514 sOm.full[ij, kl] = (self.full[map[ij], map[kl]] +
515 self.full[map[ij], nkss + map[kl]])
516 tOm.full[ij, kl] = (self.full[map[ij], map[kl]] -
517 self.full[map[ij], nkss + map[kl]])
518 return sOm, tOm
520 def timestring(self, t):
521 ti = int(t + 0.5)
522 td = ti // 86400
523 st = ''
524 if td > 0:
525 st += '%d' % td + 'd'
526 ti -= td * 86400
527 th = ti // 3600
528 if th > 0:
529 st += '%d' % th + 'h'
530 ti -= th * 3600
531 tm = ti // 60
532 if tm > 0:
533 st += '%d' % tm + 'm'
534 ti -= tm * 60
535 st += '%d' % ti + 's'
536 return st
538 def time_left(self, timer, t0, ij, nij):
539 t = timer.get_time(ij) # time for nij-ij calculations
540 t = .5 * t * (nij - ij) # estimated time for n*(n+1)/2, n=nij-(ij+1)
541 return self.timestring(t0 * (nij - ij - 1) + t)
543 def get_map(self, restrict=None):
544 """Return the reduction map for the given requirements
546 Returns
547 -------
548 map - list of original indices
549 kss - reduced KSSingles object
550 """
551 self.log('# diagonalize: %d transitions original'
552 % len(self.fullkss))
554 map = []
556 rst_dict = self.fullkss.restrict.values
557 if restrict is not None:
558 rst_dict.update(restrict)
559 kss = KSSingles(restrict=rst_dict)
560 kss.dtype = self.fullkss.dtype
562 for ij, k in zip(range(len(self.fullkss)), self.fullkss):
563 if kss.restrict.is_good(k):
564 kss.append(k)
565 map.append(ij)
566 kss.update()
567 self.log('# diagonalize: %d transitions now' % len(kss))
569 return map, kss
571 def diagonalize(self, restrict=None):
572 """Evaluate Eigenvectors and Eigenvalues:"""
573 map, kss = self.get_map(restrict)
575 nij = len(kss)
576 if map is None:
577 evec = self.full.copy()
578 else:
579 evec = np.zeros((nij, nij))
580 for ij in range(nij):
581 for kq in range(nij):
582 evec[ij, kq] = self.full[map[ij], map[kq]]
583 assert len(evec) > 0
585 self.eigenvalues, v = eigh(evec)
586 self.eigenvectors = v.T
587 self.kss = kss
589 @property
590 def kss(self):
591 if hasattr(self, '_kss'):
592 return self._kss
593 return self.fullkss
595 @kss.setter
596 def kss(self, kss):
597 """Set current (restricted) KSSingles object"""
598 self._kss = kss
600 def read(self, filename=None, fh=None):
601 """Read myself from a file"""
602 if fh is None:
603 f = open(filename)
604 else:
605 f = fh
607 f.readline()
608 nij = int(f.readline())
609 full = np.zeros((nij, nij))
610 for ij in range(nij):
611 l = [float(x) for x in f.readline().split()]
612 full[ij, ij:] = l
613 full[ij:, ij] = l
614 self.full = full
616 if fh is None:
617 f.close()
619 def write(self, filename=None, fh=None):
620 """Write current state to a file."""
621 try:
622 rank = self.paw.wfs.world.rank
623 except AttributeError:
624 rank = mpi.world.rank
625 if rank == 0:
626 if fh is None:
627 f = open(filename, 'w')
628 else:
629 f = fh
631 f.write('# OmegaMatrix\n')
632 nij = len(self.fullkss)
633 f.write('%d\n' % nij)
634 for ij in range(nij):
635 for kq in range(ij, nij):
636 f.write(' %g' % self.full[ij, kq])
637 f.write('\n')
639 if fh is None:
640 f.close()
642 def weight_Kijkq(self, ij, kq):
643 """weight for the coupling matrix terms"""
644 kss = self.fullkss
645 return 2. * np.sqrt(kss[ij].get_energy() *
646 kss[kq].get_energy() *
647 kss[ij].get_weight() *
648 kss[kq].get_weight())
650 def __str__(self):
651 str = '<OmegaMatrix> '
652 if hasattr(self, 'eigenvalues'):
653 str += 'dimension ' + ('%d' % len(self.eigenvalues))
654 str += '\neigenvalues: '
655 for ev in self.eigenvalues:
656 str += ' ' + ('%f' % (np.sqrt(ev) * Hartree))
657 return str