Coverage for gpaw/lrtddft/dielectric.py: 92%

1import numpy as np 

2 

3import gpaw 

4from ase.units import Hartree, Bohr 

5from gpaw.utilities.folder import Folder 

6 

7 

8def get_dielectric(exlist, volume, 

9 emin=0, 

10 emax=None, 

11 de=None, 

12 energyunit='eV', 

13 width=0.08, # width in energyunit 

14 comment=None, 

15 form='v' 

16 ): 

17 """Evaluate the dielectric function 

18 

19 Parameters: 

20 =============== =================================================== 

21 ``exlist`` ExcitationList 

22 ``volume`` Unit cell volume in Angstrom**3 

23 ``emin`` min. energy, set to zero if not given 

24 ``emax`` max. energy, set to cover all energies if not given 

25 ``de`` energy spacing 

26 ``energyunit`` Energy unit 'eV' or 'nm', default 'eV' 

27 ``width`` folding width in terms of the chosen energyunit 

28 =============== =================================================== 

29 all energies in [eV] 

30 

31 Returns 

32 ======= 

33 energy: energy values 

34 eps1: real part of the dielectric function 

35 eps2: imaginary part of the dielectric function 

36 N: real part of the refractive index 

37 K: imaginary part of the refractive index 

38 R: reflection coefficient 

39 """ 

40 

41 x = [] 

42 y = [] 

43 for ex in exlist: 

44 x.append(ex.get_energy() * Hartree) 

45 y.append(ex.get_oscillator_strength(form)) 

46 

47 if energyunit != 'eV': 

48 raise RuntimeError('currently only eV is supported') 

49 

50 pre = 4. * np.pi / volume * Bohr**3 * Hartree**2 

51 energies, values = Folder(width, 'RealLorentzPole').fold( 

52 x, y, de, emin, emax) 

53 eps1 = 1 + pre * values.T[0] 

54 energies, values = Folder(width, 'ImaginaryLorentzPole').fold( 

55 x, y, de, emin, emax) 

56 eps2 = pre * values.T[0] 

57 

58 N = np.sqrt(0.5 * (np.sqrt(eps1**2 + eps2**2) + eps1)) 

59 K = np.sqrt(0.5 * (np.sqrt(eps1**2 + eps2**2) - eps1)) 

60 R = ((N - 1)**2 + K**2) / ((N + 1)**2 + K**2) 

61 

62 return energies, eps1, eps2, N, K, R 

63 

64 

65def dielectric(exlist, 

66 volume, 

67 filename=None, 

68 emin=0, 

69 emax=None, 

70 de=None, 

71 energyunit='eV', 

72 width=0.08, # width in energyunit 

73 comment=None, 

74 form='v' 

75 ): 

76 """Write out the dielectric function 

77 

78 Parameters: 

79 =============== =================================================== 

80 ``exlist`` ExcitationList 

81 ``volume`` Unit cell volume in Angstrom**3 

82 ``filename`` File name for the output file, STDOUT if not given 

83 ``emin`` min. energy, set to zero if not given 

84 ``emax`` max. energy, set to cover all energies if not given 

85 ``de`` energy spacing 

86 ``energyunit`` Energy unit 'eV' or 'nm', default 'eV' 

87 ``width`` folding width in terms of the chosen energyunit 

88 =============== =================================================== 

89 all energies in [eV] 

90 """ 

91 

92 with open(filename, 'w') as out: 

93 if comment: 

94 print('#', comment, file=out) 

95 

96 print('# Dielec function', file=out) 

97 print('# GPAW version:', gpaw.__version__, file=out) 

98 print(f'# width={width} [{energyunit}]', file=out) 

99 if form == 'r': 

100 print('# length form', file=out) 

101 else: 

102 assert form == 'v' 

103 print('# velocity form', file=out) 

104 print( 

105 '# om [{}] eps1/eps0 eps2/eps0 n k R' 

106 .format(energyunit), file=out) 

107 

108 energies, eps1, eps2, N, K, R = get_dielectric( 

109 exlist=exlist, volume=volume, emin=emin, emax=emax, de=de, 

110 energyunit=energyunit, width=width, form=form) 

111 

112 for e, e1, e2, n, k, r in zip(energies, eps1, eps2, N, K, R): 

113 print('%10.5f %12.7e %12.7e %12.7e %12.7e %12.7e' % 

114 (e, e1, e2, n, k, r), file=out) 

115 

116 if filename is not None: 

117 out.close()