Where should we put our particles

[Luthaf]

Currently, the particles are stored contiguously in a vector, and the molecules have an usize field pointing to the first particle of the molecule in this vector.

struct System {
    particles: Vec<Particle>
    molecules: Vec<Molecule>,
}

struct Molecule {
    first: usize
}

Pro of this design

• good cache locality for particles
• easy to iterate over all particles in the system. A double loop over pairs is:

for (i, pi) in system.iter().enumerate() {
    for pj in system.iter().skip(i + 1) {
        // Do work
    }
}

Cons of this design

• hard to iterate over particles in a molecule, and corresponding double loops are confusing.

for molecule in system.molecules() {
    for i in &molecule {
        for j in &molecule {
            if j <= i {continue;}
            let pi = system[i];
            let pj = system[j];
            // Do work
        }
    }
}

Making the particles owned by molecules is not a solution, because it make it very hard to iterate over pairs in the whole system, and add an extra indirection layer.

A possible solution, with a POC implementation here would be to store a pointer to the first atom in the molecule, together with the molecule size. This would make it easier to work with molecular code, but add some unsafe usage.

In particular, the pointer of all molecules must be updated when adding a new particle to the system, because the vector may allocate additional memory.

So is it worth it to add complexity in the library code in order to simplify the user code?

[Luthaf]

Another solution would be to have two separated types: an "owning" molecule and a reference molecule.

The owning molecule (called Molecule) would contain a Vec<Particle>, and can be manipulated without a System. A reference molecule (MoleculeRef<'a>) would contain a slice of particles &'a [Particle] borrowed from a system, and thus can not be manipulated without a System. The owned molecule would then implement Deref<Target=MoleculeRef>, and most of the utility functions for molecules would be implemented on MoleculeRef.

This allow users to manipulate particles through the molecules, without loosing cache locality.

[g-bauer]

I am not sure if I got that. Could you elaborate about cache locality? I can see however that the current implementation is a bit uncomfortable to use. Are there problems with something like this (using horrible naming):

// method of System

/// Returns the particles of the molecule with index `id`
#[inline] 
pub fn particles_in_molecule(&self, id: usize) -> &[Particle] {
    let range = self.molecules[id].first() .. self.molecules[id].size();
    &self.particles[range] 
}

or using a &Molecule as input

// method of System

/// Returns the particles of the molecule `molecule`
#[inline] 
pub fn particles_in_molecule(&self, molecule: &Molecule) -> &[Particle] {
    let range = molecule.first() .. molecule.size();
    &self.particles[range] 
}

[Luthaf]

I had written a long comment on this issue some days ago, but forgot to post it. Here it is again !

---

The owned molecule would then implement Deref<Target=MoleculeRef>

This can not be done, because Deref targets must be references. They can not be new struct created out of thin air.

But this gave me another idea: having a ParticleGroup trait, that we can use to implement properties (see #20) and abstract the origin of the particles.

trait ParticleGroupe {
     fn particles(&self) -> &[Particle];
     fn bonds(&self) -> &[Bonds};
     // maybe some other methods dealing with topology
}

// A molecule owning its atoms
struct Molecule {
    particles: Vec<Particle>,
    bonds: // ...
}

impl ParticleGroup for Molecule {
    // ...
}

// The current system
struct System {
    particles: Vec<Particle>,
    molecules: // ...

}

impl ParticleGroup for System {
    // ...
}

// A molecule with a reference to its atoms
struct MoleculeRef<'a> {
    particles: &'a [Particle],
    bonds: // ...

}

impl<'a> ParticleGroup for MoleculeRef<'a> {
    // ...
}

impl System {
    fn molecule<'a>(&'a self, i: usize) -> MoleculeRef<'a> {
        // ...
    }
}

Then, we can use this trait in Compute algorithms:

struct Energy;

impl Compute for Energy {
    type Output = f64;
    fn compute<P: ParticleGroup>(&self, particles: P) -> f64 {
        ...
    }
}


struct CenterOfMass;

impl Compute for CenterOfMass {
    type Output = Vector3D;
    fn compute<P: ParticleGroup>(&self, particles: P) -> Vector3D {
       //  ...
    }
}

