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



    











In [120]:



    


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



    


%load_ext autoreload
%autoreload 2



    











In [6]:



    


%who
%matplotlib inline



    


















    






DWaveModel	 K_to_meV	 LA	 MCMCDriver	 PF	 SWaveModel	 TbModel	 TbParams	 array	 
cst	 curve_fit	 func	 func_short	 math	 meV_to_K	 np	 pi	 pickle	 
plt	 sys	 

















In [7]:



    


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



    


















    






Dominique Last updated: 28/03/2015 09:41:28 

CPython 2.7.8
IPython 3.0.0

compiler   : GCC 4.9.1
system     : Linux
release    : 3.16.0-31-generic
machine    : x86_64
processor  : x86_64
CPU cores  : 4
interpreter: 64bit
host name  : Everest
Git hash   : aab789eb0c93256e31988777dbee4523b69d0a31

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



    


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



    


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

d Wave

Instantiation





In [10]:



    


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



    


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



    


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"



    


















    






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: 20
Chemical potential: 0.0 eV
Hopping constant: 0.25 eV
Number of sites: 
400
temp:  21.866480965 meV

MC Driver

Instantiation





In [13]:



    


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



    


















    






False

















In [14]:



    


MC_Params = {"seed": 222315, "intervals": 100,
             "target_snapshots": 15, "observable_list":["correlation_length"]}
MY_DRIVER = MCMCDriver(MY_DWAVE_MODEL, MC_Params)

Modification





In [15]:



    


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



    


















    






False

















In [16]:



    


print MY_DRIVER
print MY_DRIVER.params



    


















    






Class: <class 'MCMC.MCMCDriver'>
Seed: 222315
Intervals: 1
Algorithm: cluster
Target snapshots: 25
Observable list: ['correlation_length']

MC Object:
Class: <class 'phase_fluctuations.DWaveModel'>
Seed: 1234567890
Uniform phase: False
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: 20
Chemical potential: 0.0 eV
Hopping constant: 0.25 eV
Number of sites: 
400

Results: None
{'intervals': 5, 'algorithm': 'cluster', 'seed': 222315, 'target_snapshots': 25, 'observable_list': ['correlation_length']}

















In [17]:



    


#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)



    











In [18]:



    


MY_DRIVER.execute()



    











In [19]:



    


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



    


















    






[   0.    1.    2.    4.    5.    8.    9.   10.   13.   16.   17.   18.
   20.   25.   26.   29.   32.   34.   36.   37.   40.   41.   45.   49.
   50.   52.   53.   58.   61.   64.   65.   68.   72.   73.   74.   80.
   81.   82.   85.   89.   90.   97.   98.  100.  101.  104.  106.  109.
  113.  116.  117.  125.  128.  130.  136.  145.  149.  162.  164.  181.
  200.]
[ 1.          0.69324935  0.63334422  0.59111875  0.57936019  0.54655568
  0.544009    0.53694147  0.51479796  0.5124851   0.51040245  0.5003111
  0.49467697  0.48429538  0.48382708  0.47891251  0.48361157  0.47139874
  0.46490092  0.46316897  0.45960648  0.46861285  0.45864352  0.45287763
  0.45085844  0.45388541  0.44508958  0.44321541  0.45309457  0.43994063
  0.44184925  0.43806486  0.44956742  0.44424431  0.44256827  0.43786972
  0.44056534  0.44022768  0.43691485  0.43104077  0.43922442  0.43290734
  0.42341455  0.42826259  0.43728858  0.44083704  0.43193061  0.43727624
  0.41788749  0.43282958  0.42731785  0.42490106  0.42069487  0.42410028
  0.42491865  0.42336991  0.42307732  0.41520815  0.42534857  0.4164745
  0.42731171]
Class: <class 'MCMC.ResultContainer'>
bcs_params: {'hopping_constant': 0.25, 'uniform_phase': False, 'use_assaad': True, 'delta': 0.25, 'width': 20, 'J_constant': 0.02808988764044944, 'g_constant': 0.25, 'chem_potential': 0.0, 'temperature': 285.46875}
mc_params: {'intervals': 5, 'algorithm': 'cluster', 'seed': 222315, 'target_snapshots': 25, 'observable_list': ['correlation_length']}
Code version: {'id': 'aab789eb0c93256e31988777dbee4523b69d0a31', 'time': 'Fri, 27 Mar 2015 22:12'}
Observable list: ['correlation_length']

















