PHYSICS 141

Winter 2004

Lecture 1: Computer Modeling

The Three Components of Computer Simulations

Computer modeling plays a very important role in science today. In the past, physical sciences were characterized by an interplay between experiment and theory. In experiment, a system is subjected to measurements, and results, expressed in numeric form, are obtained. In theory, a model of the system is constructed, usually in the form of a set of mathematical equations. The model is then validated by its ability to describe the system behavior in a few selected cases, simple enough to allow a solution to be computed from the equations. In many cases, this implies a considerable amount of simplification in order to eliminate all the complexities invariably associated with real world problems, and make the problem solvable.

In the past, theoretical models could be easily tested only in a few simple ``special circumstances''. So, for instance, in condensed matter physics a model for intermolecular forces in a specific material could be verified in a diatomic molecule, or in a perfect, infinite crystal. Even then, approximations were often required to carry out the calculation. Unfortunately, many physical problems of extreme interest (both academic and practical) fall outside the realm of these ``special circumstances''. Among them, one could mention astrophysical problems, like star formation, collisions of galaxies, and the large scale structure of the universe, the theory of the strong force which confines quarks permanently within the nucleon, termonuclear fusion, turbulunce, the physics and chemistry of defects, surfaces, clusters of atoms, organic molecules, involving a large amount of degrees of freedom; an accurate treatment of temperature effects, including anharmonicities and phase transitions; disordered systems in general, where symmetry is of no help to simplify the treatment; and this list could be continued almost indefinitely. Each computer simulation project outside the ``special circumstances'' requires three major components:

  • The physical problem and its theoretical model

    In computer modeling you are first and foremost a theorist. It is necessary to master the underlying physics of the proble we attack on the computer and to understand analytically as much as possible. We only want to ask the computer to do the things which cannot be done otherwise. The more analytic input we use, the more robust and powerful the computer simulation becomes. Needless to say, the development of computer experiments altered substantially the traditional relationship between theory and experiment. On one side, computer simulations increased the demand for accuracy of the models. For instance, a molecular dynamics simulation allows to evaluate the melting temperature of a material, modelled by means of a certain interaction law. This is a difficult test for the theoretical model to pass--and a test which has not been available in the past. Therefore, simulation ``brings to life'' the models, disclosing critical areas and providing suggestions to improve them. On the other side, simulation can often come very close to experimental conditions, to the extent that computer results can sometimes be compared directly with experimental results. When this happens, simulation becomes an extremely powerful tool not only to understand and interpret the experiments at the microscopic level, but also to study regions which are not accessible experimentally, or which would imply very expensive experiments, such as under extremely high pressure. Last but not least, computer simulations allow thought experiments--things which are just impossible to do in reality, but whose outcome greatly increases our understanding of phenomena--to be realized. Fantasy and creativity are important qualities for the computer modeller!

  • Algorithm and software implementation

    The advent of high speed computers--which started to be used in the 50s--altered the picture by inserting a new element right in between experiment and theory: the computer experiment. In a computer experiment, a model is still provided by theorists, but the calculations are carried out by the machine by following a ``recipe'' (the algorithm, implemented in a suitable programming language). In this way, complexity can be introduced and more realistic systems can be investigated, opening a road towards a better understanding of real experiments.

  • Analysis and visualization

    Analysis involves the physical interpretation of the data generated by the computer simulation. Visualization tools are indispensable in the interpretation of the results. For illustration we will discuss briefly the three components on the main theme of this course which is the computer modeling of N-body problems with particular applications to the collision of two galaxies and molecular dynamics.


    • Example: N-Body Modeling

    This course is an introduction to N-body simulations. We will briefly survey here the three main components of simulation projects on our chosen topic of the course.

    1. Physics and the theoretical model

    This part will illustrate the physical problem and its theoretical model as the first component of computer modeling in the general topic of N-body simulations. We will start with the underlying theoretical framework which is general enough to capture a very broad range of physical applications.

  • N-body equations of motion

    Any gravitational system may be idealized as a collection of N point-sized bodies, each with mass m_i, position r_i, and velocity v_i. The hamiltonian for such a system is 

    
                       ----              ----
                     1 \           2   G \      m_i m_j
(1)         H(r,v) = -  |   m_i v_i  - -  |   ----------- ,
              ~ ~    2 /        ~      2 /    |r_i - r_j|
                       ----              ----  ~     ~
                        i                i,j#i
			
			
	    G_Newton = 6.673(10) x 10^-11 m^3/kg/s^2 ,
    where H depends on all body positions and velocities, the first sum runs over all N bodies, the second runs over all pairs of bodies (twice, hence the factor of 1/2), and G is the gravitational constant. Then the equations of motion are 
                                             ----
            d(r_i)             d(v_i)    \    m_j (r_j - r_i)
	      ~                  ~                 ~     ~
