Utilities
*********

Presented here is a small set of useful utility tools that can come in
handy in various contexts in your scripts. They are simple, standalone
objects always present in the main namespace.

* Periodic Table

* Units

* Geometry tools

* File format conversion tools

* Reaction energies

* Plotting tools

* Postprocess results

What is characteristic for the "PeriodicTable" and "Units" classes
described below is that they are meant to be used in a bit different
way than all other PLAMS classes. Usually one takes a class (like
"DiracJob"), creates an instance of it ("myjob = DiracJob(...)") and
executes some of its methods ("r = myjob.run()"). In contrast, utility
classes are designed in a way similar to the so-called singleton
design pattern. That means it is not possible to create instances of
these classes. The class itself serves for "one and only instance" and
all methods should be called using the class as the calling object:

   >>> x = PeriodicTable()
   PTError: Instances of PeriodicTable cannot be created
   >>> s = PeriodicTable.get_symbol(20)
   >>> print(s)
   Ca


Periodic Table
==============

Import path:

   scm.plams.tools.periodic_table

class PeriodicTable

   A singleton class for the periodic table of elements.

   For each element the following properties are stored: atomic
   symbol, atomic mass, atomic radius and number of connectors.

   Atomic mass is, strictly speaking, atomic weight, as present in
   Mathematica's ElementData function.

   Atomic radius and number of connectors are used by "guess_bonds()".
   Note that values of radii are neither atomic radii nor covalent
   radii. They are somewhat "emprically optimized" for the bond
   guessing algorithm.

   Note:

     This class is visible in the main namespace as both
     "PeriodicTable" and "PT".

   __init__()

      Return type:
         None

   classmethod get_atomic_number(symbol)

      Convert atomic symbol to atomic number.

      Parameters:
         **symbol** (*str*) --

      Return type:
         int

   classmethod get_symbol(atnum)

      Convert atomic number to atomic symbol.

      Parameters:
         **atnum** (*int*) --

      Return type:
         str

   classmethod get_mass(arg)

      Convert atomic symbol or atomic number to atomic mass.

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         float

   classmethod get_radius(arg)

      Convert atomic symbol or atomic number to radius.

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         float

   classmethod get_connectors(arg)

      Convert atomic symbol or atomic number to number of connectors.

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         int

   classmethod get_metallic(arg)

      Convert atomic symbol or atomic number to is_metallic.

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         int

   classmethod get_electronegative(arg)

      Convert atomic symbol or atomic number to is_electronegative.

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         int

   classmethod set_mass(element, value)

      Set the mass of *element* to *value*.

      Parameters:
         * **element** (*str*) --

         * **value** (*float*) --

      Return type:
         None

   classmethod set_radius(element, value)

      Set the radius of *element* to *value*.

      Parameters:
         * **element** (*str*) --

         * **value** (*float*) --

      Return type:
         None

   classmethod set_connectors(element, value)

      Set the number of connectors of *element* to *value*.

      Parameters:
         * **element** (*str*) --

         * **value** (*int*) --

      Return type:
         None

   classmethod _get_property(arg: str | int, prop: Literal[1, 2]) -> float
   classmethod _get_property(arg: str | int, prop: Literal[3, 4, 5]) -> int
   classmethod _get_property(arg: str | int, prop: Literal[6, 7]) -> float | None

      Get property of element described by either symbol or atomic
      number. Skeleton method for "get_radius()", "get_mass()" and
      "get_connectors()".

   classmethod get_electron_affinity(arg)

      Get the electron affinity of the metal element

      Note: Returns None if data is not available

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         float | None

   classmethod get_ionization_energy(arg)

      Get the electron affinity of the metal element

      Note: Returns None if data is not available

      Parameters:
         **arg** (*str** | **int*) --

      Return type:
         float | None


Units
=====

Import path:

   scm.plams.tools.units

