How to choose between Molecular Dynamics and Monte Carlo when beginning to simulate either equilibrium or non-equilibrium systems?

Question

Answers to Are MD and MC both able to study both equilibrium and non-equilibrium systems? include lots of examples pointing to "yes".

I'm currently faced with a task of modelings small islands of 2D materials on a crystal surface - how they rotate, translate, grow, and settle down into certain rotational and translational equilibrium orientations. Right now I use Hooke's law for bonds and a 2D periodic surface potential.

For me, I can simply pop in any solver I want; for Molecular Dynamics I use Python's solve_ivp and have a Python implementation of FIRE and will be coding a FIRE method in C.

For Monte Carlo I have an ugly, slow Python implementation and there are some (presumably faster) packages out there I could try.

My question is more general, but you can use that context if it helps an answer.

Question: What are the differentiating factors that tell us how to choose between Molecular Dynamics and Monte Carlo when beginning to simulate either equilibrium or non-equilibrium systems?

Answer 1

The main difference between MD and MC is, MD is a deterministic method, and MC is a probabilistic method.

So, MC is unable to simulate something which needs total accuracy. For instance, if you want to simulate the movement of a spacecraft, MC cannot do that. Or, rather it would be highly unreliable.

Another thing to consider is, it is very hard (research is ongoing) to parallelize MC algorithms. Therefore, if you want to run MC in GPU or in a distributed system, it would be a very very difficult task.

On the other hand, MD cannot simulate on-lattice models. I.e., the atoms need to be freely movable, e.g., gas or liquid.

Finally, the decision depends on the nature of the problem. If you are simulating, say a protein folding, and you are following the all-atom method rather than Coarse Graining, MD would be very very expensive to run. Because it will take huge computing resources.

Answer 2

The shortest advice I can give you is to search the literature. You must look at the approaches that have been previously tried for your system, especially because these only tend to be published when they work reasonably well. Without doing that you are just following the advice of strangers on the Internet, and while that's perfectly acceptable when deciding where to eat for lunch, that's very difficult to justify in trying to seriously produce scientific knowledge.

---

Almost exactly a year ago I put this picture on MMSE, from Allen and Tildesley's essential textbook:

It is very tempting to start "at the top" of this diagram and say "Let's simulate surface hopping in vapour deposition!" or whatever your phenomenon is. But you need to start at the bottom and ask yourself:

• Do I have experimental results to compare to a model of my phenomenon of interest?
• Do I have a theoretical framework to explain what I see in a model of my phenomenon of interest?

Ask in addition what the precision of those results and frameworks are,¹ because that in turn determines the amount of simulation data you need -- the longer (shorter) you can run your simulation, the more (less) precise your model estimate will be.

You also face the problem of how to efficiently use your computational time. You need to search the literature! There are packages out there that are thousands of times more efficient than a Python solve_ivp script. You could spend a day writing up a Python implementation and then find that it takes a year to run your simulations -- or you could spend a month learning a package like LAMMPS or GROMACS and then run your simulations in a week. Which sounds better?

But most importantly, you must make sure your model can efficiently calculate theoretical details that can be compared with your experiments. Take, again, surface hopping during vapour deposition. You could just dump a bunch of particles on a surface and run a NumPy integrator and get a diffusivity. But when your experimental value turns out to be a hundred times smaller, what are you going to change? The model masses? The bond strength? What if, it turns out, your model was already as accurate as it could be for its simplicity, and you spend months tweaking parameters only to find out that every tweak makes things worse? And what if it turns out you had the simplest of typos (while converting units, let's say -- the bane of every molecular modeller), and you had the right answer all along?

That's why you must search the literature, and if you are new to the field, you should start by trying to replicate a known study. This also familiarises you with the process of comparing theory to model to experiment.

For what it's worth, in your situation, I'd immediately search for "kinetic Monte Carlo simulation of surface deposition", for which Google returns me an excellent-looking review. But a kMC simulation (as they're known for short) requires you to input sensible reaction rates. You'd probably need to get those from a quantum calculation of some sort, probably a nudged-elastic band calculation of the individual reaction energy profiles using density functional theory. Note that the latter does not fall in your "MC or MD" question, and the former might not have fallen in your umbrella MC category either (most kMC people I know think of themselves as doing something more specialised than regular molecular MC, using something like the GOMC package, and they're not wrong).

But again -- you shouldn't be basing your scientific research on what a stranger says on MMSE. Go look through the literature and see what the field is doing. There is no shortcut for that.

---

¹ If I want to explain why gravitational acceleration is 9.8 meters per second squared, I just need to know Newton's law of gravity and the mass and radius of the earth. If I want to explain its value to nine decimal places, suddenly I have to learn general relativity too.