Contributing a new method

A key goal of NQCDynamics.jl is to provide an accessible toolkit for implementing new nonadiabatic dynamics methods. This page details the steps you must take in order to create a new dynamics method. The existing methods are stored inside the DynamicsMethods submodule, and this is where new methods should be implemented. Technically, it is possible to implement new methods completely separately from the package by importing and extending the relevant functions but if you would like to include your method in the package, it should be added within this submodule.

Basic implementation

To implement a new dynamics method, the necessary steps are:

Create a new subtype of DynamicsMethods.Method.

This Method acts as an extra parameter inside the simulation that allows us to specify any extra information needed for our dynamics method. This can be a good place to include any temporary arrays and parameters for the simulation. Refer to the Julia manual section on Composite Types to learn how this types are created. Here, our type is called MyMethod and we have included an a parameter that will influence our dynamics:

using NQCDynamics

struct MyMethod <: DynamicsMethods.Method
    a::Float64
end

Implement DynamicsVariables

This function returns an AbstractArray of the variables to be used as the initial condition for the simulation. The array should contain all of the variables that will change during the dynamics.

For classical dynamics, this would include only the positions and velocities, however, for FSSH we must also include the continuous electronic variables and the discrete state.

For our new method, MyMethod, we implement the DynamicsVariables function and return a ComponentVector containing the velocities, positions and extra variables x. Inside this function we are free to take any inputs and manipulate them before returning the result. As an example, suppose that our x variables are randomly generated each time we run the dynamics, this could be done as follows:

using ComponentArrays

function DynamicsMethods.DynamicsVariables(sim::Simulation{<:MyMethod}, v, r, k)
    return ComponentVector(v=v, r=r, x=rand()*k)
end

Here, we take the velocities v, positions r and assign them to the output as we would for classical dynamics, but we also generate a random between 0 and k, where k was given as input.

Implement motion!(du, u, sim, t)

This function should fill du with the time-derivative of the dynamics variables u in the usual way expected by DifferentialEquations.jl. We use the in-place version, where each element of du is filled with the time derivative of the correponding element in u.

Inside the DynamicsUtils submodule there are some useful functions like velocity! and divide_by_mass! which can handle some of the common parts of the motion! function. You are free to perform whatever manipulations you like inside this function, but note that motion! is a performance critical function, called numerous times during the simulation, so you should attempt to minimise allocations inside this function.

function DynamicsMethods.motion!(du, u, sim::Simulation{MyMethod}, t)

    DynamicsUtils.velocity!(du.r, u.v, u.r, sim, t) # Set du.r equal to the velocity

    # Set the acceleration of the particles
    du.v .= -sim.method.a .* u.r # Use the `a` parameter we stored in the `method`.
    DynamicsUtils.divide_by_mass!(du.v, sim.atoms.masses) # Divide du.v by the mass

    du.x .= 1 ./ u.x # Set time derivative of `x`.

    return nothing # The return of this function is not used so the return is unimportant
end

Here we have set the time derivative of the positions equal to the velocity, the time derivative of the velocities equal to the acceleration where the force involves the parameter a. Finally, the time derivative of the extra x variable is also set.

Solve a trajectory

To perform a simulation with our new method, we can write a script in the usual format and run the dynamics. In this script, we have a single atom with a mass of 1 with a single degree of freedom. We match this by initialising the positions and velocities equal to random 1x1 matrices. The a parameter of the method has been set equal to 2.0, and the initial value of x has been set to 0.5.

sim = Simulation(Atoms(1), Free(), MyMethod(2.0))
u = DynamicsVariables(sim, rand(1,1), rand(1,1), [0.5])

sol = run_dynamics(sim, (0.0, 5.0), u, output=(OutputPosition, OutputVelocity, OutputDynamicsVariables))

4-element Dictionaries.Dictionary{Symbol, Any}
                    :Time │ [0.0, 0.05496382758024514, 0.08043646561728238, 0.1…
          :OutputPosition │ [[0.5049040961991769;;], [0.5275445981903814;;], [0…
          :OutputVelocity │ [[0.4400970698257821;;], [0.3833211646560538;;], [0…
 :OutputDynamicsVariables │ ComponentArrays.ComponentVector{Float64, Vector{Flo…

To visualise the result we can plot each of the quantities from the output table:

using Plots

plot(sol, :OutputPosition, label="Position")
plot!(sol, :OutputVelocity, label="Velocity")
plot!(sol, :OutputDynamicsVariables, label="Dynamics variables", legend=true)
ylabel!("Value(t)")

Advanced tips

Is there a custom algorithm you can implement?

Some dynamics methods have special algorithms that are tailored to the specific problem and achieve better performance than the general algorithms include in DifferentialEquations.jl. For example, ring polymer methods typically use a symplectic scheme to solve for the internal modes of the ring polymer, allowing much larger timesteps. The DifferentialEquations.jl framework provides a simple interface for adding new algorithms, check out the developer documentation to learn how it works. You can also find some examples of custom algorithms in the DynamicsMethods.IntegrationAlgorithms module.