ex11 solution

exercises/ex11_solution/Makefile

+include $(ALPS_ROOT)/share/alps/include.mk
+
+ising: ising.cpp
+	$(CXX) -Wall $(CPPFLAGS) $(CXXFLAGS) $(LDFLAGS) $(LIBS) -o ising ising.cpp

exercises/ex11_solution/data/run.py

+import numpy as np
+from math import *
+
+gamma = 1
+delta = 0.05
+J = np.linspace(0.5,2)
+nsweeps = {8: 1000000, 16: 1000000, 24: 1000000,
+          32: 1000000, 40: 1000000, 48: 1000000}
+
+i = 0
+n = 1
+for Ji in J:
+  if True:
+    basename = 'start_ML1_%s.sh'%n
+    fp = open(basename, 'w')
+    print 'bsub -W 2:00 sh '+basename
+    n += 1
+  i += 1
+  for ML in [1]:
+    for L in reversed([16,24,32,40,48]):
+     Jx = delta*Ji
+     Jy = -.5*log(delta*gamma)
+     fp.write('../ising %s %s %s %s %s\n' % (ML*L, L, Jy, Jx, nsweeps[L]))
+

exercises/ex11_solution/plot.py

+import sys, os
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pyalps
+from pyalps.plot import plot
+
+files = pyalps.getResultFiles(dirname='data')
+data = pyalps.loadMeasurements(files , ['|m|','m^2', 'Connected Susceptibility', 'Binder Cumulant U2'])
+
+for d in pyalps.flatten(data):
+    d.props['M/L'] = d.props['M'] / d.props['L']
+
+m = pyalps.collectXY(data, 'Jx', '|m|', foreach=['L', 'M'])
+chi = pyalps.collectXY(data, 'Jx', 'Connected Susceptibility', foreach=['L', 'M'])
+binder = pyalps.collectXY(data, 'Jx', 'Binder Cumulant U2', foreach=['L', 'M'])
+
+
+for d in pyalps.flatten(m):
+    d.x = np.exp(2.*d.props['Jy'])*d.x
+plt.figure()
+plot(m)
+plt.xlabel('$J/\\Gamma$')
+plt.ylabel('magnetization')
+plt.legend(loc='best', frameon=False)
+
+
+for d in pyalps.flatten(chi):
+    d.x = np.exp(2.*d.props['Jy'])*d.x
+plt.figure()
+plot(chi)
+plt.xlabel('$J/\\Gamma$')
+plt.ylabel('Connected Susceptibility, $\chi = [<m^2> -<|m|>^2] / L / M$')
+plt.legend(loc='best', frameon=False)
+
+
+for d in pyalps.flatten(binder):
+    d.x = np.exp(2.*d.props['Jy'])*d.x
+plt.figure()
+plot(binder)
+plt.axvspan(0.96, 1.02, color=(.4,.4,.4, .3)) # by zooming into the figure
+plt.axvline(0.99, color=(1,.2,.2))
+plt.xlabel('$J/\\Gamma$')
+plt.ylabel('Binder Cumulant, $U_2 = <m^2> / <|m|>^2$')
+plt.legend(loc='best', frameon=False)
+
+
+plt.show()