Implementing Body Forces within a Pair Style

Hello everyone,

I am currently working on my own Magnetoelastic pair potential within the SPIN package. I’m utilizing the a modified version of the magnetoelastic energy equation from Kittel (see here) and an atomistic strain tenson developed by Zimmerman (see here).

In my current implementation of the magnetoelastic energy equation, my magnetic effective field force is calculated through the relative displacement of my neighboring atoms (hence the use of the pair style) while the corresponding mechanical force is due to the orientation of the current atom’s spin vector only.

I’ve been looking for a way to implement my mechancial force as a body force akin to what is done in ‘fix addforce’ within the pair_style framework. However I’m having trouble on how to properly implement it (which version of ev_tally to use, ensuring external magnetic field energy is properly encapsulated, etc).

Any help is greatly appreciated, thank you.

Why not in a fix? You can request a neighor list from a fix just as well and then do a calculation over neighbors just as in a pair style?

In general, if you want to add properties computed in a pair style to other pair style, your pair style has to be used with pair_style hybrid/overlay to be added to another pair style.

Hello Axel,

Thank you for the quick response. I had not known that was possible to request neighbor lists from within a fix, I apologies. This is my first time attempting developer work on LAMMPS. I had originally used the Spin/Exchange pair style as a reference when implementing my energy equation, I’ll look into converting it into a fix as that could be the solution to my mechanical force implementation issue.

To you second statement. I have been using a combination of the Spin/Exchange, MEAM, and my magnetoelastic potential as a pair style through the use of hybrid/overlay, along with the Zeeman and Cubic options in fix precession/spin (with fix_modify energy yes enabled). I was running into issue that my energy was not being calculated as evflag would stay at zero (I’m thinking this is due to my energy being based on a reference configuration). If I had wanted to use the energy generated from the Zeeman term of the fix precession/spin due to the application of an external magnetic field, could I still call it within my magnetoelastic pair_style?

Thank you again for your help.

eflag and vflag are bitmasks that tell whether any compute or fix requires global or per-atom energy or stress data. For efficiency reasons those calculations are skipped (and hence the flags zero) if not required.

I cannot comment on the SPIN package fixes currently, since I have not looked at their implementation properly before and don’t have the time to do it now (or soon).

The situation is the following:

  • if eflag_global is set, global energy is added to either Pair::eng_vdwl or Pair::eng_coul by calling Pair::ev_tally or one of its variants. Which variant is needed depends on the structure of the pair style and the nature of the model it implements. Pairwise additive interactions use the minimal version ev_tally().
  • similarly vflag_global triggers tallying the global stress, but there is a special wrinkle here: there are two ways to do that, and if the pair style supports it, the stress is computed by summing over the F dot r values per-atom rather than computing the stress components and tally them for each pair of atoms during ev_tally(), which is slower (sometimes quite noticeable).

You have to look at the sequence of steps here: 4.6. How a timestep works — LAMMPS documentation
Fixes can be called at multiple points in that flow depending on the functions they have implemented and the bits set in the Fix::setmask() function. If the order in which the calculations you want to do matters, then hybrid/overlay can be a problem. A fix can produce a global scalar or vector or a per-atom vector or array (there are Fix::compute_scalar() and Fix::compute_vector() for global data and Fix::compute_array() for per-atom data). Single fix can have a global scalar and either a vector or an array, but a global vector and a per-atom vector/array would be ambiguous (from the script interface).

I cannot comment in more detail, since I don’t have the time to work through your descriptions and the publications you are referring to. In particular not, because this is not my are of research, so I would need to do a lot more reading to update my background.

Hello again Axel,

The words and detail you’ve provided here have been more than valuable towards rectifying my issue. I will look into both extending my pair potential as a fix (as it makes the most sense to implement my body force term) and be more colonisent on how both eflag and vflag are utilized.

Thank you again for the help.

Hello again,

I’ve gone ahead and begun to reconstruct my pair_style as a fix per your recommendation below:

Do you know any fixes that request and use a neighborlist in this manner (looping over 1st nearest neighbors akin to a pair style) that I could use a reference to building my own? It would be greatly appreciated.

$ grep add_request src/*/fix*.cpp src/fix*.cpp
src/ATC/fix_atc.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL);
src/DIELECTRIC/fix_polarize_functional.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
src/DPD-REACT/fix_rx.cpp:  neighbor->add_request(this);
src/DPD-REACT/fix_shardlow.cpp:  neighbor->add_request(this, NeighConst::REQ_GHOST|NeighConst::REQ_NEWTON_ON|NeighConst::REQ_SSA);
src/ELECTRODE/fix_electrode_conp.cpp:    auto Req = neighbor->add_request(this);
src/ELECTRODE/fix_electrode_conp.cpp:    auto matReq = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
src/ELECTRODE/fix_electrode_conp.cpp:  auto vecReq = neighbor->add_request(this);
src/EXTRA-FIX/fix_electron_stopping.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
src/KOKKOS/fix_rx_kokkos.cpp:  auto request = neighbor->add_request(this);
src/MC/fix_bond_create.cpp:  neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
src/MC/fix_bond_swap.cpp:  neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
src/ORIENT/fix_orient_bcc.cpp:  neighbor->add_request(this,NeighConst::REQ_FULL);
src/ORIENT/fix_orient_eco.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL);
src/ORIENT/fix_orient_fcc.cpp:  neighbor->add_request(this,NeighConst::REQ_FULL);
src/PERI/fix_peri_neigh.cpp:  neighbor->add_request(this,NeighConst::REQ_FULL|NeighConst::REQ_OCCASIONAL);
src/QEQ/fix_qeq_dynamic.cpp:  neighbor->add_request(this);
src/QEQ/fix_qeq_fire.cpp:  neighbor->add_request(this);
src/QEQ/fix_qeq_point.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL);
src/QEQ/fix_qeq_shielded.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL);
src/QEQ/fix_qeq_slater.cpp:  neighbor->add_request(this, NeighConst::REQ_FULL);
src/REACTION/fix_bond_react.cpp:  neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
src/REAXFF/fix_qeq_reaxff.cpp:  neighbor->add_request(this, NeighConst::REQ_NEWTON_OFF);
src/REPLICA/fix_hyper_global.cpp:  neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);
src/REPLICA/fix_hyper_local.cpp:  auto req = neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL);
src/REPLICA/fix_hyper_local.cpp:  req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL);

Thank you!