# development question

Hi everyone,

We have a force field that has three different energy term types that depend on the simultaneous positions of all the nearest neighbors to an atom. The energy can be expressed as a gradient for each atom in our simulation.

We are considering using the following procedure to identify and calculate our values:

1. All of our energy terms need to employ neighbor lists. As we understand it, LAMMPS requires a pair style to generate neighbor lists, so we will program a custom pair style. Even though it is multi-body, our first energy term can be mathematically transformed to treat it in a manner consistent with LAMMPS pair styles.

2. We will use fix atom/property to declare custom variables for the gradients corresponding to this potential energy function. (We might call them d_vv_grad_x, d_vv_grad_y, and d_vv_grad_z).

3. The other two energy terms cannot be defined in a pair style, but they need the neighbor lists. Therefore, we will write a compute function to calculate the gradients and assign them to the atom properties.

This is where we have some questions.

1. Can we then use fix set force to apply these forces to our atoms in the simulation?

2. Could we do this before invoking the pair so they would add to the forces from our pair style? or would this overwrite them?

3. If this is wrong, what approach would you recommend we undertake?

thank you,

Matthew and Barry

Hi everyone,

We have a force field that has three different energy term types that depend on the simultaneous positions of all the nearest neighbors to an atom. The energy can be expressed as a gradient for each atom in our simulation.

sorry, this sentence doesn't make sense, first you say that the energy
terms depend on positions and then that they depend on gradients.
which one is true? gradients (with respect to positions) usually come
into play when determining the forces due to a given (position
dependent) energy term...

We are considering using the following procedure to identify and calculate our values:

1. All of our energy terms need to employ neighbor lists. As we understand it, LAMMPS requires a pair style to generate neighbor lists, so we will program a custom pair style. Even though it is multi-body, our first energy term can be mathematically transformed to treat it in a manner consistent with LAMMPS pair styles.

i don't understand this statement. there is no problem *at all* to
compute a multi-body interaction from a list of pairs. a triple of
atoms is just *two* pairs of neighbors of the central atom. so all you
need to do to compute this, is to add one additional level of looping:
i.e. loop over all central atoms (i), use the neighbor list to loop
over the first neighbors and then use *the same* neighbor list to loop
over the second neighbor. the only requirement for this to work well,
is that you request a *full* neighbor list, that is a neighbor list
where *all* pairs of a atoms are listed twice (i,j and j,i). the
default in LAMMPS are "half" neighbor lists, where each unique pair is
listed only once, since the pairwise force (e.g. with lj/cut) of the
i,j pair is the same than that if the j,i pair.

the only potential problem with this kind of approach is efficiency,
since neighbor list for (non-bonded) pairwise interactions will have
to reach to a cutoff way beyond the nearest neighbors (often 10 \AA or
more), but many-body terms are often much shorter ranged, doing the
nested loop on the "long" neighborlist requires a lot of needless
looping in the nested loops. this can be improved by building a
specific "short range" neighbor list from the already existing regular
neighbor list (several manybody pair styles employ this, e.g. airebo
or comb/comb3).

2. We will use fix atom/property to declare custom variables for the gradients corresponding to this potential energy function. (We might call them d_vv_grad_x, d_vv_grad_y, and d_vv_grad_z).

