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