Writing a Monte Carlo Simulation for a Lennard Jones Fluid
Today in lecture, we learned how to use Monte Carlo to evaluate statistical mechanics and to get approximations for the thermodynamic properties of a chemical system.
These properties are expressed as high dimensional integrals:
$$ ⟨Q⟩ = \int_V Q(r^N)ρ(r^N)dr^N $$
In this equation, $ Q $ is the potential energy of the system, and $ ρ(r^N) $ stands for probability density.
Why Use Monte Carlo?
According to the idea of importance sampling, rather than evaluating the integral for every possible state, it is more efficient to evaluate the integral for situations that are more likely to actually exist.
Therefore, we will be basing these calculations on the equilibrium probability density which will help us do just that.
This is where the Metrapolis Monte Carlo simulation comes in.
In 1953, Nicholas Metropolis proposed the Metrapolis-Hasting algorithm, and W.K. Hasting developed on Metrapolis' idea in 1970. Hence we have the Metrapolis-Hasting algorithm:
$$ P_{acc} (m->n) = min [1, e^{βΔU}] $$
$ P_{acc} $ stands for the probability of accepting a move signified by a coordinate change.
$$ β = \frac{1}{T} $$
Now, Monte Carlo will let us generate the energy changes from randomly moving an atom a certain distance.
In our simulation, moves that decrease energy will always be accepted, while moves that are positive or zero are only sometimes accepted.
Radial Distribution Function (RDF)
Radial Distribution Function refers to the probability of a particle being a certain distance (r) away from a given particle.
Solids have the highest RDF peaks as the atoms of solids are packed tighly together and are in lattice positions.
On the other hand, gases are free flowing and do not have a strong constant structure like solids, so the RDF of gas does not peak as much as it does for solids. In addition, the RDF of gases eventually levels out to 1. This is because the probability of finding two particles a certain distance apart becomes constant, and this constant is based on the bulk density, which happens to be 1.
Liquids, like gas, do not stay in a tight constant structure. However, liquids do have RDF peaks like solids though to a lesser degree, with the RDF peak of the closest distance being the highest.
For molecular simulations, each individual RDF is usually averaged. As the RDF may be different at different areas of the system due to factors such as temperature, average each individual fram allows for a value close to the overall equilabrium of a system to be found.