Kinetica.jl API - ASE Interface
Optimisation and Property Calculation
Basic Species Properties
Kinetica.get_mult
— Functionget_mult(sd::SpeciesData, sid)
Calculates the spin multiplicity of the species in sd
at ID sid
.
Kinetica.get_mult!
— Functionget_mult!(sd::SpeciesData, sid)
Caches the spin multiplicity of the species in sd
at ID sid
in sd.cache[:mult]
.
Kinetica.get_charge
— Functionget_charge(sd::SpeciesData, sid)
Calculates the charge of the species in sd
at ID sid
.
Kinetica.get_charge!
— Functionget_charge!(sd::SpeciesData, sid)
Caches the charge of the species in sd
at ID sid
in sd.cache[:charge]
.
Kinetica.get_formal_charges
— Functionget_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.
Kinetica.get_formal_charges!
— Functionget_formal_charges!(sd::SpeciesData, sid)
Caches the formal charges of the species in sd
at ID sid
in sd.cache[:formal_charges]
.
Kinetica.get_initial_magmoms
— Functionget_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.
Kinetica.get_initial_magmoms!
— Functionget_initial_magmoms!(sd::SpeciesData, sid)
Caches the initial magnetic moments of the species in sd
at ID sid
in sd.cache[:initial_magmoms]
.
Kinetica.correct_magmoms_for_mult!
— Functioncorrect_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.
Kinetica.get_hydrogen_idxs
— Functionget_hydrogen_idxs(amsmi::String)
Returns indices of hydrogen atoms in atom-mapped SMILES amsmi
.
Geometry Optimisation
Kinetica.geomopt!
— Functiongeomopt!(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.
Kinetica.kabsch_fit!
— Functionkabsch_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'.
Kinetica.permute_hydrogens!
— Functionpermute_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.
NEB
Kinetica.get_initial_sys_mult
— Functionget_initial_sys_mult(mults)
Returns the combined spin multiplicity of a system of molecules.
Kinetica.get_rxn_mult
— Functionget_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.
Kinetica.neb
— Functionneb(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.
Kinetica.highest_energy_frame
— Functionhighest_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.
Vibrational Analysis
Kinetica.calc_species_vibrations!
— Functioncalc_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
.
Kinetica.calc_ts_vibrations!
— Functioncalc_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
.