17

Double precision is fairly common on newer GPUs. For instance I own a NVIDIA GTX560 Ti (fairly low end when it comes to computing) that has no issue running ViennaCL in double precision. From here (section 4) it appears all NVIDIA cards from GTX4xx onward support double precision natively. I would guess that the GROMACS information is simply outdated.


13

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


8

Since the size of each type of atom is fixed, for a given level of accuracy the asymptotic cost is dominated by far field electrostatic interactions. These are $O(n)$ using multigrid and $O(n \log n)$ using FFTs. Thus the optimal complexity is $O(n)$ per time step as a function of the number of atoms simulated. The time step is also asymptotically ...


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

This isn't a multigrid issue, it's a problem formulation issue. Consider the symmetric system $A x = b$ and suppose that $A e = 0$ for some nonzero vector $e$ (the constant when $A$ is the Laplacian with periodic or Neumann boundary conditions). This system is singular and has no solution if you choose $b$ such that $e^T b \ne 0$. Indeed, the exact periodic ...


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


5

The area you are asking about is known formally as ab initio molecular dynamics. Marx and Hütter a few years ago wrote a comprehensive monograph on recent updates in the method, and can be recommended for providing a useful summary of the (almost) current state-of-the-art in the field. However, one thing that should definitely be considered before doing any ...


5

The key is to take the differences $\Delta x$, $\Delta y$, and $\Delta z$ separately before beginning. Given the edge vectors ${\bf a}$, ${\bf b}$, and ${\bf c}$ that define a unit cell that has one corner at the origin, then your tilt factors are $b_1$, $c_1$, and $c_2$, where $n_m$ defines the $m^{th}$ component of vector ${\bf n}$. Now, for each ...


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


5

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


4

Did you check the centre-of-mass velocity of your particles, i.e. did you set it to zero at the beginning of your simulation? The centre-of-mass velocity should be preserved throughout the simulation, and should be zero to get correct temperatures (see Flying Ice Cube).


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

This is a re-formulation of my comments above as an answer. You need to enforce the boundary conditions when computing dx and dy to account for particles that interact over the periodic boundary. Imagine two particles at opposite ends of the domain: In your current setup they won't interact, but as soon as one of them moves slightly and crosses the boundary,...


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

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

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

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

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


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