Coverage for gpaw/test/gllb/test_diamond.py: 100%

1"""This calculation has the following structure. 

2 

31) Calculate the ground state of diamond. 

42) Calculate the band structure of diamond in order to obtain 

5 accurate KS band gap. 

63) Calculate the potential discontinuity using the ground state calculator 

7 and the accurate band gap. 

84) Apply the discontinuity to CBM using the band structure calculator. 

9 

10Compare to reference. 

11""" 

12import pytest 

13from ase.build import bulk 

14from gpaw import GPAW, Davidson, Mixer 

15from gpaw.mpi import world 

16from gpaw.calculator import DeprecatedParameterWarning 

17 

18 

19@pytest.mark.gllb 

20@pytest.mark.libxc 

21@pytest.mark.parametrize('deprecated_syntax', [False, True]) 

22def test_gllb_diamond(in_tmp_dir, deprecated_syntax): 

23 xc = 'GLLBSC' 

24 KS_gap_ref = 4.180237125868162 

25 QP_gap_ref = 5.469387490357182 

26 KS_gap_direct_ref = 5.682571428480924 

27 QP_gap_direct_ref = 7.452331561454596 

28 dxc_ave_ref = 1.33769712051036 

29 # M. Kuisma et. al, https://doi.org/10.1103/PhysRevB.82.115106 

30 # C: KS gap 4.14 eV, QP gap 5.41eV, expt. 5.48 eV 

31 

32 # Calculate ground state 

33 atoms = bulk('C', 'diamond', a=3.567) 

34 # We want sufficiently many grid points that the calculator 

35 # can use wfs.world for the finegd, to test that part of the code. 

36 calc = GPAW(mode='fd', 

37 h=0.2, 

38 kpts=(4, 4, 4), 

39 xc=xc, 

40 nbands=8, 

41 mixer=Mixer(0.5, 5, 10.0), 

42 parallel=dict(domain=min(world.size, 2), band=1), 

43 eigensolver=Davidson(niter=2), 

44 txt='gs.out') 

45 atoms.calc = calc 

46 atoms.get_potential_energy() 

47 

48 if deprecated_syntax: 

49 with pytest.warns(DeprecationWarning): 

50 # Calculate the discontinuity directly 

51 response = calc.hamiltonian.xc.xcs['RESPONSE'] 

52 response.calculate_delta_xc() 

53 KS_gap, dxc = response.calculate_delta_xc_perturbation() 

54 else: 

55 # Calculate the discontinuity directly 

56 homo, lumo = calc.get_homo_lumo() 

57 response = calc.hamiltonian.xc.response 

58 dxc_pot = response.calculate_discontinuity_potential(homo, lumo) 

59 KS_gap, dxc = response.calculate_discontinuity(dxc_pot) 

60 

61 QP_gap = KS_gap + dxc 

62 print('KS gap direct', KS_gap) 

63 print('QP gap direct', QP_gap) 

64 assert KS_gap == pytest.approx(KS_gap_direct_ref, abs=1e-4) 

65 assert QP_gap == pytest.approx(QP_gap_direct_ref, abs=1e-4) 

66 

67 if deprecated_syntax: 

68 with pytest.warns(DeprecatedParameterWarning): 

69 from gpaw import restart 

70 

71 calc.write('gs.gpw') 

72 # Calculate accurate KS-band gap from band structure 

73 atoms, bs_calc = restart('gs.gpw', 

74 fixdensity=True, 

75 kpts={'path': 'GX', 'npoints': 12}, 

76 symmetry='off', 

77 nbands=8, 

78 convergence={'bands': 6}, 

79 eigensolver=Davidson(niter=4), 

80 txt='bs.out') 

81 atoms.get_potential_energy() 

82 

83 # Get the accurate KS-band gap 

84 homolumo = bs_calc.get_homo_lumo() 

85 homo, lumo = homolumo 

86 else: 

87 # Calculate accurate KS-band gap from band structure 

88 bs_calc = calc.fixed_density(kpts={'path': 'GX', 'npoints': 12}, 

89 symmetry='off', 

90 nbands=8, 

91 convergence={'bands': 6}, 

92 eigensolver=Davidson(niter=4), 

93 txt='bs.out') 

94 

95 # Get the accurate KS-band gap 

96 homo, lumo = bs_calc.get_homo_lumo() 

97 

98 if deprecated_syntax: 

99 with pytest.warns(DeprecationWarning): 

100 from ase.units import Ha 

101 from gpaw.xc.gllb.c_response import ResponsePotential 

102 

103 # Calculate the discontinuity potential with accurate band gap 

104 response = calc.hamiltonian.xc.xcs['RESPONSE'] 

105 response.calculate_delta_xc(homolumo=homolumo / Ha) 

106 

107 # Calculate the discontinuity using the band structure calculator 

108 bs_response = bs_calc.hamiltonian.xc.xcs['RESPONSE'] 

109 # Copy response potential to avoid recalculation 

110 pot = ResponsePotential(response, 

111 response.Dxc_vt_sG, 

112 None, 

113 response.Dxc_D_asp, 

114 response.Dxc_Dresp_asp) 

115 pot.redistribute(bs_response) 

116 bs_response.Dxc_vt_sG = pot.vt_sG 

117 bs_response.Dxc_D_asp = pot.D_asp 

118 bs_response.Dxc_Dresp_asp = pot.Dresp_asp 

119 KS_gap, dxc = bs_response.calculate_delta_xc_perturbation() 

120 else: 

121 # Calculate the discontinuity potential with accurate band gap 

122 response = calc.hamiltonian.xc.response 

123 dxc_pot = response.calculate_discontinuity_potential(homo, lumo) 

124 

125 # Calculate the discontinuity using the band structure calculator 

126 bs_response = bs_calc.hamiltonian.xc.response 

127 KS_gap, dxc = bs_response.calculate_discontinuity(dxc_pot) 

128 

129 # Test alternative (testing) implementation 

130 if world.size == 1: 

131 with pytest.warns(Warning): 

132 dxc_ave = bs_response.calculate_discontinuity_from_average( 

133 dxc_pot, 0, True) 

134 assert dxc_ave == pytest.approx(dxc_ave_ref, abs=1e-4) 

135 

136 QP_gap = KS_gap + dxc 

137 print('KS gap', KS_gap) 

138 print('QP gap', QP_gap) 

139 assert KS_gap == pytest.approx(lumo - homo, abs=1e-10) 

140 assert KS_gap == pytest.approx(KS_gap_ref, abs=1e-4) 

141 assert QP_gap == pytest.approx(QP_gap_ref, abs=1e-4)