notebook.community

Edit and run



In [1]:

    
!ipython locate profile









    



/home/dominique/.ipython/profile_default



In [2]:

    
%install_ext https://raw.githubusercontent.com/rasbt/python_reference/master/ipython_magic/watermark.py
%load_ext watermark









    



Installed watermark.py. To use it, type:
  %load_ext watermark



In [3]:

    
import numpy as np
from numpy import linalg as LA
from numpy import array
from numpy import pi
import matplotlib.pyplot as plt



In [4]:

    
fig_width = 3.39
golden_mean = (np.sqrt(5)-1.0)/2.0    # Aesthetic ratio
fig_height = fig_width*golden_mean # height in inches
MAX_HEIGHT_INCHES = 8.0
if fig_height > MAX_HEIGHT_INCHES:
    print("WARNING: fig_height too large:" + fig_height + 
          "so will reduce to" + MAX_HEIGHT_INCHES + "inches.")
    fig_height = MAX_HEIGHT_INCHES
fig, ax = plt.subplots(figsize = (fig_width, fig_height), dpi=400, frameon=True)
small_tick_size = 8
small_label_size = 8



In [5]:

    
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
def func_short(x, a, b):
    return a * np.exp(-b * x)
def func_combined(x, a, b, d):
    return a * np.exp(-b * x) / np.power(x, d)
def func_power(x, a, b):
    return a * np.power(x, b)
def func_exponent(x, a, b):
    #return a * np.power(x - 290.112982064, b)
    return a * np.exp(-np.power(x - 290.112982064, b))
def func_full(x, a, b, c):
    return a * np.exp(-np.abs(b) * x) / np.power(x, c)



In [6]:

    
%load_ext autoreload
%autoreload 2



In [7]:

    
%who
%matplotlib inline









    



DWaveModel	 K_to_meV	 LA	 MAX_HEIGHT_INCHES	 MCMCDriver	 PF	 SWaveModel	 TbModel	 TbParams	 
array	 ax	 cst	 curve_fit	 fig	 fig_height	 fig_width	 func	 func_combined	 
func_exponent	 func_full	 func_power	 func_short	 golden_mean	 math	 meV_to_K	 np	 pi	 
pickle	 plt	 small_label_size	 small_tick_size	 sys



In [8]:

    
%watermark -a "Dominique" -d -t -u -v -h -m -g









    



Dominique Last updated: 04/04/2015 22:43:46 

CPython 2.7.8
IPython 3.0.0

compiler   : GCC 4.9.1
system     : Linux
release    : 3.16.0-33-generic
machine    : x86_64
processor  : x86_64
CPU cores  : 4
interpreter: 64bit
host name  : Olympe
Git hash   : bdd56493d7f0731ac9f7ce6fe81fcd415e283d4d

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



In [12]:

    
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"]









    



18.2283362633
1.12359550562
(1.3806488e-23, 'J K^-1', 1.3e-29)



In [13]:

    
print '$T_c^{MF} = $', Tc_mf, "K"
T_KT =  meV_to_K(0.1*250)
print r"$T_{KT} = $", T_KT, "K"









    



$T_c^{MF} = $ 1450.56491032 K
$T_{KT} = $ 290.112982064 K

Instantiation



In [14]:

    
T_CST = 0.25
MY_PARAMS = {"width":10, "chem_potential": 0.0,
                 "hopping_constant": T_CST, "J_constant": 0.1 * T_CST / 0.89,
                 "delta": 1.0 * T_CST,
                 "use_assaad": False, "broadening_delta": 0.01 * T_CST}
MY_MODEL = TbModel(MY_PARAMS)

TB Model

Modification



In [15]:

    
TB_PARAMS = {"width":14, "use_assaad": True}
MY_MODEL.set_params(TB_PARAMS)
print MY_MODEL









    



Class: <class 'phase_fluctuations.TbModel'>
broadening delta: 0.0025 eV
use Assad: True
is up to date: False

Lattice: 
Class: <class 'phase_fluctuations.TbParams'>
Width: 14
Chemical potential: 0.0 eV
Hopping constant: 0.25 eV
Number of sites: 
196

A(k, $\omega$) Computation



In [13]:

    
all_test_vectors = []
nb_samples = MY_MODEL.lattice.width
for elem in np.arange(nb_samples):
    print "handling", elem
    k_in = elem * 1.0 * np.pi * np.array([1.0, 1.0]) / nb_samples
    test2_vector = MY_MODEL.get_vector_spectral_function_point(k_in,
                                                               MY_MODEL.lattice.omega_mesh)
    all_test_vectors.append(test2_vector)









    