[Luthaf]

Could you elaborate about cache locality?

CPU have multiple level of memory available: slow memory with big capacity (the RAM), and fast memory with small capacity (the CPU L1, L2 and L3 caches). Accessing L1 cache is 4 CPU cycles, and accessing RAM is ~800 CPU cycles. But the RAM can be as big a TB of data, and the CPU L1 cache is 32 KB on Intel Haswell. The CPU also pre-fetch some RAM in the cache when accessing memory, to improve performance of further memory access.

For more information concerning CPU cache and how to use them, see:

• http://stackoverflow.com/a/16699282/4692076
• https://people.freebsd.org/~lstewart/articles/cpumemory.pdf

Here the idea is that we are often accessing our Particle sequentially. So it is better to have them stored in a Vec that will be loaded in cache all at once. By position, having Vec<Molecule> where Molecule contains a Vec<Particle> means that the particles data must be fetched from RAM (as different Vec will point to different RAM area) for every molecule. Also, as the cache will already contain the molecules data, fetching particle may replace this data, and require to fetch the molecules again at the next iteration step.

I hope I've been clear, if not I'll try again!

Are there problems with something like this (using horrible naming):

This could work, but will only give access to the particles. My proposition above is a bit more abstract, and gives access to both particles and molecular connectivity (bonds, angles, ...)

The other thing I would like to improve here is working with molecules out of a System. For example, the molecule_type algorithm, and the read_molecule function should use a single struct, not have a Molecule and some &[Particle] separated.

using a &Molecule as input

This is a bit dangerous, as there is no way to check if the molecule does come from this system or not.

[g-bauer]

Thank you! So in essence, in the most performance critical parts you want to have all information that you need very close in memory?

Your idea with ParticleGroup is to have an interface where every implementation specifies the way to access the particles (and corresponding interactions)? I think that is a good idea, since implementations in system start to get a bit messy (read: there are a lot of implementations).

// A molecule owning its atoms
struct Molecule {
    particles: Vec<Particle>,
    bonds: // ...
}

impl ParticleGroup for Molecule {
    // ...
}

// The current system
struct System {
    particles: Vec<Particle>,
    molecules: // ...
}

Is this how the "old way" would look like? After what you told me about the cache, I don't understand why there are particles inside of Molecule.

MoleculeRef can not be used separate of System, right? The particle reference would point to a slice of the particles vector in system, right?

[Luthaf]

Thank you! So in essence, in the most performance critical parts you want to have all information that you need very close in memory?

Yes.

Your idea with ParticleGroup is to have an interface where every implementation specifies the way to access the particles (and corresponding interactions)?

The initial idea was to abstract the particles and the topology of the system. Also having interactions would be nice for computing forces and energy, but I think it will prevent us from using the ParticleGroup trait for other structs than a System: if you need all the particles, the topology and the interactions, you would be better using a System containing only the particles you want.

Is this how the "old way" would look like? After what you told me about the cache, I don't understand why there are particles inside of Molecule.

What I call Molecule in this example is a new kind of molecule, let's says OwningMolecule living outside of a System. Such kind of molecule is useful when building systems, or for GCMC moves when you insert a molecule inside a system.

MoleculeRef can not be used separate of System, right? The particle reference would point to a slice of the particles vector in system, right?

Exactly. In this example, there are both a reference molecule type (like &str) and an owning molecule type (like String)

[g-bauer]

I see. The challenge will be to provide these functionalities in a structured -- in the best case, easy to use -- manner. To be honest, I have to look up the interfaces of a Systems methods every time I want to use them at the moment (.molecule(), .molecules(), .molecules[id], .molid(), ...). With that in mind

So is it worth it to add complexity in the library code in order to simplify the user code?

I'd say yes, if the "user code" involves writing stuff like new moves, propagators etc.

[Luthaf]

The challenge will be to provide these functionalities in a structured -- in the best case, easy to use -- manner.

Definitively! The initial idea would have been very hard to write, but I feel that we can manage to do something with the more recent ideas.

[g-bauer]

After reading your proposal again, I am not sure I fully understand Molecule and MoleculeRef. So my current understanding is : our "future" system will look like:

