Fix between pair computations?

Hello,

I am writing a custom fix. From what I understand, the sequence in verlet.cpp requires fixes to occur pre-force or post-force. However, my simulation requires the following sequence:

-Fix function 1
-Pair styles 1 and 2, which depend on fix function 1
-Fix function 2, which depends on pair styles 1 and 2
-Pair styles 3 and 4, which depend on fix function 2

I understand that some options could include changing verlet.cpp or moving functions from the pair style to the fix. However, I was wondering if there are any pointers on how to achieve this alternating pair and fix sequence in a less hackish manner.

I believe the fix could call a function from the pair style or the pair style could call a function from the fix, but I am 1. Not sure of the syntax and 2. Not sure if running things in this order woud cause complications with the energies and forces computed by the involved pair styles.

Thank you,
Anne

you have to explain a bit more about what you want to achieve and why you need this flow of control.

also, which parts of the computation workflow are standard LAMMPS styles and which would be custom? are all those pair styles part of a hybrid pair style?

axel.

Hi Axel,

I am using hybrid/overlay with a mix of standard and custom pair styles. The need for this sequence is mostly due to modification of and dependence on charge values as follows:

  • Fix function 1 modifies some of the charges

  • Pair styles 1 and 2 are custom and depend on these new charges. They output a small subset of short range energies.

  • Fix function 2 depends on both charges from fix function 1 and short range energies from pair styles 1-2. It updates charges on a different set of atoms from fix function 1.

  • Pair styles 3-4 are standard long-range coul (depends on all updated charges) and lj (independent).

Thank you,
Anne

Also, the custom pair styles still contribute to the forces on the atoms.
Thank you,
Anne

