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.


Table of Contents

Click here for the paper in Microsoft Word format.


Abstract

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.

Problem Description

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 friction­like 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.

Program Design

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 n second-order differentials (the Lagrangian formulation) [5] or as a system of 6 n first-order differentials (the Hamiltonian formulation). [6]

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]

F1 = (K /[ 12]2e12

— where F1 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 e12 is the unit vector from particle 2 to particle 1. The force acting on particle 2 is just the negative: F2 = -F1 .

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 F2 = -F1: 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 fj+1 = fj + h (df/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 fj+1 = fj-1 + 2h (df/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: fj­1 fj - h (df/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.

Results

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

Summary

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.

Availability of Source Code

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

Acknowledgments

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.


Notes

  1. R. Stephen Berry, Stuart A. Rice, and John Ross, Physical Chemistry (John Wiley & Sons, 1980), p. 381.
  2. 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.
  3. 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.
  4. An example of such a simplified coordinate system is found in "atomic units," as discussed in Berry, Rice, and Ross, p. 180.
  5. Herbert Goldstein, Classical Mechanics (2nd edition; Addison-Wesley, 1980), Chapter 1-2, pp. 1-69.
  6. Goldstein, Chapter 8, pp. 339-77.
  7. Richard P. Feynman, The Feynman Lectures on Physics, Vol. II: Mainly Electromagnetism and Matter (Addison-Wesley, 1964), pp. 4-2 ff.
  8. 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.
  9. Forman S. Acton, Numerical Methods That Work (Harper & Row, 1970), Ch. 5, "Ordinary Differential Equations — Initial Conditions," pp. 129-56.
  10. 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.
  11. K. F. Riley, Mathematical Methods for the Physical Sciences (Cambridge University Press, 1974), pp. 319 ff.
  12. Goldstein, pp. 143 ff, 606 ff.
  13. Wilson, Decius, and Cross, pp. 285 f.
  14. 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.
  15. Feynman, Vol I, pp. 12-3 ff.
  16. 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.
  17. Floyd H. Johnson, personal communication of 27 February 1996