Coverage for gpaw/test/exx/test_derivs.py: 78%

1import numpy as np 

2import pytest 

3from ase import Atoms 

4 

5from gpaw.grid_descriptor import GridDescriptor 

6from gpaw.hybrids.kpts import PWKPoint 

7from gpaw.kpt_descriptor import KPointDescriptor 

8from gpaw.mpi import world 

9from gpaw.projections import Projections 

10from gpaw.pw.descriptor import PWDescriptor 

11from gpaw.pw.lfc import PWLFC 

12from gpaw.hybrids.wstc import WignerSeitzTruncatedCoulomb as WSTC 

13from gpaw.spline import Spline 

14from gpaw.symmetry import Symmetry 

15from gpaw.wavefunctions.arrays import PlaneWaveExpansionWaveFunctions 

16 

17 

18class AP: 

19 my_indices = [0] 

20 comm = world 

21 rank_a = [0] 

22 

23 

24r2 = np.linspace(0, 1, 51)**2 

25 

26 

27class Setup: 

28 Delta_iiL = np.zeros((1, 1, 1)) + 0.1 

29 X_p = np.zeros(1) + 0.3 

30 ExxC = -10.0 

31 ghat_l = [Spline.from_data(0, 1.0, 1 - r2 * (1 - 2 * r2))] 

32 

33 

34@pytest.mark.xfail # XXX: reason=... 

35def test_exx_derivs(): 

36 if world.size > 1: 

37 from unittest import SkipTest 

38 raise SkipTest 

39 

40 N = 20 

41 L = 2.5 

42 nb = 2 

43 spos_ac = np.zeros((1, 3)) + 0.25 

44 

45 gd = GridDescriptor([N, N, N], np.eye(3) * L) 

46 sym = Symmetry([], gd.cell_cv) 

47 kd = KPointDescriptor(None) 

48 kd.set_symmetry(Atoms(pbc=True), sym) 

49 coulomb = WSTC(gd.cell_cv, kd.N_c) 

50 pd = PWDescriptor(10, gd, complex, kd) 

51 

52 data = pd.zeros(nb) 

53 data[0, 1] = 3.0 

54 data[1, 2] = -2.5 

55 psit = PlaneWaveExpansionWaveFunctions(nb, pd, data=data) 

56 

57 proj = Projections(nb, [1], AP(), world, dtype=complex) 

58 

59 pt = PWLFC([[Spline.from_data(0, 1.0, 1 - r2 * (1 - 2 * r2))]], pd) 

60 pt.set_positions(spos_ac) 

61 

62 f_n = np.array([1.0, 0.5]) 

63 kpt = PWKPoint(psit, proj, f_n, np.array([0.0, 0.0, 0.0]), 1.0) 

64 

65 # xx = EXX(kd, [Setup()], pt, coulomb, spos_ac) 

66 xx = (kd, [Setup()], pt, coulomb, spos_ac) 

67 

68 psit.matrix_elements(pt, out=proj) 

69 C = 0.79 

70 VV_aii = {a: np.einsum('n, ni, nj -> ij', f_n, P_ni, P_ni.conj()) * C 

71 for a, P_ni in proj.items()} 

72 

73 x = xx.calculate([kpt], [kpt], VV_aii, derivatives=True) 

74 v = x[3][0][0, 1] 

75 

76 eps = 0.00001 

77 

78 data[0, 1] = 3 + eps 

79 psit.matrix_elements(pt, out=proj) 

80 VV_aii = {a: np.einsum('n, ni, nj -> ij', f_n, P_ni, P_ni.conj()) * C 

81 for a, P_ni in proj.items()} 

82 xp = xx.calculate([kpt], [kpt], VV_aii) 

83 

84 data[0, 1] = 3 - eps 

85 psit.matrix_elements(pt, out=proj) 

86 VV_aii = {a: np.einsum('n, ni, nj -> ij', f_n, P_ni, P_ni.conj()) * C 

87 for a, P_ni in proj.items()} 

88 xm = xx.calculate([kpt], [kpt], VV_aii) 

89 

90 d = (xp[0] + xp[1] - xm[0] - xm[1]) / (2 * eps) * N**6 / L**3 / 2 

91 assert abs(v - d) < 1e-10, (v, d)