Getting started

To get started with the package we can identify the necessary ingredients to perform a simple classical dynamics simulation and walk through how to set up the simulation.

Atoms

First, we must define the particles in the simulation. For this purpose we provide the Atoms type which will contain the symbols, atomic numbers and masses for our atoms. Technically these need not be actual atoms and be a generic particle.

If using real atoms, then they can be constructed using the chemical symbols as a Vector of Julia's Symbol types, a Vector{Symbol}:

julia> Atoms([:H, :C])Atoms{Float64}([:H, :C], [1, 6], [1837.4715941070515, 21894.713607956142])

You can see that this contains two atoms labelled by their atomic numbers with their masses in atomic units.

Alternatively, if not using real atoms, Atoms can be created using a Vector{<:Real} where the provided numbers are the masses of the particles.

julia> Atoms([1, 2, 3, 4, 5, 6])Atoms{Float64}([:X, :X, :X, :X, :X, :X], [0, 0, 0, 0, 0, 0], [1.0, 2.0, 3.0, 4.0, 5.0, 6.0])

A more detailed look into the Atoms type along with a description of how to save and load structures can be found here.

Representing atomic positions and velocities

This package chooses to separate the dynamical variables from the static atomic parameters included in the Atoms type. This allows us to easily interface with other numerical packages like DifferentialEquations.jl and AdvancedMH.jl. As such, both positions and velocities are represented using Julia's standard Array type, specifically as an Array{T,2} or the Matrix{T} type, which are equivalent. If you are new to Julia, you can find a description of the Array here. The first dimension contains each atomic degree of freedom, and the second dimension contains each atom. For example, a 3D system with two atoms would have positions:

julia> using Symbolics
julia> @variables x1, y1, z1, x2, y2, z26-element Vector{Symbolics.Num}:
 x1
 y1
 z1
 x2
 y2
 z2
julia> r = [x1 x2;
            y1 y2;
            z1 z2]3×2 Matrix{Symbolics.Num}:
 x1  x2
 y1  y2
 z1  z2

For a 1D system it would be necessary to create a 1x1 matrix:

julia> r = fill(x1, (1,1))1×1 Matrix{Symbolics.Num}:
 x1

Velocities are handled in the same way as positions and the data structures are the same. Usually manual initialisation like this will only be necessary for small model systems, whereas full dimensional model system will be read from a file instead. This is explored in the Atoms documentation.

Models

The next ingredient required to set up the simulation is the Model, i.e., the potentials in which the system evolves. These Models are provided by NQCModels.jl, which is a convenient infrastructure for defining different kinds of models for adiabatic and nonadiabatic dynamics. These models can range from simple analytic potentials all the way up to multi-dimensional ab initio potentials. Refer to the NQCModels.jl page for information on the available models and a description of how to implement further models.

For now we can look at an AdiabaticModel which provides a simple harmonic potential energy function.

julia> model = Harmonic()Harmonic{Float64, Float64, Float64}
  m: Float64 1.0
  ω: Float64 1.0
  r₀: Float64 0.0
  dofs: Int64 1

Here, the four parameters (m, ω, r₀ and dofs) for this model are shown along with their types and default values. These values can be modified by specifying a new value in the constructor. For example for m:

julia> model = Harmonic(m=0.4)Harmonic{Float64, Float64, Float64}
  m: Float64 0.4
  ω: Float64 1.0
  r₀: Float64 0.0
  dofs: Int64 1

Adiabatic models implement two functions to calculate the total energy and the forces, respectively: potential(model, R) and derivative(model, R).

Let's try these out and take a look at the results:

julia> potential(model, hcat(25.0))125.0
julia> derivative(model, hcat(25.0))1×1 Matrix{Float64}:
 10.0

To make sure the model is what we expect, we can plot the potential and derivative using a custom plotting recipe. This looks pretty harmonic to me!

using Plots

plot(-5:0.1:5, model)

Simulation

To control all simulation parameters in one environment, we use the Simulation type which will contain both the Atoms and Models explained above, along with any extra information required for the simulation.

julia> sim = Simulation{Classical}(Atoms(:H), model)Simulation{Classical}:
  Atoms{Float64}([:H], [1], [1837.4715941070515])
  Harmonic{Float64, Float64, Float64}
  m: Float64 0.4
  ω: Float64 1.0
  r₀: Float64 0.0
  dofs: Int64 1

Here, we have specified that each atom has a single degree of freedom and have not provided a simulation cell. Classical is a type parameter, and specifies the dynamics method that we want to use. Check out Dynamics methods to learn about the other kinds of dynamics available.

Dynamics variables

The final ingredient before we can perform our simulation is the initial positions and velocities of our particles. For each dynamics type, the method DynamicsVariables is implemented and creates the dynamics variables for us. For classical dynamics we must provide a Matrix of velocities and of positions. These should have size = (dofs, natoms), matching the arguments of the potential and derivative functions. Usually the initial coordinates would have some physical significance, perhaps sampled from a relevant distribution, but here we use random numbers for simplicity.