class Units

   A singleton class for unit converter.

   All values are based on 2014 CODATA recommended values.

   The following constants and units are supported:

   * constants:

     * "speed_of_light" (also "c")

     * "electron_charge" (also "e")

     * "Avogadro_constant" (also "NA")

     * "Bohr_radius"

   * distance:

     * "Angstrom", "A", "Ang"

     * "Bohr", "au", "a.u."

     * "nm"

     * "pm"

   * reciprocal distance:

     * "1/Angstrom", "1/A", "Angstrom^-1", "A^-1",

     * "1/Bohr", "Bohr^-1"

   * angle:

     * "degree", "deg"

     * "radian", "rad"

     * "grad"

     * "circle"

   * charge:

     * "coulomb", "C"

     * "e"

   * energy:

     * "au", "a.u.", "Hartree"

     * "eV"

     * "kcal/mol"

     * "kJ/mol"

     * "cm^-1", "cm-1"

     * "K", "Kelvin"

     * "Hz", "Hertz"

     * "THz"

   * dipole moment:

     * "au", "a.u.", "e*bohr"

     * "Debye", "D"

     * All charge units multiplied by distance units, for example

     * "eA",  "e*A"

     * "Cm",  "C*m"

   * molecular polarizability:

     * "au", "a.u.", "(e*bohr)^2/hartree"

     * "e*A^2/V"

     * "C*m^2/V"

     * "cm^3"

     * "bohr^3"

     * "A^3", "angstrom^3", "Ang^3"

   * forces:

     * All energy units divided by angstrom or bohr, for example

     * "eV/angstrom"

     * "hartree/bohr"

   * hessian:

     * All energy units divided by angstrom^2 or bohr^2, for example

     * "eV/angstrom^2"

     * "hartree/bohr^2"

   * pressure:

     * All energy units divided by angstrom^3 or bohr^3, for example

     * "eV/angstrom^3"

     * "hartree/bohr^3"

     * And some more:

     * "Pa"

     * "GPa"

     * "bar"

     * "atm"

   Example:

      >>> print(Units.constants['speed_of_light'])
      299792458
      >>> print(Units.constants['e'])
      1.6021766208e-19
      >>> print(Units.convert(123, 'angstrom', 'bohr'))
      232.436313431
      >>> print(Units.convert([23.32, 145.0, -34.7], 'kJ/mol', 'kcal/mol'))
      [5.573613766730401, 34.655831739961755, -8.293499043977056]
      >>> print(Units.conversion_ratio('kcal/mol', 'kJ/mol'))
      4.184

   __init__()

      Return type:
         None

   classmethod conversion_ratio(inp, out)

      Return conversion ratio from unit *inp* to *out*.

      Parameters:
         * **inp** (*str*) --

         * **out** (*str*) --

      Return type:
         float

   classmethod convert(value, inp, out)

      Convert *value* from unit *inp* to *out*.

      *value* can be a single number or a container (list, tuple,
      numpy.array etc.). In the latter case a container of the same
      type and length is returned. Conversion happens recursively, so
      this method can be used to convert, for example, a list of lists
      of numbers, or any other hierarchical container structure.
      Conversion is applied on all levels, to all values that are
      numbers (also numpy number types). All other values (strings,
      bools etc.) remain unchanged. Note that type hints can be wrong
      for this method: *value* types with leaf type of int/np.int will
      be converted to float/np.float but annotated as int. *value*
      with containers with different depth will not give a useful
      typehint (but are converted correctly).

      Parameters:
         * **value** (*T*) --

         * **inp** (*str*) --

         * **out** (*str*) --

      Return type:
         *T*

   classmethod ascii2unicode(string)

      Converts '^2' to '²' etc., for prettier printing of units.

      Parameters:
         **string** (*str*) --

      Return type:
         str


Geometry tools
==============

A small module with simple functions related to 3D geometry
operations.

Import path:

   scm.plams.tools.geometry

rotation_matrix(vec1, vec2)

   Calculate the rotation matrix rotating *vec1* to *vec2*. Vectors
   can be any containers with 3 numerical values. They don't need to
   be normalized. Returns 3x3 numpy array.

   Parameters:
      * **vec1** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **vec2** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

   Return type:
      *ndarray*

axis_rotation_matrix(vector, angle, unit='radian')

   Calculate the rotation matrix rotating along the *vector* by
   *angle* expressed in *unit*.

   *vector* can be any container with 3 numerical values. They don't
   need to be normalized. A positive angle denotes counterclockwise
   rotation, when looking along *vector*. Returns 3x3 numpy array.

   Parameters:
      * **vector** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **angle** (*float*) --

      * **unit** (*str*) --

   Return type:
      *ndarray*

distance_array(array1, array2)

   Calculates distance between each pair of points in *array1* and
   *array2*. Returns 2D "numpy" array.

   Uses fast "cdist" function if "scipy" is present, otherwise falls
   back to slightly slower "numpy" loop. Arguments should be
   2-dimensional "numpy" arrays with the same second dimension. If
   *array1* is A x N and *array2* is B x N, the returned array is A x
   B.

   Parameters:
      * **array1** (*ndarray*) --

      * **array2** (*ndarray*) --

   Return type:
      *ndarray*

angle(vec1, vec2, result_unit='radian')

   Calculate an angle between vectors *vec1* and *vec2*.

   *vec1* and *vec2* should be iterable containers of length 3 (for
   example: tuple, list, numpy array). Values stored in them are
   expressed in Angstrom. Returned value is expressed in
   *result_unit*.

   This method requires all atomic coordinates to be numerical values,
   "TypeError" is raised otherwise.

   Parameters:
      * **vec1** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **vec2** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **result_unit** (*str*) --

   Return type:
      float

dihedral(p1, p2, p3, p4, unit='radian')

   Calculate the value of diherdal angle formed by points *p1*, *p2*,
   *p3* and *p4* in a 3D space. Arguments can be any containers with 3
   numerical values, also instances of "Atom". Returned value is
   always non-negative, measures the angle clockwise (looking along
   *p2-p3* vector) and is expressed in *unit*.

   Parameters:
      * **p1** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **p2** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **p3** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **p4** (*ndarray** | **Tuple**[**float**, **float**,
        **float**] **| **Sequence**[**float**]*) --

      * **unit** (*str*) --

   Return type:
      *ndarray*

