Coverage for gpaw/test/response/test_gw_si.py: 85%

1"""Test GW band-gaps for Si.""" 

2 

3import pytest 

4from ase.build import bulk 

5 

6from gpaw import GPAW 

7from gpaw.mpi import world 

8from gpaw.response.g0w0 import G0W0 

9from gpaw.response.screened_interaction import GammaIntegrationMode 

10 

11 

12def generate_si_systems(): 

13 a = 5.43 

14 si1 = bulk('Si', 'diamond', a=a) 

15 si2 = si1.copy() 

16 si2.positions -= a / 8 

17 

18 return [si1, si2] 

19 

20 

21def run(gpw_filename, nblocks, integrate_gamma, qpt=False): 

22 # This tests checks the actual numerical accuracy which is asserted below 

23 calc = GPAW(gpw_filename) 

24 e = calc.get_potential_energy() 

25 

26 integrate_gamma = GammaIntegrationMode(integrate_gamma) 

27 # The numerical integration default is too slow, so overriding 

28 integrate_gamma._N = 30 

29 

30 kwargs = dict(nbands=8, integrate_gamma=integrate_gamma, 

31 kpts=[(0, 0, 0), (0.5, 0.5, 0)], # Gamma, X 

32 ecut=40, nblocks=nblocks, 

33 frequencies={'type': 'nonlinear', 

34 'domega0': 0.1, 'omegamax': None}, 

35 eta=0.2, relbands=(-1, 2)) 

36 

37 if qpt: 

38 # This part of the code is testing for separate calculation of qpoints 

39 # which would help in trivial parallelization of GW 

40 gw = G0W0(gpw_filename, 'gw_None', **kwargs) 

41 for q in range(gw.nqpts): 

42 gw.calculate(qpoints=[q]) 

43 

44 gw = G0W0(gpw_filename, 'gw_None', **kwargs) 

45 results = gw.calculate() 

46 

47 G, X = results['eps'][0] 

48 output = [e, G[0], G[1] - G[0], X[1] - G[0], X[2] - X[1]] 

49 G, X = results['qp'][0] 

50 output += [G[0], G[1] - G[0], X[1] - G[0], X[2] - X[1]] 

51 return output 

52 

53 

54reference = {'sphere': pytest.approx([-9.253, 5.442, 2.389, 0.403, 0.000, 

55 6.261, 3.570, 1.323, 0.001], abs=0.0035), 

56 'WS': pytest.approx([-9.253, 5.442, 2.389, 0.403, 0.000, 

57 6.284, 3.551, 1.285, 0.001], abs=0.0035), 

58 '1BZ': pytest.approx([-9.252, 5.441, 2.389, 0.403, 0.000, 

59 6.337, 3.450, 1.193, 0.002], abs=0.0035), 

60 'reciprocal': pytest.approx([-9.252, 5.441, 2.389, 0.403, 0.000, 

61 6.110, 3.86, 1.624, 0.002], 

62 abs=0.0035)} 

63 

64# The systems are not 2D, thus, the reciprocal2D will yield same results as 

65# reciprocal. This is tested in test_integrate_gamma_modes. 

66reference['reciprocal2D'] = reference['reciprocal'] 

67reference['1BZ2D'] = reference['1BZ'] 

68 

69 

70@pytest.mark.response 

71@pytest.mark.slow 

72@pytest.mark.parametrize('si', [0, 1]) 

73@pytest.mark.parametrize('integrate_gamma', ['sphere', 'WS']) 

74@pytest.mark.parametrize('symm', ['all', 'no', 'tr', 'pg']) 

75@pytest.mark.parametrize('nblocks', 

76 [x for x in [1, 2, 4, 8] if x <= world.size]) 

77def test_response_gwsi(in_tmp_dir, si, symm, nblocks, integrate_gamma, 

78 scalapack, gpw_files): 

79 filename = gpw_files[f'si_gw_a{si}_{symm}'] 

80 assert run(filename, nblocks, integrate_gamma) ==\ 

81 reference[integrate_gamma] 

82 

83 

84@pytest.mark.parametrize('integrate_gamma', ['sphere', 'WS', '1BZ', 

85 'reciprocal', 

86 'reciprocal2D', 

87 '1BZ2D']) 

88@pytest.mark.response 

89def test_integrate_gamma_modes(in_tmp_dir, integrate_gamma, gpw_files): 

90 assert run(gpw_files['si_gw_a0_all'], 1, integrate_gamma) == \ 

91 reference[integrate_gamma] 

92 

93 

94@pytest.mark.response 

95@pytest.mark.ci 

96@pytest.mark.parametrize('si', [0, 1]) 

97@pytest.mark.parametrize('symm', ['all']) 

98def test_small_response_gwsi(in_tmp_dir, si, symm, scalapack, 

99 gpw_files): 

100 filename = gpw_files[f'si_gw_a{si}_{symm}'] 

101 assert run(filename, 1, 'sphere') == reference['sphere'] 

102 

103 

104@pytest.mark.response 

105@pytest.mark.ci 

106def test_few_freq_response_gwsi(in_tmp_dir, scalapack, 

107 gpw_files): 

108 if world.size > 1: 

109 nblocks = 2 

110 else: 

111 nblocks = 1 

112 

113 # This test has very few frequencies and tests that the code doesn't crash. 

114 filename = gpw_files['si_gw_a0_all'] 

115 gw = G0W0(filename, 'gw_0.2', 

116 nbands=8, integrate_gamma='sphere', 

117 kpts=[(0, 0, 0), (0.5, 0.5, 0)], # Gamma, X 

118 ecut=40, nblocks=nblocks, 

119 frequencies={'type': 'nonlinear', 

120 'domega0': 0.1, 'omegamax': 0.2}, 

121 eta=0.2, relbands=(-1, 2)) 

122 gw.calculate()