Coverage for gpaw/test/radial/test_two_phi_plw_integrals.py: 100%

1import pytest 

2 

3 

4@pytest.mark.response 

5def test_pair_density_paw_correction(): 

6 import numpy as np 

7 from gpaw.lfc import LocalizedFunctionsCollection as LFC 

8 from gpaw.grid_descriptor import GridDescriptor 

9 from gpaw.atom.radialgd import EquidistantRadialGridDescriptor 

10 from gpaw.spline import Spline 

11 from gpaw.response.paw import (calculate_pair_density_correction, 

12 LeanPAWDataset) 

13 # Initialize s, p, d (9 in total) wave and put them on grid 

14 rc = 2.0 

15 a = 2.5 * rc 

16 n = 64 

17 lmax = 2 

18 b = 8.0 

19 m = (lmax + 1)**2 

20 gd = GridDescriptor([n, n, n], [a, a, a]) 

21 r = np.linspace(0, rc, 200) 

22 g = np.exp(-(r / rc * b)**2) 

23 splines = [Spline.from_data(l=l, rmax=rc, f_g=g) for l in range(lmax + 1)] 

24 c = LFC(gd, [splines]) 

25 c.set_positions([(0.5, 0.5, 0.5)]) 

26 psi = gd.zeros(m) 

27 d0 = c.dict(m) 

28 if 0 in d0: 

29 d0[0] = np.identity(m) 

30 c.add(psi, d0) 

31 

32 # Calculate on 3d-grid < phi_i | e**(-ik.r) | phi_j > 

33 R_a = np.array([a / 2, a / 2, a / 2]) 

34 rr = gd.get_grid_point_coordinates() 

35 for dim in range(3): 

36 rr[dim] -= R_a[dim] 

37 

38 k_G = np.array([[0.0, 0.0, 0.0], [1., 0.2, 0.1], [10., 0., 10.]]) 

39 nkpt = k_G.shape[0] 

40 

41 d0 = np.zeros((nkpt, m, m), dtype=complex) 

42 for i in range(m): 

43 for j in range(m): 

44 for ik in range(nkpt): 

45 k = k_G[ik] 

46 # kk = np.sqrt(np.inner(k, k)) 

47 kr = np.inner(k, rr.T).T 

48 expkr = np.exp(-1j * kr) 

49 d0[ik, i, j] = gd.integrate(psi[i] * psi[j] * expkr) 

50 

51 # Calculate on 1d-grid < phi_i | e**(-ik.r) | phi_j > 

52 rgd = EquidistantRadialGridDescriptor(r[1], len(r)) 

53 g = [np.exp(-(r / rc * b)**2) * r**l for l in range(lmax + 1)] 

54 l_j = range(lmax + 1) 

55 rcut_j = [rc] * (lmax + 1) 

56 d1 = calculate_pair_density_correction(k_G, pawdata=LeanPAWDataset( 

57 rgd=rgd, phi_jg=g, phit_jg=np.zeros_like(g), l_j=l_j, rcut_j=rcut_j)) 

58 

59 assert d0 == pytest.approx(d1, abs=1e-8)