**Particles On a Sphere
— A Specimen Multidisciplinary Problem**^{*}

Timothy J. Rolfe

Associate Professor, Computer Science

Dakota State University

Madison SD 57042-1799

(605) 256-5166

Electronic mail: `rolfe@alpha.dsu.edu`

WWW home page: `http://www.dsu.edu/~rolfe/`

__Current information__

Timothy J. Rolfe

Professor, Computer Science

Eastern Washington University

Cheney, WA 99004-2412

(509) 359-6162

Electronic mail: `TRolfe@ewu.edu`

WWW home page: `http://penguin.ewu.edu/~trolfe/`

^{*} Copyright 1996 by the Small College Computing
Symposium

Permission to make printed or digital copies of all or part of this material for educational or personal use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies include this notice on the first page.

Appears in **SCCS Proceedings: 29th Annual Small College
Computing Symposium** (SCCS, 1996), pp. 456-466.

- Abstract
- Problem Description
- Program Design
- Results
- Summary
- Availability of Source Code
- Acknowledgements
- Notes

Click here for the paper in Microsoft Word format.

Computer professionals require knowledge of fields outside of their narrowly defined field. The example developed here shows how a problem within Chemistry (maximum separation of ligands around a central atom) requires for its solution not just background in Chemistry and programming, but also significant mathematical background as well as Physics (electrostatics and classical mechanics), numerical analysis (for solutions of the differential equations), and computer graphics.

Computing ultimately cannot operate in a vacuum: the computer must be used for some problem, generally coming from outside of Computer Science. Often, however, the problems addressed end up requiring background in a number of fields rather just in programming and the application area. A specimen problem taken from Chemistry ends up requiring in addition classical mechanics, graphics, mathematics, and numerical analysis.

Within Chemistry, the structure of a molecule can be influenced by steric hindrance — the physical packing of components around the atom to which they are mutually bound. For instance, the preferred configuration for cyclohexane is the trans or "chair" configuration rather than the cis or "boat" configuration. [1] So it is of interest to discover the arrangement that will have all components as far apart as possible.

Besides steric hindrance, another method of estimating bond angles within a molecule is to model the electrons participating in the bonding as mutually repulsive, so that the bonds are positioned with maximum separation around the atom. This can be treated as the problem of positioning particles on a sphere so as to maximize the separation of the particles. [2]

