14

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 use an extreme example – just rotates the entire system by some angle every second. These solutions obviously conserve energy, but they are blatantly incorrect. ...


11

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 generates the same bits as a massively parallel run. The techniques include Integer summation: Although each force term is computed in floating point, the total ...


7

Have you had a look at VMD? I used it ages ago to produce movies from simulation snapshots. Way back then, it could read a sequence of PDB files, render them (or generate POV-Ray scripts to raytrace them), and store them as individual images. I then used mencoder to generate MPEG-4 files out of the stills. Those were the days. I haven't used VMD since, but ...


6

Any numerical differences between A and B will become exponentially larger with time (i.e. the Lyapunov instability, as discussed in Frenkel and Smit). Even a small difference due to basis set size could result in dramatic differences in the trajectory over time. So I'm not sure that a comparison between individual trajectories will be meaningful. It may be ...


6

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 number generator would have several advantages for you: Higher quality of randomness (Mersenne Twister fails several tests for randomness) Smaller state, so ...


6

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 molecular systems!) The basic premise underlying molecular dynamics is the ergodic theorem, which states that, in the limit of long times, the time average of a ...


6

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 property of the integration, it's not a stepwise property. Because of this, local error estimates and local changes of $\Delta t$ do not necessarily preserve the ...


5

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 $m_i$ is the mass of the $i^{\textrm{th}}$ particle. Thus, the problem boils down to generate random numbers from a gaussian distribution using uniformly ...


4

By "Domain decomposition is a better choice only when linear system size considerably exceeds the range of interaction, which is seldom the case in molecular dynamics" the authors of that (very old) GROMACS paper mean that if the spatial size of the neighbour list is of the order of 1 nm, and the simulation cell is only several nanometers, then the overhead ...


4

It's possible that all the memory being used by your code is on one socket and that up to 6 cores all the tasks are running on that socket. When you get to 7+ sockets, then there are transfers between sockets to get at the memory. You may need to investigate memory affinity options for your threads. The default policy in Linux is first-touch (I think), so if ...


4

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 original box, of course you want to avoid self interactions $i = j$ but just in that one box. In other words, do not discard the electrostatic interactions between ...


4

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 the density you want Set the particles up on a regular lattice such as face centred or body centred cubic with some vacancies if required by the number of ...


4

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 systems dynamics. Molecular Physics, 87(5), 1117–1157. Retrieved from http://www.tandfonline.com/doi/abs/10.1080/00268979600100761 Here is a short example of the ...


4

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. Having said that, my advise would be to take any book where neighbour lists are explained (for instance the Frenkel & Smit book, which I guess it's what you ...


4

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$. However, note that this expression assumes that we are dealing with conservative forces. But if that's the case, the energy is conserved by definition, which means ...


3

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 it seems (1) AMBER99 is the default force field, though CHARMM, PARSE and TYL06 are supported. Gromacs' editconf reads force field from topology file generated ...


3

Periodic boundary conditions (PBC) are a common technique to get closer to the thermodynamic limit while keeping the computation feasible. You could set up an experiment simulating a box of argon with and without PBC and study how the radial distribution function g(r) changes in each case, you should observe a difference in g(r) for large values of r. ...


3

Pedro's answer would have been mine. If you want to manage balls and sticks on your own, though, VTK is a modern library for such things.


3

A common test consists of computing some ensemble averages that you can check analytically and repeating the experiment a number of times while decreasing the time step length. You should observe that, when you plot the time step length against the error in log-log scale, the data points can be easily fitted to a line with slope equal to the order of the ...


3

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 energy is not roughly conserved then something is wrong. If two particles get too close, no reasonable timestep will make the method stable, and accurately computing the forces between particles ...


3

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-research-group/spce-water-reference-calculations-10%C3%A5-cutoff See paragraph 6. Reference Calculations of Intermolecular Energy for SPC/E that includes the usual ...


3

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 solved by an integrator for 1st order ODEs or some symplectic method. Details shown using DifferentialEquations.jl As a rather thorough example, let's look at ...


3

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{r}_j-\textbf{r}_i\vert)\textbf{n}_{ji}$, where the interparticle force $f(\vert\textbf{r}_j-\textbf{r}_i\vert)$ depends only on the distance between particles ...


3

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 find useful remarks in any book of molecular simulation my personal reference being "Computational Simulation of Liquids" from Allen and Tildesley and the "...


3

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 am no expert on this, so I can only provide this a starting point for further research. In particular I do not know whether this translates to non-smooth ...


3

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 dimensionless code uses correct constants / normalization. So let's have a look. For this question it helps to write down the equations, you want to solve. Here I get ...


3

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+1}$ Draw a random momentum $\vec{p}_*$ from a Gaussian distribution. Numerically solve Hamilton's equations of motion for some time (i.e., perform some ...


3

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 function and composing the method step and the error from accumulating the single updates of size $O(h)$ to the integration result of size $O(1)$. So the total ...


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