Performing efficient Widom insertions

[lequieu]

Hi All,

I have an existing workflow to compute the chemical potential of a system using Widom insertions. My current workflow first writes the trajectory of the system with N particles to a file and then performs the Widom insertions as post processing step (this was suggested in a prior post: https://groups.google.com/g/hoomd-users/c/dtm2jikp1k0/m/qjvlWvqwBwAJ). For each frame in the trajectory, I update the simulation state to contain N+1 particles and run for a single step with dt=0 to compute the new potential energy. This information is then used to get \Delta U(N -> N+1), which is then averaged over many frames/insertions to get the chemical potential.

My question is how to do this most efficiently. I am currently performing the insertions using state.set_snapshot() which can be quite slow and inefficient for large systems, especially on the GPU. Is there a better way to do this using the existing API? Clearly, generating the Widom insertions directly on the GPU would be superior, is there a simple existing way to do this? If not, any tips for how this could be efficiently implemented?

Many thanks!

[joaander]

Are these MD or HPMC simulations? What potentials do you need to evaluate for the inserted particles?

[lequieu]

These are MD simulations. Currently, particles interact via pair.Gaussian but I'd like to do this generally for other potentials if possible.

[joaander]

To get the fastest possible evaluation (insertions/second) with pairwise energies, I would follow the same general format as HPMC's FreeVolume computation on the GPU. This generates test particle positions (1 per GPU thread) randomly and then evaluates overlaps (energies in your case) with the nearby particles in the system found in the cell list.

However, my suggestion requires writing a significant amount of GPU code. Follow @mphoward's suggestion to implement something that works with the current Python API and will complete in a very reasonable amount of time.

I don't have any suggestions for arbitrary potentials. The specific details of the potential in question determine how you can (or cannot) evaluate many insertions in parallel.

[lequieu]

Ok, thanks for the tips. For now, I'll try @mphoward 's suggestion -- I think this should be sufficient for my problem at hand. Thanks!

[mphoward]

Looks like I wrote the comment you are referring to. :)

When I did this with HOOMD 2 for LJ-type particles, I was able to do the insertion quite efficiently by doing multiple insertions in parallel. I think I made the test particle its own type, and I disabled the interactions between test particles by setting r_cut=False for those pairs. For each frame, I added Ntest particles at random positions using the snapshot interface, and I computed the potential energy of the configuration. For pair potentials, the change in energy due to the insertion of a test particle should be 2x the energy assigned to the test particle (since half the pair energy goes to each particle in a pair). You can then use this information to compute the ensemble average you need over test particles and frames.

By doing multiple trials in parallel, the cost of the snapshot setting and transfer to the GPU is amortized, so I think it was pretty cheap in the end. The number of snapshots and trials we used should be in this paper:

https://doi.org/10.1063/1.5031789

[lequieu]

Neat! That's a pretty clever approach. Thanks!