NQCModels.jl

To perform nonadiabatic molecular dynamics simulations, it is necessary to define the system Hamiltonian. For simple models, this often comes in the form of a small matrix in the diabatic representation but equally the electronic Hamiltonian could be obtained directly from ab initio electronic structure theory.

NQCModels.jl is a package that aims to provide a common interface for defining these models that is flexible enough to allow for a wide range of specifications and requirements. NQCDynamics.jl uses this interface to obtain the potentials and couplings necessary to perform the dynamics simulations. Along with the minimal interface, NQCModels.jl also provides a library of popular models often used in the field of nonadiabatic dynamics.

Depending on the quantities provided by the Model, we use Julia's abstract type system to group models that provide the same quantities. Currently, there are two top-level abstract types: ClassicalModel and QuantumModel. The ClassicalModel is used for adiabatic dynamics, providing only a classical potential in the form of an analytical function and its derivative, the force used in classical mechanics. The QuantumModel is used for nonadiabatic dynamics, where the potential is no longer an analytical function but instead a Hermitian matrix.

Using Models

In the Getting started section we briefly touched on how the ClassicalModel works and introduced one of the included models. Here let's take a look at a QuantumModel, which is more appropriate for nonadiabatic dynamics.

The DoubleWell is a one-dimensional model with two electronic states, where each electronic state is harmonic with linear coupling to the single degree of freedom.

using NQCModels

model = DoubleWell()

DoubleWell{Int64, Int64, Int64, Int64}
  mass: Int64 1
  ω: Int64 1
  γ: Int64 1
  Δ: Int64 1

Our DoubleWell implements the functions potential and derivative, their in-place versions, nstates and ndofs that return the potential, the derivative of the potential, the number of states, and the number of degrees of freedom, respectively.

julia> p = potential(model, hcat(0.2))2×2 LinearAlgebra.Hermitian{Float64, Matrix{Float64}}:
 0.302843   0.5
 0.5       -0.262843
julia> d = derivative(model, hcat(0.2))1×1 Matrix{LinearAlgebra.Hermitian{Float64, Matrix{Float64}}}:
 [1.614213562373095 0.0; 0.0 -1.2142135623730952]
julia> potential!(model, p, hcat(0.5))2×2 Matrix{Float64}:
 0.832107   0.5
 0.5       -0.582107
julia> derivative!(model, d, hcat(0.5))
julia> nstates(model)2
julia> ndofs(model)1

In general the position argument that appears in the derivative and the potential will need to be provided as an AbstractMatrix with dimensions dofsxnatoms. In the case of 1D models of a single atom, potential and derivative accept a position of type Real, but for the most consistent results it's recommended that the user wraps real valued positions in 1x1 matrices using the hcat function.

To understand how this can extend to another dimension, we can take a quick look at the GatesHollowayElbow model which is another two state diabatic model, but this one uses two dimensions to model a diatomic molecule interacting with a surface. The two degrees of freedom are the molecular bond length and the distance from the surface respectively and so the potential and derivative functions accept a position argument with two values, one for each degree of freedom.

julia> model = GatesHollowayElbow()GatesHollowayElbow
  λ₁: Float64 3.5
  λ₂: Float64 3.5
  z₀: Float64 1.4
  x₀: Float64 0.6
  α: Float64 1.028
  d: Float64 0.17161933455967182
  z12: Float64 0.5
  c: Float64 0.018374661087759297
  γ: Float64 1.0
julia> p = potential(model, [1.0 1.0])2×2 LinearAlgebra.Hermitian{Float64, Matrix{Float64}}:
 0.0710215  0.0111448
 0.0111448  0.0744945
julia> d = derivative(model, [1.0 1.0])1×2 Matrix{LinearAlgebra.Hermitian{Float64, Matrix{Float64}}}:
 [0.0810696 0.0; 0.0 -0.0129425]  …  [-0.000787036 -0.0111448; -0.0111448 0.0810696]
julia> potential!(model, p, [1.2 0.5])2×2 Matrix{Float64}:
 0.0875049  0.0183747
 0.0183747  0.0295571
julia> derivative!(model, d, [1.2 0.5])1×2 Matrix{LinearAlgebra.Hermitian{Float64, Matrix{Float64}}}:
 [0.0728352 0.0; 0.0 -0.00642707]  …  [-0.00452908 -0.0183747; -0.0183747 0.0848168]
julia> nstates(model)2
julia> ndofs(model)2

Here we see how the derivative now becomes a Matrix with size matching our input, but each entry is a Hermitian containing the elementwise derivative of the potential with respect to each degree of freedom. In this case, the Matrix has size = (1, 2), and in general the dimensions of the Matrix that wraps the Hermitian derivatives should match the dimensions of the position Matrix a given model takes as input.

Included models

NQCModels.jl includes both analytical models as well as interfaces to connect to full-dimensional potential energy surfaces. An overview of the implemented analytical models can be seen in the Analytic model library and many shall return later when we investigate the dynamics methods. An overview of interfaces to full-dimensional potential energy surfaces can be seen in the Full dimensional model library.