handling 0
handling 1
handling 2
handling 3
handling 4
handling 5
handling 6
handling 7
handling 8
handling 9
handling 10
handling 11
handling 12
handling 13

parallel code version



In [23]:

    
def get_spectral_function(x):
    import numpy as np
    import os
    print "process", os.getpid(), "handling", x
    nb_samples = MY_MODEL.lattice.width
    k_in = x * np.pi * np.array([1.0, 1.0]) / nb_samples
    return MY_MODEL.get_vector_spectral_function_point(k_in, MY_MODEL.lattice.omega_mesh)



In [ ]:

    
import multiprocessing
nb_samples = MY_MODEL.lattice.width
multiprocessing.cpu_count()
pool = multiprocessing.Pool(processes=3)
all_test_vectors = pool.map(get_spectral_function, np.arange(nb_samples))



In [19]:

    
array_coords = []
z_values = []
for elem in np.arange(nb_samples):
    k_in = elem * np.pi * np.array([1.0, 0.0]) / nb_samples
    to_add = [[k_in[0], energy/T_CST] for energy in MY_MODEL.lattice.omega_mesh]
    array_coords.extend(to_add)
    z_values.extend(all_test_vectors[elem])# / np.amax(all_test_vectors[elem]))
    
array_coords = np.array(array_coords)
z_values = np.array(z_values)



In [20]:

    
from scipy.interpolate import griddata
visu_factor = 3
fig, ax = plt.subplots(figsize = (fig_width*visu_factor, fig_height*visu_factor), dpi=400, frameon=False)

xmin = 0.0
xmax = np.pi
ymin = -5.0
ymax = 5.0
#grid_x = np.linspace(xmin, xmax, num = 200)
#grid_y = np.linspace(ymin, ymax, num = 2000)
grid_x, grid_y = np.mgrid[xmin:xmax:100j, ymin:ymax:200j]
grid_z1 = griddata(array_coords, z_values, (grid_x, grid_y), method='nearest', fill_value=0.0)
p = ax.pcolormesh(grid_x, grid_y, grid_z1, cmap=plt.get_cmap('Greys'),rasterized=True)

#for index, theta in np.ndenumerate(thetas):
#    ax.plot(theta, resonance_omega_values[index], marker='o', color='k')

cb = fig.colorbar(p, ax=ax)
ax.plot(grid_x, -2.0 * 2.0 * T_CST * np.cos(grid_x) / T_CST)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
plt.savefig("spectral_function.pdf")



In [39]:

    
import contextlib
@contextlib.contextmanager
def printoptions(*args, **kwargs):
    """Define some temporary pretty print options
    """
    original = np.get_printoptions()
    np.set_printoptions(*args, **kwargs)
    yield
    np.set_printoptions(**original)



In [40]:

    
with printoptions(precision=4, suppress=True, threshold=np.nan):
    print MY_MODEL._eigen_values
    for i in np.arange(MY_MODEL.lattice.width):
        print MY_MODEL._unitary_matrix[i]









    



[-0.5 -0.5 -0.  -0.   0.   0.   0.5  0.5]
[-0.5000+0.j  0.0000+0.j  0.7071+0.j  0.0000+0.j  0.0000+0.j  0.0000+0.j
 -0.5000+0.j  0.0000+0.j]
[-0.5000+0.j  0.0000+0.j  0.0000+0.j  0.0000+0.j  0.7071+0.j  0.0000+0.j
  0.5000+0.j  0.0000+0.j]



In [89]:

    
k_in = 1.0 * np.pi * np.array([1.0, 1.0])
test = MY_MODEL.get_vector_spectral_function_point(k_in,
                                                   MY_MODEL.lattice.omega_mesh)



In [90]:

    
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])
l_xticks = np.linspace(int(xmin), int(xmax), 9,endpoint=True)
ax.set_xticks(l_xticks)
ax.set_xticklabels(['%1.1f'  %elem for elem in l_xticks])
#ax.set_xticklabels([r'$\psi$' for elem in l_xticks])
#ax.set_yticks(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(MY_MODEL.lattice.omega_mesh/T_CST, test)
#ax.plot(MY_MODEL.lattice.omega_mesh/T_CST, test_elem)









    Out[90]:





[<matplotlib.lines.Line2D at 0x7f181ca85810>]

d Wave

Instantiation



In [18]:

    
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 [19]:

    
print MY_DWAVE_MODEL









    



Class: <class 'phase_fluctuations.DWaveModel'>
Seed: 1234567890
Uniform phase: True
Temperature: 100 K

Base class:
Class: <class 'phase_fluctuations.DWaveModel'>
broadening delta: 0.0025 eV
use Assad: True
is up to date: False