julia> v = rand(3, 3);
julia> r = rand(3, 3);
julia> DynamicsVariables(sim, v, r)ComponentVector{Float64}(v = [0.7277317953613319 0.04824897817757423 0.40924034366121953; 0.08290567051933806 0.08688717663882295 0.06102669702332608; 0.6004318020844844 0.9413022025909804 0.48106060841060816], r = [0.011434033128359222 0.9075505187571333 0.7074184691304342; 0.9808948922440597 0.23659946209668048 0.16607480522185025; 0.22719741361836676 0.8047771979672188 0.8957794560820922])

Bringing it all together

We have now covered all the parts necessary to perform our first classical dynamics simulation. Let's quickly set up our simulation parameters using what we've learned. Here we'll have two atoms in a harmonic potential, each with a single degree of freedom.

julia> atoms = Atoms([:H, :C])Atoms{Float64}([:H, :C], [1, 6], [1837.4715941070515, 21894.713607956142])
julia> sim = Simulation{Classical}(atoms, Harmonic(ω=50.0))Simulation{Classical}:
  Atoms{Float64}([:H, :C], [1, 6], [1837.4715941070515, 21894.713607956142])
  Harmonic{Float64, Float64, Float64}
  m: Float64 1.0
  ω: Float64 50.0
  r₀: Float64 0.0
  dofs: Int64 1
julia> z = DynamicsVariables(sim, randn(size(sim)), randn(size(sim)))ComponentVector{Float64}(v = [-0.5743454183642588 -0.49473532674148296], r = [0.8643469945629969 0.4357225485794407])

Now, we can finally run the trajectory using the run_dynamics function. This takes three positional arguments: the simulation parameters sim, the time span we want to solve fortspan, and the dynamics variablesz. For classical dynamics we also provide a timestepdtsince we're using theVelocityVerlet` algorithm by default.

The final keyword argument output is used to specify the quantities we want to save during the dynamics. A list of the available quantities can be found here.

tspan = (0.0, 50.0)
solution = run_dynamics(sim, (0.0, 50.0), z;
                                   dt=0.1, output=(OutputPosition, OutputVelocity))

3-element Dictionaries.Dictionary{Symbol, Any}
           :Time │ [0.0, 0.1, 0.2, 0.30000000000000004, 0.4, 0.5, 0.6, 0.7, 0.7…
 :OutputPosition │ [[0.8643469945629969 0.4357225485794407], [0.801032450300036…
 :OutputVelocity │ [[-0.5743454183642588 -0.49473532674148296], [-0.68763828845…

Here you can see the output containing the time steps and the output quantities we specified. These can be accessed directly as shown here:

julia> solution[:Time]501-element Vector{Float64}:
  0.0
  0.1
  0.2
  0.30000000000000004
  0.4
  0.5
  0.6
  0.7
  0.7999999999999999
  0.8999999999999999
  ⋮
 49.20000000000043
 49.30000000000043
 49.40000000000043
 49.50000000000043
 49.600000000000435
 49.700000000000436
 49.80000000000044
 49.90000000000044
 50.0
julia> solution[:OutputPosition]501-element Vector{Matrix{Float64}}:
 [0.8643469945629969 0.4357225485794407]
 [0.8010324503000366 0.38600025577283453]
 [0.7268193368718573 0.3358372170141637]
 [0.642717372111189 0.2852907098724078]
 [0.5498708181256797 0.2344184497715349]
 [0.4495429128652832 0.18327852408940687]
 [0.3430986829456781 0.13192932583197692]
 [0.23198637157047353 0.08042948695851257]
 [0.11771773423743546 0.02883781143397496]
 [0.0018474703185093887 -0.02278679191500276]
 ⋮
 [0.17800650970048928 0.8985404815410549]
 [0.06278139486082421 0.9397734982094069]
 [-0.05329790182119111 0.9799334550572049]
 [-0.16865204577170134 1.0189744963158294]
 [-0.2817115685929301 1.0568520438259612]
 [-0.3909382217048202 1.0935228479381776]
 [-0.4948459052206822 1.128945036896619]
 [-0.5920208873058432 1.1630781646493382]
 [-0.6811410389213527 1.1958832570305988]

As with the models, we provide custom plotting recipes to quickly visualise the results before performing further analysis by manually accessing the fields of the solution table. To use these recipes, simply provide the solution to the plot function from Plots.jl and give the name of the output quantity as the second argument. This will only work if this quantity was specified in run_dynamics.

plot(solution, :OutputPosition)
plot!(solution, :OutputVelocity)

Ensemble simulations

We have shown how to perform a single trajectory, but usually we are interested in performing many and calculating observables using statistical methods. Running more trajectories is as simple as providing the trajectories keyword to run_dynamics, but we'll go through this in more detail in the Ensemble simulations section.

What's next?

Now that we've covered the basics of classical dynamics, we're ready to explore the world of nonadiabatic dynamics. All the dynamics methods follow these patterns and anything you find elsewhere in the documentation should now seem relatively familiar.