<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 TB Model
We pick the following parameters
Step1:  d Wave
Instantiation
Step2:  Modification
Step3:  MC Driver
Instantiation
Step4:  Modification
<ASSISTANT_TASK:>
Python Code:
Tc_mf = meV_to_K(0.5*250)
print meV_to_K(pi/2.0)
print 1.0/0.89
print cst.physical_constants["Boltzmann constant"]
print '$T_c^{MF} = $', Tc_mf, "K"
T_KT =  meV_to_K(0.1*250)
print r"$T_{KT} = $", T_KT, "K"
Explanation: TB Model
We pick the following parameters:
+ hopping constant $ t= 250$ meV
+ $\Delta = 1.0 t$ so that $T_c^{MF} = 0.5 t$, and so that $\xi_0 \simeq a_0$
+ $g = -0.25$, unitless, so as to match the article's formalism, not the thesis'
+ $J = \dfrac{0.1 t}{0.89}$ so as to set $T_{KT} = 0.1 t$.
This means that we have the following physical properties
End of explanation
T_CST = 0.25
BCS_PARAMS = {"width":4, "chem_potential": 0.0,
              "hopping_constant": T_CST, "J_constant":  0.1 * T_CST / 0.89,
              "g_constant": 0.25, "delta": 1.0 * T_CST, "use_assaad": True,
              "uniform_phase": True, "temperature": 100}
MY_DWAVE_MODEL = DWaveModel(BCS_PARAMS)
print MY_DWAVE_MODEL
Explanation: d Wave
Instantiation
End of explanation
BCS_PARAMS = {"width":20, "use_assaad": True,
              "uniform_phase": True,  "temperature": 1.75*145.0}
MY_DWAVE_MODEL.set_params(BCS_PARAMS)
print MY_DWAVE_MODEL
print "temp: ", K_to_meV(MY_DWAVE_MODEL.temperature), "meV"
Explanation: Modification
End of explanation
BCS_PARAMS = {"width":20, "use_assaad": True,
              "uniform_phase": False,  "temperature": 1.75*145.0}
MY_DWAVE_MODEL.set_params(BCS_PARAMS)
print MY_DWAVE_MODEL._uniform_phase
MC_Params = {"seed": 222315, "intervals": 100,
             "target_snapshots": 15, "observable_list":["correlation_length"]}
MY_DRIVER = MCMCDriver(MY_DWAVE_MODEL, MC_Params)
Explanation: MC Driver
Instantiation
End of explanation
MC_PARAMS_MP = {"intervals": BCS_PARAMS["width"]**2 / 2,
             "target_snapshots": 25,
               "algorithm":"metropolis"}
MC_PARAMS_CLUSTER = {"intervals": 5,
             "target_snapshots": 25,
            "algorithm":"cluster"}
MY_DRIVER.set_params(MC_PARAMS_CLUSTER)
print MY_DWAVE_MODEL._uniform_phase
print MY_DRIVER
print MY_DRIVER.params
#MY_DRIVER.mc_object.set_params({"temperature": 2.0 * 145.0})
#MY_DRIVER.thermalize(20000)
MY_DRIVER.mc_object.set_params({"temperature": 1.05 * 290.0 * 1.05 / 1.12})
MY_DRIVER.thermalize(50)
MY_DRIVER.execute()
result = MY_DRIVER.result
data = result.observable_results["correlation_length"]
print data["length_values"]
print data["correlation_values"]
print result
x_data = np.sqrt(data["length_values"])
y_data = data["correlation_values"]
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
ymin = 0.0
ymax = 1.0
ax.plot(x_data, y_data)
ax.set_ylim([ymin, ymax])
popt, pcov = curve_fit(func_short, x_data[1:], y_data[1:])
print popt
ax.plot(x_data, func_short(x_data, 2.0*popt[0], 2.0*popt[1]))
print "corr length:", 1.0/popt[1]
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
ax.plot(x_data[1:], np.log(y_data[1:]))
ax.plot(x_data[1:], np.log(popt[0]) - popt[1] * x_data[1:])
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
plt.imshow(MY_DRIVER.mc_object.xy_lattice, cmap=plt.cm.hot, interpolation='none')
plt.colorbar() 
#Cf http://matplotlib.org/examples/pylab_examples/quiver_demo.html
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
plt.quiver(np.cos(MY_DRIVER.mc_object.xy_lattice), np.sin(MY_DRIVER.mc_object.xy_lattice))
MY_DRIVER.mc_object.make_wolff_step(np.pi/2.0)
dimension = MY_DRIVER.mc_object.xy_lattice.shape[0]
cluster = np.reshape(MY_DRIVER.mc_object.cluster, (dimension, dimension))
plt.imshow(cluster, cmap=plt.cm.hot, interpolation='none')
plt.colorbar()
print pi
print 5.65226755763 / (2.0 * pi), 3.77251040313 / (2.0 * pi)
neigh = MY_DRIVER.mc_object.lattice.get_neighbors()
results = pickle.load( open( "data_new.txt", "rb" ) )
data = results[0].observable_results["correlation_length"]
datas = {}
temps =np.array([])
for elem in results:
    temps = np.append(temps, elem.bcs_params['temperature'])
