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] += dzfpair;
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();