if those properties are recomputed in every step, this is not
necessary. you can maintain storage for that within the pair style. as
an example you can look at the EAM pair style, which requires a
density term due to its neighboring atoms and from that computes an
embedding term contribution, both of which are collected and stored in
per atom data (and simply require a forward and a reverse
communication to be properly distributed. other potentials have such
terms as well, but EAM is about the easiest to read and to understand
the principle. again, i am confused by the use of the term "gradient"
here. normally, i would expect that derivatives of the energy terms to
be computed analytically (with the help of an algebra software, if

3. The other two energy terms cannot be defined in a pair style, but they need the neighbor lists. Therefore, we will write a compute function to calculate the gradients and assign them to the atom properties.

you have to explain in more detail what you mean by that. as far as i
know, there is nothing that a compute would do that a pair style
cannot do. in the context of the rest of your e-mail, it seems that
you are mixing something up here.

This is where we have some questions.

1. Can we then use fix set force to apply these forces to our atoms in the simulation?

if you absolutely need to compute contributions to the force outside
the pair style (and i don't yet see a reason why this would be needed
here. the only case that i am aware of where this is used is to
compute CHARMM compatible cross terms, which are interactions on a
sequence of 5 linearly bonded backbone atoms, where that topology
would have to be entered in a data file as well), you want to
*implement* this as a fix, not a compute.

2. Could we do this before invoking the pair so they would add to the forces from our pair style? or would this overwrite them?

to add and not override (or set) forces you need to indeed *add* them
and not *set* them (hint, hint!). again, i believe this question is
caused by a more fundamental misunderstanding and thus using either
fix addforce or fix setforce should not be needed.

3. If this is wrong, what approach would you recommend we undertake?

i don't see anything in your *very* vague description of your model
that prohibits implementing it as a pair style in any way. you have to
provide proper and convincing evidence to the contrary. right now it
looks a lot like you get confused by nomenclature rather than looking
at the existing examples and understanding their flow of computation.

axel.

I'm probably not the best person to reply to your email but I decided
to reply first because I suspect other people have the same questions

We have a force field that has three different energy term types that depend on the simultaneous positions of all the nearest neighbors to an atom. The energy can be expressed as a gradient for each atom in our simulation.

Can you ellaborate a little bit on this?
Are the particles interacting with an external field of some kind?

We are considering using the following procedure to identify and calculate our values:

1. All of our energy terms need to employ neighbor lists. As we understand it, LAMMPS requires a pair style to generate neighbor lists, so we will program a custom pair style. Even though it is multi-body, our first energy term can be mathematically transformed to treat it in a manner consistent with LAMMPS pair styles.

2. We will use fix atom/property to declare custom variables for the gradients corresponding to this potential energy function. (We might call them d_vv_grad_x, d_vv_grad_y, and d_vv_grad_z).

3. The other two energy terms cannot be defined in a pair style, but they need the neighbor lists. Therefore, we will write a compute function to calculate the gradients and assign them to the atom properties.

I've always had this vague feeling that that any kind of local
(many-body) interaction between a collection of sufficiently nearby
atoms could be implemented as a pair_style in LAMMPS. Am I wrong
pair_style forces.)

This is where we have some questions.

1. Can we then use fix set force to apply these forces to our atoms in the simulation?

my guess would be no.
If you can not do this as a pair_style, then I suspect you will have
to create your own fix that modifies the forces and energies.
But I am not the one who should be replying to this email. Writing
your own fix to add forces to the particles is definitely doable, and
there are several fix_XXX.cpp files which implement this already. But
I'm not yet convinced that a fix is necessary. Why you can't
implement this as a new pair_style?

2. Could we do this before invoking the pair so they would add to the forces from our pair style? or would this overwrite them?

3. If this is wrong, what approach would you recommend we undertake?

Hopefully somebody with more knowledge will reply too.

Andrew

[...]

I've always had this vague feeling that that any kind of local
(many-body) interaction between a collection of sufficiently nearby
atoms could be implemented as a pair_style in LAMMPS. Am I wrong
pair_style forces.)

point out that even "long" pairs can be implemented as a pair style
(in fact, it is more of a problem to implement them as a bond style,
because bonds are limited to the communication cutoff, i.e. both atoms
of a bond have to be within the list of local and ghost atoms of a
sub-domain. the exception is implemented in pair style "list", which
is a bit of an anomaly (e.g. it is the only interaction in LAMMPS i
know of that *is* restricted by having to use minimum image
conventions). the curiousity of this is that implementing this as a
pair style turned out the only way to make such long distance directed
interactions work *at all* in LAMMPS. so i would say that rather than
stating that pair styles are limited, pair styles are the *least*
limited way to compute interactions and forces (of course pair style
list is rather inefficient, especially when running in parallel, but
that is the price to pay for being able to do something that would
otherwise be impossible to do at all).

This is where we have some questions.

1. Can we then use fix set force to apply these forces to our atoms in the simulation?

my guess would be no.

technically, the answer would be yes. for starters fix setforce *sets*
forces. it would be more natural to use fix addforce which *adds*
forces.
but it would in both cases require employing an atom style variable
that makes the output of that compute accessible and in the case of
fix setforce, also the existing force so far would have to be added.
but this whole discussion is moot. even though things are technically
possible, it is not necessarily the right way to do it.

axel.

To add a bit of context, In LAMMPS, a pair style is simply a black box that
is called once per Verlet timestep to compute the force
field. Originally that was only pairwise forces, but
now, as Axel said, there are lots of many-body force
fields implemented as well.

The black box takes atoms (owned and ghost) and a neighbor
list as input and returns force, energy, viral.
Inside the black box you can do that however you wish.
If you need to compute auxiliary quantities like gradients
you can allocate memory. You can communicate that
info to other procs, etc. Presumably the calculation is
only short-range, since that is what the neighbor list implicitly
encodes. Though there are pair styles which do charge equilibration
which is effectively a global operation.

The only reason to use a fix with a pair style is if there is
some new per-atom info it creates that needs to persist between time steps.
Most force fields are not defined that way. One exception
are granular force fields which have a history term for
shear friction. So the granular pair styles create a fix
to store that. The only reason to use a compute from
a pair style would be if you want to also access its calculated
info outside the pair style, e.g. as diagnostic output.
I canâ€™t think of any current pair style that does that.

Again as Axel said, nothing in your description sounds different
that what other pair styles in the code do by themselves. So I doubt
that you need an extra fix or computes or the fix atom/property command.
The latter is for per-atom values that persist from timestep to
timestep and are not defined already by the atom_style. And that
need to be accessed in multiple places, not just in a pair style.

HTH,
Steve