I am a first year PhD student and I am learning to use LAMMPS to do part of my simulations. I want to simulate active particles inside a fluid, and because the fluid interactions are DPD, it is desirable that momentum is conserved and the system is torque-free. My approach was to once active force acts on active particles, the same but opposite force should be applied to the fluid. So my problem comes here, on how to identify the two nearest particles of a specific particle.
I thought that maybe with the “bond/create” I could define a certain set of bonds that combined with the “bond_style zero” and “compute bond/local dist” this could be managed. Also combining dynamic groups, where each particle has a set of neighbor particles that is being updated each time step. More than a specific answer for this problem, I would like to know more about how this type of scenarios are typically treated in LAMMPS and if there is some usual procedure for this.
Anything of this kind using the LAMMPS scripting capabilities would be prohibitively slow since the scripting commands are executed globally and not in parallel and require communication for each step.
This is a surprisingly costly problem, since you first need to sort all neighbors by distance to know the N closest neighbors. So in comparison to finding all atoms within a cutoff distance from another, you would need an additional step and additional storage. In addition, there is the issue that in many cases, LAMMPS will use Newton’s third law to reduce the number of pairwise computations by half and thus you won’t find atom A in the neighbor list of atom B and at the same time atom B in the neighbor list of atom A. So additional communication would be needed or also a full neighbor list.
This would have to be handled in the C++ code of the force computation of a pair style. I don’t recall anybody asking about this exact issue. Since DPD is a rather approximate method and you also have random noise as part of your model, I would assume that people have determined for their applications that the need to enforce the kind of balance you want to employ is overshadowed by the dissipative contributions to the interactions in the DPD model and thus not relevant.
Adding bonded interactions has a bunch of implications. It will remove so-called special neighbors from the pairwise neighbor list (unless the special_bonds command is used to change the default settings) and bonds must not only be formed, but also broken, since atoms would diffuse away and thus no longer qualify. Also the way that LAMMPS computes forces requires them to be “local”. That is, when atoms with bonds can freely move around, you will eventually get “bond atom missing” errors.
Okay, thanks for the extensive and detailed answer. I understand that is something maybe too specific and that my initial approach was probably too vague. I will look if there is a way I can implement this from the source code and also look for other people treating similar systems.