Coverage for /builds/debichem-team/python-ase/ase/calculators/vasp/vasp_auxiliary.py: 29.74%

1import os 

2import re 

3 

4import numpy as np 

5 

6 

7def get_vasp_version(string): 

8 """Extract version number from header of stdout. 

9 

10 Example:: 

11 

12 >>> get_vasp_version('potato vasp.6.1.2 bumblebee') 

13 '6.1.2' 

14 

15 """ 

16 match = re.search(r'vasp\.(\S+)', string, re.M) 

17 return match.group(1) 

18 

19 

20class VaspChargeDensity: 

21 """Class for representing VASP charge density. 

22 

23 Filename is normally CHG.""" 

24 # Can the filename be CHGCAR? There's a povray tutorial 

25 # in doc/tutorials where it's CHGCAR as of January 2021. --askhl 

26 

27 def __init__(self, filename): 

28 # Instance variables 

29 self.atoms = [] # List of Atoms objects 

30 self.chg = [] # Charge density 

31 self.chgdiff = [] # Charge density difference, if spin polarized 

32 self.aug = '' # Augmentation charges, not parsed just a big string 

33 self.augdiff = '' # Augmentation charge differece, is spin polarized 

34 

35 # Note that the augmentation charge is not a list, since they 

36 # are needed only for CHGCAR files which store only a single 

37 # image. 

38 if filename is not None: 

39 self.read(filename) 

40 

41 def is_spin_polarized(self): 

42 if len(self.chgdiff) > 0: 

43 return True 

44 return False 

45 

46 def _read_chg(self, fobj, chg, volume): 

47 """Read charge from file object 

48 

49 Utility method for reading the actual charge density (or 

50 charge density difference) from a file object. On input, the 

51 file object must be at the beginning of the charge block, on 

52 output the file position will be left at the end of the 

53 block. The chg array must be of the correct dimensions. 

54 

55 """ 

56 # VASP writes charge density as 

57 # WRITE(IU,FORM) (((C(NX,NY,NZ),NX=1,NGXC),NY=1,NGYZ),NZ=1,NGZC) 

58 # Fortran nested implied do loops; innermost index fastest 

59 # First, just read it in 

60 for zz in range(chg.shape[2]): 

61 for yy in range(chg.shape[1]): 

62 chg[:, yy, zz] = np.fromfile(fobj, count=chg.shape[0], sep=' ') 

63 chg /= volume 

64 

65 def read(self, filename): 

66 """Read CHG or CHGCAR file. 

67 

68 If CHG contains charge density from multiple steps all the 

69 steps are read and stored in the object. By default VASP 

70 writes out the charge density every 10 steps. 

71 

72 chgdiff is the difference between the spin up charge density 

73 and the spin down charge density and is thus only read for a 

74 spin-polarized calculation. 

75 

76 aug is the PAW augmentation charges found in CHGCAR. These are 

77 not parsed, they are just stored as a string so that they can 

78 be written again to a CHGCAR format file. 

79 

80 """ 

81 import ase.io.vasp as aiv 

82 with open(filename) as fd: 

83 self.atoms = [] 

84 self.chg = [] 

85 self.chgdiff = [] 

86 self.aug = '' 

87 self.augdiff = '' 

88 while True: 

89 try: 

90 atoms = aiv.read_vasp(fd) 

91 except (KeyError, RuntimeError, ValueError): 

92 # Probably an empty line, or we tried to read the 

93 # augmentation occupancies in CHGCAR 

94 break 

95 fd.readline() 

96 ngr = fd.readline().split() 

97 ng = (int(ngr[0]), int(ngr[1]), int(ngr[2])) 

98 chg = np.empty(ng) 

99 self._read_chg(fd, chg, atoms.get_volume()) 

100 self.chg.append(chg) 

101 self.atoms.append(atoms) 

102 # Check if the file has a spin-polarized charge density part, 

103 # and if so, read it in. 

104 fl = fd.tell() 

105 # First check if the file has an augmentation charge part 

106 # (CHGCAR file.) 

107 line1 = fd.readline() 

108 if line1 == '': 

109 break 

110 elif line1.find('augmentation') != -1: 

111 augs = [line1] 

112 while True: 

113 line2 = fd.readline() 

114 if line2.split() == ngr: 

115 self.aug = ''.join(augs) 

116 augs = [] 

117 chgdiff = np.empty(ng) 

118 self._read_chg(fd, chgdiff, atoms.get_volume()) 

119 self.chgdiff.append(chgdiff) 

120 elif line2 == '': 

121 break 

122 else: 

123 augs.append(line2) 

124 if len(self.aug) == 0: 

125 self.aug = ''.join(augs) 

126 augs = [] 

127 else: 

128 self.augdiff = ''.join(augs) 

129 augs = [] 

