1import os 2import sys 3from time import time 4 5import numpy as np 6from scipy.special import sici 7from scipy.special.orthogonal import p_roots 8import ase.io.ulm as ulm 9from ase.units import Ha, Bohr 10from ase.utils.timing import timer 11 12import gpaw.mpi as mpi 13from gpaw.blacs import BlacsGrid, Redistributor 14from gpaw.fd_operators import Gradient 15from gpaw.kpt_descriptor import KPointDescriptor 16from gpaw.utilities.blas import gemmdot, axpy 17from gpaw.wavefunctions.pw import PWDescriptor 18from gpaw.xc.rpa import RPACorrelation 19 20 21class FXCCorrelation(RPACorrelation): 22 def __init__(self, 23 calc, 24 xc='RPA', 25 filename=None, 26 skip_gamma=False, 27 qsym=True, 28 nlambda=8, 29 nfrequencies=16, 30 frequency_max=800.0, 31 frequency_scale=2.0, 32 frequencies=None, 33 weights=None, 34 density_cut=1.e-6, 35 world=mpi.world, 36 nblocks=1, 37 unit_cells=None, 38 tag=None, 39 txt=sys.stdout, 40 range_rc=1.0, 41 av_scheme=None, 42 Eg=None): 43 44 RPACorrelation.__init__(self, 45 calc, 46 xc=xc, 47 filename=filename, 48 skip_gamma=skip_gamma, 49 qsym=qsym, 50 nfrequencies=nfrequencies, 51 nlambda=nlambda, 52 frequency_max=frequency_max, 53 frequency_scale=frequency_scale, 54 frequencies=frequencies, 55 weights=weights, 56 world=world, 57 nblocks=nblocks, 58 txt=txt) 59 60 self.l_l, self.weight_l = p_roots(nlambda) 61 self.l_l = (self.l_l + 1.0) * 0.5 62 self.weight_l *= 0.5 63 self.xc = xc 64 self.density_cut = density_cut 65 if unit_cells is None: 66 unit_cells = self.calc.wfs.kd.N_c 67 self.unit_cells = unit_cells 68 self.range_rc = range_rc # Range separation parameter in Bohr 69 self.av_scheme = av_scheme # Either 'density' or 'wavevector' 70 self.Eg = Eg # Band gap in eV 71 72 set_flags(self) 73 74 if tag is None: 75 76 tag = self.calc.atoms.get_chemical_formula(mode='hill') 77 78 if self.av_scheme is not None: 79 80 tag += '_' + self.av_scheme 81 82 self.tag = tag 83 84 @timer('FXC') 85 def calculate(self, ecut, nbands=None): 86 87 if self.xc not in ('RPA', 'range_RPA'): 88 # kernel not required for RPA/range_sep RPA 89 90 if isinstance(ecut, (float, int)): 91 self.ecut_max = ecut 92 else: 93 self.ecut_max = max(ecut) 94 95 # Find the first q vector to calculate kernel for 96 # (density averaging scheme always calculates all q points anyway) 97 98 q_empty = None 99 100 for iq in reversed(range(len(self.ibzq_qc))): 101 102 if not os.path.isfile('fhxc_%s_%s_%s_%s.ulm' % 103 (self.tag, self.xc, self.ecut_max, iq)): 104 q_empty = iq 105 106 if q_empty is not None: 107 108 if self.av_scheme == 'wavevector': 109 110 print('Calculating %s kernel starting from q point %s' % 111 (self.xc, q_empty), 112 file=self.fd) 113 print(file=self.fd) 114 115 if self.linear_kernel: 116 kernel = KernelWave(self.calc, self.xc, self.ibzq_qc, 117 self.fd, None, q_empty, None, 118 self.Eg, self.ecut_max, self.tag, 119 self.timer) 120 121 elif not self.dyn_kernel: 122 kernel = KernelWave(self.calc, self.xc, self.ibzq_qc, 123 self.fd, self.l_l, q_empty, None, 124 self.Eg, self.ecut_max, self.tag, 125 self.timer) 126 127 else: 128 kernel = KernelWave(self.calc, self.xc, self.ibzq_qc, 129 self.fd, self.l_l, q_empty, 130 self.omega_w, self.Eg, 131 self.ecut_max, self.tag, 132 self.timer) 133 134 else: 135 136 kernel = KernelDens(self.calc, self.xc, self.ibzq_qc, 137 self.fd, self.unit_cells, 138 self.density_cut, self.ecut_max, 139 self.tag, self.timer) 140 141 kernel.calculate_fhxc() 142 del kernel 143 144 else: 145 print('%s kernel already calculated' % self.xc, file=self.fd) 146 print(file=self.fd) 147 148 if self.xc in ('range_RPA', 'range_rALDA'): 149 150 shortrange = range_separated(self.calc, self.fd, self.omega_w, 151 self.weight_w, self.l_l, 152 self.weight_l, self.range_rc, self.xc) 153 154 self.shortrange = shortrange.calculate() 155 156 if self.calc.wfs.nspins == 1: 157 spin = False 158 else: 159 spin = True 160 161 e = RPACorrelation.calculate(self, ecut, spin=spin, nbands=nbands) 162 163 return e 164 165 @timer('Chi0(q)') 166 def calculate_q(self, chi0, pd, chi0_swGG, chi0_swxvG, chi0_swvv, m1, m2, 167 cut_G, A2_x): 168 if chi0_swxvG is None: 169 chi0_swxvG = range(2) # Not used 170 chi0_swvv = range(2) # Not used 171 chi0._calculate(pd, 172 chi0_swGG[0], 173 chi0_swxvG[0], 174 chi0_swvv[0], 175 m1, 176 m2, [0], 177 extend_head=False) 178 if len(chi0_swGG) == 2: 179 chi0._calculate(pd, 180 chi0_swGG[1], 181 chi0_swxvG[1], 182 chi0_swvv[1], 183 m1, 184 m2, [1], 185 extend_head=False) 186 print('E_c(q) = ', end='', file=self.fd) 187 188 if self.nblocks > 1: 189 if len(chi0_swGG) == 2: 190 chi0_0wGG = chi0.redistribute(chi0_swGG[0], A2_x) 191 nredist = np.product(chi0_0wGG.shape) 192 chi0.redistribute(chi0_swGG[1], A2_x[nredist:]) 193 chi0_swGG = A2_x[:2 * nredist].reshape((2, ) + chi0_0wGG.shape) 194 chi0_swGG = np.swapaxes(chi0_swGG, 2, 3) 195 else: 196 chi0_0wGG = chi0.redistribute(chi0_swGG[0], A2_x) 197 nredist = np.product(chi0_0wGG.shape) 198 chi0_swGG = A2_x[:1 * nredist].reshape((1, ) + chi0_0wGG.shape) 199 chi0_swGG = np.swapaxes(chi0_swGG, 2, 3) 200 201 if not pd.kd.gamma: 202 e = self.calculate_energy(pd, chi0_swGG, cut_G) 203 print('%.3f eV' % (e * Ha), file=self.fd) 204 self.fd.flush() 205 else: 206 nw = len(self.omega_w) 207 mynw = nw // self.nblocks 208 w1 = self.blockcomm.rank * mynw 209 w2 = w1 + mynw 210 e = 0.0 211 for v in range(3): 212 chi0_swGG[:, :, 0] = chi0_swxvG[:, w1:w2, 0, v] 213 chi0_swGG[:, :, :, 0] = chi0_swxvG[:, w1:w2, 1, v] 214 chi0_swGG[:, :, 0, 0] = chi0_swvv[:, w1:w2, v, v] 215 ev = self.calculate_energy(pd, chi0_swGG, cut_G) 216 e += ev 217 print('%.3f' % (ev * Ha), end='', file=self.fd) 218 if v < 2: 219 print('/', end='', file=self.fd) 220 else: 221 print('eV', file=self.fd) 222 self.fd.flush() 223 e /= 3 224 225 return e 226 227 @timer('Energy') 228 def calculate_energy(self, pd, chi0_swGG, cut_G): 229 """Evaluate correlation energy from chi0 and the kernel fhxc""" 230 231 ibzq2_q = [ 232 np.dot(self.ibzq_qc[i] - pd.kd.bzk_kc[0], 233 self.ibzq_qc[i] - pd.kd.bzk_kc[0]) 234 for i in range(len(self.ibzq_qc)) 235 ] 236 237 qi = np.argsort(ibzq2_q)[0] 238 239 G_G = pd.G2_qG[0]**0.5 # |G+q| 240 241 if cut_G is not None: 242 G_G = G_G[cut_G] 243 244 nG = len(G_G) 245 ns = len(chi0_swGG) 246 247 # There are three options to calculate the 248 # energy depending on kernel and/or averaging scheme. 249 250 # Option (1) - Spin-polarized form of kernel exists 251 # e.g. rALDA, rAPBE. 252 # Then, solve block diagonal form of Dyson 253 # equation (dimensions (ns*nG) * (ns*nG)) 254 # (note this does not necessarily mean that 255 # the calculation is spin-polarized!) 256 257 if self.spin_kernel: 258 r = ulm.open('fhxc_%s_%s_%s_%s.ulm' % 259 (self.tag, self.xc, self.ecut_max, qi)) 260 fv = r.fhxc_sGsG 261 262 if cut_G is not None: 263 cut_sG = np.tile(cut_G, ns) 264 cut_sG[len(cut_G):] += len(fv) // ns 265 fv = fv.take(cut_sG, 0).take(cut_sG, 1) 266 267 # the spin-polarized kernel constructed from wavevector average 268 # is already multiplied by |q+G| |q+G'|/4pi, and doesn't require 269 # special treatment of the head and wings. However not true for 270 # density average: 271 272 if self.av_scheme == 'density': 273 for s1 in range(ns): 274 for s2 in range(ns): 275 m1 = s1 * nG 276 n1 = (s1 + 1) * nG 277 m2 = s2 * nG 278 n2 = (s2 + 1) * nG 279 fv[m1:n1, 280 m2:n2] *= (G_G * G_G[:, np.newaxis] / (4 * np.pi)) 281 282 if np.prod(self.unit_cells) > 1 and pd.kd.gamma: 283 m1 = s1 * nG 284 n1 = (s1 + 1) * nG 285 m2 = s2 * nG 286 n2 = (s2 + 1) * nG 287 fv[m1, m2:n2] = 0.0 288 fv[m1:n1, m2] = 0.0 289 fv[m1, m2] = 1.0 290 291 if pd.kd.gamma: 292 G_G[0] = 1.0 293 294 e_w = [] 295 296 # Loop over frequencies 297 for chi0_sGG in np.swapaxes(chi0_swGG, 0, 1): 298 if cut_G is not None: 299 chi0_sGG = chi0_sGG.take(cut_G, 1).take(cut_G, 2) 300 chi0v = np.zeros((ns * nG, ns * nG), dtype=complex) 301 for s in range(ns): 302 m = s * nG 303 n = (s + 1) * nG 304 chi0v[m:n, m:n] = chi0_sGG[s] / G_G / G_G[:, np.newaxis] 305 chi0v *= 4 * np.pi 306 307 del chi0_sGG 308 309 e = 0.0 310 311 for l, weight in zip(self.l_l, self.weight_l): 312 chiv = np.linalg.solve( 313 np.eye(nG * ns) - l * np.dot(chi0v, fv), 314 chi0v).real # this is SO slow 315 for s1 in range(ns): 316 for s2 in range(ns): 317 m1 = s1 * nG 318 n1 = (s1 + 1) * nG 319 m2 = s2 * nG 320 n2 = (s2 + 1) * nG 321 chiv_s1s2 = chiv[m1:n1, m2:n2] 322 e -= np.trace(chiv_s1s2) * weight 323 324 e += np.trace(chi0v.real) 325 326 e_w.append(e) 327 328 else: 329 # Or, if kernel does not have a spin polarized form, 330 # 331 # Option (2) kernel does not scale linearly with lambda, 332 # so we solve nG*nG Dyson equation at each value 333 # of l. Requires kernel to be constructed 334 # at individual values of lambda 335 # 336 # Option (3) Divide correlation energy into 337 # long range part which can be integrated 338 # analytically w.r.t. lambda, and a short 339 # range part which again requires 340 # solving Dyson equation (hence no speedup, 341 # but the maths looks nice and 342 # fits with range-separated RPA) 343 # 344 # 345 # Construct/read kernels 346 347 if self.xc == 'RPA': 348 349 fv = np.eye(nG) 350 351 elif self.xc == 'range_RPA': 352 353 fv = np.exp(-0.25 * (G_G * self.range_rc)**2.0) 354 355 elif self.linear_kernel: 356 r = ulm.open('fhxc_%s_%s_%s_%s.ulm' % 357 (self.tag, self.xc, self.ecut_max, qi)) 358 fv = r.fhxc_sGsG 359 360 if cut_G is not None: 361 fv = fv.take(cut_G, 0).take(cut_G, 1) 362 363 elif not self.dyn_kernel: 364 # static kernel which does not scale with lambda 365 366 r = ulm.open('fhxc_%s_%s_%s_%s.ulm' % 367 (self.tag, self.xc, self.ecut_max, qi)) 368 fv = r.fhxc_lGG 369 370 if cut_G is not None: 371 fv = fv.take(cut_G, 1).take(cut_G, 2) 372 373 else: # dynamical kernel 374 r = ulm.open('fhxc_%s_%s_%s_%s.ulm' % 375 (self.tag, self.xc, self.ecut_max, qi)) 376 fv = r.fhxc_lwGG 377 378 if cut_G is not None: 379 fv = fv.take(cut_G, 2).take(cut_G, 3) 380 381 if pd.kd.gamma: 382 G_G[0] = 1.0 383 384 # Loop over frequencies; since the kernel has no spin, 385 # we work with spin-summed response function 386 e_w = [] 387 iw = 0 388 389 for chi0_sGG in np.swapaxes(chi0_swGG, 0, 1): 390 if cut_G is not None: 391 chi0_sGG = chi0_sGG.take(cut_G, 1).take(cut_G, 2) 392 chi0v = np.zeros((nG, nG), dtype=complex) 393 for s in range(ns): 394 chi0v += chi0_sGG[s] / G_G / G_G[:, np.newaxis] 395 chi0v *= 4 * np.pi 396 del chi0_sGG 397 398 e = 0.0 399 400 if not self.linear_kernel: 401 402 il = 0 403 for l, weight in zip(self.l_l, self.weight_l): 404 405 if not self.dyn_kernel: 406 chiv = np.linalg.solve( 407 np.eye(nG) - np.dot(chi0v, fv[il]), chi0v).real 408 else: 409 chiv = np.linalg.solve( 410 np.eye(nG) - np.dot(chi0v, fv[il][iw]), 411 chi0v).real 412 e -= np.trace(chiv) * weight 413 il += 1 414 415 e += np.trace(chi0v.real) 416 417 e_w.append(e) 418 419 iw += 1 420 421 else: 422 423 # Coupling constant integration 424 # for long-range part 425 # Do this analytically, except for the RPA 426 # simply since the analytical method is already 427 # implemented in rpa.py 428 if self.xc == 'range_RPA': 429 # way faster than np.dot for diagonal kernels 430 e_GG = np.eye(nG) - chi0v * fv 431 elif self.xc != 'RPA': 432 e_GG = np.eye(nG) - np.dot(chi0v, fv) 433 434 if self.xc == 'RPA': 435 # numerical RPA 436 elong = 0.0 437 for l, weight in zip(self.l_l, self.weight_l): 438 439 chiv = np.linalg.solve( 440 np.eye(nG) - l * np.dot(chi0v, fv), chi0v).real 441 442 elong -= np.trace(chiv) * weight 443 444 elong += np.trace(chi0v.real) 445 446 else: 447 # analytic everything else 448 elong = (np.log(np.linalg.det(e_GG)) + nG - 449 np.trace(e_GG)).real 450 # Numerical integration for short-range part 451 eshort = 0.0 452 if self.xc not in ('RPA', 'range_RPA', 'range_rALDA'): 453 fxcv = fv - np.eye(nG) # Subtract Hartree contribution 454 455 for l, weight in zip(self.l_l, self.weight_l): 456 457 chiv = np.linalg.solve( 458 np.eye(nG) - l * np.dot(chi0v, fv), chi0v) 459 eshort += (np.trace(np.dot(chiv, fxcv)).real * 460 weight) 461 462 eshort -= np.trace(np.dot(chi0v, fxcv)).real 463 464 elif self.xc in ('range_RPA', 'range_rALDA'): 465 eshort = (2 * np.pi * self.shortrange / 466 np.sum(self.weight_w)) 467 468 e = eshort + elong 469 e_w.append(e) 470 471 E_w = np.zeros_like(self.omega_w) 472 self.blockcomm.all_gather(np.array(e_w), E_w) 473 energy = np.dot(E_w, self.weight_w) / (2 * np.pi) 474 return energy 475 476 477class KernelWave: 478 def __init__(self, calc, xc, ibzq_qc, fd, l_l, q_empty, omega_w, Eg, ecut, 479 tag, timer): 480 481 self.calc = calc 482 self.gd = calc.density.gd 483 self.xc = xc 484 self.ibzq_qc = ibzq_qc 485 self.fd = fd 486 self.l_l = l_l 487 self.ns = calc.wfs.nspins 488 self.q_empty = q_empty 489 self.omega_w = omega_w 490 self.Eg = Eg 491 self.ecut = ecut 492 self.tag = tag 493 self.timer = timer 494 495 if l_l is None: # -> Kernel is linear in coupling strength 496 self.l_l = [1.0] 497 498 # Density grid 499 self.n_g = calc.get_all_electron_density(gridrefinement=2) 500 self.n_g = self.n_g.flatten() * Bohr**3 501 # Density now in units electrons/cubic Bohr 502 # For atoms with large vacuum regions 503 # this apparently can take negative values! 504 mindens = np.amin(self.n_g) 505 506 if mindens < 0: 507 print('Negative densities found! (magnitude %s)' % np.abs(mindens), 508 file=self.fd) 509 print('These will be reset to 1E-12 elec/bohr^3)', file=self.fd) 510 self.n_g[np.where(self.n_g < 0.0)] = 1.0E-12 511 512 r_g = self.gd.refine().get_grid_point_coordinates() # already in Bohr 513 self.x_g = 1.0 * r_g[0].flatten() 514 self.y_g = 1.0 * r_g[1].flatten() 515 self.z_g = 1.0 * r_g[2].flatten() 516 self.gridsize = len(self.x_g) 517 assert len(self.n_g) == self.gridsize 518 519 if self.omega_w is not None: 520 print('Calculating dynamical kernel at %s frequencies' % 521 len(self.omega_w), 522 file=self.fd) 523 524 if self.Eg is not None: 525 print('Band gap of %s eV used to evaluate kernel' % (self.Eg * Ha), 526 file=self.fd) 527 528 # Enhancement factor for GGA 529 if self.xc == 'rAPBE' or self.xc == 'rAPBEns': 530 nf_g = calc.get_all_electron_density(gridrefinement=4) 531 nf_g *= Bohr**3 532 gdf = self.gd.refine().refine() 533 grad_v = [Gradient(gdf, v, n=1).apply for v in range(3)] 534 gradnf_vg = gdf.empty(3) 535 536 for v in range(3): 537 grad_v[v](nf_g, gradnf_vg[v]) 538 539 self.s2_g = np.sqrt(np.sum(gradnf_vg[:, ::2, ::2, ::2]**2.0, 540 0)).flatten() # |\nabla\rho| 541 self.s2_g *= 1.0 / (2.0 * (3.0 * np.pi**2.0)**(1.0 / 3.0) * 542 self.n_g**(4.0 / 3.0)) 543 # |\nabla\rho|/(2kF\rho) = s 544 self.s2_g = self.s2_g**2 # s^2 545 assert len(self.n_g) == len(self.s2_g) 546 547 # Now we find all the regions where the 548 # APBE kernel wants to be positive, and hack s to = 0, 549 # so that we are really using the ALDA kernel 550 # at these points 551 apbe_g = self.get_PBE_fxc(self.n_g, self.s2_g) 552 poskern_ind = np.where(apbe_g >= 0.0) 553 if len(poskern_ind[0]) > 0: 554 print('The APBE kernel takes positive values at ' + 555 '%s grid points out of a total of %s (%3.2f%%).' % 556 (len(poskern_ind[0]), self.gridsize, 557 100.0 * len(poskern_ind[0]) / self.gridsize), 558 file=self.fd) 559 print('The ALDA kernel will be used at these points', 560 file=self.fd) 561 self.s2_g[poskern_ind] = 0.0 562 563 def calculate_fhxc(self): 564 565 print('Calculating %s kernel at %d eV cutoff' % (self.xc, self.ecut), 566 file=self.fd) 567 self.fd.flush() 568 569 for iq, q_c in enumerate(self.ibzq_qc): 570 571 if iq < self.q_empty: # don't recalculate q vectors 572 continue 573 574 thisqd = KPointDescriptor([q_c]) 575 pd = PWDescriptor(self.ecut / Ha, self.gd, complex, thisqd) 576 577 nG = pd.ngmax 578 G_G = pd.G2_qG[0]**0.5 # |G+q| 579 Gv_G = pd.get_reciprocal_vectors(q=0, add_q=False) 580 # G as a vector (note we are at a specific q point here so set q=0) 581 582 # Distribute G vectors among processors 583 # Later we calculate for iG' > iG, 584 # so stagger allocation in order to balance load 585 local_Gvec_grid_size = nG // mpi.world.size 586 my_Gints = (mpi.world.rank + np.arange( 587 0, local_Gvec_grid_size * mpi.world.size, mpi.world.size)) 588 589 if (mpi.world.rank + (local_Gvec_grid_size) * mpi.world.size) < nG: 590 my_Gints = np.append( 591 my_Gints, 592 [mpi.world.rank + local_Gvec_grid_size * mpi.world.size]) 593 594 my_Gv_G = Gv_G[my_Gints] 595 596 if (self.ns == 2) and (self.xc == 'rALDA' or self.xc == 'rAPBE'): 597 598 assert len(self.l_l) == 1 599 600 # Form spin-dependent kernel according to 601 # PRB 88, 115131 (2013) equation 20 602 # (note typo, should be \tilde{f^rALDA}) 603 # spincorr is just the ALDA exchange kernel 604 # with a step function (\equiv \tilde{f^rALDA}) 605 # fHxc^{up up} = fHxc^{down down} = fv_nospin + fv_spincorr 606 # fHxc^{up down} = fHxc^{down up} = fv_nospin - fv_spincorr 607 608 calc_spincorr = True 609 fv_spincorr = np.zeros((nG, nG), dtype=complex) 610 611 else: 612 613 calc_spincorr = False 614 615 if self.omega_w is None: 616 fv_nospin = np.zeros((len(self.l_l), nG, nG), dtype=complex) 617 else: 618 fv_nospin = np.zeros( 619 (len(self.l_l), len(self.omega_w), nG, nG), dtype=complex) 620 621 for il, l in enumerate(self.l_l): # loop over coupling constant 622 623 for iG, Gv in zip(my_Gints, my_Gv_G): # loop over G vecs 624 625 # For all kernels except JGM we 626 # treat head and wings analytically 627 if G_G[iG] > 1.0E-5 or self.xc == 'JGMs': 628 629 # Symmetrised |q+G||q+G'|, where iG' >= iG 630 mod_Gpq = np.sqrt(G_G[iG] * G_G[iG:]) 631 632 # Phase factor \vec{G}-\vec{G'} 633 deltaGv = Gv - Gv_G[iG:] 634 635 if (self.xc in ('rALDA', 'range_rALDA', 'rALDAns')): 636 637 # rALDA trick: the Hartree-XC kernel is exactly 638 # zero for densities below rho_min = 639 # min_Gpq^3/(24*pi^2), 640 # so we don't need to include these contributions 641 # in the Fourier transform 642 643 min_Gpq = np.amin(mod_Gpq) 644 rho_min = min_Gpq**3.0 / (24.0 * np.pi**2.0) 645 small_ind = np.where(self.n_g >= rho_min) 646 647 elif (self.xc == 'rAPBE' or self.xc == 'rAPBEns'): 648 649 # rAPBE trick: the Hartree-XC kernel 650 # is exactly zero at grid points where 651 # min_Gpq > cutoff wavevector 652 653 min_Gpq = np.amin(mod_Gpq) 654 small_ind = np.where(min_Gpq <= np.sqrt( 655 -4.0 * np.pi / 656 self.get_PBE_fxc(self.n_g, self.s2_g))) 657 658 else: 659 660 small_ind = np.arange(self.gridsize) 661 662 phase_Gpq = np.exp( 663 -1.0j * 664 (deltaGv[:, 0, np.newaxis] * self.x_g[small_ind] + 665 deltaGv[:, 1, np.newaxis] * self.y_g[small_ind] + 666 deltaGv[:, 2, np.newaxis] * self.z_g[small_ind])) 667 if self.omega_w is None: 668 fv_nospin[il, iG, iG:] = self.get_scaled_fHxc_q( 669 q=mod_Gpq, 670 sel_points=small_ind, 671 Gphase=phase_Gpq, 672 l=l, 673 spincorr=False, 674 w=None) 675 else: 676 for iw, omega in enumerate(self.omega_w): 677 fv_nospin[il, iw, iG, iG:] = \ 678 self.get_scaled_fHxc_q( 679 q=mod_Gpq, 680 sel_points=small_ind, 681 Gphase=phase_Gpq, 682 l=l, 683 spincorr=False, 684 w=omega) 685 686 if calc_spincorr: 687 fv_spincorr[iG, iG:] = self.get_scaled_fHxc_q( 688 q=mod_Gpq, 689 sel_points=small_ind, 690 Gphase=phase_Gpq, 691 l=1.0, 692 spincorr=True, 693 w=None) 694 695 else: 696 # head and wings of q=0 are dominated by 697 # 1/q^2 divergence of scaled Coulomb interaction 698 699 assert iG == 0 700 701 if self.omega_w is None: 702 fv_nospin[il, 0, 0] = l 703 fv_nospin[il, 0, 1:] = 0.0 704 else: 705 fv_nospin[il, :, 0, 0] = l 706 fv_nospin[il, :, 0, 1:] = 0.0 707 708 if calc_spincorr: 709 fv_spincorr[0, :] = 0.0 710 711 # End loop over G vectors 712 713 mpi.world.sum(fv_nospin[il]) 714 715 if self.omega_w is None: 716 # We've only got half the matrix here, 717 # so add the hermitian conjugate: 718 fv_nospin[il] += np.conj(fv_nospin[il].T) 719 # but now the diagonal's been doubled, 720 # so we multiply these elements by 0.5 721 fv_nospin[il][np.diag_indices(nG)] *= 0.5 722 723 else: # same procedure for dynamical kernels 724 for iw in range(len(self.omega_w)): 725 fv_nospin[il][iw] += np.conj(fv_nospin[il][iw].T) 726 fv_nospin[il][iw][np.diag_indices(nG)] *= 0.5 727 728 # End of loop over coupling constant 729 730 if calc_spincorr: 731 mpi.world.sum(fv_spincorr) 732 fv_spincorr += np.conj(fv_spincorr.T) 733 fv_spincorr[np.diag_indices(nG)] *= 0.5 734 735 # Write to disk 736 if mpi.rank == 0: 737 738 w = ulm.open( 739 'fhxc_%s_%s_%s_%s.ulm' % 740 (self.tag, self.xc, self.ecut, iq), 'w') 741 742 if calc_spincorr: 743 # Form the block matrix kernel 744 fv_full = np.empty((2 * nG, 2 * nG), dtype=complex) 745 fv_full[:nG, :nG] = fv_nospin[0] + fv_spincorr 746 fv_full[:nG, nG:] = fv_nospin[0] - fv_spincorr 747 fv_full[nG:, :nG] = fv_nospin[0] - fv_spincorr 748 fv_full[nG:, nG:] = fv_nospin[0] + fv_spincorr 749 w.write(fhxc_sGsG=fv_full) 750 751 elif len(self.l_l) == 1: 752 w.write(fhxc_sGsG=fv_nospin[0]) 753 754 elif self.omega_w is None: 755 w.write(fhxc_lGG=fv_nospin) 756 757 else: 758 w.write(fhxc_lwGG=fv_nospin) 759 w.close() 760 761 print('q point %s complete' % iq, file=self.fd) 762 self.fd.flush() 763 764 mpi.world.barrier() 765 766 def get_scaled_fHxc_q(self, q, sel_points, Gphase, l, spincorr, w): 767 # Given a coupling constant l, construct the Hartree-XC 768 # kernel in q space a la Lein, Gross and Perdew, 769 # Phys. Rev. B 61, 13431 (2000): 770 # 771 # f_{Hxc}^\lambda(q,\omega,r_s) = \frac{4\pi \lambda }{q^2} + 772 # \frac{1}{\lambda} f_{xc}(q/\lambda,\omega/\lambda^2,\lambda r_s) 773 # 774 # divided by the unscaled Coulomb interaction!! 775 # 776 # i.e. this subroutine returns f_{Hxc}^\lambda(q,\omega,r_s) 777 # * \frac{q^2}{4\pi} 778 # = \lambda * [\frac{(q/lambda)^2}{4\pi} 779 # f_{Hxc}(q/\lambda,\omega/\lambda^2,\lambda r_s)] 780 # = \lambda * [1/scaled_coulomb * fHxc computed with scaled quantities] 781 782 # Apply scaling 783 rho = self.n_g[sel_points] 784 785 # GGA enhancement factor s is lambda independent, 786 # but we might want to truncate it 787 if self.xc == 'rAPBE' or self.xc == 'rAPBEns': 788 s2_g = self.s2_g[sel_points] 789 else: 790 s2_g = None 791 792 scaled_q = q / l 793 scaled_rho = rho / l**3.0 794 scaled_rs = (3.0 / (4.0 * np.pi * scaled_rho))**(1.0 / 3.0 795 ) # Wigner radius 796 if w is not None: 797 scaled_w = w / (l**2.0) 798 else: 799 scaled_w = None 800 801 if self.Eg is not None: 802 scaled_Eg = self.Eg / (l**1.5) 803 else: 804 scaled_Eg = None 805 806 if not spincorr: 807 scaled_kernel = l * self.get_fHxc_q(scaled_rs, scaled_q, Gphase, 808 s2_g, scaled_w, scaled_Eg) 809 else: 810 scaled_kernel = l * self.get_spinfHxc_q(scaled_rs, scaled_q, 811 Gphase, s2_g) 812 813 return (scaled_kernel) 814 815 def get_fHxc_q(self, rs, q, Gphase, s2_g, w, scaled_Eg): 816 # Construct fHxc(q,G,:), divided by scaled Coulomb interaction 817 818 qF = (9.0 * np.pi / 4.0)**(1.0 / 3.0) * 1.0 / rs 819 820 if self.xc in ('rALDA', 'rALDAns', 'range_rALDA'): 821 # rALDA (exchange only) kernel 822 # Olsen and Thygesen, Phys. Rev. B 88, 115131 (2013) 823 # ALDA up to 2*qF, -vc for q >2qF (such that fHxc vanishes) 824 825 rxalda_A = 0.25 826 rxalda_qcut = qF * np.sqrt(1.0 / rxalda_A) 827 828 # construct fHxc(k,r) 829 fHxc_Gr = (0.5 + 0.0j) * ( 830 (1.0 + np.sign(rxalda_qcut - q[:, np.newaxis])) * 831 (1.0 + (-1.0) * rxalda_A * (q[:, np.newaxis] / qF)**2.0)) 832 833 elif self.xc == 'rALDAc': 834 # rALDA exchange+correlation kernel 835 # Olsen and Thygesen, Phys. Rev. B 88, 115131 (2013) 836 # Cutoff wavevector ensures kernel is continuous 837 838 rxcalda_A = self.get_heg_A(rs) 839 rxcalda_qcut = qF * np.sqrt(1.0 / rxcalda_A) 840 841 # construct fHxc(k,r) 842 fHxc_Gr = (0.5 + 0.0j) * ( 843 (1.0 + np.sign(rxcalda_qcut - q[:, np.newaxis])) * 844 (1.0 + (-1.0) * rxcalda_A * (q[:, np.newaxis] / qF)**2.0)) 845 846 elif self.xc == 'rAPBE' or self.xc == 'rAPBEns': 847 848 # Olsen and Thygesen, Phys. Rev. Lett. 112, 203001 (2014) 849 # Exchange only part of the PBE XC kernel, neglecting the terms 850 # arising from the variation of the density gradient 851 # i.e. second functional derivative 852 # d2/drho^2 -> \partial^2/\partial rho^2 at fixed \nabla \rho 853 854 fxc_PBE = self.get_PBE_fxc(pbe_rho=3.0 / (4.0 * np.pi * rs**3.0), 855 pbe_s2_g=s2_g) 856 rxapbe_qcut = np.sqrt(-4.0 * np.pi / fxc_PBE) 857 858 fHxc_Gr = (0.5 + 0.0j) * ( 859 (1.0 + np.sign(rxapbe_qcut - q[:, np.newaxis])) * 860 (1.0 + fxc_PBE / (4.0 * np.pi) * (q[:, np.newaxis])**2.0)) 861 862 elif self.xc == 'CP': 863 # Constantin & Pitarke, Phys. Rev. B 75, 245127 (2007) 864 # equation 4, omega = 0 865 # and equation 9 [note that fxc(q=0) = fxc^ALDA = -4piA/qf^2] 866 # The simplest, static error-function kernel. 867 # Produces correct q = 0 limit, but not q->oo 868 869 cp_kappa = self.get_heg_A(rs) / (qF**2.0) 870 fHxc_Gr = (1.0 + 0.0j) * np.exp(-cp_kappa * q[:, np.newaxis]**2.0) 871 872 elif self.xc == 'CP_dyn': 873 # CP kernel with frequency dependence 874 875 cp_A = self.get_heg_A(rs) 876 cp_D = self.get_heg_D(rs) 877 cp_c = cp_D / cp_A 878 cp_a = 6.0 * np.sqrt(cp_c) 879 880 cp_kappa_w = cp_A / (qF**2.0) 881 cp_kappa_w *= (1.0 + cp_a * w + cp_c * w**2.0) / (1.0 + w**2.0) 882 883 fHxc_Gr = (1.0 + 0.0j) * np.exp( 884 -cp_kappa_w * q[:, np.newaxis]**2.0) 885 886 elif self.xc == 'CDOP' or self.xc == 'CDOPs': 887 # Corradini, Del Sole, Onida and Palummo (CDOP), 888 # Phys. Rev. B 57, 14569 (1998) 889 # Reproduces exact q=0 and q -> oo limits 890 891 cdop_Q = q[:, np.newaxis] / qF 892 cdop_A = self.get_heg_A(rs) 893 cdop_B = self.get_heg_B(rs) 894 if self.xc == 'CDOPs': 895 cdop_C = np.zeros(np.shape(cdop_B)) 896 else: 897 cdop_C = self.get_heg_C(rs) 898 cdop_g = cdop_B / (cdop_A - cdop_C) 899 cdop_alpha = 1.5 / (rs**0.25) * cdop_A / (cdop_B * cdop_g) 900 cdop_beta = 1.2 / (cdop_B * cdop_g) 901 902 fHxc_Gr = -cdop_C * cdop_Q**2.0 * (1.0 + 0.0j) 903 fHxc_Gr += -cdop_B * cdop_Q**2.0 * 1.0 / (cdop_g + cdop_Q**2.0) 904 fHxc_Gr += -(cdop_alpha * cdop_Q**4.0 * 905 np.exp(-1.0 * cdop_beta * cdop_Q**2.0)) 906 907 # Hartree 908 fHxc_Gr += 1.0 909 910 elif self.xc == 'JGMs': 911 # Trevisanutto et al., 912 # Phys. Rev. B 87, 205143 (2013) 913 # equation 4, 914 # so-called "jellium with gap" model, but simplified 915 # so that it is similar to CP. 916 917 jgm_kappa = self.get_heg_A(rs) / (qF**2.0) 918 jgm_n = 3.0 / (4.0 * np.pi * rs**3.0) 919 fHxc_Gr = (1.0 + 920 0.0j) * (np.exp(-jgm_kappa * q[:, np.newaxis]**2.0) * 921 np.exp(-scaled_Eg**2.0 / 922 (4.0 * np.pi * jgm_n))) 923 924 # Integrate over r with phase 925 fHxc_Gr *= Gphase 926 fHxc_GG = np.sum(fHxc_Gr, 1) / self.gridsize 927 return (fHxc_GG) 928 929 def get_spinfHxc_q(self, rs, q, Gphase, s2_g): 930 931 qF = (9.0 * np.pi / 4.0)**(1.0 / 3.0) * 1.0 / rs 932 933 if self.xc == 'rALDA': 934 935 rxalda_A = 0.25 936 rxalda_qcut = qF * np.sqrt(1.0 / rxalda_A) 937 938 fspinHxc_Gr = ((0.5 + 0.0j) * 939 (1.0 + np.sign(rxalda_qcut - q[:, np.newaxis])) * 940 (-1.0) * rxalda_A * (q[:, np.newaxis] / qF)**2.0) 941 942 elif self.xc == 'rAPBE': 943 944 fxc_PBE = self.get_PBE_fxc(pbe_rho=(3.0 / (4.0 * np.pi * rs**3.0)), 945 pbe_s2_g=s2_g) 946 rxapbe_qcut = np.sqrt(-4.0 * np.pi / fxc_PBE) 947 948 fspinHxc_Gr = ((0.5 + 0.0j) * 949 (1.0 + np.sign(rxapbe_qcut - q[:, np.newaxis])) * 950 fxc_PBE / (4.0 * np.pi) * q[:, np.newaxis]**2.0) 951 952 fspinHxc_Gr *= Gphase 953 fspinHxc_GG = np.sum(fspinHxc_Gr, 1) / self.gridsize 954 return (fspinHxc_GG) 955 956 def get_PBE_fxc(self, pbe_rho, pbe_s2_g): 957 958 pbe_kappa = 0.804 959 pbe_mu = 0.2195149727645171 960 961 pbe_denom_g = 1.0 + pbe_mu * pbe_s2_g / pbe_kappa 962 963 F_g = 1.0 + pbe_kappa - pbe_kappa / pbe_denom_g 964 Fn_g = -8.0 / 3.0 * pbe_mu * pbe_s2_g / pbe_rho / pbe_denom_g**2.0 965 Fnn_g = (-11.0 / 3.0 / pbe_rho * Fn_g - 966 2.0 / pbe_kappa * Fn_g**2.0 * pbe_denom_g) 967 968 e_g = -3.0 / 4.0 * (3.0 / np.pi)**(1.0 / 3.0) * pbe_rho**(4.0 / 3.0) 969 v_g = 4.0 / 3.0 * e_g / pbe_rho 970 f_g = 1.0 / 3.0 * v_g / pbe_rho 971 972 pbe_f_g = f_g * F_g + 2.0 * v_g * Fn_g + e_g * Fnn_g 973 974 return (pbe_f_g) 975 976 def get_heg_A(self, rs): 977 # Returns the A coefficient, where the 978 # q ->0 limiting value of static fxc 979 # of the HEG = -\frac{4\pi A }{q_F^2} = f_{xc}^{ALDA}. 980 # We need correlation energy per electron and first and second derivs 981 # w.r.t. rs 982 # See for instance Moroni, Ceperley and Senatore, 983 # Phys. Rev. Lett. 75, 689 (1995) 984 # (and also Kohn and Sham, Phys. Rev. 140, A1133 (1965) equation 2.7) 985 986 # Exchange contribution 987 heg_A = 0.25 988 989 # Correlation contribution 990 A_ec, A_dec, A_d2ec = self.get_pw_lda(rs) 991 heg_A += (1.0 / 27.0 * rs**2.0 * (9.0 * np.pi / 4.0)**(2.0 / 3.0) * 992 (2 * A_dec - rs * A_d2ec)) 993 994 return (heg_A) 995 996 def get_heg_B(self, rs): 997 # Returns the B coefficient, where the 998 # q -> oo limiting behaviour of static fxc 999 # of the HEG is -\frac{4\pi B}{q^2} - \frac{4\pi C}{q_F^2}. 1000 # Use the parametrisation of Moroni, Ceperley and Senatore, 1001 # Phys. Rev. Lett. 75, 689 (1995) 1002 1003 mcs_xs = np.sqrt(rs) 1004 1005 mcs_a = (1.0, 2.15, 0.0, 0.435) 1006 mcs_b = (3.0, 1.57, 0.0, 0.409) 1007 1008 mcs_num = 0 1009 1010 for mcs_j, mcs_coeff in enumerate(mcs_a): 1011 mcs_num += mcs_coeff * mcs_xs**mcs_j 1012 1013 mcs_denom = 0 1014 1015 for mcs_j, mcs_coeff in enumerate(mcs_b): 1016 mcs_denom += mcs_coeff * mcs_xs**mcs_j 1017 1018 heg_B = mcs_num / mcs_denom 1019 return (heg_B) 1020 1021 def get_heg_C(self, rs): 1022 # Returns the C coefficient, where the 1023 # q -> oo limiting behaviour of static fxc 1024 # of the HEG is -\frac{4\pi B}{q^2} - \frac{4\pi C}{q_F^2}. 1025 # Again see Moroni, Ceperley and Senatore, 1026 # Phys. Rev. Lett. 75, 689 (1995) 1027 1028 C_ec, C_dec, Cd2ec = self.get_pw_lda(rs) 1029 1030 heg_C = ((-1.0) * np.pi**(2.0 / 3.0) * (1.0 / 18.0)**(1.0 / 3.0) * 1031 (rs * C_ec + rs**2.0 * C_dec)) 1032 1033 return (heg_C) 1034 1035 def get_heg_D(self, rs): 1036 # Returns a 'D' coefficient, where the 1037 # q->0 omega -> oo limiting behaviour 1038 # of the frequency dependent fxc is -\frac{4\pi D}{q_F^2} 1039 # see Constantin & Pitarke Phys. Rev. B 75, 245127 (2007) equation 7 1040 1041 D_ec, D_dec, D_d2ec = self.get_pw_lda(rs) 1042 1043 # Exchange contribution 1044 heg_D = 0.15 1045 # Correlation contribution 1046 heg_D += ((9.0 * np.pi / 4.0)**(2.0 / 3.0) * rs / 3.0 * 1047 (22.0 / 15.0 * D_ec + 26.0 / 15.0 * rs * D_dec)) 1048 return (heg_D) 1049 1050 def get_pw_lda(self, rs): 1051 # Returns LDA correlation energy and its first and second 1052 # derivatives with respect to rs according to the parametrisation 1053 # of Perdew and Wang, Phys. Rev. B 45, 13244 (1992) 1054 1055 pw_A = 0.031091 1056 pw_alp = 0.21370 1057 pw_beta = (7.5957, 3.5876, 1.6382, 0.49294) 1058 1059 pw_logdenom = 2.0 * pw_A * ( 1060 pw_beta[0] * rs**0.5 + pw_beta[1] * rs**1.0 + 1061 pw_beta[2] * rs**1.5 + pw_beta[3] * rs**2.0) 1062 1063 pw_dlogdenom = 2.0 * pw_A * (0.5 * pw_beta[0] * rs**(-0.5) + 1064 1.0 * pw_beta[1] + 1065 1.5 * pw_beta[2] * rs**0.5 + 1066 2.0 * pw_beta[3] * rs) 1067 1068 pw_d2logdenom = 2.0 * pw_A * (-0.25 * pw_beta[0] * rs**(-1.5) + 1069 0.75 * pw_beta[2] * rs**(-0.5) + 1070 2.0 * pw_beta[3]) 1071 1072 pw_logarg = 1.0 + 1.0 / pw_logdenom 1073 pw_dlogarg = (-1.0) / (pw_logdenom**2.0) * pw_dlogdenom 1074 pw_d2logarg = 2.0 / (pw_logdenom**3.0) * (pw_dlogdenom**2.0) 1075 pw_d2logarg += (-1.0) / (pw_logdenom**2.0) * pw_d2logdenom 1076 1077 # pw_ec = the correlation energy (per electron) 1078 pw_ec = -2.0 * pw_A * (1 + pw_alp * rs) * np.log(pw_logarg) 1079 1080 # pw_dec = first derivative 1081 1082 pw_dec = -2.0 * pw_A * (1 + pw_alp * rs) * pw_dlogarg / pw_logarg 1083 pw_dec += (-2.0 * pw_A * pw_alp) * np.log(pw_logarg) 1084 1085 # pw_d2ec = second derivative 1086 1087 pw_d2ec = (-2.0) * pw_A * pw_alp * pw_dlogarg / pw_logarg 1088 pw_d2ec += (-2.0) * pw_A * ((1 + pw_alp * rs) * 1089 (pw_d2logarg / pw_logarg - 1090 (pw_dlogarg**2.0) / (pw_logarg**2.0))) 1091 pw_d2ec += (-2.0 * pw_A * pw_alp) * pw_dlogarg / pw_logarg 1092 1093 return ((pw_ec, pw_dec, pw_d2ec)) 1094 1095 1096class range_separated: 1097 def __init__(self, calc, fd, frequencies, freqweights, l_l, lweights, 1098 range_rc, xc): 1099 1100 self.calc = calc 1101 self.fd = fd 1102 self.frequencies = frequencies 1103 self.freqweights = freqweights 1104 self.l_l = l_l 1105 self.lweights = lweights 1106 self.range_rc = range_rc 1107 self.xc = xc 1108 1109 self.cutoff_rs = 36.278317 1110 1111 if self.xc == 'range_RPA': 1112 print( 1113 'Using range-separated RPA approach, with parameter %s Bohr' % 1114 (self.range_rc), 1115 file=self.fd) 1116 1117 nval_g = calc.get_all_electron_density( 1118 gridrefinement=4, skip_core=True).flatten() * Bohr**3 1119 self.dv = self.calc.density.gd.dv / 64.0 # 64 = gridrefinement^3 1120 1121 density_cut = 3.0 / (4.0 * np.pi * self.cutoff_rs**3.0) 1122 if (nval_g < 0.0).any(): 1123 print('Warning, negative densities found! (Magnitude %s)' % 1124 np.abs(np.amin(nval_g)), 1125 file=self.fd) 1126 print('These will be ignored', file=self.fd) 1127 if (nval_g < density_cut).any(): 1128 nval_g = nval_g[np.where(nval_g > density_cut)] 1129 print('Not calculating correlation energy ', 1130 'contribution for densities < %3.2e elecs/Bohr ^ 3' % 1131 (density_cut), 1132 file=self.fd) 1133 1134 densitysum = np.sum(nval_g * self.dv) 1135 valence = self.calc.wfs.setups.nvalence 1136 1137 print('Density integrates to %s electrons' % (densitysum), 1138 file=self.fd) 1139 1140 print('Renormalized to %s electrons' % (valence), file=self.fd) 1141 1142 nval_g *= valence / densitysum 1143 self.rs_g = (3.0 / (4.0 * np.pi * nval_g))**(1.0 / 3.0) 1144 1145 self.rsmin = np.amin(self.rs_g) 1146 self.rsmax = np.amax(self.rs_g) 1147 1148 def calculate(self): 1149 1150 print('Generating tables of electron gas energies...', file=self.fd) 1151 1152 table_SR = self.generate_tables() 1153 1154 print('...done', file=self.fd) 1155 # Now interpolate the table to calculate local density terms 1156 E_SR = np.sum(np.interp(self.rs_g, table_SR[:, 0], 1157 table_SR[:, 1])) * self.dv 1158 1159 # RPA energy minus long range correlation 1160 print('Short range correlation energy/unit cell = %5.4f eV \n' % 1161 ((E_SR) * Ha), 1162 file=self.fd) 1163 e = E_SR 1164 1165 return (e) 1166 1167 def generate_tables(self): 1168 1169 # Finite difference steps for density and k vec 1170 rs_step = 0.01 1171 k_step = 0.01 1172 1173 rs_r = np.arange(self.rsmin - rs_step, self.rsmax + rs_step, rs_step) 1174 1175 table_SR = np.empty((len(rs_r), 2)) 1176 table_SR[:, 0] = rs_r 1177 for iR, Rs in enumerate(rs_r): 1178 1179 qF = (9.0 * np.pi / 4.0)**(1.0 / 3.0) * 1.0 / Rs 1180 1181 q_k = np.arange(k_step, 10.0 * qF, k_step) 1182 1183 if self.xc == 'range_RPA': 1184 # Correlation energy per electron, in Hartree, per k 1185 Eeff_k, Erpa_k = self.RPA_corr_hole(q_k, Rs) 1186 ESR_k = Erpa_k - Eeff_k 1187 elif self.xc == 'range_rALDA': 1188 ESR_k = self.rALDA_corr_hole(q_k, Rs) 1189 1190 # Integrate over k 1191 table_SR[iR, 1] = k_step * np.sum(ESR_k) 1192 1193 return table_SR 1194 1195 def RPA_corr_hole(self, q, rs): 1196 1197 # Integrating this quantity over q, gives 1198 # correlation energy per unit volume 1199 # calcuated with a Coulomb-like effective 1200 # interaction, in Hartree 1201 # = 1/(2\pi) * \sum_{\vec{q}} \int_0^infty ds 1202 # * [ ln (1 - v_eff \chi_0) + v_eff \chi_0] 1203 # = 1/(2\pi) * \int 4 \pi q^2 dq /((2\pi)^3) 1204 # * \int_0^infty ds [ ln (1 - v_eff \chi_0) + v_eff \chi_0] 1205 # = 1/(4\pi^3) * \int q^2 dq \int_0^infty ds 1206 # * [ ln (1 - v_eff \chi_0) + v_eff \chi_0] 1207 1208 veff = 4.0 * np.pi / (q * q) * np.exp( 1209 -0.25 * q * q * self.range_rc * self.range_rc) 1210 vc = 4.0 * np.pi / (q * q) 1211 1212 # Do the integral over frequency using Gauss-Legendre 1213 1214 eeff_q = np.zeros(len(q)) 1215 erpa_q = np.zeros(len(q)) 1216 1217 for u, freqweight in zip(self.frequencies, self.freqweights): 1218 1219 chi0 = self.calc_lindhard(q, u, rs) 1220 1221 eff_integrand = np.log(np.ones(len(q)) - veff * chi0) + veff * chi0 1222 eeff_q += eff_integrand * freqweight 1223 1224 rpa_integrand = np.log(np.ones(len(q)) - vc * chi0) + vc * chi0 1225 erpa_q += rpa_integrand * freqweight 1226 1227 # Per unit volume 1228 1229 eeff_q *= 1.0 / (4.0 * np.pi**3.0) * q * q 1230 erpa_q *= 1.0 / (4.0 * np.pi**3.0) * q * q 1231 1232 return (eeff_q, erpa_q) 1233 1234 def rALDA_corr_hole(self, q, rs): 1235 qF = (9.0 * np.pi / 4.0)**(1.0 / 3.0) * 1.0 / rs 1236 1237 veff = 4.0 * np.pi / (q * q) * ((1.0 - 0.25 * q * q / 1238 (qF * qF)) * 0.5 * 1239 (1.0 + np.sign(2.0 * qF - q))) 1240 fxc = veff - 4.0 * np.pi / (q * q) 1241 1242 esr_q = np.zeros(len(q)) 1243 1244 for u, freqweight in zip(self.frequencies, self.freqweights): 1245 1246 chi0 = self.calc_lindhard(q, u, rs) 1247 1248 esr_u = np.zeros(len(q)) 1249 1250 for l, lweight in zip(self.l_l, self.lweights): 1251 1252 chil = chi0 / (1.0 - l * fxc * chi0) 1253 esr_u += lweight * (-fxc) * (chil - chi0) 1254 1255 esr_q += freqweight * esr_u 1256 1257 esr_q *= 1.0 / (4.0 * np.pi**3.0) * q * q 1258 1259 return (esr_q) 1260 1261 def calc_lindhard(self, q, u, rs): 1262 1263 # Calculate Lindhard function at imaginary frequency u 1264 1265 qF = (9.0 * np.pi / 4.0)**(1.0 / 3.0) * 1.0 / rs 1266 1267 Q = q / 2.0 / qF 1268 U = u / q / qF 1269 lchi = ((Q * Q - U * U - 1.0) / 4.0 / Q * np.log( 1270 (U * U + (Q + 1.0) * (Q + 1.0)) / (U * U + (Q - 1.0) * (Q - 1.0)))) 1271 lchi += -1.0 + U * np.arctan((1.0 + Q) / U) + U * np.arctan( 1272 (1.0 - Q) / U) 1273 lchi *= qF / (2.0 * np.pi * np.pi) 1274 return (lchi) 1275 1276 1277class KernelDens: 1278 def __init__(self, calc, xc, ibzq_qc, fd, unit_cells, density_cut, ecut, 1279 tag, timer): 1280 1281 self.calc = calc 1282 self.gd = calc.density.gd 1283 self.xc = xc 1284 self.ibzq_qc = ibzq_qc 1285 self.fd = fd 1286 self.unit_cells = unit_cells 1287 self.density_cut = density_cut 1288 self.ecut = ecut 1289 self.tag = tag 1290 self.timer = timer 1291 1292 self.A_x = -(3 / 4.) * (3 / np.pi)**(1 / 3.) 1293 1294 self.n_g = calc.get_all_electron_density(gridrefinement=1) 1295 self.n_g *= Bohr**3 1296 1297 if xc[-3:] == 'PBE': 1298 nf_g = calc.get_all_electron_density(gridrefinement=2) 1299 nf_g *= Bohr**3 1300 gdf = self.gd.refine() 1301 grad_v = [Gradient(gdf, v, n=1).apply for v in range(3)] 1302 gradnf_vg = gdf.empty(3) 1303 for v in range(3): 1304 grad_v[v](nf_g, gradnf_vg[v]) 1305 self.gradn_vg = gradnf_vg[:, ::2, ::2, ::2] 1306 1307 qd = KPointDescriptor(self.ibzq_qc) 1308 self.pd = PWDescriptor(ecut / Ha, self.gd, complex, qd) 1309 1310 @timer('FHXC') 1311 def calculate_fhxc(self): 1312 1313 print('Calculating %s kernel at %d eV cutoff' % (self.xc, self.ecut), 1314 file=self.fd) 1315 if self.xc[0] == 'r': 1316 self.calculate_rkernel() 1317 else: 1318 assert self.xc[0] == 'A' 1319 self.calculate_local_kernel() 1320 1321 def calculate_rkernel(self): 1322 1323 gd = self.gd 1324 ng_c = gd.N_c 1325 cell_cv = gd.cell_cv 1326 icell_cv = 2 * np.pi * np.linalg.inv(cell_cv) 1327 vol = np.linalg.det(cell_cv) 1328 1329 ns = self.calc.wfs.nspins 1330 n_g = self.n_g # density on rough grid 1331 1332 fx_g = ns * self.get_fxc_g(n_g) # local exchange kernel 1333 qc_g = (-4 * np.pi * ns / fx_g)**0.5 # cutoff functional 1334 flocal_g = qc_g**3 * fx_g / (6 * np.pi**2) # ren. x-kernel for r=r' 1335 Vlocal_g = 2 * qc_g / np.pi # ren. Hartree kernel for r=r' 1336 1337 ng = np.prod(ng_c) # number of grid points 1338 r_vg = gd.get_grid_point_coordinates() 1339 rx_g = r_vg[0].flatten() 1340 ry_g = r_vg[1].flatten() 1341 rz_g = r_vg[2].flatten() 1342 1343 print(' %d grid points and %d plane waves at the Gamma point' % 1344 (ng, self.pd.ngmax), 1345 file=self.fd) 1346 1347 # Unit cells 1348 R_Rv = [] 1349 weight_R = [] 1350 nR_v = self.unit_cells 1351 nR = np.prod(nR_v) 1352 for i in range(-nR_v[0] + 1, nR_v[0]): 1353 for j in range(-nR_v[1] + 1, nR_v[1]): 1354 for h in range(-nR_v[2] + 1, nR_v[2]): 1355 R_Rv.append(i * cell_cv[0] + j * cell_cv[1] + 1356 h * cell_cv[2]) 1357 weight_R.append((nR_v[0] - abs(i)) * (nR_v[1] - abs(j)) * 1358 (nR_v[2] - abs(h)) / float(nR)) 1359 if nR > 1: 1360 # with more than one unit cell only the exchange kernel is 1361 # calculated on the grid. The bare Coulomb kernel is added 1362 # in PW basis and Vlocal_g only the exchange part 1363 dv = self.calc.density.gd.dv 1364 gc = (3 * dv / 4 / np.pi)**(1 / 3.) 1365 Vlocal_g -= 2 * np.pi * gc**2 / dv 1366 print(' Lattice point sampling: ' + '(%s x %s x %s)^2 ' % 1367 (nR_v[0], nR_v[1], nR_v[2]) + 1368 ' Reduced to %s lattice points' % len(R_Rv), 1369 file=self.fd) 1370 1371 l_g_size = -(-ng // mpi.world.size) 1372 l_g_range = range(mpi.world.rank * l_g_size, 1373 min((mpi.world.rank + 1) * l_g_size, ng)) 1374 1375 fhxc_qsGr = {} 1376 for iq in range(len(self.ibzq_qc)): 1377 fhxc_qsGr[iq] = np.zeros( 1378 (ns, len(self.pd.G2_qG[iq]), len(l_g_range)), dtype=complex) 1379 1380 inv_error = np.seterr() 1381 np.seterr(invalid='ignore') 1382 np.seterr(divide='ignore') 1383 1384 t0 = time() 1385 # Loop over Lattice points 1386 for i, R_v in enumerate(R_Rv): 1387 # Loop over r'. f_rr and V_rr are functions of r (dim. as r_vg[0]) 1388 if i == 1: 1389 print( 1390 ' Finished 1 cell in %s seconds' % int(time() - t0) + 1391 ' - estimated %s seconds left' % int( 1392 (len(R_Rv) - 1) * (time() - t0)), 1393 file=self.fd) 1394 self.fd.flush() 1395 if len(R_Rv) > 5: 1396 if (i + 1) % (len(R_Rv) / 5 + 1) == 0: 1397 print(' Finished %s cells in %s seconds' % 1398 (i, int(time() - t0)) + 1399 ' - estimated %s seconds left' % int( 1400 (len(R_Rv) - i) * (time() - t0) / i), 1401 file=self.fd) 1402 self.fd.flush() 1403 for g in l_g_range: 1404 rx = rx_g[g] + R_v[0] 1405 ry = ry_g[g] + R_v[1] 1406 rz = rz_g[g] + R_v[2] 1407 1408 # |r-r'-R_i| 1409 rr = ((r_vg[0] - rx)**2 + (r_vg[1] - ry)**2 + 1410 (r_vg[2] - rz)**2)**0.5 1411 1412 n_av = (n_g + n_g.flatten()[g]) / 2. 1413 fx_g = ns * self.get_fxc_g(n_av, index=g) 1414 qc_g = (-4 * np.pi * ns / fx_g)**0.5 1415 x = qc_g * rr 1416 osc_x = np.sin(x) - x * np.cos(x) 1417 f_rr = fx_g * osc_x / (2 * np.pi**2 * rr**3) 1418 if nR > 1: # include only exchange part of the kernel here 1419 V_rr = (sici(x)[0] * 2 / np.pi - 1) / rr 1420 else: # include the full kernel (also hartree part) 1421 V_rr = (sici(x)[0] * 2 / np.pi) / rr 1422 1423 # Terms with r = r' 1424 if (np.abs(R_v) < 0.001).all(): 1425 tmp_flat = f_rr.flatten() 1426 tmp_flat[g] = flocal_g.flatten()[g] 1427 f_rr = tmp_flat.reshape(ng_c) 1428 tmp_flat = V_rr.flatten() 1429 tmp_flat[g] = Vlocal_g.flatten()[g] 1430 V_rr = tmp_flat.reshape(ng_c) 1431 del tmp_flat 1432 1433 f_rr[np.where(n_av < self.density_cut)] = 0.0 1434 V_rr[np.where(n_av < self.density_cut)] = 0.0 1435 1436 f_rr *= weight_R[i] 1437 V_rr *= weight_R[i] 1438 1439 # r-r'-R_i 1440 r_r = np.array([r_vg[0] - rx, r_vg[1] - ry, r_vg[2] - rz]) 1441 1442 # Fourier transform of r 1443 for iq, q in enumerate(self.ibzq_qc): 1444 q_v = np.dot(q, icell_cv) 1445 e_q = np.exp(-1j * gemmdot(q_v, r_r, beta=0.0)) 1446 f_q = self.pd.fft((f_rr + V_rr) * e_q, iq) * vol / ng 1447 fhxc_qsGr[iq][0, :, g - l_g_range[0]] += f_q 1448 if ns == 2: 1449 f_q = self.pd.fft(V_rr * e_q, iq) * vol / ng 1450 fhxc_qsGr[iq][1, :, g - l_g_range[0]] += f_q 1451 1452 mpi.world.barrier() 1453 1454 np.seterr(**inv_error) 1455 1456 for iq, q in enumerate(self.ibzq_qc): 1457 npw = len(self.pd.G2_qG[iq]) 1458 fhxc_sGsG = np.zeros((ns * npw, ns * npw), complex) 1459 l_pw_size = -(-npw // mpi.world.size) # parallelize over PW below 1460 l_pw_range = range(mpi.world.rank * l_pw_size, 1461 min((mpi.world.rank + 1) * l_pw_size, npw)) 1462 1463 if mpi.world.size > 1: 1464 # redistribute grid and plane waves in fhxc_qsGr[iq] 1465 bg1 = BlacsGrid(mpi.world, 1, mpi.world.size) 1466 bg2 = BlacsGrid(mpi.world, mpi.world.size, 1) 1467 bd1 = bg1.new_descriptor(npw, ng, npw, 1468 -(-ng // mpi.world.size)) 1469 bd2 = bg2.new_descriptor(npw, ng, -(-npw // mpi.world.size), 1470 ng) 1471 1472 fhxc_Glr = np.zeros((len(l_pw_range), ng), dtype=complex) 1473 if ns == 2: 1474 Koff_Glr = np.zeros((len(l_pw_range), ng), dtype=complex) 1475 1476 r = Redistributor(bg1.comm, bd1, bd2) 1477 r.redistribute(fhxc_qsGr[iq][0], fhxc_Glr, npw, ng) 1478 if ns == 2: 1479 r.redistribute(fhxc_qsGr[iq][1], Koff_Glr, npw, ng) 1480 else: 1481 fhxc_Glr = fhxc_qsGr[iq][0] 1482 if ns == 2: 1483 Koff_Glr = fhxc_qsGr[iq][1] 1484 1485 # Fourier transform of r' 1486 for iG in range(len(l_pw_range)): 1487 f_g = fhxc_Glr[iG].reshape(ng_c) 1488 f_G = self.pd.fft(f_g.conj(), iq) * vol / ng 1489 fhxc_sGsG[l_pw_range[0] + iG, :npw] = f_G.conj() 1490 if ns == 2: 1491 v_g = Koff_Glr[iG].reshape(ng_c) 1492 v_G = self.pd.fft(v_g.conj(), iq) * vol / ng 1493 fhxc_sGsG[npw + l_pw_range[0] + iG, :npw] = v_G.conj() 1494 1495 if ns == 2: # f_00 = f_11 and f_01 = f_10 1496 fhxc_sGsG[:npw, npw:] = fhxc_sGsG[npw:, :npw] 1497 fhxc_sGsG[npw:, npw:] = fhxc_sGsG[:npw, :npw] 1498 1499 mpi.world.sum(fhxc_sGsG) 1500 fhxc_sGsG /= vol 1501 1502 if mpi.rank == 0: 1503 w = ulm.open( 1504 'fhxc_%s_%s_%s_%s.ulm' % 1505 (self.tag, self.xc, self.ecut, iq), 'w') 1506 if nR > 1: # add Hartree kernel evaluated in PW basis 1507 Gq2_G = self.pd.G2_qG[iq] 1508 if (q == 0).all(): 1509 Gq2_G = Gq2_G.copy() 1510 Gq2_G[0] = 1. 1511 vq_G = 4 * np.pi / Gq2_G 1512 fhxc_sGsG += np.tile(np.eye(npw) * vq_G, (ns, ns)) 1513 w.write(fhxc_sGsG=fhxc_sGsG) 1514 w.close() 1515 mpi.world.barrier() 1516 print(file=self.fd) 1517 1518 def calculate_local_kernel(self): 1519 # Standard ALDA exchange kernel 1520 # Use with care. Results are very difficult to converge 1521 # Sensitive to density_cut 1522 ns = self.calc.wfs.nspins 1523 gd = self.gd 1524 pd = self.pd 1525 cell_cv = gd.cell_cv 1526 icell_cv = 2 * np.pi * np.linalg.inv(cell_cv) 1527 vol = np.linalg.det(cell_cv) 1528 1529 fxc_sg = ns * self.get_fxc_g(ns * self.n_g) 1530 fxc_sg[np.where(self.n_g < self.density_cut)] = 0.0 1531 1532 r_vg = gd.get_grid_point_coordinates() 1533 1534 for iq in range(len(self.ibzq_qc)): 1535 Gvec_Gc = np.dot(pd.get_reciprocal_vectors(q=iq, add_q=False), 1536 cell_cv / (2 * np.pi)) 1537 npw = len(Gvec_Gc) 1538 l_pw_size = -(-npw // mpi.world.size) 1539 l_pw_range = range(mpi.world.rank * l_pw_size, 1540 min((mpi.world.rank + 1) * l_pw_size, npw)) 1541 fhxc_sGsG = np.zeros((ns * npw, ns * npw), dtype=complex) 1542 for s in range(ns): 1543 for iG in l_pw_range: 1544 for jG in range(npw): 1545 fxc = fxc_sg[s].copy() 1546 dG_c = Gvec_Gc[iG] - Gvec_Gc[jG] 1547 dG_v = np.dot(dG_c, icell_cv) 1548 dGr_g = gemmdot(dG_v, r_vg, beta=0.0) 1549 ft_fxc = gd.integrate(np.exp(-1j * dGr_g) * fxc) 1550 fhxc_sGsG[s * npw + iG, s * npw + jG] = ft_fxc 1551 1552 mpi.world.sum(fhxc_sGsG) 1553 fhxc_sGsG /= vol 1554 1555 Gq2_G = self.pd.G2_qG[iq] 1556 if (self.ibzq_qc[iq] == 0).all(): 1557 Gq2_G[0] = 1. 1558 vq_G = 4 * np.pi / Gq2_G 1559 fhxc_sGsG += np.tile(np.eye(npw) * vq_G, (ns, ns)) 1560 1561 if mpi.rank == 0: 1562 w = ulm.open( 1563 'fhxc_%s_%s_%s_%s.ulm' % 1564 (self.tag, self.xc, self.ecut, iq), 'w') 1565 w.write(fhxc_sGsG=fhxc_sGsG) 1566 w.close() 1567 mpi.world.barrier() 1568 print(file=self.fd) 1569 1570 def get_fxc_g(self, n_g, index=None): 1571 if self.xc[-3:] == 'LDA': 1572 return self.get_lda_g(n_g) 1573 elif self.xc[-3:] == 'PBE': 1574 return self.get_pbe_g(n_g, index=index) 1575 else: 1576 raise '%s kernel not recognized' % self.xc 1577 1578 def get_lda_g(self, n_g): 1579 return (4. / 9.) * self.A_x * n_g**(-2. / 3.) 1580 1581 def get_pbe_g(self, n_g, index=None): 1582 if index is None: 1583 gradn_vg = self.gradn_vg 1584 else: 1585 gradn_vg = self.calc.density.gd.empty(3) 1586 for v in range(3): 1587 gradn_vg[v] = (self.gradn_vg[v] + 1588 self.gradn_vg[v].flatten()[index]) / 2 1589 1590 kf_g = (3. * np.pi**2 * n_g)**(1 / 3.) 1591 s2_g = np.zeros_like(n_g) 1592 for v in range(3): 1593 axpy(1.0, gradn_vg[v]**2, s2_g) 1594 s2_g /= 4 * kf_g**2 * n_g**2 1595 1596 e_g = self.A_x * n_g**(4 / 3.) 1597 v_g = (4 / 3.) * e_g / n_g 1598 f_g = (1 / 3.) * v_g / n_g 1599 1600 kappa = 0.804 1601 mu = 0.2195149727645171 1602 1603 denom_g = (1 + mu * s2_g / kappa) 1604 F_g = 1. + kappa - kappa / denom_g 1605 Fn_g = -mu / denom_g**2 * 8 * s2_g / (3 * n_g) 1606 Fnn_g = -11 * Fn_g / (3 * n_g) - 2 * Fn_g**2 / kappa 1607 1608 fxc_g = f_g * F_g 1609 fxc_g += 2 * v_g * Fn_g 1610 fxc_g += e_g * Fnn_g 1611 """ 1612 # Contributions from varying the gradient 1613 #Fgrad_vg = np.zeros_like(gradn_vg) 1614 #Fngrad_vg = np.zeros_like(gradn_vg) 1615 #for v in range(3): 1616 # axpy(1.0, mu / den_g**2 * gradn_vg[v] / (2 * kf_g**2 * n_g**2), 1617 # Fgrad_vg[v]) 1618 # axpy(-8.0, Fgrad_vg[v] / (3 * n_g), Fngrad_vg[v]) 1619 # axpy(-2.0, Fgrad_vg[v] * Fn_g / kappa, Fngrad_vg[v]) 1620 1621 #tmp = np.zeros_like(fxc_g) 1622 #tmp1 = np.zeros_like(fxc_g) 1623 1624 #for v in range(3): 1625 #self.grad_v[v](Fgrad_vg[v], tmp) 1626 #axpy(-2.0, tmp * v_g, fxc_g) 1627 #for u in range(3): 1628 #self.grad_v[u](Fgrad_vg[u] * tmp, tmp1) 1629 #axpy(-4.0/kappa, tmp1 * e_g, fxc_g) 1630 #self.grad_v[v](Fngrad_vg[v], tmp) 1631 #axpy(-2.0, tmp * e_g, fxc_g) 1632 #self.laplace(mu / den_g**2 / (2 * kf_g**2 * n_g**2), tmp) 1633 #axpy(1.0, tmp * e_g, fxc_g) 1634 """ 1635 1636 return fxc_g 1637 1638 1639""" 1640 def get_fxc_libxc_g(self, n_g): 1641 ### NOT USED AT THE MOMENT 1642 gd = self.calc.density.gd.refine() 1643 1644 xc = XC('GGA_X_' + self.xc[2:]) 1645 #xc = XC('LDA_X') 1646 #sigma = np.zeros_like(n_g).flat[:] 1647 xc.set_grid_descriptor(gd) 1648 sigma_xg, gradn_svg = xc.calculate_sigma(np.array([n_g])) 1649 1650 dedsigma_xg = np.zeros_like(sigma_xg) 1651 e_g = np.zeros_like(n_g) 1652 v_sg = np.array([np.zeros_like(n_g)]) 1653 1654 xc.calculate_gga(e_g, np.array([n_g]), v_sg, sigma_xg, dedsigma_xg) 1655 1656 sigma = sigma_xg[0].flat[:] 1657 gradn_vg = gradn_svg[0] 1658 dedsigma_g = dedsigma_xg[0] 1659 1660 libxc = LibXC('GGA_X_' + self.xc[2:]) 1661 #libxc = LibXC('LDA_X') 1662 libxc.initialize(1) 1663 libxc_fxc = libxc.xc.calculate_fxc_spinpaired 1664 1665 fxc_g = np.zeros_like(n_g).flat[:] 1666 d2edndsigma_g = np.zeros_like(n_g).flat[:] 1667 d2ed2sigma_g = np.zeros_like(n_g).flat[:] 1668 1669 libxc_fxc(n_g.flat[:], fxc_g, sigma, d2edndsigma_g, d2ed2sigma_g) 1670 fxc_g = fxc_g.reshape(np.shape(n_g)) 1671 d2edndsigma_g = d2edndsigma_g.reshape(np.shape(n_g)) 1672 d2ed2sigma_g = d2ed2sigma_g.reshape(np.shape(n_g)) 1673 1674 tmp = np.zeros_like(fxc_g) 1675 tmp1 = np.zeros_like(fxc_g) 1676 1677 #for v in range(3): 1678 #self.grad_v[v](d2edndsigma_g * gradn_vg[v], tmp) 1679 #axpy(-4.0, tmp, fxc_g) 1680 1681 #for u in range(3): 1682 #for v in range(3): 1683 #self.grad_v[v](d2ed2sigma_g * gradn_vg[u] * gradn_vg[v], tmp) 1684 #self.grad_v[u](tmp, tmp1) 1685 #axpy(4.0, tmp1, fxc_g) 1686 1687 #self.laplace(dedsigma_g, tmp) 1688 #axpy(2.0, tmp, fxc_g) 1689 1690 return fxc_g[::2, ::2, ::2] 1691 1692 def get_numerical_fxc_sg(self, n_sg): 1693 ### NOT USED AT THE MOMENT 1694 gd = self.calc.density.gd.refine() 1695 delta = 1.e-4 1696 1697 if self.xc[2:] == 'LDA': 1698 xc = XC('LDA_X') 1699 v1xc_sg = np.zeros_like(n_sg) 1700 v2xc_sg = np.zeros_like(n_sg) 1701 xc.calculate(gd, (1 + delta) * n_sg, v1xc_sg) 1702 xc.calculate(gd, (1 - delta) * n_sg, v2xc_sg) 1703 fxc_sg = (v1xc_sg - v2xc_sg) / (2 * delta * n_sg) 1704 else: 1705 fxc_sg = np.zeros_like(n_sg) 1706 xc = XC('GGA_X_' + self.xc[2:]) 1707 vxc_sg = np.zeros_like(n_sg) 1708 xc.calculate(gd, n_sg, vxc_sg) 1709 for s in range(len(n_sg)): 1710 for x in range(len(n_sg[0])): 1711 for y in range(len(n_sg[0, 0])): 1712 for z in range(len(n_sg[0, 0, 0])): 1713 v1xc_sg = np.zeros_like(n_sg) 1714 n1_sg = n_sg.copy() 1715 n1_sg[s, x, y, z] *= (1 + delta) 1716 xc.calculate(gd, n1_sg, v1xc_sg) 1717 num = v1xc_sg[s, x, y, z] - vxc_sg[s, x, y, z] 1718 den = delta * n_sg[s, x, y, z] 1719 fxc_sg[s, x, y, z] = num / den 1720 1721 return fxc_sg[:, ::2, ::2, ::2] 1722""" 1723 1724 1725def set_flags(self): 1726 """ Based on chosen fxc and av. scheme set up true-false flags """ 1727 1728 if self.xc not in ( 1729 'RPA', 1730 'range_RPA', # range separated RPA a la Bruneval 1731 'rALDA', # renormalized kernels 1732 'rAPBE', 1733 'range_rALDA', 1734 'rALDAns', # no spin (ns) 1735 'rAPBEns', 1736 'rALDAc', # rALDA + correlation 1737 'CP', # Constantin Pitarke 1738 'CP_dyn', # Dynamical form of CP 1739 'CDOP', # Corradini et al 1740 'CDOPs', # CDOP without local term 1741 'JGMs', # simplified jellium-with-gap kernel 1742 'JGMsx', # simplified jellium-with-gap kernel, 1743 # constructed with exchange part only 1744 # so that it scales linearly with l 1745 'ALDA' # standard ALDA 1746 ): 1747 raise RuntimeError('%s kernel not recognized' % self.xc) 1748 1749 if (self.xc == 'rALDA' or self.xc == 'rAPBE' or self.xc == 'ALDA'): 1750 1751 if self.av_scheme is None: 1752 self.av_scheme = 'density' 1753 # Two-point scheme default for rALDA and rAPBE 1754 1755 self.spin_kernel = True 1756 # rALDA/rAPBE are the only kernels which have spin-dependent forms 1757 1758 else: 1759 self.spin_kernel = False 1760 1761 if self.av_scheme == 'density': 1762 assert (self.xc == 'rALDA' or self.xc == 'rAPBE' 1763 or self.xc == 'ALDA'), ('Two-point density average ' + 1764 'only implemented for rALDA and rAPBE') 1765 1766 elif self.xc not in ('RPA', 'range_RPA'): 1767 self.av_scheme = 'wavevector' 1768 else: 1769 self.av_scheme = None 1770 1771 if self.xc in ('rALDAns', 'rAPBEns', 'range_RPA', 'JGMsx', 'RPA', 'rALDA', 1772 'rAPBE', 'range_rALDA', 'ALDA'): 1773 self.linear_kernel = True # Scales linearly with coupling constant 1774 else: 1775 self.linear_kernel = False 1776 1777 if self.xc == 'CP_dyn': 1778 self.dyn_kernel = True 1779 else: 1780 self.dyn_kernel = False 1781 1782 if self.xc == 'JGMs' or self.xc == 'JGMsx': 1783 assert (self.Eg is not None), 'JGMs kernel requires a band gap!' 1784 self.Eg /= Ha # Convert from eV 1785 else: 1786 self.Eg = None 1787