16
votes
Accepted
Why is leapfrog integration symplectic and RK4 not, if the latter is more accurate?
TL;DR: It depends on what kind of accuracy you need.
Energy conservation does not automatically equal accuracy. Suppose, you want to simulate the solar system, and you are using a solver that – to ...
- 2,127
7
votes
Accepted
MD Simulation: Reference for the Neighbor's List Method
I'd recommend "The Art of Molecular Dynamics Simulation" by D. C. Rapaport. The code samples are written in C. I'm not a huge fan of the programming style of the book, but at least it's not ...
- 687
7
votes
Accepted
The velocity Verlet method and variable time steps
The second formula is just the velocity verlet, and it's correct but if you adapt time steps then it's not symplectic. In a separate answer I describe in quite detail that symplecticness is a global ...
- 12.4k
6
votes
What good are hard-sphere event-driven molecular dynamics simulations in the face of chaos?
Yes, it is possible to show that the statistical behavior of the approximate system will reach that of the "exact" system. (This is true even though hard-sphere dynamics do not accurately describe ...
- 3,513
6
votes
Accepted
Understanding leapfrog integration algorithm
In the second code the full time stepping is given by three lines in main()
...
- 485
6
votes
Starting configuration for Molecular Dynamics
Usually one needs to employ periodic boundary conditions (at least in the horizontal directions). Any atoms which fly outside of the box will be mapped to the opposite side. This also has to be ...
- 3,005
5
votes
Fortran code for Ewald summation
The National Institute of Standards and Technology (NIST) proposes different inputs and outputs to test your implementation. Have a look there:
https://www.nist.gov/mml/csd/chemical-informatics-...
- 207
5
votes
Understanding leapfrog integration algorithm
The first code is wrong if you have multiple particles that interact. Due to the structure of the loop, the interaction forces with particles at the start of the list are computed with the updated ...
- 6,149
4
votes
How to integrate numerically Nosé Hoover equation?
You can find explicit pseudo-code for Nosé-Hoover Chains in Section 4 of:
Martyna, G. J., Tuckerman, M. E., Tobias, D. J., & Klein, M. L. (1996). Explicit reversible integrators for extended ...
- 4,004
4
votes
Calculation of Mean Square Displacement for Brownian dynamics system with Lennard Jones interactions in python3
This will not an answer to your problem, more an excessive comment and few things you might consider, when writing such code even for self educational purpose.
Constants
You asked whether your ...
- 1,305
4
votes
Calculate forces on atoms from potential energy of system and position of atoms
Generally speaking, the force acting on particle $i$ is just $-\nabla E(\vec{r}_i)$ (note the minus sign), where $\nabla$ is the gradient operator and $\vec{r}_i$ is the position of particle $i$. ...
- 687
3
votes
What good are hard-sphere event-driven molecular dynamics simulations in the face of chaos?
Something that is at least very close to what your are asking for is the Shadowing Lemma, which roughly states that near any numerical solution of a dynamical system, you will find a real solution.
I ...
- 2,127
3
votes
Accepted
How to compute forces in multi-particle MD
I will only be developing the computation of forces on a given atom (part 1 of BalazsToth answer), for this point is actually a bit tricky and I think it deserves to be a bit more explained. You can ...
- 323
3
votes
How to compute forces in multi-particle MD
I have the following remarks:
1. Force calculation
The equation of motion of the $i^\text{th}$ atom with mass $m_i$ reads as
$\textbf{a}_i=\frac{d\textbf{v}_i}{dt}=\frac{1}{m_i}\cdot f(\vert\textbf{...
- 430
3
votes
How to integrate numerically Nosé Hoover equation?
tl;dr
The Nose-Hoover equations are normally defined by their Hamiltonian. The two partial derivatives of the Hamiltonian define a partitioned set of ODEs. From there the partitioned ODEs are either ...
- 12.4k
3
votes
Accepted
Which are the right configurations in the Markov chain of a Hamiltonian Monte Carlo algorithm?
The algorithm for perfoming a single HMC step is as follows:
Input: Some initial configuration $\vec{y}_i$ and momentum $\vec{p}_i$.
Output: Next configuration $\vec{y}_{i+1}$ and momentum $\vec{p}_{i+...
- 451
3
votes
Accepted
How error scales with numerical precision in molecular dynamics?
In general in a correctly implemented fixed-step ODE solver method you have 3 sources for numerical errors:
the theoretical method truncation error,
the floating point error from evaluating the ODE ...
- 6,149
3
votes
How can I write a unit test for this routine?
Divide and conquer
refactor your code, so that you have smaller methods with well defined purposes. Give these methods short descriptive names of what they are doing. You make the code more testable ...
- 3,005
3
votes
Accepted
Computing LJ force from LJ potential, or not?
Which one is computationally more accutate [sic] and why?
The first one is more accurate. It uses the exact derivative of the LJ potential.
By contrast, the second one not only uses a numerical ...
- 198
3
votes
Starting configuration for Molecular Dynamics
The biggest issue in your subsequent simulation won't be the periodicity of the box -- it will be trying to equilibrate your system down from what looks like a tangle of lipids into an orderly ...
- 131
2
votes
Difficulty with possibly the simplest MD simulation
Two thoughts
A) A certain degree of energy dissipation in MD simulations is considered normal and can be counterbalanced with so called thermostat algorithms e.g. this. In your particular case I can ...
- 3,005
2
votes
Oil/Water interface simulation using GROMACS
A simple solution to begin with would be to create two systems: one with decane in a box of width 1.295nm and another one with water in a box of width 2nm (depth and height should be as in the picture ...
- 4,004
2
votes
Why are Hamiltonian dynamics used in MCMC?
I'm a little late to reply but I found this review to be really informative. One particularly nice feature of Hamiltonian flows is that they preserve volume in phase space. If we take an infinitesimal ...
- 10.4k
2
votes
Fortran code for Ewald summation
I don't think your formula is correct for the force. You can rewrite your energy expression in an easier form, and it makes taking the derivative easier as well.
The way you have it is
\begin{equation}...
- 121
2
votes
Accepted
Can LINCS algorithm be used for colliding molecules?
Constrained simulation methods such as LINCS are there to keep certain degrees of freedom fixed throughout the simulation, so you can use them to make each of your molecules behave like a rigid body.
...
- 4,004
2
votes
Time Reversibility of Velocity Verlet Algorithm
I can share my perspective on this topic and can try to prove that the Leapfrog (velocity Verlet) is time-reversible, according to an appropriate definition of time reversibility (very nice property ...
- 621
2
votes
Time Reversibility of Velocity Verlet Algorithm
I think for time reversibility we need to see if we can predict velocity and position at time $t$ if the values at $t+\Delta t$ is known.
Verlet Algorithm
The integrator equation is given as:
$r_{t+\...
- 21
2
votes
Who can explain the minimum image convention in molecular dynamic simulations?
The question as stated does make sense in the context of molecular dynamics simulations but is very broad.
I recommend that you should check the manual of the MD package your are using about the image ...
- 1,305
Only top scored, non community-wiki answers of a minimum length are eligible