Nonadiabatic ring polymer molecular dynamics (NRPMD)

Theory

Nonadiabatic ring polymer molecular dynamics (NRPMD) is a method that uses the ring polymer formalism to include quantum effects in the nuclear dynamics and mapping variables for the electronic degrees of freedom. ([9], [10], [11]) This results in a classical dynamics in an extended phasespace of the ring polymer with each bead coupled to a set of classical mapping variables. Originally, this method was proposed as a simple combination of the Meyer-Miller-Stock-Thoss mapping formalism with RPMD but has since been rigorously derived from nonadiabatic Matsubara dynamics ([12]).

The classical Hamiltonian conserved by NRPMD is given by

\[H_N = \sum_{\alpha=1}^N \left[ \frac{P_\alpha^2}{2M} + V_0(R_\alpha) + \frac{M}{2\beta_N^2\hbar^2} (R_\alpha - R_{\alpha-1})^2 + \frac{1}{\hbar} \sum_{nm}V_{nm}(R_\alpha) ([q_\alpha]_n[q_\alpha]_m + [p_\alpha]_n[p_\alpha]_m - \delta_{nm}\hbar) \right]\]

which contains $N$ replicas with positions $R_\alpha$ and momenta $P_\alpha$ joined by harmonic springs. $M$ is the mass, $\beta_N = \beta / N$ is the inverse temperature scaled by the number of beads. $V_0(R_\alpha)$ is the state independent potential. $V_{nm}$ are the matrix elements of the diabatic potential. The sum runs over all pairs of states. Each replica has a set of mapping variables $[q_\alpha]_n$ and $[p_\alpha]_n$ that interact only within the set associated with a single replica. The consequence of this is that the electronic dynamics is not contaminated by interbead coupling.

The equations of motion obtained from this Hamiltonian are

\[\begin{aligned} \dot{R}_\alpha &= \frac{P_\alpha}{M} \\ \dot{P}_\alpha &= - \frac{M}{\beta_N^2 \hbar^2}(2 R_\alpha - R_{\alpha+1} - R_{\alpha-1}) - \nabla_{R_\alpha} V_0(R_\alpha) - \frac{1}{2\hbar}\sum_{nm} \nabla_{R_\alpha} V_{nm}(R_\alpha) ([q_\alpha]_n[q_\alpha]_m + [p_\alpha]_n[p_\alpha]_m - \delta_{nm}\hbar) \\ [\dot{q}_\alpha]_n &= \frac{1}{\hbar} \sum_m V_{nm}(R_\alpha)[p_\alpha]_m \\ [\dot{p}_\alpha]_n &= -\frac{1}{\hbar} \sum_m V_{nm}(R_\alpha)[q_\alpha]_m \end{aligned}\]

Implementation details

Solving the differential equations

For mapping variable methods of this type, a symplectic algorithm ([13]) exists. This algorithm has the advantage of long time stability and can be easily combined with the standard algorithms for ring polymer time-evolution. For NRPMD we have implemented this algorithm using the Cayley modified ring polymer propagator ([14]) and obtain accurate and efficient dynamics. For few beads, similar performance to the OrdinaryDiffEq.jl algorithms is obtained, but as the number of beads increases this algorithm becomes more effective.

Generating the initial distribution

Currently, we provide this functionality only for nonequilibrium simulations where the nuclear part of the distribution is separable from the electronic part. Typically, the nuclear distribution will be sampled using Langevin dynamics or Monte Carlo sampling and the electronic variables are confined to a single electronic state. This is appropriate for modelling photoexcitation dynamics but is not yet suitable for equilibrium simulations. Equilibrium dynamics would require also sampling a thermal distribution for the mapping variables.

Form of the Hamiltonian

The diabatic models defined in NQCModels.jl are of the appropriate form for this method though they provide the potential as a single matrix, rather than separating the state-dependent and independent parts. It has been suggested that defining the Hamiltonian such that the lowest eigenvalue of the diabatic matrix is zero everywhere leads to improved convergence in the sampling ([9]). However, here we have not done this for simplicity when defining the models.

Example

Using NRPMD we can reproduce the Fig. 3a in the 2019 paper of Chowdhury and Huo ([11]).

First we generate a thermal ring polymer distribution in a harmonic potential. A simple way to do this is to use Monte Carlo sampling for the positions and obtain velocities from a Boltzmann distribution.

using NQCDynamics

atom = Atoms(1)

sim = RingPolymerSimulation(atom, Harmonic(dofs=1), 4; temperature=1/16)

r0 = zeros(size(sim))
steps = 5e3 # Number of Monte Carlo steps
step_size = Dict(:X=>1.0) # Monte Carlo step size for species :X
output = InitialConditions.ThermalMonteCarlo.run_advancedmh_sampling(sim, r0, steps, step_size)
velocities = VelocityBoltzmann(1/16, masses(sim), (1,1))

distribution = DynamicalDistribution(velocities, output, size(sim)) * PureState(1)

