Thank you for porting the MOLC package to the latest dev, it seems correct to me.
The method you suggest also works fine, and does not require any new compute. To summarise, the inter-molecular (aka cohesion) energy for a molecular system with long-range electrostatics can be obtained from the equation:
E_\text{inter}=E_\text{pair}-E_\text{intra}
where E_\text{pair}=\sum E_\text{vdwl}+E_\text{coul}+E_\text{long} and E_\text{intra} is computed as follows:
compute intra all group/group all pair yes kspace yes molecule intra
variable enb equal evdwl+ecoul+elong
variable inter equal v_enb-c_intra
compute inter2 all inter kspace yes # the new compute.
compute inter3 all group/group all pair yes kspace yes molecule inter
The intermolecular energy computed in this way is numerically identical to that computed with the new compute inter command:
v_inter c_inter2 c_inter3 v_enb c_intra
-12353.99 -12353.937 -21524.214 -8025.8405 4328.1496
To me this result means that the intermolecular energy computed with computed group/group is not correct, as it is equivalent to a rerun with neigh_modify exclude molecule/intra all. Another way to see that something does not add up is by summing c_intra+c_inter3, which does not return v_enb.
@stamoor, can we take this offline? Or, better, in the Open Review? I still need to publish this stuff, and your revision is definitely the best I can hope for.