Coverage for gpaw/test/xc/test_gga_atom.py: 100%

1import numpy as np 

2import numpy.random as ra 

3import pytest 

4from gpaw.setup import create_setup 

5from gpaw.grid_descriptor import GridDescriptor 

6from gpaw.lfc import LFC 

7from gpaw.spline import Spline 

8from gpaw.xc import XC 

9from gpaw.utilities import pack_density 

10from gpaw.mpi import serial_comm 

11 

12 

13def test_xc_gga_atom(): 

14 rng = ra.default_rng(8) 

15 for name in ['LDA', 'PBE']: 

16 xc = XC(name) 

17 s = create_setup('N', xc) 

18 ni = s.ni 

19 nao = s.nao 

20 wt0_j = s.basis_functions_J 

21 

22 rcut = s.xc_correction.rgd.r_g[-1] 

23 

24 wt_j = [] 

25 for wt0 in wt0_j: 

26 data = [wt0(r) for r in np.arange(121) * rcut / 100] 

27 data[-1] = 0.0 

28 l = wt0.get_angular_momentum_number() 

29 wt_j.append(Spline.from_data(l, 1.2 * rcut, data)) 

30 

31 a = rcut * 1.2 * 2 + 1.0 

32 n = 70 

33 n = 90 

34 gd = GridDescriptor((n, n, n), (a, a, a), comm=serial_comm) 

35 pr = LFC(gd, [wt_j]) 

36 pr.set_positions([(0.5, 0.5, 0.5)]) 

37 

38 psit_ig = np.zeros((nao, n, n, n)) 

39 pr.add(psit_ig, {0: np.eye(nao)}) 

40 

41 nii = ni * (ni + 1) // 2 

42 D_p = np.zeros(nii) 

43 

44 e_g = np.zeros((n, n, n)) 

45 n_g = np.zeros((1, n, n, n)) 

46 v_g = np.zeros((1, n, n, n)) 

47 

48 P_ni = 0.2 * rng.random((20, ni)) 

49 P_ni[:, nao:] = 0.0 

50 D_ii = np.dot(np.transpose(P_ni), P_ni) 

51 D_p = pack_density(D_ii) 

52 p = 0 

53 for i1 in range(nao): 

54 for i2 in range(i1, nao): 

55 n_g += D_p[p] * psit_ig[i1] * psit_ig[i2] 

56 p += 1 

57 p += ni - nao 

58 

59 p = LFC(gd, [[s.nct]]) 

60 p.set_positions([(0.5, 0.5, 0.5)]) 

61 p.add(n_g[0], 1.0) 

62 e_g = gd.zeros() 

63 xc.calculate(gd, n_g, v_g, e_g) 

64 

65 r2_g = np.sum((np.indices((n, n, n)) - n / 2)**2, axis=0) 

66 dv_g = gd.dv * np.less(r2_g, (rcut / a * n)**2) 

67 

68 E2 = -np.dot(e_g.ravel(), dv_g.ravel()) 

69 

70 s.xc_correction.n_qg[:] = 0.0 

71 s.xc_correction.nc_g[:] = 0.0 

72 E1 = (xc.calculate_paw_correction(s, D_p.reshape(1, -1)) 

73 + s.xc_correction.e_xc0) 

74 

75 print(name, E1, E2, E1 - E2) 

76 assert E1 == pytest.approx(E2, abs=0.0013)