8
$\begingroup$

I've been looking through the other answers, and I haven't found a good answer yet. I'll describe the simulation I am running, and then the problems. The program simulates particles undergoing Brownian dynamic motion (thermal motion) alternating with stochastic processes (kinetic Monte-Carlo). When I move to OpenMP code, I gain performance (finally), but I lose the determinism in the system. Therefore, I have two problems.

  1. I can only repeat the same 'experiment' twice using a single threaded version. Otherwise, the order in which the threads operate on the particles destroys any determinism.
  2. Profiling how much faster the OpenMP code runs as a function of threads, or versus the serial version, requires running 10x sims at every point since I test for a particular ending configuration, and the speed of the program is very sensitive to the configuration of the particles. Basically, the closer more particles are, the longer it takes to run, so different ending configurations can take drastically different amounts of time.

I know that this is a normal problem found in MD, but I am seeing if there are any ways around it rather than running the single threaded version. If that is the case, then so be it, I was just hoping for something else mostly so that I could optimize the OpenMP code and environment variables.

Improve this question
$\endgroup$

4 Answers 4

Reset to default
12
$\begingroup$

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

  1. Integer summation: Although each force term is computed in floating point, the total force on each atom is summed in fixed point. As a result, the sums are both commutative and associative, and therefore trivial to parallelize without losing determinism. All state (positions and momenta) are maintained in fixed point as well.

  2. Any required random numbers are generated using counter mode techniques, as described in http://www.thesalmons.org/john/random123/papers/random123sc11.pdf. A counter mode generator amounts to encrypting the integer sequence 0, 1, 2, $\ldots$ using a sufficiently strong cypher. See the paper for details.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I will take a look at this when I have a chance. Right now I am trying to get results ASAP, so that takes precedence over figuring this out, even if this would save me tons of time. $\endgroup$
    – NuclearAlchemist
    Commented May 30, 2014 at 16:02
  • 1
    $\begingroup$ This sounds very interesting. Can you comment whether there is any performance loss associated to these techniques? When doing BD, you needs lots of random numbers (basically three per particle per step, if I'm not mistaken), so performance is of considerable importance. $\endgroup$
    – olenz
    Commented Jun 2, 2014 at 8:40
  • 2
    $\begingroup$ Several of the counter mode generators are faster than Mersenne twister. It also occasionally helps that they have much less state: 256 bits vs. lots of Mersenne twister. The period is correspondingly lower, but this isn't a problem unless you have a machine that can perform $O(2^{128})$ operations. All the details are in the linked paper. $\endgroup$
    – Geoffrey Irving
    Commented Jun 2, 2014 at 16:55
  • 2
    $\begingroup$ The performance implications of integer summation are small as long as computing a given force takes more than a couple flops. The only extra operation is a float to int conversion per computed force. On BlueGene this is a bit worse due to the lack of integer vector instructions. Summing in fixed point does require you to work out the right scale, but also gives you additional precision (no bits wasted on exponents). $\endgroup$
    – Geoffrey Irving
    Commented Jun 2, 2014 at 16:59
2
$\begingroup$

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 any variable (e.g. any particle position) in two runs of a simulation is different, the trajectories will diverge after some time.

A typical source of such an error is the summation order of the forces in the main MD loop: when you sum up the forces in a different order, this might already lead to differences in the least significant bits of the sum. To get full determinism, you would have to make sure that the forces are always summed up in the same order.

Unfortunately, it is pretty much impossible to ensure the order of summation. When the compiler optimizes code, one of the first things it does is to change the order. This means that when you compile the same code with different compilers (or maybe even different versions of the same compiler, or different optimization levels), you can not expect the results to be identical on a binary level. Furthermore, some modern architectures do instruction reordering even on the hardware level, or they provide faster instructions where the least significant bit is undetermined. So even when you execute the same binary on the same machine, it is not necessarily guaranteed that both runs create the same result. This problem is even enhanced when you parallelize the code in any efficient way, as this will for sure modify the summation order.

Note that this means that even when you manage to write a program that creates binary reproducible results on your machine, you cannot expect that anybody else will be able to reproduce your results unless he has the same version of your compiler, the same version of your OS and the same hardware.

Consequently, I have to say that for MD, we simply cannot gain binary reproducibility. On the other hand, the results from an MD usually are statistical mean values anyway. So, what we should aim for when doing MD simulations is statistical reproducibility rather than binary reproducibility.

Improve this answer
$\endgroup$
2
  • $\begingroup$ See my DESRES answer for why "impossible" clearly isn't true. $\endgroup$
    – Geoffrey Irving
    Commented May 29, 2014 at 19:45
  • $\begingroup$ I stated "close to impossible". The techniques you describe in your posting are pretty involved to implement, and may cost considerable performance. I'm keen to hear on the performance aspect. Still, I wonder: is binary determinism worth the hassle, as we are after statistical results, anyway? $\endgroup$
    – olenz
    Commented Jun 2, 2014 at 8:46
1
$\begingroup$

Do you know where exactly the non-determinism is stemming from? I don't do MD, but I can think of two scenarios:

  • The particle updates (time stepping) are using the most recently available particle data, so another OpenMP thread may have already updated particle that wasn't updated in the serial version. Fixing this would require all threads to only update based on the previous time step's values. (I would think this would be required anyway to maintain time accuracy though.)

  • You mentioned "stochastic processes." How are you seeding the random number generator used here? If you force all threads to use the same seed, that should remove any non-determinism from the RNG.

Improve this answer
$\endgroup$
4
  • $\begingroup$ The problem is the latter. I give each thread it's own RNG with a known seed, but the mapping from thread to particle isn't always the same, since it is determined by the openmp runtime. So we could get different random numbers, and therefore forces, for the same particle depending on which thread operates on it first. $\endgroup$
    – NuclearAlchemist
    Commented May 28, 2014 at 16:42
  • 1
    $\begingroup$ Is this just done with a simple OpenMP pragma on a loop? Specifying STATIC scheduling (with # of particles an integer multiple of number of procs) should make it deterministic then for a fixed number of threads. Alternatively, you could serialize the RNG part of the code: have only thread 0 compute all the random numbers from a fixed seed and store them in a vector, and then execute the parallel loops. If you're doing significant work on the updates and with a relatively small number of cores, it shouldn't impact performance too much. $\endgroup$
    – Aurelius
    Commented May 28, 2014 at 17:45
  • $\begingroup$ I was going to avoid the static clause since it significantly impacts performance for some of the configurations I wind up in. I don't know how much of a problem serializing the RNG would be, that is probably the way to go, since it should be somewhere less than 3k random numbers per iteration. $\endgroup$
    – NuclearAlchemist
    Commented May 28, 2014 at 19:12
  • $\begingroup$ There is never a need to serialize random number generation. Use thesalmons.org/john/random123/papers/random123sc11.pdf. $\endgroup$
    – Geoffrey Irving
    Commented May 29, 2014 at 19:37
1
$\begingroup$

You need to make sure that every particle update uses the same random number, whether you run on one of many threads. As an extreme case, consider that you have one RNG associated with every particle. Then, when you do the $n$th update for particle $m$, you will get the $n$th random number from RNG $m$, regardless of the thread you do this update on.

This may be a bit extreme, but you could have one RNG per, say, 1000 particles and each thread can only update blocks of 1000 particles at once, using the associated RNG. Again, this way it makes no difference which thread makes the update since it draws random numbers deterministically.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Just use one RNG per particle using counter mode generators. $\endgroup$
    – Geoffrey Irving
    Commented May 30, 2014 at 3:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.