# Simulating Physical Systems

Posted: Friday, February 29, 2008

Due dates:

• Part 1, Ideal gas simulation, Friday, March 7, 10pm;
• Part 2, Ising spin simulation, Friday March 21, 10pm

## Introduction

You have already seen a number of examples using randomization to solve problems and to generate experimental results. In the percolation project, grids were generated according to a specified probability distribution to experimentally determine the percolation probability. This approach is an example of a Monte Carlo method, an algorithm that relies on repeated random sampling to compute results. Monte Carlo methods are often used for simulating physical and mathematical systems when other methods are impossible or computationally infeasible.

## Objective and Background

The objective of this project is to use the Monte Carlo “demon algorithm” to estimate parameters in physical systems. We will examine two systems, an ideal gas with moving molecules and a magnetic grid, also called lattice, of magnetized particles having spins. We could simulate such systems by computing the pairwise interactions among all molecules/particles at every time step in the simulation, but that approach would be computationally intractable. Moreover, in these systems we're not typically interested in the specific behavior of each individual molecule/particle, but rather in the overall system behavior. For example, we want to know the expected distribution of molecule velocities in a gas (e.g., to understand its temperature), not which molecules have collided with which other molecules.

The key concept behind the demon algorithm is that any simulation must both preserve the total energy of the system (energy is conserved) and insure that no particle winds up with negative kinetic energy. The algorithm proceeds by choosing a molecule/particle at random and proposing to make a small random change to its energy (e.g., to simulate the effect of a collision). To compensate for the change in energy of the molecule/particle, the opposite change is made to a special “demon” molecule/particle, thus conserving the total energy of the system.

In the ideal gas simulation we change molecule energy by adjusting its velocity and in the magnetic grid we change particle spins. Computationally, the project uses many of the concepts you have used before: loops, lists, conditionals, functions, and VPython and MatPlotLib for visualization.

Before describing the demon algorithm in detail, we give a brief description of the two systems to be simulated. In the ideal gas, we want to study the velocity of molecules in the gas. This average velocity gives a measure of the temperature of the gas. Consider the state of the gas in an insulated container at an instant of time. Each of its molecules has some velocity that remains constant until the next collision involving the molecules. Because the molecules move quickly, on average, many of them will have experienced collisions during a short time span. Consequently, most of the molecules will have substantially different velocities than they originally did. This process can, in principle, be simulated directly. In practice, doing so is computationally too expensive or is infeasible (e.g., one step of the situation could mean solving n differential equations, with each equation having n² terms).

The second system considered is the Ising Model, which simulates an idealized 2-D lattice arrangement of particles that have a spin; i.e., their state is either +1 or -1. Such a system models, for example, magnetism and allows one to estimate average magnetization at varying energy levels. The Ising model has been studied extensively and exact solutions are computationally not feasible. For an n x n grid, exact solutions would need to consider 2^(n²) states; for a 16×16 grid, this corresponds to ~10^77 states contributing to averages.

## Sketch of Demon Algorithm

The demon algorithm is a common technique for simulating selected parameters of a complex system by maintaining a given macroscopic energy, while randomly changing the microstate of the system. If a trial change would result in a reduction of energy in the system, additional energy is given to a demon (fictitious) particle. If the change would result in additional system energy, the difference is taken from the demon (assuming it has sufficient energy). In our two simulations, making a trial change will refer to changing a particle's velocity or magnetic spin, respectively.

1. Initialize the system state. The sum of the system energy and the demon energy should be equal to a specified total energy. The demon energy must always be non-negative.
2. Execute one step of the simulation. One simulation step consists of making a specified number of trials, with each trial doing the following: Make a trial change to the state of the system and record the change in system energy as `deltaEnergy`. If the change is legitimate, i.e., if the demon energy resulting from making the change is non-negative, the change is accepted and we set
• `demonEnergy = demonEnergy - deltaEnergy`
• `systemEnergy = systemEnergy + deltaEnergy`
• Update any parameters influenced by making one trial.
3. Repeat making trials until one simulation step is complete.
4. After one step, update any parameters influenced by one step. Execute as many steps as specified by the simulation.

