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

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 List, Optional, TextIO, Tuple 

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 for i, num in enumerate(numofatoms): 

247 numofatoms[i] = int(num) 

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

249 

250 # Check if Selective dynamics is switched on 

251 sdyn = fd.readline() 

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

253 

254 # Check if atom coordinates are cartesian or direct 

255 if selective_dynamics: 

256 ac_type = fd.readline() 

257 else: 

258 ac_type = sdyn 

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

260 tot_natoms = sum(numofatoms) 

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

262 if selective_dynamics: 

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

264 for atom in range(tot_natoms): 

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

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

267 if selective_dynamics: 

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

269 

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

271 if cartesian: 

272 atoms_pos *= scale 

273 atoms.set_positions(atoms_pos) 

274 else: 

275 atoms.set_scaled_positions(atoms_pos) 

276 if selective_dynamics: 

277 set_constraints(atoms, selective_flags) 

278 

279 return atoms 

280 

281 

282def read_lattice_velocities(fd): 

283 """ 

284 Read lattice velocities and vectors from POSCAR/CONTCAR. 

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

286 throws away these lines. 

287 """ 

288 fd.readline() # initialization state 

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

290 fd.readline() 

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

292 fd.readline() 

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

294 

295 

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

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

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

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

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

301 return None 

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

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

304 for atom in range(1, natoms): 

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

306 assert len(words) == 3 

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

308 

309 # unit conversion from Angstrom/fs to ASE units 

310 return atoms_vel * (Ang / fs) 

311 

312 

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

314 """Set constraints based on selective_flags""" 

315 from ase.constraints import FixAtoms, FixConstraint, FixScaled 

316 

317 constraints: List[FixConstraint] = [] 

318 indices = [] 

319 for ind, sflags in enumerate(selective_flags): 

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

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

322 elif sflags.all(): 

323 indices.append(ind) 

324 if indices: 

325 constraints.append(FixAtoms(indices)) 

326 if constraints: 

327 atoms.set_constraint(constraints) 

328 

329 

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

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

332 it = ImageIterator(vop.outcarchunks) 

333 return it(filename, index=index) 

334 

335 

336@reader 

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

338 """Import OUTCAR type file. 

339 

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

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

342 """ 

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

344 g = iread_vasp_out(filename, index=index) 

345 # Code borrowed from formats.py:read 

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

347 # Return list of atoms 

348 return list(g) 

349 else: 

350 # Return single atoms object 

351 return next(g) 

352 

353 

354@reader 

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

356 """Import XDATCAR file. 

357 

358 Parameters 

359 ---------- 

360 index : int or slice or str 

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

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

363 

364 Notes 

365 ----- 

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

367 retrieved from the XDATCAR will not have constraints. 

368 """ 

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

370 images = [] 

371 

372 cell = np.eye(3) 

373 atomic_formula = '' 

374 

375 while True: 

376 comment_line = fd.readline() 

377 if "Direct configuration=" not in comment_line: 

378 try: 

379 lattice_constant = float(fd.readline()) 

380 except Exception: 

381 # XXX: When would this happen? 

382 break 

383 

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

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

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

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

388 

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

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

391 total = sum(numbers) 

392 

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

394 for n, sym in enumerate(symbols)) 

395 

396 fd.readline() 

397 

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

399 

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

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

402 images.append(image) 

403 

404 if index is None: 

405 index = -1 

406 

407 if isinstance(index, str): 

408 index = string2index(index) 

409 

410 return images[index] 

411 

412 

413def __get_xml_parameter(par): 

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

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

416 handling. 

417 

418 """ 

419 def to_bool(b): 

420 if b == 'T': 

421 return True 

422 else: 

423 return False 

424 

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

426 

427 text = par.text 

428 if text is None: 

429 text = '' 

430 

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

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

433 

434 try: 

435 if par.tag == 'v': 

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

437 else: 

438 return var_type(text.strip()) 

439 except ValueError: 

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

441 return None 

442 

443 

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

445 """Parse vasprun.xml file. 