the way i see this can work would be by implementing both fixes and pair styles 1 and 2 as a single fix and run it at the PRE_FORCE stage (i.e. right after the forces are cleared but before any of the standard pair styles.

you can request perpetual neighbor lists also from fixes and then just do the whole thing as multiple steps.
you will need to do multiple forward communications (to send the modified charges from local atoms to ghost atoms) and reverse communication (to collect forces on ghost atoms to their local atoms) and clear out the forces on the ghost atoms after the reverse communication for that to work properly.

fixes can contribute to forces and virial. but you may need to do the latter manually (if they can use F dot r) or tally them with individual forces.

another approach would be a custom time integration style, but that could require complex changes to much of the inner workings of LAMMPS.

axel.

Hi Axel,
Thank you very much for the feedback. I will try that out.
Anne

Hello, just following up on this question. I was able to figure out neighbor list operations.

However, I am stuck on how to call ev_tally and virial_fdotr_compute from the fix instead of the pair. Is there a straightforward way to do this?

Thank you,
Anne

Hello, just following up on this question. I was able to figure out neighbor list operations.

However, I am stuck on how to call ev_tally and virial_fdotr_compute from the fix instead of the pair. Is there a straightforward way to do this?

there is an ev_init() / ev_tally() API in the fix base class, but it works differently from pair styles. I suspect, It doesn’t even support delaying the global virial computation through F dot r instead of tallying it for each individual pairwise force computation. You need to read the source to be certain.
Most importantly, It does not add contributions to the pair style accumulators, but there is a separate accumulator for fixes.

You have to be very careful with the F dot r computation, because you can only use it if no other forces have yet been added, since this computation needs to be done by looping across local and ghost atoms and thus before forces are collected with a reverse communication (in the default setting of force->newton_pair == 1). This may have implications for the pair styles you use as well; they may not be allowed to use the F dot r virial computation, if the forces are already “tainted”.

Axel.

Hi Axel,

Thank you for the pointers. To work around what was mentioned, my idea to implement pair_style functionality in a pre_force fix is shown as follows. I wanted to check if there are any issues that I may be overlooking if I take this approach:

  • Create this custom pre_force fix, along with a corresponding custom pair style. The fix accesses the pair style using a pointer as shown below.

  • For energy accumulation/tallying, the fix uses the pair style’s variable (eng_vdwl) and function (ev_tally)

  • Similarly, for the virial, the fix uses the pair style’s function (virial_fdotr_compute)

  • To not interfere with virial computation, the fix stores computed forces on special f_tmp arrays local to the fix. (Relevant forward/reverse communcation to handle the local/ghost atoms is done during the preparation of f_tmp.)

  • When the custom pair style runs later, it extracts f_tmp from the custom fix and adds to the regular force arrays. (Forces are not urgently needed at the pre_force stage). No energy or virial functions are performed in the pair style, since they were performed at the pre_force fix stage already.

Relevant code excerpts from the custom fix:

// Get pointer to access pair style functionality from the fix
int inum,jnum;
int *ilist,*jlist;
int *numneigh;
int **firstneigh;
int nsub,itmp;
for(nsub=0;nsub<3;nsub++){
paircustom = (Pair *) force->pair_match(“pair/custom”,1,nsub);
if(paircustom){
inum = paircustom->list->inum;
ilist = paircustom->list->ilist;
numneigh = paircustom->list-> numneigh;
firstneigh = paircustom->list-> firstneigh;
break;
}
}
if(!paircustom) error->all(FLERR,“Pair style not found”);

[… pair style computations and loop structure … ]

// save force computations for the pair to arrays local to the fix
// (perform relevant fwd/rev communications and add an extract function)

f_tmp[i][0] += dx*fpair;

f_tmp[i][1] += dyfpair;
f_tmp[i][2] += dz
fpair;
if (j < nlocal) {
f_tmp[j][0] -= dx*fpair;

f_tmp[j][1] -= dy*fpair;

f_tmp[j][2] -= dz*fpair;
}

// Energy commands for the pair

// Typical energy commands from pair_style implemented in fix
if (eflag) evdwl = [ energy computed for pair];

paircustom->eng_vdwl += evdwl;

if (evflag)

paircustom->ev_tally(i,j,nlocal,newton_pair,

evdwl,0.0,fpair,dx,dy,dz);

// Call virial function for the pair

// (Forces are saved to f_tmp at this stage to not intefere with virial calculation)

if (vflag_fdotr) paircustom->virial_fdotr_compute();

that seems overly complex and I am not convinced it will work as expected. It would be cleaner to just copy the body of your custom pair styles into the fix and adapt that. But if you want to go down this route, I would start implementing the whole thing as a “mega” pair style and take inspiration from pair style hybrid/overlay. That one shows you the logic and effort required to maintain and run multiple pair styles from within some class and how to make certain all settings are consistent.

I also would not want to mess with even trying to get the virial from F dot r but explicitly disable this (after creating the sub-pair styles, you can just set no_virial_fdotr_compute to 1 for each instance) and then the virial is correctly computed for global and per-atom contribution via ev_tally(). whether F dot r really has a significant performance impact is something that I would look into after everything works. you can selectively decide which data to communicate during forward and reverse communication by using a flag (see pair style eam/cd for an example).

Axel

Thanks Axel,

Copying the pair style functionality to the fix is basically what I am doing. I will take your advice on the viral_fdotr_compute, and therefore I will handle force accumulation in the fix now.

The only thing remaining is how to handle eng_vdwl and ev_tally. Would the previously mentioned approach of calling them using a paircustom-> pointer work?

As long as this does not present issues, I feel it would be not too complex.

Thank you,
Anne

I am also considering the mega pair style idea. Would it allow for running a fix in between pair styles? (This was my original issue)

if you compute/accumulate energies and virial contributions inside they fix, they should be accumulated into the fix accumulators.
the accumulation is a bit more complex as you need to factor in that you have either global or per-atom energies and stress contributions and that the global stress contribution of a pair style may or may not be accumulated via ev_tally depending on whether the pair style support F dot r and whether per-atom stress is requested.

Axel.

I am also considering the mega pair style idea. Would it allow for running a fix in between pair styles? (This was my original issue)

why does it have to be a fix? pair styles like eam have multiple “stages” and are thus similar in spirit to what I understand you want to do. pair styles that need to be used in combination with fixes (e.g. reax/c, comb, comb3, coul/streitz) usually do the fix operation either completely before or completely after the entire force computation.

my recommendation is to either do everything from inside a pair style (and then have sub-pair styles and do the “fix work” in between) or do the entire work inside the fix until you can call the existing conventional pair styles.

if you want to alternate between pair computation and fix calls, you need to write your own integration style and have one or more classes like the Force class that manage the pair (and other) styles computing forces. we’ve been over this before.

at this point, i am getting tired of offering the same advice multiple times.

bottom line: LAMMPS is open source, so you are free to implement and change things any which way you like, but at the same time, you cannot expect much enthusiasm to help you, if you ask for advice and then decide to go against it. Also please keep in mind that I don’t have a crystal ball and don’t know many details of what you intend to be doing, so I cannot give very specific advice and definitely no guarantees about what is going to work (well) and what not and what would be the best approach.

Axel.

Hi Axel, I better understand the mega pair style idea now and it sounds like a more straightforward route as opposed to the fix. Thank you for your time in assisting with this.
Anne