See comparison of two results, one obtained with compute heat_fix and
another with this script on fig LJ_90K.jpg
on linux cluster using 8 processors.
It looks like it should work with any type of potentials.
Reference Schelling PR B65, 144306 - it's not important how you treat
three-body term.

However, and it's confusing me, both methods demonstrate dependence of
the number of processors (see Fig LJ.jpg)

Dear German,
Thanks for keeping this alive. I have a quick question. What is the cutoff you are using? For the case of 8 processors that the each processor has 3X3X3 box size. I am thinking that some interactions might be getting calculating twice (perhaps) due to the size constraints. Is it possible for you to run a larger box may be 12X12X12 for 1 and 8 processors and observe the differences.

So you've added this diagnostic to a run to compute the flux
yourself. Are you getting the same run on any number
of procs (independent of the diagnostic)? If so, and the
diagnostic just depends on quantities like per-atom eng
and virial (which also should not change as a function
of the number of procs), then I don't see how you would
get different answers as P changes.

So you've added this diagnostic to a run to compute the flux
yourself. Are you getting the same run on any number
of procs (independent of the diagnostic)? If so, and the
diagnostic just depends on quantities like per-atom eng
and virial (which also should not change as a function
of the number of procs),

Yes they are different for different number of processors.
I checked it on linux cluster - compiled with intel_mkl and openmpi,
and also on cray xt4 computer - compiled with "storm".
It has this dependence of number of processors.

Dear Vikas,
I've followed Stive's recommendation to add loop geom in velocity command and
started Equilibration in nve for 100000 iterations at the same volume.
Previously
I did Equilibration and thermalisation in NPT for 100000 iterations
before Equilibration.
Now results for different number of processors (fig LJ.jpg) and
different type of
heat flux calculations (fig LJ_90K.jpg) looks reasonable.

Dear German,
I think it should work. You can check it this way. Do couple of test simulations.

With no overlay and single pair 1 2 interaction

With overlay and both pairs of 1 2 interaction but make the appropriate coefficients equal to zero in the second pair (which you did not use in case 1).

With overlay and all coefficients non-zero

Case 1 and 2 should give you same answers while case 3 will be different than both. You can try to do this with both pairs seperately just to be sure. Let us know what you find if you do this.