I think I found a bug in compute group/group. Say I define two groups
of atoms A and B that sum to the total number of molecules. I then
use compute group/group for the different subsets of atoms:
compute AB A group/group B
compute BA B group/group A
compute AA A group/group A
compute BB B group/group B
compute All all group/group all
If I then sum the force as (AB + BA + AA + BB = All), the total
force=0 as it should:
variable dum equal (c_AB[1]+c_BA[1]+c_AA[1]+c_BB[1])
thermo_style custom c_All[1] v_dum
However, if I can only get the total energy to match if I use (AB +
AA + BB = All):
variable dum equal (0.5*c_AB+0.5*c_BA+c_AA+c_BB)
thermo_style custom pe c_All v_dum
This seems to suggest that the energy AB interactions are
double-counted in compute group/group.
Stan
Nice analysis - I'll take a look.
Steve
So imagine group A and group B are 2 simply
2 different single
atoms. If you use compute AB, then I think
you would expect to get the force between atoms
A and B, and the energy of the A-B interaction,
not half the energy.
So if you sum compute AB and compute BA,
you will get zero force (b/c it's signed)
and 2x the energy of the A-B interaction.
We could add the 1/2 factor to everything, but
I don't think that's what the user expects when
they just use compute AB, and for your combined
calculation it seems like just a bookkeeping thing.
Steve
Steve,
I see what you are saying. I've almost finished a kspace (Ewald/PPPM)
version of compute group/group but I'm not sure how to double the AB
interactions in kspace if you have some atoms only in one group and
some that are in both groups. It's easier if you separate the AB and
BA interactions explicitly. I'll send you some code soon so you can
take a look.
Stan
We could enforce an error condition with compute group/group
that the 2 sets of atoms either need to be the same group,
or disjoint, i.e. no atom in both groups. Doesn't seem logical
to use it when some atoms (but not all) are in both groups.
Steve
Steve,
I think I got it to work--I perform the group-group calcs again a
second time to get the BA energy (only). Enforcing the error
conditions would also work too.