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(0x8c73ab8a269bf427, 0x6a1e8ec61e0f8ca5, 0x2e962b97b372d542, 0x23b4ff8950ab1328, 0xadb12e9436b64a78), 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.248670699304881
 9.238455928832426
 9.24411508416325
 9.238443181945959
 9.25376445502399
 9.251970971138366
 9.24385013456313
 9.245434988656488
 9.246848471373104
 9.254286975941088

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