Coverage for gpaw/cli/dos.py: 91%

1"""CLI-code for dos-subcommand.""" 

2from pathlib import Path 

3from typing import Union, List, Tuple, Optional 

4 

5import numpy as np 

6 

7from gpaw import GPAW 

8from gpaw.setup import Setup 

9from gpaw.spherical_harmonics import names as ylmnames 

10from gpaw.dos import DOSCalculator 

11 

12 

13class CLICommand: 

14 """Calculate (projected) density of states from gpw-file.""" 

15 

16 @staticmethod 

17 def add_arguments(parser): 

18 add = parser.add_argument 

19 add('gpw', metavar='gpw-file') 

20 add('-p', '--plot', action='store_true', 

21 help='Plot the DOS.') 

22 add('-i', '--integrated', action='store_true', 

23 help='Calculate integrated DOS.') 

24 add('-w', '--width', type=float, default=0.1, 

25 help='Width of Gaussian. Use 0.0 to use linear tetrahedron ' 

26 'interpolation.') 

27 add('-a', '--atom', help='Project onto atoms: "Cu-spd,H-s" or use ' 

28 'atom indices "12-spdf". Particular m-values can be obtained ' 

29 'like this: "N-p0,N-p1,N-p2. For p-orbitals, m=0,1,2 translates ' 

30 'to y, z and x. For d-orbitals, m=0,1,2,3,4 translates ' 

31 'to xy, yz, 3z2-r2, zx and x2-y2.') 

32 add('-t', '--total', action='store_true', 

33 help='Show both PDOS and total DOS.') 

34 add('-r', '--range', nargs=2, metavar=('emin', 'emax'), 

35 help='Energy range in eV relative to Fermi level.') 

36 add('-n', '--points', type=int, default=400, help='Number of points.') 

37 add('--soc', action='store_true', 

38 help='Include spin-orbit coupling.') 

39 

40 @staticmethod 

41 def run(args): 

42 if args.range is None: 

43 emin = None 

44 emax = None 

45 else: 

46 emin, emax = (float(e) for e in args.range) 

47 dos(args.gpw, 

48 plot=args.plot, 

49 width=args.width, 

50 integrated=args.integrated, 

51 projection=args.atom, 

52 soc=args.soc, 

53 emin=emin, 

54 emax=emax, 

55 npoints=args.points, 

56 show_total=args.total) 

57 

58 

59def parse_projection_string(projection: str, 

60 symbols: List[str], 

61 setups: List[Setup] 

62 ) -> List[Tuple[str, 

63 List[Tuple[int, 

64 int, 

65 Union[None, int]]]]]: 

66 """Create labels and lists of (a, l, m)-tuples. 

67 

68 Example: 

69 

70 >>> from gpaw.setup import create_setup 

71 >>> setup = create_setup('Li') 

72 >>> s, py = parse_projection_string('Li-sp0', 

73 ... ['Li', 'Li'], 

74 ... [setup, setup]) 

75 >>> s 

76 ('Li-s', [(0, 0, None), (1, 0, None)]) 

77 >>> py 

78 ('Li-p(y)', [(0, 1, 0), (1, 1, 0)]) 

79 

80 * "Li-s" will have contributions from l=0 and atoms 0 and 1 

81 * "Li-p(y)" will have contributions from l=1, m=0 and atoms 0 and 1 

82 

83 """ 

84 result: List[Tuple[str, List[Tuple[int, int, Union[None, int]]]]] = [] 

85 for proj in projection.split(','): 

86 s, ll = proj.split('-') 

87 if s.isdigit(): 

88 A = [int(s)] 

89 s = '#' + s 

90 else: 

91 A = [a for a, symbol in 

92 enumerate(symbols) 

93 if symbol == s] 

94 if not A: 

95 raise ValueError('No such atom: ' + s) 

96 for l, m in parse_lm_string(ll): 

97 label = s + '-' + 'spdfg'[l] 

