Coverage for gpaw/nlopt/linear.py: 98%

47 statements  

« prev     ^ index     » next       coverage.py v7.7.1, created at 2025-07-14 00:18 +0000

1import numpy as np 

2from ase.units import Ha 

3 

4from gpaw.mpi import world 

5 

6 

7def get_chi_tensor( 

8 nlodata, 

9 freqs=[1.0], eta=0.05, 

10 ftol=1e-4, Etol=1e-6, eshift=0.0, 

11 band_n=None, out_name=None): 

12 """ 

13 Calculate full linear susceptibility tensor for nonmagnetic semiconductors; 

14 array will be saved to disk if out_name is given. 

15 

16 Parameters 

17 ---------- 

18 nlodata 

19 Data object of type NLOData. 

20 freqs 

21 Excitation frequency array (a numpy array or list). 

22 eta 

23 Broadening, a number or an array (default 0.05 eV). 

24 Etol, ftol 

25 Tolerance in energy and occupancy to consider degeneracy. 

26 eshift 

27 Bandgap correction. 

28 band_n 

29 List of bands in the sum (default 0 to nb). 

30 out_name 

31 If it is given: output filename. 

32 

33 Returns 

34 ------- 

35 np.ndarray: 

36 Full linear susceptibility tensor (3, 3, nw). 

37 

38 """ 

39 

40 # Covert inputs in eV to Ha 

41 freqs = np.array(freqs) 

42 nw = len(freqs) 

43 w_lc = (freqs + 1j * eta) / Ha 

44 Etol /= Ha 

45 eshift /= Ha 

46 

47 # Load the required data 

48 k_info = nlodata.distribute() 

49 if k_info: 

50 tmp = list(k_info.values())[0] 

51 nb = len(tmp[1]) 

52 if band_n is None: 

53 band_n = list(range(nb)) 

54 

55 # Initialize the outputs 

56 sum_vvl = np.zeros((3, 3, nw), complex) 

57 

58 # Do the calculations 

59 for _, (we, f_n, E_n, p_vnn) in k_info.items(): 

60 tmp = np.zeros((3, 3, nw), complex) 

61 for v1 in range(3): 

62 for v2 in range(v1, 3): 

63 sum_l = calc_chi( 

64 w_lc, f_n, E_n, p_vnn, [v1, v2], 

65 band_n, ftol, Etol, eshift=eshift) 

66 tmp[v1, v2] = sum_l 

67 tmp[v2, v1] = sum_l 

68 # Add it to previous with a weight 

69 sum_vvl += tmp * we 

70 

71 world.sum(sum_vvl) 

72 

73 # Multiply prefactors (4pi from eps0, 8pi^3 from BZ) 

74 prefactor = 4 * np.pi / (2 * np.pi)**3 

75 chi_vvl = prefactor * sum_vvl 

76 

77 # Save it to the file 

78 if world.rank == 0 and out_name is not None: 

79 tmp = chi_vvl.reshape(9, nw) 

80 lin = np.vstack((freqs, tmp)) 

81 np.save(out_name, lin) 

82 

83 return chi_vvl 

84 

85 

86def calc_chi( 

87 w_l, f_n, E_n, p_vnn, pol_v, 

88 band_n=None, ftol=1e-4, Etol=1e-6, eshift=0): 

89 """ 

90 Loop over bands for computing the response 

91 

92 Input: 

93 w_l Complex frequency array 

94 f_n Fermi levels 

95 E_n Energies 

96 p_vnn Momentum matrix elements 

97 pol_v Tensor element 

98 band_n Band list 

99 Etol, ftol Tol. in energy and fermi to consider degeneracy 

100 eshift Bandgap correction 

101 Output: 

102 sum_l Output sum value (array with length of w_l) 

103 """ 

104 

105 # Initialize variables 

106 nb = len(f_n) 

107 if band_n is None: 

108 band_n = list(range(nb)) 

109 sum_l = np.zeros(w_l.size, complex) 

110 

111 # Loop over bands 

112 for nni in band_n: 

113 for mmi in band_n: 

114 # Use TRS to reduce 

115 if mmi <= nni: 

116 continue 

117 fnm = f_n[nni] - f_n[mmi] 

118 Emn = E_n[mmi] - E_n[nni] + fnm * eshift 

119 if np.abs(fnm) < ftol or np.abs(Emn) < Etol: 

120 continue 

121 # *2 for real, /2 for TRS, *2 for m<n 

122 sum_l += 2 * fnm * np.real( 

123 p_vnn[pol_v[0], nni, mmi] * p_vnn[pol_v[1], mmi, nni]) \ 

124 / (Emn * (Emn**2 - w_l**2)) 

125 

126 return sum_l