Chemistry · Lesson

Molecular dynamics, Monte Carlo and sampling

Read time integration, equilibrium sampling, state weights and correlated observations under their declared models.

By the end, you can

Apply one velocity-Verlet step and convert model time.
Include reverse proposals in a Metropolis–Hastings acceptance account.
Combine microscopic energies and degeneracies without double counting.
Calculate finite mean variance under a stated covariance model.

Molecular dynamics follows a specified time-evolution model

In a classical molecular-dynamics model, positions and velocities evolve under assigned forces and integration rules. Initial conditions, boundaries and any temperature or pressure control are part of the simulation. A time step must be appropriate to the represented motions and numerical method.

A stable-looking trajectory does not automatically establish accurate dynamics or representative equilibrium sampling. Conservation or other expected behaviour should be checked under the actual ensemble and algorithm.

A time step and a sampling step have different meanings

Velocity Verlet updates position from the current velocity and acceleration, then uses both endpoint accelerations for the velocity. Its time increment belongs to the dynamical model. Changing the increment changes numerical error as well as the time associated with a fixed number of steps. The single oscillator step below does not certify long trajectories.

A Monte Carlo proposal instead changes a state through an algorithm. In the two-state example the proposals are asymmetric, so the reverse proposal probability enters the acceptance ratio. Detailed balance checks that the assigned target is stationary. Suitable access to the relevant states and adequate sampling remain separate requirements.

Rejected proposals keep the current state. Removing them from an average changes the sampling account. A number of proposals has no supplied conversion into physical molecular time in this example.

Monte Carlo can sample states without representing physical time

A Monte Carlo procedure proposes changes and accepts or rejects them according to a specified rule. A valid equilibrium scheme targets a defined distribution under its conditions. Proposal probabilities and acceptance rules must fit together.

The number of Monte Carlo steps is generally different from elapsed physical time. A rejected proposal is part of the sampling algorithm, rather than proof that a real molecular event did not occur. Check that the procedure reaches the relevant states instead of only one region.

A population and its estimated mean need the correct weights

An energy row can stand for one microscopic state or for a group. In the assigned two-group example, four equivalent B states contribute four equal Boltzmann terms. Once those terms are included in a macrostate free energy, multiplying by the same degeneracy again would double count.

An observable average then needs a sampling-uncertainty account. Positive correlation between nearby observations increases the variance of the mean relative to independent observations with the same marginal variance. The finite sum below includes the actual number of observation pairs at each separation.

The variance-equivalent effective count summarizes one assigned covariance law. A long series can still fail to sample a relevant state, and a variance formula conditional on stationarity does not prove stationarity. Preserve the initial conditions, equilibration choices, boundaries, observations and transformations used to assess a real simulation.

Correlation increases the standard error of the mean

Standard error / kJ mol⁻¹

Number of observations N

Assigned correlated observations
Independent reference
— Finite covariance-model result

Underlying data

Number of observations N	Standard error / kJ mol⁻¹
20	0.54884552846152 · Assigned correlated observations
80	0.287174838304302 · Assigned correlated observations
160	0.20452782855411442 · Assigned correlated observations
20	0.2683281572999747 · Independent reference
80	0.13416407864998736 · Independent reference
160	0.09486832980505137 · Independent reference

Original exact stationary covariance model with σ=1.20 kJ mol⁻¹ and ρ=0.650. Circles and dashed curve use the finite correlated-mean variance; squares use independent observations with the same marginal variance. No actual trajectory or interval coverage is represented.

An ensemble average needs the correct state weights

A population calculation needs the state energies or appropriate free energies, degeneracies and temperature convention. A higher-energy state can still contribute, and several states can collectively outweigh one lower-energy state.

Distinguish a microscopic energy from a macroscopic free energy that already includes a state’s internal multiplicity. Counting degeneracy again can double-count that contribution. State what each row of a model represents.

Many saved frames are not necessarily many independent samples

Nearby trajectory frames can be correlated. Longer sampling or more saved frames can characterize a trace without proportionally increasing the independent information about an equilibrium average. Rare transitions and disconnected regions can limit coverage.

Use suitable equilibration, convergence, correlation and replicate checks for the intended estimate. The original sampling examples below contain invented distributions and observations; they do not prove convergence for an actual molecular simulation.

Worked examples

One integration step and its time account

Use one-dimensional velocity Verlet for a symbolic neutral oscillator: m=4.00×10⁻²⁶ kg, k=8.00 N m⁻¹, x₀=1.00×10⁻¹⁰ m, v₀=−120 m s⁻¹ and Δt=2.00 fs. No thermostat or barostat is present. Calculate one step and the time assigned to 60,000 such steps.

a₀=−kx₀/m=−2.00×10¹⁶ m s⁻². x₁=x₀+Δtv₀+Δt²a₀/2=9.972×10⁻¹¹ m=0.09972 nm.
a₁=−kx₁/m. v₁=v₀+Δt(a₀+a₁)/2=−159.944 m s⁻¹. Initial total energy mv₀²/2+kx₀²/2 is 4.028800×10⁻²⁰ J; one-step energy is 4.028795526272×10⁻²⁰ J.
The change is −4.473728×10⁻²⁶ J for this step. It checks a numerical result in this assigned model; it proves neither long-time stability nor accurate molecular dynamics.
60,000×2.00 fs=120 ps of model time. That account does not depend on how frequently frames are saved, and it does not assign physical time to an unrelated Monte Carlo proposal count.

