# Tutorials:Code-02 C++

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'])
m2=pyalps.collectXY(data,x='BETA',y='m^2',foreach=['L'])
m4=pyalps.collectXY(data,x='BETA',y='m^4',foreach=['L'])

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


and plot the Binder cumulant using the commands:

plt.figure()
pyalps.pyplot.plot(u)
plt.xlabel('Inverse Temperature $\beta$')
plt.ylabel('Binder Cumulant U4 $g$')
plt.title('2D Ising model')
plt.show()