cell_shape(lattice)

   Converts lattice vectors to lengths and angles (in radians) Sets
   internal cell size data, based on set of cell vectors.

   *lattice* is list containing three cell vectors (a 3x3 matrix)

   Parameters:
      **lattice** (*ndarray** |
      **Sequence**[**Sequence**[**float**]**]*) --

   Return type:
      *List*[float] | None

cellvectors_from_shape(box)

   Converts lengths and angles (in radians) of lattice vectors to the
   lattice vectors.

   *box* should be an iterable of length 3, containing the lengths a,
   b, c, or an iterable of length 6 additionally containing the angles
   alpha, beta, gamma

   Parameters:
      **box** (*Sequence**[**float**]*) --

   Return type:
      *List*[*List*[float]]


File format conversion tools
============================

A small module for converting VASP output to AMS-like output, and for
converting ASE .traj trajectory files to the .rkf format.

Import path:

   scm.plams.tools.converters

traj_to_rkf(trajfile, rkftrajectoryfile, task=None, timestep=0.25)

   Convert ase .traj file to .rkf file. NOTE: The order of atoms (or
   the number of atoms) cannot change between frames!

   trajfile : str
      path to a .traj file

   rkftrajectoryfile : str
      path to the output .rkf file (will be created)

   task : str
      Which task to write. If None it is auto-determined.

   timestep: float
      Which timestep to write when task == 'moleculardynamics'

   Returns : 2-tuple (coords, cell)
      The final coordinates and cell in angstrom

   Parameters:
      * **trajfile** (*str*) --

      * **rkftrajectoryfile** (*str*) --

      * **task** (*str** | **None*) --

      * **timestep** (*float*) --

   Return type:
      *Tuple*[ndarray | None, Cell | None]

vasp_output_to_ams(vasp_folder, wdir=None, overwrite=False, write_engine_rkf=True, task=None, timestep=0.25)

   Converts VASP output (OUTCAR, ...) to AMS output (ams.rkf,
   vasp.rkf)

   Returns: a string containing the directory where ams.rkf was
   written

   vasp_folder : str
      path to a directory with an OUTCAR, INCAR, POTCAR etc. files

   wdir : str or None
      directory in which to write the ams.rkf and vasp.rkf files If
      None, a subdirectory "AMSJob" of vasp_folder will be created

   overwrite : bool
      if False, first check if wdir already contains ams.rkf and
      vasp.rkf, in which case do nothing if True, overwrite if exists

   write_engine_rkf : bool
      If True, also write vasp.rkf alongside ams.rkf. The vasp.rkf
      file will only contain an AMSResults section (energy, gradients,
      stress tensor). It will not contain the DOS or the band
      structure.

   task : str
      Which task to write to ams.rkf. If None it is auto-determined
      (probably set to 'geometryoptimization')

   timestep : float
      If task='moleculardynamics', which timestep (in fs) between
      frames to write

   Parameters:
      * **vasp_folder** (*str*) --

      * **wdir** (*str** | **None*) --

      * **overwrite** (*bool*) --

      * **write_engine_rkf** (*bool*) --

      * **task** (*str** | **None*) --

      * **timestep** (*float*) --

   Return type:
      str

qe_output_to_ams(qe_outfile, wdir=None, overwrite=False, write_engine_rkf=True)

   Converts a qe .out file to ams.rkf and qe.rkf.

   Returns: a string containing the directory where ams.rkf was
   written

   If the filename ends in .out, check if a .results directory exists.
   In that case, place the AMSJob subdirectory in the .results
   directory.

   Otherwise, create a new directory called filename.AMSJob

   qe_outfile : str
      path to the qe output file

   Parameters:
      * **qe_outfile** (*str*) --

      * **wdir** (*str** | **None*) --

      * **overwrite** (*bool*) --

      * **write_engine_rkf** (*bool*) --

   Return type:
      str

gaussian_output_to_ams(outfile, wdir=None, overwrite=False, write_engine_rkf=True)

   Converts a Gaussian .out file to ams.rkf and gaussian.rkf.

   Returns: a string containing the directory where ams.rkf was
   written

   If the filename ends in .out, check if a .results directory exists.
   In that case, place the AMSJob subdirectory in the .results
   directory.

   Otherwise, create a new directory called filename.AMSJob

   outfile : str
      path to the gaussian output file

   Parameters:
      * **outfile** (*str*) --

      * **wdir** (*str** | **None*) --

      * **overwrite** (*bool*) --

      * **write_engine_rkf** (*bool*) --

   Return type:
      str

rkf_to_ase_traj(rkf_file, out_file, get_results=True)

   Convert an ams.rkf trajectory to a different trajectory format
   (.xyz, .traj, anything supported by ASE)

   rkf_file: str
      Path to an ams.rkf file

   out_file: str
      Path to the .traj or .xyz file that will be created. If the file
      exists it will be overwritten. If a .xyz file is specified it
      will use the normal ASE format (not the AMS format).

   get_results: bool
      Whether to include results like energy, forces, and stress in
      the trajectory.

   Parameters:
      * **rkf_file** (*str*) --

      * **out_file** (*str*) --

      * **get_results** (*bool*) --

   Return type:
      *List*[Atoms]

