Thermal Hamiltonian Monte Carlo

Our implementation of Hamiltonian Monte Carlo (HMC) is a light wrapper around the AdvancedHMC.jl package. If you want to learn about the HMC theory, refer to the references and documentation provided with AdvancedHMC.jl.

Currently, our implementation works for systems with classical nuclei only (i.e. Simulation but not RingPolymerSimulation).

Example

In this example we use Hamiltonian Monte Carlo to sample the canonical distribution of a 3 dimensional harmonic oscillator potential containing 4 atoms.

using NQCDynamics
using Unitful
using UnitfulAtomic

sim = Simulation(Atoms([:H, :H, :C, :C]), Harmonic(dofs=3); temperature=300u"K")
r0 = randn(size(sim))
chain, stats = InitialConditions.ThermalMonteCarlo.run_advancedhmc_sampling(sim, r0, 1e4)

┌ Info: Finished 1000 adapation steps
│   adaptor =
│    StanHMCAdaptor(
│        pc=WelfordVar,
│        ssa=NesterovDualAveraging(γ=0.05, t_0=10.0, κ=0.75, δ=0.5, state.ϵ=1.2047514421543333),
│        init_buffer=75, term_buffer=50, window_size=25,
│        state=window(76, 950), window_splits(100, 150, 250, 450, 950)
│    )
│   κ.τ.integrator = Leapfrog(ϵ=1.2)
└   h.metric = DiagEuclideanMetric([0.0009206900730121456, 0.0 ...])
┌ Info: Finished 10000 sampling steps for 1 chains in 1.055385028 (s)
│   h = Hamiltonian(metric=DiagEuclideanMetric([0.0009206900730121456, 0.0 ...]), kinetic=AdvancedHMC.GaussianKinetic())
│   κ = AdvancedHMC.HMCKernel{AdvancedHMC.FullMomentumRefreshment, AdvancedHMC.Trajectory{AdvancedHMC.MultinomialTS, AdvancedHMC.Leapfrog{Float64}, AdvancedHMC.GeneralisedNoUTurn{Float64}}}(AdvancedHMC.FullMomentumRefreshment(), Trajectory{AdvancedHMC.MultinomialTS}(integrator=Leapfrog(ϵ=1.2), tc=AdvancedHMC.GeneralisedNoUTurn{Float64}(10, 1000.0)))
│   EBFMI_est = 0.5009181913667998
└   average_acceptance_rate = 0.5449400458233716

The Monte Carlo chain contains the nuclear configurations that we have sampled:

chain

