Coverage for gpaw/response/site_paw.py: 100%

1import numpy as np 

2 

3from gpaw.gaunt import gaunt 

4from gpaw.sphere.rshe import RealSphericalHarmonicsExpansion 

5from gpaw.sphere.integrate import radial_truncation_function 

6 

7 

8def calculate_site_matrix_element_correction( 

9 pawdata, rshe: RealSphericalHarmonicsExpansion, 

10 rcut_p, drcut, lambd_p): 

11 r"""Calculate the PAW correction to a site matrix element. 

12 

13 For the site matrix element 

14 

15 f^ap_nn' = <ψ_n|Θ(r∊Ω_ap)f(r)|ψ_n'> 

16 

17 the PAW correction for each pair of partial waves i and i' is given by 

18 l 

19 __ __ 

20 ap \ \ m,mi,mi' / a a ̰ a ̰ a 

21 F = / / g | r^2 dr θ(r<rc) f (r) [φ(r) φ(r) - φ(r) φ(r)] 

22 ii' ‾‾ ‾‾ l,li,li' / p lm i i' i i' 

23 l m=-l 

24 

25 where g refer to the Gaunt coefficients and f_lm(r) are the real 

26 spherical harmonics expansion coefficients of the function f(r). 

27 

28 Here, we evaluate the correction F_ii'^ap for various smooth truncation 

29 functions θ_p(r<rc), parametrized by rc, Δrc and λ. 

30 """ 

31 rgd = rshe.rgd 

32 assert rgd is pawdata.xc_correction.rgd 

33 ni = pawdata.ni # Number of partial waves 

34 l_j = pawdata.l_j # l-index for each radial function index j 

35 lmax = max(l_j) 

36 G_LLL = gaunt(lmax) 

37 assert max(rshe.l_M) <= lmax * 2 

38 # (Real) radial functions for the partial waves 

39 phi_jg = pawdata.phi_jg 

40 phit_jg = pawdata.phit_jg 

41 # Truncate the radial functions to span only the radial grid coordinates 

42 # which need correction 

43 assert np.allclose(rgd.r_g - pawdata.rgd.r_g[:rgd.N], 0.) 

44 phi_jg = np.array(phi_jg)[:, :rgd.N] 

45 phit_jg = np.array(phit_jg)[:, :rgd.N] 

46 

47 # Calculate smooth truncation functions and allocate array 

48 Np = len(rcut_p) 

49 assert len(lambd_p) == Np 

50 theta_pg = [radial_truncation_function(rgd.r_g, rcut, drcut, lambd) 

51 for rcut, lambd in zip(rcut_p, lambd_p)] 

52 F_pii = np.zeros((Np, ni, ni), dtype=float) 

53 

54 # Loop of radial function indices for partial waves i and i' 

55 i1_counter = 0 

56 for j1, l1 in enumerate(l_j): 

57 i2_counter = 0 

58 for j2, l2 in enumerate(l_j): 

59 # Calculate the radial partial wave correction 

60 dn_g = phi_jg[j1] * phi_jg[j2] - phit_jg[j1] * phit_jg[j2] 

61 

62 # Generate m-indices for each radial function 

63 for m1 in range(2 * l1 + 1): 

64 for m2 in range(2 * l2 + 1): 

65 # Set up the i=(l,m) index for each partial wave 

66 i1 = i1_counter + m1 

67 i2 = i2_counter + m2 

68 

69 # Loop through the real spherical harmonics of the local 

70 # function f(r) 

71 for L, f_g in zip(rshe.L_M, rshe.f_gM.T): 

72 # Angular integral 

73 gaunt_coeff = G_LLL[l1**2 + m1, l2**2 + m2, L] 

74 if gaunt_coeff == 0: 

75 continue 

76 # Radial integral 

77 for p, theta_g in enumerate(theta_pg): 

78 F_pii[p, i1, i2] += \ 

79 gaunt_coeff * rgd.integrate_trapz( 

80 theta_g * f_g * dn_g) 

81 

82 # Add to i and i' counters 

83 i2_counter += 2 * l2 + 1 

84 i1_counter += 2 * l1 + 1 

85 return F_pii