1
$\begingroup$

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.

|improve this question
$\endgroup$
  • 1
    $\begingroup$ The Outer[] version is quite fast (almost idiomatic); you should maybe look into making Coulomb[] more efficient. $\endgroup$ – J. M. will be back soon♦ Aug 3 '12 at 6:40
  • 2
    $\begingroup$ Are you sure about this definition? Norm[v, p] is the $p$-norm of the vector v. Maybe you wanted EuclideanDistance[]? $\endgroup$ – J. M. will be back soon♦ Aug 3 '12 at 7:04
  • 6
    $\begingroup$ A tip: when you ask for help with things like this (fastest, cleanest, etc), strip your domain specific details from the problem. No one cares if it is Coulomb potential or something else. These details only make it harder for people to focus on your problem, which is why even a list-manipulation question hasn't received an answer in 7 hrs (usually you get 3-4 in the first hr). Assume some simple 2 argument function f that's close enough in time complexity to your Coulomb, and use res = RandomReal[...] (use a seed). Give us something concrete to objectively measure improvement $\endgroup$ – rm -rf♦ Aug 3 '12 at 14:01
  • 1
    $\begingroup$ @R.M Specific details were added because the form of Coulomb[] was suggested to be important for performance. Why speculate on complexity of Coulomb if I can simply post it. I'm a physicist and I think in terms of physics, and it seems unfair to say no one cares. I don't care if you don't care. $\endgroup$ – Pragabhava Aug 3 '12 at 20:40
  • 1
    $\begingroup$ @Manuel It's fine to post Coulomb's definition... I was asking you to strip all the body-body, residue-residue details which are only confusing. If you understand it, good for you! But keep in mind that you're seeking help from a Q&A site, where most people might not be physicists like you. I'm certainly not. In such cases, it is helpful if you kept the question solely on the list-manipulation part of it without other details. Moving on, am I right in assuming that your Caa, Crr and Cbb work as you expect (i.e., give you the right answer) and all you want to do now is improve speed? $\endgroup$ – rm -rf♦ Aug 3 '12 at 20:50
4
$\begingroup$

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.

|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for your response. I need some help understanding the code. In {res1x,_Real,2},{res1q,_Real,1}, what does the 1 and 2 stand for? $\endgroup$ – Pragabhava Aug 6 '12 at 17:26
  • $\begingroup$ J.M. pointed out that Outer was as fast as possible, therefore I'm accepting this reply as an answer (along with J.M.'s comment), due to significant speed improvement. Thanks Simon. $\endgroup$ – Pragabhava Aug 8 '12 at 2:49
  • $\begingroup$ @Manuel, the 1 and 2 are the ranks of the arrays. In this case, res1x is a 2D array of reals and res1q1 is a 1D array of reals. $\endgroup$ – Simon Woods Aug 8 '12 at 11:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.