Coverage for gpaw/response/jdos.py: 100%

1import numpy as np 

2 

3from gpaw.response import ResponseContext 

4from gpaw.response.frequencies import ComplexFrequencyDescriptor 

5from gpaw.response.pair_integrator import (PairFunctionIntegrator, 

6 DynamicPairFunction) 

7from gpaw.response.chiks import get_temporal_part 

8 

9 

10class JDOS(DynamicPairFunction): 

11 

12 def __init__(self, spincomponent, q_c, zd): 

13 self.spincomponent = spincomponent 

14 super().__init__(q_c, zd) 

15 

16 def zeros(self): 

17 nz = len(self.zd) 

18 return np.zeros(nz, float) 

19 

20 

21class JDOSCalculator(PairFunctionIntegrator): 

22 r"""Joint density of states calculator. 

23 

24 Here, the joint density of states of collinear systems is defined as the 

25 spectral part of the four-component Kohn-Sham susceptibility, 

26 see [PRB 103, 245110 (2021)]: 

27 

28 __ __ 

29 1 \ \ / 

30 g^μν(q,ω) = ‾ / / | σ^μ_ss' σ^ν_s's (f_nks - f_n'k+qs') 

31 V ‾‾ ‾‾ \ 

32 k t \ 

33 x δ(ħω - [ε_n'k's'-ε_nks]) | 

34 / 

35 

36 where t is a composite band and spin transition index: (n, s) -> (n', s'). 

37 """ 

38 

39 def __init__(self, gs, context=None, 

40 nbands=None, bandsummation='pairwise', 

41 **kwargs): 

42 """Contruct the JDOSCalculator 

43 

44 Parameters 

45 ---------- 

46 gs : ResponseGroundStateAdapter 

47 context : ResponseContext 

48 nbands : int 

49 Number of bands to include in the sum over states 

50 bandsummation : str 

51 Band summation strategy (does not change the result, but can affect 

52 the run-time). 

53 'pairwise': sum over pairs of bands 

54 'double': double sum over band indices. 

55 kwargs : see gpaw.response.pair_integrator.PairFunctionIntegrator 

56 """ 

57 if context is None: 

58 context = ResponseContext() 

59 assert isinstance(context, ResponseContext) 

60 

61 super().__init__(gs, context, **kwargs) 

62 

63 self.nbands = nbands 

64 self.bandsummation = bandsummation 

65 

66 def calculate(self, spincomponent, q_c, zd): 

67 """Calculate g^μν(q,ω) using a lorentzian broadening of the δ-function 

68 

69 Parameters 

70 ---------- 

71 spincomponent : str 

72 Spin component (μν) of the joint density of states. 

73 Currently, '00', 'uu', 'dd', '+-' and '-+' are implemented. 

74 q_c : list or np.array 

75 Wave vector in relative coordinates 

76 zd : ComplexFrequencyDescriptor 

77 Complex frequencies ħz to evaluate g^μν(q,ω) at, where z = ω + iη 

78 and η > 0 yields the HWHM of the resulting lorentzian broadening. 

79 """ 

80 assert isinstance(zd, ComplexFrequencyDescriptor) 

81 assert zd.upper_half_plane 

82 

83 # Prepare to sum over bands and spins 

84 transitions = self.get_band_and_spin_transitions( 

85 spincomponent, nbands=self.nbands, 

86 bandsummation=self.bandsummation) 

87 

88 self.context.print(self.get_info_string( 

89 q_c, len(zd), spincomponent, self.nbands, len(transitions))) 

90 

91 # Set up output data structure 

92 jdos = JDOS(spincomponent, q_c, zd) 

93 

94 # Perform actual in-place integration 

95 self.context.print('Integrating the joint density of states:') 

96 self._integrate(jdos, transitions) 

97 

98 return jdos 

99 

100 def add_integrand(self, kptpair, weight, jdos): 

101 r"""Add the g^μν(q,ω) integrand of the outer k-point integral: 

102 __ 

103 -1 \ σ^μ_ss' σ^ν_s's (f_nks - f_n'k's') 

104 (...)_k = ‾‾ Im / ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾ 

105 π ‾‾ ħω - (ε_n'k's' - ε_nks) + iħη 

106 t 

107 

108 NB: Since the implemented spin matrices are real, the dissipative part 

109 is equal to the imaginary part (up to a factor of π) of the full 

110 integrand. 

111 """ 

112 # Construct jdos integrand via the imaginary part of the frequency 

113 # dependence in χ_KS^μν(q,z) 

114 if jdos.spincomponent == '00' and self.gs.nspins == 1: 

115 weight = 2 * weight 

116 x_mytz = get_temporal_part(jdos.spincomponent, jdos.zd.hz_z, 

117 kptpair, self.bandsummation) 

118 integrand_mytz = -x_mytz.imag / np.pi 

119 

120 with self.context.timer('Perform sum over t-transitions'): 

121 # Sum over local transitions 

122 integrand_z = np.sum(integrand_mytz, axis=0) 

123 # Sum over distributed transitions 

124 kptpair.tblocks.blockcomm.sum(integrand_z) 

125 jdos.array[:] += weight * integrand_z 

126 

127 def get_info_string(self, q_c, nz, spincomponent, nbands, nt): 

128 """Get information about the joint density of states calculation""" 

129 info_list = ['', 

130 'Calculating the joint density of states with:', 

131 f' q_c: [{q_c[0]}, {q_c[1]}, {q_c[2]}]', 

132 f' Number of frequency points: {nz}', 

133 f' Spin component: {spincomponent}', 

134 self.get_band_and_transitions_info_string(nbands, nt), 

135 '', 

136 self.get_basic_info_string()] 

137 return '\n'.join(info_list)