Based on steric hindrance, the problem of calculating the
positions could be approached purely mathematically, as an *N*-parameter
maximization problem. As such it presents a significant challenge,
either for a purely mathematical approach or as a minimization
problem using successive approximations: for** ***n***
**ligands, there are 3 ** n** Cartesian
coordinates, of which we are allowed to discard 6 as independent
of the interior of the system itself (the center of mass position
and the system's orientation in space). [3] The problem can,
however, be approached in a multidisciplinary fashion as a
computer simulation of a system operating under the constraints
of classical mechanics.

The model used here comes directly from the bonding electron model: the particles are treated as mutually repulsive, where the forces acting are the electrostatic repulsions (no gravitational forces), subject to the constraint that the particles remain on the surface of a sphere. Within Physics, classical mechanics provides the formulas for calculating the forces acting on each of the particles. Within Computer Science, the field of numerical analysis provides the tools for the simulation as the solution of the differential equations coming from the Physics.

The problem of finding the *maximum* separation
translates into the problem of finding the minimum in the *N*-dimensional
potential energy surface. While the Physics will have the
particles moving towards the potential energy minimum, the
particles also acquire kinetic energy; if the Physics are
correctly implemented, the total system is conservative, with
energy exchanged between potential and kinetic forms. The kinetic
energy can be bled away in two ways: (A) adding a frictionlike
force, or (B) intermittently freezing the system (removal of
all motion).

The user interface also provides some challenges in graphics, mechanics and computing. There is the obvious problem of displaying the propagating system, but there is also the problem of finding and implementing the viewing position appropriate for displaying the symmetries of the final system when the potential energy minimum has been found. The symmetries themselves provide a tie in with the mathematical area of Group Theory and are also central to the area of Inorganic Chemistry.

The first consideration is the coordinate system to use for the model. While we might view our problem as positioning electrons about an atom, we restrict the problem of positioning points of equal mass all at the same distance from the center. As such, we are positioning points on the surface of a sphere, and we choose a simplified coordinate system, [4] setting the sphere radius and the particle mass to unity; we also treat the particles as all being of equal charge. Electrostatic repulsion (effectively, the magnitude of that charge) will be determined by a force factor entered by the program user.

The next consideration is which model of Newtonian Physics to
use in setting up the differential equations. Given ** n**
particles within three-dimensional Cartesian space, one can view
the system of differential equations as a system of 3

Using the Hamiltonian formulation, we represent each point in the developing system by six coordinates: three Cartesian coordinates for position and three further Cartesian coordinates for momentum (since in the simplified system the particle mass is unity, the momentum and velocity components are identical). The momentum components provide the first derivative information for the positions, and the electrostatic forces provide the first derivative information for the momenta.

The electrostatic forces are determined by Coulomb's Law: given points with identical charge, we can collect that charge information with other constants involved in Coulomb's Law into a constant K to generate the following equation: [7]

*F*_{1} = (K /[ *r *_{12}]^{2}) *e*_{12}

— where *F*_{1} is the (vector) force
acting on particle 1, [*r *_{12}]^{2} is the (scalar) square of the
distance between particle 1 and particle 2, and *e*_{12}
is the unit vector from particle 2 to particle 1. The force
acting on particle 2 is just the negative: *F*_{2}
= -*F*_{1 }.

The total force acting on a particle in the system is the sum
of all of the pair-wise forces with the other particles in the
system. For programming purposes, we can take advantage of the
fact that *F*_{2} = -*F*_{1}:
the force calculation will involve a nested loop:

for (p1 = 0; p1 < Size-1, p1++) for (p2 = p1 + 1; p2 < Size; p2++) { Body of the calculation, summing force components on the particles. }

A more difficult problem is the constraint that the particles remain on the surface of a unit sphere. Without such a constraint our collection of charges will fly off towards infinity. The way we accomplish this is to force the total force vector to remain perpendicular to the radius vector from the center to the particle; that is, for the force to act only along the sphere's surface. In addition, after a time step each particle has its Cartesian position adjusted back onto the unit sphere (since even an infinitesimal step on the surface tangent to the sphere moves away from the spherical surface).

Once we calculate the necessary derivatives, we are ready to
solve the differential equations: the state of the system at time*
T*+*h* depends upon the state of the system at time* T*
(and possibly earlier time points). [8, 9, 10]

Based on the forward differences approximation to the first
derivative one can obtain the iteration equation *f*_{j+1}
= *f*_{j} + *h* (d*f*/dx)_{j }and
the technique is called Euler's method. Based on the central
differences approximation to the first derivative one can obtain
the iteration equation *f*_{j+1} = *f*_{j-1}
+ 2*h* (d*f*/dx)_{j }and the technique is
called Milne's method. [11]

The program as developed uses Milne's method; hence it is
necessary to retain both the present system configuration and the
previous one. In addition there is a start-up problem: the
initial system does *not* have a history, so that there is
no configuration for* T* = -*h*. The answer
is to generate one by back-propagation based on Euler's method: *f*_{j1 }= *f*_{j}
- *h* (d*f*/dx)_{j}.

At this point it is time to face implementation issues for
constraining the particles to remain on the unit sphere. While
one could use cross products to determine the tangent surface at
an arbitrary point on the sphere and to discard force components
normal to that surface, it is appreciably easier to rotate every
point onto the *z* axis. Once each point takes its time step,
it is rotated back to its original position.

To move a point onto the *z* axis we will need just two
of the three Eulerian angles for performing an arbitrary rotation
of a body in space. [12, 13] Specifically, we will need the two
angles for the transformation from Cartesian into spherical polar
coordinates. Within the physical sciences (using a convention
exactly opposite that of mathematicians), the angle theta is the
angle from the *z* axis to the point (and can have values
between 0 and 180). The angle phi is the angle from the *x*
axis to the projection of the point onto the *xy* plane (and
can have values between 0 and 360). Hence the point can be
rotated to the *z* axis by a rotation of -phi about the *z*
axis (that is, into the *xz* plane), followed by a rotation
of -theta about the *y* axis. Once the point (position and
momentum information for the present and previous configurations
along with the derivatives) has been rotate onto the *z* axis,
removal of force components off the unit sphere is simply a
matter of zeroing out the *z* components. After the step has
been taken, the constraint of remaining on the unit sphere can be
handled by replacing the position vector with the unit vector in
that direction. The point is moved back to its position by a
rotation of +theta about the *y* axis, followed by a
rotation of +phi about the *z* axis.

Newtonian Physics provides a validity check on system propagation: if the forces acting on the particles are correctly calculated, and the time step is correctly implemented, the total energy of the system will remain constant, with energy flowing back and forth between kinetic and potential energy [14] — unless a pathway is provided to bleed energy from the total system.

Since we are interested in finding the maximum separation, we need to provide a pathway for removing kinetic energy from the system. The program will then find the minimum in potential energy, the state of the system with the particles at maximum separation.

One pathway is provided by adding a friction-like force: a force exactly opposed to the momentum vector and proportional to the particle's speed. [15] Implementation in the calculation of the forces is straight-forward, given momentum components and the constant of proportionality (a friction coefficient).

A second pathway is provided by totally removing kinetic
energy, effectively suddenly freezing the system. The optimal
time for freezing the system will be at the point at which energy
begins to flow from kinetic energy back into potential energy.
The only tricky point is that a fictitious history is needed for
the system at time *T*-*h* consistent with zero momenta
at time* T*; the same back-propagation used for initializing
the system can also be used here. (Otherwise the Milne method
used for propagating the total system will generate anomalous
results due to the inconsistency of the system at time* T*-*h*
and time* T*.)

Based on a "simplest first" design philosophy, the
initial system display module uses ANSI character graphics on a
24-row by 80-column video screen (viewing the projection of the
system onto the *yz* plane, using *x* axis
information for brightness to distinguish points above the plane
from points below it). This allows portability across computers
without porting a complex graphical display library. At a later
point, the display module can be customized to a more accurate
and pleasing form specific to the target machine. (By the time of
the presentation at this conference, there will be a generation
of the program specific to the Borland C++ environment for the
IBM PC family displaying the points as circles of varying color
and size.) For display of a trajectory not in real time, it is a
simple matter to collect coordinate information at individual
steps and then use a more complex graphical rendering system.

Once we have a display module to examine the system, we will
probably want to rotate the total system the better to view
symmetries in the system. [16] Since our force calculation
requires modules for rotations about the *z* axis and the *y*
axis, the complete suite only requires one more module for
rotation about the *x* axis. The logic has already been
discussed for moving a given point onto the *z* axis. In
addition, one can use three points to determine a plane and then
set that plane perpendicular to the *z* axis: the cross
product of two difference vectors based on those three points
will generate a vector normal to that surface; hence one can use
earlier logic to move that surface normal vector onto the *z* axis.

The initial configuration for a trajectory can be specified
within a file or it can be generated randomly. Space
considerations preclude discussing the hazards of randomly
choosing angular spherical polar coordinates for the distribution.
The technique used in its place is to randomly choose points
within the cubical box with *x*, *y*, and *z* each
in the open range (1, +1); triples found to lie outside of
the unit sphere are rejected. The accepted triple is then
transformed into a unit vector, giving a point on the unit sphere.

Once an initial configuration is obtained, one begins the trajectory by specifying the duration (through the final value for time), the granularity (through the step size), and the magnitude of the repulsive force between particles. (For the friction pathway, the coefficient of friction is also entered here.)

The trajectory will be in general a complex interaction among
all of the particles in the simulation. Consequently it would be
highly impractical to attempt to reproduce an entire trajectory
on paper. We can, however, show the __start__ of a trajectory.

The specimen will be a random configuration of four points,
initially at rest on the unit sphere and interacting under an
electrostatic repulsive force. In Figure 1, the initial position
is shown by a symbol O, where the size
of the symbol represents position information relative to the
plane of the paper. (Specifically, the page represents the *yz* plane,
with the *x* axis normal to the plane of the paper. The
position of the symbol represents the projection onto the *yz* plane
with perspective approximated by the size of the symbol.) The
final position is shown by a symbol X;
again, the size of the symbol represents position relative to the
plane of the paper. Within the arbitrary units of the simulation,
the integration time step (*h*) is 0.001, the force constant
is 0.005, and the time lapse is from *T*=0 to *T*=5,
and then from *T*=5 to *T*=10.

*T* = 0 to *T* = 5

*T* = 5 to *T* = 10

**Figure 1.**

Since we are attempting to represent three dimensions in the
plane, the instances in the *T* = 0 to *T* =
5 display which appear to have the O
and X overlapping do not indicate
unmoving points, but rather ones that move around the sphere. In
T = 5 to T = 10 display we see appreciable position changes as
velocities build under the repulsive forces. In this step the two
paths that appear to cross probably cross opposite sides of the
sphere since they actually end on opposite sides. The line drawn
is NOT the path taken, just a visual queue to the initial and
final positions for a given point.

For the trajectories discussed below, we will show only the single initial configuration and then several final configurations (using filled circles whose sizes represent position relative to the plane of the paper).

Figure 2 shows the initial and final configurations for a
random run with six particles from *T*=0 to *T*=50, in
which there is no pathway to drain kinetic energy from the system.
The same initial configuration will also be used for specimen
runs with pathways draining kinetic energy.

Initial configuration used for N=6 runs

Final configuration with no kinetic energy drain

**Figure 2**

Figure 3 shows the behavior of (from top to bottom in the plot) total energy, potential energy, and kinetic energy.

**Figure 3**

Total energy remains appreciable constant throughout the course of the trajectory, suggesting that we have indeed correctly implemented the Physics for this system.

The next test is to generate a trajectory with a pathway for
kinetic energy loss. Testing first the friction pathway, from the
common initial configuration we obtain the final configuration in
Figure 4 (shown both as obtained directly from the run and also
as rotated so that one point is forced to the *z* axis
and another point is forced to lie in the *yz* plane).

Final configuration with friction as kinetic
energy drain

Final friction configuration, rotated

**Figure 4**

Figure 5 shows the behavior of total, potential, and kinetic energies.

**Figure 5**

The final configuration is moving toward the minimum energy configuration, with the particles at the vertices of an octahedron; that is, with them sitting on the three Cartesian axes at ±1. There are four points on an equator in the rotated view, with two points moving towards the poles.

The final specimen trajectory displays the results of the alternative pathway to remove kinetic energy: at any trajectory step at which energy moves from kinetic energy back into potential energy (that is, when the trajectory has passed a local minimum in potential energy), remove all kinetic energy by removing all velocity. Figure 6 shows again the final configuration (starting from the common initial configuration) and that configuration rotated to ease viewing the symmetries.

Final configuration with momentum freezing as
kinetic energy drain

Final freeze configuration, rotated

**Figure 6**

Behavior of total energy, potential energy, and kinetic energy is shown in Figure 7.

**Figure 7**

The final results are appreciably closer to the octahedral
configuration. The final configuration at *T*=50 is
approximately looking at an octahedron down a vector normal to
one of the triangular faces of the octahedron. The rotated
configuration is approximately looking down a vector from one
vertex to the opposite vertex (the far vertex is obscured by the
near one).

If either trajectory with energy drain is allowed to continue further, the points do settle exactly onto the expected configuration that can be rotated to have points at ±1 on each of the Cartesian axes.

The presentation of this paper at the SCCS conference will
include animation of a trajectory of five particles passing
through a metastable configuration on the way to the final
minimum: specifically, through placement at the vertices of a
pyramid with four faces and a base on the way to the trigonal
bipyramid — two pyramids of three faces and base joined at the
base. Figure 8 shows the *T*=100 and *T*=200
configurations, slightly rotated so that no foreground point
hides a background point. (In the display below, the bases of the
two pyramids are approximately in the plane of the paper; one
projects above the paper; the other, below.)

*T* = 100

*T* = 200

**Figure 8**

As Floyd Johnson of Northwestern College in Orange City, IA, commented in an electronic mail exchange about this paper, [17] Computer Science ends up tying into every other subject, even municipal garbage collection.

This paper has discussed one particular application in which a computing solution to a problem in one field has required knowledge of a number of other fields. This interrelatedness of disciplines is not peculiar to computing within the physical sciences. Computing the solution to any problem of significant size will probably also require pulling in ideas for a number of other disciplines, whether it be the propagation of galaxies in astrophysics, or word counts and stylistic analysis in a study of the authorship of the Pauline Epistles in the New Testament.

Programs developed by the author for this project are
available through the"Programs"
link on author's WorldWide Web home page: `http://penguin.ewu.edu/~rolfe/`.
Direct link to the code page

I wish to thank my colleagues Dr. Brian Carlson and Dr. Mary P. Freier of Dakota State University for fruitful discussions about this project, as well as Floyd H. Johnson for fruitful discussions in the electronic mail exchange mentioned above. I especially wish to thank Dr. Freier for her offer to provide an animation (to be shown at the SCCS conference itself) of a specimen trajectory using some of the graphical animation tools available.

I first encountered this problem and undertook it as a recreational diversion in 1981 while I was a graduate student in theoretical physical chemistry under Stuart A. Rice at the University of Chicago. Through the obscuring mists of the elapsed 15 years I would like to thank fellow graduate students Dr. Alan Belch, Dr. R. Michael Townsend, and Dr. Ralph Weber for fruitful discussions back then, and I am grateful to Stuart Rice for the idle cycles on his research group's PDP-11/34 used for that diversion. For the purposes of developing this paper, though, I have developed the program from scratch without making appreciable reference to the Fortran version from 1981-2.

- R. Stephen Berry, Stuart A. Rice, and John Ross,
*Physical Chemistry*(John Wiley & Sons, 1980), p. 381. - Claxton and Benson ("Stereochemistry and Seven
Coordination,"
**Canadian Journal of Chemistry***44*(1966), 157-63) include a reference to a paper on the arragements of particles interacting with a 1/4 interparticle potential: L. Foeppel,**J. Reine Angew. Math.***141*(1912), 251. - E. Bright Wilson, Jr., J. C. Decius, and Paul C. Cross,
*Molecular Vibrations*(Dover Publications, 1980: unabridged and corrected re-publication of the 1955 publication by McGraw-Hill), pp. 6-27, 54 ff. - An example of such a simplified coordinate system is found in "atomic units," as discussed in Berry, Rice, and Ross, p. 180.
- Herbert Goldstein,
*Classical Mechanics*(2nd edition; Addison-Wesley, 1980), Chapter 1-2, pp. 1-69. - Goldstein, Chapter 8, pp. 339-77.
- Richard P. Feynman,
*The Feynman Lectures on Physics, Vol. II: Mainly Electromagnetism and Matter*(Addison-Wesley, 1964), pp. 4-2 ff. - Richard L. Burden and J. Douglas Faires,
*Numerical Analysis*(3rd edition; PWS-Kent, 1985), Ch. 5, "Initial-Value Problems for Ordinary Differential Equations," pp. 199-289. - Forman S. Acton,
*Numerical Methods That Work*(Harper & Row, 1970), Ch. 5, "Ordinary Differential Equations — Initial Conditions," pp. 129-56. - William H. Press, Brian P. Plannery, Saul A. Teukolsky,
and William T. Vetterling,
*Numerical Recipes: The Art of Scientific Computing*(Cambridge University Press, 1986), Ch. 15, "Integration of Ordinary Differential Equations," pp. 547-77. - K. F. Riley,
*Mathematical Methods for the Physical Sciences*(Cambridge University Press, 1974), pp. 319 ff. - Goldstein, pp. 143 ff, 606 ff.
- Wilson, Decius, and Cross, pp. 285 f.
- Richard P. Feynman,
*The Feynman Lectures on Physics, Vol. I: Mainly Mechanics, Radiation, and Heat*(Addison-Wesley, 1964), Ch. 4, "Conservation of Energy", pp. 4-1 ff. - Feynman, Vol I, pp. 12-3 ff.
- For a discussion of rotation of axes in three-space
within computer graphics, see Cornel Pokorny,
*Comptuer Graphics: An Object-Oriented Approach to the Art and Science*(Franklin, Beedle & Associates, 1994), pp. 152 ff. - Floyd H. Johnson, personal communication of 27 February 1996