Coverage for gpaw/xc/ri/ribasis.py: 15%

1from dataclasses import replace 

2from gpaw.basis_data import BasisFunction 

3from collections import defaultdict 

4 

5 

6def generate_ri_basis(basis, accuracy): 

7 lmax = 2 

8 

9 # TODO: Hartree 

10 def poisson(n_g, l): 

11 return Hartree(basis.rgd, n_g, l) 

12 

13 # Auxiliary basis functions per angular momentum channel 

14 auxt_lng = defaultdict(lambda: []) 

15 

16 ribf_j = [] 

17 

18 def add(aux_g, l, rc=None): 

19 ribf = BasisFunction(n=None, l=l, rc=rc, phit_g=aux_g, 

20 type='auxiliary') 

21 ribf_j.append(ribf) 

22 

23 def basisloop(): 

24 for j, bf in enumerate(basis.bf_j): 

25 yield j, bf.l, bf.rc, bf.phit_g 

26 

27 # Double basis function loop to create product orbitals 

28 for j1, l1, rc1, phit1_g in basisloop(): 

29 for j2, l2, rc2, phit2_g in basisloop(): 

30 # Loop only over ordered pairs 

31 if j1 > j2: 

32 continue 

33 

34 # Loop over all possible angular momentum states what the 

35 # product l1 x l2 creates. 

36 for l in range((l1 + l2) % 2, l1 + l2 + 1, 2): 

37 if l > lmax: 

38 continue 

39 

40 # min: The support of basis function is the intersection 

41 # of the individual supports. 

42 add(phit1_g * phit2_g, l, rc=min(rc1, rc2)) 

43 

44 for l, auxt_ng in auxt_lng.items(): 

45 print(l, auxt_ng) 

46 print(f' l={l}') 

47 for n, auxt_g in enumerate(auxt_ng): 

48 print(f' {n}') 

49 

50 return replace(basis, ribf_j=ribf_j) 

51 

52 

53def Hartree(rgd, n_g, l): 

54 v_g = rgd.poisson(n_g, l) 

55 v_g[1:] /= rgd.r_g[1:] 

56 v_g[0] = v_g[1] 

57 return v_g