Why inner, middle, and outer keyword in respa does not support eam potentials

Dear Lammps users,

I’m trying to repeat some simulations from the literature work using respa with eam potentials. But it seems that the force decomposition algorithm (with inner, middle, and outer keyword in respa) has not been implemented for eam potential in lammps, but implemented for LJ potential. I wonder why? Is it because technically it is very difficult or even impossible to decompose forces from eam potential due to its manybody feature, where the atoms in one layer will also affect the force contribution from atoms in another layer?

Thank you,

Zhiliang

Dear Lammps users,

I’m trying to repeat some simulations from the literature work using respa with eam potentials. But it seems that the

was this simulation performed with LAMMPS? what are the specific
directions how r-RESPA was applied, what cutoff ranges and what
timesteps?

force decomposition algorithm (with inner, middle, and outer keyword in respa) has not been implemented for eam
potential in lammps, but implemented for LJ potential. I wonder why? Is it because technically it is very difficult or even impossible to decompose forces from eam potential due to its manybody feature, where the atoms in one layer will also affect the force contribution from atoms in another layer?

it is technically difficult, and the benefit is limited.

1) for 3-tier inner/middle/outer pairwise interactions to work
efficiently, you need a significant separation and geometry changes
between the distance ranges, and that requires using a rather long
outer cutoff. EAM potentials, have a quite short cutoff, so that won't
make sense. you need some overlap region to smoothly switch between
the distance ranges so you retain sufficient energy conservation.

2) for 2-tier inner/outer decomposition, it is more likely to have a
(small) performance benefit from r-RESPA. in fact, i find this often
to be more efficient than the three level setup even for lj/cut, but
again, somebody needs to program the switching function.

3) the eam force computation is a 2-step process. you have to walk the
neighbor list twice. you first compute the embedding energy
contribution from the neighboring atoms, make a reverse communication
(i.e. reduction) and then compute the embedding term contribution and
do a forward communication before computing forces and energies from a
pairwise term and the embedding term. so you need that first step in
any case, or figure out a way to decompose and store those
contributions for each r-RESPA level in addition to the forces and
avoid having to recompute them. LAMMPS is not set up for that
currently.

in short, it would be a lot of work, but it is not obvious, if there
would be a sufficient benefit for typical EAM potentials the way
LAMMPS is currently programmed. in today's HPC environment, it is much
easier to obtain access to additional computational resources, and
computers are much more capable, than they were when r-RESPA was a big
benefit.
nowadays, you can probably gain almost as much performance, if not
more by simply using MPI+OpenMP or MPI+GPU MPI+KOKKOS or
MPI+USER-INTEL with plain verlet time integration.

axel.

Thank you Axel for the detailed information. It helped me a lot.

Zhiliang

Thank you Axel for the detailed information. It helped me a lot.

then why don't you repay the favor and answer *my* questions?

axel.

Sorry Axel,

I looked directly into your answer at the end of my message and didn’t notice your question in the middle.

The simulation was performed on lammps, that’s why I asked my question here. I was actually trying to reproduce some simulation results from a literature work since their results is different with mine. I have to address this issue (brought by the reviewer) to get our paper published. The system is nanostructured Cu modeled with EAM potential. The authors claimed that they used multiple timestep on EAM potentials with their own MD code. The time step is 2 and 6 fs, but they did not mention cutoff at all. I tried respa in lammps, but did not reproduce their results. I think it should be that their MD code support force decomposition with EAM potentials, but lammps does not.

Zhiliang

Sorry Axel,

I looked directly into your answer at the end of my message and didn’t notice your question in the middle.

The simulation was performed on lammps, that’s why I asked my question here. I was actually trying to reproduce some simulation results from a literature work since their results is different with mine. I have to address this issue (brought by the reviewer) to get our paper published. The system is nanostructured Cu modeled with EAM potential. The authors claimed that they used multiple timestep on EAM potentials with their own MD code. The time step is 2 and 6 fs, but they did not mention cutoff at all. I tried respa in lammps, but did not reproduce their results. I think it should be that their MD code support force decomposition with EAM potentials, but lammps does not.

if the statistics are sufficiently well converged for either
simulation and the simulation parameters adequate, it should not
matter whether you use r-RESPA or (velocity) verlet for time
integration, or whether you have domain or force decomposition.

FWIW, the cutoff for EAM potentials is embedded in the potential file.
it is the 5th datum on the 3rd line.

i am very curious how you would be able to do a meaningful
multi-timestepping with EAM potentials, as their cutoffs are usually
too short. to have a significant separation of time scales. i suggest
you have a closer look at the reference publication and contact the
authors, if needed. using an in-house MD code always has a much higher
risk to have undetected bugs in the code than in community maintained
codes. i am constantly surprised that we or others sometimes find bugs
in parts of LAMMPS have been used for a long time and by quite a few
people. that said, not every unconfirmed simulation result is due to
bugs in code, incomplete sampling and bad statistics or unfortunate
choice of simulation settings can cause this just as well.

axel.