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 [5].

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(0x48fe30df6de1a29e, 0x3a091825261ec9aa, 0x2ae04641a2573fb0, 0xa980471ecddf69e9, 0x6ffcca6ead0e236b)), 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.780254804052891
 9.780254793490254
 9.780254817855637
 9.780282580684883
 9.780255800342381
 9.780255533654351
 9.780254794623353
 9.780254792553725
 9.780254794417026
 9.780254807883898

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