temps = np.unique(temps)
for temp in temps:
    datas[temp] = np.array([elem for elem in results if elem.bcs_params['temperature']==temp])
x_datas = {}
y_datas = {}
for temp in temps:
    x_datas[temp] = np.sqrt(datas[temp][0].observable_results["correlation_length"]["length_values"])
    y_datas[temp] = np.zeros((x_datas[temp].size))
    total_sum = 0
    for elem in datas[temp]:
        y_datas[temp] +=\
            elem.observable_results["correlation_length"]["correlation_values"]
    y_datas[temp] /= datas[temp].size
fig, ax = plt.subplots(figsize = (14, 12), dpi=100, frameon=False)
corr_lens = {}
chosen_fun = func_power
for temp in temps[:1]:
    x_data = x_datas[temp]
    y_data = y_datas[temp]
    print temp
    ax.plot(x_data, y_data, label=str(temp))
    #popt, pcov = curve_fit(func, x_data[0:], y_data[0:])
    popt, pcov = curve_fit(chosen_fun, x_data[1:], y_data[1:])
    print "temp: ", temp, "params: ", popt, r"$\eta$: ", -popt[1]
    corr_lens[temp] = 1.0/popt[1]
    ax.plot(x_data, chosen_fun(x_data, popt[0], popt[1]))
chosen_fun = func_full
for temp in temps[1:]:
    x_data = x_datas[temp]
    y_data = y_datas[temp]
    print temp
    ax.plot(x_data, y_data, label=str(temp))
    #popt, pcov = curve_fit(func, x_data[0:], y_data[0:])
    popt, pcov = curve_fit(chosen_fun, x_data[1:], y_data[1:])
    print "temp: ", temp, "params: ", popt, "length: ", 1.0/popt[1], 1.0/popt[-1]
    corr_lens[temp] = 1.0/popt[1]
    ax.plot(x_data, chosen_fun(x_data, popt[0], popt[1], popt[2]))
    
ax.legend()
plt.savefig("transition.pdf")
fig, ax = plt.subplots(figsize = (14, 12), dpi=100, frameon=False)
x_es = np.sort(np.array(corr_lens.keys()))
y_es = np.array([corr_lens[elem] for elem in x_es])
ax.plot(x_es[1:], y_es[1:])
ax.grid(True)
plt.savefig("corr_length.pdf")
popt, pcov = curve_fit(func_exponent, x_es[1:], y_es[1:])
nu = popt[1]
print "nu", popt[1]
fig, ax = plt.subplots(figsize = (14, 12), dpi=100, frameon=False)
ax.plot(x_es[1:], y_es[1:])
ax.plot(x_es[1:-1], func_exponent(x_es[1:-1], popt[0], popt[1]))
ax.grid(True)
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
temp = 275.0
x_data = x_datas[temp]
y_data = y_datas[temp]
print x_data.shape
print y_data.shape
print x_es.shape
print data["length_values"].shape
l_values = y_data
popt, pcov = curve_fit(func, x_data, y_data)
ax.plot(np.sqrt(data["length_values"]), np.log(l_values - popt[2]))
ax.plot(np.sqrt(data["length_values"]), np.log(popt[0]) - popt[1] * np.sqrt(data["length_values"]))
results = pickle.load( open( "dos_alltemps.txt", "rb" ) )
datas = {}
temps =np.array([])
for elem in results:
    temps = np.append(temps, elem.bcs_params['temperature'])
temps = np.unique(temps)
for temp in temps:
    datas[temp] = np.array([elem for elem in results if elem.bcs_params['temperature']==temp])
print temps
print datas[275.0][0]
x_datas = {}
y_datas = {}
for temp in temps:
    x_datas[temp] = datas[temp][0].observable_results["DOS"]["omega_mesh"]
    y_datas[temp] = np.zeros((x_datas[temp].size))
    total_sum = 0
    for elem in datas[temp]:
        y_datas[temp] +=\
            elem.observable_results["DOS"]["DOS_values"]
    y_datas[temp] /= datas[temp].size
fig, ax = plt.subplots(figsize = (8, 14), dpi=100, frameon=False)
#for i in range(len(temps)):
selected_temps = [0,1,  2, 3, 4, 6, 8]
for i in range(len(selected_temps)):
    temp = temps[selected_temps[i]]
    x_data = x_datas[temp]
    y_data = y_datas[temp] + i * 0.7
    ax.plot(x_data, y_data, label=(r'T={:3.2f}$T_{{KT}}$').format(temp/T_KT))
    
ax.legend()
plt.savefig("all_dos.pdf")
Explanation: Modification
End of explanation
<END_TASK>