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