Separate System and Interactions

[antoinewdg]

If we want to do any parallelization, we will need the System object to be Send and Sync, and it is not Sync because of the RefCells used in Interaction.
The good news is that we actually do not need to have the interactions passed around threads, so making the system (here the configuration of the particles) and the interactions two separate entities would solve the problem. The bad news is that it would require some non trivial refactoring. I'll work on a draft for this and open a PR.

As a side note : according to the rust doc, it would make more sense if GlobalPotential::energy (and others) were immutable and the RefCell was actually inside the struct implementing GlobalPotential if it were to use caching, the System should not be aware of this mutability. Is there a particular reason that it is done the other way around ?

[antoinewdg]

Would having a separation like that be OK ?

#[derive(Clone)]
pub struct SystemConfiguration {
    /// Unit cell of the system
    cell: UnitCell,
    /// List of particles in the system
    particles: Vec<Particle>,
    /// Molecules in the system
    molecules: Vec<Molecule>,
    /// Molecules indexes for all the particles
    molids: Vec<usize>,
}
/// The System type hold all the data about a simulated system.
///
/// This data contains:
///
///   - an unit cell, containing the system;
///   - a list of particles in the system;
///   - a list of molecules in the system;
///   - a list of interactions, associating particles kinds and potentials
///   - a hash map associating particles names and particles kinds.
///
/// In the implementation, the particles contained in a molecule are guaranteed
/// to be contiguous in memory. This allow for faster access when iterating over
/// molecules, and easier molecule removal in the system.
#[derive(Clone)]
pub struct System {
    configuration: SystemConfiguration,
    /// Interactions manages the associations between particles and potentials
    interactions: Interactions,
    /// Current step of the simulation
    step: u64,
    /// Externally managed temperature for the system, if the propagation
    /// algorithm does not update the velocities.
    external_temperature: Option<f64>
}

[g-bauer]

In this context maybe it is good to talk about #18 again? Also @Luthaf and I had a conversation about how to represent molecules, how to get bonding information from the input and how to apply a move for a certain molecule. While this is not directly connected to parallelization, it would also change the structure of System and maybe Interactions.

@antoinewdg What would be the first step? Using something like rayon with the energy and force loops?

[antoinewdg]

In this context maybe it is good to talk about #18 again?

Probably. My take on this is that no matter what we do, the system should keep this internally and expose a clear interface, and that it we are getting closer and closer to a moment where we would have to sit down and think about it.

What would be the first step? Using something like rayon with the energy and force loops?

Exactly. Energy and virial computation would be simpler because they are just a simple reduction, while forces are a little subtler (avoiding computation of the same force twice as we can do sequentially does not seem easy).

[g-bauer]

(avoiding computation of the same force twice as we can do sequentially does not seem easy)

I'm not sure I understand this? We also have this issue with energies where the outer loop is over all particles minus one (with index i) and the inner loop over all particles starting at j = i+1. Or do I misunderstand you here?

[antoinewdg]

The issue is actually not computing the values forces for each pair, it's the aggregation that causes problems: energy and virial are just a big reduce operation over all pairs, while forces is a weird scatter.

[g-bauer]

Ah, I see what you mean 😄

[Luthaf]

Global interactions are inside a RefCell because it was faster for Ewald summation to just check the borrow state once. As you say, it should be moved to something like Ewald = RefCell<InnerEwald> or so, to handle caching borrowing. I believe a simple way to fix the Sync issue would be to use a Mutex or a RwLock instead of the RefCell

The issue I have with a separation of System and Interactions is that they are tightly coupled right now: the system setup the particles kinds from the interactions, and will only work with one interaction instance. We could decouple them if we don't get a performance penalty, maybe by changing how we store and represent interactions.

And yes, we should definitively think again about #18, with requirements from the parallelism (we may want to use multiple subsystems for MPI style parallelisation); possibility to refer to molecules in the system by molecular species; possibility to have SoA style optimization; and offering a nice API. I don't know how much these are coupled together, and how many trade-offs there will be.

Maybe we could have a video call at some point to actually sit down and discuss all of this?

[antoinewdg]

The issue I have with a separation of System and Interactions is that they are tightly coupled right now

Indeed this is a problem I did not see. Mutex or a RwLock would be a solution but we really would not want them in the middle of parallel loops (badly used they would slow things down, and we do not need them).
What we could also do is restrict the interactions stored in the BtreeMaps to be Sync (they would not be allowed to cache) and have the system exposes the maps (maybe behind a trait to hide implementation details) so that the compute functions could pass them around threads along with the SystemConfiguration.

Maybe we could have a video call at some point to actually sit down and discuss all of this?

That is definitely a good idea !

[Luthaf]

Sorry, I might be too tired, but I really don't get your solution. The idea is to have only Sync Interactions, and pass the (read-only) interactions maps aside with the systems? Then how would one implement caching with this?

[antoinewdg]

Sorry, I might be too tired, but I really don't get your solution.

No problem, I'll elaborate a bit:

• the SystemConfiguration I suggested above is totally Sync, which is good because it is needed by almost every computation and will be passed around
• GlobalPotentials (including Coulomb) are not Sync, which is no problem : how I see it (I may be wrong) they are called sequentially by the system and implement parallelism internally if needed
• pair, bond, angle and dihedral potentials, on the other hand, need to be Sync: the parallelism is implemented by the system (in the various Compute loops) and can not be allowed to cache

[g-bauer]

Maybe we could have a video call at some point to actually sit down and discuss all of this?

👍

[antoinewdg]

Can we organize this call in the not too far future ? Parallelizing stuff is my original goal, and I can't work on it as long a we do not make these changes.
I have of quite flexible timetable so maybe you two try to find a good time and I'll adapt. If you think that is useful I can make a more elaborate proposal and we could use it as a base to iterate on during the call.

[Luthaf]

Can we organize this call in the not too far future ?

Yes, sure. I'll send an email about it.

Parallelizing stuff is my original goal, and I can't work on it as long a we do not make these changes.

I disagree here =). We can go for incremental changes, and just start with the non-optimal way of using RwLock and iterate from here. I think refactoring the system is needed, but it will take a lot of time, and it does not prevent you from working on parallelisation in parallel!

[antoinewdg]

We can go for incremental changes, and just start with the non-optimal way of using RwLock and iterate from here.

That's totally right, thanks for pointing it out !

[Luthaf]

Changing the issue to track the System/Interactions separation now that the system is Sync (#132).

---

As I said during the video call, this whole issue comes from the need to cache Ewald computation in MC moves. But the cost and update methods of the cache takes &System and can not mutate the Interactions inside the System. Hence the RefCell and now RwLock.

We can resolve this by separating the System in two as proposed in #126 (comment), and then change the EnergyCache function to take separated arguments for configuration and interactions.

The main issue with this way of doing this is the atomic names <=> atomic kind association. We need to be sure a specific atomic name correspond to the right atomic kind, as both are stored in the Atom, while the canonical source for this information is in the Interactions.