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

1from math import pi 

2 

3import numpy as np 

4import pytest 

5 

6from gpaw.atom.radialgd import EquidistantRadialGridDescriptor 

7from gpaw.grid_descriptor import GridDescriptor 

8from gpaw.xc import XC 

9 

10 

11@pytest.mark.ci 

12def test_xc_XC2(): 

13 rgd = EquidistantRadialGridDescriptor(0.01, 100) 

14 

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

16 xc = XC(name) 

17 for nspins in [1, 2]: 

18 n = rgd.zeros(nspins) 

19 v = rgd.zeros(nspins) 

20 n[:] = np.exp(-rgd.r_g**2) 

21 n[-1] *= 2 

22 E = xc.calculate_spherical(rgd, n, v) 

23 i = 23 

24 x = v[-1, i] * rgd.dv_g[i] 

25 n[-1, i] += 0.000001 

26 Ep = xc.calculate_spherical(rgd, n, v) 

27 n[-1, i] -= 0.000002 

28 Em = xc.calculate_spherical(rgd, n, v) 

29 x2 = (Ep - Em) / 0.000002 

30 print(name, nspins, E, x, x2, x - x2) 

31 assert x == pytest.approx(x2, abs=1e-9) 

32 n[-1, i] += 0.000001 

33 if nspins == 1: 

34 ns = rgd.empty(2) 

35 ns[:] = n / 2 

36 Es = xc.calculate_spherical(rgd, ns, 0 * ns) 

37 assert E == pytest.approx(Es, abs=1e-13) 

38 

39 N = 20 

40 a = 1.0 

41 gd = GridDescriptor((N, N, N), (a, a, a)) 

42 

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

44 xc = XC(name) 

45 for nspins in [1, 2]: 

46 n = gd.empty(nspins) 

47 n.fill(0.03) 

48 z = np.arange(gd.beg_c[2], gd.end_c[2]) * a / N 

49 n[:] += 0.01 * np.sin(2 * pi * z / a) 

50 if nspins == 2: 

51 n[1] += 0.01 * np.cos(2 * pi * z / a) 

52 n /= nspins 

53 

54 v = 0.0 * n 

55 E = xc.calculate(gd, n, v) 

56 

57 here = (gd.beg_c[0] <= 1 < gd.end_c[0] and 

58 gd.beg_c[1] <= 2 < gd.end_c[1] and 

59 gd.beg_c[2] <= 3 < gd.end_c[2]) 

60 if here: 

61 x = v[-1, 1, 2, 3] * gd.dv 

62 n[-1, 1, 2, 3] += 0.000001 

63 Ep = xc.calculate(gd, n, v) 

64 if here: 

65 n[-1, 1, 2, 3] -= 0.000002 

66 Em = xc.calculate(gd, n, v) 

67 x2 = (Ep - Em) / 0.000002 

68 if here: 

69 print(name, nspins, E, x, x2, x - x2) 

70 assert x == pytest.approx(x2, abs=1e-11) 

71 n[-1, 1, 2, 3] += 0.000001 

72 

73 if nspins == 1: 

74 ns = gd.empty(2) 

75 ns[:] = n / 2 

76 Es = xc.calculate(gd, ns, 0 * ns) 

77 assert E == pytest.approx(Es, abs=1e-13)