Coverage for ase / io / vasp.py: 90.70%

1# fmt: off 

2 

3""" 

4This module contains functionality for reading and writing an ASE 

5Atoms object in VASP POSCAR format. 

6 

7""" 

8from __future__ import annotations 

9 

10import re 

11from pathlib import Path 

12from typing import TextIO 

13 

14import numpy as np 

15 

16from ase import Atoms 

17from ase.constraints import FixAtoms, FixedLine, FixedPlane, FixScaled 

18from ase.io import ParseError 

19from ase.io.formats import string2index 

20from ase.io.utils import ImageIterator 

21from ase.symbols import Symbols 

22from ase.units import Ang, fs 

23from ase.utils import reader, writer 

24 

25from .vasp_parsers import vasp_outcar_parsers as vop 

26 

27__all__ = [ 

28 'read_vasp', 'read_vasp_out', 'iread_vasp_out', 'read_vasp_xdatcar', 

29 'read_vasp_xml', 'write_vasp', 'write_vasp_xdatcar' 

30] 

31 

32EVTOJ = 1.60217733E-19 # VASP constant.F 

33 

34 

35def parse_poscar_scaling_factor(line: str) -> np.ndarray: 

36 """Parse scaling factor(s) in the second line in POSCAR/CONTCAR. 

37 

38 This can also be one negative number or three positive numbers. 

39 

40 https://www.vasp.at/wiki/index.php/POSCAR#Full_format_specification 

41 

42 """ 

43 scale = [] 

44 for _ in line.split()[:3]: 

45 try: 

46 scale.append(float(_)) 

47 except ValueError: 

48 break 

49 if len(scale) not in {1, 3}: 

50 raise RuntimeError('The number of scaling factors must be 1 or 3.') 

51 if len(scale) == 3 and any(_ < 0.0 for _ in scale): 

52 raise RuntimeError('All three scaling factors must be positive.') 

53 return np.array(scale) 

54 

55 

56def get_atomtypes(fname): 

57 """Given a file name, get the atomic symbols. 

58 

59 The function can get this information from OUTCAR and POTCAR 

60 format files. The files can also be compressed with gzip or 

61 bzip2. 

62 

63 """ 

64 fpath = Path(fname) 

65 

66 atomtypes = [] 

67 atomtypes_alt = [] 

68 if fpath.suffix == '.gz': 

69 import gzip 

70 opener = gzip.open 

71 elif fpath.suffix == '.bz2': 

72 import bz2 

73 opener = bz2.BZ2File 

74 else: 

75 opener = open 

76 with opener(fpath) as fd: 

77 for line in fd: 

78 if 'TITEL' in line: 

79 atomtypes.append(line.split()[3].split('_')[0].split('.')[0]) 

80 elif 'POTCAR:' in line: 

81 atomtypes_alt.append( 

82 line.split()[2].split('_')[0].split('.')[0]) 

83 

84 if len(atomtypes) == 0 and len(atomtypes_alt) > 0: 

85 # old VASP doesn't echo TITEL, but all versions print out species 

86 # lines preceded by "POTCAR:", twice 

87 if len(atomtypes_alt) % 2 != 0: 

88 raise ParseError( 

89 f'Tried to get atom types from {len(atomtypes_alt)}' 

90 '"POTCAR": lines in OUTCAR, but expected an even number' 

91 ) 

