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 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 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 -------- 

452 >>> import ase.io 

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

454 """ 

455 

456 import xml.etree.ElementTree as ET 

457 from collections import OrderedDict 

458 

459 from ase.calculators.singlepoint import ( 

460 SinglePointDFTCalculator, 

461 SinglePointKPoint, 

462 ) 

463 from ase.units import GPa 

464 

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

466 

467 atoms_init = None 

468 calculation = [] 

469 ibz_kpts = None 

470 kpt_weights = None 

471 parameters = OrderedDict() 

472 

473 try: 

474 for event, elem in tree: 

475 

476 if event == 'end': 

477 if elem.tag == 'kpoints': 

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

479 kpts_params = OrderedDict() 

480 parameters['kpoints_generation'] = kpts_params 

481 for par in subelem.iter(): 

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

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

484 kpts_params[parname] = __get_xml_parameter(par) 

485 

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

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

488 

489 for i, kpt in enumerate(kpts): 

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

491 

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

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

494 

495 elif elem.tag == 'parameters': 

496 for par in elem.iter(): 

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

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

499 parameters[parname] = __get_xml_parameter(par) 

500 

501 elif elem.tag == 'atominfo': 

502 species = [] 

503 

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

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

506 

507 natoms = len(species) 

508 

509 elif (elem.tag == 'structure' 

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

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

512 

513 for i, v in enumerate( 

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

515 cell_init[i] = np.array( 

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

517 

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

519 

520 for i, v in enumerate( 

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

522 scpos_init[i] = np.array( 

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

524 

525 constraints = [] 

526 fixed_indices = [] 

527 

528 for i, entry in enumerate( 

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

530 flags = (np.array( 

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

532 if flags.all(): 

533 fixed_indices.append(i) 

534 elif flags.any(): 

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

536 

537 if fixed_indices: 

538 constraints.append(FixAtoms(fixed_indices)) 

539 

540 atoms_init = Atoms(species, 

541 cell=cell_init, 

542 scaled_positions=scpos_init, 

543 constraint=constraints, 

544 pbc=True) 

545 

546 elif elem.tag == 'dipole': 

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

548 if dblock is not None: 

549 dipole = np.array( 

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

551 

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

553 calculation.append(elem) 

554 

555 except ET.ParseError as parse_error: 

556 if atoms_init is None: 

557 raise parse_error 

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

559 calculation = calculation[:-1] 

560 if not calculation: 

561 yield atoms_init 

562 

563 if calculation: 

564 if isinstance(index, int): 

565 steps = [calculation[index]] 

566 else: 

567 steps = calculation[index] 

568 else: 

569 steps = [] 

570 

571 for step in steps: 

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

573 for i, vector in enumerate( 

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

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

576 volume = np.linalg.det(cell) 

577 

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

579 

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

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

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

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

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

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

586 free_energy -= pressure * volume 

587 

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

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

590 # include classical VDW corrections. So, first calculate 

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

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

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

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

595 if dipoles: 

596 lastdipole = dipoles[-1] 

597 else: 

598 lastdipole = None 

599 

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

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

602 

603 energy = free_energy + de 

604 

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

606 for i, vector in enumerate( 

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

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

609 

610 forces = None 

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

612 if fblocks is not None: 

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

614 for i, vector in enumerate(fblocks): 

615 forces[i] = np.array( 

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

617 

618 stress = None 

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

620 if sblocks is not None: 

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

622 for i, vector in enumerate(sblocks): 

623 stress[i] = np.array( 

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

625 stress *= -0.1 * GPa 

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

627 

628 dipole = None 

629 if lastdipole is not None: 

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

631 if dblock is not None: 

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

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

634 

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

636 if dblock is not None: 

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

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

639 

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

641 if efermi is not None: 

642 efermi = float(efermi.text) 

643 

644 kpoints = [] 

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

646 kblocks = step.findall( 

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

648 if kblocks is not None: 

649 for spin, kpoint in enumerate(kblocks): 

650 eigenvals = kpoint.findall('r') 

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

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

653 for j, val in enumerate(eigenvals): 

654 val = val.text.split() 

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

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

657 if len(kblocks) == 1: 

658 f_n *= 2 

659 kpoints.append( 

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

661 eps_n, f_n)) 

662 if len(kpoints) == 0: 

663 kpoints = None 

664 

665 # DFPT properties 

666 # dielectric tensor 

667 dielectric_tensor = None 

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

669 if sblocks is not None: 

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

671 for ii, vector in enumerate(sblocks): 

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

673 

674 # Born effective charges 

675 born_charges = None 

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

677 if fblocks is not None: 

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

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

680 for jj, vector in enumerate(block): 

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

682 

683 atoms = atoms_init.copy() 

684 atoms.set_cell(cell) 

685 atoms.set_scaled_positions(scpos) 

686 atoms.calc = SinglePointDFTCalculator( 

687 atoms, 

688 energy=energy, 

689 forces=forces, 

690 stress=stress, 

691 free_energy=free_energy, 

692 ibzkpts=ibz_kpts, 

693 efermi=efermi, 

694 dipole=dipole, 

695 dielectric_tensor=dielectric_tensor, 

696 born_effective_charges=born_charges 

697 ) 

698 atoms.calc.name = 'vasp' 

699 atoms.calc.kpts = kpoints 

700 atoms.calc.parameters = parameters 

701 yield atoms 

702 

703 

704@writer 

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

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

707 

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

709 

710 Args: 

711 fd (str, fp): Output file 

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

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

714 checked. 

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

716 of elements. 

717 

718 """ 

719 

720 images = iter(images) 

721 image = next(images) 

722 

723 if not isinstance(image, Atoms): 

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

725 

726 _write_xdatcar_header(fd, image, label) 

727 _write_xdatcar_config(fd, image, index=1) 

728 for i, image in enumerate(images): 

729 _write_xdatcar_header(fd, image, label) 

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

731 # and we already wrote the first block. 

732 _write_xdatcar_config(fd, image, i + 2) 

733 

734 

735def _write_xdatcar_header(fd, atoms, label): 

736 symbol_count = _symbol_count_from_symbols(atoms.get_chemical_symbols()) 

737 

738 if label is None: 

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

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

741 

742 # Not using lattice constants, set it to 1 

743 fd.write(' 1\n') 

744 

745 float_string = '{:11.6f}' 

746 for row_i in range(3): 

747 fd.write(' ') 

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

749 fd.write('\n') 

750 

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

752 

753 

754def _write_xdatcar_config(fd, atoms, index): 

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

756 

757 Args: 

758 fd (fd): writeable Python file descriptor 

759 atoms (ase.Atoms): Atoms to write 

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

761 

762 """ 

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

764 float_string = '{:11.8f}' 

765 scaled_positions = atoms.get_scaled_positions() 

766 for row in scaled_positions: 

767 fd.write(' ') 

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

769 fd.write('\n') 

770 

771 

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

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

774 

775 Args: 

776 symbols (iterable of str) 

777 

778 Returns 

779 ------- 

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

781 

782 Example: 

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

784 >>> _symbols_count_from_symbols(s) 

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

786 """ 

787 sc = [] 

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