rkf_to_ase_atoms(rkf_file, get_results=True)

   Convert an ams.rkf trajectory to a list of ASE atoms

   rkf_file: str
      Path to an ams.rkf file

   get_results: bool
      Whether to include results like energy, forces, and stress in
      the trajectory.

   Returns: a list of all the ASE Atoms objects.

   Parameters:
      * **rkf_file** (*str*) --

      * **get_results** (*bool*) --

   Return type:
      *List*[Atoms]

file_to_traj(outfile, trajfile)

   outfile : str
      path to existing file (OUTCAR, qe.out, etc.)

   trajfile : str
      will be created

   Parameters:
      * **outfile** (*str*) --

      * **trajfile** (*str*) --

   Return type:
      str


Reaction energies
=================

*New in AMS2026*: The "balance" function is new in AMS2026. For usage,
see the Balance reaction equations example.

Import paths:

   scm.plams.tools.reaction
   scm.plams.tools.reaction_energies

balance(reactants, products, method='plams')

   Method to find balanced reaction converting reactants into
   products.

   reactants: sequence of str | Species
      The reactants. If string, should be the summary formula like
      "C2H6O". Note: an element may only appear once in the string,
      "HCOOH" is not allowed.

      To set charges or min_coeff, use a Species.

   products: sequence of str | Species
      The products.

   method: "plams" or "sympy" (default "plams")
      Which method to use. "sympy" requires that sympy is installed.

   Returns: ReactionEquation
      A balanced reaction equation.

   Example:

   >>> from scm.plams.tools.reaction import balance, Species
   >>> balance(reactants=["CH4", "O2"], products=["CO2", "H2O"])
   >>> balance(reactants=[Species("H", charge=1), Species("OH", charge=-1)], products=[Species("H2O")])

   Parameters:
      * **reactants** (*Sequence**[**str** | **Species**]*) --

      * **products** (*Sequence**[**str** | **Species**]*) --

      * **method** (*Literal**[**'plams'**, **'sympy'**]*) --

   Return type:
      *ReactionEquation*

class ReactionEquation(reactants, products, method='plams')

   Class representing a reaction, which can then be balanced to any
   selected set of molecules

   Parameters:
      * **reactants** (*Sequence**[**Molecule** | **str**]*) --

      * **products** (*Sequence**[**Molecule** | **str**]*) --

      * **method** (*str*) --

   __init__(reactants, products, method='plams')

      Initiate the reaction

      * "reactants" -- List of PLAMS molecules or strings representing
        a molecular formula

      * "products" -- List of PLAMS molecules or strings representing
        a molecular formula

      Parameters:
         * **reactants** (*Sequence**[**Molecule** | **str**]*) --

         * **products** (*Sequence**[**Molecule** | **str**]*) --

         * **method** (*str*) --

   prepare_state()

      Do the time consuming stuff that needs to be done before we call
      balance

      Return type:
         None

   balance(min_coeffs=None)

      Balance the equation to any set of molecules specified in
      min_coeffs

      * "min_coeffs" -- Vector representing the minimal allowable
        (integer) coefficient of each molecule
           [0, 0, 1, 0, 0]

      Parameters:
         **min_coeffs** (*Iterable**[**int**] **| **None*) --

      Return type:
         *ndarray*[*Any*, *dtype*[*ScalarType*]] | None

   set_minimum_coefficients(min_coeffs=None)

      Set the minimum coefficient values

      * "min_coeffs" -- Vector representing the minimal allowable
        (integer) coefficient of each molecule
           [0, 0, 1, 0, 0]

      Parameters:
         **min_coeffs** (*ndarray**[**Any**,
         **dtype**[**ScalarType**]**] **| **None*) --

      Return type:
         None

   get_nullspace_basis()

      Here we try to get the nullspace vectors

      Return type:
         *ndarray*[*Any*, *dtype*[*ScalarType*]]

   get_reduced_basis()

      Select only the rows in basis that share the same block with the
      main product

      Return type:
         *ndarray*[*Any*, *dtype*[*ScalarType*]]

   setup_optimizer(basis)

      Solve the problem using ILP.

      min(sum(x)), with basis@y = x, and x >= min_coeffs

      * "basis" -- m x n matrix, representing the null space vectors
        as rows

      Parameters:
         **basis** (*ndarray**[**Any**, **dtype**[**ScalarType**]**]*)
         --

      Return type:
         None

   optimize_coefficients()

      Solve the problem using ILP.

      min(sum(x)), with basis@y = x, and x >= min_coeffs

      Return type:
         *ndarray*[*Any*, *dtype*[*ScalarType*]]

   refine_optimization()

      Find all solutions with the same cost, and select the one with
      the lowest indices

      Note: This is done by repressing the highest coefficient in the
      previous solution,
         and optimizing again, untill no improvement is made.

      Return type:
         *ndarray*[*Any*, *dtype*[*ScalarType*]]

   property reactants: List[str]

      Get the list of reactant formulas

   property products: List[str]

      Get the list of product formulas

   set_charges(reactant_charges, product_charges)

      Set the charges for reactant and product molecules

      * "reactant_charges" -- List of charges of the reactant
        molecules

      * "product_charges" -- List of charges of the product molecules

      Parameters:
         * **reactant_charges** (*List**[**float**]*) --

         * **product_charges** (*List**[**float**]*) --

      Return type:
         None

   set_coefficients(icoeffs)

      Set one of the entries in self.equivalent_coefficients as the
      final result

      Parameters:
         **icoeffs** (*int*) --

      Return type:
         None

   __format__(fmt=None)

      Formats the chemical reaction.

      fmt: str
         If 'smiles', will print the reaction using SMILES if
         possible.

      Example:

      >>> print(f"{reaction:smiles}")

      Parameters:
         **fmt** (*str** | **None*) --

      Return type:
         str

   __str__()

      Write the balanced reaction

      Return type:
         str

   property reaction_charge: int | str

      Write the balanced reaction

   matrix_as_string(mat=None, space=8)

      Print a numpy matrix in nice format

      Parameters:
         * **mat** (*ndarray**[**Any**, **dtype**[**ScalarType**]**]
           **| **None*) --

         * **space** (*int*) --

      Return type:
         str

   _set_elements()

      Set the element data from the formulas

      Return type:
         None

   _set_matrix()

      Get the matrix from the molecular formulas

      Return type:
         None

   static _elements_from_formula(formula)

      Get all elements from the formula

      Parameters:
         **formula** (*str*) --

      Return type:
         *Dict*[str, int]

   __weakref__

      list of weak references to the object (if defined)

