New Feature: Support for bond-type interactions

So from talking to Rick yesterday, its clear that once we have a basic DPD implementation, the next thing that will be needed to do "short-range" only is the ability to create molecules, so we need to be able to handle bonds/angles/dihedrals/etc.

For a traditional implementation, this can basically just be a list of pairs/triplets/quartets, and can then be evaluated simply with:

for(int i = 0; i < nr_bonds; i++){
  xpart = &parts[bonds[i][0]];
  ypart = &parts[bonds[i][1][;
  evaluate_bonds(xpart, ypart);
}

For our setup this is not so straightforward, as we can't just access particles at random. However, what we can maybe do is have the bonds as a seperate region and field space, i.e.:

fspace bond( part_reg : region(part)){
  ipart : ptr(part, part_reg),
  jpart : ptr(part, part_reg)
}

task create_bond( parts : region(part), i : ptr(part, parts), j : ptr(part, parts)
  return [bond(parts)]{ ipart = i, jpart = j}
end

The difficulty is I have no idea how this will perform (or if it will do what I expect), but definitely a good starting point to try.

Initial ideas then:

[ ] Test if the recursive field space idea can work for some naive implementation.
[ ] If so, then begin building the infrastructure for how to build up bonds, and potentially run them in parallel.
[ ] Add a new invoke function for bonds, allowing user-defined functions on bonded pairs.
[ ] IO modules used for DPD/other methods need a way to build up these bonds automatically.

Once this works for bonds, we can extend it to triplets/quartets as desired.

Ok, so I struggled with the recursive field space idea for a while, and am not sure it is straightforward. I was thinking I could just use int1d to simulate the pointers to the particles. However, either of those ways is not sufficient.

The difficulty is due to the neighbour search algorithm for pairwise interactions being able to freely move the particles around the array (as required to avoid repartitioning). This allows parallleism and is essential for short-range forces, but looks to be making the bonded interactions way more difficult.

The big difficulty I can see, is if we have a mapping of:

bond.ipart => particle index

and the particle index changes, there's no straightforward way to update the "pointer" to the particle's new index.

Implementing a double sided link is also not straightforward, as particles appear multiple times in the bonded list potentially.

We could add additional functionality to create a globally available list of

particle ID => particle index

and keep this updated whenever particles are moved. However since this would use write I'm not sure its straightforward to do in parallel alongside the push/pull section of the code, so it would have to be an additional task, with partitioning etc. to enable this. It might not even be doable in parallel (which would be bad).

I think I'll have to implement a proof of concept for the list, and see if it works or if its horribly slow. Its an O(n) operation so in principle it should be cheapish, it just needs a good way to be parallelised.

The other concern is that the bonded tasks do not parallelise amazingly. It would almost be better if we could keep track of bonds as they are inside/between cells, and run them similarly to pairwise interactions. This would take way more bookkeeping but maybe lead to cleaner code, but doesn't help with the previous problem realistically.

Ok, so I looked back at old an old-task based MD code I was involved in, and the last thing is how it worked.

It had a structure called an engine_set, which contained:

  struct bond **bonds;
  struct angle **angles;
  struct dihedral **dihedrals;
  ...

These were of dimensions n_sets by nr_bonds.

After each timestep, there was a call to engine_bonded_reshuffle, which loops over the set of all bonds, and moves them to the appropriate engine_set (or sets, if accross a cell boundary) according to the cellpair containing the bond. The bonds are then done "asymmetrically" - in actuality it checks if both particles belong to the current engine_set and only updates the particles that do.

I think we can do a similar thing, though with a bit more infrastructure.

We start of by initialising some bonded array, of double the size or the number of bonds we expect to use. When creating a bond, we create two copies of it, one with each particle as part i and the other as part j. We then assign the cell the bond belongs to to part_i, and cell2 the bond belongs to to part_j.

We then create a cell partition of the bonds, and also as with the tradequeues for neighbour search, we allow for space to move bonds between cells, which we do in the same way.

The difficulty is then how to compute the bonds in a cell or pair of cells efficiently. The only way I can think to do it is to loop through the cell, and for each particle loop through the bonds to find all the relevant bonded interactions. The issue with this is the double loop ordering which is a bit horrible, but for now I don't have a better idea.

stfc / RegentParticleDSL

New Feature: Support for bond-type interactions #81