446 

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

448 from vasprun.xml file 

449 

450 Examples: 

451 >>> import ase.io 

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

453 """ 

454 

455 import xml.etree.ElementTree as ET 

456 from collections import OrderedDict 

457 

458 from ase.calculators.singlepoint import ( 

459 SinglePointDFTCalculator, 

460 SinglePointKPoint, 

461 ) 

462 from ase.units import GPa 

463 

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

465 

466 atoms_init = None 

467 calculation = [] 

468 ibz_kpts = None 

469 kpt_weights = None 

470 parameters = OrderedDict() 

471 

472 try: 

473 for event, elem in tree: 

474 

475 if event == 'end': 

476 if elem.tag == 'kpoints': 

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

478 kpts_params = OrderedDict() 

479 parameters['kpoints_generation'] = kpts_params 

480 for par in subelem.iter(): 

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

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

483 kpts_params[parname] = __get_xml_parameter(par) 

484 

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

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

487 

488 for i, kpt in enumerate(kpts): 

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

490 

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

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

493 

494 elif elem.tag == 'parameters': 

495 for par in elem.iter(): 

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

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

498 parameters[parname] = __get_xml_parameter(par) 

499 

500 elif elem.tag == 'atominfo': 

501 species = [] 

502 

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

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

505 

506 natoms = len(species) 

507 

508 elif (elem.tag == 'structure' 

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

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

511 

512 for i, v in enumerate( 

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

514 cell_init[i] = np.array( 

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

516 

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

518 

519 for i, v in enumerate( 

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

521 scpos_init[i] = np.array( 

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

523 

524 constraints = [] 

525 fixed_indices = [] 

526 

527 for i, entry in enumerate( 

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

529 flags = (np.array( 

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

531 if flags.all(): 

532 fixed_indices.append(i) 

533 elif flags.any(): 

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

535 

536 if fixed_indices: 

537 constraints.append(FixAtoms(fixed_indices)) 

538 

539 atoms_init = Atoms(species, 

540 cell=cell_init, 

541 scaled_positions=scpos_init, 

542 constraint=constraints, 

543 pbc=True) 

544 

545 elif elem.tag == 'dipole': 

546 dblock = elem.find('v[@name="dipole"]') 

547 if dblock is not None: 

548 dipole = np.array( 

549 [float(val) for val in dblock.text.split()]) 

550 

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

552 calculation.append(elem) 

553 

554 except ET.ParseError as parse_error: 

555 if atoms_init is None: 

556 raise parse_error 

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

558 calculation = calculation[:-1] 

559 if not calculation: 

560 yield atoms_init 

561 

562 if calculation: 

563 if isinstance(index, int): 

564 steps = [calculation[index]] 

565 else: 

566 steps = calculation[index] 

567 else: 

568 steps = [] 

569 

570 for step in steps: 

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

572 for i, vector in enumerate( 

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

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

575 volume = np.linalg.det(cell) 

576 

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

578 

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

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

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

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

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

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

585 free_energy -= pressure * volume 

586 

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

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

589 # include classical VDW corrections. So, first calculate 

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

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

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

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

594 if dipoles: 

595 lastdipole = dipoles[-1] 

596 else: 

597 lastdipole = None 

598 

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

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

601 

602 energy = free_energy + de 

603 

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

605 for i, vector in enumerate( 

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

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

608 

609 forces = None 

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

611 if fblocks is not None: 

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

613 for i, vector in enumerate(fblocks): 

614 forces[i] = np.array( 

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

616 

617 stress = None 

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

619 if sblocks is not None: 

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

621 for i, vector in enumerate(sblocks): 

622 stress[i] = np.array( 

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

624 stress *= -0.1 * GPa 

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

626 

627 dipole = None 

628 if lastdipole is not None: 

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

630 if dblock is not None: 

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

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

633 

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

635 if dblock is not None: 

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

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

638 

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

640 if efermi is not None: 

641 efermi = float(efermi.text) 

642 

643 kpoints = [] 

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

645 kblocks = step.findall( 

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

647 if kblocks is not None: 

648 for spin, kpoint in enumerate(kblocks): 

649 eigenvals = kpoint.findall('r') 

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

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

652 for j, val in enumerate(eigenvals): 

653 val = val.text.split() 

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

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

656 if len(kblocks) == 1: 

657 f_n *= 2 

658 kpoints.append( 

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

660 eps_n, f_n)) 

661 if len(kpoints) == 0: 

662 kpoints = None 

663 

664 # DFPT properties 

665 # dielectric tensor 

666 dielectric_tensor = None 

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

668 if sblocks is not None: 

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

670 for ii, vector in enumerate(sblocks): 

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

672 

673 # Born effective charges 

674 born_charges = None 

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

676 if fblocks is not None: 

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

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

679 for jj, vector in enumerate(block): 

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

681 

682 atoms = atoms_init.copy() 

683 atoms.set_cell(cell) 

684 atoms.set_scaled_positions(scpos) 

685 atoms.calc = SinglePointDFTCalculator( 

686 atoms, 

687 energy=energy, 

688 forces=forces, 

689 stress=stress, 

690 free_energy=free_energy, 

691 ibzkpts=ibz_kpts, 

692 efermi=efermi, 

693 dipole=dipole, 

694 dielectric_tensor=dielectric_tensor, 

695 born_effective_charges=born_charges 

696 ) 

697 atoms.calc.name = 'vasp' 

698 atoms.calc.kpts = kpoints 

699 atoms.calc.parameters = parameters 

700 yield atoms 

701 

702 

703@writer 

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

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

706 

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

708 

709 Args: 

710 fd (str, fp): Output file 

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

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

713 checked. 

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

715 of elements. 

716 

717 """ 

