Intro to Molecular Dynamics: Coding MD From Scratch

This is a brief introduction to how MD simulations work: essentially numerically solving Newton’s equations for a bunch of interacting objects. It is a bit different from my usual videos because it is just something I was screwing around with for fun and thought it might be interesting to post. If you notice any issues with the code or the math please let me know. As I said, I was just screwing around so it’s not as extensively checked as publication-quality research.

Correction: When running the water simulations, there was an error in the epsilon values I was using. The videos I included are not too far off but I thought I should mention it.

Correction: As pointed out by MCMC271, simply rescaling the velocities like I did does not produce the correct probability distributions for the NVT ensemble and is thus not a good way to keep the temperature constant. I should have stated this in the video and pointed the viewer to better thermostats.

References:
(Still not sure if to include the code so it may or may not be public)

In case I do not get around to making a video on using available MD software, I recommend looking into LAMMPS, NAMD or GROMACS:
LAMMPS:
I particularly recommend the .pdf files numbered SoR 01, 05 and 11
NAMD:
GROMACS:
If you feel like paying for software, you should also check out Materials Studio and/or COMSOL.

Chapters:
0:00 Hello
0:39 Newton's equations
1:21 Code
2:30 Visualization (matplotlib)
3:10 Boundary conditions (periodic)
4:35 BCs (reflecting)
5:55 Visualization (OVITO)
7:45 Lennard-Jones interactions
13:14 Periodic BC interaction discussion
14:48 Particle types
15:12 Microcanonical (NVE) ensemble
16:33 Canonical ensemble (fixing T)
19:26 Bond potentials
22:17 Bond angles
26:32 Dihedral angles
27:11 Electrostatics
28:28 Combining potentials
29:33 Polymers
30:08 Potential cutoff
30:53 Gravity
31:28 Summary

[uoy]

Extremely high quality video. Great job!

[niconeuman]

Thank you very much! I started a python MD code from an older video that had no follow up and yours contains everything I need to follow! This is very useful for me! Thanks again!

[PunmasterSTP]

As far as tutorials go, this was incredible, and this video (and channel) deserves a much wider audience. Thanks for putting all this thought and effort into making it!

[MCMC]

This is a very helpful video for understanding MD simulation from scratch, even with wonderful visualization! I would like to appreciate your efforts. I'd like add a minor comment that velocity rescaling is not a suitable method for maintaining the probability function of the canonical ensemble (NVT). To ensure the system follows the probability function of the canonical ensemble, one should use other thermostats like Langevin or Nose-Hoover thermostats.

[thehungmai]

A good video. I'm freshman at University. I'm so grateful video because It helped in my study. Thank you and hopefully that you will continue to launch video like this

[Rexd]

I have no idea how to do any of this math but I am very interested after watching this and attempting to work with NAMD and VMD

[gabrielgoetten]

Hi Russel, excellent video! Thank you for going into the basics, and the way you did is very unique. However I do have a question as to why not share the code? Are you afraid of someone stealing your original idea ?

[uoy]

A few questions if you don't mind:
When you mean a better numerical scheme, the verlet algorithm would be the standard choice am I correct?
I've gotten into MD just this week hoping to use it to see the interaction between the anion and cation of an Ionic liquid, these are fairly complex molecules, do you think it's worth building an MD from scratch for this or would it be better to use GROMACS or LAMMPS? Especially if i'm only looking for a qualitative understanding and accuracy is not a priority

[bomcimtube]

Hello, what is the software do you use in your research to build and model polymeric systems?

[Черезтерникзвёздам-кр]

Hello. May i ask a question? When you write formula gradPij=-24e[...]*(xij * x" yij * y" zij * z"] . xij means xij=( xi-xj ) etc. Will it be more right to use module xij=| xi-xj | ?

[eskicia]

hiya, do you have any platform I can maybe contact you on? I am quite stuck on a project relating to MD for school

[GewoonFinn-]

I am working on this, however your NVE slide is a bit odd not well explained to my understanding, the pseudo code is not really understandable by me. lambda is defined two times different, besides if you start your simulation how do you choose KEold, you set it first equal to KEnew? Lambda might be very big at first, and how is now delta E defined? Probably it is determined by your starting velocity. And I guess starting with zero velocity is not the nicest thing to start with as dividing by zero is pretty unpleasant interesting stuff.