Coverage for /builds/hweiske/ase/ase/calculators/vasp/vasp_auxiliary.py: 13.33%
195 statements
« prev ^ index » next coverage.py v7.2.7, created at 2024-04-22 11:22 +0000
1import os
2import re
4import numpy as np
7def get_vasp_version(string):
8 """Extract version number from header of stdout.
10 Example::
12 >>> get_vasp_version('potato vasp.6.1.2 bumblebee')
13 '6.1.2'
15 """
16 match = re.search(r'vasp\.(\S+)', string, re.M)
17 return match.group(1)
20class VaspChargeDensity:
21 """Class for representing VASP charge density.
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
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
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)
41 def is_spin_polarized(self):
42 if len(self.chgdiff) > 0:
43 return True
44 return False
46 def _read_chg(self, fobj, chg, volume):
47 """Read charge from file object
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.
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
65 def read(self, filename):
66 """Read CHG or CHGCAR file.
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.
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.
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.
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 (OSError, ValueError, IndexError):
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)
137 def _write_chg(self, fobj, chg, volume, format='chg'):
138 """Write charge density
140 Utility function similar to _read_chg but for writing.
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')
185 def write(self, filename, format=None):
186 """Write VASP charge density in CHG format.
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.
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 long_format=False)
213 fd.write('\n')
214 for dim in chg.shape:
215 fd.write(' %4i' % dim)
216 fd.write('\n')
217 vol = self.atoms[ii].get_volume()
218 self._write_chg(fd, chg, vol, format)
219 if format == 'chgcar':
220 fd.write(self.aug)
221 if self.is_spin_polarized():
222 if format == 'chg':
223 fd.write('\n')
224 for dim in chg.shape:
225 fd.write(' %4i' % dim)
226 fd.write('\n') # a new line after dim is required
227 self._write_chg(fd, self.chgdiff[ii], vol, format)
228 if format == 'chgcar':
229 # a new line is always provided self._write_chg
230 fd.write(self.augdiff)
231 if format == 'chg' and len(self.chg) > 1:
232 fd.write('\n')
235class VaspDos:
236 """Class for representing density-of-states produced by VASP
238 The energies are in property self.energy
240 Site-projected DOS is accesible via the self.site_dos method.
242 Total and integrated DOS is accessible as numpy.ndarray's in the
243 properties self.dos and self.integrated_dos. If the calculation is
244 spin polarized, the arrays will be of shape (2, NDOS), else (1,
245 NDOS).
247 The self.efermi property contains the currently set Fermi
248 level. Changing this value shifts the energies.
250 """
252 def __init__(self, doscar='DOSCAR', efermi=0.0):
253 """Initialize"""
254 self._efermi = 0.0
255 self.read_doscar(doscar)
256 self.efermi = efermi
258 # we have determine the resort to correctly map ase atom index to the
259 # POSCAR.
260 self.sort = []
261 self.resort = []
262 if os.path.isfile('ase-sort.dat'):
263 with open('ase-sort.dat') as file:
264 lines = file.readlines()
265 for line in lines:
266 data = line.split()
267 self.sort.append(int(data[0]))
268 self.resort.append(int(data[1]))
270 def _set_efermi(self, efermi):
271 """Set the Fermi level."""
272 ef = efermi - self._efermi
273 self._efermi = efermi
274 self._total_dos[0, :] = self._total_dos[0, :] - ef
275 try:
276 self._site_dos[:, 0, :] = self._site_dos[:, 0, :] - ef
277 except IndexError:
278 pass
280 def _get_efermi(self):
281 return self._efermi
283 efermi = property(_get_efermi, _set_efermi, None, "Fermi energy.")
285 def _get_energy(self):
286 """Return the array with the energies."""
287 return self._total_dos[0, :]
289 energy = property(_get_energy, None, None, "Array of energies")
291 def site_dos(self, atom, orbital):
292 """Return an NDOSx1 array with dos for the chosen atom and orbital.
294 atom: int
295 Atom index
296 orbital: int or str
297 Which orbital to plot
299 If the orbital is given as an integer:
300 If spin-unpolarized calculation, no phase factors:
301 s = 0, p = 1, d = 2
302 Spin-polarized, no phase factors:
303 s-up = 0, s-down = 1, p-up = 2, p-down = 3, d-up = 4, d-down = 5
304 If phase factors have been calculated, orbitals are
305 s, py, pz, px, dxy, dyz, dz2, dxz, dx2
306 double in the above fashion if spin polarized.
308 """
309 # Correct atom index for resorting if we need to. This happens when the
310 # ase-sort.dat file exists, and self.resort is not empty.
311 if self.resort:
312 atom = self.resort[atom]
314 # Integer indexing for orbitals starts from 1 in the _site_dos array
315 # since the 0th column contains the energies
316 if isinstance(orbital, int):
317 return self._site_dos[atom, orbital + 1, :]
318 n = self._site_dos.shape[1]
320 from .vasp_data import PDOS_orbital_names_and_DOSCAR_column
321 norb = PDOS_orbital_names_and_DOSCAR_column[n]
323 return self._site_dos[atom, norb[orbital.lower()], :]
325 def _get_dos(self):
326 if self._total_dos.shape[0] == 3:
327 return self._total_dos[1, :]
328 elif self._total_dos.shape[0] == 5:
329 return self._total_dos[1:3, :]
331 dos = property(_get_dos, None, None, 'Average DOS in cell')
333 def _get_integrated_dos(self):
334 if self._total_dos.shape[0] == 3:
335 return self._total_dos[2, :]
336 elif self._total_dos.shape[0] == 5:
337 return self._total_dos[3:5, :]
339 integrated_dos = property(_get_integrated_dos, None, None,
340 'Integrated average DOS in cell')
342 def read_doscar(self, fname="DOSCAR"):
343 """Read a VASP DOSCAR file"""
344 with open(fname) as fd:
345 natoms = int(fd.readline().split()[0])
346 [fd.readline() for _ in range(4)]
347 # First we have a block with total and total integrated DOS
348 ndos = int(fd.readline().split()[2])
349 dos = []
350 for _ in range(ndos):
351 dos.append(np.array([float(x) for x in fd.readline().split()]))
352 self._total_dos = np.array(dos).T
353 # Next we have one block per atom, if INCAR contains the stuff
354 # necessary for generating site-projected DOS
355 dos = []
356 for _ in range(natoms):
357 line = fd.readline()
358 if line == '':
359 # No site-projected DOS
360 break
361 ndos = int(line.split()[2])
362 line = fd.readline().split()
363 cdos = np.empty((ndos, len(line)))
364 cdos[0] = np.array(line)
365 for nd in range(1, ndos):
366 line = fd.readline().split()
367 cdos[nd] = np.array([float(x) for x in line])
368 dos.append(cdos.T)
369 self._site_dos = np.array(dos)