Coverage for gpaw/test/xc/test_lxc_xcatom.py: 96%

1import pytest 

2import numpy as np 

3import numpy.random as ra 

4from gpaw.setup import create_setup 

5from gpaw.xc import XC 

6from gpaw.test import gen 

7 

8 

9@pytest.mark.libxc 

10def test_xc_lxc_xcatom(in_tmp_dir): 

11 setups = {} 

12 for functional in [ 

13 'LDA_X', 'LDA_X+LDA_C_PW', 'LDA_X+LDA_C_VWN', 'LDA_X+LDA_C_PZ', 

14 'GGA_X_PBE+GGA_C_PBE', 'GGA_X_PBE_R+GGA_C_PBE', 

15 'GGA_X_B88+GGA_C_P86', 'GGA_X_B88+GGA_C_LYP', 

16 'GGA_X_FT97_A+GGA_C_LYP']: 

17 s = gen('N', xcname=functional) 

18 setups[functional] = s 

19 

20 tolerance = 0.000005 # libxc must reproduce old gpaw energies 

21 # zero Kelvin: in Hartree 

22 

23 reference = { # version 0.9.1 

24 'LDA_X+LDA_C_PW': 2.28836113207, # 'LDA' 

25 'GGA_X_PBE+GGA_C_PBE': 2.3366049993, # 'PBE' 

26 'GGA_X_PBE_R+GGA_C_PBE': 2.34496288319} # 'revPBE' 

27 

28 reference_libxc = { # svnversion 5252 

29 'LDA_X': 1.95030600807, 

30 'LDA_X+LDA_C_PW': 2.23194461135, 

31 'LDA_X+LDA_C_VWN': 2.23297429824, 

32 'LDA_X+LDA_C_PZ': 2.23146045547, 

33 'GGA_X_PBE+GGA_C_PBE': 2.28208665019, 

34 'GGA_X_PBE_R+GGA_C_PBE': 2.29201920843, 

35 'GGA_X_B88+GGA_C_P86': 2.30508027546, 

36 'GGA_X_B88+GGA_C_LYP': 2.28183010548, 

37 'GGA_X_FT97_A+GGA_C_LYP': 2.26846048873} 

38 

39 libxc_set = [ 

40 'LDA_X', 'LDA_X+LDA_C_PW', 'LDA_X+LDA_C_VWN', 'LDA_X+LDA_C_PZ', 

41 'GGA_X_PBE+GGA_C_PBE', 'GGA_X_PBE_R+GGA_C_PBE', 

42 'GGA_X_B88+GGA_C_P86', 'GGA_X_B88+GGA_C_LYP', 

43 'GGA_X_FT97_A+GGA_C_LYP'] 

44 

45 x = 0.000001 

46 for xcname in libxc_set: 

47 # note: using ra.default_rng() not compatible with legacy results 

48 rng = ra.RandomState(8) 

49 xc = XC(xcname) 

50 s = create_setup('N', xc, setupdata=setups[xcname]) 

51 ni = s.ni 

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

53 D_p = 0.1 * rng.random(nii) + 0.4 

54 H_p = np.zeros(nii) 

55 

56 E1 = xc.calculate_paw_correction(s, 

57 D_p.reshape(1, -1), 

58 H_p.reshape(1, -1)) 

59 dD_p = x * rng.random(nii) 

60 D_p += dD_p 

61 dE = np.dot(H_p, dD_p) / x 

62 E2 = xc.calculate_paw_correction(s, D_p.reshape(1, -1)) 

63 print(xcname, dE, (E2 - E1) / x) 

64 assert dE == pytest.approx((E2 - E1) / x, abs=0.003) 

65 

66 E2s = xc.calculate_paw_correction( 

67 s, 

68 np.array([0.5 * D_p, 0.5 * D_p]), 

69 np.array([H_p, H_p])) 

70 print(E2, E2s) 

71 assert E2 == pytest.approx(E2s, abs=1.0e-12) 

72 

73 if xcname in reference: # compare with old gpaw 

74 print('A:', E2, reference[xcname]) 

75 assert E2 == pytest.approx(reference[xcname], abs=tolerance) 

76 

77 if xc in reference_libxc: # compare with reference libxc 

78 print('B:', E2, reference_libxc[xcname]) 

79 assert E2 == pytest.approx(reference_libxc[xcname], abs=tolerance) 

80 

81 D_sp = 0.1 * rng.random((2, nii)) + 0.2 

82 H_sp = np.zeros((2, nii)) 

83 

84 E1 = xc.calculate_paw_correction(s, D_sp, H_sp) 

85 dD_sp = x * rng.random((2, nii)) 

86 D_sp += dD_sp 

87 dE = np.dot(H_sp.ravel(), dD_sp.ravel()) / x 

88 E2 = xc.calculate_paw_correction(s, D_sp, H_sp) 

89 print(dE, (E2 - E1) / x) 

90 assert dE == pytest.approx((E2 - E1) / x, abs=0.005)