Introduction

Performing dynamics simulations is at the core of this package's functionality (as you might have guessed from the name). This section of the documentation will describe how to perform dynamics simulations, building on the introduction from Getting started.

Since we use DifferentialEquations to perform the dynamics, it is most natural to split up the system parameters from the dynamics variables. This manifests itself as two separate data types: the Simulation, and the DynamicsVariables.

The Simulation holds all the static information about the system: the atoms, the model, the temperature, the cell and the dynamics method.

atoms = Atoms(2000) # Single atom with mass = 2000 a.u.
sim = Simulation{Ehrenfest}(atoms, TullyModelOne(); temperature=0, cell=InfiniteCell())

Simulation{Ehrenfest{Float64}}:
  Atoms{Float64}([:X], [0], [2000.0])
  TullyModelOne{Float64, Float64, Float64, Float64}
  a: Float64 0.01
  b: Float64 1.6
  c: Float64 0.005
  d: Float64 1.0

Here we have initialised the simulation parameters, including the default temperature and cell explicitly. sim takes the place of the p parameter seen throughout DifferentialEquations.

For DifferentialEquations to allow for a wide variety of solvers, the input arrays DynamicsVariables must be AbstractArrays. In nonadiabatic dynamics simulations, we usually have different groups of variables that behave in particular ways. For example: for mapping variable methods we have positions, velocities, and two sets of mapping variables representing the electronic degrees of freedom.

For this purpose we use the ComponentVector, which allows us to arbitrarily partition the variables into their subgroups. This allows us to keep all the variables in a single array as required by DifferentialEquations, whilst still having them partitioned for convenient computation and readable code.

v0 = hcat(10) / 2000
r0 = hcat(-5)
u0 = DynamicsVariables(sim, v0, r0, PureState(1))

ComponentVector{Float64}(v = [0.005;;], r = [-5.0;;], σreal = [1.0 0.0; 0.0 0.0], σimag = [0.0 0.0; 0.0 0.0])

Since each dynamics method has a different set of variables, each method implements DynamicsVariables(sim, args...), which will convert the input into the correct structure. This helps to ensure each method follows a similar workflow, making it easy to switch between different methods. The output of this function takes the place of the u argument seen throughout DifferentialEquations.

With both the Simulation and DynamicsVariables in hand, the central function is run_dynamics which allows us to perform a single dynamics trajectory. run_dynamics takes the simulation parameters sim and the initial conditions u0, along with a time span tspan that the trajectory will cover.

tspan = (0.0, 2000.0)
run_dynamics(sim, tspan, u0; output=OutputDynamicsVariables, dt=1.0)

2-element Dictionaries.Dictionary{Symbol, Any}
                    :Time │ [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …
 :OutputDynamicsVariables │ ComponentArrays.ComponentVector{Float64, Vector{Flo…

The output is a dictionary containing entries for :Time and our requested output quantity. Output is a required keyword and the code will error unless at least one quantity is specified. By passing a Tuple to the output keyword argument we can ask for multiple quantities.

out = run_dynamics(sim, tspan, u0; output=(OutputPosition, OutputAdiabaticPopulation), dt=1.0)

3-element Dictionaries.Dictionary{Symbol, Any}
                      :Time │ [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0…
            :OutputPosition │ [[-5.0;;], [-4.99500000134185;;], [-4.99000000538…
 :OutputAdiabaticPopulation │ [[1.0, 0.0], [1.0, 8.315490106702257e-27], [1.0, …

The quantities that are available are listed here. More quantities can be added by defining new functions with the signature f(sol, i). The first argument is the DifferentialEquations.jl solution object and the second is the trajectory index.

This time we can see that the output contains the two quantities that we asked for.

using Plots
plot(out, :OutputPosition)

plot(out, :OutputAdiabaticPopulation)

All of the dynamics methods work in a similar way. For details on a specific method along with examples, please see the method specific page in the sidebar under Dynamics methods.