NQCDistributions.jl

Storing and sampling distributions

In order to perform ensembles of trajectories, it is useful to have a convenient way to generate distributions of velocities and positions which can be sampled to initialise trajectories. The NQCDistributions.jl package contains the types and functions that seek to address this requirement as painlessly as possible.

For quantum classical nonadiabatic dynamics simulations, the initial distributions contain both nuclear and electronic degrees of freedom.

This page describes the types that can be used to represent nuclear and electronic distributions and demonstrates how they can be combined into a product distribution.

Nuclear Distributions using DynamicalDistribution

When handling distributions for the nuclear degrees of freedom, the DynamicalDistribution type can be used to store initial velocities and positions:

d = DynamicalDistribution(10, 5, (3, 2))

Here, we have created a delta distribution with fixed velocities and positions, the final argument specifies the size of each sample. The (3, 2) case shown here would be appropriate when using 2 atoms each with 3 degrees of freedom.

julia> rand(d)ComponentVector{Int64}(v = [10 10; 10 10; 10 10], r = [5 5; 5 5; 5 5])

DynamicalDistribution is flexible and each of the first two arguments can be Real, Vector or Sampleable.

Each of these options can be composed in any combination, let's see the case where we have an analytic distribution of positions and a preset collection of velocities:

using Distributions

velocity = [[1.0 1.0;1.0 1.0], [2.0 2.0; 2.0 2.0], [3.0 3.0; 3.0 3.0]]
position = Normal()
d = DynamicalDistribution(velocity, position, (2, 2))
rand(d)

ComponentVector{Float64}(v = [3.0 3.0; 3.0 3.0], r = [-0.6637635515551271 0.7174105899241786; 1.3632084599735597 0.6787430624233933])

This has generated normally distributed positions along with one of the three velocities we provided.

Sampling the nuclear distribution

To learn how to generate configurations to use with the DynamicalDistribution, read on to the next sections about the included sampling methods.

VelocityBoltzmann

When performing equilibrium simulations it is often desirable to initialise trajectories with thermal velocities, e.g. in combination with positions obtained from Monte Carlo sampling. These can be obtained for each atom from a gaussian distribution of the appropriate width, or alternatively, using the VelocityBoltzmann distribution which simplifies the process. This takes the temperature, masses and size of the system and ensures the samples you obtain are of the correct shape:

using NQCDynamics
using Unitful

velocity = VelocityBoltzmann(300u"K", rand(10), (3, 10))
rand(velocity)

3×10 Matrix{Float64}:
 -0.0289863  0.0070886    0.844678   …  -0.0271776   -0.0596305  -0.2134
 -0.0288267  0.00754776   0.0764973      0.0828133    0.0105265  -0.0552835
  0.0132569  0.0434114   -0.12115        0.00928299  -0.021202    0.0364977

This can be handed directly to the DynamicalDistribution when Boltzmann velocities are required.

distribution = DynamicalDistribution(velocity, 1, (3, 10))
rand(distribution)

ComponentVector{Float64}(v = [0.01045141631449122 -0.04103561924858895 … -0.009574444812328197 0.09911929979647943; 0.044282957435553154 0.0015077015538464314 … -0.05889224258003931 -0.22990558804002448; -0.05655743861896711 -0.02074325824216273 … 0.013502184413835624 0.021595630553922114], r = [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0])

Wigner distributions

For harmonic oscillator systems, we have implemented the analytic Wigner distribution. These are just mormal distributions of the appropriate width but can be accessed easily as in the following:

julia> using NQCDistributions
julia> omega = 1.0;
julia> beta = 1e-3;
julia> mass = 10;
julia> dist = PositionHarmonicWigner(omega, beta, mass, centre=0.0)Distributions.Normal{Float64}(μ=0.0, σ=10.000000416666651)
julia> rand(dist)16.243185980500716
julia> dist = VelocityHarmonicWigner(omega, beta, mass)Distributions.Normal{Float64}(μ=0.0, σ=10.000000416666651)
julia> rand(dist)6.655784390086933

These can also be given to the DynamicalDistribution since they are just univariate normal distributions.

Nuclear distributions for Ring polymer simulations

The components making up a DynamicalDistribution can all be adapted for use in Ring Polymer simulations.

For samplable components based on nuclear configurations, simply use three-dimensional arrays instead of two-dimensional ones.

Pre-defined distribution functions such as VelocityBoltzmann can be turned into a three-dimensional version using the RingPolymerWrapper:

velocity = VelocityBoltzmann(300u"K", rand(10), (3, 10))
velocity_ring_polymer = RingPolymerWrapper(velocity, n_beads, Int[]) # RingPolymerWrapper(Distribution, number of ring-polymer beads, indices of atoms to treat classically)

Electronic distributions

For nonadiabatic dynamics, the initial electronic variables must also be sampled. For this, we can use an ElectronicDistribution which will tell our simulation how we want to sample the initial variables. Currently, two of these are provided, the PureState and the MixedState. The PureState is used for nonequilibrium simulations when the population is confined to a single state, whereas MixedState allows for a mixed state distribution.

julia> using NQCDistributions
julia> PureState(1, Diabatic())PureState{Diabatic}(1, Diabatic())
julia> PureState(2, Adiabatic())PureState{Adiabatic}(2, Adiabatic())
julia> MixedState([1, 2], Diabatic())MixedState{Int64, Diabatic}([1, 2], Diabatic())

These structs contain only the minimal information about the distributions, whereas the sampling of the distribution is handled separately by each of the different methods.

Product distributions - Combining Nuclear and electronic distributions

The initial nuclear and electronic states of a system can be combined in a product distribution, which can be used in place of a purely nuclear distribution for methods considering electronic dynamics.

nuclear_dist = DynamicalDistribution(velocity, position, (2, 2))
electronic_dist = PureState(2, Adiabatic())

total_dist = nuclear_dist * electronic_dist

rand(total_dist.nuclear) # Returns a random nuclear configuration

total_dist.electronic.state # Returns the chosen electronic state.