I think a few clarifications are necessary here:
1. The linked-cell code can mimic the naive O(N^2) algorithm if you make
the following changes to md_main(...) (In main.jl)
n_cell = [ 1 for dim = 1:DIM ]
l_cell = [ 2r_cutoff for dim = 1:DIM ]
With this modification the code is pretty much like a direct sum algorithm.
In fact, this would be slightly faster than using the 2 x 2 x 2 grid for
the the current benchmark problem since it is cheaper here (just 125
particles!) to do an explicit search than to do all the bookkeeping about
which particle is in which cell. The real advantage of the linked-cell
method would show up as you increase the domain size and the number of
particles.
2. The linked-cell method is perfectly valid as a serial algorithm (though
the size problems you can tackle would accordingly be restricted). The
basic idea behind the method is that if the forces between the particles
are short-ranged, the calculation of the inter-particle forces is best
restricted to the set of all particles within a cut-off radius of a given
particle. To accomplish this the physical domain is represented as a
disjoint union smaller cells whose physical dimensions are of the order of
the cut-off radius for the inter-particle potential - this ensures that
each cell handshakes only with its nearest neighbor cells (in the geometric
sense). An elementary description of the method can be found
here: http://en.wikipedia.org/wiki/Cell_lists. Note that this method would
not work in the presence of long-range forces. Regarding parallelization,
it is not necessary that each cell lives in a different processor (and
moreover, that is likely to be quite inefficient). Parallelizing this
method is essentially about exploiting the geometric layout of the cells
and mapping groups of them into various processors so as to minimize the
talking time between the cells on different processors. I'm currently
working on this and will have something interesting to share soon,
hopefully.
3. I had slightly modified the initial positions to make sure that all the
particles are inside [0,10]x[0,10]x[0,10] initially (you can find the
modified He-LJ.xyz file in the link I gave earlier). Running the code with
the original code .xyz file would likely discard the points (just 4 I
think) outside the domain, and this could be one reason why it runs faster
on your machine. Currently, periodic boundary conditions are enforced only
during computation, but not during initialization. Or maybe my laptop has
grown old :)
Cheers,
Amuthan
Post by LuthafThat's already way better than the simple and naive version ! On my
computer, this gives me something 4x slower than serial LAMMPS code.
I'll try to integrate this into the package. It will need changes in the
data structures, but this is worth it !
I'll have to think a bit about how I can organise the data for
parallelism. The basic idea is to have every cell on a different processor,
isn't it ?
Bests,
Guillaume
Hi Guillaume: I've added the serial version of the linked-cell MD code
here: https://gist.github.com/Amuthan
The code is pretty basic and needless to say, a lot more functionality
needs to be added (though I should admit that I'm more interested in
extending this code to a parallel version than to develop it as a
full-fledged MD package). I've made a few changes to make it work for your
benchmark; I didn't spend too much time optimizing it, but the current
version is just 7 times slower than the LAMMPS code. I guess someone with a
better understanding of Julia can optimize it further.
Let me know if you have any comments or questions or suggestions for
improvement.
Cheers,
Amuthan
Post by LuthafHi Amuthan,
The current code uses the direct sum because this was the easiest
algorithm to implement, and your understanding is totally correct.
Adding other algorithms is a piece of cake: it is as easy as subtyping
`BaseForceComputer` (I could not came with a better name), and providing a
`call` method.
So I am planning to add at least (serial) Verlet neighbor lists for the
forces computations, and maybe some parallel code.
But I don't know a lot about parallel methods in MD, and I still have a
lot to learn about it. It would definitively be very nice to have it !
I think that only a few changes are needed to use julia parallel
framework in the package. All the data are stored in objects of type
Array3D, which could be modified to use DArray instead.
I am open to change the internal data structure if this can improve the
performances.
Is your neighbor list code available somewhere ? I am interested in
reading it.
Cheers,
Guillaume
This looks good -- thanks for the package! I have a basic question
regarding the force computation: I had a quick look at the code and it
appears to me that the current code uses an O(N^2) algorithm (direct sum)
for this. Is my understanding correct? Software like LAMMPS use very
efficient algorithms for short and long range potentials and this might
possibly account for the fact that the benchmark is 80x slower.
I had developed a serial MD code using neighbor lists last year and I'm
currently in the process of parallelizing it to study the suitability of
Julia for developing (reasonably) parallel research codes. Is
parallelization something you're looking into? I would be very interested
in hearing about this. Thanks!
Cheers,
Amuthan
Post by LuthafHi julians !
I am very happy to announce the release of Jumos, a Julia package for
molecular simulations.
You can find the code here: https://github.com/Luthaf/Jumos.jl, and the
http://jumos.readthedocs.org/en/latest/.
This package aims at being as extensible as possible. Every algorithm is
a type that can be subtyped and changed at runtime, during the simulation.
This makes it easy to write new algorithms, experiment with them, and
combine algorithms in all the fancy ways you can imagine.
Presently, only molecular dynamic is implemented, but I am planning to
add energy minimisation, and trajectory replay from files (maybe
monte-carlo simulations if I find some time). A handful of tools are
- Velocity-Verlet and Verlet integrator;
- Lennard-Jones and Harmonic potential;
- User-defined potentials;
- RDF and density profile analysis;
- Berendsen thermostat;
- Velocity-rescale thermostat;
- Temperature and energy computations.
And a lot more are planned: https://github.com/Luthaf/Jumos.jl/issues.
If you want to help, contributions are very welcome ! Just pick up an
issue or create a new one as a feature request.
The first benchmark (against LAMMPS) gives me a speed 80 times slower
than LAMMPS. This is already great for a dynamic language, but there is
still room for improvement !
Any feedback on the implementation is also very welcome.
Regards
Guillaume
PS: What does the ANN in the title mean ? It seems to be everywhere,
but I don't have any clue about its meaning