15
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,022
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.2k
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,523
6
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 ...
- 677
6
votes
Fortran code for Ewald summation
Notice that the summation in $U^{(r)}$ is incorrect. You want to sum over all the copies of the atoms in the lattice of periodic boxes, not just those whose indices satisfy $i > j$. In the ...
- 3,994
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 ...
- 2,627
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
Accepted
Understanding leapfrog integration algorithm
In the second code the full time stepping is given by three lines in main()
...
- 465
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 ...
- 3,994
4
votes
Creating dense random configuration in for molecular dynamics
As I am assuming by the word random you mean "liquid or gas like" is there any reason you can't use your MD program itself?
Pick the number of particles and the volume of the unit cell so
as to get ...
- 491
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$. ...
- 677
4
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 ...
- 5,519
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.2k
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
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,022
3
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,275
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+...
- 401
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 ...
- 5,519
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 ...
- 2,627
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
periodic boundary conditions for triclinic box
Well there's always the general approach. Let's say we have a particle at a position $\bf r$ and the lattice vectors, i.e. the vectors that forming the sides of the simulation cell, are $\bf a_i$. We ...
- 491
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 ...
- 2,627
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
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.1k
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 ...
- 3,994
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.
...
- 3,994
Only top scored, non community-wiki answers of a minimum length are eligible