718 

719 images = iter(images) 

720 image = next(images) 

721 

722 if not isinstance(image, Atoms): 

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

724 

725 _write_xdatcar_header(fd, image, label) 

726 _write_xdatcar_config(fd, image, index=1) 

727 for i, image in enumerate(images): 

728 _write_xdatcar_header(fd, image, label) 

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

730 # and we already wrote the first block. 

731 _write_xdatcar_config(fd, image, i + 2) 

732 

733 

734def _write_xdatcar_header(fd, atoms, label): 

735 symbol_count = _symbol_count_from_symbols(atoms.get_chemical_symbols()) 

736 

737 if label is None: 

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

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

740 

741 # Not using lattice constants, set it to 1 

742 fd.write(' 1\n') 

743 

744 float_string = '{:11.6f}' 

745 for row_i in range(3): 

746 fd.write(' ') 

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

748 fd.write('\n') 

749 

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

751 

752 

753def _write_xdatcar_config(fd, atoms, index): 

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

755 

756 Args: 

757 fd (fd): writeable Python file descriptor 

758 atoms (ase.Atoms): Atoms to write 

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

760 

761 """ 

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

763 float_string = '{:11.8f}' 

764 scaled_positions = atoms.get_scaled_positions() 

765 for row in scaled_positions: 

766 fd.write(' ') 

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

768 fd.write('\n') 

769 

770 

771def _symbol_count_from_symbols(symbols: Symbols) -> List[Tuple[str, int]]: 

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

773 

774 Args: 

775 symbols (iterable of str) 

776 

777 Returns: 

778 list of pairs [(el1, c1), (el2, c2), ...] 

779 

780 Example: 

781 >>> s = Atoms('Ar3NeHe2ArNe').symbols 

782 >>> _symbols_count_from_symbols(s) 

783 [('Ar', 3), ('Ne', 1), ('He', 2), ('Ar', 1), ('Ne', 1)] 

784 """ 

785 sc = [] 

786 psym = str(symbols[0]) # we cast to str to appease mypy