98 if m is not None: 

99 name = ylmnames[l][m] 

100 label += f'({name})' 

101 result.append((label, [(a, l, m) for a in A])) 

102 

103 return result 

104 

105 

106def parse_lm_string(s: str) -> List[Tuple[int, Union[None, int]]]: 

107 """Parse 'spdf' kind of string to numbers. 

108 

109 Return list of (l, m) tuples with m=None if not specified: 

110 

111 >>> parse_lm_string('sp') 

112 [(0, None), (1, None)] 

113 >>> parse_lm_string('p0p1p2') 

114 [(1, 0), (1, 1), (1, 2)] 

115 """ 

116 result = [] 

117 while s: 

118 l = 'spdfg'.index(s[0]) 

119 m: Union[None, int] 

120 if s[1:2].isnumeric(): 

121 m = int(s[1:2]) 

122 s = s[2:] 

123 else: 

124 m = None 

125 s = s[1:] 

126 result.append((l, m)) 

127 return result 

128 

129 

130def dos(filename: Union[Path, str], 

131 *, 

132 plot=False, 

133 width=0.1, 

134 integrated=False, 

135 projection=None, 

136 soc=False, 

137 emin=None, 

138 emax=None, 

139 npoints=200, 

140 show_total=None): 

141 """Calculate density of states. 

142 

143 filename: str 

144 Name of restart-file. 

145 plot: bool 

146 Show a plot. 

147 width: float 

148 Width of Gaussians. 

149 integrated: bool 

150 Calculate integrated DOS. 

151 

152 """ 

153 from ase.spectrum.dosdata import GridDOSData 

154 from ase.spectrum.doscollection import GridDOSCollection 

155 

156 calc = GPAW(filename) 

157 

158 doscalc = DOSCalculator.from_calculator(calc, soc) 

159 energies = doscalc.get_energies(emin, emax, npoints) 

160 nspins = doscalc.nspins 

161 spinlabels = [''] if nspins == 1 else [' up', ' dn'] 

162 spins: List[Optional[int]] = [None] if nspins == 1 else [0, 1] 

163 

164 dosobjs = GridDOSCollection([], energies) 

165 

166 # Note: ignoring wrong GridDOSCollection.__add__ hint below. 

167 

168 if projection is None or show_total: 

169 for spin, label in zip(spins, spinlabels): 

170 dos = doscalc.raw_dos(energies, spin=spin, width=width) 

171 dosobjs += GridDOSData(energies, dos, # type: ignore 

172 {'label': 'DOS' + label}) 

173 

174 if projection is not None: 

175 symbols = calc.atoms.get_chemical_symbols() 

176 projs = parse_projection_string( 

177 projection, symbols, calc.setups) 

178 for label, contributions in projs: 

179 for spin, spinlabel in zip(spins, spinlabels): 

180 dos = np.zeros_like(energies) 

181 for a, l, m in contributions: 

182 dos += doscalc.raw_pdos(energies, 

183 a, l, m, 

184 spin=spin, 

185 width=width) 

186 dosobjs += GridDOSData(energies, dos, # type: ignore 

187 {'label': label + spinlabel}) 

188 

189 if integrated: 

190 de = energies[1] - energies[0] 

191 energies = energies + de / 2 

192 dosobjs = GridDOSCollection( 

193 [GridDOSData(energies, 

194 np.cumsum(obj.get_weights()) * de, 

195 obj.info) 

196 for obj in dosobjs]) 

197 ylabel = 'iDOS [electrons]' 

198 else: 

199 ylabel = 'DOS [electrons/eV]' 

200 

201 if plot: 

202 ax = dosobjs.plot() 

203 ax.set_xlabel(r'$\epsilon-\epsilon_F$ [eV]') 

204 ax.set_ylabel(ylabel) 

205 import matplotlib.pyplot as plt 

206 plt.show() 

207 

208 return dosobjs