Lattice: 
Class: <class 'phase_fluctuations.PairingParams'>
J constant: 0.0280898876404 eV
U constant: -0.825 eV
g constant: 0.25 eV
delta: 0.25 eV

Base class:
Class: <class 'phase_fluctuations.PairingParams'>
Width: 4
Chemical potential: 0.0 eV
Hopping constant: 0.25 eV
Number of sites: 
16

Modification



In [20]:

    
BCS_PARAMS = {"width":12, "use_assaad": True,
              "uniform_phase": True,  "temperature": 1.75*145.0,
              "delta": 1.0 * T_CST}
MY_DWAVE_MODEL.set_params(BCS_PARAMS)
print MY_DWAVE_MODEL
print "temp: ", K_to_meV(MY_DWAVE_MODEL.temperature), "meV"









    



Class: <class 'phase_fluctuations.DWaveModel'>
Seed: 1234567890
Uniform phase: True
Temperature: 253.75 K

Base class:
Class: <class 'phase_fluctuations.DWaveModel'>
broadening delta: 0.0025 eV
use Assad: True
is up to date: False

Lattice: 
Class: <class 'phase_fluctuations.PairingParams'>
J constant: 0.0280898876404 eV
U constant: -0.825 eV
g constant: 0.25 eV
delta: 0.25 eV

Base class:
Class: <class 'phase_fluctuations.PairingParams'>
Width: 12
Chemical potential: 0.0 eV
Hopping constant: 0.25 eV
Number of sites: 
144
temp:  21.866480965 meV



In [129]:

    
k_in = 9.0 * np.pi * np.array([1.0, 1.0]) / MY_DWAVE_MODEL.lattice.width
test = MY_DWAVE_MODEL.get_vector_spectral_function_point(k_in,
                                                         MY_DWAVE_MODEL.lattice.omega_mesh)



In [130]:

    
test_elem = []
for omega in MY_DWAVE_MODEL.lattice.omega_mesh:
    test_elem.append(MY_DWAVE_MODEL.get_spectral_function_point(k_in, omega))

debug



In [128]:

    
test_index = 105
omega_test = MY_DWAVE_MODEL.lattice.omega_mesh[test_index]
omega = MY_DWAVE_MODEL.lattice.omega_mesh[test_index]

a = MY_DWAVE_MODEL.get_vector_spectral_function_point(k_in, np.array([omega_test]))
c = MY_DWAVE_MODEL.get_spectral_function_point(k_in, omega_test)
print np.allclose(a, c)









    



True



In [131]:

    
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])
l_xticks = np.linspace(int(xmin), int(xmax), 9,endpoint=True)
ax.set_xticks(l_xticks)
ax.set_xticklabels(['%1.1f'  %elem for elem in l_xticks])
#ax.set_xticklabels([r'$\psi$' for elem in l_xticks])
#ax.set_yticks(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(MY_DWAVE_MODEL.lattice.omega_mesh/T_CST, test, label = "vector")
ax.plot(MY_DWAVE_MODEL.lattice.omega_mesh/T_CST, test_elem, label = "single")
ax.legend()









    Out[131]:





<matplotlib.legend.Legend at 0x7f4514d87510>



In [132]:

    
test2 = []
for omega in MY_DWAVE_MODEL.lattice.omega_mesh:
    test2.append(MY_DWAVE_MODEL.get_spectral_function_point(np.pi * np.array([1.0, 0.0]),
                                                            omega))



In [133]:

    
k_in = np.pi * np.array([1.0, 0.0])
test2_vector = MY_DWAVE_MODEL.get_vector_spectral_function_point(k_in,
                                                         MY_DWAVE_MODEL.lattice.omega_mesh)



In [134]:

    
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])
l_xticks = np.linspace(int(xmin), int(xmax), 9,endpoint=True)
ax.set_xticks(l_xticks)
ax.set_xticklabels(['%1.1f'  %elem for elem in l_xticks])
#ax.set_xticklabels([r'$\psi$' for elem in l_xticks])
#ax.set_yticks(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(MY_DWAVE_MODEL.lattice.omega_mesh/T_CST, test2_vector)
ax.plot(MY_DWAVE_MODEL.lattice.omega_mesh/T_CST, test2)









    Out[134]:





[<matplotlib.lines.Line2D at 0x7f4514649c10>]

Scan BZ



In [21]:

    
all_test_vectors = []
nb_samples = MY_DWAVE_MODEL.lattice.width
for elem in np.arange(nb_samples):
    print "handling", elem
    k_in = elem * np.pi * np.array([1.0, 1.0]) / nb_samples
    test2_vector = MY_DWAVE_MODEL.get_vector_spectral_function_point(k_in,
                                                                     MY_DWAVE_MODEL.lattice.omega_mesh)
    all_test_vectors.append(test2_vector)









    



