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