Older functions:

balance_equation(reactants, products, normalization='r0', normalization_value=1.0)

   Calculate stoichiometric coefficients This only works if *
   number_of_chemical_elements == len(reactants)+len(products), OR *
   number_of_chemical_elements == len(reactants)+len(products)-1

   Returns: a 2-tuple (coeffs_reactants, coeffs_products)
      coeffs_reactants is a list with length == len(reactants)
      coeffs_products is a list with length == len(products)

   reactants: a list of amsjobs, or a list of paths to ams.results
   folders or ams.rkf files or .xyz files, or a list of Molecules, or
   a list of stoichiometry dicts, or a list of Molecules
      The reactants

   products: a list of amsjobs, or a list of paths to ams.results
   folders or ams.rkf files, or a list of Molecules or .xyz files, or
   a list of stoichiometry dicts, or a list of Molecules
      The products

   normalization: str
      'r0' for the first reactant, 'r1' for the second reactant, etc.
      'p0' for the first product, 'p1' for the second product, etc.
      This normalizes the chemical equation such that the coefficient
      in front of the specified species is normalization_value

   normalization_value : float
      The coefficient to normalize to

   EXAMPLE:

      balance_equation(
          reactants=[
              {'N': 2, 'H': 8, 'Cr': 2, 'O': 7}
          ],
          products=[
              {'Cr': 2, 'O': 3},
              {'N': 2},
              {'H': 2, 'O': 1}
          ])

   The above returns a tuple "([1.0], [1.0, 1.0, 1.0, 4.0])"

   Parameters:
      * **reactants** (*Sequence**[**AMSJob** | **Molecule** | **str**
        | **Dict**[**str**, **int**]**]*) --

      * **products** (*Sequence**[**AMSJob** | **Molecule** | **str**
        | **Dict**[**str**, **int**]**]*) --

      * **normalization** (*str*) --

      * **normalization_value** (*float*) --

   Return type:
      *Tuple*[*List*[float], *List*[float]]

reaction_energy(reactants, products, normalization='r0', unit='hartree')

   Calculates a reaction energy from an unbalanced chemical equation
   (the equation is first balanced)

   reactants: a list of amsjobs or paths to ams results folders,
      The recatnts

   products: a list of amsjobs or paths to ams results folders
      The products

   normalization: str
      normalize the chemical equation by setting the corresponding
      coefficient to 1.

      * 'r0': first reactant

      * 'r1': second reactant, ...

      * 'p0: first product,

      * 'p1': second product, ...

   unit: str
      Unit of the reaction energy

   Returns: a 3-tuple (coeffs_reactants, coeffs_products,
   reaction_energy)

   Parameters:
      * **reactants** (*Sequence**[**AMSJob** | **str**]*) --

      * **products** (*Sequence**[**AMSJob** | **str**]*) --

      * **normalization** (*str*) --

      * **unit** (*str*) --

   Return type:
      *Tuple*[*List*[float], *List*[float], float]


Plotting tools
==============

See also:

  * Band structure with BAND

Tools for creating plots with matplotlib.

Import path:

   scm.plams.tools.plot

The "scm.plams.tools.plot" module also contains small reusable helpers
for common analysis tasks. For example, "linear_fit_extrapolate_to_0"
performs a linear regression and returns the fitted line extended to
"x = 0".

Example:

   from scm.plams.tools.plot import linear_fit_extrapolate_to_0

   fit_x, fit_y, slope, intercept = linear_fit_extrapolate_to_0(
       [1.0, 2.0, 3.0],
       [3.0, 5.0, 7.0],
   )