130 elif line1.split() == ngr: 

131 chgdiff = np.empty(ng) 

132 self._read_chg(fd, chgdiff, atoms.get_volume()) 

133 self.chgdiff.append(chgdiff) 

134 else: 

135 fd.seek(fl) 

136 

137 def _write_chg(self, fobj, chg, volume, format='chg'): 

138 """Write charge density 

139 

140 Utility function similar to _read_chg but for writing. 

141 

142 """ 

143 # Make a 1D copy of chg, must take transpose to get ordering right 

144 chgtmp = chg.T.ravel() 

145 # Multiply by volume 

146 chgtmp = chgtmp * volume 

147 # Must be a tuple to pass to string conversion 

148 chgtmp = tuple(chgtmp) 

149 # CHG format - 10 columns 

150 if format.lower() == 'chg': 

151 # Write all but the last row 

152 for ii in range((len(chgtmp) - 1) // 10): 

153 fobj.write(' %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G\ 

154 %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G\n' % chgtmp[ii * 10:(ii + 1) * 10]) 

155 # If the last row contains 10 values then write them without a 

156 # newline 

157 if len(chgtmp) % 10 == 0: 

158 fobj.write(' %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G' 

159 ' %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G' % 

160 chgtmp[len(chgtmp) - 10:len(chgtmp)]) 

161 # Otherwise write fewer columns without a newline 

162 else: 

163 for ii in range(len(chgtmp) % 10): 

164 fobj.write((' %#11.5G') % 

165 chgtmp[len(chgtmp) - len(chgtmp) % 10 + ii]) 

166 # Other formats - 5 columns 

167 else: 

168 # Write all but the last row 

169 for ii in range((len(chgtmp) - 1) // 5): 

170 fobj.write(' %17.10E %17.10E %17.10E %17.10E %17.10E\n' % 

171 chgtmp[ii * 5:(ii + 1) * 5]) 

172 # If the last row contains 5 values then write them without a 

173 # newline 

174 if len(chgtmp) % 5 == 0: 

175 fobj.write(' %17.10E %17.10E %17.10E %17.10E %17.10E' % 

176 chgtmp[len(chgtmp) - 5:len(chgtmp)]) 

177 # Otherwise write fewer columns without a newline 

178 else: 

179 for ii in range(len(chgtmp) % 5): 

180 fobj.write((' %17.10E') % 

181 chgtmp[len(chgtmp) - len(chgtmp) % 5 + ii]) 

182 # Write a newline whatever format it is 

183 fobj.write('\n') 

184 

185 def write(self, filename, format=None): 

186 """Write VASP charge density in CHG format. 

187 

188 filename: str 

189 Name of file to write to. 

190 format: str 

191 String specifying whether to write in CHGCAR or CHG 

192 format. 

193 

194 """ 

195 import ase.io.vasp as aiv 

196 if format is None: 

197 if filename.lower().find('chgcar') != -1: 

198 format = 'chgcar' 

199 elif filename.lower().find('chg') != -1: 

200 format = 'chg' 

201 elif len(self.chg) == 1: 

202 format = 'chgcar' 

203 else: 

204 format = 'chg' 

205 with open(filename, 'w') as fd: 

206 for ii, chg in enumerate(self.chg): 

207 if format == 'chgcar' and ii != len(self.chg) - 1: 

208 continue # Write only the last image for CHGCAR 

209 aiv.write_vasp(fd, 

210 self.atoms[ii], 

211 direct=True) 

212 fd.write('\n') 

213 for dim in chg.shape: 

214 fd.write(' %4i' % dim) 

215 fd.write('\n') 

216 vol = self.atoms[ii].get_volume() 

217 self._write_chg(fd, chg, vol, format) 

218 if format == 'chgcar': 

219 fd.write(self.aug) 

220 if self.is_spin_polarized(): 

221 if format == 'chg': 

222 fd.write('\n') 

223 for dim in chg.shape: 

224 fd.write(' %4i' % dim) 

225 fd.write('\n') # a new line after dim is required 

226 self._write_chg(fd, self.chgdiff[ii], vol, format) 

227 if format == 'chgcar': 

228 # a new line is always provided self._write_chg 

229 fd.write(self.augdiff) 

230 if format == 'chg' and len(self.chg) > 1: 

231 fd.write('\n') 

232 

233 

234class VaspDos: 

235 """Class for representing density-of-states produced by VASP 

236 

237 The energies are in property self.energy 

238 

239 Site-projected DOS is accesible via the self.site_dos method. 

240 

241 Total and integrated DOS is accessible as numpy.ndarray's in the 

242 properties self.dos and self.integrated_dos. If the calculation is 

243 spin polarized, the arrays will be of shape (2, NDOS), else (1, 

244 NDOS). 

245 

246 The self.efermi property contains the currently set Fermi 

247 level. Changing this value shifts the energies. 

248 

249 """ 

250 

251 def __init__(self, doscar='DOSCAR', efermi=0.0): 

252 """Initialize""" 

253 self._efermi = 0.0 

254 self.read_doscar(doscar) 

255 self.efermi = efermi 

256 

257 # we have determine the resort to correctly map ase atom index to the 

258 # POSCAR. 

259 self.sort = [] 

260 self.resort = [] 

261 if os.path.isfile('ase-sort.dat'): 

262 with open('ase-sort.dat') as file: 

263 lines = file.readlines() 

264 for line in lines: 

265 data = line.split() 

266 self.sort.append(int(data[0])) 

267 self.resort.append(int(data[1])) 

268 

269 def _set_efermi(self, efermi): 

270 """Set the Fermi level.""" 

271 ef = efermi - self._efermi 

272 self._efermi = efermi 

273 self._total_dos[0, :] = self._total_dos[0, :] - ef 

274 try: 

275 self._site_dos[:, 0, :] = self._site_dos[:, 0, :] - ef 

276 except IndexError: 

277 pass 

278 

279 def _get_efermi(self): 

280 return self._efermi 

281 

282 efermi = property(_get_efermi, _set_efermi, None, "Fermi energy.") 

283 

284 def _get_energy(self): 

285 """Return the array with the energies.""" 

286 return self._total_dos[0, :] 

287 

288 energy = property(_get_energy, None, None, "Array of energies") 

289 

290 def site_dos(self, atom, orbital): 

291 """Return an NDOSx1 array with dos for the chosen atom and orbital. 

292 

293 atom: int 

294 Atom index 

295 orbital: int or str 

296 Which orbital to plot 

297 

298 If the orbital is given as an integer: 

299 If spin-unpolarized calculation, no phase factors: 

300 s = 0, p = 1, d = 2 

301 Spin-polarized, no phase factors: 

302 s-up = 0, s-down = 1, p-up = 2, p-down = 3, d-up = 4, d-down = 5 

303 If phase factors have been calculated, orbitals are 

304 s, py, pz, px, dxy, dyz, dz2, dxz, dx2 

305 double in the above fashion if spin polarized. 

306 

307 """ 

308 # Correct atom index for resorting if we need to. This happens when the 

309 # ase-sort.dat file exists, and self.resort is not empty. 

310 if self.resort: 

311 atom = self.resort[atom] 

312 

313 # Integer indexing for orbitals starts from 1 in the _site_dos array 

314 # since the 0th column contains the energies 

315 if isinstance(orbital, int): 

316 return self._site_dos[atom, orbital + 1, :] 

317 n = self._site_dos.shape[1] 

318 

319 from .vasp_data import PDOS_orbital_names_and_DOSCAR_column 

320 norb = PDOS_orbital_names_and_DOSCAR_column[n] 

321 

322 return self._site_dos[atom, norb[orbital.lower()], :] 

323 

324 def _get_dos(self): 

325 if self._total_dos.shape[0] == 3: 

326 return self._total_dos[1, :] 

327 elif self._total_dos.shape[0] == 5: 

328 return self._total_dos[1:3, :] 

329 

330 dos = property(_get_dos, None, None, 'Average DOS in cell') 

331 

332 def _get_integrated_dos(self): 

333 if self._total_dos.shape[0] == 3: 

334 return self._total_dos[2, :] 

335 elif self._total_dos.shape[0] == 5: 

336 return self._total_dos[3:5, :] 

337 

338 integrated_dos = property(_get_integrated_dos, None, None, 

339 'Integrated average DOS in cell') 

340 

341 def read_doscar(self, fname="DOSCAR"): 

342 """Read a VASP DOSCAR file""" 

343 with open(fname) as fd: 

344 natoms = int(fd.readline().split()[0]) 

345 [fd.readline() for _ in range(4)] 

346 # First we have a block with total and total integrated DOS 

347 ndos = int(fd.readline().split()[2]) 

348 dos = [] 

349 for _ in range(ndos): 

350 dos.append(np.array([float(x) for x in fd.readline().split()])) 

351 self._total_dos = np.array(dos).T 

352 # Next we have one block per atom, if INCAR contains the stuff 

353 # necessary for generating site-projected DOS 

354 dos = [] 

355 for _ in range(natoms): 

356 line = fd.readline() 

357 if line == '': 

358 # No site-projected DOS 

359 break 

360 ndos = int(line.split()[2]) 

361 line = fd.readline().split() 

362 cdos = np.empty((ndos, len(line))) 

363 cdos[0] = np.array(line) 

364 for nd in range(1, ndos): 

365 line = fd.readline().split() 

366 cdos[nd] = np.array([float(x) for x in line]) 

367 dos.append(cdos.T) 

368 self._site_dos = np.array(dos)