handling 0
handling 1
handling 2
handling 3
handling 4
handling 5
handling 6
handling 7
handling 8
handling 9
handling 10
handling 11



In [22]:

    
array_coords = []
z_values = []
for elem in np.arange(nb_samples):
    k_in = elem * np.pi * np.array([1.0, 0.0]) / nb_samples
    to_add = [[k_in[0], energy/T_CST] for energy in MY_DWAVE_MODEL.lattice.omega_mesh]
    array_coords.extend(to_add)
    z_values.extend(all_test_vectors[elem])# / np.amax(all_test_vectors[elem]))
    
array_coords = np.array(array_coords)
z_values = np.array(z_values)



In [23]:

    
from scipy.interpolate import griddata
visu_factor = 3
fig, ax = plt.subplots(figsize = (fig_width*visu_factor, fig_height*visu_factor), dpi=400, frameon=False)

xmin = 0.0
xmax = np.pi
ymin = -5.0
ymax = 5.0
#grid_x = np.linspace(xmin, xmax, num = 200)
#grid_y = np.linspace(ymin, ymax, num = 2000)
grid_x, grid_y = np.mgrid[xmin:xmax:100j, ymin:ymax:200j]
grid_z1 = griddata(array_coords, z_values, (grid_x, grid_y), method='nearest', fill_value=0.0)
p = ax.pcolormesh(grid_x, grid_y, grid_z1, cmap=plt.get_cmap('Greys'),rasterized=True)

#for index, theta in np.ndenumerate(thetas):
#    ax.plot(theta, resonance_omega_values[index], marker='o', color='k')

cb = fig.colorbar(p, ax=ax)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])









    Out[23]:





(-5.0, 5.0)



In [24]:

    
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])
l_xticks = np.linspace(int(xmin), int(xmax), 9,endpoint=True)
ax.set_xticks(l_xticks)
ax.set_xticklabels(['%1.1f'  %elem for elem in l_xticks])
#ax.set_xticklabels([r'$\psi$' for elem in l_xticks])
#ax.set_yticks(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(MY_DWAVE_MODEL.lattice.omega_mesh/T_CST,
        MY_DWAVE_MODEL.get_vector_spectral_function_point(np.array([0.6, 0.6]),
                                                         MY_DWAVE_MODEL.lattice.omega_mesh))









    Out[24]:





[<matplotlib.lines.Line2D at 0x7fadcdcf6590>]

Scan BZ: gap



In [25]:

    
all_test_vectors = []
nb_samples = MY_DWAVE_MODEL.lattice.width
for elem in np.arange(nb_samples):
    print "handling", elem
    k_in = elem * np.pi * np.array([1.0, 0.0]) / nb_samples
    test2_vector = MY_DWAVE_MODEL.get_vector_spectral_function_point(k_in,
                                                                     MY_DWAVE_MODEL.lattice.omega_mesh)
    all_test_vectors.append(test2_vector)









    



handling 0
handling 1
handling 2
handling 3
handling 4
handling 5
handling 6
handling 7
handling 8
handling 9
handling 10
handling 11



In [26]:

    
array_coords = []
z_values = []
for elem in np.arange(nb_samples):
    k_in = elem * np.pi * np.array([1.0, 0.0]) / nb_samples
    to_add = [[k_in[0], energy/T_CST] for energy in MY_DWAVE_MODEL.lattice.omega_mesh]
    array_coords.extend(to_add)
    z_values.extend(all_test_vectors[elem])# / np.amax(all_test_vectors[elem]))
    
array_coords = np.array(array_coords)
z_values = np.array(z_values)



In [27]:

    
from scipy.interpolate import griddata
visu_factor = 3
fig, ax = plt.subplots(figsize = (fig_width*visu_factor, fig_height*visu_factor), dpi=400, frameon=False)

xmin = 0.0
xmax = np.pi
ymin = -5.0
ymax = 5.0
#grid_x = np.linspace(xmin, xmax, num = 200)
#grid_y = np.linspace(ymin, ymax, num = 2000)
grid_x, grid_y = np.mgrid[xmin:xmax:100j, ymin:ymax:200j]
grid_z1 = griddata(array_coords, z_values, (grid_x, grid_y), method='nearest', fill_value=0.0)
p = ax.pcolormesh(grid_x, grid_y, grid_z1, cmap=plt.get_cmap('Greys'),rasterized=True)

#for index, theta in np.ndenumerate(thetas):
#    ax.plot(theta, resonance_omega_values[index], marker='o', color='k')

cb = fig.colorbar(p, ax=ax)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])









    Out[27]:





(-5.0, 5.0)



In [ ]: