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.1716243001640527),
│        init_buffer=75, term_buffer=50, window_size=25,
│        state=window(76, 950), window_splits(100, 150, 250, 450, 950)
│    )
│   κ.τ.integrator = Leapfrog(ϵ=1.17)
└   h.metric = DiagEuclideanMetric([0.000971159518946784, 0.00 ...])
┌ Info: Finished 10000 sampling steps for 1 chains in 1.196925815 (s)
│   h = Hamiltonian(metric=DiagEuclideanMetric([0.000971159518946784, 0.00 ...]), 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.17), tc=AdvancedHMC.GeneralisedNoUTurn{Float64}(10, 1000.0)))
│   EBFMI_est = 0.4556994867987458
└   average_acceptance_rate = 0.5721180849410592

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

chain

10000-element Vector{Matrix{Float64}}:
 [-0.12673634671146405 0.05885557729528944 0.1939908036153582 0.30211506496555385; -0.23079273774043607 -0.25668085644950556 -0.23173040019943725 0.2718333113325193; 0.7572152220802026 0.49909804739942953 0.4802643207938082 -0.16280807211258316]
 [-0.12673634671146405 0.05885557729528944 0.1939908036153582 0.30211506496555385; -0.23079273774043607 -0.25668085644950556 -0.23173040019943725 0.2718333113325193; 0.7572152220802026 0.49909804739942953 0.4802643207938082 -0.16280807211258316]
 [-0.12673634671146405 0.05885557729528944 0.1939908036153582 0.30211506496555385; -0.23079273774043607 -0.25668085644950556 -0.23173040019943725 0.2718333113325193; 0.7572152220802026 0.49909804739942953 0.4802643207938082 -0.16280807211258316]
 [-0.12673634671146405 0.05885557729528944 0.1939908036153582 0.30211506496555385; -0.23079273774043607 -0.25668085644950556 -0.23173040019943725 0.2718333113325193; 0.7572152220802026 0.49909804739942953 0.4802643207938082 -0.16280807211258316]
 [0.026942686060849944 -0.08823987329516632 -0.08060265600675856 -0.07282225649765892; 0.06969751134334226 0.0762873867707482 0.13925510000892846 -0.07508131343548405; -0.30318942836396945 -0.1862126313240074 -0.17756789773921833 0.053408357312542704]
 [0.026942686060849944 -0.08823987329516632 -0.08060265600675856 -0.07282225649765892; 0.06969751134334226 0.0762873867707482 0.13925510000892846 -0.07508131343548405; -0.30318942836396945 -0.1862126313240074 -0.17756789773921833 0.053408357312542704]
 [0.02339535846176699 -0.034796924259236535 -0.028358861212473482 -0.047549763548871984; -0.01855547162951464 0.02797287179223134 -0.007855169462258421 0.025854534548282383; -0.040152681773775484 0.030683169894240983 0.03676131272724288 0.0077380986742362]
 [0.02339535846176699 -0.034796924259236535 -0.028358861212473482 -0.047549763548871984; -0.01855547162951464 0.02797287179223134 -0.007855169462258421 0.025854534548282383; -0.040152681773775484 0.030683169894240983 0.03676131272724288 0.0077380986742362]
 [0.0165755527694333 -0.007368264154717304 -0.04296071321544532 -0.02540743418609091; -0.008071927915912346 0.04223350560651044 0.000765571084669573 0.01224968331533875; -0.011960117071897758 0.025197578742755244 0.061090156287509995 0.032574367921005]
 [0.0165755527694333 -0.007368264154717304 -0.04296071321544532 -0.02540743418609091; -0.008071927915912346 0.04223350560651044 0.000765571084669573 0.01224968331533875; -0.011960117071897758 0.025197578742755244 0.061090156287509995 0.032574367921005]
 ⋮
 [0.026592368908189997 0.011717336931192028 -0.012319608105480682 0.031464142325002956; 0.025884200214275503 0.010205116126045917 0.053870379073900215 -0.031976808234493705; -0.05642563588251998 0.021949354328078294 0.03406808292682186 -0.010512953442351348]
 [0.007205165985835718 -0.02143067799904936 -0.04704163004286905 -0.010912408956891294; -0.032226890310870816 -0.0004677276669845364 0.0002584991767331096 0.029470232519479778; 0.017990585320177686 -0.010739682769971944 -0.04066373247543281 0.025299042585685376]
 [0.009119546247320589 0.05006204864706457 -0.008992362510350384 0.005663358478823061; -0.009106615059138082 -0.008463001060229347 0.039531221907103384 0.0586423898411144; -0.012515724900399077 0.014081164976074426 0.010842166526094654 0.07587175285983216]
 [-0.06283630443532423 0.03370560633704679 0.011740365835803804 -0.06269523066672893; 0.017349111214293472 -0.023295868915768216 0.017849503972381206 0.015774262686575814; 0.020363611359301582 0.018497683164407273 -0.027469733631352564 0.011368444831670793]
 [0.045818550834164404 -0.04619243539006194 -0.013094979940475033 0.06566599223342126; -0.01768269159590705 0.00835135186100662 0.012970916388572132 0.013521052798452754; -0.005598619143901799 -0.027813789250128694 0.02606607882867374 0.013397150441315082]
 [0.04173607904442254 -0.016983795715604533 0.009796787370845795 -0.05793217622890634; 0.02995464524488115 0.0044423626767740665 0.039529743328171135 0.019974107472083426; -0.003762930310015182 -0.023061061246159205 0.0030899196983427794 0.07897911596344803]
 [-0.041215778241813515 0.03997790405995584 0.03959991105941456 -0.025902601794583688; 0.06650612350607747 0.044936123493064584 0.0017254441814970808 0.025176224036341666; 0.03183326488390651 0.00017614462182924395 -0.014115543723297333 0.012556776071490486]
 [-0.041215778241813515 0.03997790405995584 0.03959991105941456 -0.025902601794583688; 0.06650612350607747 0.044936123493064584 0.0017254441814970808 0.025176224036341666; 0.03183326488390651 0.00017614462182924395 -0.014115543723297333 0.012556776071490486]
 [-0.011729800273325278 0.010103249894489134 0.0621115646724079 0.03506994335105567; 0.019352876211531325 0.006242000309435926 0.03159266533535275 -0.001191328023692595; -0.037971760140191706 -0.03901093434477773 0.0005004104690609131 -0.02592915333811737]

