Tutorials:Code-02 C++

Jump to: navigation, search

This tutorial shows how to write a Monte Carlo simulation in C++ using the ALPS alea library for the evaluation of observables. In a second step we will write the measurements into a standard HDF5 file format which allows us to use tools from the ALPS suite for further data analysis and plotting.

As a simple example, we will write a simulation of the classical 2D Ising model with local updates. The file ising-skeleton.cpp contains a skeleton code which already has all the infrastructure we will need: First it includes all needed headers, then it initializes a random number generator and three alps::RealObservable objects. Then it sets up a square lattice of Ising spins. It also provides a table of probabilities that can be used for Metropolis updates. The interface is the same as in the python script you implemented in the previous tutorial .

Your job is again to complete the methods step() and measure(): step() should choose a random spin from the lattice and flip it with the Metropolis probability p_{accept} = min(1,e^{-\beta \Delta E}) where \Delta E is the energy change the spin flip would cause. measure() determines the energy and magnetization of a spin configuration and adds this sample to the observable objects.

After replacing all ellipses with code, you can compile the simulation with this Makefile: Save the Makefile to the same directory as the .cpp file, edit the second line to point to your ALPS installation (if you haven't already set the environment variable ALPS_ROOT) and type make. This will produce an executable ising. Run it and you will see a scan over different values of \beta = 1/k_B T.

You can reuse your Binder cumulant python script from the previous tutorial in exactly the same way:

data = pyalps.loadMeasurements(pyalps.getResultFiles(pattern='ising.L*'),['m^2', 'm^4'])

for i in range(len(m2)):
    d = pyalps.DataSet()
    d.props = m2[i].props
    d.x= m2[i].x
    d.y = m4[i].y/m2[i].y/m2[i].y

and plot the Binder cumulant using the commands:

plt.xlabel('Inverse Temperature $\beta$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')