Coverage for gpaw/test/response/test_paw.py: 100%

1import numpy as np 

2import pytest 

3from functools import partial 

4 

5from ase.units import Ha 

6from ase.data import chemical_symbols 

7 

8from gpaw.response import ResponseGroundStateAdapter 

9from gpaw.response.qpd import SingleQPWDescriptor 

10from gpaw.response.groundstate import ResponsePAWDataset 

11from gpaw.response.paw import (calculate_pair_density_correction, 

12 calculate_matrix_element_correction) 

13from gpaw.response.site_paw import calculate_site_matrix_element_correction 

14from gpaw.response.localft import add_LSDA_trans_fxc 

15 

16from gpaw.setup import create_setup 

17from gpaw.sphere.rshe import calculate_reduced_rshe 

18 

19 

20def setups(): 

21 for symbol in chemical_symbols: 

22 try: 

23 setup = create_setup(symbol) 

24 except FileNotFoundError: 

25 pass 

26 else: 

27 yield setup 

28 

29 

30@pytest.mark.response 

31@pytest.mark.serial 

32@pytest.mark.parametrize('setup', setups()) 

33def test_paw_corrections(setup): 

34 radial_points = 2**10 

35 if setup.symbol in {'I', 'Hg', 'Pb'}: 

36 # More points where needed, for performance. 

37 # https://gitlab.com/gpaw/gpaw/-/issues/984 

38 radial_points *= 4 

39 

40 G_Gv = np.zeros((5, 3)) 

41 G_Gv[:, 0] = np.linspace(0, 20, 5) 

42 pawdata = ResponsePAWDataset(setup, radial_points=radial_points) 

43 calculate_pair_density_correction(G_Gv, pawdata=pawdata) 

44 

45 

46@pytest.mark.response 

47def test_paw_correction_consistency(gpw_files): 

48 """Test consistency of the pair density PAW corrections.""" 

49 gs = ResponseGroundStateAdapter.from_gpw_file(gpw_files['fe_pw']) 

50 

51 # Set up plane-wave basis 

52 ecut = 50 # eV 

53 qpd = SingleQPWDescriptor.from_q([0., 0., 0.5], # N-point 

54 ecut / Ha, gs.gd) 

55 qG_Gv = qpd.get_reciprocal_vectors(add_q=True) 

56 

57 # Calculate ordinary pair density corrections 

58 pawdata = gs.pawdatasets.by_atom[0] 

59 Q1_Gii = calculate_pair_density_correction(qG_Gv, pawdata=pawdata) 

60 

61 # Calculate pair density as a generalized matrix element 

62 # Expand unity in real spherical harmonics 

63 rgd = pawdata.xc_correction.rgd 

64 Y_nL = pawdata.xc_correction.Y_nL 

65 f_ng = np.ones((Y_nL.shape[0], rgd.N)) 

66 rshe, _ = calculate_reduced_rshe(rgd, f_ng, Y_nL, lmax=0) 

67 # Calculate correction 

68 Q2_Gii = calculate_matrix_element_correction(qG_Gv, pawdata, rshe) 

69 

70 assert Q2_Gii == pytest.approx(Q1_Gii, rel=1e-3, abs=1e-5) 

71 

72 

73@pytest.mark.response 

74@pytest.mark.serial 

75def test_site_paw_correction_consistency(gpw_files): 

76 """Test consistency of generalized matrix elements.""" 

77 gs = ResponseGroundStateAdapter.from_gpw_file(gpw_files['fe_pw']) 

78 

79 # Expand the LDA fxc kernel in real spherical harmonics 

80 pawdata = gs.pawdatasets.by_atom[0] 

81 micro_setup = gs.micro_setups[0] 

82 add_fxc = partial(add_LSDA_trans_fxc, fxc='ALDA') 

83 rshe, _ = micro_setup.expand_function(add_fxc, wmin=1e-8) 

84 

85 # Calculate PAW correction with G + q = 0 

86 qG_Gv = np.zeros((1, 3)) 

87 nF_Gii = calculate_matrix_element_correction(qG_Gv, pawdata, rshe) 

88 

89 # Calculate PAW correction with site cutoff exceeding the augmentation 

90 # sphere radius 

91 augr = gs.get_aug_radii()[0] 

92 rcut_p = [augr * 1.5] 

93 drcut = augr * 0.25 

94 lambd_p = [0.5] 

95 nF_pii = calculate_site_matrix_element_correction(pawdata, rshe, 

96 rcut_p, drcut, lambd_p) 

97 

98 assert nF_Gii[0] == pytest.approx(nF_pii[0])