Utilities¶

Periodic Table¶

class PeriodicTable[source]¶

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__()[source]¶

Initialize self. See help(type(self)) for accurate signature.

classmethod get_atomic_number(symbol)[source]¶

Convert atomic symbol to atomic number.

classmethod get_symbol(atnum)[source]¶

Convert atomic number to atomic symbol.

classmethod get_mass(arg)[source]¶

Convert atomic symbol or atomic number to atomic mass.

classmethod get_radius(arg)[source]¶

Convert atomic symbol or atomic number to radius.

classmethod get_connectors(arg)[source]¶

Convert atomic symbol or atomic number to number of connectors.

classmethod get_metallic(arg)[source]¶

Convert atomic symbol or atomic number to number of connectors.

classmethod get_electronegative(arg)[source]¶

Convert atomic symbol or atomic number to number of connectors.

classmethod set_mass(element, value)[source]¶

Set the mass of element to value.

classmethod set_radius(element, value)[source]¶

Set the radius of element to value.

classmethod set_connectors(element, value)[source]¶

Set the number of connectors of element to value.

classmethod _get_property(arg, prop)[source]¶

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

Units¶

class Units[source]¶

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

• 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__()[source]¶

Initialize self. See help(type(self)) for accurate signature.

classmethod conversion_ratio(inp, out)[source]¶

Return conversion ratio from unit inp to out.

classmethod convert(value, inp, out)[source]¶

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.

classmethod ascii2unicode(string)[source]¶

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

Geometry tools¶

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

rotation_matrix(vec1, vec2)[source]¶

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.

axis_rotation_matrix(vector, angle, unit='radian')[source]¶

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.

distance_array(array1, array2)[source]¶

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.

angle(vec1, vec2, result_unit='radian')[source]¶

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.

dihedral(p1, p2, p3, p4, unit='radian')[source]¶

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.

cell_shape(lattice)[source]¶

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

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

cellvectors_from_shape(box)[source]¶

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

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.

traj_to_rkf(trajfile, rkftrajectoryfile, task=None, timestep=0.25)[source]¶

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

trajfilestr

path to a .traj file

rkftrajectoryfilestr

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

taskstr

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

timestep: float

Which timestep to write when task == ‘moleculardynamics’

Returns2-tuple (coords, cell)

The final coordinates and cell in angstrom

vasp_output_to_ams(vasp_folder, wdir=None, overwrite=False, write_engine_rkf=True, task=None, timestep=0.25)[source]¶

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

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

vasp_folderstr

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

wdirstr 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

overwritebool

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_rkfbool

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.

taskstr

Which task to write to ams.rkf. If None it is auto-determined (probably set to ‘geometryoptimization’)

timestepfloat

If task=’moleculardynamics’, which timestep (in fs) between frames to write

qe_output_to_ams(qe_outfile, wdir=None, overwrite=False, write_engine_rkf=True)[source]¶

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_outfilestr

path to the qe output file

gaussian_output_to_ams(outfile, wdir=None, overwrite=False, write_engine_rkf=True)[source]¶

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

outfilestr

path to the gaussian output file

rkf_to_ase_traj(rkf_file, out_file, get_results=True)[source]¶

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.

rkf_to_ase_atoms(rkf_file, get_results=True)[source]¶

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.

file_to_traj(outfile, trajfile)[source]¶

outfilestr

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

trajfilestr

will be created

Reaction energies¶

balance_equation(reactants, products, normalization='r0', normalization_value=1.0)[source]¶

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_valuefloat

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])

reaction_energy(reactants, products, normalization='r0', unit='hartree')[source]¶

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)

Plotting tools¶

Tools for creating plots with matplotlib.

plot_band_structure(x, y_spin_up, y_spin_down=None, labels=None, fermi_energy=None, zero=None, show=False)[source]¶

Plots an electronic band structure from DFTB or BAND 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’, ‘vbmax’, ‘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.

show: bool

If True, call plt.show() at the end

plot_molecule(molecule, figsize=None, ax=None, keep_axis=False, **kwargs)[source]¶

Show a molecule in a Jupyter notebook

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)[source]¶

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_variablestr

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

Postprocess results¶

Tools for postprocessing the results.

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

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 hight =1, defaults to None

Returns

the x and y arrays

Return type

Tuple[np.ndarray, np.ndarray]