x₁=0.09972 nm, v₁=−159.944 m s⁻¹; 60,000 steps correspond to 120 ps in this time-evolution model.

An asymmetric proposal needs its reverse probability

For two symbolic configurations of the same inventory, assign target probabilities π_A=0.750 and π_B=0.250. Propose A→B with probability 0.600 and B→A with 0.200; the remaining proposal probability stays in each state. Apply Metropolis–Hastings.

Acceptance α_AB=min[1,π_Bq_BA/(π_Aq_AB)]=min[1,0.250×0.200/(0.750×0.600)]=1/9. α_BA=1.
The actual off-diagonal transition probabilities are P_AB=0.600/9=0.0666667 and P_BA=0.200. Staying probabilities are 0.9333333 and 0.800, including rejected proposals.
Equilibrium fluxes are π_AP_AB=0.0500 and π_BP_BA=0.0500. Their equality checks detailed balance and the assigned invariant target in this connected finite chain.
Ignoring the unequal reverse proposals would target a different distribution. Rejections retain the current state and belong in the sampling average. This demonstration assigns no physical time or real molecular convergence.

Acceptances 1/9 and one; both target-weighted transition fluxes are 0.0500 per step.

Several higher-energy states can outweigh one lower-energy state

At 310.0 K assign macrostate A one microscopic state at 0 kJ mol⁻¹ and macrostate B four states, each at 3.40 kJ mol⁻¹. No other states exist in this toy account. Find the B population.

Use the molar gas constant R=N_Ak_B=8.314462618… J mol⁻¹ K⁻¹, with Avogadro constant N_A and Boltzmann constant k_B. The relative weights are w_A=1 and w_B=4exp[−3.40/(0.008314462618…×310.0)]=1.069486053.
P_B=w_B/(1+w_B)=0.516788239 and P_A=0.483211761. Mean assigned molar energy is 3.40P_B=1.75708001 kJ mol⁻¹.
For these state groups, ΔG_macro=3.40−RTln4=−0.173150719 kJ mol⁻¹. Each B microstate is higher in energy, yet the whole B group is slightly more populated.
If an input free energy already included this internal multiplicity, multiplying by four again would count it twice. These assigned states are not measured molecular conformers.

P_B=0.516788239; the fourfold multiplicity changes the macrostate comparison.

Correlated frames change the uncertainty of a mean

Assign a stationary observable with N=80 equally spaced observations, marginal standard deviation σ=1.20 kJ mol⁻¹ and covariance Cov(Y_i,Y_j)=σ²ρ^|i−j|, ρ=0.650. Calculate the standard error of the sample mean under this exact stipulated covariance law.

For a mean of N observations, Var(mean)=σ²/N[1+2Σ_(k=1)^(N−1)(1−k/N)ρ^k]. The finite factor is 4.581632653.
SE=1.20√(4.581632653/80)=0.287174838 kJ mol⁻¹. If the 80 observations were independent, SE would be 1.20/√80=0.134164079 kJ mol⁻¹.
Define a variance-equivalent effective count N_eff=N/factor=17.4610245. It is fractional and describes this mean-variance comparison; it is not a guaranteed count of independent frames.
The covariance is assigned rather than estimated. Its exact variance does not establish equilibration, rare-state coverage, interval coverage or convergence of an actual trajectory.

Model-based SE=0.287174838 kJ mol⁻¹ and variance-equivalent N_eff=17.4610245.

Try it yourself

Use dt = 4.00 fs with the same initial state and 60000 steps for the elapsed-time account. Find the one-step state and elapsed time.
Show the answer and reasoning
x₁=0.09936 nm, v₁=−199.744 m s⁻¹; energy change −4.0828928×10⁻²⁵ J; 60,000 steps give 240 ps.
The integration state and numerical energy error change with dt; elapsed model time is the number of steps times dt.
Change only proposal B -> A to 0.500, with selfB = 0.500. Find both acceptances, off-diagonal transitions and equilibrium fluxes.
Show the answer and reasoning
α_AB=5/18, α_BA=1; P_AB=1/6, P_BA=1/2; both fluxes 0.125.
The reverse proposal probability belongs in the acceptance ratio. Detailed balance is retained after the adjusted acceptance.
Change only B degeneracy from 4 to 2. Find P_B and the macrostate free-energy difference at 310.0 K.
Show the answer and reasoning
P_B=0.348425122; ΔG_macro=+1.61342464 kJ mol⁻¹.
The energy of each microscopic B state stays the same while the number of states changes.
Double only the number of observations to 160, keeping the same assigned covariance law and marginal sigma. Find N_eff and the standard error.
Show the answer and reasoning
N_eff=34.4237102; SE=0.204527829 kJ mol⁻¹.
More observations reduce mean variance, but correlated observations remain less informative than the same number of independent draws.

What to carry forward

Integration time follows the stated dynamical time step.
Equilibrium sampling requires compatible proposal and acceptance probabilities.
Multiplicity must be included once on the defined state basis.
Correlated observations and incomplete state coverage limit different aspects of a simulation.

References

Teaching resources and reference material supporting this lesson’s explanations and examples.