Ring polymer surface hopping (RPSH)

Ring polymer surface hopping was one of the early attempts to extend RPMD to the realm of nonadiabatic dynamics [17]. On the surface, the concept is reasonably simple. Since RPMD proceeds on a single adiabatic surface, it should be possible to directly combine the FSSH scheme with ring polymer dynamics to approximately include nuclear quantum effects in the surface hopping dynamics. However, there are some ambiguities surrounding the exact implementation when considering how to couple the electronic equations to the nuclear equations and how the velocity rescaling should be implemented.

Originally, two varieties were proposed: the bead approximation and the centroid approximation. The centroid approximation is the simpler of the two and involves directly replacing the classical particle in the FSSH algorithm with the ring polymer centroid. This means that the nonadiabatic couplings evaluated at the centroid and the centroid velocity are used to propagate the electronic equations, and the kinetic energy is conserved on the centroid level. This is the version that is implemented here.

The bead approximation involves evaluating the nonadiabatic couplings for every bead and using these contributions from every bead to propagate the electronics. This version acts to conserve the kinetic energy for the entire ring polymer. [17] describes both the centroid and bead approximations, [18] uses the centroid approximation.

Example

In this example we can apply RPSH to the ThreeStateMorse model as shown in the supporting info of [18]. This model has a single particle with mass of 20000 a.u. and we use 4 beads for the ring polymer.

using NQCDynamics
using Unitful

atoms = Atoms(20000)
model = ThreeStateMorse()
sim = RingPolymerSimulation{FSSH}(atoms, model, 4; temperature=300u"K")

For our initial conditions let's use a position distribution centred at 2.1 a.u. with Boltzmann velocities at 300 K.

using Distributions: Normal

position = Normal(2.1, 1 / sqrt(20000 * 0.005))
velocity = VelocityBoltzmann(300u"K" * nbeads(sim), masses(sim), (1,1))
distribution = DynamicalDistribution(velocity, position, size(sim)) * PureState(1)

Now let's run an ensemble of trajectories that sample from this distribution. For the output we will receive the diabatic population at intervals of t=50 and it will be averaged over all trajectories by the :mean keyword.

solution = run_dynamics(sim, (0.0, 3000.0), distribution;
    saveat=50, trajectories=50, dt=1,
    output=TimeCorrelationFunctions.PopulationCorrelationFunction(sim, Diabatic()),
    reduction=MeanReduction())

2-element Dictionaries.Dictionary{Symbol, Vector}
                          :Time │ [0.0, 50.0, 100.0, 150.0, 200.0, 250.0, 300.0…
 :PopulationCorrelationFunction │ [[1.0 -6.446665200863662e-45 -3.6946187132410…

Here we plot diabatic population of each state as a function of time. The result closely resembles the literature reference ([18]). The small discrepancy that occurs at around t=2000 is due to our use of an alternative method to calculate the diabatic populations. A discussion on this topic is available from [19].

using Plots

plot(0:50:3000, [p[1,1] for p in solution[:PopulationCorrelationFunction]], label="State 1")
plot!(0:50:3000, [p[1,2] for p in solution[:PopulationCorrelationFunction]], label="State 2")
plot!(0:50:3000, [p[1,3] for p in solution[:PopulationCorrelationFunction]], label="State 3")
xlabel!("Time /a.u.")
ylabel!("Population")

For our simulation we are using a Normal distribution to initialise our ring polymer configuration. Since ring polymer surface hopping has not been rigorously derived, this choice is somewhat arbitrary and it is possible that better results could be achieved using a free ring polymer distribution instead. [20] provides a theoretical description of how nonequilibrium simulations using RPMD should be performed. This techniques here should likely be applied to RPSH too.