compute intermolecular energy with compute command

Andrew,

Attached is the experimental kspace compute group/group code, which works for both Ewald and PPPM. I just updated it to compile with the latest version of LAMMPS, though it may still break some of the USER add-on packages. I also included a brief 3 page summary of the equations involved (really this is just an extension of the same method used to get per-atom energy for Ewald/PPPM). For Ewald, the compute is very efficient because the needed per-atom structure factors are already calculated and stored previously in Ewald::compute(), so it’s just a matter of summing them over atoms in the different groups. For PPPM, the compute requires 2 extra forward FFTs each time (one for each group), which is about half the cost of a regular PPPM::compute(), which is also not too expensive. I haven’t extended the compute to the Ewald/n package yet but it could easily be done. I’d welcome any feedback.

Thanks,

Stan

P. S. Here is a summary of the tests I have performed so far for SPC/E. I created two groups, A and B, that correspond to the oxygen and hydrogen atoms respectively, as well as a group of a single atom:

group A type 1 # (oxygen atoms)
group B type 2 # (hydrogen atoms)
group 1atom id 1 # (one atom)

compute AB A group/group B
compute BA B group/group A
compute All all group/group all
compute AA A group/group A
compute BB B group/group B
compute one 1atom group/group all

I then tested that all of the following conditions were met:

Energy:
All = pe
AB + AA + BB = All = pe
AB = BA
one = per-atom energy for that atom

Force (for all 3 components):
All = 0
AB + BA + AA + BB = All = 0
AB = - BA
one = per-atom force for that atom

kspace_group_group.tar.gz (158 KB)

Just added this as a 21May12 patch with some options
in the compute group/group command to control whether
the pair and/or kspace contribution gets added in.

This is a nice piece of work by Stan. I didn't know it was
possible to do these calculations in a group-wise manner.

Stan - I added your eqs PDF to the doc pages as well
with links from compute group/group and compute pe/atom.
I think you should add a brief intro paragraph to the doc
and put your name/email on it, so people will know who
to credit.

Axel, Mike, Christian - you can look at Ewald and PPPM to
see see if your Kspace variants can support this. If you add
it, just set a group_group_enable flag = 1 at the top of your
constructors, else it should be off by default.

Steve