Coverage for /builds/alexhroom/ase/ase/io/aims.py: 93.09%
811 statements
« prev ^ index » next coverage.py v7.5.3, created at 2024-08-05 14:37 +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 == "aims_command":
545 continue
546 elif key == "plus_u":
547 continue
548 elif key == "smearing":
549 name = parameters["smearing"][0].lower()
550 if name == "fermi-dirac":
551 name = "fermi"
552 width = parameters["smearing"][1]
553 if name == "methfessel-paxton":
554 order = parameters["smearing"][2]
555 order = " %d" % order
556 else:
557 order = ""
559 fd.write(
560 format_aims_control_parameter(
561 "occupation_type", (name, width, order), "%s %f%s"
562 )
563 )
564 elif key == "output":
565 for output_type in value:
566 fd.write(format_aims_control_parameter(key, output_type, "%s"))
567 elif key == "vdw_correction_hirshfeld" and value:
568 fd.write(format_aims_control_parameter(key, "", "%s"))
569 elif isinstance(value, bool):
570 fd.write(
571 format_aims_control_parameter(
572 key, str(value).lower(), ".%s."))
573 elif isinstance(value, (tuple, list)):
574 fd.write(
575 format_aims_control_parameter(
576 key, " ".join([str(x) for x in value]), "%s"
577 )
578 )
579 elif isinstance(value, str):
580 fd.write(format_aims_control_parameter(key, value, "%s"))
581 else:
582 fd.write(format_aims_control_parameter(key, value, "%r"))
584 if cubes:
585 cubes.write(fd)
587 fd.write(lim + "\n\n")
589 # Get the species directory
590 species_dir = get_species_directory
591 # dicts are ordered as of python 3.7
592 species_array = np.array(list(dict.fromkeys(atoms.symbols)))
593 # Grab the tier specification from the parameters. THis may either
594 # be None, meaning the default should be used for all species, or a
595 # list of integers/None values giving a specific basis set size
596 # for each species in the calculation.
597 tier = parameters.pop("tier", None)
598 tier_array = np.full(len(species_array), tier)
599 # Path to species files for FHI-aims. In this files are specifications
600 # for the basis set sizes depending on which basis set tier is used.
601 species_dir = get_species_directory(parameters.get("species_dir"))
602 # Parse the species files for each species present in the calculation
603 # according to the tier of each species.
604 species_basis_dict = parse_species_path(
605 species_array=species_array, tier_array=tier_array,
606 species_dir=species_dir)
607 # Write the basis functions to be included for each species in the
608 # calculation into the control.in file (fd).
609 write_species(fd, species_basis_dict, parameters)
612def get_species_directory(species_dir=None):
613 """Get the directory where the basis set information is stored
615 If the requested directory does not exist then raise an Error
617 Parameters
618 ----------
619 species_dir: str
620 Requested directory to find the basis set info from. E.g.
621 `~/aims2022/FHIaims/species_defaults/defaults_2020/light`.
623 Returns
624 -------
625 Path
626 The Path to the requested or default species directory.
628 Raises
629 ------
630 RuntimeError
631 If both the requested directory and the default one is not defined
632 or does not exit.
633 """
634 if species_dir is None:
635 species_dir = os.environ.get("AIMS_SPECIES_DIR")
637 if species_dir is None:
638 raise RuntimeError(
639 "Missing species directory! Use species_dir "
640 + "parameter or set $AIMS_SPECIES_DIR environment variable."
641 )
643 species_path = Path(species_dir)
644 if not species_path.exists():
645 raise RuntimeError(
646 f"The requested species_dir {species_dir} does not exist")
648 return species_path
651def write_species(control_file_descriptor, species_basis_dict, parameters):
652 """Write species for the calculation depending on basis set size.
654 The calculation should include certain basis set size function depending
655 on the numerical settings (light, tight, really tight) and the basis set
656 size (minimal, tier1, tier2, tier3, tier4). If the basis set size is not
657 given then a 'standard' basis set size is used for each numerical setting.
658 The species files are defined according to these standard basis set sizes
659 for the numerical settings in the FHI-aims repository.
661 Note, for FHI-aims in ASE, we don't explicitly give the numerical setting.
662 Instead we include the numerical setting in the species path: e.g.
663 `~/aims2022/FHIaims/species_defaults/defaults_2020/light` this path has
664 `light`, the numerical setting, as the last folder in the path.
666 Example - a basis function might be commented in the standard basis set size
667 such as "# hydro 4 f 7.4" and this basis function should be
668 uncommented for another basis set size such as tier4.
670 Args:
671 control_file_descriptor: File descriptor for the control.in file into
672 which we need to write relevant basis functions to be included for
673 the calculation.
674 species_basis_dict: Dictionary where keys as the species symbols and
675 each value is a single string containing all the basis functions
676 to be included in the caclculation.
677 parameters: Calculation parameters as a dict.
678 """
679 # Now for every species (key) in the species_basis_dict, save the
680 # relevant basis functions (values) from the species_basis_dict, by
681 # writing to the file handle (species_file_descriptor) given to this
682 # function.
683 for species_symbol, basis_set_text in species_basis_dict.items():
684 control_file_descriptor.write(basis_set_text)
685 if parameters.get("plus_u") is not None:
686 if species_symbol in parameters.plus_u:
687 control_file_descriptor.write(
688 f"plus_u {parameters.plus_u[species_symbol]} \n")
691def parse_species_path(species_array, tier_array, species_dir):
692 """Parse the species files for each species according to the tier given.
694 Args:
695 species_array: An array of species/element symbols present in the unit
696 cell (e.g. ['C', 'H'].)
697 tier_array: An array of None/integer values which define which basis
698 set size to use for each species/element in the calcualtion.
699 species_dir: Directory containing FHI-aims species files.
701 Returns:
702 Dictionary containing species as keys and the basis set specification
703 for each species as text as the value for the key.
704 """
705 if len(species_array) != len(tier_array):
706 raise ValueError(
707 f"The species array length: {len(species_array)}, "
708 f"is not the same as the tier_array length: {len(tier_array)}")
710 species_basis_dict = {}
712 for symbol, tier in zip(species_array, tier_array):
713 path = species_dir / f"{atomic_numbers[symbol]:02}_{symbol}_default"
714 # Open the species file:
715 with open(path, encoding="utf8") as species_file_handle:
716 # Read the species file into a string.
717 species_file_str = species_file_handle.read()
718 species_basis_dict[symbol] = manipulate_tiers(
719 species_file_str, tier)
720 return species_basis_dict
723def manipulate_tiers(species_string: str, tier: Union[None, int] = 1):
724 """Adds basis set functions based on the tier value.
726 This function takes in the species file as a string, it then searches
727 for relevant basis functions based on the tier value to include in a new
728 string that is returned.
730 Args:
731 species_string: species file (default) for a given numerical setting
732 (light, tight, really tight) given as a string.
733 tier: The basis set size. This will dictate which basis set functions
734 are included in the returned string.
736 Returns:
737 Basis set functions defined by the tier as a string.
738 """
739 if tier is None: # Then we use the default species file untouched.
740 return species_string
741 tier_pattern = r"(# \".* tier\" .*|# +Further.*)"
742 top, *tiers = re.split(tier_pattern, species_string)
743 tier_comments = tiers[::2]
744 tier_basis = tiers[1::2]
745 assert len(
746 tier_comments) == len(tier_basis), "Something wrong with splitting"
747 n_tiers = len(tier_comments)
748 assert tier <= n_tiers, f"Only {n_tiers} tiers available, you choose {tier}"
749 string_new = top
750 for i, (c, b) in enumerate(zip(tier_comments, tier_basis)):
751 b = re.sub(r"\n( *for_aux| *hydro| *ionic| *confined)", r"\n#\g<1>", b)
752 if i < tier:
753 b = re.sub(
754 r"\n#( *for_aux| *hydro| *ionic| *confined)", r"\n\g<1>", b)
755 string_new += c + b
756 return string_new
759# Read aims.out files
760scalar_property_to_line_key = {
761 "free_energy": ["| Electronic free energy"],
762 "number_of_iterations": ["| Number of self-consistency cycles"],
763 "magnetic_moment": ["N_up - N_down"],
764 "n_atoms": ["| Number of atoms"],
765 "n_bands": [
766 "Number of Kohn-Sham states",
767 "Reducing total number of Kohn-Sham states",
768 "Reducing total number of Kohn-Sham states",
769 ],
770 "n_electrons": ["The structure contains"],
771 "n_kpts": ["| Number of k-points"],
772 "n_spins": ["| Number of spin channels"],
773 "electronic_temp": ["Occupation type:"],
774 "fermi_energy": ["| Chemical potential (Fermi level)"],
775}
778class AimsOutChunk:
779 """Base class for AimsOutChunks"""
781 def __init__(self, lines):
782 """Constructor
784 Parameters
785 ----------
786 lines: list of str
787 The set of lines from the output file the encompasses either
788 a single structure within a trajectory or
789 general information about the calculation (header)
790 """
791 self.lines = lines
793 def reverse_search_for(self, keys, line_start=0):
794 """Find the last time one of the keys appears in self.lines
796 Parameters
797 ----------
798 keys: list of str
799 The key strings to search for in self.lines
800 line_start: int
801 The lowest index to search for in self.lines
803 Returns
804 -------
805 int
806 The last time one of the keys appears in self.lines
807 """
808 for ll, line in enumerate(self.lines[line_start:][::-1]):
809 if any(key in line for key in keys):
810 return len(self.lines) - ll - 1
812 return LINE_NOT_FOUND
814 def search_for_all(self, key, line_start=0, line_end=-1):
815 """Find the all times the key appears in self.lines
817 Parameters
818 ----------
819 key: str
820 The key string to search for in self.lines
821 line_start: int
822 The first line to start the search from
823 line_end: int
824 The last line to end the search at
826 Returns
827 -------
828 list of ints
829 All times the key appears in the lines
830 """
831 line_index = []
832 for ll, line in enumerate(self.lines[line_start:line_end]):
833 if key in line:
834 line_index.append(ll + line_start)
835 return line_index
837 def parse_scalar(self, property):
838 """Parse a scalar property from the chunk
840 Parameters
841 ----------
842 property: str
843 The property key to parse
845 Returns
846 -------
847 float
848 The scalar value of the property
849 """
850 line_start = self.reverse_search_for(
851 scalar_property_to_line_key[property])
853 if line_start == LINE_NOT_FOUND:
854 return None
856 line = self.lines[line_start]
857 return float(line.split(":")[-1].strip().split()[0])
860class AimsOutHeaderChunk(AimsOutChunk):
861 """The header of the aims.out file containint general information"""
863 def __init__(self, lines):
864 """Constructor
866 Parameters
867 ----------
868 lines: list of str
869 The lines inside the aims.out header
870 """
871 super().__init__(lines)
872 self._k_points = None
873 self._k_point_weights = None
875 @lazyproperty
876 def constraints(self):
877 """Parse the constraints from the aims.out file
879 Constraints for the lattice vectors are not supported.
880 """
882 line_inds = self.search_for_all("Found relaxation constraint for atom")
883 if len(line_inds) == 0:
884 return []
886 fix = []
887 fix_cart = []
888 for ll in line_inds:
889 line = self.lines[ll]
890 xyz = [0, 0, 0]
891 ind = int(line.split()[5][:-1]) - 1
892 if "All coordinates fixed" in line:
893 if ind not in fix:
894 fix.append(ind)
895 if "coordinate fixed" in line:
896 coord = line.split()[6]
897 if coord == "x":
898 xyz[0] = 1
899 elif coord == "y":
900 xyz[1] = 1
901 elif coord == "z":
902 xyz[2] = 1
903 keep = True
904 for n, c in enumerate(fix_cart):
905 if ind == c.index:
906 keep = False
907 break
908 if keep:
909 fix_cart.append(FixCartesian(ind, xyz))
910 else:
911 fix_cart[n].mask[xyz.index(1)] = 1
912 if len(fix) > 0:
913 fix_cart.append(FixAtoms(indices=fix))
915 return fix_cart
917 @lazyproperty
918 def initial_cell(self):
919 """Parse the initial cell from the aims.out file"""
920 line_start = self.reverse_search_for(["| Unit cell:"])
921 if line_start == LINE_NOT_FOUND:
922 return None
924 return [
925 [float(inp) for inp in line.split()[-3:]]
926 for line in self.lines[line_start + 1:line_start + 4]
927 ]
929 @lazyproperty
930 def initial_atoms(self):
931 """Create an atoms object for the initial geometry.in structure
932 from the aims.out file"""
933 line_start = self.reverse_search_for(["Atomic structure:"])
934 if line_start == LINE_NOT_FOUND:
935 raise AimsParseError(
936 "No information about the structure in the chunk.")
938 line_start += 2
940 cell = self.initial_cell
941 positions = np.zeros((self.n_atoms, 3))
942 symbols = [""] * self.n_atoms
943 for ll, line in enumerate(
944 self.lines[line_start:line_start + self.n_atoms]):
945 inp = line.split()
946 positions[ll, :] = [float(pos) for pos in inp[4:7]]
947 symbols[ll] = inp[3]
949 atoms = Atoms(symbols=symbols, positions=positions)
951 if cell:
952 atoms.set_cell(cell)
953 atoms.set_pbc([True, True, True])
954 atoms.set_constraint(self.constraints)
956 return atoms
958 @lazyproperty
959 def is_md(self):
960 """Determine if calculation is a molecular dynamics calculation"""
961 return LINE_NOT_FOUND != self.reverse_search_for(
962 ["Complete information for previous time-step:"]
963 )
965 @lazyproperty
966 def is_relaxation(self):
967 """Determine if the calculation is a geometry optimization or not"""
968 return LINE_NOT_FOUND != self.reverse_search_for(
969 ["Geometry relaxation:"])
971 @lazymethod
972 def _parse_k_points(self):
973 """Get the list of k-points used in the calculation"""
974 n_kpts = self.parse_scalar("n_kpts")
975 if n_kpts is None:
976 return {
977 "k_points": None,
978 "k_point_weights": None,
979 }
980 n_kpts = int(n_kpts)
982 line_start = self.reverse_search_for(["| K-points in task"])
983 line_end = self.reverse_search_for(["| k-point:"])
984 if (
985 (line_start == LINE_NOT_FOUND)
986 or (line_end == LINE_NOT_FOUND)
987 or (line_end - line_start != n_kpts)
988 ):
989 return {
990 "k_points": None,
991 "k_point_weights": None,
992 }
994 k_points = np.zeros((n_kpts, 3))
995 k_point_weights = np.zeros(n_kpts)
996 for kk, line in enumerate(self.lines[line_start + 1:line_end + 1]):
997 k_points[kk] = [float(inp) for inp in line.split()[4:7]]
998 k_point_weights[kk] = float(line.split()[-1])
1000 return {
1001 "k_points": k_points,
1002 "k_point_weights": k_point_weights,
1003 }
1005 @lazyproperty
1006 def n_atoms(self):
1007 """The number of atoms for the material"""
1008 n_atoms = self.parse_scalar("n_atoms")
1009 if n_atoms is None:
1010 raise AimsParseError(
1011 "No information about the number of atoms in the header."
1012 )
1013 return int(n_atoms)
1015 @lazyproperty
1016 def n_bands(self):
1017 """The number of Kohn-Sham states for the chunk"""
1018 line_start = self.reverse_search_for(
1019 scalar_property_to_line_key["n_bands"])
1021 if line_start == LINE_NOT_FOUND:
1022 raise AimsParseError(
1023 "No information about the number of Kohn-Sham states "
1024 "in the header.")
1026 line = self.lines[line_start]
1027 if "| Number of Kohn-Sham states" in line:
1028 return int(line.split(":")[-1].strip().split()[0])
1030 return int(line.split()[-1].strip()[:-1])
1032 @lazyproperty
1033 def n_electrons(self):
1034 """The number of electrons for the chunk"""
1035 line_start = self.reverse_search_for(
1036 scalar_property_to_line_key["n_electrons"])
1038 if line_start == LINE_NOT_FOUND:
1039 raise AimsParseError(
1040 "No information about the number of electrons in the header."
1041 )
1043 line = self.lines[line_start]
1044 return int(float(line.split()[-2]))
1046 @lazyproperty
1047 def n_k_points(self):
1048 """The number of k_ppoints for the calculation"""
1049 n_kpts = self.parse_scalar("n_kpts")
1050 if n_kpts is None:
1051 return None
1053 return int(n_kpts)
1055 @lazyproperty
1056 def n_spins(self):
1057 """The number of spin channels for the chunk"""
1058 n_spins = self.parse_scalar("n_spins")
1059 if n_spins is None:
1060 raise AimsParseError(
1061 "No information about the number of spin "
1062 "channels in the header.")
1063 return int(n_spins)
1065 @lazyproperty
1066 def electronic_temperature(self):
1067 """The electronic temperature for the chunk"""
1068 line_start = self.reverse_search_for(
1069 scalar_property_to_line_key["electronic_temp"]
1070 )
1071 if line_start == LINE_NOT_FOUND:
1072 return 0.10
1074 line = self.lines[line_start]
1075 return float(line.split("=")[-1].strip().split()[0])
1077 @lazyproperty
1078 def k_points(self):
1079 """All k-points listed in the calculation"""
1080 return self._parse_k_points()["k_points"]
1082 @lazyproperty
1083 def k_point_weights(self):
1084 """The k-point weights for the calculation"""
1085 return self._parse_k_points()["k_point_weights"]
1087 @lazyproperty
1088 def header_summary(self):
1089 """Dictionary summarizing the information inside the header"""
1090 return {
1091 "initial_atoms": self.initial_atoms,
1092 "initial_cell": self.initial_cell,
1093 "constraints": self.constraints,
1094 "is_relaxation": self.is_relaxation,
1095 "is_md": self.is_md,
1096 "n_atoms": self.n_atoms,
1097 "n_bands": self.n_bands,
1098 "n_electrons": self.n_electrons,
1099 "n_spins": self.n_spins,
1100 "electronic_temperature": self.electronic_temperature,
1101 "n_k_points": self.n_k_points,
1102 "k_points": self.k_points,
1103 "k_point_weights": self.k_point_weights,
1104 }
1107class AimsOutCalcChunk(AimsOutChunk):
1108 """A part of the aims.out file correponding to a single structure"""
1110 def __init__(self, lines, header):
1111 """Constructor
1113 Parameters
1114 ----------
1115 lines: list of str
1116 The lines used for the structure
1117 header: dict
1118 A summary of the relevant information from the aims.out header
1119 """
1120 super().__init__(lines)
1121 self._header = header.header_summary
1123 @lazymethod
1124 def _parse_atoms(self):
1125 """Create an atoms object for the subsequent structures
1126 calculated in the aims.out file"""
1127 start_keys = [
1128 "Atomic structure (and velocities) as used in the preceding "
1129 "time step",
1130 "Updated atomic structure",
1131 "Atomic structure that was used in the preceding time step of "
1132 "the wrapper",
1133 ]
1134 line_start = self.reverse_search_for(start_keys)
1135 if line_start == LINE_NOT_FOUND:
1136 return self.initial_atoms
1138 line_start += 1
1140 line_end = self.reverse_search_for(
1141 [
1142 'Next atomic structure:',
1143 'Writing the current geometry to file "geometry.in.next_step"'
1144 ],
1145 line_start
1146 )
1147 if line_end == LINE_NOT_FOUND:
1148 line_end = len(self.lines)
1150 cell = []
1151 velocities = []
1152 atoms = Atoms()
1153 for line in self.lines[line_start:line_end]:
1154 if "lattice_vector " in line:
1155 cell.append([float(inp) for inp in line.split()[1:]])
1156 elif "atom " in line:
1157 line_split = line.split()
1158 atoms.append(Atom(line_split[4], tuple(
1159 float(inp) for inp in line_split[1:4])))
1160 elif "velocity " in line:
1161 velocities.append([float(inp) for inp in line.split()[1:]])
1163 assert len(atoms) == self.n_atoms
1164 assert (len(velocities) == self.n_atoms) or (len(velocities) == 0)
1165 if len(cell) == 3:
1166 atoms.set_cell(np.array(cell))
1167 atoms.set_pbc([True, True, True])
1168 elif len(cell) != 0:
1169 raise AimsParseError(
1170 "Parsed geometry has incorrect number of lattice vectors."
1171 )
1173 if len(velocities) > 0:
1174 atoms.set_velocities(np.array(velocities))
1175 atoms.set_constraint(self.constraints)
1177 return atoms
1179 @lazyproperty
1180 def forces(self):
1181 """Parse the forces from the aims.out file"""
1182 line_start = self.reverse_search_for(["Total atomic forces"])
1183 if line_start == LINE_NOT_FOUND:
1184 return None
1186 line_start += 1
1188 return np.array(
1189 [
1190 [float(inp) for inp in line.split()[-3:]]
1191 for line in self.lines[line_start:line_start + self.n_atoms]
1192 ]
1193 )
1195 @lazyproperty
1196 def stresses(self):
1197 """Parse the stresses from the aims.out file"""
1198 line_start = self.reverse_search_for(
1199 ["Per atom stress (eV) used for heat flux calculation"]
1200 )
1201 if line_start == LINE_NOT_FOUND:
1202 return None
1203 line_start += 3
1204 stresses = []
1205 for line in self.lines[line_start:line_start + self.n_atoms]:
1206 xx, yy, zz, xy, xz, yz = (float(d) for d in line.split()[2:8])
1207 stresses.append([xx, yy, zz, yz, xz, xy])
1209 return np.array(stresses)
1211 @lazyproperty
1212 def stress(self):
1213 """Parse the stress from the aims.out file"""
1214 from ase.stress import full_3x3_to_voigt_6_stress
1216 line_start = self.reverse_search_for(
1217 [
1218 "Analytical stress tensor - Symmetrized",
1219 "Numerical stress tensor",
1220 ]
1222 ) # Offest to relevant lines
1223 if line_start == LINE_NOT_FOUND:
1224 return None
1226 stress = [
1227 [float(inp) for inp in line.split()[2:5]]
1228 for line in self.lines[line_start + 5:line_start + 8]
1229 ]
1230 return full_3x3_to_voigt_6_stress(stress)
1232 @lazyproperty
1233 def is_metallic(self):
1234 """Checks the outputfile to see if the chunk corresponds
1235 to a metallic system"""
1236 line_start = self.reverse_search_for(
1237 ["material is metallic within the approximate finite "
1238 "broadening function (occupation_type)"])
1239 return line_start != LINE_NOT_FOUND
1241 @lazyproperty
1242 def total_energy(self):
1243 """Parse the energy from the aims.out file"""
1244 atoms = self._parse_atoms()
1246 if np.all(atoms.pbc) and self.is_metallic:
1247 line_ind = self.reverse_search_for(["Total energy corrected"])
1248 else:
1249 line_ind = self.reverse_search_for(["Total energy uncorrected"])
1250 if line_ind == LINE_NOT_FOUND:
1251 raise AimsParseError("No energy is associated with the structure.")
1253 return float(self.lines[line_ind].split()[5])
1255 @lazyproperty
1256 def dipole(self):
1257 """Parse the electric dipole moment from the aims.out file."""
1258 line_start = self.reverse_search_for(["Total dipole moment [eAng]"])
1259 if line_start == LINE_NOT_FOUND:
1260 return None
1262 line = self.lines[line_start]
1263 return np.array([float(inp) for inp in line.split()[6:9]])
1265 @lazyproperty
1266 def dielectric_tensor(self):
1267 """Parse the dielectric tensor from the aims.out file"""
1268 line_start = self.reverse_search_for(["PARSE DFPT_dielectric_tensor"])
1269 if line_start == LINE_NOT_FOUND:
1270 return None
1272 # we should find the tensor in the next three lines:
1273 lines = self.lines[line_start + 1:line_start + 4]
1275 # make ndarray and return
1276 return np.array([np.fromstring(line, sep=' ') for line in lines])
1278 @lazyproperty
1279 def polarization(self):
1280 """ Parse the polarization vector from the aims.out file"""
1281 line_start = self.reverse_search_for(["| Cartesian Polarization"])
1282 if line_start == LINE_NOT_FOUND:
1283 return None
1284 line = self.lines[line_start]
1285 return np.array([float(s) for s in line.split()[-3:]])
1287 @lazymethod
1288 def _parse_hirshfeld(self):
1289 """Parse the Hirshfled charges volumes, and dipole moments from the
1290 ouput"""
1291 atoms = self._parse_atoms()
1293 line_start = self.reverse_search_for(
1294 ["Performing Hirshfeld analysis of fragment charges and moments."]
1295 )
1296 if line_start == LINE_NOT_FOUND:
1297 return {
1298 "charges": None,
1299 "volumes": None,
1300 "atomic_dipoles": None,
1301 "dipole": None,
1302 }
1304 line_inds = self.search_for_all("Hirshfeld charge", line_start, -1)
1305 hirshfeld_charges = np.array(
1306 [float(self.lines[ind].split(":")[1]) for ind in line_inds]
1307 )
1309 line_inds = self.search_for_all("Hirshfeld volume", line_start, -1)
1310 hirshfeld_volumes = np.array(
1311 [float(self.lines[ind].split(":")[1]) for ind in line_inds]
1312 )
1314 line_inds = self.search_for_all(
1315 "Hirshfeld dipole vector", line_start, -1)
1316 hirshfeld_atomic_dipoles = np.array(
1317 [
1318 [float(inp) for inp in self.lines[ind].split(":")[1].split()]
1319 for ind in line_inds
1320 ]
1321 )
1323 if not np.any(atoms.pbc):
1324 hirshfeld_dipole = np.sum(
1325 hirshfeld_charges.reshape((-1, 1)) * atoms.get_positions(),
1326 axis=1,
1327 )
1328 else:
1329 hirshfeld_dipole = None
1330 return {
1331 "charges": hirshfeld_charges,
1332 "volumes": hirshfeld_volumes,
1333 "atomic_dipoles": hirshfeld_atomic_dipoles,
1334 "dipole": hirshfeld_dipole,
1335 }
1337 @lazymethod
1338 def _parse_eigenvalues(self):
1339 """Parse the eigenvalues and occupancies of the system. If eigenvalue
1340 for a particular k-point is not present in the output file
1341 then set it to np.nan
1342 """
1344 atoms = self._parse_atoms()
1346 line_start = self.reverse_search_for(["Writing Kohn-Sham eigenvalues."])
1347 if line_start == LINE_NOT_FOUND:
1348 return {"eigenvalues": None, "occupancies": None}
1350 line_end_1 = self.reverse_search_for(
1351 ["Self-consistency cycle converged."], line_start
1352 )
1353 line_end_2 = self.reverse_search_for(
1354 [
1355 "What follows are estimated values for band gap, "
1356 "HOMO, LUMO, etc.",
1357 "Current spin moment of the entire structure :",
1358 "Highest occupied state (VBM)"
1359 ],
1360 line_start,
1361 )
1362 if line_end_1 == LINE_NOT_FOUND:
1363 line_end = line_end_2
1364 elif line_end_2 == LINE_NOT_FOUND:
1365 line_end = line_end_1
1366 else:
1367 line_end = min(line_end_1, line_end_2)
1369 n_kpts = self.n_k_points if np.all(atoms.pbc) else 1
1370 if n_kpts is None:
1371 return {"eigenvalues": None, "occupancies": None}
1373 eigenvalues = np.full((n_kpts, self.n_bands, self.n_spins), np.nan)
1374 occupancies = np.full((n_kpts, self.n_bands, self.n_spins), np.nan)
1376 occupation_block_start = self.search_for_all(
1377 "State Occupation Eigenvalue [Ha] Eigenvalue [eV]",
1378 line_start,
1379 line_end,
1380 )
1381 kpt_def = self.search_for_all("K-point: ", line_start, line_end)
1383 if len(kpt_def) > 0:
1384 kpt_inds = [int(self.lines[ll].split()[1]) - 1 for ll in kpt_def]
1385 elif (self.n_k_points is None) or (self.n_k_points == 1):
1386 kpt_inds = [0]
1387 else:
1388 raise ParseError("Cannot find k-point definitions")
1390 assert len(kpt_inds) == len(occupation_block_start)
1391 spins = [0] * len(occupation_block_start)
1393 if self.n_spins == 2:
1394 spin_def = self.search_for_all("Spin-", line_start, line_end)
1395 assert len(spin_def) == len(occupation_block_start)
1397 spins = [int("Spin-down eigenvalues:" in self.lines[ll])
1398 for ll in spin_def]
1400 for occ_start, kpt_ind, spin in zip(
1401 occupation_block_start, kpt_inds, spins):
1402 for ll, line in enumerate(
1403 self.lines[occ_start + 1:occ_start + self.n_bands + 1]
1404 ):
1405 if "***" in line:
1406 warn_msg = f"The {ll+1}th eigenvalue for the "
1407 "{kpt_ind+1}th k-point and {spin}th channels could "
1408 "not be read (likely too large to be printed "
1409 "in the output file)"
1410 warnings.warn(warn_msg)
1411 continue
1412 split_line = line.split()
1413 eigenvalues[kpt_ind, ll, spin] = float(split_line[3])
1414 occupancies[kpt_ind, ll, spin] = float(split_line[1])
1415 return {"eigenvalues": eigenvalues, "occupancies": occupancies}
1417 @lazyproperty
1418 def atoms(self):
1419 """Convert AimsOutChunk to Atoms object and add all non-standard
1420outputs to atoms.info"""
1421 atoms = self._parse_atoms()
1423 atoms.calc = SinglePointDFTCalculator(
1424 atoms,
1425 energy=self.free_energy,
1426 free_energy=self.free_energy,
1427 forces=self.forces,
1428 stress=self.stress,
1429 stresses=self.stresses,
1430 magmom=self.magmom,
1431 dipole=self.dipole,
1432 dielectric_tensor=self.dielectric_tensor,
1433 polarization=self.polarization,
1434 )
1435 return atoms
1437 @property
1438 def results(self):
1439 """Convert an AimsOutChunk to a Results Dictionary"""
1440 results = {
1441 "energy": self.free_energy,
1442 "free_energy": self.free_energy,
1443 "total_energy": self.total_energy,
1444 "forces": self.forces,
1445 "stress": self.stress,
1446 "stresses": self.stresses,
1447 "magmom": self.magmom,
1448 "dipole": self.dipole,
1449 "fermi_energy": self.E_f,
1450 "n_iter": self.n_iter,
1451 "hirshfeld_charges": self.hirshfeld_charges,
1452 "hirshfeld_dipole": self.hirshfeld_dipole,
1453 "hirshfeld_volumes": self.hirshfeld_volumes,
1454 "hirshfeld_atomic_dipoles": self.hirshfeld_atomic_dipoles,
1455 "eigenvalues": self.eigenvalues,
1456 "occupancies": self.occupancies,
1457 "dielectric_tensor": self.dielectric_tensor,
1458 "polarization": self.polarization,
1459 }
1461 return {
1462 key: value for key,
1463 value in results.items() if value is not None}
1465 # Properties from the aims.out header
1466 @lazyproperty
1467 def initial_atoms(self):
1468 """The initial structure defined in the geoemtry.in file"""
1469 return self._header["initial_atoms"]
1471 @lazyproperty
1472 def initial_cell(self):
1473 """The initial lattice vectors defined in the geoemtry.in file"""
1474 return self._header["initial_cell"]
1476 @lazyproperty
1477 def constraints(self):
1478 """The relaxation constraints for the calculation"""
1479 return self._header["constraints"]
1481 @lazyproperty
1482 def n_atoms(self):
1483 """The number of atoms for the material"""
1484 return self._header["n_atoms"]
1486 @lazyproperty
1487 def n_bands(self):
1488 """The number of Kohn-Sham states for the chunk"""
1489 return self._header["n_bands"]
1491 @lazyproperty
1492 def n_electrons(self):
1493 """The number of electrons for the chunk"""
1494 return self._header["n_electrons"]
1496 @lazyproperty
1497 def n_spins(self):
1498 """The number of spin channels for the chunk"""
1499 return self._header["n_spins"]
1501 @lazyproperty
1502 def electronic_temperature(self):
1503 """The electronic temperature for the chunk"""
1504 return self._header["electronic_temperature"]
1506 @lazyproperty
1507 def n_k_points(self):
1508 """The number of electrons for the chunk"""
1509 return self._header["n_k_points"]
1511 @lazyproperty
1512 def k_points(self):
1513 """The number of spin channels for the chunk"""
1514 return self._header["k_points"]
1516 @lazyproperty
1517 def k_point_weights(self):
1518 """k_point_weights electronic temperature for the chunk"""
1519 return self._header["k_point_weights"]
1521 @lazyproperty
1522 def free_energy(self):
1523 """The free energy for the chunk"""
1524 return self.parse_scalar("free_energy")
1526 @lazyproperty
1527 def n_iter(self):
1528 """The number of SCF iterations needed to converge the SCF cycle for
1529the chunk"""
1530 return self.parse_scalar("number_of_iterations")
1532 @lazyproperty
1533 def magmom(self):
1534 """The magnetic moment for the chunk"""
1535 return self.parse_scalar("magnetic_moment")
1537 @lazyproperty
1538 def E_f(self):
1539 """The Fermi energy for the chunk"""
1540 return self.parse_scalar("fermi_energy")
1542 @lazyproperty
1543 def converged(self):
1544 """True if the chunk is a fully converged final structure"""
1545 return (len(self.lines) > 0) and ("Have a nice day." in self.lines[-5:])
1547 @lazyproperty
1548 def hirshfeld_charges(self):
1549 """The Hirshfeld charges for the chunk"""
1550 return self._parse_hirshfeld()["charges"]
1552 @lazyproperty
1553 def hirshfeld_atomic_dipoles(self):
1554 """The Hirshfeld atomic dipole moments for the chunk"""
1555 return self._parse_hirshfeld()["atomic_dipoles"]
1557 @lazyproperty
1558 def hirshfeld_volumes(self):
1559 """The Hirshfeld volume for the chunk"""
1560 return self._parse_hirshfeld()["volumes"]
1562 @lazyproperty
1563 def hirshfeld_dipole(self):
1564 """The Hirshfeld systematic dipole moment for the chunk"""
1565 atoms = self._parse_atoms()
1567 if not np.any(atoms.pbc):
1568 return self._parse_hirshfeld()["dipole"]
1570 return None
1572 @lazyproperty
1573 def eigenvalues(self):
1574 """All outputted eigenvalues for the system"""
1575 return self._parse_eigenvalues()["eigenvalues"]
1577 @lazyproperty
1578 def occupancies(self):
1579 """All outputted occupancies for the system"""
1580 return self._parse_eigenvalues()["occupancies"]
1583def get_header_chunk(fd):
1584 """Returns the header information from the aims.out file"""
1585 header = []
1586 line = ""
1588 # Stop the header once the first SCF cycle begins
1589 while (
1590 "Convergence: q app. | density | eigen (eV) | Etot (eV)"
1591 not in line
1592 and "Convergence: q app. | density, spin | eigen (eV) |"
1593 not in line
1594 and "Begin self-consistency iteration #" not in line
1595 ):
1596 try:
1597 line = next(fd).strip() # Raises StopIteration on empty file
1598 except StopIteration:
1599 raise ParseError(
1600 "No SCF steps present, calculation failed at setup."
1601 )
1603 header.append(line)
1604 return AimsOutHeaderChunk(header)
1607def get_aims_out_chunks(fd, header_chunk):
1608 """Yield unprocessed chunks (header, lines) for each AimsOutChunk image."""
1609 try:
1610 line = next(fd).strip() # Raises StopIteration on empty file
1611 except StopIteration:
1612 return
1614 # If the calculation is relaxation the updated structural information
1615 # occurs before the re-initialization
1616 if header_chunk.is_relaxation:
1617 chunk_end_line = (
1618 "Geometry optimization: Attempting to predict improved coordinates."
1619 )
1620 else:
1621 chunk_end_line = "Begin self-consistency loop: Re-initialization"
1623 # If SCF is not converged then do not treat the next chunk_end_line as a
1624 # new chunk until after the SCF is re-initialized
1625 ignore_chunk_end_line = False
1626 while True:
1627 try:
1628 line = next(fd).strip() # Raises StopIteration on empty file
1629 except StopIteration:
1630 break
1632 lines = []
1633 while chunk_end_line not in line or ignore_chunk_end_line:
1634 lines.append(line)
1635 # If SCF cycle not converged or numerical stresses are requested,
1636 # don't end chunk on next Re-initialization
1637 patterns = [
1638 (
1639 "Self-consistency cycle not yet converged -"
1640 " restarting mixer to attempt better convergence."
1641 ),
1642 (
1643 "Components of the stress tensor (for mathematical "
1644 "background see comments in numerical_stress.f90)."
1645 ),
1646 "Calculation of numerical stress completed",
1647 ]
1648 if any(pattern in line for pattern in patterns):
1649 ignore_chunk_end_line = True
1650 elif "Begin self-consistency loop: Re-initialization" in line:
1651 ignore_chunk_end_line = False
1653 try:
1654 line = next(fd).strip()
1655 except StopIteration:
1656 break
1658 yield AimsOutCalcChunk(lines, header_chunk)
1661def check_convergence(chunks, non_convergence_ok=False):
1662 """Check if the aims output file is for a converged calculation
1664 Parameters
1665 ----------
1666 chunks: list of AimsOutChunks
1667 The list of chunks for the aims calculations
1668 non_convergence_ok: bool
1669 True if it is okay for the calculation to not be converged
1671 Returns
1672 -------
1673 bool
1674 True if the calculation is converged
1675 """
1676 if not non_convergence_ok and not chunks[-1].converged:
1677 raise ParseError("The calculation did not complete successfully")
1678 return True
1681@reader
1682def read_aims_output(fd, index=-1, non_convergence_ok=False):
1683 """Import FHI-aims output files with all data available, i.e.
1684 relaxations, MD information, force information etc etc etc."""
1685 header_chunk = get_header_chunk(fd)
1686 chunks = list(get_aims_out_chunks(fd, header_chunk))
1687 check_convergence(chunks, non_convergence_ok)
1689 # Relaxations have an additional footer chunk due to how it is split
1690 if header_chunk.is_relaxation:
1691 images = [chunk.atoms for chunk in chunks[:-1]]
1692 else:
1693 images = [chunk.atoms for chunk in chunks]
1694 return images[index]
1697@reader
1698def read_aims_results(fd, index=-1, non_convergence_ok=False):
1699 """Import FHI-aims output files and summarize all relevant information
1700 into a dictionary"""
1701 header_chunk = get_header_chunk(fd)
1702 chunks = list(get_aims_out_chunks(fd, header_chunk))
1703 check_convergence(chunks, non_convergence_ok)
1705 # Relaxations have an additional footer chunk due to how it is split
1706 if header_chunk.is_relaxation and (index == -1):
1707 return chunks[-2].results
1709 return chunks[index].results