Correlation Functions for the Ising Model

Correlation functions are a major subject in statistical mechanics. Many experiments (i.e., scattering probes like X-rays) measure the correlation functions. Many scaling properties, both near critical points and in systems out of equilibrium, are best studied using correlation functions. Today we will study the spin-spin correlation functions C(r) = <S(i,j) S(i+r,j)> in the two-dimensional Ising model.

There are two parts to this assignment. We will provide a function CorrelationObserverSlow, which computes the correlation function in the most brute-force way possible. The first part of the assignment is to write a much faster version CorrelationObserverFFT, based on the one we provide, which uses fast Fourier transforms to compute the same function much faster.

Recognizing that some of the class needs some time to finish SWIGging their existing code, and that for some it may be hard to motivate writing a fast version of an observer whose usefulness isn’t yet clear, I’m encouraging some of the class to just use CorrelationObserverSlow, and skip the first part of the assignment. (You can always come back to it!)

 

Writing CorrelationFunctionFFT

  1. Copy CorrelationObserverSlow.h and CorrelationObserverSlow.cpp from the directory LMCSolution\LMCFiles to your LMCFiles directory.
  2. Change the default size of your simulation (set in the constructor to IsingSimulation) to 128x128. Our FFT package only works on vectors of length a power of two (although CorrelationObserverSlow works for any length vector).
  3. Construct a correlationObserver of type CorrelationFunctionSlow, and attach it to sim. Make sure it runs. See how much it slows down the simulation, especially for larger systems (512x512, for example).
  4. We have installed the BLAS routines and high-performance FFT routines onto the system. To link to them, set Project\Settings\Link, Object/library Modules, to
  5. F:\"Program Files"\Intel\PLSuite\lib\Pentium.II\Intel\mkl_s.lib.

  6. Copy the two CorrelationObserverSlow files into new files named CorrelationObserverFFT. Edit all the names "Slow" to "FFT" in the two files, using "Replace".
  7. Add the files to the project, change correlationObserver to type FFT, and make sure it runs.
  8. We’ll be using the C routine zfft1dc (z for double-precision, 1d, C language). At the top of the .h file, inside a #ifndef SWIG line, add the line extern "C" void zfft1dc(double*,double*,int,int,double*);
  9. The new Correlation(r) will automatically add r and width-r: you should remove width-r from the definition.
  10. You’ll need three protected vectors (double *): FTreal, FTimag, and wsave. Set them to NULL in the constructor. (We don’t know yet how big the width is.) Delete them if they exist in the destructor.
  11. In Update, when we test to see if the width is right, we should allocate the three new vectors. The FT vectors should have "width" components, and the wsave vector should have length equal to 3*width.
  12. At the same time (when width changes) call zfft1dc(FTreal,FTimag,width,0,wsave). The zero tells the routine to set up a table of coefficients: it doesn’t alter the FT arrays.
  13. The correlation function is given by the inverse Fourier transform of the absolute square of the Fourier transform of the data. Looping over rows i up to height,
  14. Fill the array FTreal[j] with S->Get(i,j);

    Zero Ftimag

    Call zfft1dc(FTreal, FTimag, width, -1, wsave); the minus gives the forward transform.

    Fill the array FTreal[j] with the absolute square, and FTimag with zeros.

    Call zfft1dc(FTreal, FTimag, width, 1, wsave); the 1 gives inverse

    Add FTreal[j] to correlationFunction[j], and increment numberAveraged.

  15. Compile, and run. Check to make sure that it gives the same answers, up to rounding errors, as did the slow version.
  16. Check to see if it slows down the simulation less than the old version. The FFT takes L log(L) steps, so the new correlation function takes time L2 logL, where the old correlation function scaled like L3.

The fast Fourier transform is so efficient we use it whenever we can.

SWIGging the Code

I’ve modified Makefile.win and Makefile.msc to pass the library path to the Intel routines. If you have skipped the first part, you can continue SWIGging with the old Makefiles. Otherwise, you should add

MYLIB = F:\PROGRA~1\Intel\PLSuite\lib\Pentium.II\Intel\mkl_s.lib

to your Makefile.msc, and copy my new version of Makefile.win.

Statistical Mechanics: Entropy, Order Parameters, and Complexity, now available at Oxford University Press (USA, Europe).