Subject: Summary: new square root improved molecular dynamics

From: John Marelius <[email-masked]>
Date: Fri, 6 Sep 1996 21:41:30 +0200 (MDT)

This posting summarises the results of the discussion on the usenet
newsgroup comp.sys.intel on square-root calculations in molecular dynamics
(MD) applications. Using a novel root algorithm and general optimisations
of the MD code, performance has increased by up to 55%, depending on the
simulated system.

Background:
The computational cost in molecular dynamics simulations (with
water represented as individual molecules) is dominated by
calculations of non-bonded (long-range) interaction energies and forces
between pairs of atoms. In particular, the cost scales with the square of
the number of atoms. The energy function has the following form:

         A(i,j)   B(i,j)   q(i)*q(j)
E(i,j) = ------ - ------ + --------- 
          r^12     r^6         r
where r is the distance between atoms i and j, A and B are constants
depending on the kinds of atoms and the q:s are the (partial) electric
charges of the atoms.

If it wasn't for Coulomb, who introduced an odd power of r (r^-1) into his
law, this could be evaluatued without a square root calculation. (In
cartesian space r^2 is easily obtained from the co-ordinates.)

The core of the problem is that square root calculations and also
divisions take up to two orders of magnitude longer than multiplications 
and additions to execute on most computers.

One could of course consider working in another co-ordinate
system but, unfortunately, other aspects of the potential
energy fuctions (such as bond angles) then become more demanding.

Attacking the problem:
Knowing the above, plus the fact that the MD program I use spends nearly
50% of its time with calculations of sqrt(1/r^2), I decided to try to
improve the overall speed by finding a faster way to calculate square
roots.

The initial discussion on comp.sys.intel and comp.lang.asm.x86 taught me
that some types of CPUs have floating-point unit (FPU) instructions
for square root, while others have not (a function library call is
necessary). The speed difference between the two ways is not as
big as I first expected, and most modern processors can do the square root
in  hardware.

Approaches by me and others involving integer-math, look-up tables,
Taylor polynominals or Newton-Raphson iterations ended up being
either slower or less accurate than the standard function.

One person who became interested in the problem was Terje Mathisen, he
started thinking and coding and has overwhelmed me with a stream of
ever-improving codes for fast and accurate functions to calculate
sqrt(1/x).

Please refer to his posting (coming soon), describing the details of the
algorithm! 

Most modern processors have more than one execution unit for
floating-point calculations (parallellism) and each of these can start
working on new instructions before previous instructions are finished
(pipelining). This opens for another type of optimisation, namely
rewriting the progam so that it can make maximum benefit from this.

Results:
I incorporated the new sqrt(1/x) algorithm in the MD program. I also
modified the critical sub-routines to do two or three independent
calculations in parallell to avoid 'pipeline stalls'.

To benchmark the new program it was tested in very short simulations of 
two typical systems:
1) A system dominated by protein (ca 2000 atoms + ca 1000 water atoms
(=ca 333 water molecules)) 100 dynamics steps were done.
2) A system dominated by water (one 12-atom solute molecule + 1200 water
atoms), where the non-bonded interactions are even more dominant in terms
of CPU time. 250 dynamics steps were done.
The performace tests of the new program are indeed very
encouraging!
The tests were run on a 400 Mhz DEC AlphaStation 500.

			system 1		system 2
			time	rel. speed	time	relative speed
			(s)	(%)		(s)	(%)
original program	37.0	100		37.4	100	
new program		30.2	123		24.1	155

Initiation overheads have been subtracted.

Conclusion:
My week-long simulation jobs will now complete in 4.5 to 5.7 days instead!

I am very grateful to Terje Mathisen for the new square root algorithm.
I also acknowledge the many people all over the internet who contributed
valuable information.

John Marelius
graduate student in molecular biotechnology

+-----------------------------------------------------------------------------+
| 	   John Marelius                                                      |
| 	   Dept. of Molecular Biology, Uppsala University                     |
| E-mail:  [email-masked]   						      |
| www:     http://www.ibg.uu.se/~john					      |
| fax:     +46 18 536971						      |
| address: Box 590, S-751 24 Uppsala, Sweden	     			      |
+-----------------------------------------------------------------------------+