and stats contains extra information about the sampling procedure:

stats

10000-element Vector{NamedTuple}:
 (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -776.1796117878215, hamiltonian_energy = 2851.654388871193, hamiltonian_energy_error = -3973.482213938826, max_hamiltonian_energy_error = -4069.6111309484127, tree_depth = 1, 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 = -776.1796117878215, hamiltonian_energy = 780.9867355872484, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 2.07179694607044e11, 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 = -776.1796117878215, hamiltonian_energy = 779.6866245045163, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 8.757894855079418e8, 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 = -776.1796117878215, hamiltonian_energy = 781.7367172594874, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 238527.81894299408, tree_depth = 0, numerical_error = true, step_size = 0.13192866018819532, nom_step_size = 0.13192866018819532, is_adapt = true)
 (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -114.20558073884924, hamiltonian_energy = 637.9202457567856, hamiltonian_energy_error = -143.64639769584528, max_hamiltonian_energy_error = -143.64639769584528, 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 = -114.20558073884924, hamiltonian_energy = 119.9468871855698, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 12093.77132004534, 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 = -5.602527087369944, hamiltonian_energy = 103.27887628972464, hamiltonian_energy_error = -15.648722179556515, max_hamiltonian_energy_error = -15.648722179556515, 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 = 5.811854917211381e-287, log_density = -5.602527087369944, hamiltonian_energy = 10.438469548863415, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 659.0820219065217, tree_depth = 1, numerical_error = false, step_size = 0.10545494475829686, nom_step_size = 0.10545494475829686, is_adapt = true)
 (n_steps = 7, is_accept = true, acceptance_rate = 1.0, log_density = -5.4685527920125345, hamiltonian_energy = 9.534742771905078, hamiltonian_energy_error = -0.016422764319424132, max_hamiltonian_energy_error = -0.1441230322093361, 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 = 5.412702454668533e-268, log_density = -5.4685527920125345, hamiltonian_energy = 10.391568981365335, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 615.4040564247556, tree_depth = 1, numerical_error = false, step_size = 0.1030960141257045, nom_step_size = 0.1030960141257045, is_adapt = true)
 ⋮
 (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -6.116387655252432, hamiltonian_energy = 9.7933765215547, hamiltonian_energy_error = -0.1748326685688255, max_hamiltonian_energy_error = -0.499435505247126, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.5702870760224644, log_density = -3.9383036367454203, hamiltonian_energy = 12.219811086813925, hamiltonian_energy_error = -0.7001957757386528, max_hamiltonian_energy_error = 2.092480183653141, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.4012389657469034, log_density = -7.414139331965236, hamiltonian_energy = 10.341091677963956, hamiltonian_energy_error = 1.3329859892207043, max_hamiltonian_energy_error = 1.4071247496017616, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.9654005750202046, log_density = -6.423278337944673, hamiltonian_energy = 11.304394759288325, hamiltonian_energy_error = -0.3792746184110545, max_hamiltonian_energy_error = -0.3792746184110545, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 1.0, log_density = -5.849176274901759, hamiltonian_energy = 9.811923316304139, hamiltonian_energy_error = -0.18991478177656518, max_hamiltonian_energy_error = -0.43925703590216614, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.443689539695955, log_density = -7.975564107063307, hamiltonian_energy = 12.027468180803707, hamiltonian_energy_error = 0.8441526331130973, max_hamiltonian_energy_error = 1.2090182971460717, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.877663622118931, log_density = -7.3604626700566955, hamiltonian_energy = 12.775353966349666, hamiltonian_energy_error = -0.4873985887312493, max_hamiltonian_energy_error = -0.4873985887312493, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.11759945228502912, log_density = -7.3604626700566955, hamiltonian_energy = 18.532406610464925, hamiltonian_energy_error = 0.0, max_hamiltonian_energy_error = 3.6238839408870334, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, is_adapt = false)
 (n_steps = 3, is_accept = true, acceptance_rate = 0.9503837996190093, log_density = -5.461172422554094, hamiltonian_energy = 10.699720321452574, hamiltonian_energy_error = -0.632224469221212, max_hamiltonian_energy_error = -0.632224469221212, tree_depth = 2, numerical_error = false, step_size = 1.1716243001640527, nom_step_size = 1.1716243001640527, 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.00598271217822692
julia> austrip(sim.temperature) * 3 * 4 / 20.005700260814198944