In [20]:



    


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]



    


















    






[ 0.59743421  0.03282085]
corr length: 30.4684374262










    
























In [21]:



    


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:])



    


















    
Out[21]:








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










    
























In [22]:



    


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()



    


















    
Out[22]:








<matplotlib.colorbar.Colorbar instance at 0x7f511534cb48>










    
























In [23]:



    


#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))



    


















    
Out[23]:








<matplotlib.quiver.Quiver at 0x7f512c100b50>










    
























In [24]:



    


MY_DRIVER.mc_object.make_wolff_step(np.pi/2.0)



    











In [25]:



    


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()



    


















    
Out[25]:








<matplotlib.colorbar.Colorbar instance at 0x7f51150e4050>










    
























In [26]:



    


print pi
print 5.65226755763 / (2.0 * pi), 3.77251040313 / (2.0 * pi)



    


















    






3.14159265359
0.899586321475 0.600413678524

















In [27]:



    


neigh = MY_DRIVER.mc_object.lattice.get_neighbors()



    











In [40]:



    


results = pickle.load( open( "data_new.txt", "rb" ) )
data = results[0].observable_results["correlation_length"]



    











In [31]:



    


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



    











In [33]:



    


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



    











In [127]:



    


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



    


















    






275.0
temp:  275.0 params:  [ 0.70093587 -0.15788423] $\eta$:  0.157884234645
320.0
temp: 









    






/usr/local/lib/python2.7/dist-packages/IPython/kernel/__main__.py:25: RuntimeWarning: divide by zero encountered in power
/usr/local/lib/python2.7/dist-packages/IPython/kernel/__main__.py:30: RuntimeWarning: divide by zero encountered in divide










    






 320.0 params:  [ 0.52874106  0.02796279  0.55519556] length:  35.7618114405 1.80116714994
360.0
temp:  360.0 params:  [ 0.4456464   0.35189884  0.31269504] length:  2.84172576118 3.19800407285
395.0
temp:  395.0 params:  [ 0.42949439  0.50909813  0.26239391] length:  1.96425784109 3.81106406145
430.0
temp:  430.0 params:  [ 0.39402915  0.59545075  0.26372577] length:  1.67940002376 3.79181759901
470.0
temp:  470.0 params:  [ 0.34413265  0.57988514  0.399743  ] length:  1.72447944138 2.50160729759
490.0
temp:  490.0 params:  [ 0.33145845  0.66841071  0.37406881] length:  1.49608614668 2.67330494891
550.0
temp:  550.0 params:  [ 0.33158425  0.86752044  0.32816207] length:  1.15271059461 3.04727477154
575.0
temp:  575.0 params:  [ 0.3680106   1.03963076  0.2309973 ] length:  0.961879967852 4.32905491086










    
























In [126]:



    


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



    


















    
























In [125]:



    


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)



    


















    






nu 1.0










    
























In [93]:



    


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



    


















    






(135,)
(135,)
(9,)
(135,)










    






/usr/local/lib/python2.7/dist-packages/IPython/kernel/__main__.py:13: RuntimeWarning: invalid value encountered in log










    
Out[93]:








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










    
























In [128]:



    


results = pickle.load( open( "dos_alltemps.txt", "rb" ) )



    











In [129]:



    


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



    











In [132]:



    


print temps
print datas[275.0][0]



    


















    






[ 275.  320.  360.  395.  430.  470.  490.  550.  575.  625.]
Class: <class 'MCMC.ResultContainer'>
bcs_params: {'hopping_constant': 0.25, 'uniform_phase': False, 'use_assaad': True, 'delta': 0.25, 'width': 32, 'seed': 123456789, 'J_constant': 0.02808988764044944, 'g_constant': 0.25, 'chem_potential': 0.0, 'temperature': 275.0}
mc_params: {'intervals': 5, 'seed': 234567, 'target_snapshots': 32, 'algorithm': 'cluster', 'observable_list': ['correlation_length', 'DOS']}
Code version: NA
Observable list: ['correlation_length', 'DOS']

















In [135]:



    


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



    











In [183]:



    


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



    


















    
























In [ ]:

Content source: dombrno/PG

Similar notebooks:

notebook.community | gallery | about