linear_fit_extrapolate_to_0(x, y)

   Perform a linear regression on "x" and "y" and return the fit
   extended to "x = 0".

   x: sequence of float
      X values for the linear regression.

   y: sequence of float
      Y values for the linear regression.

   Returns: tuple
      "fit_x", "fit_y", "slope", "intercept".

   If "0" is already present in "x", it is not appended a second time.

   Parameters:
      * **x** (*Sequence**[**float**]*) --

      * **y** (*Sequence**[**float**]*) --

   Return type:
      *Tuple*[*ndarray*, *ndarray*, float, float]

plot_band_structure(x, y_spin_up, y_spin_down=None, labels=None, fermi_energy=None, zero=None, ax=None)

   Plots an electronic band structure from DFTB, BAND, or
   QuantumEspresso engines with matplotlib.

   To control the appearance of the plot you need to call
   "plt.ylim(bottom, top)", "plt.title(title)", etc. manually outside
   this function.

   x: list of float
      Returned by AMSResults.get_band_structure()

   y_spin_up: 2D numpy array of float
      Returned by AMSResults.get_band_structure()

   y_spin_down: 2D numpy array of float. If None, the spin down bands
   are not plotted.
      Returned by AMSResults.get_band_structure()

   labels: list of str
      Returned by AMSResults.get_band_structure()

   fermi_energy: float
      Returned by AMSResults.get_band_structure(). Should have the
      same unit as "y".

   zero: None or float or one of 'fermi', 'vbm', 'vbmax', 'cbm',
   'cbmin'
      Shift the curves so that y=0 is at the specified value. If None,
      no shift is performed. 'fermi', 'vbmax', and 'cbmin' require
      that the "fermi_energy" is not None. Note: 'vbmax' and 'cbmin'
      calculate the zero as the highest (lowest) eigenvalue smaller
      (greater) than or equal to "fermi_energy". This is NOT
      necessarily equal to the valence band maximum or conduction band
      minimum as calculated by the compute engine.

   Additional parameters:

   "ax": matplotlib axis
      The axis. If None, one will be created

   Parameters:
      * **x** (*List**[**float**]*) --

      * **y_spin_up** (*ndarray*) --

      * **y_spin_down** (*ndarray** | **None*) --

      * **labels** (*List**[**str**] **| **None*) --

      * **fermi_energy** (*float** | **None*) --

      * **zero** (*Literal**[**'fermi'**, **'vbm'**, **'vbmax'**,
        **'cbm'**, **'cbmin'**] **| **float** | **None*) --

      * **ax** (*plt.Axes** | **None*) --

   Return type:
      plt.Axes

plot_phonons_band_structure(x, y, labels=None, zero=None, ax=None)

   Plots a phonons band structure from DFTB, BAND or QuantumEspresso
   engines with matplotlib.

   To control the appearance of the plot you need to call
   "plt.ylim(bottom, top)", "plt.title(title)", etc. manually outside
   this function.

   x: list of float
      Returned by AMSResults.get_phonons_band_structure()

   y: 2D numpy array of float
      Returned by AMSResults.get_phonons_band_structure()

   labels: list of str
      Returned by AMSResults.get_phonons_band_structure()

   zero: None or float
      Shift the curves so that y=0 is at the specified value. If None,
      no shift is performed.

   Additional parameters:

   "ax": matplotlib axis
      The axis. If None, one will be created

   Parameters:
      * **x** (*List**[**float**]*) --

      * **y** (*ndarray*) --

      * **labels** (*List**[**str**] **| **None*) --

      * **zero** (*float** | **None*) --

      * **ax** (*plt.Axes** | **None*) --

   Return type:
      plt.Axes

plot_phonons_dos(energy, total_dos, dos_per_species, dos_per_atom, dos_type='total', ax=None)

   Plots the phonons DOS from DFTB, BAND or QuantumEspresso engines
   with matplotlib.

   To control the appearance of the plot you need to call
   "plt.ylim(bottom, top)", "plt.title(title)", etc. manually outside
   this function.

   energy: list of float
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   total_dos: list of float
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   dos_per_species: dictionary of list of float
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   dos_per_atom: dictionary of list of float
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   dos_type: str
      Specifies the kind of plot to show. Possible options:
         * "total": Total DOS.

         * "species": DOS decomposed by species.

         * "atom": DOS decomposed by atom.

   Additional parameters:

   "ax": matplotlib axis
      The axis. If None, one will be created

   Parameters:
      * **energy** (*List**[**float**]*) --

      * **total_dos** (*List**[**float**]*) --

      * **dos_per_species** (*Dict**[**str**, **List**[**float**]**]*)
        --

      * **dos_per_atom** (*Dict**[**str**, **List**[**float**]**]*) --

      * **dos_type** (*str*) --

      * **ax** (*plt.Axes** | **None*) --

   Return type:
      plt.Axes

