Posted: Monday, April 20, 2009

Due date: Friday, Mai 1, 2009, 10pm

In part I of the project on the Monte Carlo “demon algorithm” you used the demon algorithm to simulate a simple, uniform physical system representing ideal gas with moving molecules inside a closed container. The extension of this project considers 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.

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). See the March 11 lecture for additional detail and more background.

The energy of an Ising lattice is the sum of contributions representing the interactions between each spin and its nearest neighbors. 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 simulation uses the demon algorithm as described in part I. The core of the Ising simulation is the function ising_model(N, totalEnergy, steps, vis=True). Its arguments are as follows:

1. N: the side of the square grid
2. totalEnergy: The total amount of energy. That is, the energy of the N×N lattice + the demon energy.
3. steps: The number of steps to execute. Each step consists of N^2 demon trials, as described below.
4. visuals: Whether or not to draw visuals for the simulation. Similar to the corresponding parameter in ideal_gas.

Function ising_model runs the simulation on a grid of size N×N, performing steps iterations, with each iteration performing N² trials (through flips). 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 8 (when changing from -4 to +4, or vice versa). 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.

For the Ising simulations, choose as the initial state the state with the lowest energy (i.e., all spins are in the same direction). This state has a system energy of -2N². The total energy will vary from -2N² to 0. Note that when the total energy is -2N², the demon energy is 0. For this particular state, the system cannot change from its initial state, because there is no energy for the transition. However, it is an easy state to start the simulation in. Once the system energy -2N² + C, the demon has an energy of C, and now the system can move out from its initial state by taking energy from the demon. Your experimental results will show that the system will quickly reach an equilibrium independent of the initial state.

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 returns three values:

• M_aver, the average magnetization
• M2_aver, the average square magnetizations
• the variance associated with M2_aver

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 to compute the variance.

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

For the experimental portion of the Ising spin simulation, use 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)).

Finally, describe the results of your experiment. Make sure to run a number of simulations (using other parameters) so you have seen the a number of 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?
  • Why does the average magnetization approach zero at higher energies?

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.