Coverage for gpaw/test/vdw/test_potential.py: 100%

1"""Test correctness of vdW-DF potential.""" 

2import pytest 

3from math import pi 

4from gpaw.grid_descriptor import GridDescriptor 

5import numpy as np 

6from gpaw.xc import XC 

7from gpaw.mpi import world 

8 

9 

10@pytest.mark.libxc 

11def test_vdw_potential(): 

12 N = 10 

13 a = 2.0 

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

15 

16 def paired(): 

17 xc = XC('vdW-DF') 

18 n = 0.3 * np.ones((1, N, N, N)) 

19 n += 0.01 * np.cos(np.arange(N) * 2 * pi / N) 

20 v = 0.0 * n 

21 xc.calculate(gd, n, v) 

22 n2 = 1.0 * n 

23 i = 1 

24 n2[0, i, i, i] += 0.00002 

25 x = v[0, i, i, i] * gd.dv 

26 E2 = xc.calculate(gd, n2, v) 

27 n2[0, i, i, i] -= 0.00004 

28 E2 -= xc.calculate(gd, n2, v) 

29 x2 = E2 / 0.00004 

30 print(i, x, x2, x - x2, x / x2) 

31 assert x == pytest.approx(x2, abs=2e-11) 

32 

33 def polarized(): 

34 xc = XC('vdW-DF') 

35 n = 0.04 * np.ones((2, N, N, N)) 

36 n[1] = 0.3 

37 n[0] += 0.02 * np.sin(np.arange(N) * 2 * pi / N) 

38 n[1] += 0.2 * np.cos(np.arange(N) * 2 * pi / N) 

39 v = 0.0 * n 

40 xc.calculate(gd, n, v) 

41 n2 = 1.0 * n 

42 i = 1 

43 n2[0, i, i, i] += 0.00002 

44 x = v[0, i, i, i] * gd.dv 

45 E2 = xc.calculate(gd, n2, v) 

46 n2[0, i, i, i] -= 0.00004 

47 E2 -= xc.calculate(gd, n2, v) 

48 x2 = E2 / 0.00004 

49 print(i, x, x2, x - x2, x / x2) 

50 assert x == pytest.approx(x2, abs=2e-10) 

51 

52 if world.size == 1: 

53 polarized() 

54 paired()