It is not obvious that this process can effectively simulate a physical system. In fact, there is no rigorous proof that it does. In practice, however, this algorithm works very well. The first part of the project uses the demon algorithm in an ideal gas simulation.

## Part I: Ideal gas simulation

An ideal gas is a simple, uniform system: Molecules with energy are bouncing around inside a closed container. The internal molecular dynamics is simple – free motion with occasional intermolecular collisions. When using the demon algorithm for a simulation, we view the result of a collision as a random change of the molecule's velocity. Each molecule has a particular momentum, and collisions exchange part of that momentum.

Consider the velocity of a single molecule over time. Each collision can be thought of as a random change in velocity. Since we are interested in the average velocity of all molecules, the exact number and time of the collisions is not important. Rather, what matters is the total change in velocity over the time period we are considering. The demon algorithm takes advantage of this idea. We choose a total energy for the system (view it as the temperature). One step of the simulation consists of a number of trials. One trial perturbs the velocity of a random molecule by a random amount. The change is accepted if system and demon energy are less than the total energy. The difference between the system energy and the total energy is held by the algorithm’s namesake, the demon.

The core of the simulation is a function ideal_gas:

#### ideal_gas(N, totalEnergy, steps, hist=True)

The function simulates N molecules in a container with system energy bounded by totalEnergy. It executes steps steps, with each step consisting of N demon trials. The fourth parameter controls whether or not histograms are generated, which is described in more detail in the experimental work section.

The energy of the system state is fully described by the molecular velocities. Recall that the energy of a moving molecule is ½mv², where m is the mass of the molecule and v is the velocity. For our simulation, assume unit mass (i.e., m=1). Steps 1 and 2 of the demon algorithm now look like this:

1. Initialize the system state by setting an initial velocity for each molecule.
• Set the velocity of each molecule such that each accounts for an equal portion of the totalEnergy. The velocity of each molecule in the initial state is thus `sqrt(2*totalEnergy/N)`.
• This setting implies that `systemEnergy` = `totalEnergy` and leaves no energy for the demon (`demonEnergy=0`) since the `systemEnergy` + `demonEnergy` must not exceed `totalEnergy`.
2. Execute `steps` iterations, where one step consists of N trials. Making a trial means choosing a random molecule and changing its velocity by a random value in (`-Dv`, `Dv`), hence changing its energy. In all simulations, choose Dv to correspond to 1/10th of the totalEnergy. If the old velocity was `vOld` and the new velocity is `vNew`, then the change in system energy is `deltaEnergy = vNew²/2 - vOld²/2`. Use the description of step 2 given in the demon sketch to determine whether to accept or reject the change.

Note that since random choices determine which molecules are selected and not all trial changes are allowed, it is not the case that every molecule changes in a step. The entire simulation will execute up to N×steps velocity changes, resulting in an equal number of demon energy changes. Function `ideal_gas` returns one value, the average system energy, where the average is taken over the system energy values after each simulation step.

We will sample the demon energy after every simulation step and store these demon energy values in a list of size step. If histograms are to be generated, the step demon values will be used in histogram Demon_Energy. The final velocities of the N molecules will be reported in histogram Final_Molecule_Velocity. This is described in more deatil in the next section.

### Experiments for the Ideal Gas Simulation

Your experimental work should consider systems with `N`= 50 to 500, with an increment of 50. For each N, run function ideal_gas with steps = 3000 and totalEnergy = 500. You should use small values of N and steps for testing your code. The simulation should generate one graph and two types of histograms. You can use MatPlotLib or VPython to plot the graph and the histograms (all plots can be generated after the simulation).

The graph N_versus_Energy uses the value returned by function ideal_gas, which represents the average system energy of a simulation with particular value of N. You should graph the average system energy on the y-axis and N on the x-axis. In the same plot, show the total system energy in a different color (in the required experiment it is 500 for every value of N). You can generate the plot after the entire simulation or generate it incrementally after ideal_gas returns a value.

