Kinetica.jl API - ASE Interface

get_mult(sd::SpeciesData, sid)

Calculates the spin multiplicity of the species in sd at ID sid.

get_mult!(sd::SpeciesData, sid)

Caches the spin multiplicity of the species in sd at ID sid in sd.cache[:mult].

get_charge(sd::SpeciesData, sid)

Calculates the charge of the species in sd at ID sid.

get_charge!(sd::SpeciesData, sid)

Caches the charge of the species in sd at ID sid in sd.cache[:charge].

get_formal_charges(amsmi::String)
get_formal_charges(sd::SpeciesData, sid)

Calculates formal charges on each atom of a species.

Can be given an atom-mapped SMILES amsmi to calculate from, or a species in sd at ID sid, in which case the atom-mapped SMILES is calculated on-the-fly.

get_formal_charges!(sd::SpeciesData, sid)

Caches the formal charges of the species in sd at ID sid in sd.cache[:formal_charges].

get_initial_magmoms(amsmi::String)
get_initial_magmoms(sd::SpeciesData, sid)

Calculates initial magnetic moments on each atom of a species.

Can be given an atom-mapped SMILES amsmi to calculate from, or a species in sd at ID sid, in which case the atom-mapped SMILES is calculated on-the-fly.

get_initial_magmoms!(sd::SpeciesData, sid)

Caches the initial magnetic moments of the species in sd at ID sid in sd.cache[:initial_magmoms].

correct_magmoms_for_mult(reac_magmoms::Vector{Float64}, prod_magmoms::Vector{Float64}, mult::Int)

Identifies a set of initial magnetic moments that are consistent with spin multiplicity mult across a whole reaction.

As single-reference electronic structure methods cannot handle switching of electronic states along a NEB path, a consistent multiplicity must be used. In cases where initial magnetic moments are needed, these must be consistent with the mult for both the reactant and product systems.

Attempts to correct existing magmoms by initially flipping lone radical electrons, e.g. by converting a dissociated system of two monoradicals to be a spin up-spin down pair with total mult of 1 to match with a singlet-state reactant. If this cannot be achieved, resorts to flipping half of an electron pair, e.g. by converting singlet carbene to triplet carbene.

get_hydrogen_idxs(amsmi::String)

Returns indices of hydrogen atoms in atom-mapped SMILES amsmi.

geomopt!(sd::SpeciesData, i, calc_builder[, calcdir::String="./", optimiser="LBFGSLineSearch", 
         fmax=0.01, maxiters=nothing, check_isomorphic=true, kwargs...])
geomopt!(frame::Dict{String, Any}, calc_builder[, calcdir::String="./", mult::Int=1, 
         chg::Int=0, formal_charges=nothing, initial_magmoms=nothing, 
         optimiser="BFGSLineSearch", fmax=0.01, maxiters=1000, check_isomorphic=true, 
         kwargs...])

Runs an ASE-driven geometry optimisation of the species in frame.

Can be run directly from a frame, or a frame can be extracted from sd.xyz[i] in theh case of the second method. With this method, formal charges, total charge and spin multiplicity are assumed to have been calculated and cached in sd.cache. When running directly from a frame, this information must be passed manually.

calc_builder should be a struct with a functor that returns a correctly constructed ASE calculator for the system at hand. While ASE can handle many system-specific calculator details from an Atoms object, quantities such as spin multiplicity and sometimes charge must be input separately. For this reason, the calc_builder functor must take mult::Int and charge::Int as its first two arguments. Any other arguments can be passed via this method's kwargs.

Some ASE calculators handle charged species at the Atoms level. These require an array of formal charges on each atom to be given during Atoms construction. If formal_charges is provided, this will occur. If not, all formal charges will be assumed to be zero.

IMPORTANTLY, any formal_charges array MUST match the atom ordering of the provided frame. This can by calling get_formal_charges on an atom-mapped SMILES from atom_map_smiles.

This optimisation always uses one of ASE's optimisers out of "BFGSLineSearch", "fire", "bfgs" or "lbfgs". The maximum force fmax and maximum number of iterations maxiters that this optimiser converges with can also be controlled.

Directly modifies the atomic positions and energy of the passed in frame. Energies returned are in eV. Returns a boolean for whether the optimisation was a conv.

kabsch_fit!(frame1::Dict{String, Any}, frame2::Dict{String, Any})

Computes the maximum overlap of atoms within frame1 to frame2.

Modifies the atomic positions in frame1 to be as close as possible to those in frame2 through a combination of translation and rotation. Uses the Kabsch algorithm, as implemented in the Python package 'rmsd'.

