[BUG] Compute RDF unphysical with neighbor multi

Summary

compute rdf produces incorrect results with neighbor multi, the curve decays to zero. My scenario is a WCA liquid for which I want an RDF for comparably long distances. The issue does not occur with neighbor bin.

LAMMPS Version and Platform

7 Feb 2024 - Development - patch_7Feb2024_update1-147-g4d89741d8c

Expected Behavior

neighbor multi should lead to the same RDF as neighbor bin.

Actual Behavior

The RDF produced with neighbor multi decays to zero as shown in the plot.

Steps to Reproduce

A simplified input script is attached (in.spheres.txt). Equilibrate a system of WCA spheres. Next, make sure to use neighbor bin and measure the RDF with compute rdf. Subsequently, turn to neighbor multi and measure it again. The discrepancies can be spotted directly in the output files without plotting if desired.

Further Information, Files, and Links

Longer force cutoffs can solve the issue. Note that if I calculate the RDF in post-processing, it’s fine in both cases.

Thank you for the clear bug report and simple example. I was able to reproduce the bug in the current development branch of LAMMPS.

I narrowed down the issue to the fact that custom cutoffs in neighbor requests, NeighRequest->cutoff, are not propagated to the multi (or multi/old) stencils since they don't use the NStencil->cutneighmaxsq variable. There appears to be a note alluding to this in nstencil.cpp:

lammps/src/nstencil.cpp

Lines 167 to 168 in 4d89741

    
           // overwrite Neighbor cutoff with custom value set by requestor 
        
           // only works for style = BIN (checked by Neighbor class)

@sjplimp Do you know whether there is an (old) error message this comment is referencing or if there was a fundamental reason only BIN was supported? It looks like the original multi/old was not supported either.

I don't see an obvious reason why we can't update the stencil cutoffs (NStencil->cutcollectionsq for multi and NStencil->cuttypesq for multi/old) to account for a custom cutoff request (other than performance), but maybe there's a gotcha.

I'll do some digging into older versions of neighbor.cpp in the meantime.

@jtclemm - I think this was not supported simply b/c no one had ever asked for it - based on that comment I think it probably threw an error, or possibly built the cutoff-requested neighbor list with bin, instead of multi

Given that a multi system is likely to have multiple type-dependent cutoffs, what is the desired behavior with compute rdf, which just allows specifying a single (presumbably longer) cutoff ? If that cutoff is applied to all type pairs, won't the neighbor list potentially be huge ? Or should the specified value be like an extra skin distance applied to the cutoff of each type pair ?

Given that a multi system is likely to have multiple type-dependent cutoffs, what is the desired behavior with compute rdf, which just allows specifying a single (presumbably longer) cutoff ? If that cutoff is applied to all type pairs, won't the neighbor list potentially be huge ? Or should the specified value be like an extra skin distance applied to the cutoff of each type pair ?

If we added a skin/scaled cutoffs, how would we handle RDFs calculated using all distribution atom types? E.g. in a binary system, if we want the full RDF for large particles with a scale factor $\lambda$. With cutoffs $\lambda r_{LS}$ and $\lambda r_{LL}$, the computed RDF would include both types up to $\lambda r_{LS}$ then only include large particles up to $\lambda r_{LL}$. Would we forcibly truncate the RDF at the minimum interaction distance $\lambda r_{LS}$ (or $r_{LS} + \epsilon$ for a skin $\epsilon$)?

For a fixed cutoff, can we assume the user picks a reasonable value? And if they don't, then the neighbor list would just overflow and they'd need to adjust parameters?

For a fixed cutoff, can we assume the user picks a reasonable value?

We do this now. People can and have tried using insanely large cutoffs for compute rdf and were complaining about the slowdown and memory requirements. There is no stopping people making bad choices. You can do the same with many settings. So for compute rdf (and some other computes that use neighbor lists and may have explicit cutoffs with their neighbor list requests, I know there are some modifications along those lines in preparation) if you specify a cutoff explicitly that should produce a "bin" style neighbor list with that cutoff and that is that. Anything else conflicts with the principle of least surprise.

What is done without specifying an explicit cutoff is an open question as far as I am concerned. Perhaps, the compute rdf code should have a requirement that either "bin" or "nsq" neighbor lists are used globally or use an explicit cutoff.

if you specify a cutoff explicitly that should produce a "bin" style neighbor list with that cutoff and that is that.