10000-element Vector{Matrix{Float64}}:
 [-0.18702369151065656 -0.1429664086847366 -0.1362092782389347 -0.7443053669778412; -0.28236355809407754 0.04621241489844985 0.09495473110531549 -0.019242249857139127; -0.002262898366013011 -0.028272418984512993 -0.08293215911708007 -0.12095048228444102]
 [-0.18702369151065656 -0.1429664086847366 -0.1362092782389347 -0.7443053669778412; -0.28236355809407754 0.04621241489844985 0.09495473110531549 -0.019242249857139127; -0.002262898366013011 -0.028272418984512993 -0.08293215911708007 -0.12095048228444102]
 [-0.18702369151065656 -0.1429664086847366 -0.1362092782389347 -0.7443053669778412; -0.28236355809407754 0.04621241489844985 0.09495473110531549 -0.019242249857139127; -0.002262898366013011 -0.028272418984512993 -0.08293215911708007 -0.12095048228444102]
 [-0.18702369151065656 -0.1429664086847366 -0.1362092782389347 -0.7443053669778412; -0.28236355809407754 0.04621241489844985 0.09495473110531549 -0.019242249857139127; -0.002262898366013011 -0.028272418984512993 -0.08293215911708007 -0.12095048228444102]
 [0.05411772237761993 0.07761846990189805 0.04725199421439377 0.282698577452358; 0.10602824123910376 -0.05024911671492067 0.014168215044135263 0.029087114457002627; 0.015761316771472945 -0.05608814462814801 0.021401815636497082 0.07749758928214279]
 [0.05411772237761993 0.07761846990189805 0.04725199421439377 0.282698577452358; 0.10602824123910376 -0.05024911671492067 0.014168215044135263 0.029087114457002627; 0.015761316771472945 -0.05608814462814801 0.021401815636497082 0.07749758928214279]
 [-0.010847706907335348 0.0031361445906835095 0.03248546228578895 0.0349869805049067; 0.041291661409957936 -0.03816742080531686 -0.05806999884707099 0.042427181915870534; 0.0482793330560119 0.04203858886727082 -0.04819451984959677 -0.0028352797007443406]
 [-0.010847706907335348 0.0031361445906835095 0.03248546228578895 0.0349869805049067; 0.041291661409957936 -0.03816742080531686 -0.05806999884707099 0.042427181915870534; 0.0482793330560119 0.04203858886727082 -0.04819451984959677 -0.0028352797007443406]
 [-0.03318207912950866 -0.01417239672031793 0.05849562296315093 -0.021356006589175778; 0.029578329866626965 -0.06441472599325589 -0.017358858777341057 -0.08240129473859376; -0.03112167738017768 -0.01196706747843436 0.013556026885373523 0.04226188572988374]
 [-0.03318207912950866 -0.01417239672031793 0.05849562296315093 -0.021356006589175778; 0.029578329866626965 -0.06441472599325589 -0.017358858777341057 -0.08240129473859376; -0.03112167738017768 -0.01196706747843436 0.013556026885373523 0.04226188572988374]
 ⋮
 [0.03731739143286175 -0.01979398303706387 0.03242861300931135 -0.03842771617092325; 0.03787939767476857 0.08199414511944919 0.03477142042299815 -0.04638042884707205; -0.048374806170426984 0.01070278078536862 0.012357774794438808 -0.010002510938107619]
 [-0.0053862733812210725 0.04023610797709845 0.03622087612683379 0.0588770889995055; 0.035858986454445754 0.03427719303151072 0.01867901868603431 -0.020257157247901565; -0.0008203090023072152 0.015218581843829545 0.007758316316040532 0.01731457373666642]
 [0.0363560914139223 -0.00984942535382146 0.00569033790003054 -0.018161613880326874; 0.0248633374429187 -0.0060975753572064945 0.029351513887321923 -0.04804011363259912; 0.013768419516317878 -0.021928255904137328 -0.03402619475583456 -0.0038771315858318195]
 [-0.0258471901354908 0.029902848769395043 -0.02245383675717488 0.00862939614571068; -0.007660914625753218 0.009943228166429941 -0.018627733075385476 0.05777608004877373; 0.01827789196201937 -0.017210869088939408 0.028005826769716465 -0.020804944585612148]
 [0.046669588000757885 -0.030763903080171667 -0.00250005067740576 -0.017501322663089795; -0.003537954067175754 0.018737602872320826 -0.011522827932484526 -0.04444370268351591; -0.018105672174349743 -0.016056375878169903 -0.029156459610815387 -0.02989761071108038]
 [0.046669588000757885 -0.030763903080171667 -0.00250005067740576 -0.017501322663089795; -0.003537954067175754 0.018737602872320826 -0.011522827932484526 -0.04444370268351591; -0.018105672174349743 -0.016056375878169903 -0.029156459610815387 -0.02989761071108038]
 [-0.012432116664398689 0.051429973079082994 0.0016178699301378089 -0.0008155541399164079; -0.009780880055441839 0.022605456803000072 0.018831024714077528 -0.03488243778189068; -0.05133986561366238 0.017797535757866592 0.040691133717691626 0.009851093917678563]
 [-0.007704314526807598 -0.03297974576307204 0.008138547376680544 -0.03025915837964807; 0.001205224787018705 -0.03593106158304516 -0.01834586362249781 0.03095304141792815; 0.05503404853494512 -0.005143470667760986 0.002154813561756022 -0.014275976103803486]
 [0.011030193478164368 0.02024675179806119 0.016247179360897555 -0.03498667206646926; 0.004928048034646061 0.029918554047772983 0.004046415301190612 -0.04720163184911465; -0.03274082844662785 -0.022322078301015486 0.001905698891515881 -0.012239400865161794]

and stats contains extra information about the sampling procedure:

stats