struct System {
    particles: Vec<Particle>,
    molecules: Vec<MoleculeRef>
}

where Particles are only "physically" present inside of particles and we store a reference to a slice inside the MoleculeRef. The sole purpose of a Molecule is then to have a self-contained (i.e. including its particles) representation, that we can use to build molecules outside of the system. Then, there will be a functionality to "consume" a Molecule to a MoleculeRef into the system I guess?

Anyways, I really like the idea of a trait as an interface to get the particles of a group (we should work on that asap, it would make things cleaner and more elegant). Maybe we could even add a function that gives access to positions directly (which from my understanding could not be a slice but a vec)?

[Luthaf]

So my current understanding is : our "future" system will look like:

Roughly. We can not have MoleculeRef directly inside the system because Rust does not allow to have both the owner and the borrower in the same struct. I was thinking to extract the Bonding informations and use it like this:

struct Bonding {
    range: Range<usize>,
    bonds: Vec<(usize, usize)>,
    angles: Vec<(usize, usize, usize)>,
    // ...
}

struct Molecule {
    particles: Vec<Particle>,
    bonding: Bonding
}

struct System {
    particles: Vec<Particle>,
    bondings: Vec<Bonding>
}

struct MoleculeRef<'a> {
    particles: &'a [Particles],
    bonding: &'a Bonding
}

impl System {
    fn add_molecule(&mut self, molecule: Molecule {
        self.particles.extend(molecule.particles);
        self.bonding.push(molecule.bonding);
    }

    fn molecule(self, id: usize) -> MoleculeRef {
        let bonding = &self.bondings[id];
        MoleculeRef {
            particles: &[self.particles[bonding.range],
            bonding: bonding
        }
    }
}

where Particles are only "physically" present inside of particles and we store a reference to a slice inside the MoleculeRef. The sole purpose of a Molecule is then to have a self-contained (i.e. including its particles) representation, that we can use to build molecules outside of the system. Then, there will be a functionality to "consume" a Molecule to a MoleculeRef into the system I guess?

That is right. But this design can be improved if anyone has more idea/suggestions.

Maybe we could even add a function that gives access to positions directly (which from my understanding could not be a slice but a vec)?

I don't know if this would be worth it, as this would allocate new memory and make it slower. If the users need to collect the positions, they can do it by themselves.

[g-bauer]

I though about the above again. From a "physical" point of view it makes a lot of sense to put all information of a particle together in a single struct, including position and velocities.

Coming from Fortran (and a lot of legacy codes) I am used to see everything stored separately in a "struct of arrays" style, i.e. there is a vector with all particles' positions, one with all velocities, one with all masses, etc. As far as I understood what you wrote about the cache, this is a very effective way to store the data. For example, a great part of the time is spend on computing distances, where a lot of extra information is not needed (masses, particle names, IDs, velocities ...).

I was wondering, did you think about using such a "struct of arrays" instead of "arrays of structs" approach when you designed sys? Is this a typical case of "clarity and easy-to-use over pure speed"? What made you chose the current way? What is your opinion on using rust's nice abstraction methods to use a different (more efficient**?** - would that be the case?) way to store data while still offering nice and easy-to-use APIs similar to what we currently do with molecules (which are effectively a "mask" to access certain particles out of an array of all particles)?

[Luthaf]

It is fun that you are asking about this, because I have a proof-of-concept of using the upcoming macros 1.1 feature (stable in rustc 1.15) to generate the interface code to have all the cache niceties of using a struct of array, while retaining the API niceties of the current array of structs approach.

I was inspired to look at this by this blog post, and it looks actually feasible. I will send the code to a branch here (it will only build using the beta compiler for a few more weeks), so that we can discuss the exact API.

[g-bauer]

=) I also read that blog post, after I came across a vlog of Jonathan Blow where he discussed SoA vs AoS in his implementation of a programming language for games.

because I have a proof-of-concept of using the upcoming macros 1.1 feature (stable in rustc 1.15) to generate the interface code to have all the cache niceties of using a struct of array

That's awesome! I am looking forward to see what you came up with.