The ideal gas simulation should generate two types of histograms. Whether and when the histograms are generated is determined by the fourth parameter of ideal_gas. For example, `ideal_gas(N, 500, 3000, N==500)` will plot the histograms only when N=500. Note that omitting the 4th argument will plot the histograms every time ideal_gas is called. The histograms are generated within function ideal_gas.

Histogram Final_Molecule_Velocity uses the values of the N molecule velocities after steps iterations of ideal_gas (i.e., the values just before ideal_gas returns). Plot the velocities on x-axis and their frequencies on y-axis. You should use the signed velocities and 0 is thus in the center of the x-axis. You should experiment with the most appropriate number of bins.

Histogram Demon_Energy plots the observed demon energy. The plot is generated after the steps steps of ideal_gas are completed and it records the demon energy after each step. This generates an array with steps demon energies. Plot the demon energies on the x-axis and their frequencies on the y-axis. Again, select an appropriate number of bins.

If you are using VPython, you can show a histogram plot as soon you have data needed. If you are using Matplotlib, examples shown in class generate the plots before the simulation terminates. Here is an example on how to show Matplotlib plots during the computation.

You need to hand in a discussion of your experimental results. The following questions should be addressed:

• In graph N_versus_Energy, how close is the actual system energy to the total energy? How does it change as N changes?
• Using histogram Final_Particle_Velocity, what kind of distribution do you judge the velocities to come from? How close is the molecule velocity to the initial velocity of `sqrt(2*totalEnergy/N)`?
• Using histogram Demon_Energy, what kind of distribution do you judge the demon energies to come from?

## Part II: The Ising Model

The Ising model is a physical system defined on a grid of spins (in the literature, this grid is often called a lattice). Consider a 2-dimensional grid. Every element in the grid represents a magnetized particle that either points up (value +1) or down (value -1). The figure below shows two views of a 3×4 grid, with spins represented by arrows as well as colors (red is +1, white is -1). Particles interact with their neighbors: two adjacent particles pointing in the same direction contribute one unit of negative energy, while particles pointing in opposite directions contribute one unit of positive energy. We assume periodic boundary conditions, i.e., the left-most particle is adjacent to the right-most particle in the same row, and the top-most particle is adjacent to the bottom-most particle in the same column. This topology is that of a torus (but we draw it as a grid).

As an example, the left 3×4 grid has a total energy of 0. The contribution of each interaction for the circled spin is labeled. The result of the interactions when flipping that spin is shown in the right grid–the initial energy contribution of the circled particle is +2 and it changes to -2. The total energy of the lattice changes from 0 to -4.

The core of the simulation is a function ising_model:

#### ising_model(N, totalEnergy, steps, vis=True)

The simulation determines the average magnetization of the grid for a particular system energy. The basics of the demon algorithm remain the same.

• Choose as initial state one with the lowest energy (i.e., all spins are in the same direction). This state has a system energy of -2N². Your experimental results will show that the system will quickly reach an equilibrium independent of the initial state.
• The simulation runs for steps iterations, with each step consisting of N² trials. One trial corresponds to perturbing the energy of the system by flipping the magnetization of a random particle. Calculating the change in energy after a flip needs to consider the interactions with particle's four neighbors (the torus gives every particle four neighbors). For one flip, the energy could increase or decrease by at most 4. The change in energy is to be “absorbed” by the demon particle. A proposed change is rejected if it would make the demon's energy negative.

In summary, function ising_model runs the simulation on a grid of size `N×N`, performing steps iterations, with each iteration performing `N²` attempted changes (through flips). The fourth parameter of ising_model controls whether visuals are generated, similar to the fourth parameter in ideal_gas.

Function ising_model returns three values: M_aver, the average magnetization, M2_aver, the average square magnetizations, and the associated variance. The magnetization of a lattice is the sum of the magnetization of its particles (the change of one particle thus changes the magnetization by 2). Function ising_model records the magnetization and the square of the magnetization after each trial (this means storing steps*N*N values). M_aver, M2_aver, and the variance are computed just before ising_model returns. You can use any the of numerically stable methods seen in Lab 8 to compute the variance.

