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. ...


12

The conjugate gradient method is good for finding the minimum of a strictly convex functional. This is typical when you reformulate a nonlinear elliptic PDE as an optimization problem. If you want to learn about it, I recommend you read about the CG method for linear systems first, for which An Introduction to the Conjugate Gradient Method Without the ...


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

To ask for a progression year-on-year or even more fine-grained is asking for a bit too much detail, I think. Available information is more coarse-grained. This paper reports a nine orders of magnitude progress in simulated timescales over three decades of molecular dynamics research on proteins in water.


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

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

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

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

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

Since you are specifically looking for a library, you might be interested in OpenMM. Two popular visualization packages are PyMol and VMD.


3

This link points to a thesis that contains an otherwise unpublished chapter (chapter 6) on the simulation of a droplet of water on graphite using a force field that includes a linear response model for the dipole moment. A lot of the approach is based on the following publication: Abdolnour Toukmaji, Celeste Sagui, John Board, and Tom Darden. Efficient ...


3

You can also look at CMAES. It essentially boils down to CG for convex functions, yet represents global and robust optimizer for other types of functions (including non-convex functions with multiple minima). I have not seen, however, its application to anything larger than a couple of hundreds unknowns. Note also that CG can be applied in combination with ...


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 ...


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