diff --git a/exercises/ex04_solution/exchange_energies_L6.png b/exercises/ex04_solution/exchange_energies_L6.png
new file mode 100644
index 0000000000000000000000000000000000000000..1a3534f0b3fb4c265d9b0c4ea15cbede6c70c594
Binary files /dev/null and b/exercises/ex04_solution/exchange_energies_L6.png differ
diff --git a/exercises/ex04_solution/exchange_energies_L7.png b/exercises/ex04_solution/exchange_energies_L7.png
new file mode 100644
index 0000000000000000000000000000000000000000..18e042aa76f8426ee9e3f4d14ede4abbdb14fa2b
Binary files /dev/null and b/exercises/ex04_solution/exchange_energies_L7.png differ
diff --git a/exercises/ex04_solution/heised.ipynb b/exercises/ex04_solution/heised.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..2e201362b215fd7f1d17e991b60e3e66634a348e
--- /dev/null
+++ b/exercises/ex04_solution/heised.ipynb
@@ -0,0 +1,120 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Exact Diagonalization\n",
+ "\n",
+ "A C++ code for exact diagonalization is attached to this solution sheet, and can be compiled by following the instructions provided with the skeleton code. Numeric results are available in the openchain.dat and periodicchain.dat for an open and periodic boundary conditions, respectively. The figures below can be reproduced with the `plot.py` script."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "<center>Figure 1: Ground state energy converging to value of −0.44325 for the thermodynamic limit.</center>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<img src=\"exchange_energies_L6.png\" width=\"300\"/>"
+ ],
+ "text/plain": [
+ "<IPython.core.display.Image object>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "text/html": [
+ "<img src=\"exchange_energies_L7.png\" width=\"300\"/>"
+ ],
+ "text/plain": [
+ "<IPython.core.display.Image object>"
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "from IPython.display import Image\n",
+ "from IPython.display import display\n",
+ "display(Image(url=\"exchange_energies_L6.png\", width=300), Image(url=\"exchange_energies_L7.png\", width=300))"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Figure 2: Strength of the local contribution of each bond to the Heisenberg Hamiltonian\n",
+ "for even and odd open chains."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "In Figure 1 we observe the oscillation of the ground state energy between even and odd number of sites in both lattice topologies. In the case of periodic boundary conditions this is easily understood by the appearance of a frustrated spin in an odd chain, i.e. one spin can align in order to favor the bond with the previous spin, or with the next spin, but not both at the same time.\n",
+ "\n",
+ "For an open chain there is no apparent frustration, but the increase in energy per site can still be understood from the local contributions to the Heisenberg Hamiltonian as shown in Figure 2. The major contribution to the ground state comes from the dimerized state, i.e. a chain of singlets: this is only possible in an even chain (see left panel). In the odd chain (right panel) the system tries to build singlets from both boundaries, which is again giving rise to frustration in the center where the two incompatible dimer patterns meet. In the thermodynamic limit L → ∞ the difference between even and odd length vanishes."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "\n",
+ "<center>Figure 3: Extrapolation of the energy gap. Extrapolation is performed only on the last\n",
+ "four lattice sizes.</center>\n",
+ "\n",
+ "By targeting the two lowest eigenenergies in the Lanczos solver (or a similar iterative\n",
+ "eigensolver) for each quantum number sector we identify the smallest excitation from the\n",
+ "ground state. Figure 3 shows how the energy gap nicely extrapolates to zero ∼ 1/L in\n",
+ "the thermodynamic limit."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.4.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/exercises/ex04_solution/heised_cpp/.gitignore b/exercises/ex04_solution/heised_cpp/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..15611dcce53a3f349eb894689871ce655a932081
--- /dev/null
+++ b/exercises/ex04_solution/heised_cpp/.gitignore
@@ -0,0 +1 @@
+alps-2.2*
diff --git a/exercises/ex04_solution/heised_cpp/Makefile b/exercises/ex04_solution/heised_cpp/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..f4e2fe584b1edc27f56a1f0a3244dbac8b5fd575
--- /dev/null
+++ b/exercises/ex04_solution/heised_cpp/Makefile
@@ -0,0 +1,12 @@
+CXX = g++
+CPPFLAGS = -std=c++14 -I${ALPSDIR}/boost -I${ALPSDIR}/alps/src
+
+ALPSDIR = alps-2.2.b4-src-with-boost
+
+heised: heised.cpp
+ ${CXX} ${CPPFLAGS} -o $@ $^ -llapack -lblas
+
+prepare_alps:
+ mkdir -p ${ALPSDIR}/build
+ -cd ${ALPSDIR}/build; cmake -D Boost_ROOT_DIR:PATH=../boost ../alps/
+ cp ${ALPSDIR}/build/src/alps/config.h ${ALPSDIR}/alps/src/alps/
diff --git a/exercises/ex04_solution/heised_cpp/energy_gap.png b/exercises/ex04_solution/heised_cpp/energy_gap.png
new file mode 100644
index 0000000000000000000000000000000000000000..8928c5d034caf6a5f386efec40f81b94a888422b
Binary files /dev/null and b/exercises/ex04_solution/heised_cpp/energy_gap.png differ
diff --git a/exercises/ex04_solution/heised_cpp/groundstate_energy.png b/exercises/ex04_solution/heised_cpp/groundstate_energy.png
new file mode 100644
index 0000000000000000000000000000000000000000..978505a5335b620084599a07eb4b9092222f315a
Binary files /dev/null and b/exercises/ex04_solution/heised_cpp/groundstate_energy.png differ
diff --git a/exercises/ex04_solution/heised_cpp/heised.cpp b/exercises/ex04_solution/heised_cpp/heised.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..deaf68726ac0fdd03b4b1e60ed50aff75628d14c
--- /dev/null
+++ b/exercises/ex04_solution/heised_cpp/heised.cpp
@@ -0,0 +1,214 @@
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <cassert>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/random.hpp>
+#include <boost/tuple/tuple.hpp>
+
+// FROM ALPS: #include <alps/utility/bitops.hpp>
+namespace alps
+{
+ inline static int BX_(long x)
+ { return ((x) - (((x)>>1)&0x77777777) - (((x)>>2)&0x33333333) - (((x)>>3)&0x11111111)); }
+
+ inline long popcnt(uint32_t x)
+ { return (((BX_(x)+(BX_(x)>>4)) & 0x0F0F0F0F) % 255); }
+}
+
+class SpinHalfBasis
+{
+public:
+ enum { NOT_FOUND=-1 };
+ typedef unsigned State;
+ typedef unsigned Index;
+
+ SpinHalfBasis(unsigned l, unsigned nups);
+
+ State state( Index i ) const { return states_[i]; }
+ Index index( State s ) const { return index_[s]; }
+ Index dimension() const { return states_.size(); }
+
+private:
+ std::vector<State> states_;
+ std::vector<Index> index_;
+};
+
+SpinHalfBasis::SpinHalfBasis(unsigned l, unsigned nups) // l: chain length; nups: spins-ups (fixes the symmetry sector).
+: index_(State(1<<l), NOT_FOUND) // initialize index_ to have 2^l entries each of which is the NOT_FOUND flag.
+{
+ // find all states with [nups] up spins
+ for( State s = 0; s < index_.size(); ++s )
+ {
+ if( alps::popcnt(s) == nups ) // if the state is in this sector keep it.
+ {
+ index_[s] = states_.size();
+ states_.push_back(s);
+ }
+ }
+}
+
+class HeisenbergHamiltonian : public SpinHalfBasis
+{
+public:
+ typedef double Scalar;
+ typedef boost::numeric::ublas::vector<Scalar> Vector;
+
+ HeisenbergHamiltonian(unsigned l, unsigned nups, bool periodic, double j);
+
+ void multiply(const Vector& x, Vector& y) const;
+
+private:
+ unsigned l_; // length of the chain
+ bool periodic_; // periodic boundary conditions off/on
+ double j_; // coupling J
+};
+
+HeisenbergHamiltonian::HeisenbergHamiltonian(unsigned l, unsigned nups, bool periodic, double j)
+: SpinHalfBasis(l,nups)
+, l_(l)
+, periodic_(periodic)
+, j_(j)
+{
+}
+
+/// Calculate y = H x
+void HeisenbergHamiltonian::multiply(const Vector& x, Vector& y) const
+{
+ // check dimensions
+ assert( x.size() == dimension() );
+ assert( y.size() == dimension() );
+
+ // diagonal part: +J/4 for parallel, -J/4 for antiparallel neighbors
+ State mask = (1<<(l_-1)) - 1; // contains a single 0 at the m.s.b. and 1s everywhere else
+ for( Index i = 0; i < dimension(); ++i )
+ {
+ State s = state(i);
+ int p = alps::popcnt(mask & ( s ^ (~s>>1) )); // number of parallel pairs. s ^ (~s>>1): implements an EQUAL over neighbour spins.
+ y[i] = .25*j_*( 2.*p - l_ + 1 )*x[i]; // # par - # antipar = 2 * # par - # bonds.
+ }
+
+ // off-diagonal part: {01,10} -> J/2 {10,01}
+ for( Index i = 0; i < dimension(); ++i )
+ {
+ State s = state(i);
+ for( int r = 0; r < l_-1; ++r )
+ {
+ State sflip = s ^ (3<<r); // flip spins (r,r+1). (3<<r): bit pattern with all zeros except at r and r+1.
+ Index j = index(sflip);
+ if( j != NOT_FOUND ) // we need to filter out events like e.g. (1,1) flip-> (0,0), since this kicks us out of the symmetry sector.
+ y[j] += .5*j_*x[i];
+ }
+ }
+
+ // periodic boundaries
+ if( periodic_ )
+ {
+ for( Index i = 0; i < dimension(); ++i )
+ {
+ State s = state(i);
+ // diagonal
+ int p = 1 & ( s ^ (~s>>(l_-1)) );
+ y[i] += .25*j_*( 2.*p - 1 )*x[i];
+ // off-diagonal
+ State sflip = s ^ ( 1 | (1<<(l_-1)) ); // ( 1 | (1<<(l_-1)) ): bit pattern with 1s at least and most sig. bits. the rest are zero.
+ Index j = index(sflip);
+ if( j != NOT_FOUND )
+ y[j] += .5*j_*x[i];
+ }
+ }
+}
+
+namespace ietl
+{
+ // overload mult function inside ietl for use with our HeisenbergHamiltonian objects.
+ inline void mult( const HeisenbergHamiltonian& h, const HeisenbergHamiltonian::Vector& x, HeisenbergHamiltonian::Vector& y )
+ {
+ h.multiply( x, y );
+ }
+}
+
+// ietl::mult() needs to be declared before including these
+#include <ietl/interface/ublas.h>
+#include <ietl/lanczos.h>
+
+std::pair< std::vector<double>, std::vector<int> >
+diagonalize( const HeisenbergHamiltonian& h, unsigned nvals=1, unsigned maxiter=1000)
+{
+ typedef ietl::vectorspace<HeisenbergHamiltonian::Vector> Vectorspace;
+ typedef ietl::lanczos<HeisenbergHamiltonian,Vectorspace> Lanczos;
+ typedef ietl::lanczos_iteration_nlowest<double> Iteration;
+ typedef boost::mt19937 Generator;
+
+ Vectorspace vspace(h.dimension());
+ Lanczos solver(h,vspace);
+ Iteration iter(maxiter,nvals);
+
+ solver.calculate_eigenvalues(iter,Generator()); // call the solver for the eigenvalues.
+
+ if( iter.iterations() == maxiter )
+ std::cerr << "Lanczos did not converge in " << iter.iterations() << " iterations." << std::endl;
+ else
+ std::cout << " Lanczos converged after " << iter.iterations() << " iterations." << std::endl;
+
+ return std::make_pair(solver.eigenvalues(),solver.multiplicities());
+}
+
+void do_chain(int l, bool periodic, double j, unsigned nstates, std::ostream& datfile)
+{
+ std::cout << "+++ Diagonalizing S=1/2 Heisenberg " << (periodic ? "periodic" : "open")
+ << " chain with L=" << l << ", J=" << j << std::endl;
+
+ std::vector<double> all_energies;
+
+ for( unsigned n = 0; n <= l; ++n ) // iterate over # up spins. (sectors)
+ {
+ double sz = n-l/2.;
+ std::cout << "--- Sector Sz=" << sz << ": ";
+ HeisenbergHamiltonian ham(l,n,periodic,j);
+ std::cout << ham.dimension() << " basis states" << std::endl;
+
+ unsigned nvals = std::min(nstates,ham.dimension());
+ std::vector<double> evals;
+ std::vector<int> mults;
+ boost::tie(evals,mults) = diagonalize(ham,nvals);
+
+ for( unsigned i = 0; i < nvals; ++i )
+ std::cout << n-l/2. << '\t' << evals[i] << '\t' << mults[i] << '\n';
+ all_energies.insert(all_energies.end(),evals.begin(),evals.begin()+nvals);
+ }
+
+ // print data to file
+ std::sort(all_energies.begin(),all_energies.end());
+ datfile << l;
+ for( unsigned i = 0; i < nstates; ++i ) datfile << '\t' << all_energies[i];
+ datfile << std::endl;
+}
+
+void write_dat_header(std::ostream& datfile,unsigned nstates)
+{
+ datfile << "# Sz";
+ for( unsigned i = 0; i < nstates; ++i )
+ datfile << "\tE" << i;
+ datfile << "\n";
+ datfile.precision(10);
+}
+
+int main()
+{
+ int lmin = 2;
+ int lmax = 20;
+ double j = 1.;
+ unsigned nstates = 2; // How many of the lowest eigenvalues (energies) do you want to calculate? ( nstates = 1 ... just the ground state)
+
+ std::ofstream opendata("openchain.dat");
+ write_dat_header(opendata,nstates);
+ for( int l = lmin; l <= lmax; ++l )
+ do_chain(l,false,j,nstates,opendata);
+
+ std::ofstream periodicdata("periodicchain.dat");
+ write_dat_header(periodicdata,nstates);
+ for( int l = lmin; l <= lmax; ++l )
+ do_chain(l,true,j,nstates,periodicdata);
+}
diff --git a/exercises/ex04_solution/heised_cpp/openchain.dat b/exercises/ex04_solution/heised_cpp/openchain.dat
new file mode 100644
index 0000000000000000000000000000000000000000..243b04d106dc6ae8fddc4df48c171c58441f8ee1
--- /dev/null
+++ b/exercises/ex04_solution/heised_cpp/openchain.dat
@@ -0,0 +1,20 @@
+# Sz E0 E1
+2 -0.75 0.25
+3 -1 -1
+4 -1.616025404 -0.9571067812
+5 -1.927886253 -1.927886253
+6 -2.493577134 -2.001995357
+7 -2.836239681 -2.836239681
+8 -3.374932599 -2.982240488
+9 -3.736321706 -3.736321706
+10 -4.258035207 -3.93067359
+11 -4.632093302 -4.632093302
+12 -5.142090633 -4.861147937
+13 -5.525322097 -5.525322097
+14 -6.026724662 -5.780492604
+15 -6.416920492 -6.416920492
+16 -6.911737146 -6.692460429
+17 -7.307408708 -7.307408708
+18 -7.797011069 -7.599284884
+19 -8.197105742 -8.197105742
+20 -8.682473334 -8.502378698
diff --git a/exercises/ex04_solution/heised_cpp/periodicchain.dat b/exercises/ex04_solution/heised_cpp/periodicchain.dat
new file mode 100644
index 0000000000000000000000000000000000000000..d59ce225e727a17ab7d7fdbd9d859da1f56c5b88
--- /dev/null
+++ b/exercises/ex04_solution/heised_cpp/periodicchain.dat
@@ -0,0 +1,20 @@
+# Sz E0 E1
+2 -1.5 0.5
+3 -0.75 -0.75
+4 -2 -1
+5 -1.868033989 -1.868033989
+6 -2.802775638 -2.118033989
+7 -2.855179257 -2.855179257
+8 -3.651093409 -3.128419064
+9 -3.797299784 -3.797299784
+10 -4.515446354 -4.092207347
+11 -4.718936363 -4.718936363
+12 -5.387390917 -5.031543404
+13 -5.62958433 -5.62958433
+14 -6.263549534 -5.956443824
+15 -6.533667572 -6.533667572
+16 -7.142296361 -6.872106678
+17 -7.433520165 -7.433520165
+18 -8.022749087 -7.781499637
+19 -8.330488454 -8.330488454
+20 -8.90438653 -8.686440986
diff --git a/exercises/ex04_solution/heised_cpp/plot.py b/exercises/ex04_solution/heised_cpp/plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..4da32ad8c0966611995696f4a2469713ed8786ae
--- /dev/null
+++ b/exercises/ex04_solution/heised_cpp/plot.py
@@ -0,0 +1,48 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+tasks = [('open','openchain.dat'),('periodic','periodicchain.dat')]
+
+# create two figures
+plt.figure(1)
+plt.title('Ground state energy per site')
+plt.xlabel('$1/L$')
+plt.ylabel('$E_0/L$')
+plt.figure(2)
+plt.title('Gap')
+plt.xlabel('$1/L$')
+plt.ylabel('$E_1-E_0$')
+plt.xlim(0,0.3)
+plt.ylim(0,0.8)
+
+for t in tasks:
+ # load data
+ name = t[0]
+ data = np.loadtxt(t[1])
+ lengths = data[:,0]
+ invlen = 1./lengths
+ e0 = data[:,1]
+ e1 = data[:,2]
+
+ # plot ground state energy per site
+ gs = e0/lengths
+ plt.figure(1)
+ plt.plot(invlen,gs,label=name)
+ plt.legend(loc='upper left')
+
+ # restrict to chains of even length
+ even = lengths%2 == 0
+ invlen = invlen[even]
+ gap = (e1-e0)[even]
+
+ # fit to form y = m*x = m/L for the 4 largest chains
+ numpoints = 4
+ m = np.linalg.lstsq(np.vstack(invlen[-numpoints:]),gap[-numpoints:])[0]
+
+ # plot gap and 1/L fit
+ plt.figure(2)
+ plt.plot(invlen,gap,'.',label=name)
+ invlen = np.concatenate((invlen,[0.]))
+ plt.plot(invlen,m*invlen)
+ plt.legend(loc='upper left')
+plt.show()
\ No newline at end of file