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.
Thanks for the comment. You're right - I should have said that in the video. I just wanted to illustrate a simple example of a thermostat, but I didn't make it clear that this isn't a good one and doesn't give the correct distributions.
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!
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!
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
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
The verlet algorithm is a good choice. It's a popular choice among the commercial MD simulators. Building your own MD code from scratch is a nice way to learn about MD and get a feel for how it works, but I would not recommend it as a practical method for getting a good MD simulator. Analogy: If you want to understand how cars work, building your own might be fun and informative, but it is not a practical way to get to the supermarket, especially if you can borrow someone else's car for free. I can't speak for you, but GROMACS and LAMMPS are faster, more bug-free, more capable, and generally better than anything I, personally, could write from scratch myself. You asked a month ago, and I only just noticed. I hope your MD journey is going well.
Hey. It depends what I'm doing. I mostly do SCFT calculations and I write my own code for that. I also do coarse grained particle simulations with my own code and other simulations with gromacs or lammps (not as common for me but good options if you want to get into particle simulations)
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 ?
Thanks! The basic ideas are not original and everything one needs to recreate the code is in the video so I'm not afraid of anyone stealing it. I'm just not really used to sharing code publicly so a bit insecure about putting it out there lol. But that's not a good reason so I'll go make it public. Disclaimer that I've put everywhere already: I have not checked it as carefully as I would if, for example, I was collecting data or trying to calculate precise values to publish, so watch out for bugs and/or edge cases.
@@PolymerTheory Alright thanks for the explanation. I think those who will test it will probably use it as a study guide, and will follow your advise, so you do not need to worry....
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 | ?
If you replace xij with |xij| etc. the force always points in the same direction e.g. two interacting particles will have the same force vector on eachother, as opposed to equal and opposite forces. Assuming you mean the grad P equation at about 9:17.
Feel free to email me at the contact email given on my papers - I don't want to post my email address in a RUclips comment because it's more likely to be scraped and hence even more spam than usual.
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.
Typically you start the simulation with a kinetic energy that reflects the temperature that you want it. This, combined with the potential energy defined by the distribution, gives the total energy of the system. E_old and E_new are just the energies before and after you update the kinetic energy, at any given step. Delta E is the difference between the energy of the system and what you want it to be. You can't change the potential energy without moving particles, but you can change the kinetic energy, so in practice it is the difference between the kinetic energy and what you want it to be. That changes at every step. I'm not sure I follow your questions so may not have answered them well.
The code just simulates the dynamics, which only depends on equations of motion and thus derivatives of the potential, not the actual potential itself. If you want to calculate the internal energy then LJPot itself is required. The internal energy (and thus the actual potential) is required for a bunch of stuff you might want to know about the system, like the heat capacity or if you want to do thermodynamic integration or lots of other things. I wanted to include the function/subroutine but since the code just evolves the system in time, it's not useful to calculate LJPot for anything in the scope of the video.
It's easiest to put it right after the particle position is updated. In the example code, I put it at the end of the 'update' routine but it can also be after, depending how you implement it.
Extremely high quality video. Great job!
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
I'm very glad to have inspired you! I hope it goes well 😊
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.
Thanks for the comment. You're right - I should have said that in the video. I just wanted to illustrate a simple example of a thermostat, but I didn't make it clear that this isn't a good one and doesn't give the correct distributions.
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!
My pleasure! I know you like trying out simulation techniques - let me know if you play with the code :)
@@PolymerTheory I will for sure.
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!
Very glad it was useful!
Extremely helpful
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
I'm glad I could be helpful :)
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
The verlet algorithm is a good choice. It's a popular choice among the commercial MD simulators.
Building your own MD code from scratch is a nice way to learn about MD and get a feel for how it works, but I would not recommend it as a practical method for getting a good MD simulator. Analogy: If you want to understand how cars work, building your own might be fun and informative, but it is not a practical way to get to the supermarket, especially if you can borrow someone else's car for free. I can't speak for you, but GROMACS and LAMMPS are faster, more bug-free, more capable, and generally better than anything I, personally, could write from scratch myself.
You asked a month ago, and I only just noticed. I hope your MD journey is going well.
Hello, what is the software do you use in your research to build and model polymeric systems?
Hey. It depends what I'm doing. I mostly do SCFT calculations and I write my own code for that. I also do coarse grained particle simulations with my own code and other simulations with gromacs or lammps (not as common for me but good options if you want to get into particle simulations)
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 ?
Thanks! The basic ideas are not original and everything one needs to recreate the code is in the video so I'm not afraid of anyone stealing it. I'm just not really used to sharing code publicly so a bit insecure about putting it out there lol. But that's not a good reason so I'll go make it public.
Disclaimer that I've put everywhere already: I have not checked it as carefully as I would if, for example, I was collecting data or trying to calculate precise values to publish, so watch out for bugs and/or edge cases.
@@PolymerTheory Alright thanks for the explanation. I think those who will test it will probably use it as a study guide, and will follow your advise, so you do not need to worry....
@@gabrielgoetten Yeah, I think you're right. Just being paranoid.
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 | ?
If you replace xij with |xij| etc. the force always points in the same direction e.g. two interacting particles will have the same force vector on eachother, as opposed to equal and opposite forces. Assuming you mean the grad P equation at about 9:17.
hiya, do you have any platform I can maybe contact you on? I am quite stuck on a project relating to MD for school
Feel free to email me at the contact email given on my papers - I don't want to post my email address in a RUclips comment because it's more likely to be scraped and hence even more spam than usual.
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.
Typically you start the simulation with a kinetic energy that reflects the temperature that you want it. This, combined with the potential energy defined by the distribution, gives the total energy of the system. E_old and E_new are just the energies before and after you update the kinetic energy, at any given step. Delta E is the difference between the energy of the system and what you want it to be. You can't change the potential energy without moving particles, but you can change the kinetic energy, so in practice it is the difference between the kinetic energy and what you want it to be. That changes at every step.
I'm not sure I follow your questions so may not have answered them well.
Why is LJPot never used?
The code just simulates the dynamics, which only depends on equations of motion and thus derivatives of the potential, not the actual potential itself. If you want to calculate the internal energy then LJPot itself is required. The internal energy (and thus the actual potential) is required for a bunch of stuff you might want to know about the system, like the heat capacity or if you want to do thermodynamic integration or lots of other things. I wanted to include the function/subroutine but since the code just evolves the system in time, it's not useful to calculate LJPot for anything in the scope of the video.
Where does reflectBC() go?
It's easiest to put it right after the particle position is updated. In the example code, I put it at the end of the 'update' routine but it can also be after, depending how you implement it.