Splitting long-range and bonded forces

Dear LAMMPS users and developers,

I am trying to use the python interface to lammps to evaluate atomic forces for my system, and I wish to split the force values into two terms: the long-range interaction involving pairwise and k_space contributions a well as short-range interactions involving 2, 3 and 4-body terms. I have tried two methods:

  1. Evaluate the total forces using

lammps_pot.lmp.command(‘run 0’)
f = lammps_pot.lmp.gather_atoms(“f”, 1, 3)

and then evaluate the short range interactions via:

lammps_pot.lmp.command(‘kspace_modify compute no’)
lammps_pot.lmp.command(‘pair_modify compute no tail no’)

lammps_pot.lmp.command(‘run 0’)

fb = lammps_pot.lmp.gather_atoms(“f”, 1, 3)

Obtaining the long-range contributions by total_force - short_range_force

  1. Evaluate total forces by

lammps_pot.lmp.command(‘run 0’)
f = lammps_pot.lmp.gather_atoms(“f”, 1, 3)

and then finding out the long-range interactions using:

lammps_pot.lmp.command(‘bond_coeff * 0 0 0’)
lammps_pot.lmp.command(‘angle_coeff * 0 0’)
lammps_pot.lmp.command(‘dihedral_coeff * 0 1 1’)
lammps_pot.lmp.command(‘improper_coeff * 0 0 0’)

lammps_pot.lmp.command(‘run 0’)

flr = lammps_pot.lmp.gather_atoms(“f”, 1, 3)

Obtaining the short-range forces using total_force - long_range_force.

It would seem that both methods should give me the correct answer, but they produce different results which lead me to believe that there is a flaw in my approach. Could you let me know which would be the best way to obtain those split-contribution forces and how should I modify my methods to get it? I think I would prefer something along the lines of options one since it should be slightly faster as there is no need to re-set forcefield parameters nor re-calculate pairwise interactions which should be the most expensive part of the force evaluation.

I am sorry if a similar question has been already answerd, I have tried browsing through the mailing list but came empty handed.
All the best,
Jacek Golebiowski

Dear LAMMPS users and developers,

I am trying to use the python interface to lammps to evaluate atomic
forces for my system, and I wish to split the force values into two terms:
the long-range interaction involving pairwise and k_space contributions a
well as short-range interactions involving 2, 3 and 4-body terms. I have
tried two methods:

do you need this information only for analysis purposes or during the run?​

in the first case, you can record a trajectory and then use the rerun
command with corresponding settings to only compute the bonded forces
(using pair style zero with a suitable cutoff and no kspace) and subtract
them from the total forces.

in the second case, i would consider instrumenting run style respa. that
one already keeps copies for forces of the various respa levels, and you
can also define levels with a factor of 1, which will recover a velocity
verlet integration scheme.

axel.

Dear Axel,

I need to access the forces throughout the run, however, only on selected timesteps. Would running the simulation with RESPA style with the options:

run_stype respa 2 1 bond 1 pair 2

be just as fast as the standard verlet run_style (accomplishes the same, I imagine), and if not, would the difference be significant?

In case the speed is comparable, would you know how is it possible to access the forces from two different RESPA levels separately using the python interface?

Thank you for the suggestions,

all the best,

Jacek

Dear Axel,

I need to access the forces throughout the run, however, only on selected
timesteps. Would running the simulation with RESPA style with the options:

run_stype respa 2 1 bond 1 pair 2

be just as fast as the standard verlet run_style (accomplishes the same, I
imagine), and if not, would the difference be significant?

​why don't you just make a test? it is trivial to do, won't take long, and
is authoritative (and not just me guessing). it will *definitely* beat your
very inefficient and badly parallelizing python hacks.

In case the speed is comparable, would you know how is it possible to
access the forces from two different RESPA levels separately using the
python interface?

​that would have to be programmed, first into the library interface and
then into the python wrapper.
however, if you worry about speed (and parallel efficiency), you don't want
to do that, but rather access this information internally.​

​in general, it is extremely difficult to give good advice without context.
e.g. a more parallel approach to solve the same issue would be to set up a
multi-partition run, where you do an on-the-fly equivalent of rerun for two
different settings based on the last timestep coordinates while advancing
the full interaction system concurrently.

axel.

Dear Axel,

I guess I was not clear enough, the majority of my simulation are timesteps when I do not need to know the separate force contributions. As a result, I am looking for a way of evaluating the long-range/short-range forces that will be easy to implement rather than fast but won’t slow down the rest of the standard run.

I agree that I should probably run the tests myself on my setup and see if this works for me.

Again, thank you for the suggestion, I will try to extend the library interface and work from there.

All the best,

Jacek

Axel’s earlier suggestion of using the rerun command sounds
best for what you are explaining. Just save dump snapshots
of the steps on which you want the detailed info. Then
do 2 reruns on just those snapshots to extract long vs short range
info.

Steve

Axel's earlier suggestion of using the rerun command sounds
best for what you are explaining. Just save dump snapshots
of the steps on which you want the detailed info. Then
do 2 reruns on just those snapshots to extract long vs short range
info.

​from the python interface, this could be done even more effectively:
a) create 2 (two!) instances of LAMMPS. one with the regular force field,
the second with pair style zero and no kspace style
b) run a chunk of trajectory with the full instance
c) extract the forces
d) save the positions through calling write_dump
e) call read_dump from the _other_ LAMMPS instance and do a run 0 to
compute forces
f) extract forces for the bonds and compute the forces for the non-bonded
interactions through subtracting the bond forces from the total forces
g) do whatever needs to be done with these two separate sets of forces
h) continue from point b)​

through using two LAMMPS instances, the costly, complicated and error prone
overhead of turning interactions on and off can be avoided. and you really
only need to recompute the info about one part of the two subsystems, since
you already know the sum. computing the second subsystem would only be
required as a test for consistency.

​axel.​

Dear Axel,