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 ...
user avatar
  • 1,819
11 votes
Accepted

Removing non-determinism from molecular dynamics code

Both the standard cluster and custom supercomputer (Anton) versions of molecular dynamics at D. E. Shaw Research are both deterministic and parallel invariant. That is, a test run on a single core ...
user avatar
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 ...
user avatar
6 votes
Accepted

Generating random numbers for long molecular dynamics simulations

I suspect the reason why generating the random numbers on the fly is slower for you is due to the rather large state of the Mersenne Twister. Switching to something like the PCG or XorShift+ random ...
user avatar
  • 2,463
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 ...
user avatar
  • 3,483
6 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 ...
user avatar
5 votes
Accepted

How do I generate Maxwell-Boltzmann variates using a uniform distribution random number generator?

The initial velocities are drawn from a Gaussian distribution with variance $$\sigma_i^2=\frac{k_{\textrm{B}}T}{m_i},$$ where $k_{\textrm{B}}$ denotes Boltzmann's constant, $T$ is the temperature and $...
user avatar
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-...
user avatar
  • 207
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 ...
user avatar
  • 436
4 votes

How can I make velocity verlet algorithm more stable?

You can try using an adaptive timestep, where dt is adjusted at each step if the velocities get too large. However, the Verlet method is symplectic, so if the total ...
user avatar
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 ...
user avatar
4 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 FORTRAN. ...
user avatar
  • 657
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$. ...
user avatar
  • 657
3 votes
Accepted

PQR files make from pdb2pqr are different to GROMACS

The Van der Walls radius - last column in the output - is calculated from the force field. Probably pdb2pqr and editconf uses different force fields, hence different radius. I don't use pdb2pqr, but ...
user avatar
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 ...
user avatar
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{...
user avatar
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 ...
user avatar
  • 323
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 ...
user avatar
  • 1,819
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 ...
user avatar
  • 1,285
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+...
user avatar
  • 391
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 ...
user avatar
  • 3,511
2 votes

Removing non-determinism from molecular dynamics code

I believe that it close to impossible to get binary reproducibility (or "determinsm") in a stochastic MD simulation. Such simulations are fundamentally chaotic - even when the least significant bit of ...
user avatar
  • 411
2 votes

Molecular Dynamics: Diffusion with PBC

You need to unfold the positions after the simulation is done. I use the following subroutine: ...
user avatar
2 votes

Unexpected behavior, values tend to converge instead of fluctuate. (MD)

Your algorithm is a first-order Euler and not the Verlet algorithm. What is the actual box size when you run the code? If you made a mistake and your particles end up being well separated (i.e. ...
user avatar
2 votes

Is Langevin thermostat/equation correct when trying to model time-dependent behaviour of a molecule?

If $H = H(q, p)$ is the Hamiltonian of your mechanical system, $T$ is the temperature and $\gamma$ is the friction, then Langevin dynamics (LD), \begin{equation*} \left\{ \begin{aligned} \...
user avatar
2 votes

periodic boundary conditions for triclinic box

I have never implemented periodicity for this symmetry class, although I have done it for hexagonal crystals before. Let us start with a (2D) box like the one depicted in the image In this case, we ...
user avatar
  • 7,902
2 votes

Proper handling of 1-2 1-3 1-4 bonded neighbors in long range electrostatic solver

You have to include the bonded interactions in the list of so-called masked pairs. If you consult [1], this is the first term in equation (2.5). That will subtract the contribution from the long-...
user avatar
2 votes

energy drift in molecular dynamics

You've measured the relative deviation from the initial energy, whereas "energy drift" in MD is only usefully measured as a rate of absolute deviation, and only after the system has equilibrated. (...
user avatar
  • 287

Only top scored, non community-wiki answers of a minimum length are eligible