Ehrenfest molecular dynamics

The Ehrenfest method is a mixed quantum-classical dynamics method in which the total wavefunction is factorized into slow (nuclear) variables, which are treated classically, and fast ones (electrons) which remain quantum-mechanical. In the Ehrenfest method, nuclei move according to classical mechanics on a potential energy surface given by the expectation value of the electronic Hamiltonian.

The time dependence of the electronic wavefunction is expanded into an adiabatic basis and follows the time-dependent Schr\"odinger equation.

\[i\hbar \dot{c}_i(t) = V_i(\mathbf{R}) c_i (t) - i\hbar \sum_j \dot{\mathbf{R}} \cdot \mathbf{d}_{ij}(\mathbf{R})c_j(t)\]

Example

Below the example of the Ehrenfest implementation is presented, using model from [12].

At the start, we assign atoms and initialise the simulation using the mass and model from NQCModels.jl.

using NQCDynamics

atoms = Atoms(1980)
sim = Simulation{Ehrenfest}(atoms, AnanthModelOne())

Simulation{Ehrenfest{Float64}}:
  Atoms{Float64}([:X], [0], [1980.0])
  AnanthModelOne{Float64, Float64, Float64, Float64}
  a: Float64 0.01
  b: Float64 1.6
  c: Float64 0.005
  d: Float64 1.0

Next, the initial distribution is created:

using Distributions
e = 0.03
k = sqrt(e*2*atoms.masses[1])
r = Normal(-5, 1/sqrt(0.25))
v = k / atoms.masses[1]
distribution = DynamicalDistribution(v, r, size(sim))* PureState(1, Adiabatic())

NQCDistributions.ProductDistribution{NQCDistributions.FixedFill{Float64}, NQCDistributions.UnivariateFill{Distributions.Normal{Float64}}, PureState{Adiabatic}}(DynamicalDistribution{NQCDistributions.FixedFill{Float64}, NQCDistributions.UnivariateFill{Distributions.Normal{Float64}}}(NQCDistributions.FixedFill{Float64}(0.005504818825631803, (1, 1)), NQCDistributions.UnivariateFill{Distributions.Normal{Float64}}(Distributions.Normal{Float64}(μ=-5.0, σ=2.0), (1, 1)), Random.Xoshiro(0x904f9892cd710e68, 0xef8589e74732188a, 0x0aaf50294f00f64e, 0xa424eb6193512b8c, 0x114f36c235445860), Int64[]), PureState{Adiabatic}(1, Adiabatic()))

To run an ensemble simulation we additionally choose number of trajectories n_traj and timespan tspan and we pass all the established settings to the run_dynamics function. In this example we output velocities by specifying output=OutputVelocity and store the final values in the final_velocities array. Following that, we calculate final momenta.

n_traj = 10
tspan = (0.0, 3000.0)
solution = run_dynamics(sim, tspan, distribution;
    trajectories=n_traj, output=OutputVelocity, dt=1.0)
final_velocities = [r[:OutputVelocity][end] for r in solution]
momenta = reduce(vcat, final_velocities*atoms.masses[1])

10×1 Matrix{Float64}:
 9.67388264375856
 9.678171850484937
 9.677024614663175
 9.67577333469034
 9.674967728499983
 9.675082632477238
 9.67644021809818
 9.695798931621782
 9.675538768173183
 9.676468130294541

using Plots
histogram(momenta)
xlims!(-20,20)