Monte Carlo Simulations

The mc.jl function implements the Monte Carlo method.

Monte Carlo simulations have a completely different principle from dynamics simulations, but they are expected to sample the same set of configurations if the thermodynamic conditions are the same. Here we will perform a Monte Carlo simulation and verify how similar it is to molecular dynamics simulations.

Unlike MD, MC has no time. There is a generation of consecutive random positions, which are accepted or rejected according to the Metropolis criterion:

If $V(\vec{x}_j) \leqslant V(\vec{x}_i)$, then $P(i\to j) = 1$

If $V(\vec{x}_j) > V(\vec{x}_i)$, then $P(i\to j) = e^{-(V_j-V_i)/kT}$

The second criterion is numerically satisfied by comparing the result of $e^{-(V_j-V_i)/kT}$ with a random number drawn between 0 and 1. In our examples, $kT=0.6$.

This procedure generates a sequence of configurations, which in practice are correlated because new configurations are generally generated by perturbations of the previous configurations. The perturbations must be chosen to minimize correlation, while keeping the acceptance rate reasonable. Acceptance rates of around 20 to 30\% are considered ideal.

7.1. Running

Run the mc function. We will vary the magnitude of the perturbations. The positional perturbations are Gaussian, and the input magnitude is the standard deviation. The number of steps corresponds to the number of new structures, not necessarily accepted, that are generated:

julia> out = mc(sys,Options(alpha=0.05,nsteps=50_000));

For 50,000 steps, try different perturbations, until an acceptance rate of around 30\% is obtained. (Something close to $0.08\textrm{\AA}$).

Once the perturbation is chosen, run the program with 200,000 steps, which means approximately 60,000 steps will be accepted (for a 30\% acceptance rate).

Observe the evolution of the potential energy by making plots with:

plot(out,ylim=[-100,100], label="Potential", xlabel="step")

In this case it does not make sense to show the kinetic energy, which is not defined since the particles do not have actual velocities. Save the plot to a pdf file for later comparison, using:

julia> savefig("./mc.pdf")

Observe the trajectory (with the same commands as before).

7.2. Complete summary code

using FundamentosDMC, Plots
sys = System(n=100,sides=[100,100])
minimize!(sys)
out = mc(sys,Options(alpha=0.05,nsteps=50_000))
plot(out,ylim=[-100,100], label="Potential", xlabel="step")