Compute group/group returns zero when called on the same group

LAMMPS LAMMPS Development
,
1

Hi,

I am trying to monitor the contribution of pair interactions between and within groups during a free energy calculation that uses ti/spring. In particular, I would like to understand if the pair interactions between a lambda-coupled group and the rest of the system are switched off during the simulation (the documentation for ti/spring doesn’t state explicitly if are only the pair interactions within groups to be switched off, or the ones between groups too).

I have used compute group/group to calculate the energy within a solid substrate, between the substrate and a surroundings liquid, and within the liquid:

compute  gg_11 substrate group/group substrate
compute  gg_22 liquid group/group liquid
compute  gg_12 substrate group/group liquid

I am using a EAM potential (pair_style eam/alloy); the liquid is molten Al, the solid is BCC W (T=1200K, constant volume).

Upon printing c_gg_ij[1] (here stored in pe_ij) I get always zero for intra-group interactions:

Step          Temp          PotEng        v_pe_11        v_pe_12        v_pe_22
52000   1224.561      -89537.731      0              4.834049       0
52100   1209.435      -89143.851      0             -14.194153      0
52200   1228.4934     -89122.292      0             -3.983723       0
52300   1215.5963     -88610.611      0              2.5011471      0
52400   1213.2597     -88849.873      0              6.1760416      0
52500   1213.8632     -89317.644      0             -9.5075325      0
52600   1209.2473     -88988.876      0             -15.825294      0
52700   1219.1943     -89140.501      0              2.2555602      0
52800   1202.2384     -89182.093      0              4.1723518      0
52900   1206.1197     -88746.281      0             -5.0029934      0
53000   1230.0787     -88926.291      0             -28.176244      0
53100   1211.95       -88362.79       0              6.1961017      0
53200   1207.41       -88747.516      0              21.171291      0
53300   1216.4681     -89341.206      0             -8.095013       0
53400   1215.7294     -89691.73       0              6.2598819      0
…

I am not sure about the behaviour of the 11 and 22 terms. The documentation for compute group/group states:

Define a computation that calculates the total energy and force interaction between two groups of atoms: the compute group and the specified group2. The two groups can be the same.

So I would not expect the intra-group energy to be trivially zero. Or is it? I am going to prevent probable follow-up questions: setting pair no doesn’t have any effect. Also, there’s no molecule definition and no long-range interactions requiring kspace calculations (it’s EAM).

I don’t think it relates to my problem, but I don’t fully get what this comment means:

EAM potentials will re-use previously(?) computed embedding term contributions, so the computed pairwise forces and energies are based on the whole system and not valid if particles have been moved since(?).

I guess this could matter in the molten phase, but still doesn’t explain the behaviour of gg_11.

2

Ok, I think I found the “problem”: It’s not a 4-dimensional vector (1 component for energy + 3 for forces), but a scalar + a 3 dimensional vector.

So this:

variable        pe_11 equal c_gg_11
variable        pe_22 equal c_gg_22
variable        pe_12 equal c_gg_12
thermo_style    custom step temp pe v_pe_11 v_pe_12 v_pe_22

works.

Sorry

3

Could you try running a test where you replace EAM with lj/cut?

Generally speaking, the compute group/group functionality assumes that interactions are pairwise-additive and thus you can loop over all atoms in a group and then compute interactions selectively looping over all pairs of atoms formed with the atoms in the second group.
Internally this is done by calling the single() function for each pair instead of compute() which loops over all pairs of neighbors in the case of a pair style more efficiently.

For many-body potentials this presents a problem, since those have not only contributions that are solely determine by the individual pairs, but also by larger tuples. Thus, almost all many-body do not support this single() function and thus not compute group/group. EAM is an exception because its computation is almost pair-wise. It is effectively a two step process, when in the first step you accumulate the amount of (electron) density from its neighbors at the point of each atom. This is then summed and communicated across the whole system and then you have a second loop over all atoms where you compute the pair-wise term of the EAM potential and the embedding term based on the accumulated density information. Since the pair style does not deleted that density information, you can later call the single() function and use the same information to recompute the interaction and you should obtain the same total force and energy if you loop over the full system as it is done with compute().

It would require some more detailed look and some reproducer inputs to determine the cause of the inconsistency or whether this is an inconsistency. Thus my request to test this with a proper pair-wise additive potential (you could also use Morse or something else, if that would be easier for you) to see whether to first look at how you set up the compute or whether the EAM pair style still does what the documentation says it does.

4

No worries. I am relieved that there is a simple resolution.
I was concerned that something rather subtle had been broken (it has happened before).

5

Great! And thanks for the explanation.