### Experiments for the Ising Spin Simulation

Make sure your program runs correctly and you have tested it with smaller values!

For the experimental portion of the Ising spin simulation, use first N=20 and steps=400 (this implies 160,000 trials in ising_model). The simulation will consider total energies from -2N² to 0 (-800 to 0), with a step size of 25. This makes 32 call to ising_model.

Run the simulation once with N=40 and steps=300. The simulation will consider total energies from -2N² to 0, which is from -3200 to 0. Choose the step size so that again 32 calls to ising_model are made. In general, the step size should be 2N*N/32. This allows you to run larger grid sizes in a reasonable time (but you will experience a slowdown).

The following visuals should be generated during the simulation:

1. Visualize the NxN grid during the entire simulation after every particle flip. Use two contrasting colors to show the grid after every change using VPython (similar to the VPython visualization done in lab 4 and the percolation project).
2. Histogram Demon_Energy plots the observed demon energy. The histogram is only generated when visualization is turned on and it generates one plot for each call to ising_model.
• The plot is generated after the steps steps of function ising_model are completed and it records the demon energy after each step. These calculations generate an array with steps demon energies.
• Plot the demon energies on the x-axis and their frequencies on the y-axis. You should be able to reuse the histogram from the ideal_gas simulation.
3. Graph Magnetization_versus_Energy_and_SimulationSteps shows the magnetization of the grid (on the y-axis) as a function of both total energy (z-axis) and simulation steps (x-axis). It is only shown when visualization is turned on.
• One data point is generated after one step of the simulation. Note that the steps values generated in one call to ising_model represent one “curve” in the plot.
• A single data point is a triple: index of simulation step, total energy (which varies from -2*N² to 0), and current magnetization of the grid.
• This plot is drawn in 3D. See the provided VPython code skeleton to get started.
4. Graph Magnetization_versus_Energy illustrates statistical properties of magnetization. This graph is generated using the values returned by ising_model (recall that in the simulation specified, ising_model is called 32 times).
• The graph shows the total system energy in the x-axis and the average magnetization M_aver and the average square magnetization M2_aver on the y-axis.
• To achieve a more uniform plot, the values of M_aver and M2_aver should be scaled to be less than 1 (i.e., use M_aver/N*N and M2_aver/N*N).
• The plot should also show the standard deviation of magnetization as functions of total energy (i.e., sqrt(variance)).

Before you describe your experiment, make sure you have run a number of simulations and have seen at least ten Magnetization_versus_Energy graphs. You should address the following questions.

• Demon_Energy histogram
• Is the Demon_Energy histogram stable across simulations?
• How does it change with varying energies?
• Is the distribution the same as for the ideal gas model?
• Magnetization_versus_Energy graph
• Which curves are most stable?
• What aspects of the grid behavior are revealed in this graph?
• What does the graph mean, in terms of the system being modeled (magnetization)?
• Magnetization_versus_Energy_and_SimulationSteps
• Is this graph more or less useful than the Magnetization_versus_Energy plot?
• Does this graph reveal more behaviors of the lattice?
• Why does the average magnetization approach zero at higher energies?

## Submission instructions

#### Part I

For part I, submit function `ideal_gas` and the main program driving the simulation in `ideal_gas.py`. Your main program will contain the loop for the values of N considered. The rest of the code will differ significantly between students depending how the visualization is handled (.e., the graphs N_versus_Energy, Final_Particle_Velocity, and Demon_Energy). You also need to include a text file (.txt, .doc or .rtf) that includes the discussion of your experiment and answers the questions listed. These files should all be placed into one folder and submitted. You will have the option to revise your discussion of the experiment. However, in order to be able to do so you need to submit an initial discussion and a submit your revision with part 2.

#### Part II

For part II, submit function `ising_model` and the main program driving the simulation in `ising_model.py`. Your main program will contain the loop for the different total energies considered. The rest of the code will differ significantly between students depending how the visualization is handled. You also need to include a file (.txt, .doc or .rtf) that includes the discussion of your experiment and answers the issues listed. These files should all be placed into one folder and submitted.