Fastest way to sum pairwise potentials

Question

I've been writing a function to sum a pairwise potential on two lists, i.e. two charged bodies, each containing N >> 1 and M >> 1 atoms respectively.

I need to be able to calculate the potential of atom-atom, residue-residue or body-body potential.

In order to do so, I've defined the atom-atom, residue-residue, and body-body potentials as

Caa[atom1_,atom2_]:= Coulomb[atom1,atom2];
Crr[res1_,res2_]:= Total[Outer[Caa,res1,res2,1],2]
Cbb[body1_,body2_]:= Total[Outer[Crr,body1,body2,1],2]

where

body1 = {res11,res12,res13,...}
body2 = {res21,res22,res32,...}

res11 = {atom111,atom112,atom113,...}
res12 = {atom121,atom122,atom123,...}
...
res21 = {atom211,atom212,atom213,...}
res22 = {atom221,atom222,atom223,...}
...

atom111 = {{x111,y111,z111},q111}
atom112 = {{x112,y112,z112},q112}
...
atom211 = {{x211,y211,z211},q211}
atom212 = {{x212,y212,z212},q212}
...

i.e, two lists (bodies) of residues, which in turn are lists of atoms, all properly indexed.

My question is, what is the most efficient way to define Crrand Cbb?

I've tried with loops, Sum, Table, and a the Total[Outer[...],2] definition shown, but all seem to take almost the same (very long) time when doing a Timing check.

Needless to say, Coulomb is a radial function, i.e., it only depends on the distance between atom1and atom2, and atom = {{x,y,z},q}, where q is the charge.

--EDIT 1--

Coulomb[atom1_,atom2_]:= atom1[[2]] atom2[[2]]/Norm[atom1[[1]]-atom2[[1]]]

--EDIT 2--

In http://pastebin.com/Yf9TKSDx, you can find a small example of body1 and body2.

Here, body1[[i]], will give you the ith residue of the first body, body1[[i,j]], the jth atom of the ith residue of body1.

The atom-atom potential is calculated by doing

Caa[body1[[i,j]],body2[[k,l]]]

the residue-residue potential

Crr[body1[[i]],body2[[j]]]

and the body-body potential

Cbb[body1,body2]

hope this clarifies the problem.

Answer 1

Assuming your Coloumb function should be:

Coulomb[atom1_,atom2_]:= atom1[[2]] atom2[[2]] / Norm[atom1[[1]] - atom2[[1]]] 

then the calculation of Cbb[body1,body2] takes about 240 ms (on my PC) with your example data.

The best scope for speed-up is probably to compile Crr. Here I define a compiled function:

Crrcomp = Compile[{{res1x, _Real, 2}, {res1q, _Real, 1}, {res2x, _Real, 2}, {res2q, _Real, 1}},
  Plus @@ Flatten[Outer[Times, res1q, res2q]/
     Map[Sqrt[Dot[#.#]] &, Outer[Plus, res1x, -res2x, 1], {2}]]]

and redefine Crr to call this compiled function:

Crr[res1_, res2_] := Crrcomp[res1[[All, 1]], res1[[All, 2]], res2[[All, 1]], res2[[All, 2]]]

With this change, Cbb[body1,body2] now takes about 7 ms.

If you are using version 8 and have a suitable C compiler installed, you can add the options CompilationTarget -> "C" and RuntimeOptions -> "Speed" to the Compile function and this brings the timing down to about 4.5 ms.