This might get complicated. For instance, if someone is using the multi comm style then neighbors won't be communicated as expected for a "bin" style neighbor list. In fact, this is a complication for Steve's skin suggestion as well...

I'm not sure there is an obvious solution that wouldn't be a major development task touching many parts of the code. In the meantime, what if we just

Error if RDF is used with multi
Update NeighRequest->choose_stencil() to not match req->cut with style multi or multi_old
Add a "volunteer needed" feature request to add multi support to compute RDF

Alternatively for (1), we could just error if a custom cutoff is requested with multi and add a warning to the code/documentation explaining that one has to be careful using compute RDF with multi since the cutoffs of the specified types may not correspond to the cutoffs of the collections.

Since the RDF doesn't need to be done frequently, it would be quite feasible to compute the RDF using rerun with "bin" neighbor lists after running the simulation with "multi" and keeping a trajectory. So I would go for a hard error.

@sjplimp @jtclemm I have added a change to PR #4094 that will disallow using a custom cutoff in compute rdf with "multi" neighbor lists in commit 45e8ee1. I am currently looking into other computes and fixes that request neighbor lists with a custom cutoff, that would be affected as well; for example compute adf. If I cannot tell for certain, I will add a warning instead of an error. @jtclemm can you have a look at fix nonaffine/displacement?

@ellio167 @tadmor This issue could also affect pair style kim. It will probably not affect the correctness but could negatively affect performance. Have you ever tested some models with neighbor style multi?

@jtclemm can you have a look at fix nonaffine/displacement?

Hey @akohlmey, I made the necessary changes in this branch: https://github.com/jtclemm/lammps/tree/updating-nonaffine-disp. Looking through the fix, I made some additional changes to clean it up (e.g. using less ambiguous variable names and enabling some unnecessarily disabled options). Should I merge the changes in your PR or just create a separate PR?

In addition to these error checks in different fixes/computes/pair styles, would it be good to add a fail-safe to neighbor->choose_stencil() to prevent req->cut from matching with anything but bin styles? With a clarifying comment explaining why, something to the effect of:

"Custom cutoffs currently only impact bin-style nstencils by overwriting cutneighmaxsq. To apply to multi and multi/old stencils, cutcollectionsq and cuttypesq need to be updated, respectively, and additional checks need to be performed to ensure all atoms are communicated appropriately with the multi/reduce comm style. Given the nature of multi, an alternative solution (e.g. scaling cutoffs by a fixed factor) maybe more useful."

@jtclemm thanks, I just pulled them into my branch.

I don't think additional changes are needed. This is a very unusual case:

$ grep -l set_cutoff *.cpp */*.cpp
compute_rdf.cpp
neigh_request.cpp
pair_hybrid.cpp
EXTRA-COMPUTE/compute_adf.cpp
EXTRA-COMPUTE/compute_ave_sphere_atom.cpp
EXTRA-COMPUTE/compute_composition_atom.cpp
EXTRA-COMPUTE/compute_efield_wolf_atom.cpp
EXTRA-FIX/fix_nonaffine_displacement.cpp
GPU/pair_lj_cut_tip4p_long_gpu.cpp
KIM/pair_kim.cpp
REPLICA/fix_hyper_local.cpp

We have taken care of almost all cases. The GPU package does not allow "multi" neighbor lists, for pair style kim, "multi" should not make a difference depending on what I see in the code (it will only be slower) and fix hyper/local is such a special case, that I don't consider it to happen in our lifetime. The most likely cases are the computes.

@jtclemm BTW, you should be able to check out most pull requests, modify and commit and then push them back. This applies to all branches submitted to a pull requests for LAMMPS members with write privilege and where the submitter has allowed pushes (either directly or through rules for the group account).
This is best done with the github CLI: https://cli.github.com/
After authenticating with "gh login" (only needs to be done once per checkout) you can check out and modify pull requests

gh co 4094
[...]
git commit
git push

It often saves a lot of time doing small changes like this than having to explain and request changes and then re-check.

his applies to all branches submitted to a pull requests for LAMMPS members with write privilege and where the submitter has allowed pushes

Ah, that's very helpful to know. Thanks for the tip, I'll try it out next time.

	// overwrite Neighbor cutoff with custom value set by requestor
	// only works for style = BIN (checked by Neighbor class)