Most convenient way to compute custom pair-level calculations without saving trajectories?

[CreditDefaultSwap]

Hello. I would like to compute some custom pair-based quantity, say,
V=sum_(i>j) v(rij)
where rij is the pair distance and v(.) is some arbitrary function. I am wondering what would be the easiest and fastest way to do this in hoomd directly at the simulation step, rather than saving trajectory files? The reason being that I would be simulating very large systems for a long time and would rather save single scalar values of what I need, rather than entire trajectories for later analysis.

Would doing a custom logging action and computing v(.) from the neighbor lists be the best way? Or is there some other easy, built in way of doing this that I somehow missed in the documentation.

Thank you for the tips and help

[joaander]

Would doing a custom logging action and computing v(.) from the neighbor lists be the best way?

Yes. This is one of the main reasons why HOOMD-blue has a Python interface.

Or is there some other easy, built in way of doing this that I somehow missed in the documentation.

No. Users are expected to implement custom operations. HOOMD-blue implements only the methods necessary to perform MD and MC simulations. All the built in optimizations and data structures are heavily optimized to provide the best performance in these cases.

[joaander]

self is the class instance, as describe in the Python documentation. You would pass an instance of Cell to your action's __init__ which sets self.cpu_local_nlist_arrays.

[CreditDefaultSwap]

Thanks again for your answer. As I realize, this task is a more non-trivial than it seems and information in the docs is a bit challenging at times. I have figured out how to pass the cell list and compute it locally, at the MPI rank level, as follows:

 with self._state.cpu_local_snapshot as snap:
            with self.nlist.cpu_local_nlist_arrays as arrays: #returns neighbors lists in local index, including ghosts particles
                nlist_iter = zip(arrays.head_list, arrays.n_neigh)
                for i, (head, nn) in enumerate(nlist_iter):
                    pos_i = snap.particles.position[i]
                    for j_idx in range(head, head + nn):
                        j = arrays.nlist[j_idx] 
                        pos_j = snap.particles.position_with_ghost[j]; 
                        #further operations to calculate arbitrary v(rij)
 

What I noticed is that, because some particles in the neighbor list may not be in that particular rank, I have to ensure to call position with ghost particles listed. This then allows me to compute all relevant rij quantities. Following this, however, I have the two problems I couldn't really find any info on the docs about:

1. Neighbor pairs appear twice on each rank i.e. if I summing over quantities I'd be double summing them when adding all contributions from each rank. How can I deal with this issue, or is there some other explicit way of doing this type of collective additions.
2. how to collect the results from all independent ranks and combine them into a single value. For instance this code should be able to calculate the energy, but once I've calculate whatever local contribution from each rank, the final value of energy remains isolated to each rank. Is there something like, or is there any way to do the equivalent of an MPI 'gather' in hoomd following this procedure?

I much appreciate all the help so far

[joaander]

Please understand that HOOMD-blue is built for high performance first and for a clean/extendable API second. By opting into using the high performance zero-copy data structures, your Python code must now adhere to esoteric data layout (and duplication) decisions made to enable high performance simulations at the C++ and GPU kernel level. I assume that users of the zero-copy API are familiar with most of the C++ implementation details and are choosing to prototype an algorithm in Python that will later be ported to a C++ component if it proves successful. To get the most performance out of the zero-copy APIs, you should pass the arrays directly to numpy whenever possible and avoid writing for loops in Python.

If performance is less of a concern for you, you can use local_pair_list with MPI decomposition or pair_list on rank 0 only for a more Pythonic approach.

1. Neighbor pairs appear twice on each rank i.e. if I summing over quantities I'd be double summing them when adding all contributions from each rank. How can I deal with this issue, or is there some other explicit way of doing this type of collective additions.

When using the GPU, the internal neighbor list stores both i,j and j,i neighbors. It stores only i,j neighbors when using a CPU device. The simplest way to handle only unique pairs is to check i < j.

2. how to collect the results from all independent ranks and combine them into a single value. For instance this code should be able to calculate the energy, but once I've calculate whatever local contribution from each rank, the final value of energy remains isolated to each rank. Is there something like, or is there any way to do the equivalent of an MPI 'gather' in hoomd following this procedure?

Use mpi4py. You can pass the MPI4py communicator to Communicator (and then pass that communicator to Device) to ensure that HOOMD-blue's internal domain decomposition strategy communicates with the same communicator you are using.

[CreditDefaultSwap]

Just wanted to confirm that I got this to work, though it is certainly still not as fast I thought. Either way I appreciate the help and support given!

[joaander]

I'm not surprised. We implemented these APIs with the best intentions, but in practice they rarely deliver. Stay tuned to this discussion board - we are currently working on some new designs to solve these problems.