Thermal Metropolis-Hastings Monte Carlo
Metropolis-Hastings Monte Carlo is a popular method for sampling the canonical distribution for a molecular system. Our implementations uses AdvancedMH.jl from the Turing organisation.
For a classical Simulation, the algorithm involves proposing new configurations in a random walk starting from an initial configuration. These are accepted or rejected based upon the Metropolis-Hastings criteria. The result is a Markov chain that samples the canonical distribution.
!!! Legacy version
Prior to the use of [`AdvancedMH.jl`](https://github.com/TuringLang/AdvancedMH.jl),
an alternative version of the algorithm was implemented that works for both classical
and ring polymer systems: `MetropolisHastings.run_monte_carlo_sampling(sim, R0, Δ, passes)`
This is currently still included in the code but should be regarded as deprecated and
will likely be removed/combined with the [`AdvancedMH.jl`](https://github.com/TuringLang/AdvancedMH.jl)
version.Example 1
We can perform the sampling by setting up a classical simulation in the usual way and providing an appropriate initial configuration.
using NQCDynamics
sim = Simulation(Atoms([:H, :H, :H, :H, :H]), Harmonic(); temperature=15)
r0 = zeros(size(sim))1×5 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0Then we must also specify the total number of steps and the size of each step. These can be provided in a dictionary for each species to allow for different step sizes depending on the element in the simulation.
steps = 1e4
step_size = Dict(:H=>1)Dict{Symbol, Int64} with 1 entry:
:H => 1Now we can run the sampling. The extra keyword argument movement_ratio is used to specify the fraction of the system moved during each Monte Carlo step. If we attempt to move the entire system at once, we can expect a very low acceptance ratio, whereas is we move only a single atom, the sampling will take much longer. You will likely have to experiment with this parameter to achieve optimal sampling.
Keyword arguments relating to how much of the system to sample were recently changed. The existing move_ratio and internal_ratio arguments are no longer used.
Instead, there are now two options to specify how much of the system to move:
movement_ratio: Defines which fraction of the system to move. 1 moves the entire system.
stop_ratio: Defines which fraction of the system not to move. 1 stops the entire system.
movement_ratio_internal: Defines which proportion of ring polymer normal modes to perturb. 1 moves the entire system.
stop_ratio_internal: Defines which proportion of ring polymer normal modes not to perturb. 1 stops the entire system.
using NQCDynamics.InitialConditions: ThermalMonteCarlo
chain = ThermalMonteCarlo.run_advancedmh_sampling(sim, r0, steps, step_size; movement_ratio=0.5)10000-element Vector{Matrix{Float64}}:
[0.0 0.0 … 0.0 0.0]
[0.0 0.0 … 0.0 1.3648077044021925]
[0.5377529809081826 0.0 … 0.0 1.3648077044021925]
[0.5377529809081826 0.0 … 0.0 1.9337298373855603]
[-0.050866865001635175 0.0 … -0.21727795487781842 1.9337298373855603]
[2.0469557865132892 0.0 … -0.21727795487781842 1.9337298373855603]
[2.0469557865132892 0.0 … 0.30340200895367164 1.8595419658728367]
[-1.1985738316051053 0.0 … -0.19302334108688213 -0.7141587012626569]
[-1.1419722132797596 -1.8169326384132611 … -0.19302334108688213 -0.7141587012626569]
[-1.1419722132797596 -1.8169326384132611 … -0.19302334108688213 -0.7141587012626569]
⋮
[0.48250056255314533 3.2843153014442756 … -1.1471093372723562 -2.079076584655851]
[0.48250056255314533 3.2843153014442756 … -1.1471093372723562 -2.079076584655851]
[0.48250056255314533 2.002405061699317 … -1.056153225252739 -2.8817200151833804]
[0.48250056255314533 2.002405061699317 … -1.056153225252739 -2.8817200151833804]
[0.48250056255314533 0.6820787268235271 … -1.056153225252739 -2.8817200151833804]
[0.48250056255314533 -0.27498517428669944 … -1.056153225252739 -2.8817200151833804]
[0.48250056255314533 -0.27498517428669944 … -1.056153225252739 -2.442418853961172]
[0.8241432396359415 0.2248767306419267 … 0.5181826296687606 -2.051656876550303]
[1.499937891807123 0.2248767306419267 … 0.5181826296687606 -2.051656876550303]Now that our sampling is complete we can evaluate the potential energy expectation value. Here we use the @estimate macro which will evaluate the given function for every configuration inside chain and return the average. Here we can see that the energy we obtain closely matches that predicted by the equipartition theorem.
julia> Estimators.@estimate potential_energy(sim, chain)38.88040511152187julia> sim.temperature / 2 * 537.5
Example 2
Here, we obtain a thermal distribution in a simple model system with some additional tweaks to try and sample a larger configuration space.
First we set up the system in the usual way, here we're using an NO molecule with a harmonic interaction between the atoms. Notice that we use Unitful.jl to specify the temperature.
using Unitful
atoms = Atoms([:N, :O])
model = DiatomicHarmonic(1.0)
sim = Simulation{Classical}(atoms, model; temperature=300u"K")Then we have to specify the parameters for the Monte Carlo simulation and perform the sampling. Δ contains the step sizes for each of the species, R0 the initial geometry and samples the number of configurations we want to obtain.
AdvancedMH.jl provides some additional options to control the sampling process. To (hopefully) include a larger configuration space in the final results, we set a thinning of 10, meaning that we only keep every 10th proposed configuration. In addition, we discard the first 100 samples, since our initial configuration might not lie in the equilibrium. This is set with discard_initial. Further explanations of the keyword arguments can be found in the AbstractMCMC.jl documentation.
Δ = Dict(:N => 0.1, :O => 0.1)
R0 = [1.0 0.0; 0.0 0.0; 0.0 0.0]
samples = 1000
output = InitialConditions.ThermalMonteCarlo.run_advancedmh_sampling(sim, R0, samples, Δ; movement_ratio=0.5, thinning=10, discard_initial=100)Sampling 0%| | ETA: N/A
Sampling 0%|▎ | ETA: 0:00:07
Sampling 1%|▍ | ETA: 0:00:03
Sampling 1%|▋ | ETA: 0:00:03
Sampling 2%|▉ | ETA: 0:00:02
Sampling 2%|█ | ETA: 0:00:02
Sampling 3%|█▎ | ETA: 0:00:01
Sampling 3%|█▌ | ETA: 0:00:01
Sampling 4%|█▋ | ETA: 0:00:01
Sampling 4%|█▉ | ETA: 0:00:01
Sampling 5%|██▏ | ETA: 0:00:01
Sampling 5%|██▎ | ETA: 0:00:01
Sampling 6%|██▌ | ETA: 0:00:01
Sampling 6%|██▊ | ETA: 0:00:01
Sampling 7%|██▉ | ETA: 0:00:01
Sampling 7%|███▏ | ETA: 0:00:01
Sampling 8%|███▍ | ETA: 0:00:01
Sampling 8%|███▌ | ETA: 0:00:01
Sampling 9%|███▊ | ETA: 0:00:01
Sampling 9%|████ | ETA: 0:00:00
Sampling 10%|████▏ | ETA: 0:00:00
Sampling 10%|████▍ | ETA: 0:00:00
Sampling 11%|████▋ | ETA: 0:00:00
Sampling 11%|████▊ | ETA: 0:00:00
Sampling 12%|█████ | ETA: 0:00:00
Sampling 12%|█████▎ | ETA: 0:00:00
Sampling 13%|█████▍ | ETA: 0:00:00
Sampling 13%|█████▋ | ETA: 0:00:00
Sampling 14%|█████▉ | ETA: 0:00:00
Sampling 14%|██████ | ETA: 0:00:00
Sampling 15%|██████▎ | ETA: 0:00:00
Sampling 15%|██████▌ | ETA: 0:00:00
Sampling 16%|██████▋ | ETA: 0:00:00
Sampling 16%|██████▉ | ETA: 0:00:00
Sampling 17%|███████▏ | ETA: 0:00:00
Sampling 17%|███████▎ | ETA: 0:00:00
Sampling 18%|███████▌ | ETA: 0:00:00
Sampling 18%|███████▊ | ETA: 0:00:00
Sampling 19%|███████▉ | ETA: 0:00:00
Sampling 19%|████████▏ | ETA: 0:00:00
Sampling 20%|████████▍ | ETA: 0:00:00
Sampling 20%|████████▌ | ETA: 0:00:00
Sampling 21%|████████▊ | ETA: 0:00:00
Sampling 21%|█████████ | ETA: 0:00:00
Sampling 22%|█████████▏ | ETA: 0:00:00
Sampling 22%|█████████▍ | ETA: 0:00:00
Sampling 23%|█████████▋ | ETA: 0:00:00
Sampling 23%|█████████▊ | ETA: 0:00:00
Sampling 24%|██████████ | ETA: 0:00:00
Sampling 24%|██████████▎ | ETA: 0:00:00
Sampling 25%|██████████▍ | ETA: 0:00:00
Sampling 25%|██████████▋ | ETA: 0:00:00
Sampling 26%|██████████▉ | ETA: 0:00:00
Sampling 26%|███████████ | ETA: 0:00:00
Sampling 27%|███████████▎ | ETA: 0:00:00
Sampling 27%|███████████▌ | ETA: 0:00:00
Sampling 28%|███████████▋ | ETA: 0:00:00
Sampling 28%|███████████▉ | ETA: 0:00:00
Sampling 29%|████████████▏ | ETA: 0:00:00
Sampling 29%|████████████▎ | ETA: 0:00:00
Sampling 30%|████████████▌ | ETA: 0:00:00
Sampling 30%|████████████▊ | ETA: 0:00:00
Sampling 31%|████████████▉ | ETA: 0:00:00
Sampling 31%|█████████████▏ | ETA: 0:00:00
Sampling 32%|█████████████▍ | ETA: 0:00:00
Sampling 32%|█████████████▌ | ETA: 0:00:00
Sampling 33%|█████████████▊ | ETA: 0:00:00
Sampling 33%|██████████████ | ETA: 0:00:00
Sampling 34%|██████████████▏ | ETA: 0:00:00
Sampling 34%|██████████████▍ | ETA: 0:00:00
Sampling 35%|██████████████▋ | ETA: 0:00:00
Sampling 35%|██████████████▊ | ETA: 0:00:00
Sampling 36%|███████████████ | ETA: 0:00:00
Sampling 36%|███████████████▎ | ETA: 0:00:00
Sampling 37%|███████████████▍ | ETA: 0:00:00
Sampling 37%|███████████████▋ | ETA: 0:00:00
Sampling 38%|███████████████▉ | ETA: 0:00:00
Sampling 38%|████████████████ | ETA: 0:00:00
Sampling 39%|████████████████▎ | ETA: 0:00:00
Sampling 39%|████████████████▌ | ETA: 0:00:00
Sampling 40%|████████████████▋ | ETA: 0:00:00
Sampling 40%|████████████████▉ | ETA: 0:00:00
Sampling 41%|█████████████████▏ | ETA: 0:00:00
Sampling 41%|█████████████████▎ | ETA: 0:00:00
Sampling 42%|█████████████████▌ | ETA: 0:00:00
Sampling 42%|█████████████████▊ | ETA: 0:00:00
Sampling 43%|█████████████████▉ | ETA: 0:00:00
Sampling 43%|██████████████████▏ | ETA: 0:00:00
Sampling 44%|██████████████████▍ | ETA: 0:00:00
Sampling 44%|██████████████████▌ | ETA: 0:00:00
Sampling 45%|██████████████████▊ | ETA: 0:00:00
Sampling 45%|███████████████████ | ETA: 0:00:00
Sampling 46%|███████████████████▏ | ETA: 0:00:00
Sampling 46%|███████████████████▍ | ETA: 0:00:00
Sampling 47%|███████████████████▌ | ETA: 0:00:00
Sampling 47%|███████████████████▊ | ETA: 0:00:00
Sampling 48%|████████████████████ | ETA: 0:00:00
Sampling 48%|████████████████████▏ | ETA: 0:00:00
Sampling 49%|████████████████████▍ | ETA: 0:00:00
Sampling 49%|████████████████████▋ | ETA: 0:00:00
Sampling 50%|████████████████████▊ | ETA: 0:00:00
Sampling 50%|█████████████████████ | ETA: 0:00:00
Sampling 51%|█████████████████████▎ | ETA: 0:00:00
Sampling 51%|█████████████████████▍ | ETA: 0:00:00
Sampling 52%|█████████████████████▋ | ETA: 0:00:00
Sampling 52%|█████████████████████▉ | ETA: 0:00:00
Sampling 53%|██████████████████████ | ETA: 0:00:00
Sampling 53%|██████████████████████▎ | ETA: 0:00:00
Sampling 54%|██████████████████████▌ | ETA: 0:00:00
Sampling 54%|██████████████████████▋ | ETA: 0:00:00
Sampling 55%|██████████████████████▉ | ETA: 0:00:00
Sampling 55%|███████████████████████▏ | ETA: 0:00:00
Sampling 55%|███████████████████████▎ | ETA: 0:00:00
Sampling 56%|███████████████████████▌ | ETA: 0:00:00
Sampling 56%|███████████████████████▊ | ETA: 0:00:00
Sampling 57%|███████████████████████▉ | ETA: 0:00:00
Sampling 57%|████████████████████████▏ | ETA: 0:00:00
Sampling 58%|████████████████████████▍ | ETA: 0:00:00
Sampling 58%|████████████████████████▌ | ETA: 0:00:00
Sampling 59%|████████████████████████▊ | ETA: 0:00:00
Sampling 59%|█████████████████████████ | ETA: 0:00:00
Sampling 60%|█████████████████████████▏ | ETA: 0:00:00
Sampling 60%|█████████████████████████▍ | ETA: 0:00:00
Sampling 61%|█████████████████████████▋ | ETA: 0:00:00
Sampling 61%|█████████████████████████▊ | ETA: 0:00:00
Sampling 62%|██████████████████████████ | ETA: 0:00:00
Sampling 62%|██████████████████████████▎ | ETA: 0:00:00
Sampling 63%|██████████████████████████▍ | ETA: 0:00:00
Sampling 63%|██████████████████████████▋ | ETA: 0:00:00
Sampling 64%|██████████████████████████▉ | ETA: 0:00:00
Sampling 64%|███████████████████████████ | ETA: 0:00:00
Sampling 65%|███████████████████████████▎ | ETA: 0:00:00
Sampling 65%|███████████████████████████▌ | ETA: 0:00:00
Sampling 66%|███████████████████████████▋ | ETA: 0:00:00
Sampling 66%|███████████████████████████▉ | ETA: 0:00:00
Sampling 67%|████████████████████████████▏ | ETA: 0:00:00
Sampling 67%|████████████████████████████▎ | ETA: 0:00:00
Sampling 68%|████████████████████████████▌ | ETA: 0:00:00
Sampling 68%|████████████████████████████▊ | ETA: 0:00:00
Sampling 69%|████████████████████████████▉ | ETA: 0:00:00
Sampling 69%|█████████████████████████████▏ | ETA: 0:00:00
Sampling 70%|█████████████████████████████▍ | ETA: 0:00:00
Sampling 70%|█████████████████████████████▌ | ETA: 0:00:00
Sampling 71%|█████████████████████████████▊ | ETA: 0:00:00
Sampling 71%|██████████████████████████████ | ETA: 0:00:00
Sampling 72%|██████████████████████████████▏ | ETA: 0:00:00
Sampling 72%|██████████████████████████████▍ | ETA: 0:00:00
Sampling 73%|██████████████████████████████▋ | ETA: 0:00:00
Sampling 73%|██████████████████████████████▊ | ETA: 0:00:00
Sampling 74%|███████████████████████████████ | ETA: 0:00:00
Sampling 74%|███████████████████████████████▎ | ETA: 0:00:00
Sampling 75%|███████████████████████████████▍ | ETA: 0:00:00
Sampling 75%|███████████████████████████████▋ | ETA: 0:00:00
Sampling 76%|███████████████████████████████▉ | ETA: 0:00:00
Sampling 76%|████████████████████████████████ | ETA: 0:00:00
Sampling 77%|████████████████████████████████▎ | ETA: 0:00:00
Sampling 77%|████████████████████████████████▌ | ETA: 0:00:00
Sampling 78%|████████████████████████████████▋ | ETA: 0:00:00
Sampling 78%|████████████████████████████████▉ | ETA: 0:00:00
Sampling 79%|█████████████████████████████████▏ | ETA: 0:00:00
Sampling 79%|█████████████████████████████████▎ | ETA: 0:00:00
Sampling 80%|█████████████████████████████████▌ | ETA: 0:00:00
Sampling 80%|█████████████████████████████████▊ | ETA: 0:00:00
Sampling 81%|█████████████████████████████████▉ | ETA: 0:00:00
Sampling 81%|██████████████████████████████████▏ | ETA: 0:00:00
Sampling 82%|██████████████████████████████████▍ | ETA: 0:00:00
Sampling 82%|██████████████████████████████████▌ | ETA: 0:00:00
Sampling 83%|██████████████████████████████████▊ | ETA: 0:00:00
Sampling 83%|███████████████████████████████████ | ETA: 0:00:00
Sampling 84%|███████████████████████████████████▏ | ETA: 0:00:00
Sampling 84%|███████████████████████████████████▍ | ETA: 0:00:00
Sampling 85%|███████████████████████████████████▋ | ETA: 0:00:00
Sampling 85%|███████████████████████████████████▊ | ETA: 0:00:00
Sampling 86%|████████████████████████████████████ | ETA: 0:00:00
Sampling 86%|████████████████████████████████████▎ | ETA: 0:00:00
Sampling 87%|████████████████████████████████████▍ | ETA: 0:00:00
Sampling 87%|████████████████████████████████████▋ | ETA: 0:00:00
Sampling 88%|████████████████████████████████████▉ | ETA: 0:00:00
Sampling 88%|█████████████████████████████████████ | ETA: 0:00:00
Sampling 89%|█████████████████████████████████████▎ | ETA: 0:00:00
Sampling 89%|█████████████████████████████████████▌ | ETA: 0:00:00
Sampling 90%|█████████████████████████████████████▋ | ETA: 0:00:00
Sampling 90%|█████████████████████████████████████▉ | ETA: 0:00:00
Sampling 91%|██████████████████████████████████████▏ | ETA: 0:00:00
Sampling 91%|██████████████████████████████████████▎ | ETA: 0:00:00
Sampling 92%|██████████████████████████████████████▌ | ETA: 0:00:00
Sampling 92%|██████████████████████████████████████▊ | ETA: 0:00:00
Sampling 93%|██████████████████████████████████████▉ | ETA: 0:00:00
Sampling 93%|███████████████████████████████████████▏ | ETA: 0:00:00
Sampling 94%|███████████████████████████████████████▍ | ETA: 0:00:00
Sampling 94%|███████████████████████████████████████▌ | ETA: 0:00:00
Sampling 95%|███████████████████████████████████████▊ | ETA: 0:00:00
Sampling 95%|████████████████████████████████████████ | ETA: 0:00:00
Sampling 96%|████████████████████████████████████████▏ | ETA: 0:00:00
Sampling 96%|████████████████████████████████████████▍ | ETA: 0:00:00
Sampling 97%|████████████████████████████████████████▋ | ETA: 0:00:00
Sampling 97%|████████████████████████████████████████▊ | ETA: 0:00:00
Sampling 98%|█████████████████████████████████████████ | ETA: 0:00:00
Sampling 98%|█████████████████████████████████████████▎| ETA: 0:00:00
Sampling 99%|█████████████████████████████████████████▍| ETA: 0:00:00
Sampling 99%|█████████████████████████████████████████▋| ETA: 0:00:00
Sampling 100%|█████████████████████████████████████████▉| ETA: 0:00:00
Sampling 100%|██████████████████████████████████████████| Time: 0:00:00We can calculate the distance between each atom and plot the bond length throughout the sampling.
using LinearAlgebra
plot([norm(R[:,1] .- R[:,2]) for R in output])
The result of this simulation seamlessly interfaces with the DynamicalDistribution presented in the previous section and output can be readily passed to provide the position distribution. The Monte Carlo sampling does not include velocities but these can be readily obtained from the Maxwell-Boltzmann distribution.