I am wondering if there is any way to compute the pair distribution function for pairs of atoms within a system in such a way that will exclude any i,j pairs which are within a certain number of bonds from one another?
Thanks in advance for any thoughts on this,
An idea that comments to mind is to save a file containing the position of the atoms for a representative amount of microstates and write a python code that is able to do the calculation of the RDF for you. When writing the python code, you would be able to identify these specific neighbors that you are planning on excluding from the RDF.
dump 2 all atom 500 dump.pos
fix 1 all npt temp 300 300 $(100.0*dt) aniso 0.0 0.0 $(1000.0*dt)
This would save you 1 microstate every 500 timesteps (so 1000 microstates in total) and save them all in a dump.pos file.
Please have a closer look at the compute rdf command. There is a note explaining what happens, if you have bonds and use the special_bonds settings. Unless you have long-range coulomb interactions, any 1-2, 1-3, or 1-4 pairs may be already excluded from the neighbor list.
If you do have long-range coulomb interactions, you can follow the suggestion of @ceciliaalvares and record a trajectory first (for an rdf calculation it is sufficient to store positions only every 1ps) and then you can use the “rerun” command together with using pair style zero (or lj/cut, but zero should be faster) instead of the original pair style, to recompute the rdf with the special neighbors excluded.
Thank you both for the suggestions!
I did spot the note on special_bonds, however for my system the culprit pair of atoms are not within the same bond angle or dihedral - so looks like I will have to solve this with python as @ceciliaalvares suggested.