Coverage for gpaw/response/g0w0_kernels.py: 97%

1from ase.units import Ha 

2import numpy as np 

3from gpaw.xc.fxc import KernelWave, XCFlags, FXCCache 

4from gpaw.xc.rpa import GCut 

5 

6from gpaw.response.qpd import SingleQPWDescriptor 

7from gpaw.pw.descriptor import PWMapping 

8 

9 

10class G0W0Kernel: 

11 def __init__(self, xc, context, **kwargs): 

12 self.xc = xc 

13 self.context = context 

14 self.xcflags = XCFlags(xc) 

15 self._kwargs = kwargs 

16 

17 def calculate(self, qpd): 

18 if self.xc == 'RPA': 

19 return np.eye(qpd.ngmax) 

20 

21 return calculate_spinkernel( 

22 qpd=qpd, 

23 xcflags=self.xcflags, 

24 context=self.context, 

25 **self._kwargs) 

26 

27 

28def calculate_spinkernel(*, ecut, xcflags, gs, qd, qpd, context): 

29 assert xcflags.spin_kernel 

30 xc = xcflags.xc 

31 ns = gs.nspins 

32 ibzq_qc = qd.ibzk_kc 

33 iq = np.argmin(np.linalg.norm(ibzq_qc - qpd.q_c[np.newaxis], axis=1)) 

34 assert np.allclose(ibzq_qc[iq], qpd.q_c) 

35 

36 ecut_max = ecut * Ha # XXX very ugly this 

37 

38 cache = FXCCache(comm=context.comm, 

39 tag=gs.atoms.get_chemical_formula(mode='hill'), 

40 xc=xc, ecut=ecut_max) 

41 handle = cache.handle(iq) 

42 

43 if not handle.exists(): 

44 # Somehow we calculated many q even though this function 

45 # only works on one q? Very confusing. 

46 kernel = KernelWave( 

47 q_empty=iq, ibzq_qc=qd.ibzk_kc, 

48 xc=xcflags.xc, 

49 ecut=ecut_max, gs=gs, 

50 context=context) 

51 

52 # The first time we miss the cache, we calculate /all/ iq. 

53 # (Whether that's the best strategy can be discussed.) 

54 for iq_calculated, array in kernel.calculate_fhxc(): 

55 cache.handle(iq_calculated).write(array) 

56 

57 fv = handle.read() 

58 assert fv is not None 

59 

60 # If we want a reduced plane-wave description, create qpd mapping 

61 if qpd.ecut < ecut: 

62 # Recreate nonreduced plane-wave description corresponding to ecut_max 

63 qpdnr = SingleQPWDescriptor.from_q(qpd.q_c, ecut, qpd.gd, 

64 gammacentered=qpd.gammacentered) 

65 pw_map = PWMapping(qpd, qpdnr) 

66 gcut = GCut(pw_map.G2_G1) 

67 fv = gcut.spin_cut(fv, ns=ns) 

68 

69 return fv