10000-element Vector{NamedTuple}:
 (n_steps = 5, is_accept = true, acceptance_rate = 1.0, log_density = -390.25758825412584, hamiltonian_energy = 1501.2299630427765, hamiltonian_energy_error = -2121.1688089600007, max_hamiltonian_energy_error = -2180.5127923849573, tree_depth = 2, numerical_error = false, step_size = 0.05, nom_step_size = 0.05, is_adapt = true)
 (n_steps = 1, is_accept = true, acceptance_rate = 0.0, log_density = -390.25758825412584, hamiltonian_energy = 393.1070026384194, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 1.0404938890405309e11, tree_depth = 0, numerical_error = true, step_size = 1.241032542311506, nom_step_size = 1.241032542311506, is_adapt = true)
 (n_steps = 1, is_accept = true, acceptance_rate = 0.0, log_density = -390.25758825412584, hamiltonian_energy = 397.24483736171396, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 4.364525427314434e8, tree_depth = 0, numerical_error = true, step_size = 0.5, nom_step_size = 0.5, is_adapt = true)
 (n_steps = 1, is_accept = true, acceptance_rate = 0.0, log_density = -390.25758825412584, hamiltonian_energy = 392.8832460386449, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 122797.08970005889, tree_depth = 0, numerical_error = true, step_size = 0.13192866018819532, nom_step_size = 0.13192866018819532, is_adapt = true)
 (n_steps = 7, is_accept = true, acceptance_rate = 1.0, log_density = -60.93224600609105, hamiltonian_energy = 323.9591409443618, hamiltonian_energy_error = -71.46262068455957, max_hamiltonian_energy_error = -73.5520725878269, tree_depth = 2, numerical_error = false, step_size = 0.028716309633808685, nom_step_size = 0.028716309633808685, is_adapt = true)
 (n_steps = 1, is_accept = true, acceptance_rate = 0.0, log_density = -60.93224600609105, hamiltonian_energy = 67.3765712542574, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 5967.825821060741, tree_depth = 0, numerical_error = true, step_size = 0.11260612534950616, nom_step_size = 0.11260612534950616, is_adapt = true)
 (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -9.036288404483178, hamiltonian_energy = 61.01945215300068, hamiltonian_energy_error = -7.477740224098589, max_hamiltonian_energy_error = -7.477740224098589, tree_depth = 2, numerical_error = false, step_size = 0.02340023160450795, nom_step_size = 0.02340023160450795, is_adapt = true)
 (n_steps = 1, is_accept = true, acceptance_rate = 0.0, log_density = -9.036288404483178, hamiltonian_energy = 18.400044421712792, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 1092.9572281306102, tree_depth = 0, numerical_error = true, step_size = 0.10545494475829686, nom_step_size = 0.10545494475829686, is_adapt = true)
 (n_steps = 7, is_accept = true, acceptance_rate = 0.8477579905927107, log_density = -10.724122118676133, hamiltonian_energy = 19.012676570218908, hamiltonian_energy_error = 0.2068971157834838, max_hamiltonian_energy_error = 0.33061403928448385, tree_depth = 3, numerical_error = false, step_size = 0.021583114937805906, nom_step_size = 0.021583114937805906, is_adapt = true)
 (n_steps = 1, is_accept = true, acceptance_rate = 1.3405651372572878e-12, log_density = -10.724122118676133, hamiltonian_energy = 14.852974123619404, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 27.33792984665297, tree_depth = 1, numerical_error = false, step_size = 0.06374536590866547, nom_step_size = 0.06374536590866547, is_adapt = true)
 ⋮
 (n_steps = 3, is_accept = true, acceptance_rate = 0.2825045035606695, log_density = -9.756508547108144, hamiltonian_energy = 20.598945478632473, hamiltonian_energy_error = 0.4441730882679309, max_hamiltonian_energy_error = 2.5608963633596957, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -5.388561283090916, hamiltonian_energy = 11.591153170000984, hamiltonian_energy_error = -1.746824353339063, max_hamiltonian_energy_error = -1.746824353339063, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.8046577028038145, log_density = -3.9203271696881927, hamiltonian_energy = 8.00560414584591, hamiltonian_energy_error = -0.4999394338932195, max_hamiltonian_energy_error = 0.5121813622919866, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.43155980510175546, log_density = -4.121385908580034, hamiltonian_energy = 9.385840413320771, hamiltonian_energy_error = 0.0859474598379304, max_hamiltonian_energy_error = 1.89261220217284, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.44330161679399765, log_density = -4.335715670082412, hamiltonian_energy = 10.553570720391104, hamiltonian_energy_error = 0.007197779946954697, max_hamiltonian_energy_error = 2.2810592208929883, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.07811043142420909, log_density = -4.335715670082412, hamiltonian_energy = 16.829497016083845, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 5.338022214810152, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.46390194914255845, log_density = -5.097815703845807, hamiltonian_energy = 9.005633195753806, hamiltonian_energy_error = 0.3041664277484859, max_hamiltonian_energy_error = 1.2569026757983153, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.7918176893110767, log_density = -4.199637550531728, hamiltonian_energy = 8.468198411888977, hamiltonian_energy_error = -0.293665507479405, max_hamiltonian_energy_error = 0.4219595872106261, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.3674349054192214, log_density = -3.6351383487962248, hamiltonian_energy = 12.303569198695826, hamiltonian_energy_error = -0.1944807768311474, max_hamiltonian_energy_error = 3.777590005667383, tree_depth = 2, numerical_error = false, step_size = 1.2047514421543333, nom_step_size = 1.2047514421543333, is_adapt = false)

Here we should see that the energy expectation for the generated ensemble matches with the equipartition theorem:

julia> Estimators.@estimate potential_energy(sim, chain)0.0058615208982014475
julia> austrip(sim.temperature) * 3 * 4 / 20.005700260814198944