Coverage for /builds/debichem-team/python-ase/ase/io/aims.py: 93.06%
807 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-03-06 04:00 +0000
1"""Defines class/functions to write input and parse output for FHI-aims."""
2import os
3import re
4import time
5import warnings
6from functools import cached_property
7from pathlib import Path
8from typing import Any, Dict, List, Union
10import numpy as np
12from ase import Atom, Atoms
13from ase.calculators.calculator import kpts2mp
14from ase.calculators.singlepoint import SinglePointDFTCalculator
15from ase.constraints import FixAtoms, FixCartesian
16from ase.data import atomic_numbers
17from ase.io import ParseError
18from ase.units import Ang, fs
19from ase.utils import deprecated, reader, writer
21v_unit = Ang / (1000.0 * fs)
23LINE_NOT_FOUND = object()
26class AimsParseError(Exception):
27 """Exception raised if an error occurs when parsing an Aims output file"""
29 def __init__(self, message):
30 self.message = message
31 super().__init__(self.message)
34# Read aims geometry files
35@reader
36def read_aims(fd, apply_constraints=True):
37 """Import FHI-aims geometry type files.
39 Reads unitcell, atom positions and constraints from
40 a geometry.in file.
42 If geometric constraint (symmetry parameters) are in the file
43 include that information in atoms.info["symmetry_block"]
44 """
46 lines = fd.readlines()
47 return parse_geometry_lines(lines, apply_constraints=apply_constraints)
50def parse_geometry_lines(lines, apply_constraints=True):
52 from ase import Atoms
53 from ase.constraints import (
54 FixAtoms,
55 FixCartesian,
56 FixCartesianParametricRelations,
57 FixScaledParametricRelations,
58 )
60 atoms = Atoms()
62 positions = []
63 cell = []
64 symbols = []
65 velocities = []
66 magmoms = []
67 symmetry_block = []
68 charges = []
69 fix = []
70 fix_cart = []
71 xyz = np.array([0, 0, 0])
72 i = -1
73 n_periodic = -1
74 periodic = np.array([False, False, False])
75 cart_positions, scaled_positions = False, False
76 for line in lines:
77 inp = line.split()
78 if inp == []:
79 continue
80 if inp[0] in ["atom", "atom_frac"]:
82 if inp[0] == "atom":
83 cart_positions = True
84 else:
85 scaled_positions = True
87 if xyz.all():
88 fix.append(i)
89 elif xyz.any():
90 fix_cart.append(FixCartesian(i, xyz))
91 floatvect = float(inp[1]), float(inp[2]), float(inp[3])
92 positions.append(floatvect)
93 symbols.append(inp[4])
94 magmoms.append(0.0)
95 charges.append(0.0)
96 xyz = np.array([0, 0, 0])
97 i += 1
99 elif inp[0] == "lattice_vector":
100 floatvect = float(inp[1]), float(inp[2]), float(inp[3])
101 cell.append(floatvect)
102 n_periodic = n_periodic + 1
103 periodic[n_periodic] = True
105 elif inp[0] == "initial_moment":
106 magmoms[-1] = float(inp[1])
108 elif inp[0] == "initial_charge":
109 charges[-1] = float(inp[1])
111 elif inp[0] == "constrain_relaxation":
112 if inp[1] == ".true.":
113 fix.append(i)
114 elif inp[1] == "x":
115 xyz[0] = 1
116 elif inp[1] == "y":
117 xyz[1] = 1
118 elif inp[1] == "z":
119 xyz[2] = 1
121 elif inp[0] == "velocity":
122 floatvect = [v_unit * float(line) for line in inp[1:4]]
123 velocities.append(floatvect)
125 elif inp[0] in [
126 "symmetry_n_params",
127 "symmetry_params",
128 "symmetry_lv",
129 "symmetry_frac",
130 ]:
131 symmetry_block.append(" ".join(inp))
133 if xyz.all():
134 fix.append(i)
135 elif xyz.any():
136 fix_cart.append(FixCartesian(i, xyz))
138 if cart_positions and scaled_positions:
139 raise Exception(
140 "Can't specify atom positions with mixture of "
141 "Cartesian and fractional coordinates"
142 )
143 elif scaled_positions and periodic.any():
144 atoms = Atoms(
145 symbols,
146 scaled_positions=positions,
147 cell=cell,
148 pbc=periodic)
149 else:
150 atoms = Atoms(symbols, positions)
152 if len(velocities) > 0:
153 if len(velocities) != len(positions):
154 raise Exception(
155 "Number of positions and velocities have to coincide.")
156 atoms.set_velocities(velocities)
158 fix_params = []
160 if len(symmetry_block) > 5:
161 params = symmetry_block[1].split()[1:]
163 lattice_expressions = []
164 lattice_params = []
166 atomic_expressions = []
167 atomic_params = []
169 n_lat_param = int(symmetry_block[0].split(" ")[2])
171 lattice_params = params[:n_lat_param]
172 atomic_params = params[n_lat_param:]
174 for ll, line in enumerate(symmetry_block[2:]):
175 expression = " ".join(line.split(" ")[1:])
176 if ll < 3:
177 lattice_expressions += expression.split(",")
178 else:
179 atomic_expressions += expression.split(",")
181 fix_params.append(
182 FixCartesianParametricRelations.from_expressions(
183 list(range(3)),
184 lattice_params,
185 lattice_expressions,
186 use_cell=True,
187 )
188 )
190 fix_params.append(
191 FixScaledParametricRelations.from_expressions(
192 list(range(len(atoms))), atomic_params, atomic_expressions
193 )
194 )
196 if any(magmoms):
197 atoms.set_initial_magnetic_moments(magmoms)
198 if any(charges):
199 atoms.set_initial_charges(charges)
201 if periodic.any():
202 atoms.set_cell(cell)
203 atoms.set_pbc(periodic)
204 if len(fix):
205 atoms.set_constraint([FixAtoms(indices=fix)] + fix_cart + fix_params)
206 else:
207 atoms.set_constraint(fix_cart + fix_params)
209 if fix_params and apply_constraints:
210 atoms.set_positions(atoms.get_positions())
211 return atoms
214def get_aims_header():
215 """Returns the header for aims input files"""
216 lines = ["#" + "=" * 79]
217 for line in [
218 "Created using the Atomic Simulation Environment (ASE)",
219 time.asctime(),
220 ]:
221 lines.append("# " + line + "\n")
222 return lines
225def _write_velocities_alias(args: List, kwargs: Dict[str, Any]) -> bool:
226 arg_position = 5
227 if len(args) > arg_position and args[arg_position]:
228 args[arg_position - 1] = True
229 elif kwargs.get("velocities", False):
230 if len(args) < arg_position:
231 kwargs["write_velocities"] = True
232 else:
233 args[arg_position - 1] = True
234 else:
235 return False
236 return True
239# Write aims geometry files
240@deprecated(
241 "Use of `velocities` is deprecated, please use `write_velocities`",
242 category=FutureWarning,
243 callback=_write_velocities_alias,
244)
245@writer
246def write_aims(
247 fd,
248 atoms,
249 scaled=False,
250 geo_constrain=False,
251 write_velocities=False,
252 velocities=False,
253 ghosts=None,
254 info_str=None,
255 wrap=False,
256):
257 """Method to write FHI-aims geometry files.
259 Writes the atoms positions and constraints (only FixAtoms is
260 supported at the moment).
262 Args:
263 fd: file object
264 File to output structure to
265 atoms: ase.atoms.Atoms
266 structure to output to the file
267 scaled: bool
268 If True use fractional coordinates instead of Cartesian coordinates
269 symmetry_block: list of str
270 List of geometric constraints as defined in:
271 :arxiv:`1908.01610`
272 write_velocities: bool
273 If True add the atomic velocity vectors to the file
274 velocities: bool
275 NOT AN ARRAY OF VELOCITIES, but the legacy version of
276 `write_velocities`
277 ghosts: list of Atoms
278 A list of ghost atoms for the system
279 info_str: str
280 A string to be added to the header of the file
281 wrap: bool
282 Wrap atom positions to cell before writing
284 .. deprecated:: 3.23.0
285 Use of ``velocities`` is deprecated, please use ``write_velocities``.
286 """
288 if scaled and not np.all(atoms.pbc):
289 raise ValueError(
290 "Requesting scaled for a calculation where scaled=True, but "
291 "the system is not periodic")
293 if geo_constrain:
294 if not scaled and np.all(atoms.pbc):
295 warnings.warn(
296 "Setting scaled to True because a symmetry_block is detected."
297 )
298 scaled = True
299 elif not np.all(atoms.pbc):
300 warnings.warn(
301 "Parameteric constraints can only be used in periodic systems."
302 )
303 geo_constrain = False
305 for line in get_aims_header():
306 fd.write(line + "\n")
308 # If writing additional information is requested via info_str:
309 if info_str is not None:
310 fd.write("\n# Additional information:\n")
311 if isinstance(info_str, list):
312 fd.write("\n".join([f"# {s}" for s in info_str]))
313 else:
314 fd.write(f"# {info_str}")
315 fd.write("\n")
317 fd.write("#=======================================================\n")
319 i = 0
320 if atoms.get_pbc().any():
321 for n, vector in enumerate(atoms.get_cell()):
322 fd.write("lattice_vector ")
323 for i in range(3):
324 fd.write(f"{vector[i]:16.16f} ")
325 fd.write("\n")
327 fix_cart = np.zeros((len(atoms), 3), dtype=bool)
328 for constr in atoms.constraints:
329 if isinstance(constr, FixAtoms):
330 fix_cart[constr.index] = (True, True, True)
331 elif isinstance(constr, FixCartesian):
332 fix_cart[constr.index] = constr.mask
334 if ghosts is None:
335 ghosts = np.zeros(len(atoms))
336 else:
337 assert len(ghosts) == len(atoms)
339 wrap = wrap and not geo_constrain
340 scaled_positions = atoms.get_scaled_positions(wrap=wrap)
342 for i, atom in enumerate(atoms):
343 if ghosts[i] == 1:
344 atomstring = "empty "
345 elif scaled:
346 atomstring = "atom_frac "
347 else:
348 atomstring = "atom "
349 fd.write(atomstring)
350 if scaled:
351 for pos in scaled_positions[i]:
352 fd.write(f"{pos:16.16f} ")
353 else:
354 for pos in atom.position:
355 fd.write(f"{pos:16.16f} ")
356 fd.write(atom.symbol)
357 fd.write("\n")
358 # (1) all coords are constrained:
359 if fix_cart[i].all():
360 fd.write(" constrain_relaxation .true.\n")
361 # (2) some coords are constrained:
362 elif fix_cart[i].any():
363 xyz = fix_cart[i]
364 for n in range(3):
365 if xyz[n]:
366 fd.write(f" constrain_relaxation {'xyz'[n]}\n")
367 if atom.charge:
368 fd.write(f" initial_charge {atom.charge:16.6f}\n")
369 if atom.magmom:
370 fd.write(f" initial_moment {atom.magmom:16.6f}\n")
372 if write_velocities and atoms.get_velocities() is not None:
373 v = atoms.get_velocities()[i] / v_unit
374 fd.write(f" velocity {v[0]:.16f} {v[1]:.16f} {v[2]:.16f}\n")
376 if geo_constrain:
377 for line in get_sym_block(atoms):
378 fd.write(line)
381def get_sym_block(atoms):
382 """Get symmetry block for Parametric constraints in atoms.constraints"""
383 from ase.constraints import (
384 FixCartesianParametricRelations,
385 FixScaledParametricRelations,
386 )
388 # Initialize param/expressions lists
389 atomic_sym_params = []
390 lv_sym_params = []
391 atomic_param_constr = np.zeros((len(atoms),), dtype="<U100")
392 lv_param_constr = np.zeros((3,), dtype="<U100")
394 # Populate param/expressions list
395 for constr in atoms.constraints:
396 if isinstance(constr, FixScaledParametricRelations):
397 atomic_sym_params += constr.params
399 if np.any(atomic_param_constr[constr.indices] != ""):
400 warnings.warn(
401 "multiple parametric constraints defined for the same "
402 "atom, using the last one defined"
403 )
405 atomic_param_constr[constr.indices] = [
406 ", ".join(expression) for expression in constr.expressions
407 ]
408 elif isinstance(constr, FixCartesianParametricRelations):
409 lv_sym_params += constr.params
411 if np.any(lv_param_constr[constr.indices] != ""):
412 warnings.warn(
413 "multiple parametric constraints defined for the same "
414 "lattice vector, using the last one defined"
415 )
417 lv_param_constr[constr.indices] = [
418 ", ".join(expression) for expression in constr.expressions
419 ]
421 if np.all(atomic_param_constr == "") and np.all(lv_param_constr == ""):
422 return []
424 # Check Constraint Parameters
425 if len(atomic_sym_params) != len(np.unique(atomic_sym_params)):
426 warnings.warn(
427 "Some parameters were used across constraints, they will be "
428 "combined in the aims calculations"
429 )
430 atomic_sym_params = np.unique(atomic_sym_params)
432 if len(lv_sym_params) != len(np.unique(lv_sym_params)):
433 warnings.warn(
434 "Some parameters were used across constraints, they will be "
435 "combined in the aims calculations"
436 )
437 lv_sym_params = np.unique(lv_sym_params)
439 if np.any(atomic_param_constr == ""):
440 raise OSError(
441 "FHI-aims input files require all atoms have defined parametric "
442 "constraints"
443 )
445 cell_inds = np.where(lv_param_constr == "")[0]
446 for ind in cell_inds:
447 lv_param_constr[ind] = "{:.16f}, {:.16f}, {:.16f}".format(
448 *atoms.cell[ind])
450 n_atomic_params = len(atomic_sym_params)
451 n_lv_params = len(lv_sym_params)
452 n_total_params = n_atomic_params + n_lv_params
454 sym_block = []
455 if n_total_params > 0:
456 sym_block.append("#" + "=" * 55 + "\n")
457 sym_block.append("# Parametric constraints\n")
458 sym_block.append("#" + "=" * 55 + "\n")
459 sym_block.append(
460 "symmetry_n_params {:d} {:d} {:d}\n".format(
461 n_total_params, n_lv_params, n_atomic_params
462 )
463 )
464 sym_block.append(
465 "symmetry_params %s\n" % " ".join(lv_sym_params + atomic_sym_params)
466 )
468 for constr in lv_param_constr:
469 sym_block.append(f"symmetry_lv {constr:s}\n")
471 for constr in atomic_param_constr:
472 sym_block.append(f"symmetry_frac {constr:s}\n")
473 return sym_block
476def format_aims_control_parameter(key, value, format="%s"):
477 """Format a line for the aims control.in
479 Parameter
480 ---------
481 key: str
482 Name of the paramteter to format
483 value: Object
484 The value to pass to the parameter
485 format: str
486 string to format the the text as
488 Returns
489 -------
490 str
491 The properly formatted line for the aims control.in
492 """
493 return f"{key:35s}" + (format % value) + "\n"
496# Write aims control.in files
497@writer
498def write_control(fd, atoms, parameters, verbose_header=False):
499 """Write the control.in file for FHI-aims
500 Parameters
501 ----------
502 fd: str
503 The file object to write to
504 atoms: atoms.Atoms
505 The Atoms object for the requested calculation
506 parameters: dict
507 The dictionary of all paramters for the calculation
508 verbose_header: bool
509 If True then explcitly list the paramters used to generate the
510 control.in file inside the header
511 """
513 parameters = dict(parameters)
514 lim = "#" + "=" * 79
516 if parameters["xc"] == "LDA":
517 parameters["xc"] = "pw-lda"
519 cubes = parameters.pop("cubes", None)
521 for line in get_aims_header():
522 fd.write(line + "\n")
524 if verbose_header:
525 fd.write("# \n# List of parameters used to initialize the calculator:")
526 for p, v in parameters.items():
527 s = f"# {p}:{v}\n"
528 fd.write(s)
529 fd.write(lim + "\n")
531 assert "kpts" not in parameters or "k_grid" not in parameters
532 assert "smearing" not in parameters or "occupation_type" not in parameters
534 for key, value in parameters.items():
535 if key == "kpts":
536 mp = kpts2mp(atoms, parameters["kpts"])
537 dk = 0.5 - 0.5 / np.array(mp)
538 fd.write(
539 format_aims_control_parameter(
540 "k_grid",
541 tuple(mp),
542 "%d %d %d"))
543 fd.write(
544 format_aims_control_parameter(
545 "k_offset",
546 tuple(dk),
547 "%f %f %f"))
548 elif key in ("species_dir", "tier"):
549 continue
550 elif key == "aims_command":
551 continue
552 elif key == "plus_u":
553 continue
554 elif key == "smearing":
555 name = parameters["smearing"][0].lower()
556 if name == "fermi-dirac":
557 name = "fermi"
558 width = parameters["smearing"][1]
559 if name == "methfessel-paxton":
560 order = parameters["smearing"][2]
561 order = " %d" % order
562 else:
563 order = ""
565 fd.write(
566 format_aims_control_parameter(
567 "occupation_type", (name, width, order), "%s %f%s"
568 )
569 )
570 elif key == "output":
571 for output_type in value:
572 fd.write(format_aims_control_parameter(key, output_type, "%s"))
573 elif key == "vdw_correction_hirshfeld" and value:
574 fd.write(format_aims_control_parameter(key, "", "%s"))
575 elif isinstance(value, bool):
576 fd.write(
577 format_aims_control_parameter(
578 key, str(value).lower(), ".%s."))
579 elif isinstance(value, (tuple, list)):
580 fd.write(
581 format_aims_control_parameter(
582 key, " ".join([str(x) for x in value]), "%s"
583 )
584 )
585 elif isinstance(value, str):
586 fd.write(format_aims_control_parameter(key, value, "%s"))
587 else:
588 fd.write(format_aims_control_parameter(key, value, "%r"))
590 if cubes:
591 cubes.write(fd)
593 fd.write(lim + "\n\n")
595 # Get the species directory
596 species_dir = get_species_directory
597 # dicts are ordered as of python 3.7
598 species_array = np.array(list(dict.fromkeys(atoms.symbols)))
599 # Grab the tier specification from the parameters. THis may either
600 # be None, meaning the default should be used for all species, or a
601 # list of integers/None values giving a specific basis set size
602 # for each species in the calculation.
603 tier = parameters.pop("tier", None)
604 tier_array = np.full(len(species_array), tier)
605 # Path to species files for FHI-aims. In this files are specifications
606 # for the basis set sizes depending on which basis set tier is used.
607 species_dir = get_species_directory(parameters.get("species_dir"))
608 # Parse the species files for each species present in the calculation
609 # according to the tier of each species.
610 species_basis_dict = parse_species_path(
611 species_array=species_array, tier_array=tier_array,
612 species_dir=species_dir)
613 # Write the basis functions to be included for each species in the
614 # calculation into the control.in file (fd).
615 write_species(fd, species_basis_dict, parameters)
618def get_species_directory(species_dir=None):
619 """Get the directory where the basis set information is stored
621 If the requested directory does not exist then raise an Error
623 Parameters
624 ----------
625 species_dir: str
626 Requested directory to find the basis set info from. E.g.
627 `~/aims2022/FHIaims/species_defaults/defaults_2020/light`.
629 Returns
630 -------
631 Path
632 The Path to the requested or default species directory.
634 Raises
635 ------
636 RuntimeError
637 If both the requested directory and the default one is not defined
638 or does not exit.
639 """
640 if species_dir is None:
641 species_dir = os.environ.get("AIMS_SPECIES_DIR")
643 if species_dir is None:
644 raise RuntimeError(
645 "Missing species directory! Use species_dir "
646 + "parameter or set $AIMS_SPECIES_DIR environment variable."
647 )
649 species_path = Path(species_dir)
650 if not species_path.exists():
651 raise RuntimeError(
652 f"The requested species_dir {species_dir} does not exist")
654 return species_path
657def write_species(control_file_descriptor, species_basis_dict, parameters):
658 """Write species for the calculation depending on basis set size.
660 The calculation should include certain basis set size function depending
661 on the numerical settings (light, tight, really tight) and the basis set
662 size (minimal, tier1, tier2, tier3, tier4). If the basis set size is not
663 given then a 'standard' basis set size is used for each numerical setting.
664 The species files are defined according to these standard basis set sizes
665 for the numerical settings in the FHI-aims repository.
667 Note, for FHI-aims in ASE, we don't explicitly give the numerical setting.
668 Instead we include the numerical setting in the species path: e.g.
669 `~/aims2022/FHIaims/species_defaults/defaults_2020/light` this path has
670 `light`, the numerical setting, as the last folder in the path.
672 Example - a basis function might be commented in the standard basis set size
673 such as "# hydro 4 f 7.4" and this basis function should be
674 uncommented for another basis set size such as tier4.
676 Args:
677 control_file_descriptor: File descriptor for the control.in file into
678 which we need to write relevant basis functions to be included for
679 the calculation.
680 species_basis_dict: Dictionary where keys as the species symbols and
681 each value is a single string containing all the basis functions
682 to be included in the caclculation.
683 parameters: Calculation parameters as a dict.
684 """
685 # Now for every species (key) in the species_basis_dict, save the
686 # relevant basis functions (values) from the species_basis_dict, by
687 # writing to the file handle (species_file_descriptor) given to this
688 # function.
689 for species_symbol, basis_set_text in species_basis_dict.items():
690 control_file_descriptor.write(basis_set_text)
691 if parameters.get("plus_u") is not None:
692 if species_symbol in parameters.plus_u:
693 control_file_descriptor.write(
694 f"plus_u {parameters.plus_u[species_symbol]} \n")
697def parse_species_path(species_array, tier_array, species_dir):
698 """Parse the species files for each species according to the tier given.
700 Args:
701 species_array: An array of species/element symbols present in the unit
702 cell (e.g. ['C', 'H'].)
703 tier_array: An array of None/integer values which define which basis
704 set size to use for each species/element in the calcualtion.
705 species_dir: Directory containing FHI-aims species files.
707 Returns:
708 Dictionary containing species as keys and the basis set specification
709 for each species as text as the value for the key.
710 """
711 if len(species_array) != len(tier_array):
712 raise ValueError(
713 f"The species array length: {len(species_array)}, "
714 f"is not the same as the tier_array length: {len(tier_array)}")
716 species_basis_dict = {}
718 for symbol, tier in zip(species_array, tier_array):
719 path = species_dir / f"{atomic_numbers[symbol]:02}_{symbol}_default"
720 # Open the species file:
721 with open(path, encoding="utf8") as species_file_handle:
722 # Read the species file into a string.
723 species_file_str = species_file_handle.read()
724 species_basis_dict[symbol] = manipulate_tiers(
725 species_file_str, tier)
726 return species_basis_dict
729def manipulate_tiers(species_string: str, tier: Union[None, int] = 1):
730 """Adds basis set functions based on the tier value.
732 This function takes in the species file as a string, it then searches
733 for relevant basis functions based on the tier value to include in a new
734 string that is returned.
736 Args:
737 species_string: species file (default) for a given numerical setting
738 (light, tight, really tight) given as a string.
739 tier: The basis set size. This will dictate which basis set functions
740 are included in the returned string.
742 Returns:
743 Basis set functions defined by the tier as a string.
744 """
745 if tier is None: # Then we use the default species file untouched.
746 return species_string
747 tier_pattern = r"(# \".* tier\" .*|# +Further.*)"
748 top, *tiers = re.split(tier_pattern, species_string)
749 tier_comments = tiers[::2]
750 tier_basis = tiers[1::2]
751 assert len(
752 tier_comments) == len(tier_basis), "Something wrong with splitting"
753 n_tiers = len(tier_comments)
754 assert tier <= n_tiers, f"Only {n_tiers} tiers available, you choose {tier}"
755 string_new = top
756 for i, (c, b) in enumerate(zip(tier_comments, tier_basis)):
757 b = re.sub(r"\n( *for_aux| *hydro| *ionic| *confined)", r"\n#\g<1>", b)
758 if i < tier:
759 b = re.sub(
760 r"\n#( *for_aux| *hydro| *ionic| *confined)", r"\n\g<1>", b)
761 string_new += c + b
762 return string_new
765# Read aims.out files
766scalar_property_to_line_key = {
767 "free_energy": ["| Electronic free energy"],
768 "number_of_iterations": ["| Number of self-consistency cycles"],
769 "magnetic_moment": ["N_up - N_down"],
770 "n_atoms": ["| Number of atoms"],
771 "n_bands": [
772 "Number of Kohn-Sham states",
773 "Reducing total number of Kohn-Sham states",
774 "Reducing total number of Kohn-Sham states",
775 ],
776 "n_electrons": ["The structure contains"],
777 "n_kpts": ["| Number of k-points"],
778 "n_spins": ["| Number of spin channels"],
779 "electronic_temp": ["Occupation type:"],
780 "fermi_energy": ["| Chemical potential (Fermi level)"],
781}
784class AimsOutChunk:
785 """Base class for AimsOutChunks"""
787 def __init__(self, lines):
788 """Constructor
790 Parameters
791 ----------
792 lines: list of str
793 The set of lines from the output file the encompasses either
794 a single structure within a trajectory or
795 general information about the calculation (header)
796 """
797 self.lines = lines
799 def reverse_search_for(self, keys, line_start=0):
800 """Find the last time one of the keys appears in self.lines
802 Parameters
803 ----------
804 keys: list of str
805 The key strings to search for in self.lines
806 line_start: int
807 The lowest index to search for in self.lines
809 Returns
810 -------
811 int
812 The last time one of the keys appears in self.lines
813 """
814 for ll, line in enumerate(self.lines[line_start:][::-1]):
815 if any(key in line for key in keys):
816 return len(self.lines) - ll - 1
818 return LINE_NOT_FOUND
820 def search_for_all(self, key, line_start=0, line_end=-1):
821 """Find the all times the key appears in self.lines
823 Parameters
824 ----------
825 key: str
826 The key string to search for in self.lines
827 line_start: int
828 The first line to start the search from
829 line_end: int
830 The last line to end the search at
832 Returns
833 -------
834 list of ints
835 All times the key appears in the lines
836 """
837 line_index = []
838 for ll, line in enumerate(self.lines[line_start:line_end]):
839 if key in line:
840 line_index.append(ll + line_start)
841 return line_index
843 def parse_scalar(self, property):
844 """Parse a scalar property from the chunk
846 Parameters
847 ----------
848 property: str
849 The property key to parse
851 Returns
852 -------
853 float
854 The scalar value of the property
855 """
856 line_start = self.reverse_search_for(
857 scalar_property_to_line_key[property])
859 if line_start == LINE_NOT_FOUND:
860 return None
862 line = self.lines[line_start]
863 return float(line.split(":")[-1].strip().split()[0])
866class AimsOutHeaderChunk(AimsOutChunk):
867 """The header of the aims.out file containint general information"""
869 def __init__(self, lines):
870 """Constructor
872 Parameters
873 ----------
874 lines: list of str
875 The lines inside the aims.out header
876 """
877 super().__init__(lines)
879 @cached_property
880 def constraints(self):
881 """Parse the constraints from the aims.out file
883 Constraints for the lattice vectors are not supported.
884 """
886 line_inds = self.search_for_all("Found relaxation constraint for atom")
887 if len(line_inds) == 0:
888 return []
890 fix = []
891 fix_cart = []
892 for ll in line_inds:
893 line = self.lines[ll]
894 xyz = [0, 0, 0]
895 ind = int(line.split()[5][:-1]) - 1
896 if "All coordinates fixed" in line:
897 if ind not in fix:
898 fix.append(ind)
899 if "coordinate fixed" in line:
900 coord = line.split()[6]
901 if coord == "x":
902 xyz[0] = 1
903 elif coord == "y":
904 xyz[1] = 1
905 elif coord == "z":
906 xyz[2] = 1
907 keep = True
908 for n, c in enumerate(fix_cart):
909 if ind == c.index:
910 keep = False
911 break
912 if keep:
913 fix_cart.append(FixCartesian(ind, xyz))
914 else:
915 fix_cart[n].mask[xyz.index(1)] = 1
916 if len(fix) > 0:
917 fix_cart.append(FixAtoms(indices=fix))
919 return fix_cart
921 @cached_property
922 def initial_cell(self):
923 """Parse the initial cell from the aims.out file"""
924 line_start = self.reverse_search_for(["| Unit cell:"])
925 if line_start == LINE_NOT_FOUND:
926 return None
928 return [
929 [float(inp) for inp in line.split()[-3:]]
930 for line in self.lines[line_start + 1:line_start + 4]
931 ]
933 @cached_property
934 def initial_atoms(self):
935 """Create an atoms object for the initial geometry.in structure
936 from the aims.out file"""
937 line_start = self.reverse_search_for(["Atomic structure:"])
938 if line_start == LINE_NOT_FOUND:
939 raise AimsParseError(
940 "No information about the structure in the chunk.")
942 line_start += 2
944 cell = self.initial_cell
945 positions = np.zeros((self.n_atoms, 3))
946 symbols = [""] * self.n_atoms
947 for ll, line in enumerate(
948 self.lines[line_start:line_start + self.n_atoms]):
949 inp = line.split()
950 positions[ll, :] = [float(pos) for pos in inp[4:7]]
951 symbols[ll] = inp[3]
953 atoms = Atoms(symbols=symbols, positions=positions)
955 if cell:
956 atoms.set_cell(cell)
957 atoms.set_pbc([True, True, True])
958 atoms.set_constraint(self.constraints)
960 return atoms
962 @cached_property
963 def is_md(self):
964 """Determine if calculation is a molecular dynamics calculation"""
965 return LINE_NOT_FOUND != self.reverse_search_for(
966 ["Complete information for previous time-step:"]
967 )
969 @cached_property
970 def is_relaxation(self):
971 """Determine if the calculation is a geometry optimization or not"""
972 return LINE_NOT_FOUND != self.reverse_search_for(
973 ["Geometry relaxation:"])
975 @cached_property
976 def _k_points(self):
977 """Get the list of k-points used in the calculation"""
978 n_kpts = self.parse_scalar("n_kpts")
979 if n_kpts is None:
980 return {
981 "k_points": None,
982 "k_point_weights": None,
983 }
984 n_kpts = int(n_kpts)
986 line_start = self.reverse_search_for(["| K-points in task"])
987 line_end = self.reverse_search_for(["| k-point:"])
988 if (
989 (line_start == LINE_NOT_FOUND)
990 or (line_end == LINE_NOT_FOUND)
991 or (line_end - line_start != n_kpts)
992 ):
993 return {
994 "k_points": None,
995 "k_point_weights": None,
996 }
998 k_points = np.zeros((n_kpts, 3))
999 k_point_weights = np.zeros(n_kpts)
1000 for kk, line in enumerate(self.lines[line_start + 1:line_end + 1]):
1001 k_points[kk] = [float(inp) for inp in line.split()[4:7]]
1002 k_point_weights[kk] = float(line.split()[-1])
1004 return {
1005 "k_points": k_points,
1006 "k_point_weights": k_point_weights,
1007 }
1009 @cached_property
1010 def n_atoms(self):
1011 """The number of atoms for the material"""
1012 n_atoms = self.parse_scalar("n_atoms")
1013 if n_atoms is None:
1014 raise AimsParseError(
1015 "No information about the number of atoms in the header."
1016 )
1017 return int(n_atoms)
1019 @cached_property
1020 def n_bands(self):
1021 """The number of Kohn-Sham states for the chunk"""
1022 line_start = self.reverse_search_for(
1023 scalar_property_to_line_key["n_bands"])
1025 if line_start == LINE_NOT_FOUND:
1026 raise AimsParseError(
1027 "No information about the number of Kohn-Sham states "
1028 "in the header.")
1030 line = self.lines[line_start]
1031 if "| Number of Kohn-Sham states" in line:
1032 return int(line.split(":")[-1].strip().split()[0])
1034 return int(line.split()[-1].strip()[:-1])
1036 @cached_property
1037 def n_electrons(self):
1038 """The number of electrons for the chunk"""
1039 line_start = self.reverse_search_for(
1040 scalar_property_to_line_key["n_electrons"])
1042 if line_start == LINE_NOT_FOUND:
1043 raise AimsParseError(
1044 "No information about the number of electrons in the header."
1045 )
1047 line = self.lines[line_start]
1048 return int(float(line.split()[-2]))
1050 @cached_property
1051 def n_k_points(self):
1052 """The number of k_ppoints for the calculation"""
1053 n_kpts = self.parse_scalar("n_kpts")
1054 if n_kpts is None:
1055 return None
1057 return int(n_kpts)
1059 @cached_property
1060 def n_spins(self):
1061 """The number of spin channels for the chunk"""
1062 n_spins = self.parse_scalar("n_spins")
1063 if n_spins is None:
1064 raise AimsParseError(
1065 "No information about the number of spin "
1066 "channels in the header.")
1067 return int(n_spins)
1069 @cached_property
1070 def electronic_temperature(self):
1071 """The electronic temperature for the chunk"""
1072 line_start = self.reverse_search_for(
1073 scalar_property_to_line_key["electronic_temp"]
1074 )
1075 if line_start == LINE_NOT_FOUND:
1076 return 0.10
1078 line = self.lines[line_start]
1079 return float(line.split("=")[-1].strip().split()[0])
1081 @property
1082 def k_points(self):
1083 """All k-points listed in the calculation"""
1084 return self._k_points["k_points"]
1086 @property
1087 def k_point_weights(self):
1088 """The k-point weights for the calculation"""
1089 return self._k_points["k_point_weights"]
1091 @cached_property
1092 def header_summary(self):
1093 """Dictionary summarizing the information inside the header"""
1094 return {
1095 "initial_atoms": self.initial_atoms,
1096 "initial_cell": self.initial_cell,
1097 "constraints": self.constraints,
1098 "is_relaxation": self.is_relaxation,
1099 "is_md": self.is_md,
1100 "n_atoms": self.n_atoms,
1101 "n_bands": self.n_bands,
1102 "n_electrons": self.n_electrons,
1103 "n_spins": self.n_spins,
1104 "electronic_temperature": self.electronic_temperature,
1105 "n_k_points": self.n_k_points,
1106 "k_points": self.k_points,
1107 "k_point_weights": self.k_point_weights,
1108 }
1111class AimsOutCalcChunk(AimsOutChunk):
1112 """A part of the aims.out file correponding to a single structure"""
1114 def __init__(self, lines, header):
1115 """Constructor
1117 Parameters
1118 ----------
1119 lines: list of str
1120 The lines used for the structure
1121 header: dict
1122 A summary of the relevant information from the aims.out header
1123 """
1124 super().__init__(lines)
1125 self._header = header.header_summary
1127 @cached_property
1128 def _atoms(self):
1129 """Create an atoms object for the subsequent structures
1130 calculated in the aims.out file"""
1131 start_keys = [
1132 "Atomic structure (and velocities) as used in the preceding "
1133 "time step",
1134 "Updated atomic structure",
1135 "Atomic structure that was used in the preceding time step of "
1136 "the wrapper",
1137 ]
1138 line_start = self.reverse_search_for(start_keys)
1139 if line_start == LINE_NOT_FOUND:
1140 return self.initial_atoms
1142 line_start += 1
1144 line_end = self.reverse_search_for(
1145 [
1146 'Next atomic structure:',
1147 'Writing the current geometry to file "geometry.in.next_step"'
1148 ],
1149 line_start
1150 )
1151 if line_end == LINE_NOT_FOUND:
1152 line_end = len(self.lines)
1154 cell = []
1155 velocities = []
1156 atoms = Atoms()
1157 for line in self.lines[line_start:line_end]:
1158 if "lattice_vector " in line:
1159 cell.append([float(inp) for inp in line.split()[1:]])
1160 elif "atom " in line:
1161 line_split = line.split()
1162 atoms.append(Atom(line_split[4], tuple(
1163 float(inp) for inp in line_split[1:4])))
1164 elif "velocity " in line:
1165 velocities.append([float(inp) for inp in line.split()[1:]])
1167 assert len(atoms) == self.n_atoms
1168 assert (len(velocities) == self.n_atoms) or (len(velocities) == 0)
1169 if len(cell) == 3:
1170 atoms.set_cell(np.array(cell))
1171 atoms.set_pbc([True, True, True])
1172 elif len(cell) != 0:
1173 raise AimsParseError(
1174 "Parsed geometry has incorrect number of lattice vectors."
1175 )
1177 if len(velocities) > 0:
1178 atoms.set_velocities(np.array(velocities))
1179 atoms.set_constraint(self.constraints)
1181 return atoms
1183 @cached_property
1184 def forces(self):
1185 """Parse the forces from the aims.out file"""
1186 line_start = self.reverse_search_for(["Total atomic forces"])
1187 if line_start == LINE_NOT_FOUND:
1188 return None
1190 line_start += 1
1192 return np.array(
1193 [
1194 [float(inp) for inp in line.split()[-3:]]
1195 for line in self.lines[line_start:line_start + self.n_atoms]
1196 ]
1197 )
1199 @cached_property
1200 def stresses(self):
1201 """Parse the stresses from the aims.out file"""
1202 line_start = self.reverse_search_for(
1203 ["Per atom stress (eV) used for heat flux calculation"]
1204 )
1205 if line_start == LINE_NOT_FOUND:
1206 return None
1207 line_start += 3
1208 stresses = []
1209 for line in self.lines[line_start:line_start + self.n_atoms]:
1210 xx, yy, zz, xy, xz, yz = (float(d) for d in line.split()[2:8])
1211 stresses.append([xx, yy, zz, yz, xz, xy])
1213 return np.array(stresses)
1215 @cached_property
1216 def stress(self):
1217 """Parse the stress from the aims.out file"""
1218 from ase.stress import full_3x3_to_voigt_6_stress
1220 line_start = self.reverse_search_for(
1221 [
1222 "Analytical stress tensor - Symmetrized",
1223 "Numerical stress tensor",
1224 ]
1226 ) # Offest to relevant lines
1227 if line_start == LINE_NOT_FOUND:
1228 return None
1230 stress = [
1231 [float(inp) for inp in line.split()[2:5]]
1232 for line in self.lines[line_start + 5:line_start + 8]
1233 ]
1234 return full_3x3_to_voigt_6_stress(stress)
1236 @cached_property
1237 def is_metallic(self):
1238 """Checks the outputfile to see if the chunk corresponds
1239 to a metallic system"""
1240 line_start = self.reverse_search_for(
1241 ["material is metallic within the approximate finite "
1242 "broadening function (occupation_type)"])
1243 return line_start != LINE_NOT_FOUND
1245 @cached_property
1246 def total_energy(self):
1247 """Parse the energy from the aims.out file"""
1248 if np.all(self._atoms.pbc) and self.is_metallic:
1249 line_ind = self.reverse_search_for(["Total energy corrected"])
1250 else:
1251 line_ind = self.reverse_search_for(["Total energy uncorrected"])
1252 if line_ind == LINE_NOT_FOUND:
1253 raise AimsParseError("No energy is associated with the structure.")
1255 return float(self.lines[line_ind].split()[5])
1257 @cached_property
1258 def dipole(self):
1259 """Parse the electric dipole moment from the aims.out file."""
1260 line_start = self.reverse_search_for(["Total dipole moment [eAng]"])
1261 if line_start == LINE_NOT_FOUND:
1262 return None
1264 line = self.lines[line_start]
1265 return np.array([float(inp) for inp in line.split()[6:9]])
1267 @cached_property
1268 def dielectric_tensor(self):
1269 """Parse the dielectric tensor from the aims.out file"""
1270 line_start = self.reverse_search_for(
1271 ["DFPT for dielectric_constant:--->",
1272 "PARSE DFPT_dielectric_tensor"],
1273 )
1274 if line_start == LINE_NOT_FOUND:
1275 return None
1277 # we should find the tensor in the next three lines:
1278 lines = self.lines[line_start + 1:line_start + 4]
1280 # make ndarray and return
1281 return np.array([np.fromstring(line, sep=' ') for line in lines])
1283 @cached_property
1284 def polarization(self):
1285 """ Parse the polarization vector from the aims.out file"""
1286 line_start = self.reverse_search_for(["| Cartesian Polarization"])
1287 if line_start == LINE_NOT_FOUND:
1288 return None
1289 line = self.lines[line_start]
1290 return np.array([float(s) for s in line.split()[-3:]])
1292 @cached_property
1293 def _hirshfeld(self):
1294 """Parse the Hirshfled charges volumes, and dipole moments from the
1295 ouput"""
1296 line_start = self.reverse_search_for(
1297 ["Performing Hirshfeld analysis of fragment charges and moments."]
1298 )
1299 if line_start == LINE_NOT_FOUND:
1300 return {
1301 "charges": None,
1302 "volumes": None,
1303 "atomic_dipoles": None,
1304 "dipole": None,
1305 }
1307 line_inds = self.search_for_all("Hirshfeld charge", line_start, -1)
1308 hirshfeld_charges = np.array(
1309 [float(self.lines[ind].split(":")[1]) for ind in line_inds]
1310 )
1312 line_inds = self.search_for_all("Hirshfeld volume", line_start, -1)
1313 hirshfeld_volumes = np.array(
1314 [float(self.lines[ind].split(":")[1]) for ind in line_inds]
1315 )
1317 line_inds = self.search_for_all(
1318 "Hirshfeld dipole vector", line_start, -1)
1319 hirshfeld_atomic_dipoles = np.array(
1320 [
1321 [float(inp) for inp in self.lines[ind].split(":")[1].split()]
1322 for ind in line_inds
1323 ]
1324 )
1326 if not np.any(self._atoms.pbc):
1327 positions = self._atoms.get_positions()
1328 hirshfeld_dipole = np.sum(
1329 hirshfeld_charges.reshape((-1, 1)) * positions,
1330 axis=1,
1331 )
1332 else:
1333 hirshfeld_dipole = None
1334 return {
1335 "charges": hirshfeld_charges,
1336 "volumes": hirshfeld_volumes,
1337 "atomic_dipoles": hirshfeld_atomic_dipoles,
1338 "dipole": hirshfeld_dipole,
1339 }
1341 @cached_property
1342 def _eigenvalues(self):
1343 """Parse the eigenvalues and occupancies of the system. If eigenvalue
1344 for a particular k-point is not present in the output file
1345 then set it to np.nan
1346 """
1348 line_start = self.reverse_search_for(["Writing Kohn-Sham eigenvalues."])
1349 if line_start == LINE_NOT_FOUND:
1350 return {"eigenvalues": None, "occupancies": None}
1352 line_end_1 = self.reverse_search_for(
1353 ["Self-consistency cycle converged."], line_start
1354 )
1355 line_end_2 = self.reverse_search_for(
1356 [
1357 "What follows are estimated values for band gap, "
1358 "HOMO, LUMO, etc.",
1359 "Current spin moment of the entire structure :",
1360 "Highest occupied state (VBM)"
1361 ],
1362 line_start,
1363 )
1364 if line_end_1 == LINE_NOT_FOUND:
1365 line_end = line_end_2
1366 elif line_end_2 == LINE_NOT_FOUND:
1367 line_end = line_end_1
1368 else:
1369 line_end = min(line_end_1, line_end_2)
1371 n_kpts = self.n_k_points if np.all(self._atoms.pbc) else 1
1372 if n_kpts is None:
1373 return {"eigenvalues": None, "occupancies": None}
1375 eigenvalues = np.full((n_kpts, self.n_bands, self.n_spins), np.nan)
1376 occupancies = np.full((n_kpts, self.n_bands, self.n_spins), np.nan)
1378 occupation_block_start = self.search_for_all(
1379 "State Occupation Eigenvalue [Ha] Eigenvalue [eV]",
1380 line_start,
1381 line_end,
1382 )
1383 kpt_def = self.search_for_all("K-point: ", line_start, line_end)
1385 if len(kpt_def) > 0:
1386 kpt_inds = [int(self.lines[ll].split()[1]) - 1 for ll in kpt_def]
1387 elif (self.n_k_points is None) or (self.n_k_points == 1):
1388 kpt_inds = [0]
1389 else:
1390 raise ParseError("Cannot find k-point definitions")
1392 assert len(kpt_inds) == len(occupation_block_start)
1393 spins = [0] * len(occupation_block_start)
1395 if self.n_spins == 2:
1396 spin_def = self.search_for_all("Spin-", line_start, line_end)
1397 assert len(spin_def) == len(occupation_block_start)
1399 spins = [int("Spin-down eigenvalues:" in self.lines[ll])
1400 for ll in spin_def]
1402 for occ_start, kpt_ind, spin in zip(
1403 occupation_block_start, kpt_inds, spins):
1404 for ll, line in enumerate(
1405 self.lines[occ_start + 1:occ_start + self.n_bands + 1]
1406 ):
1407 if "***" in line:
1408 warn_msg = f"The {ll + 1}th eigenvalue for the "
1409 "{kpt_ind+1}th k-point and {spin}th channels could "
1410 "not be read (likely too large to be printed "
1411 "in the output file)"
1412 warnings.warn(warn_msg)
1413 continue
1414 split_line = line.split()
1415 eigenvalues[kpt_ind, ll, spin] = float(split_line[3])
1416 occupancies[kpt_ind, ll, spin] = float(split_line[1])
1417 return {"eigenvalues": eigenvalues, "occupancies": occupancies}
1419 @cached_property
1420 def atoms(self):
1421 """Convert AimsOutChunk to Atoms object and add all non-standard
1422outputs to atoms.info"""
1423 atoms = self._atoms
1425 atoms.calc = SinglePointDFTCalculator(
1426 atoms,
1427 energy=self.free_energy,
1428 free_energy=self.free_energy,
1429 forces=self.forces,
1430 stress=self.stress,
1431 stresses=self.stresses,
1432 magmom=self.magmom,
1433 dipole=self.dipole,
1434 dielectric_tensor=self.dielectric_tensor,
1435 polarization=self.polarization,
1436 )
1437 return atoms
1439 @property
1440 def results(self):
1441 """Convert an AimsOutChunk to a Results Dictionary"""
1442 results = {
1443 "energy": self.free_energy,
1444 "free_energy": self.free_energy,
1445 "total_energy": self.total_energy,
1446 "forces": self.forces,
1447 "stress": self.stress,
1448 "stresses": self.stresses,
1449 "magmom": self.magmom,
1450 "dipole": self.dipole,
1451 "fermi_energy": self.E_f,
1452 "n_iter": self.n_iter,
1453 "hirshfeld_charges": self.hirshfeld_charges,
1454 "hirshfeld_dipole": self.hirshfeld_dipole,
1455 "hirshfeld_volumes": self.hirshfeld_volumes,
1456 "hirshfeld_atomic_dipoles": self.hirshfeld_atomic_dipoles,
1457 "eigenvalues": self.eigenvalues,
1458 "occupancies": self.occupancies,
1459 "dielectric_tensor": self.dielectric_tensor,
1460 "polarization": self.polarization,
1461 }
1463 return {
1464 key: value for key,
1465 value in results.items() if value is not None}
1467 @property
1468 def initial_atoms(self):
1469 """The initial structure defined in the geoemtry.in file"""
1470 return self._header["initial_atoms"]
1472 @property
1473 def initial_cell(self):
1474 """The initial lattice vectors defined in the geoemtry.in file"""
1475 return self._header["initial_cell"]
1477 @property
1478 def constraints(self):
1479 """The relaxation constraints for the calculation"""
1480 return self._header["constraints"]
1482 @property
1483 def n_atoms(self):
1484 """The number of atoms for the material"""
1485 return self._header["n_atoms"]
1487 @property
1488 def n_bands(self):
1489 """The number of Kohn-Sham states for the chunk"""
1490 return self._header["n_bands"]
1492 @property
1493 def n_electrons(self):
1494 """The number of electrons for the chunk"""
1495 return self._header["n_electrons"]
1497 @property
1498 def n_spins(self):
1499 """The number of spin channels for the chunk"""
1500 return self._header["n_spins"]
1502 @property
1503 def electronic_temperature(self):
1504 """The electronic temperature for the chunk"""
1505 return self._header["electronic_temperature"]
1507 @property
1508 def n_k_points(self):
1509 """The number of electrons for the chunk"""
1510 return self._header["n_k_points"]
1512 @property
1513 def k_points(self):
1514 """The number of spin channels for the chunk"""
1515 return self._header["k_points"]
1517 @property
1518 def k_point_weights(self):
1519 """k_point_weights electronic temperature for the chunk"""
1520 return self._header["k_point_weights"]
1522 @cached_property
1523 def free_energy(self):
1524 """The free energy for the chunk"""
1525 return self.parse_scalar("free_energy")
1527 @cached_property
1528 def n_iter(self):
1529 """The number of SCF iterations needed to converge the SCF cycle for
1530the chunk"""
1531 return self.parse_scalar("number_of_iterations")
1533 @cached_property
1534 def magmom(self):
1535 """The magnetic moment for the chunk"""
1536 return self.parse_scalar("magnetic_moment")
1538 @cached_property
1539 def E_f(self):
1540 """The Fermi energy for the chunk"""
1541 return self.parse_scalar("fermi_energy")
1543 @cached_property
1544 def converged(self):
1545 """True if the chunk is a fully converged final structure"""
1546 return (len(self.lines) > 0) and ("Have a nice day." in self.lines[-5:])
1548 @property
1549 def hirshfeld_charges(self):
1550 """The Hirshfeld charges for the chunk"""
1551 return self._hirshfeld["charges"]
1553 @property
1554 def hirshfeld_atomic_dipoles(self):
1555 """The Hirshfeld atomic dipole moments for the chunk"""
1556 return self._hirshfeld["atomic_dipoles"]
1558 @property
1559 def hirshfeld_volumes(self):
1560 """The Hirshfeld volume for the chunk"""
1561 return self._hirshfeld["volumes"]
1563 @property
1564 def hirshfeld_dipole(self):
1565 """The Hirshfeld systematic dipole moment for the chunk"""
1566 if not np.any(self._atoms.pbc):
1567 return self._hirshfeld["dipole"]
1569 return None
1571 @property
1572 def eigenvalues(self):
1573 """All outputted eigenvalues for the system"""
1574 return self._eigenvalues["eigenvalues"]
1576 @property
1577 def occupancies(self):
1578 """All outputted occupancies for the system"""
1579 return self._eigenvalues["occupancies"]
1582def get_header_chunk(fd):
1583 """Returns the header information from the aims.out file"""
1584 header = []
1585 line = ""
1587 # Stop the header once the first SCF cycle begins
1588 while (
1589 "Convergence: q app. | density | eigen (eV) | Etot (eV)"
1590 not in line
1591 and "Convergence: q app. | density, spin | eigen (eV) |"
1592 not in line
1593 and "Begin self-consistency iteration #" not in line
1594 ):
1595 try:
1596 line = next(fd).strip() # Raises StopIteration on empty file
1597 except StopIteration:
1598 raise ParseError(
1599 "No SCF steps present, calculation failed at setup."
1600 )
1602 header.append(line)
1603 return AimsOutHeaderChunk(header)
1606def get_aims_out_chunks(fd, header_chunk):
1607 """Yield unprocessed chunks (header, lines) for each AimsOutChunk image."""
1608 try:
1609 line = next(fd).strip() # Raises StopIteration on empty file
1610 except StopIteration:
1611 return
1613 # If the calculation is relaxation the updated structural information
1614 # occurs before the re-initialization
1615 if header_chunk.is_relaxation:
1616 chunk_end_line = (
1617 "Geometry optimization: Attempting to predict improved coordinates."
1618 )
1619 else:
1620 chunk_end_line = "Begin self-consistency loop: Re-initialization"
1622 # If SCF is not converged then do not treat the next chunk_end_line as a
1623 # new chunk until after the SCF is re-initialized
1624 ignore_chunk_end_line = False
1625 while True:
1626 try:
1627 line = next(fd).strip() # Raises StopIteration on empty file
1628 except StopIteration:
1629 break
1631 lines = []
1632 while chunk_end_line not in line or ignore_chunk_end_line:
1633 lines.append(line)
1634 # If SCF cycle not converged or numerical stresses are requested,
1635 # don't end chunk on next Re-initialization
1636 patterns = [
1637 (
1638 "Self-consistency cycle not yet converged -"
1639 " restarting mixer to attempt better convergence."
1640 ),
1641 (
1642 "Components of the stress tensor (for mathematical "
1643 "background see comments in numerical_stress.f90)."
1644 ),
1645 "Calculation of numerical stress completed",
1646 ]
1647 if any(pattern in line for pattern in patterns):
1648 ignore_chunk_end_line = True
1649 elif "Begin self-consistency loop: Re-initialization" in line:
1650 ignore_chunk_end_line = False
1652 try:
1653 line = next(fd).strip()
1654 except StopIteration:
1655 break
1657 yield AimsOutCalcChunk(lines, header_chunk)
1660def check_convergence(chunks, non_convergence_ok=False):
1661 """Check if the aims output file is for a converged calculation
1663 Parameters
1664 ----------
1665 chunks: list of AimsOutChunks
1666 The list of chunks for the aims calculations
1667 non_convergence_ok: bool
1668 True if it is okay for the calculation to not be converged
1670 Returns
1671 -------
1672 bool
1673 True if the calculation is converged
1674 """
1675 if not non_convergence_ok and not chunks[-1].converged:
1676 raise ParseError("The calculation did not complete successfully")
1677 return True
1680@reader
1681def read_aims_output(fd, index=-1, non_convergence_ok=False):
1682 """Import FHI-aims output files with all data available, i.e.
1683 relaxations, MD information, force information etc etc etc."""
1684 header_chunk = get_header_chunk(fd)
1685 chunks = list(get_aims_out_chunks(fd, header_chunk))
1686 check_convergence(chunks, non_convergence_ok)
1688 # Relaxations have an additional footer chunk due to how it is split
1689 if header_chunk.is_relaxation:
1690 images = [chunk.atoms for chunk in chunks[:-1]]
1691 else:
1692 images = [chunk.atoms for chunk in chunks]
1693 return images[index]
1696@reader
1697def read_aims_results(fd, index=-1, non_convergence_ok=False):
1698 """Import FHI-aims output files and summarize all relevant information
1699 into a dictionary"""
1700 header_chunk = get_header_chunk(fd)
1701 chunks = list(get_aims_out_chunks(fd, header_chunk))
1702 check_convergence(chunks, non_convergence_ok)
1704 # Relaxations have an additional footer chunk due to how it is split
1705 if header_chunk.is_relaxation and (index == -1):
1706 return chunks[-2].results
1708 return chunks[index].results