Dear LAMMPS developers and users,
From time to time, we may introduce certain substances that have to be carefully treated with a much smaller time step than the original system, thus leading to a much smaller global time step. This can bring about a great amount of unnecessary computation in terms of the original substance. I am not sure whether this circumstatnce exists in molecular dynamics. But for users interested in mesocale simulation (or rather, smoothed particle hydrodynamics, SPH) like me, a large time step can be used for liquid but a much smaller time step has to be set for the gas. As a consequence, one has to use a very small time step in gas-liquid two-phase flow even if there is only very little amount of gas.
This prompts me to ponder whether we can set different time steps for different particles. LAMMPS has supported RESPA currently that is also a sort of multi-timescale integration. But RESPA is not applicable to my case because different time steps should be set for different atoms, rather than different kinds of forces (RESPA). It is easy to determine if a particle’s info should be updated in SPH. A peratom quantity (dstep) can be defined to describe how many time steps elapsed since the last update of each particle. If a particle is not gas, has no gaseous neighbors, and dstep < stepmax, then the particle’s info does not need to be updated in the current time step. Otherwise, update its info. I have written a new compute style update to implement this (see the attatched files).
compute_update.cpp (4.6 KB)
compute_update.h (1.2 KB)
However, I soon find there can be many challenges. I have not yet digged into the source code of RESPA. For the traditional Verlet scheme, the forces in the last time step are cleared before force computation. So I must store the forces at each timestep. Worse still, some pair styles computes not only the forces but also heat flux and density change (see pair_style sph/taitwater/morris, and pair_style rheo). Then it seems I have to store heat flux and density change too. Also, it seems I need to modify each fix/pair/bond/angle/… style I am going to use to disable the update of specified particles. Maybe a tedious work…
So I am writing here to ask if this practice is worthwhile, how to implement it in the most convenient way, and whether there can be better ways for speedup. Any suggestions or comments are appreciated!
The situation is not quite as simple as you make it, because you also have to handle the case where the fast particles are interacting with the slow particles.
There has been a precedent, where an all-atom simulation (short timestep) was combined with a coarse grain model (longer timestep) and the assignment to the different sub-systems (all-atom, coarse grain, mixed interations) was done with different pair style hybrid sub-styles. It would also be possible to do something similar with (non-overlapping) groups.
The basic idea of r-RESPA is that you compute the forces/interactions at a specific level and store them in an internal fix. The time stepping itself is done for the smallest timestep, but the forces are copied from the internal fix if they don’t need to be computed at a given level for the specific time step. So you avoid looping over the neighbors of each particle and just copy the forces from the previous step.
Whether this is sufficient for your use case is something that you may have to figure out and, if necessary, perform the relevant derivations and tests. There is detailed literature on this topic.
If you need to cache more properties, you should create a custom run style with the corresponding internal fix.
The fundamental point is that it is not about how often you update the positions, but how often you (need to) compute the interactions.
1 Like
Thank you for your quick reply.
As to this precedent, is there any publication or source code for reference?
I will spend some time figuring out how r-RESPA works.
Yes. Interaction computation is the most expensive. I plan to update the positions at each small time step, and update interactions asynchronously.
There is the “hybrid” keyword in the “run_style respa” command which is implemented in the src/respa.cpp file.
This was added about 10 years ago in the 28 Jul 2015 LAMMPS release which has this description:
Sam Genheden (U of Southampton) added the capability to rRESPA via the “run_style respa” (run_style command — LAMMPS documentation) command to use “pair_style hybrid” (pair_style hybrid command — LAMMPS documentation) where different sub-styles are invoked at different rRESPA levels. This is useful for hybrid coarse-grained/all-atom models.
I suggest you try to look up Sam’s publications or those of his then advisor, Mario Orsi, who developed the ELBA coarse-grain force field for which this was used.
1 Like
Got it. I will read them. Thanks again.