plot_phonons_thermodynamic_properties(temperature, properties, units, ax=None)

   Plots the phonons thermodynamic properties from DFTB, BAND or
   QuantumEspresso engines with matplotlib.

   To control the appearance of the plot you need to call
   "plt.ylim(bottom, top)", "plt.title(title)", etc. manually outside
   this function.

   temperature: list of float
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   properties: dictionary of list of float
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   units: dictionary of str
      Returned by AMSResults.get_phonons_thermodynamic_properties()

   Additional parameters:

   "ax": matplotlib axis
      The axis. If None, one will be created

   Parameters:
      * **temperature** (*List**[**float**]*) --

      * **properties** (*Dict**[**str**, **List**[**float**]**]*) --

      * **units** (*Dict**[**str**, **str**]*) --

      * **ax** (*plt.Axes** | **None*) --

   Return type:
      plt.Axes

plot_molecule(molecule, figsize=None, ax=None, keep_axis=False, **kwargs)

   Show a molecule in a Jupyter notebook

   Parameters:
      * **molecule** (*ase.Atoms** | **Molecule*) --

      * **figsize** (*int** | **Tuple**[**int**] **| **None*) --

      * **ax** (*plt.Axes** | **None*) --

      * **keep_axis** (*bool*) --

      * **kwargs** (*Any*) --

   Return type:
      plt.Axes

plot_image_grid(images, rows=None, cols=None, figsize=None, show_labels=True, save_path=None)

   Plot a dictionary of images in a matplotlib grid.

   Parameters:
      * **images** (*Dict**[**str**, **PilImage.Image**]*) --
        dictionary with labels as keys and images as values; iteration
        order determines image order in the grid

      * **rows** (*int** | **None*) -- number of rows in the grid; if
        "None", infer from "cols" and number of images

      * **cols** (*int** | **None*) -- number of columns in the grid;
        if "None", infer from "rows" and number of images

      * **figsize** (*Tuple**[**float**, **float**] **| **None*) --
        matplotlib figure size; if "None", uses a grid-proportional
        default

      * **show_labels** (*bool*) -- whether to show labels above
        images; labels are taken from dictionary keys

      * **save_path** (*str** | **PathLike** | **None*) -- optional
        path to save the plotted grid image using matplotlib "savefig"

   Returns:
      2D numpy array of matplotlib axes with shape "(rows, cols)"

   Return type:
      np.ndarray

plot_correlation(job1, job2, section, variable, alt_section=None, alt_variable=None, file='ams', multiplier=1.0, unit=None, save_txt=None, ax=None, show_xy=True, show_linear_fit=True, show_mad=True, show_rmsd=True, xlabel=None, ylabel=None)

   Plot a correlation plot from AMS .rkf files

   job1: AMSJob or List[AMSJob]
      Job(s) plotted on x-axis

   job2: AMSJob or List[AMSJob]
      job2: Job(s) plotted on y-axis

   section: str
      section: section to read on .rkf files

   variable: str
      variable: variable to read

   alt_section: str
      Section to read on .rkf files for job2. If not specified it will
      be the same as "section"

   alt_variable : str
      Variable to read for job2. If not specified it will be the same
      as "variable".

   file: str, optional
      file: "ams" or "engine", defaults to "ams"

   multiplier: float, optional
      multiplier: Numbers will be multiplied by this number, defaults
      to 1.0

   unit: str, optional
      unit: unit will be shown in the plot, defaults to None

   save_txt: str, optional
      save_txt: If not None, save the xy data to this text file,
      defaults to None

   ax: matplotlib axis, optional
      ax: matplotlib axis, defaults to None

   show_xy: bool, optional
      show_xy: Whether to show y=x line, defaults to True

   show_linear_fit: bool, optional
      show_linear_fit: Whether to perform and show a linear fit,
      defaults to True

   show_mad: bool, optional
      show_mad: Whether to show mean absolute deviation, defaults to
      True

   show_rmsd: bool, optional
      show_rmsd: Whether to show root-mean-square deviation, defaults
      to True

   xlabel: str, optional
      xlabel: The x-label. If not given will be a list of job names,
      defaults to None

   ylabel: str, optional
      ylabel: THe y-label. If not given will be al ist of job names,
      defaults to None

   Returns: A matplotlib axis

   Parameters:
      * **job1** (*AMSJob** | **List**[**AMSJob**]*) --

      * **job2** (*AMSJob** | **List**[**AMSJob**]*) --

      * **section** (*str*) --

      * **variable** (*str*) --

      * **alt_section** (*str** | **None*) --

      * **alt_variable** (*str** | **None*) --

      * **file** (*str*) --

      * **multiplier** (*float*) --

      * **unit** (*str** | **None*) --

      * **save_txt** (*str** | **None*) --

      * **ax** (*plt.Axes** | **None*) --

      * **show_xy** (*bool*) --

      * **show_linear_fit** (*bool*) --

      * **show_mad** (*bool*) --

      * **show_rmsd** (*bool*) --

      * **xlabel** (*str** | **None*) --

      * **ylabel** (*str** | **None*) --

   Return type:
      plt.Axes

plot_msd(job, start_time_fit_fs=None, ax=None)

   job: AMSMSDJob
      The job for which to plot the results

   start_time_fit_fs: float
      The start time (in fs) for which to perform the linear fit

   ax: matplotlib axis
      The axis. If None, one will be created

   Returns: matplotlib axis

   Parameters:
      * **job** (*AMSMSDJob*) --

      * **start_time_fit_fs** (*float** | **None*) --

      * **ax** (*plt.Axes** | **None*) --

   Return type:
      plt.Axes