(2)         ------ = v_i ,     ------ = G |   --------------- ,
                     ~
              dt                 dt      /     |r_j - r_i|^3
                                         ----   ~     ~
                                          j#i
    where the sum runs over all bodies except body i.

    N-body systems obey several basic conservation laws which are derived by manipulating Eq. 2. However, they may also be recognized directly from the form of the hamiltonian; Noether's theorem states that each symmetry of H gives rise to a corresponding conservation law. A symmetry of the hamiltonian is a transformation which leaves the physical system unchanged. For example, translation in time, t -> t + Delta t, is a symmetry of Eq. 1 because H is not an explicit function of time; consequently the total system energy E = T + U = H is conserved. Likewise, symmetry with respect to translation in space, r -> r + Delta r, implies conservation of total linear momentum, and symmetry with respect to rotation gives rise to conservation of total angular momentum.

    Eq. 2 is easy to generalize to molecular dynamics, if we replace the gravitational force by the appropriate molecular force acting between molecules in the system we are modeling.

  • Galaxy collisions

    We will illustrate the power of solving the N-body equations with two examples. The first is the simulation of two colliding galaxies on parabolic orbits as shown in the figure (sequence) . At close encounter during the collision, the simulation is a close match to what is observed on the sky for galaxy N4676 which is also known as the Mice. The observation on the sky is shown on the left and the simulation results are shown on the right:

    Josh Barnes and John Hibbard

  • Fullerene-based Nanoscale Gears

    D.H. Robertson, B.I. Dunlap, & C.T. White

    Composition of nanogears

    These nanogears are constructed totally out of Carbon. The different colors of the structures above is used only to clarify the molecular dynamics simulations to be presented later. The nanogears as depicted were constructed in two parts: the gear and the shaft. They are formed using the standard fullerene rules for forming closed surface structures. i.e.
    number_pentagons - number_heptagons - 2*number_octagons = 12
    For these gears, there are 4 pentagons at the tip of a tooth and a octagon between adjacent teeth. This equates to 24 pentagons and 6 octagons which if we place in the above equations gives 24 - 2*6 = 12. Satisfying the rull for generating a closed fullerene.

    Next these gears were joined to an appropriate shaft. In this case, part of a graphitic tubule. In terms of joining the tubule to the gear, a tubule with a similar symmetry as the gear (6-fold rotation axis) was chosen. In tubule terms this is a (6,6) tubule or serpentine tube. The tubule is capped on one end by part of C60 (containing 6 pentagons). The end joining to the gear uses 6 heptagons (allowing the negative curvature). Note that in the shaft the number of heptagons equally balance the number of pentagons -- allowing the structure to remain a closed fullerene.

    Energetics from Reactive Hydrocarbon Potential

    After forming these gears, they were minimized using the Brenner Reactive Hydrocarbon Potential (Phys. Rev. B. 42, 9458 (1990)). In terms of stability, the gears (starenes) as well as the total gear/shaft assembly is overall more stable than Buckminsterfullerene (C60) although it may have some regions with a higher strain energy (tips of the teeth).

    • 2. Algorithm and software implementation

    The 6N coupled ordinary differential equations as displayed in Eq. 2 have to be solved numerically using some finite discretization of the time variable. Variable discretization schemes are used in the literature, like the leapfrog algorithm, or the Nordsieck predictor-corrector method. The force calculation on the right side of Eq. 2 is also a challenging problem. A naive approach will require N^2 calculations, sophisticated procedures, like FFT or the Treecode will reduce the number of operations to NlogN.

    To keep software development accessible to every student, the code we will work with will always be available in Fortran and C programming languages. It will be the student's choice which language to use. Since the discussions and applications will never be completely symmetric, it is advantageous if teams have some knowledge in both languages.

    The Fortran implementation of the famous Treecode for the collision of two galaxies is about 3000 lines and organized in a very transparent fashion. We will also work in the class with its equivalent C language implementation.

    3. Analysis and visualization

    The above two examples of N-body simulations are briefly discussed here.

  • Two movies on the Mice (N4676)

    The first movie Mice is the final product of the project on galaxy collision. The second movie Mice (rotate) shows that the two tidal tails are not in the same plane which explains the observed skyview of N4676.

  • Molecular Dynamics Simulations of One Gear Driving the Other

    After the gears were minimized, molecular dynamics simulations were performed to study the dynamic behavior of these gears and whether or not these gears would make robust components of future nanomachines. The molecular dynamics (MD) simulations were performed using the Brenner hydrocarbon potential and a Nordsieck predictor-corrector algorithm. These simulations were performed at constant energy and the reactive potential allows these gears to break if appropriate.

    In these MD simulations the back half atoms of the shafts of each of the gears were constrained to only rotate. Thus the shafts were held fixed in space. To perform the simulation the constrained atoms of the red gear were given an angular velocity that ramped up over time until reaching a final angular velocity and then these atoms continued to rotate at this velocity. The white gear only acted in response to the red gear with the constraint that the back half of the shaft could only rotate. These conditions allowed the gears to be driven and held fixed in space by their shafts while still allowing part of their shafts and their gear-heads to flex and respond under the forces of the reactive potential. The simulations below were run for 100ps with two terminal angular velocities, 0.1 rev/ps and 0.5 rev/ps, and two ramp times, 50ps and 10ps.

    Simulations Results

    In the graphs below the green line is the angular velocity of the driven atoms of the red gear, the cyan line is the angular velocity of the white shaft, and the red line is the internal temperature of the two gears.
    0 to 0.1 revolution/ps in 50ps
    From the graph it can be seen that the green line rises to 0.1 rev/ps by 50ps and then levels off. The cyan line for the most part follows the green line although there are variations in its angular velocity. This is a beating or chattering as the white gear is being forced to rotate. Besides this chattering the white gear matches angular velocity and is driven by the red gear. The temperature of the system continues to rise throughout this process. This increase in temperature is results partly because energy is continuously being added to this system with no mechanism for heat dissipation (or dampening of the gear chatter). A MPEG movie of this simulation is available (click for movie). Alternately a series of snapshots is available taken every 5ps (right to left, top to bottom) is available. From these animations or images it can be seen that the gear-heads do fit together but also flex and possible slip.