permute_hydrogens!(frame1::Dict{String, Any}, hidxs::Vector{Vector{Int}}, frame2::Dict{String, Any})

Corrects erroneous RDKit atom-mapping on hydrogen atoms.

Due to the implicit handling of hydrogens in RDKit, they can sometimes receive incorrect atom mapping duing conversion of a geometry to atom-mapped SMILES. This can cause additional bond breakage or rotation over the course of a NEB calculation.

This is fixed by permuting every pair of hydrogen atoms in frame1 (usually a reactant geometry) and comparing whether the positional RMSD between it and frame2 (a product geometry) decreases. If this is the case, the permutation is made permanent.

This process repeats over all hydrogens in a system until no more permutations are accepted, indicating the atom mapping with the least likelihood to cause NEB problems has been found.

get_initial_sys_mult(mults)

Returns the combined spin multiplicity of a system of molecules.

get_rxn_mult(n_reacs::Int, rmult::Int, n_prods::Int, pmult::Int)
get_rxn_mult(reacsys::Dict{String, Any}, prodsys::Dict{String, Any})

Chooses a spin multiplicity to apply over a whole reaction given reactant and product multiplicities.

Given an addition (2->1) or a dissociation (1->2), picks the mult on the side with the fewest species, as this usually results in the species that are expected to be more stable (reactants for dissociation, products for association) being in their most stable electronic configuration.

Given a substitution (2->2) or a rearrangement (1->1), the mult should always be balanced. However, abstractions from radicals can lead to unbalanced mult, e.g. [CH2]C + [H] –> C=C + [H][H], in which case the lower mult should be taken. There are also rare cases involving carbenes that may be unbalanced, e.g. [CH]C –> C=C, in which case the lower mult should also be taken.

neb(reacsys, prodsys, calc::ASENEBCalculator[, calcdir="./", kwargs...])

Interpolates and runs (CI-)NEB for the reaction defined by endpoints reacsys and prodsys.

Given ExtXYZ frames for reactants and products reacsys and prodsys respectively, performs an interpolation to generate an initial reaction path and then optimises it with ASE's NEB implementation.

Interpolation scheme is selected by calc.interpolation, where "linear" and "idpp" are the currently available options.

NEB optimisation is performed by the optimiser defined in calc.neb_optimiser. If calc.climb=false, runs a regular NEB calculation until maximum force is less than calc.ftol. If calc.climb=true, runs a regular NEB calculation until maximum force is less than calc.climb_ftol, then enables CINEB and runs until maximum force is less than calc.ftol. Each optimisation only runs until calc.maxiters iterations have elapsed.

kwargs are passed to the Kinetica calculator builder.

Returns the optimised NEB path as a Python list of ASE Atoms objects.

highest_energy_frame(images::Py)

Finds the highest energy NEB image in images, returns as a frame.

images should be a Python list of ASE Atoms objects representing a NEB path (see the neb(reacsys, prodsys, calc) method). Locates the highest energy image, converts it to an ExtXYZ frame and returns this frame.

calc_species_vibrations(sd::SpeciesData, sid, calc_builder[, calcdir="./", refresh=false, delta=0.01, ivetol=0.1, kwargs...])

Calculates vibrational energies of a species.

Given a species with an optimised geometry in sd.xyz[sid], uses ASE to calculate vibrational energies with a finite difference approximation of its Hessian, using the ASE calculator constructed through calc_builder.

Writes the vibrational energies in eV as an array into sd.cache[:vib_energies]. Defaults to not overwriting already calculated vibrational energies for a species, instead returnong immediately. If recalculations are desired, use refresh=true.

If a species has imaginary frequencies, these can be ignored below a vibrational energy threshold ivetol, which removes them from the final vibrational energy array. All imaginary frequencies can be ignored by passing ivetol=0.0.

calc_ts_vibrations(ts_cache::Dict{Symbol, Any}, rid, calc_builder[, calcdir="./", delta=0.01, ivetol=0.1, kwargs...])

Calculates vibrational energies of a transition state.

Given a TS with an optimised geometry in ts_cache[:xyz][rid], uses ASE to calculate vibrational energies with a finite difference apporoximation of its Hessian, using the ASE calculator constructed through calc_builder.

Writes the vibrational energies in eV as an array into ts_cache[:vib_energies].

If a species has imaginary frequencies, these can be ignored below a vibrational energy threshold ivetol, which removes them from the final vibrational energy array. All imaginary frequencies can be ignored by passing ivetol=0.0.