plot_work_function(coordinate, planarAverage, macroscopicAverage, Efermi, Vbulk, Vvacuum, WF, ax=None)

   Plots an Electrostatic Potential Profile from AMS-QE with
   matplotlib.

   To control the appearance of the plot you need to call
   "plt.ylim(bottom, top)", "plt.title(title)", etc. manually outside
   this function.

   "coordinate": 1D array of float.
      Returned by AMSResults.get_work_function_results().

   "planarAverage": 1D array of float.
      Returned by AMSResults.get_work_function_results().

   "macroscopicAverage": 1D array of float.
      Returned by AMSResults.get_work_function_results(). Should have
      the same unit as "planarAverage".

   "Efermi": float.
      Returned by AMSResults.get_work_function_results(). Should have
      the same unit as "planarAverage".

   "Vbulk": float.
      Returned by AMSResults.get_work_function_results(). Should have
      the same unit as "planarAverage".

   "Vvacuum": Tuple[float,float].
      Returned by AMSResults.get_work_function_results(). Should have
      the same unit as "planarAverage".

   "WF": Tuple[float,float].
      Returned by AMSResults.get_work_function_results(). Should have
      the same unit as "planarAverage".

   Additional parameters:

   "ax": matplotlib axis
      The axis. If None, one will be created

   Returns: matplotlib axis

   Parameters:
      * **coordinate** (*ndarray*) --

      * **planarAverage** (*ndarray*) --

      * **macroscopicAverage** (*ndarray*) --

      * **Efermi** (*float*) --

      * **Vbulk** (*float*) --

      * **Vvacuum** (*Tuple**[**float**, **float**]*) --

      * **WF** (*Tuple**[**float**, **float**]*) --

      * **ax** (*plt.Axes** | **None*) --

   Return type:
      plt.Axes

plot_grid_molecules(molecules, legends=None, molsPerRow=2, subImgSize=(200, 200), ax=None, save_svg_path=None, **kwargs)

   Plot series of molecules in a grid using RDKit

   Parameters:
      * **ax** (*plt.Axes** | **None*) -- if provided molecules are
        plotted in these axes, note that the quality of the image
        might reduce, defaults to None

      * **save_svg_path** (*str** | **PathLike** | **None*) --
        pathlike of the file, with formats .svg to save it as image,
        it returns the svg string, defaults to None

      * **molecules** (*List**[**Molecule** | **ChemicalSystem**]*) --

      * **legends** (*List**[**str**] **| **None*) --

      * **molsPerRow** (*int*) --

      * **subImgSize** (*Tuple**[**int**, **int**]*) --

      * **kwargs** (*Any*) --

   Returns:
      an image of the molecules

   Return type:
      pil.Image or plt.Axes or string


Postprocess results
===================

Tools for postprocessing the results.

Import path:

   scm.plams.tools.postprocess_results

The "scm.plams.tools.postprocess_results" module contains helpers such
as "moving_average" for smoothing paired "x"/"y" data and
"broaden_results" for constructing broadened spectra.

Example:

   from scm.plams.tools.postprocess_results import moving_average

   avg_x, avg_y = moving_average([1.0, 2.0, 3.0], [3.0, 5.0, 7.0], window=2)

moving_average(x, y, window)

   Calculate a moving average of x and y.

   Parameters:
      * **x** (*Sequence**[**float**]*) -- X values

      * **y** (*Sequence**[**float**]*) -- Y values

      * **window** (*int*) -- Moving-average window size

   Returns:
      "x_moving_averaged", "y_moving_averaged"

   Return type:
      Tuple[np.ndarray, np.ndarray]

broaden_results(centers, areas, broadening_width=40, broadening_type='gaussian', x_data=(0, 4000, 0.5), post_process=None)

   convenient function for create xy curve with gaussian or lorentzian
   peaks. Used to obtain IR or Raman spectra for example.

   Parameters:
      * **centers** (*np.ndarray*) -- x positions of the centers of
        the peaks

      * **areas** (*np.ndarray*) -- areas of the peaks

      * **broadening_width** (*Union**[**float**, **np.ndarray**]**,
        **optional*) -- if is np.ndarray each peak broadening is
        assigned otherwise apply the same value for all the peaks,
        defaults to 40

      * **broadening_type** (*Literal**[**'gaussian'**,
        **'lorentzian'**]**, **optional*) -- the line shape of the
        peak, defaults to "gaussian"

      * **x_data** (*Union**[**np.ndarray**, **Tuple**[**float**,
        **float**, **float**]**]**, **optional*) -- it can be a
        np.ndarray or you can set a tuple with (min,max,step_size),
        defaults to (0, 4000, 0.5)

      * **post_process** (*Literal**[**'max_to_1'**]**, **optional*)
        -- if max_to_1 the resulted spectrum has the max peak height
        =1, defaults to None

   Returns:
      the x and y arrays

   Return type:
      Tuple[np.ndarray, np.ndarray]