92 atomtypes = atomtypes_alt[0:len(atomtypes_alt) // 2] 

93 

94 return atomtypes 

95 

96 

97def atomtypes_outpot(posfname, numsyms): 

98 """Try to retrieve chemical symbols from OUTCAR or POTCAR 

99 

100 If getting atomtypes from the first line in POSCAR/CONTCAR fails, it might 

101 be possible to find the data in OUTCAR or POTCAR, if these files exist. 

102 

103 posfname -- The filename of the POSCAR/CONTCAR file we're trying to read 

104 

105 numsyms -- The number of symbols we must find 

106 

107 """ 

108 posfpath = Path(posfname) 

109 

110 # Check files with exactly same path except POTCAR/OUTCAR instead 

111 # of POSCAR/CONTCAR. 

112 fnames = [posfpath.with_name('POTCAR'), 

113 posfpath.with_name('OUTCAR')] 

114 # Try the same but with compressed files 

115 fsc = [] 

116 for fnpath in fnames: 

117 fsc.append(fnpath.parent / (fnpath.name + '.gz')) 

118 fsc.append(fnpath.parent / (fnpath.name + '.bz2')) 

119 for f in fsc: 

120 fnames.append(f) 

121 # Code used to try anything with POTCAR or OUTCAR in the name 

122 # but this is no longer supported 

123 

124 tried = [] 

125 for fn in fnames: 

126 if fn in posfpath.parent.iterdir(): 

127 tried.append(fn) 

128 at = get_atomtypes(fn) 

129 if len(at) == numsyms: 

130 return at 

131 

132 raise ParseError('Could not determine chemical symbols. Tried files ' + 

133 str(tried)) 

134 

135 

136def get_atomtypes_from_formula(formula): 

137 """Return atom types from chemical formula (optionally prepended 

138 with and underscore). 

139 """ 

140 from ase.symbols import string2symbols 

141 symbols = string2symbols(formula.split('_')[0]) 

142 atomtypes = [symbols[0]] 

143 for s in symbols[1:]: 

144 if s != atomtypes[-1]: 

145 atomtypes.append(s) 

146 return atomtypes 

147 

148 

149@reader 

150def read_vasp(fd): 

151 """Import POSCAR/CONTCAR type file. 

152 

153 Reads unitcell, atom positions and constraints from the POSCAR/CONTCAR 

154 file and tries to read atom types from POSCAR/CONTCAR header, if this 

155 fails the atom types are read from OUTCAR or POTCAR file. 

156 """ 

157 atoms = read_vasp_configuration(fd) 

158 velocity_init_line = fd.readline() 

159 if velocity_init_line.strip() and velocity_init_line[0].lower() == 'l': 

160 read_lattice_velocities(fd) 

161 velocities = read_velocities_if_present(fd, len(atoms)) 

162 if velocities is not None: 

163 atoms.set_velocities(velocities) 

164 return atoms 

165 

166 

167def read_vasp_configuration(fd): 

168 """Read common POSCAR/CONTCAR/CHGCAR/CHG quantities and return Atoms.""" 

169 from ase.data import chemical_symbols 

170 

171 # The first line is in principle a comment line, however in VASP 

172 # 4.x a common convention is to have it contain the atom symbols, 

173 # eg. "Ag Ge" in the same order as later in the file (and POTCAR 

174 # for the full vasp run). In the VASP 5.x format this information 

175 # is found on the fifth line. Thus we save the first line and use 

176 # it in case we later detect that we're reading a VASP 4.x format 

177 # file. 

178 line1 = fd.readline() 

179 

180 scale = parse_poscar_scaling_factor(fd.readline()) 

181 

182 # Now the lattice vectors 

183 cell = np.array([fd.readline().split()[:3] for _ in range(3)], dtype=float) 

184 # Negative scaling factor corresponds to the cell volume. 

185 if scale[0] < 0.0: 

186 scale = np.cbrt(-1.0 * scale / np.linalg.det(cell)) 

187 cell *= scale # This works for both one and three scaling factors. 

188 

189 # Number of atoms. Again this must be in the same order as 

190 # in the first line 

191 # or in the POTCAR or OUTCAR file 

192 atom_symbols = [] 

193 numofatoms = fd.readline().split() 

194 # Check whether we have a VASP 4.x or 5.x format file. If the 

195 # format is 5.x, use the fifth line to provide information about 

196 # the atomic symbols. 

197 vasp5 = False 

198 try: 

199 int(numofatoms[0]) 

200 except ValueError: 

201 vasp5 = True 

202 atomtypes = numofatoms 

203 numofatoms = fd.readline().split() 

204 

205 # check for comments in numofatoms line and get rid of them if necessary 

206 commentcheck = np.array(['!' in s for s in numofatoms]) 

207 if commentcheck.any(): 

208 # only keep the elements up to the first including a '!': 

209 numofatoms = numofatoms[:np.arange(len(numofatoms))[commentcheck][0]] 

210 

211 if not vasp5: 

212 # Split the comment line (first in the file) into words and 

213 # try to compose a list of chemical symbols 

214 from ase.formula import Formula 

215 atomtypes = [] 

216 for word in line1.split(): 

217 word_without_delims = re.sub(r"-|_|,|\.|=|[0-9]|^", "", word) 

218 if len(word_without_delims) < 1: 

219 continue 

220 try: 

221 atomtypes.extend(list(Formula(word_without_delims))) 

222 except ValueError: 

223 # print(atomtype, e, 'is comment') 

224 pass 

225 # Now the list of chemical symbols atomtypes must be formed. 

226 # For example: atomtypes = ['Pd', 'C', 'O'] 

227 

228 numsyms = len(numofatoms) 

229 if len(atomtypes) < numsyms: 

230 # First line in POSCAR/CONTCAR didn't contain enough symbols. 

231 

232 # Sometimes the first line in POSCAR/CONTCAR is of the form 

233 # "CoP3_In-3.pos". Check for this case and extract atom types 

234 if len(atomtypes) == 1 and '_' in atomtypes[0]: 

235 atomtypes = get_atomtypes_from_formula(atomtypes[0]) 

236 else: 

237 atomtypes = atomtypes_outpot(fd.name, numsyms) 

238 else: 

239 try: 

240 for atype in atomtypes[:numsyms]: 

241 if atype not in chemical_symbols: 

242 raise KeyError 

243 except KeyError: 

244 atomtypes = atomtypes_outpot(fd.name, numsyms) 

245 

246 # Certain versions of VASP compiled with hdf5 support will write POTCAR 

247 # symbols with hashes to the CONTCAR (e.g. "Na_pv/6a2f546d)". 

248 for i, atomtype in enumerate(atomtypes): 

249 potcar_label = atomtype.split('/')[0] 

250 element_symbol = potcar_label.split('_')[0] 

251 atomtypes[i] = element_symbol 

252 

253 for i, num in enumerate(numofatoms): 

254 numofatoms[i] = int(num) 

255 atom_symbols.extend(numofatoms[i] * [atomtypes[i]]) 

256 

257 # Check if Selective dynamics is switched on 

258 sdyn = fd.readline() 

259 selective_dynamics = sdyn[0].lower() == 's' 

260 

261 # Check if atom coordinates are cartesian or direct 

262 if selective_dynamics: 

263 ac_type = fd.readline() 

264 else: 

265 ac_type = sdyn 

266 cartesian = ac_type[0].lower() in ['c', 'k'] 

267 tot_natoms = sum(numofatoms) 

268 atoms_pos = np.empty((tot_natoms, 3)) 

269 if selective_dynamics: 

270 selective_flags = np.empty((tot_natoms, 3), dtype=bool) 

271 for atom in range(tot_natoms): 

272 ac = fd.readline().split() 

273 atoms_pos[atom] = [float(_) for _ in ac[0:3]] 

274 if selective_dynamics: 

275 selective_flags[atom] = [_ == 'F' for _ in ac[3:6]] 

276 

277 atoms = Atoms(symbols=atom_symbols, cell=cell, pbc=True) 

278 if cartesian: 

279 atoms_pos *= scale 

280 atoms.set_positions(atoms_pos) 

281 else: 

282 atoms.set_scaled_positions(atoms_pos) 

283 if selective_dynamics: 

284 set_constraints(atoms, selective_flags) 

285 

286 return atoms 

287 

288 

289def read_lattice_velocities(fd): 

290 """ 

291 Read lattice velocities and vectors from POSCAR/CONTCAR. 

292 As lattice velocities are not yet implemented in ASE, this function just 

293 throws away these lines. 

294 """ 

295 fd.readline() # initialization state 

296 for _ in range(3): # lattice velocities 

297 fd.readline() 

298 for _ in range(3): # lattice vectors 

299 fd.readline() 

300 fd.readline() # get rid of 1 empty line below if it exists 

301 

302 

303def read_velocities_if_present(fd, natoms) -> np.ndarray | None: 

304 """Read velocities from POSCAR/CONTCAR if present, return in ASE units.""" 

305 # Check if it is the velocities block or the MD extra block 

306 words = fd.readline().split() 

307 if len(words) <= 1: # MD extra block or end of file 

308 return None 

309 atoms_vel = np.empty((natoms, 3)) 

310 atoms_vel[0] = (float(words[0]), float(words[1]), float(words[2])) 

311 for atom in range(1, natoms): 

312 words = fd.readline().split() 

313 assert len(words) == 3 

314 atoms_vel[atom] = (float(words[0]), float(words[1]), float(words[2])) 

315 

316 # unit conversion from Angstrom/fs to ASE units 

317 return atoms_vel * (Ang / fs) 

318 

319 

320def set_constraints(atoms: Atoms, selective_flags: np.ndarray): 

321 """Set constraints based on selective_flags""" 

322 from ase.constraints import FixAtoms, FixConstraint, FixScaled 

323 

324 constraints: list[FixConstraint] = [] 

325 indices = [] 

326 for ind, sflags in enumerate(selective_flags): 

327 if sflags.any() and not sflags.all(): 

328 constraints.append(FixScaled(ind, sflags, atoms.get_cell())) 

329 elif sflags.all(): 

330 indices.append(ind) 

331 if indices: 

332 constraints.append(FixAtoms(indices)) 

333 if constraints: 

334 atoms.set_constraint(constraints) 

335 

336 

337def iread_vasp_out(filename, index=-1): 

338 """Import OUTCAR type file, as a generator.""" 

339 it = ImageIterator(vop.outcarchunks) 

340 return it(filename, index=index) 

341 

342 

343@reader 

344def read_vasp_out(filename='OUTCAR', index=-1): 

345 """Import OUTCAR type file. 

346 

347 Reads unitcell, atom positions, energies, and forces from the OUTCAR file 

348 and attempts to read constraints (if any) from CONTCAR/POSCAR, if present. 

349 """ 

350 # "filename" is actually a file-descriptor thanks to @reader 

351 g = iread_vasp_out(filename, index=index) 

352 # Code borrowed from formats.py:read 

353 if isinstance(index, (slice, str)): 

354 # Return list of atoms 

355 return list(g) 

356 else: 

357 # Return single atoms object 

358 return next(g) 

359 

360 

361@reader 

362def read_vasp_xdatcar(filename='XDATCAR', index=-1): 

363 """Import XDATCAR file. 

364 

365 Parameters 

366 ---------- 

367 index : int or slice or str 

368 Which frame(s) to read. The default is -1 (last frame). 

369 See :func:`ase.io.read` for details. 

370 

371 Notes 

372 ----- 

373 Constraints ARE NOT stored in the XDATCAR, and as such, Atoms objects 

374 retrieved from the XDATCAR will not have constraints. 

375 """ 

376 fd = filename # @reader decorator ensures this is a file descriptor 

377 images = [] 

378 

379 cell = np.eye(3) 

380 atomic_formula = '' 

381 

382 while True: 

383 comment_line = fd.readline() 

384 if "Direct configuration=" not in comment_line: 

385 try: 

386 lattice_constant = float(fd.readline()) 

387 except Exception: 

388 # XXX: When would this happen? 

389 break 

390 

391 xx = [float(x) for x in fd.readline().split()] 

392 yy = [float(y) for y in fd.readline().split()] 

393 zz = [float(z) for z in fd.readline().split()] 

394 cell = np.array([xx, yy, zz]) * lattice_constant 

395 

396 symbols = fd.readline().split() 

397 numbers = [int(n) for n in fd.readline().split()] 

398 total = sum(numbers) 

399 

400 atomic_formula = ''.join(f'{sym:s}{numbers[n]:d}' 

401 for n, sym in enumerate(symbols)) 

402 

403 fd.readline() 

404 

405 coords = [np.array(fd.readline().split(), float) for _ in range(total)] 

406 

407 image = Atoms(atomic_formula, cell=cell, pbc=True) 

408 image.set_scaled_positions(np.array(coords)) 

409 images.append(image) 

410 

411 if index is None: 

412 index = -1 

413 

414 if isinstance(index, str): 

415 index = string2index(index) 

416 

417 return images[index] 

418 

419 

420def __get_xml_parameter(par): 

421 """An auxiliary function that enables convenient extraction of 

422 parameter values from a vasprun.xml file with proper type 

423 handling. 

424 

425 """ 

426 def to_bool(b): 

427 if b == 'T': 

428 return True 

429 else: 

430 return False 

431 

432 to_type = {'int': int, 'logical': to_bool, 'string': str, 'float': float} 

433 

434 text = par.text 

435 if text is None: 

436 text = '' 

437 

438 # Float parameters do not have a 'type' attrib 

439 var_type = to_type[par.attrib.get('type', 'float')] 

440 

441 try: 

442 if par.tag == 'v': 

443 return list(map(var_type, text.split())) 

444 else: 

445 return var_type(text.strip()) 

446 except ValueError: 

447 # Vasp can sometimes write "*****" due to overflow 

448 return None 

449 

450 

451def read_vasp_xml(filename='vasprun.xml', index=-1): 

452 """Parse vasprun.xml file. 

453 

454 Reads unit cell, atom positions, energies, forces, and constraints 

455 from vasprun.xml file 

456 

457 Examples 

458 -------- 

459 >>> import ase.io 

460 >>> ase.io.write("out.traj", ase.io.read("vasprun.xml", index=":")) 

461 """ 

462 

463 import xml.etree.ElementTree as ET 

464 from collections import OrderedDict 

465 

466 tree = ET.iterparse(filename, events=['start', 'end']) 

467 

468 atoms_init = None 

469 calculation = [] 

470 ibz_kpts = None 

471 kpt_weights = None 

472 parameters = OrderedDict() 

473 

474 try: 

475 for event, elem in tree: 

476 

477 if event == 'end': 

478 if elem.tag == 'kpoints': 

479 for subelem in elem.iter(tag='generation'): 

480 kpts_params = OrderedDict() 

481 parameters['kpoints_generation'] = kpts_params 

482 for par in subelem.iter(): 

483 if par.tag in ['v', 'i'] and "name" in par.attrib: 

484 parname = par.attrib['name'].lower() 

485 kpts_params[parname] = __get_xml_parameter(par) 

486 

487 kpts = elem.findall("varray[@name='kpointlist']/v") 

488 ibz_kpts = np.zeros((len(kpts), 3)) 

489 

490 for i, kpt in enumerate(kpts): 

491 ibz_kpts[i] = [float(val) for val in kpt.text.split()] 

492 

493 kpt_weights = elem.findall('varray[@name="weights"]/v') 

494 kpt_weights = [float(val.text) for val in kpt_weights] 

495 

496 elif elem.tag == 'parameters': 

497 for par in elem.iter(): 

498 if par.tag in ['v', 'i']: 

499 parname = par.attrib['name'].lower() 

500 parameters[parname] = __get_xml_parameter(par) 

501 

502 elif elem.tag == 'atominfo': 

503 species = [] 

504 

505 for entry in elem.find("array[@name='atoms']/set"): 

506 species.append(entry[0].text.strip()) 

507 

508 natoms = len(species) 

509 

510 elif (elem.tag == 'structure' 

511 and elem.attrib.get('name') == 'initialpos'): 

512 cell_init = np.zeros((3, 3), dtype=float) 

513 

514 for i, v in enumerate( 

515 elem.find("crystal/varray[@name='basis']")): 

516 cell_init[i] = np.array( 

517 [float(val) for val in v.text.split()]) 

518 

519 scpos_init = np.zeros((natoms, 3), dtype=float) 

520 

521 for i, v in enumerate( 

522 elem.find("varray[@name='positions']")): 

523 scpos_init[i] = np.array( 

524 [float(val) for val in v.text.split()]) 

525 

526 constraints = [] 

527 fixed_indices = [] 

528 

529 for i, entry in enumerate( 

530 elem.findall("varray[@name='selective']/v")): 

531 flags = (np.array( 

532 entry.text.split() == np.array(['F', 'F', 'F']))) 

533 if flags.all(): 

534 fixed_indices.append(i) 

535 elif flags.any(): 

536 constraints.append(FixScaled(i, flags, cell_init)) 

537 

538 if fixed_indices: 

539 constraints.append(FixAtoms(fixed_indices)) 

540 

541 atoms_init = Atoms(species, 

542 cell=cell_init, 

543 scaled_positions=scpos_init, 

544 constraint=constraints, 

545 pbc=True) 

546 

547 elif event == 'start' and elem.tag == 'calculation': 

548 calculation.append(elem) 

549 

550 except ET.ParseError as parse_error: 

551 if atoms_init is None: 

552 raise parse_error 

553 if calculation and calculation[-1].find("energy") is None: 

554 calculation = calculation[:-1] 

555 if not calculation: 

556 yield atoms_init 

557 

558 if calculation: 

559 if isinstance(index, int): 

560 steps = [calculation[index]] 

561 else: 

562 steps = calculation[index] 

563 else: 

564 steps = [] 

565 

566 for step in steps: 

567 yield atoms_from_step(step, ibz_kpts=ibz_kpts, kpt_weights=kpt_weights, 

568 parameters=parameters, atoms=atoms_init.copy(), 

569 natoms=natoms) 

570 

571 

572def atoms_from_step(step, ibz_kpts, kpt_weights, parameters, atoms, natoms): 

573 from ase.calculators.singlepoint import ( 

574 SinglePointDFTCalculator, 

575 SinglePointKPoint, 

576 ) 

577 from ase.units import GPa 

578 

579 cell = np.zeros((3, 3), dtype=float) 

580 for i, vector in enumerate( 

581 step.find('structure/crystal/varray[@name="basis"]')): 

582 cell[i] = np.array([float(val) for val in vector.text.split()]) 

583 volume = np.linalg.det(cell) 

584 

585 free_energy = float(step.find('energy/i[@name="e_fr_energy"]').text) 

586 

587 # https://gitlab.com/ase/ase/-/merge_requests/2685 

588 # e_fr_energy in calculation/energy is actually an enthalpy including 

589 # the PV term, unlike that in /calculation/scstep/energy or in OUTCAR, 

590 # and therefore we need to subtract the PV term. 

591 pressure = parameters.get('pstress', 0.0) 

592 pressure *= 1e-22 / EVTOJ # kbar -> eV/A3 

593 free_energy -= pressure * volume 

594 

595 # Workaround for VASP bug, e_0_energy contains the wrong value 

596 # in calculation/energy, but calculation/scstep/energy does not 

597 # include classical VDW corrections. So, first calculate 

598 # e_0_energy - e_fr_energy from calculation/scstep/energy, then 

599 # apply that correction to e_fr_energy from calculation/energy. 

600 lastscf = step.findall('scstep/energy')[-1] 

601 dipoles = step.findall('scstep/dipole') 

602 if dipoles: 

603 lastdipole = dipoles[-1] 

604 else: 

605 lastdipole = None 

606 

607 de = (float(lastscf.find('i[@name="e_0_energy"]').text) - 

608 float(lastscf.find('i[@name="e_fr_energy"]').text)) 

609 

610 energy = free_energy + de 

611 

612 scpos = np.zeros((natoms, 3), dtype=float) 

613 for i, vector in enumerate( 

614 step.find('structure/varray[@name="positions"]')): 

615 scpos[i] = np.array([float(val) for val in vector.text.split()]) 

616 

617 forces = None 

618 fblocks = step.find('varray[@name="forces"]') 

619 if fblocks is not None: 

620 forces = np.zeros((natoms, 3), dtype=float) 

621 for i, vector in enumerate(fblocks): 

622 forces[i] = np.array( 

623 [float(val) for val in vector.text.split()]) 

624 

625 stress = None 

626 sblocks = step.find('varray[@name="stress"]') 

627 if sblocks is not None: 

628 stress = np.zeros((3, 3), dtype=float) 

629 for i, vector in enumerate(sblocks): 

630 stress[i] = np.array( 

631 [float(val) for val in vector.text.split()]) 

632 stress *= -0.1 * GPa 

633 stress = stress.reshape(9)[[0, 4, 8, 5, 2, 1]] 

634 

635 dipole = None 

636 if lastdipole is not None: 

637 dblock = lastdipole.find('v[@name="dipole"]') 

638 if dblock is not None: 

639 dipole = np.zeros((1, 3), dtype=float) 

640 dipole = np.array([float(val) for val in dblock.text.split()]) 

641 

642 dblock = step.find('dipole/v[@name="dipole"]') 

643 if dblock is not None: 

644 dipole = np.zeros((1, 3), dtype=float) 

645 dipole = np.array([float(val) for val in dblock.text.split()]) 

646 

647 efermi = step.find('dos/i[@name="efermi"]') 

648 if efermi is not None: 

649 efermi = float(efermi.text) 

650 

651 kpoints = [] 

652 for ikpt in range(1, len(ibz_kpts) + 1): 

653 kblocks = step.findall( 

654 'eigenvalues/array/set/set/set[@comment="kpoint %d"]' % ikpt) 

655 if kblocks is not None: 

656 for spin, kpoint in enumerate(kblocks): 

657 eigenvals = kpoint.findall('r') 

658 eps_n = np.zeros(len(eigenvals)) 

659 f_n = np.zeros(len(eigenvals)) 

660 for j, val in enumerate(eigenvals): 

661 val = val.text.split() 

662 eps_n[j] = float(val[0]) 

663 f_n[j] = float(val[1]) 

664 if len(kblocks) == 1: 

665 f_n *= 2 

666 kpoints.append( 

667 SinglePointKPoint(kpt_weights[ikpt - 1], spin, ikpt, 

668 eps_n, f_n)) 

669 if len(kpoints) == 0: 

670 kpoints = None 

671 

672 # DFPT properties 

673 # dielectric tensor 

674 dielectric_tensor = None 

675 sblocks = step.find('varray[@name="dielectric_dft"]') 

676 if sblocks is not None: 

677 dielectric_tensor = np.zeros((3, 3), dtype=float) 

678 for ii, vector in enumerate(sblocks): 

679 dielectric_tensor[ii] = np.fromstring(vector.text, sep=' ') 

680 

681 # Born effective charges 

682 born_charges = None 

683 fblocks = step.find('array[@name="born_charges"]') 

684 if fblocks is not None: 

685 born_charges = np.zeros((natoms, 3, 3), dtype=float) 

686 for ii, block in enumerate(fblocks[1:]): # 1. element = dimension 

687 for jj, vector in enumerate(block): 

688 born_charges[ii, jj] = np.fromstring(vector.text, sep=' ') 

689 

690 atoms.set_cell(cell) 

691 atoms.set_scaled_positions(scpos) 

692 atoms.calc = SinglePointDFTCalculator( 

693 atoms, 

694 energy=energy, 

695 forces=forces, 

696 stress=stress, 

697 free_energy=free_energy, 

698 ibzkpts=ibz_kpts, 

699 efermi=efermi, 

700 dipole=dipole, 

701 dielectric_tensor=dielectric_tensor, 

702 born_effective_charges=born_charges 

703 ) 

704 atoms.calc.name = 'vasp' 

705 atoms.calc.kpts = kpoints 

706 atoms.calc.parameters = parameters 

707 return atoms 

708 

709 

710@writer 

711def write_vasp_xdatcar(fd, images, label=None): 

712 """Write VASP MD trajectory (XDATCAR) file 

713 

714 Only Vasp 5 format is supported (for consistency with read_vasp_xdatcar) 

715 

716 Args: 

717 fd (str, fp): Output file 

718 images (iterable of Atoms): Atoms images to write. These must have 

719 consistent atom order and lattice vectors - this will not be 

720 checked. 

721 label (str): Text for first line of file. If empty, default to list 

722 of elements. 

723 

724 """ 

725 

726 images = iter(images) 

727 image = next(images) 

728 

729 if not isinstance(image, Atoms): 

730 raise TypeError("images should be a sequence of Atoms objects.") 

731 

732 _write_xdatcar_header(fd, image, label) 

733 _write_xdatcar_config(fd, image, index=1) 

734 for i, image in enumerate(images): 

735 _write_xdatcar_header(fd, image, label) 

736 # Index is off by 2: 1-indexed file vs 0-indexed Python; 

737 # and we already wrote the first block. 

738 _write_xdatcar_config(fd, image, i + 2) 

739 

740 

741def _write_xdatcar_header(fd, atoms, label): 

742 symbol_count = _symbol_count_from_symbols(atoms.get_chemical_symbols()) 

743 

744 if label is None: 

745 label = ' '.join([s for s, _ in symbol_count]) 

746 fd.write(label + '\n') 

747 

748 # Not using lattice constants, set it to 1 

749 fd.write(' 1\n') 

750 

751 float_string = '{:11.6f}' 

752 for row_i in range(3): 

753 fd.write(' ') 

754 fd.write(' '.join(float_string.format(x) for x in atoms.cell[row_i])) 

755 fd.write('\n') 

756 

757 fd.write(_symbol_count_string(symbol_count, vasp5=True)) 

758 

759 

760def _write_xdatcar_config(fd, atoms, index): 

761 """Write a block of positions for XDATCAR file 

762 

763 Args: 

764 fd (fd): writeable Python file descriptor 

765 atoms (ase.Atoms): Atoms to write 

766 index (int): configuration number written to block header 

767 

768 """ 

769 fd.write(f"Direct configuration={index:6d}\n") 

770 float_string = '{:11.8f}' 

771 scaled_positions = atoms.get_scaled_positions(wrap=False) 

772 for row in scaled_positions: 

773 fd.write(' ') 

774 fd.write(' '.join([float_string.format(x) for x in row])) 

775 fd.write('\n') 

776 

777 

778def _symbol_count_from_symbols(symbols: Symbols) -> list[tuple[str, int]]: 

779 """Reduce list of chemical symbols into compact VASP notation 

780 

781 Args: 

782 symbols (iterable of str) 

783