NQCDistributions.ProductDistribution{NQCDistributions.RingPolymerWrapper{NQCDistributions.UnivariateArray{2, Distributions.Normal{Float64}}}, NQCDistributions.ConfigurationVector{Vector{Array{Float64, 3}}}, PureState{Diabatic}}(DynamicalDistribution{NQCDistributions.RingPolymerWrapper{NQCDistributions.UnivariateArray{2, Distributions.Normal{Float64}}}, NQCDistributions.ConfigurationVector{Vector{Array{Float64, 3}}}}(NQCDistributions.RingPolymerWrapper{NQCDistributions.UnivariateArray{2, Distributions.Normal{Float64}}}(NQCDistributions.UnivariateArray{2, Distributions.Normal{Float64}}(Distributions.Normal{Float64}[Distributions.Normal{Float64}(μ=0.0, σ=0.25);;]), (1, 1, 4), Int64[]), NQCDistributions.ConfigurationVector{Vector{Array{Float64, 3}}}([[0.0;;; 0.0;;; 0.0;;; 0.0], [0.0;;; 0.0;;; 0.0;;; 0.0], [0.0;;; 0.0;;; 0.0;;; 0.0], [-0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983], [-0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983], [-0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983], [-0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983;;; -0.21407429209636983], [-0.0053779611686197715;;; -0.0053779611686197715;;; -0.0053779611686197715;;; -0.0053779611686197715], [-0.0053779611686197715;;; -0.0053779611686197715;;; -0.0053779611686197715;;; -0.0053779611686197715], [0.3465254531854926;;; 0.3465254531854926;;; 0.3465254531854926;;; 0.3465254531854926]  …  [0.42993766887967627;;; 0.39484919016300746;;; 0.38898097517046926;;; -0.37526517458729336], [0.42993766887967627;;; 0.39484919016300746;;; 0.38898097517046926;;; -0.37526517458729336], [0.42993766887967627;;; 0.39484919016300746;;; 0.38898097517046926;;; -0.37526517458729336], [0.242577622927569;;; 0.38115033910887464;;; 0.2016209292183619;;; -0.7362864154373752], [0.242577622927569;;; 0.38115033910887464;;; 0.2016209292183619;;; -0.7362864154373752], [0.242577622927569;;; 0.38115033910887464;;; 0.2016209292183619;;; -0.7362864154373752], [-0.06326507943914242;;; 0.07530763674216323;;; -0.10422177314834952;;; -1.0421291178040866], [-0.06326507943914242;;; 0.07530763674216323;;; -0.10422177314834952;;; -1.0421291178040866], [-0.06326507943914242;;; 0.07530763674216323;;; -0.10422177314834952;;; -1.0421291178040866], [-0.06326507943914242;;; 0.07530763674216323;;; -0.10422177314834952;;; -1.0421291178040866]]), Random.Xoshiro(0x833d3bea9340db16, 0x8405c951be83778e, 0xf8450b3b03c6d12f, 0x585da5d8cb5793bd, 0x533758f4711261d6)), PureState{Diabatic}(1, Diabatic()))

We can check the distribution by plotting the phasespace diagram for each of the points in our distribution:

using Plots

nuclear = distribution.nuclear
flat_position = reduce(vcat, (nuclear.position[i][:] for i in 1:length(nuclear)))
flat_velocity = reduce(vcat, (rand(nuclear.velocity)[:] for _ in 1:length(nuclear)))
scatter(flat_position, flat_velocity)

The simulation method is given as the type parameter {NRPMD} and the simulation constructor is given the atoms, model, number of beads, temperature and degrees of freedom.

sim = RingPolymerSimulation{NRPMD}(atom, DoubleWell(γ=0.1), 4; temperature=1/16)

RingPolymerSimulation{NRPMD{Float64}}:

  Atoms{Float64}([:X], [0], [1.0])

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

  with 4 beads.

Next, we can use this distribution as a starting point for the dynamics simulations. This will result in each trajectory starting from a random configuration in the distribution. For NRPMD, the electronic variables are sampled from a Gaussian, independent of the initial electronic state. The electronic state is introduced in the correlation function expression when correlating the initial and final populations.

The quantities output by the ensemble simulation are specified by the output and the reduction. The output follows the DifferentialEquations format where we provide a function that determines the output of each trajectory. The reduction can be one of :mean, :append, or :sum, which will determine how the data from each trajectory is combined.

output = TimeCorrelationFunctions.PopulationCorrelationFunction(sim, Diabatic())

Finally, we can combine the parameters and run the simulation. The resulting plot shows the time dependent population difference and closely matches the figure from the paper we were attempting to reproduce. Nice!

ensemble = run_dynamics(sim, (0.0, 30.0), distribution;
    trajectories=100, output, reduction=MeanReduction(), dt=0.1)

plot(0:0.1:30, [p[1,1]-p[2,1] for p in ensemble[:PopulationCorrelationFunction]])
xlabel!("Time")
ylabel!("Population difference")