Comparison between MC and MD

We will compare the set of structures generated by molecular dynamics with a set of structures generated by a Monte Carlo simulation.

8.1. Generating a good MD sampling

Run the md_langevin program again with parameters $\lambda=0.01$, for 100 thousand steps, and save the trajectory to a file with a suitable name, such as md.xyz:

julia> md_out = md_langevin(
           sys,
           Options(lambda=0.01,nsteps=100_000,trajectory_file="md.xyz")
       );

8.2. Generating a good MC sampling

Run a long Monte Carlo simulation (300_000 steps)

julia> mc_out = mc(
           sys,
           Options(alpha=0.05,nsteps=300_000,trajectory_file="mc.xyz")
       );

8.3. Comparing potential energies

The first column of md_out contains the potential energy from the molecular dynamics simulation. Let us plot this energy:

julia> plot(md_out[:,1],label="Potential energy- MD",xlabel="step")

And let us add to the same plot the potential energy obtained from the Monte Carlo simulation, which is in mc_out (note the ! at the end of plot!, indicating that the previous plot will be modified):

julia> plot!(mc_out,label="Potential energy- MC",xlabel="step")

Note the similarity, or difference, between the two plots. Keep in mind that in neither of these simulations did we explicitly control the potential energy.

8.4. Radial distribution function

We will compare the average structure obtained using MD with the average structure obtained with MC. For this we will use the $g(r)$ function, called the ``radial distribution function''. This function measures the probability of finding a particle at a distance $r$ from another particle in the real system, relative to that same probability if there were no interactions between the particles.

In our two-dimensional case, the number of particles per unit area is $\rho=n/A$, where $n=100$ is the number of particles and $A=100^2$ is the total area of the simulated system. The expected number of particles in a distance interval between $r$ and $r+\Delta r$ from each particle is therefore $n(r)=\rho A(r)$, where $A(r)$ is the area of a circular shell with inner radius $r$ and outer radius $r+\Delta r$:

We see that $A(r)=\pi (r+\Delta r)^2 - \pi r^2 \approx 2\pi r\Delta r$. Thus, the expected average number of particles would be $n(r)=2\pi r\Delta r\rho$, if there were no interactions.

Interactions cause the number of particles at each distance to differ from a homogeneous distribution. If there are favorable interactions, for example, the probability of finding two particles close together is higher. This particle distribution is one of the most important structural parameters.

8.5. Computing $g(r)$

The program radial_distribution.jl computes, from a trajectory, the function $g(r)=n'(r)/n(r)$, where $n(r)$ is defined above, and $n'(r)$ is the average number of particles effectively observed between $r$ and $r+\Delta r$ in the simulation.

To compute the radial distribution function of the molecular dynamics simulation, run:

julia> rmd, gmd = radial_distribution(sys,"md.xyz")

rmc and rmd are the distances and relative particle densities at each distance, relative to the system density.

Next, let us obtain the $g(r)$ of the Monte Carlo simulation:

julia> rmc, gmc = radial_distribution(sys,"mc.xyz")

Create a plot of the distribution function from the molecular dynamics simulation with:

julia> plot(rmd,gmd,xlabel="r",ylabel="g(r)",label="MD")

Note the increase in local density at short distances, which results from the favorable interactions between particles. Also note that the relative density tends to 1.0 at large distances, when the positional correlation between particles is no longer significant.

Add to the same plot the distribution function obtained from the Monte Carlo simulation:

julia> plot!(rmc,gmc,xlabel="r",ylabel="g(r)",label="MC")

Compare. Did the simulations, with their completely different natures, sample the same structures?

8.2. Complete summary code

using FundamentosDMC, Plots
sys = System(n=100,sides=[100,100])
minimize!(sys)

md_out = md_langevin(
    sys,
    Options(lambda=0.01,nsteps=100_000,trajectory_file="md.xyz")
)
mc_out = mc(
    sys,
    Options(alpha=0.05,nsteps=300_000,trajectory_file="mc.xyz")
)
plot(md_out[:,1],label="Potential energy- MD",xlabel="step")
plot!(mc_out,label="Potential energy- MC",xlabel="step")
savefig("potential_energy.pdf")

rmd, gmd = radial_distribution(sys,"md.xyz")
rmc, gmc = radial_distribution(sys,"mc.xyz")
plot(rmd,gmd,xlabel="r",ylabel="g(r)",label="MD")
plot!(rmc,gmc,xlabel="r",ylabel="g(r)",label="MC")
savefig("gr.pdf")