In [1]:
%install_ext https://raw.githubusercontent.com/rasbt/python_reference/master/ipython_magic/watermark.py
%load_ext watermark
In [2]:
import numpy as np
from numpy import linalg as LA
from numpy import array
from numpy import pi
import matplotlib.pyplot as plt
import sys
import math
from scipy.optimize import curve_fit
import pickle
if "/home/dominique/Code/PG/Source" not in sys.path:
sys.path.append("/home/dominique/Code/PG/Source")
import phase_fluctuations as PF
from phase_fluctuations import TbModel, TbParams, SWaveModel, DWaveModel
from MCMC import MCMCDriver
import scipy.constants as cst
def K_to_meV(in_temp):
return cst.physical_constants["Boltzmann constant in eV/K"][0] * in_temp * 1000.0
def meV_to_K(in_temp):
return in_temp / 1000.0 / cst.physical_constants["Boltzmann constant in eV/K"][0]
def func(x, a, b, c):
return a * np.exp(-b * x) + c
In [3]:
%load_ext autoreload
%autoreload 2
In [4]:
%who
%matplotlib inline
In [5]:
%watermark -a "Dominique" -d -t -u -v -h -m -g
We pick the following parameters:
This means that we have the following physical properties
In [13]:
Tc_mf = meV_to_K(0.5*250)
print r"$T_c^{MF} = $", Tc_mf, "K"
print r"$T_{KT} = $", Tc_mf/10.0, "K"
In [6]:
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)
In [7]:
print MY_DWAVE_MODEL
In [8]:
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"
In [9]:
dos_values = np.real(MY_DWAVE_MODEL.get_dos())
In [10]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
x_values = MY_DWAVE_MODEL.lattice.omega_mesh / T_CST
xmin = np.amin(x_values)
xmax = np.amax(x_values)
ax.set_xlim([xmin, xmax])
ax.set_xticks(np.linspace(int(xmin), int(xmax), 13,endpoint=True))
ax.plot(x_values, dos_values)
Out[10]:
In [26]:
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
In [27]:
MC_Params = {"seed": 222315, "intervals": 100,
"target_snapshots": 10, "observable_list":["correlation_length"]}
MY_DRIVER = MCMCDriver(MY_DWAVE_MODEL, MC_Params)
In [28]:
print MY_DRIVER.mc_object.xy_lattice
In [31]:
MC_PARAMS = {"intervals": 100,
"target_snapshots": 25}
MY_DRIVER.set_params(MC_PARAMS)
print MY_DWAVE_MODEL._uniform_phase
In [32]:
print MY_DRIVER
print MY_DRIVER.params
In [33]:
MY_DRIVER.mc_object.set_params({"temperature": 0.0000000000001 * 145.0})
MY_DRIVER.thermalize(20000)
In [50]:
MY_DRIVER.thermalize(1000000)
print MY_DRIVER.mc_object.xy_lattice
In [51]:
MY_DRIVER.execute()
In [52]:
result = MY_DRIVER.result
data = result.observable_results["correlation_length"]
print data["length_values"].size
print data["correlation_values"]
print "\n"
print result
x_data = np.sqrt(data["length_values"])
y_data = data["correlation_values"]
In [53]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
ax.plot(x_data, y_data)
popt, pcov = curve_fit(func, x_data, y_data)
print popt
print "corr length:", 1.0/popt[1]
ax.plot(x_data, func(x_data, popt[0], popt[1], popt[2]))
Out[53]:
In [ ]:
In [ ]:
In [60]:
results = pickle.load( open( "result_corr_dwave_new.txt", "rb" ) )
results = np.append(results, pickle.load( open( "result_corr_dwave_new2.txt", "rb" ) ))
#results = np.append(results, pickle.load( open( "result_corr_dwave_4.txt", "rb" ) ))
#results = np.append(results, pickle.load( open( "result_corr_dwave_5.txt", "rb" ) ))
#results = np.append(results, pickle.load( open( "result_corr_dwave_6.txt", "rb" ) ))
data = results[0].observable_results["correlation_length"]
l_values = []
for elem in data['correlation_values']:
l_values.append(np.average(elem))
print data["length_values"].size
print len(l_values)
In [61]:
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[temps[0]].size
In [62]:
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
#np.array([np.average(zob) for zob in elem.observable_results["correlation_length"]["correlation_values"]])
In [65]:
fig, ax = plt.subplots(figsize = (14, 12), dpi=100, frameon=False)
corr_lens = {}
for temp in temps:
x_data = x_datas[temp]
y_data = y_datas[temp]
ax.plot(x_data, y_data, label=str(temp))
popt, pcov = curve_fit(func, x_data, y_data)
print "temp: ", temp, "params: ", popt, "length: ", 1.0/popt[1]
corr_lens[temp] = 1.0/popt[1]
ax.plot(x_data, func(x_data, popt[0], popt[1], popt[2]))
ax.legend()
Out[65]:
In [70]:
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, y_es)
ax.grid(True)
plt.savefig("corr_length.pdf")
In [ ]:
a = np.array([1.125, 1.25, 1.5, 1.75, 2.0, 3.0, 4.0, 5.0])
In [199]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
for temp in temps:
x_data = x_datas[temp]
y_data = y_datas[temp]
ax.plot(x_data, y_data, label=str(temp))
popt, pcov = curve_fit(func, x_data, y_data)
print "temp: ", temp, "params: ", popt, "length: ", 1.0/popt[1]
ax.plot(x_data, func(x_data, popt[0], popt[1], popt[2]))
ax.legend()
Out[199]:
In [123]:
print results
In [121]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
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"]))
Out[121]:
In [454]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
x_values = MY_DRIVER.mc_object.lattice.omega_mesh / T_CST
xmin = np.amin(x_values)
xmax = np.amax(x_values)
ax.set_xlim([xmin, xmax])
ax.set_xticks(np.linspace(int(xmin), int(xmax), 13,endpoint=True))
ax.plot(x_values, dos_values)
Out[454]:
In [562]:
print MY_DRIVER.observable_results
In [383]:
MY_DRIVER.mc_object.temperature = 1.0 * 145.0
MY_DRIVER.thermalize(20000)
In [385]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
x_axis_ampl = 4.5
xmin = -x_axis_ampl
xmax = x_axis_ampl
ax.set_xlim([xmin, xmax])
#ax.set_ylim([0.0, 1.0])
ax.set_xticks(np.linspace(int(xmin), int(xmax), 9,endpoint=True))
#ax.set_yticks(np.linspace(xmin, xmax, 5,endpoint=True))
#ax.set_xticklabels(['%1.1f' %elem for elem in np.linspace(xmin, xmax,5,endpoint=True)])
#ax.set_yticklabels(['%1.1f' %elem for elem in np.linspace(xmin, xmax,5,endpoint=True)])
MY_DRIVER.execute()
dos_values = np.real(MY_DRIVER.observable_results["DOS"])
ax.plot(MY_DRIVER.mc_object.lattice.omega_mesh / T_CST, dos_values)
Out[385]:
In [ ]:
In [ ]:
1500/ 450
In [ ]:
In [ ]:
In [582]:
#300K
temps = [15, 113, 150, 155, 300, 450, 600, 800, 1200]
names = ['./result_dwave_' + str(temp)+'.txt' for temp in temps]
input_files = [np.loadtxt(name) for name in names]
In [585]:
print len(input_files[3])
In [609]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
x_axis_ampl = 6.0 * T_CST
nb_ticks = 1001
xvalues = np.linspace(-x_axis_ampl, x_axis_ampl, nb_ticks, endpoint=True)
#ymin = 0.0
#ymax = 0.5
xmin = -x_axis_ampl
xmax = x_axis_ampl
ax.set_xlim([xmin, xmax])
#ax.set_ylim([0.0, 2.0])
ax.set_xticks(np.linspace(int(xmin), int(xmax), 9,endpoint=True))
#ax.set_yticks(np.linspace(xmin, xmax, 5,endpoint=True))
#ax.set_xticklabels(['%1.1f' %elem for elem in np.linspace(xmin, xmax,5,endpoint=True)])
#ax.set_yticklabels(['%1.1f' %elem for elem in np.linspace(xmin, xmax,5,endpoint=True)])
for i in range(len(input_files)):
ax.plot(xvalues/T_CST, input_files[i], ls = '-', label = str(temps[i]))
ax.legend(loc=2)
plt.savefig("DOS.pdf")
In [297]:
#200K
out = np.loadtxt('result_dwave_750.txt')
In [298]:
fig, ax = plt.subplots(figsize = (10, 8), dpi=100, frameon=False)
x_axis_ampl = 4.2
nb_ticks = 1001
xvalues = np.linspace(-x_axis_ampl, x_axis_ampl, nb_ticks, endpoint=True)
#ymin = 0.0
#ymax = 0.5
xmin = -x_axis_ampl
xmax = x_axis_ampl
ax.set_xlim([xmin, xmax])
#ax.set_ylim([0.0, 1.0])
ax.set_xticks(np.linspace(int(xmin), int(xmax), 9,endpoint=True))
#ax.set_yticks(np.linspace(xmin, xmax, 5,endpoint=True))
#ax.set_xticklabels(['%1.1f' %elem for elem in np.linspace(xmin, xmax,5,endpoint=True)])
#ax.set_yticklabels(['%1.1f' %elem for elem in np.linspace(xmin, xmax,5,endpoint=True)])
ax.plot(xvalues/0.25, out)
plt.savefig("DOS_dwave_T_KT.pdf")
In [227]:
150.0*1.75
Out[227]:
In [303]:
print cst.physical_constants["Boltzmann constant in eV/K"][0]
In [ ]: