Coverage for /builds/kinetik161/ase/ase/calculators/vasp/vasp_auxiliary.py: 13.13%
198 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +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 fd = open(filename)
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, and
103 # if so, read it in.
104 fl = fd.tell()
105 # First check if the file has an augmentation charge part (CHGCAR
106 # 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 fd.close()
138 def _write_chg(self, fobj, chg, volume, format='chg'):
139 """Write charge density
141 Utility function similar to _read_chg but for writing.
143 """
144 # Make a 1D copy of chg, must take transpose to get ordering right
145 chgtmp = chg.T.ravel()
146 # Multiply by volume
147 chgtmp = chgtmp * volume
148 # Must be a tuple to pass to string conversion
149 chgtmp = tuple(chgtmp)
150 # CHG format - 10 columns
151 if format.lower() == 'chg':
152 # Write all but the last row
153 for ii in range((len(chgtmp) - 1) // 10):
154 fobj.write(' %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G\
155 %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G\n' % chgtmp[ii * 10:(ii + 1) * 10])
156 # If the last row contains 10 values then write them without a
157 # newline
158 if len(chgtmp) % 10 == 0:
159 fobj.write(' %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G'
160 ' %#11.5G %#11.5G %#11.5G %#11.5G %#11.5G' %
161 chgtmp[len(chgtmp) - 10:len(chgtmp)])
162 # Otherwise write fewer columns without a newline
163 else:
164 for ii in range(len(chgtmp) % 10):
165 fobj.write((' %#11.5G') %
166 chgtmp[len(chgtmp) - len(chgtmp) % 10 + ii])
167 # Other formats - 5 columns
168 else:
169 # Write all but the last row
170 for ii in range((len(chgtmp) - 1) // 5):
171 fobj.write(' %17.10E %17.10E %17.10E %17.10E %17.10E\n' %
172 chgtmp[ii * 5:(ii + 1) * 5])
173 # If the last row contains 5 values then write them without a
174 # newline
175 if len(chgtmp) % 5 == 0:
176 fobj.write(' %17.10E %17.10E %17.10E %17.10E %17.10E' %
177 chgtmp[len(chgtmp) - 5:len(chgtmp)])
178 # Otherwise write fewer columns without a newline
179 else:
180 for ii in range(len(chgtmp) % 5):
181 fobj.write((' %17.10E') %
182 chgtmp[len(chgtmp) - len(chgtmp) % 5 + ii])
183 # Write a newline whatever format it is
184 fobj.write('\n')
186 def write(self, filename, format=None):
187 """Write VASP charge density in CHG format.
189 filename: str
190 Name of file to write to.
191 format: str
192 String specifying whether to write in CHGCAR or CHG
193 format.
195 """
196 import ase.io.vasp as aiv
197 if format is None:
198 if filename.lower().find('chgcar') != -1:
199 format = 'chgcar'
200 elif filename.lower().find('chg') != -1:
201 format = 'chg'
202 elif len(self.chg) == 1:
203 format = 'chgcar'
204 else:
205 format = 'chg'
206 with open(filename, 'w') as fd:
207 for ii, chg in enumerate(self.chg):
208 if format == 'chgcar' and ii != len(self.chg) - 1:
209 continue # Write only the last image for CHGCAR
210 aiv.write_vasp(fd,
211 self.atoms[ii],
212 direct=True,
213 long_format=False)
214 fd.write('\n')
215 for dim in chg.shape:
216 fd.write(' %4i' % dim)
217 fd.write('\n')
218 vol = self.atoms[ii].get_volume()
219 self._write_chg(fd, chg, vol, format)
220 if format == 'chgcar':
221 fd.write(self.aug)
222 if self.is_spin_polarized():
223 if format == 'chg':
224 fd.write('\n')
225 for dim in chg.shape:
226 fd.write(' %4i' % dim)
227 fd.write('\n') # a new line after dim is required
228 self._write_chg(fd, self.chgdiff[ii], vol, format)
229 if format == 'chgcar':
230 # a new line is always provided self._write_chg
231 fd.write(self.augdiff)
232 if format == 'chg' and len(self.chg) > 1:
233 fd.write('\n')
236class VaspDos:
237 """Class for representing density-of-states produced by VASP
239 The energies are in property self.energy
241 Site-projected DOS is accesible via the self.site_dos method.
243 Total and integrated DOS is accessible as numpy.ndarray's in the
244 properties self.dos and self.integrated_dos. If the calculation is
245 spin polarized, the arrays will be of shape (2, NDOS), else (1,
246 NDOS).
248 The self.efermi property contains the currently set Fermi
249 level. Changing this value shifts the energies.
251 """
253 def __init__(self, doscar='DOSCAR', efermi=0.0):
254 """Initialize"""
255 self._efermi = 0.0
256 self.read_doscar(doscar)
257 self.efermi = efermi
259 # we have determine the resort to correctly map ase atom index to the
260 # POSCAR.
261 self.sort = []
262 self.resort = []
263 if os.path.isfile('ase-sort.dat'):
264 file = open('ase-sort.dat')
265 lines = file.readlines()
266 file.close()
267 for line in lines:
268 data = line.split()
269 self.sort.append(int(data[0]))
270 self.resort.append(int(data[1]))
272 def _set_efermi(self, efermi):
273 """Set the Fermi level."""
274 ef = efermi - self._efermi
275 self._efermi = efermi
276 self._total_dos[0, :] = self._total_dos[0, :] - ef
277 try:
278 self._site_dos[:, 0, :] = self._site_dos[:, 0, :] - ef
279 except IndexError:
280 pass
282 def _get_efermi(self):
283 return self._efermi
285 efermi = property(_get_efermi, _set_efermi, None, "Fermi energy.")
287 def _get_energy(self):
288 """Return the array with the energies."""
289 return self._total_dos[0, :]
291 energy = property(_get_energy, None, None, "Array of energies")
293 def site_dos(self, atom, orbital):
294 """Return an NDOSx1 array with dos for the chosen atom and orbital.
296 atom: int
297 Atom index
298 orbital: int or str
299 Which orbital to plot
301 If the orbital is given as an integer:
302 If spin-unpolarized calculation, no phase factors:
303 s = 0, p = 1, d = 2
304 Spin-polarized, no phase factors:
305 s-up = 0, s-down = 1, p-up = 2, p-down = 3, d-up = 4, d-down = 5
306 If phase factors have been calculated, orbitals are
307 s, py, pz, px, dxy, dyz, dz2, dxz, dx2
308 double in the above fashion if spin polarized.
310 """
311 # Correct atom index for resorting if we need to. This happens when the
312 # ase-sort.dat file exists, and self.resort is not empty.
313 if self.resort:
314 atom = self.resort[atom]
316 # Integer indexing for orbitals starts from 1 in the _site_dos array
317 # since the 0th column contains the energies
318 if isinstance(orbital, int):
319 return self._site_dos[atom, orbital + 1, :]
320 n = self._site_dos.shape[1]
322 from .vasp_data import PDOS_orbital_names_and_DOSCAR_column
323 norb = PDOS_orbital_names_and_DOSCAR_column[n]
325 return self._site_dos[atom, norb[orbital.lower()], :]
327 def _get_dos(self):
328 if self._total_dos.shape[0] == 3:
329 return self._total_dos[1, :]
330 elif self._total_dos.shape[0] == 5:
331 return self._total_dos[1:3, :]
333 dos = property(_get_dos, None, None, 'Average DOS in cell')
335 def _get_integrated_dos(self):
336 if self._total_dos.shape[0] == 3:
337 return self._total_dos[2, :]
338 elif self._total_dos.shape[0] == 5:
339 return self._total_dos[3:5, :]
341 integrated_dos = property(_get_integrated_dos, None, None,
342 'Integrated average DOS in cell')
344 def read_doscar(self, fname="DOSCAR"):
345 """Read a VASP DOSCAR file"""
346 fd = open(fname)
347 natoms = int(fd.readline().split()[0])
348 [fd.readline() for nn in range(4)] # Skip next 4 lines.
349 # First we have a block with total and total integrated DOS
350 ndos = int(fd.readline().split()[2])
351 dos = []
352 for nd in range(ndos):
353 dos.append(np.array([float(x) for x in fd.readline().split()]))
354 self._total_dos = np.array(dos).T
355 # Next we have one block per atom, if INCAR contains the stuff
356 # necessary for generating site-projected DOS
357 dos = []
358 for na in range(natoms):
359 line = fd.readline()
360 if line == '':
361 # No site-projected DOS
362 break
363 ndos = int(line.split()[2])
364 line = fd.readline().split()
365 cdos = np.empty((ndos, len(line)))
366 cdos[0] = np.array(line)
367 for nd in range(1, ndos):
368 line = fd.readline().split()
369 cdos[nd] = np.array([float(x) for x in line])
370 dos.append(cdos.T)
371 self._site_dos = np.array(dos)
372 fd.close()