Added tutorial0

tutorials/README.md

@@ -1,6 +1,23 @@
 # List of tutorials
+- [Tutorial0](#tutorial0)
 - [Tutorial1](#tutorial1)
 
+## Tutorial0
+This is an example of an single-band Anderson impurity model.
+$python mk_model.py                        #generate Coulomb matrix & hybridization function & generate hopping matrix
+$mpirun -np 24 /path/to/hybmat input.ini   #run the solver with 24 MPI processes
+$python plot.py                            #plot Green's function. The plot will be written into GF.eps.
+
+At the end of the simulation, the results are written into the HDF5 file "input.out.h5".
+A HDF5 file can be read easily e.g. by using the h5py library in Python.
+In the script file "plot.py", one reads out the Green's function from the output file,
+plot the data, and save the plot as a eps file.
+Note that latex must be required to run plot.py.
+If latex is not available on your cluster, please just download input.out.h5 to your local system and run plot.py.
+The graph looks like this (simulation time was 5min with 24 MPI processes).
+![](tutorial0/GF.png)
+
+
 ## Tutorial1
 **We strongly recommend to run this tutorial at least with 24 MPI processes and 600 sec. Otherwise, the calculation will fail to converge.**
 We explain how to solve a three-orbital model with spin-orbit coupling and Slater-Kanamori interaction.

tutorials/tutorial0/gen_Uijkl.py

@@ -0,0 +1,37 @@
+import numpy as np
+import sys
+
+up = 0
+down = 1
+
+# H_int = (1/2) * \sum_{ijkl} U_{ijkl} c_i^\dagegr c_j^\dagger c_k c_l
+def generate_U_tensor_onsite(n_orb, U):
+    U_tensor = np.zeros((n_orb,2,n_orb,2,n_orb,2,n_orb,2),dtype=complex)
+
+    for iorb in xrange(n_orb):
+        U_tensor[iorb, up, iorb, down, iorb, down, iorb, up] = U
+        U_tensor[iorb, down, iorb, up, iorb, up, iorb, down] = U
+
+    return U_tensor, 2*n_orb
+
+n_site = 2
+Uval = 4.0
+
+V_mat = np.identity(2*n_site,dtype=complex)
+
+U_tensor, num_elem = generate_U_tensor_onsite(n_site, Uval)
+
+f = open("Uijkl.txt", "w")
+print >>f, num_elem
+line = 0
+for iorb1 in xrange(n_site):
+    for iorb2 in xrange(n_site):
+        for iorb3 in xrange(n_site):
+            for iorb4 in xrange(n_site):
+                for isp in xrange(2):
+                    for isp2 in xrange(2):
+                        if U_tensor[iorb1,isp,iorb2,isp2,iorb3,isp2,iorb4,isp] != 0.0:
+                            print >>f, line, "   ", 2*iorb1+isp, 2*iorb2+isp2, 2*iorb3+isp2, 2*iorb4+isp, U_tensor[iorb1,isp,iorb2,isp2,iorb3,isp2,iorb4,isp].real, U_tensor[iorb1,isp,iorb2,isp2,iorb3,isp2,iorb4,isp].imag
+                            line += 1
+
+f.close()

tutorials/tutorial0/gen_hopping.py

@@ -0,0 +1,14 @@
+import sys
+import numpy as np
+from scipy.integrate import simps
+from math import pi
+
+nf=4
+norb=nf/2
+
+f = open('hopping.txt','w')
+for iorb in xrange(nf):
+    for jorb in xrange(nf):
+        print>>f, iorb, jorb, 0.0, 0.0
+f.close()
+

tutorials/tutorial0/gen_hyb.py

@@ -0,0 +1,54 @@
+import numpy as np
+from scipy.integrate import simps
+from math import pi
+
+def ft_to_tau_hyb(ndiv_tau, beta, matsubara_freq, tau, Vek, data_n, data_tau, cutoff):
+ for it in range(ndiv_tau+1):
+     tau_tmp=tau[it]
+     if it==0:
+         tau_tmp=1E-4*(beta/ndiv_tau)
+     if it==ndiv_tau:
+         tau_tmp=beta-1E-4*(beta/ndiv_tau)
+     ztmp=0.0
+     for im in range(cutoff):
+         ztmp+=(data_n[im]+1J*Vek/matsubara_freq[im])*np.exp(-1J*matsubara_freq[im]*tau_tmp)
+     ztmp=ztmp/beta
+     data_tau[it]=2.0*ztmp.real-0.5*Vek
+
+
+vbeta=20.0
+ndiv_tau=1000
+nf=4
+
+matsubara_freq=np.zeros((ndiv_tau,),dtype=float)
+for im in range(ndiv_tau):
+    matsubara_freq[im]=((2*im+1)*np.pi)/vbeta
+
+tau=np.zeros((ndiv_tau+1,),dtype=float)
+for it in range(ndiv_tau+1):
+    tau[it]=(vbeta/ndiv_tau)*it
+
+ndiv_dos = 10000
+W=2.0
+e_smp=np.linspace(-W,W,ndiv_dos)
+dos_smp=np.sqrt(W**2-e_smp**2)/(0.5*pi*W**2)
+ek_var = simps(dos_smp*(e_smp**2), e_smp)
+
+#Bath Green's function (g_omega, g_tau)
+g_omega=np.zeros((ndiv_tau,),dtype=complex)
+g_tau=np.zeros((ndiv_tau+1,),dtype=complex)
+for im in range(ndiv_tau):
+    f_tmp = dos_smp/(1J*matsubara_freq[im]-e_smp)
+    g_omega[im]=simps(f_tmp,e_smp)
+
+ft_to_tau_hyb(ndiv_tau,vbeta,matsubara_freq,tau,1.0,g_omega,g_tau,ndiv_tau)
+
+f = open('delta.txt', 'w')
+for i in range(ndiv_tau+1):
+    for j in range(nf):
+        for k in range(nf):
+            if j==k:
+                print >>f, i, j, k, g_tau[i].real, g_tau[i].imag
+            else:
+                print >>f, i, j, k, 0.0, 0.0
+f.close()

tutorials/tutorial0/input.ini

@@ -0,0 +1,22 @@
+seed=20
+timelimit=300
+
+[model]
+sites=1
+spins=2
+# Or, you can use sites=2, spins=1.
+# Only the difference is that a global spin flip is not attempted for this alternative choice.
+coulomb_tensor_input_file="Uijkl.txt"
+hopping_matrix_input_file="hopping.txt"
+beta=20.0
+
+#Delta(tau) must involve data at the end point \tau=\beta.
+delta_input_file="delta.txt"
+
+#The number of tau points in Delta(tau) - 1
+n_tau_hyb=1000
+
+[measurement.G1]
+n_legendre=50
+n_tau=1000
+n_matsubara=500

tutorials/tutorial0/mk_model.py

@@ -0,0 +1,102 @@
+import numpy as np
+import sys
+from math import pi
+from scipy.integrate import simps
+import random
+
+up = 0
+down = 1
+
+def generate_U_tensor_onsite(n_orb, U):
+    U_tensor = np.zeros((n_orb,2,n_orb,2,n_orb,2,n_orb,2),dtype=complex)
+
+    for iorb in xrange(n_orb):
+        U_tensor[iorb, up, iorb, down, iorb, down, iorb, up] = U
+        U_tensor[iorb, down, iorb, up, iorb, up, iorb, down] = U
+
+    return U_tensor, 2*n_orb
+
+def ft_to_tau_hyb(ntau, beta, matsubara_freq, tau, Vek, data_n, data_tau, cutoff):
+ for it in range(ntau+1):
+     tau_tmp=tau[it]
+     if it==0:
+         tau_tmp=1E-4*(beta/ntau)
+     if it==ntau:
+         tau_tmp=beta-1E-4*(beta/ntau)
+     ztmp=0.0
+     for im in range(cutoff):
+         ztmp+=(data_n[im]+1J*Vek/matsubara_freq[im])*np.exp(-1J*matsubara_freq[im]*tau_tmp)
+     ztmp=ztmp/beta
+     data_tau[it]=2.0*ztmp.real-0.5*Vek
+
+
+#Parameters
+n_site = 1
+nf = n_site*2
+
+Uval = 4.0
+vbeta = 20.0
+ntau = 1000
+mu = 0.5*Uval
+
+#Generate Coulomb tensor
+U_tensor, num_elem = generate_U_tensor_onsite(n_site, Uval)
+f = open("Uijkl.txt", "w")
+print >>f, num_elem
+line = 0
+for iorb1 in xrange(n_site):
+    for iorb2 in xrange(n_site):
+        for iorb3 in xrange(n_site):
+            for iorb4 in xrange(n_site):
+                for isp in xrange(2):
+                    for isp2 in xrange(2):
+                        if U_tensor[iorb1,isp,iorb2,isp2,iorb3,isp2,iorb4,isp] != 0.0:
+                            print >>f, line, "   ", 2*iorb1+isp, 2*iorb2+isp2, 2*iorb3+isp2, 2*iorb4+isp, U_tensor[iorb1,isp,iorb2,isp2,iorb3,isp2,iorb4,isp].real, U_tensor[iorb1,isp,iorb2,isp2,iorb3,isp2,iorb4,isp].imag
+                            line += 1
+
+f.close()
+
+#Generate hopping matrix
+f = open('hopping.txt','w')
+for iorb in xrange(nf):
+    for jorb in xrange(nf):
+        if iorb == jorb:
+            print>>f, iorb, jorb, -mu, 0.0
+        else:
+            print>>f, iorb, jorb, 0.0, 0.0
+f.close()
+
+
+#Generate hybridization function
+matsubara_freq=np.zeros((ntau,),dtype=float)
+for im in range(ntau):
+    matsubara_freq[im]=((2*im+1)*np.pi)/vbeta
+
+tau=np.zeros((ntau+1,),dtype=float)
+for it in range(ntau+1):
+    tau[it]=(vbeta/ntau)*it
+
+ndiv_dos = 10000
+W=2.0
+e_smp=np.linspace(-W,W,ndiv_dos)
+dos_smp=np.sqrt(W**2-e_smp**2)/(0.5*pi*W**2)
+ek_var = simps(dos_smp*(e_smp**2), e_smp)
+
+#Bath Green's function (g_omega, g_tau)
+g_omega=np.zeros((ntau,),dtype=complex)
+g_tau=np.zeros((ntau+1,),dtype=complex)
+for im in range(ntau):
+    f_tmp = dos_smp/(1J*matsubara_freq[im]-e_smp)
+    g_omega[im]=simps(f_tmp,e_smp)
+
+ft_to_tau_hyb(ntau,vbeta,matsubara_freq,tau,1.0,g_omega,g_tau,ntau)
+
+f = open('delta.txt', 'w')
+for i in range(ntau+1):
+    for j in range(nf):
+        for k in range(nf):
+            if j==k:
+                print >>f, i, j, k, g_tau[i].real, g_tau[i].imag
+            else:
+                print >>f, i, j, k, 0.0, 0.0
+f.close()

tutorials/tutorial0/plot.py

@@ -0,0 +1,114 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import pylab
+import h5py
+
+def read_param(h5, name):
+    if '/parameters/dictionary/'+name in h5:
+        return h5['/parameters/dictionary/'+name].value
+    elif '/parameters/'+name in h5:
+        return h5['/parameters/'+name].value
+    else:
+        raise RuntimeError("Parameter "+ name + " not found") 
+
+def compute_Tnl(n_matsubara, n_legendre):
+    Tnl = np.zeros((n_matsubara, n_legendre), dtype=complex)
+    for n in xrange(n_matsubara):
+        sph_jn = scipy.special.sph_jn(n_legendre, (n+0.5)*np.pi)[0]
+        for il in xrange(n_legendre):
+            Tnl[n,il] = ((-1)**n) * ((1J)**(il+1)) * np.sqrt(2*il + 1.0) * sph_jn[il]
+    return Tnl
+
+def read_h5(p):
+    r = {}
+
+    print p+'.out.h5'
+    h5 = h5py.File(p+'.out.h5','r')
+
+    r["SITES"] = read_param(h5, 'model.sites')
+    r["BETA"] = read_param(h5, 'model.beta')
+
+    def load_g(path):
+        N = h5[path].shape[0]
+        M = h5[path].shape[1]
+        data = h5[path].value.reshape(N,M,M,2)
+        return data[:,:,:,0] + 1J*data[:,:,:,1]
+
+    r["Gtau"] = load_g('/gtau/data')
+
+    r["Gomega"] = load_g('/gf/data')
+
+    r["Sign"] = h5['/simulation/results/Sign/mean/value'].value
+
+    r["Sign_count"] = h5['/simulation/results/Sign/count'].value
+
+    return r
+
+# input: G2 in the mix basis
+# return G2 in the Matsubara freq. domain: (i,j,k,l, fermionic freq, fermionic freq, bosnic freq)
+def compute_G2_matsubara(g2_l, niw_f):
+    nl_G2 = g2_l.shape[4]
+    Tnl_G2 = compute_Tnl(niw_f, nl_G2)
+    tmp = np.tensordot(g2_l, Tnl_G2.conjugate(), axes=(5,1))
+    return np.tensordot(Tnl_G2, tmp, axes=(1,4)).transpose((1,2,3,4,0,6,5))
+
+prefix_list = ['input']
+result_list = []
+for p in prefix_list:
+    result_list.append(read_h5(p))
+
+color_list = ['r', 'g', 'b', 'y', 'k', 'm']
+params = {
+    'backend': 'ps',
+    'axes.labelsize': 24,
+    'text.fontsize': 24,
+    'legend.fontsize': 18,
+    'xtick.labelsize': 24,
+    'ytick.labelsize': 24,
+    'text.usetex': True,
+    }
+pylab.rcParams.update(params)
+plt.figure(1,figsize=(8,8))
+plt.subplot(211)
+plt.xlabel(r'$\tau/\beta$', fontname='serif')
+plt.ylabel(r'$-\mathrm{Re}G(\tau)$', fontname='serif')
+plt.yscale('log')
+
+plt.subplot(212)
+plt.xlabel(r'$\omega_n$', fontname='serif')
+plt.ylabel(r'$-\mathrm{Im}G(i\omega_n)$', fontname='serif')
+plt.xscale('log')
+plt.yscale('log')
+
+for i in range(len(result_list)):
+    norb = result_list[i]["SITES"]
+    beta = result_list[i]["BETA"]
+    nf = norb*2
+
+    sign = result_list[i]["Sign"]
+    gf_legendre = result_list[i]["Gtau"]
+    gomega_l = result_list[i]["Gomega"]
+
+    print "The number of measurements is ", result_list[i]["Sign_count"]
+
+    tau_point = np.linspace(0.0, 1.0, gf_legendre.shape[0])
+    plt.subplot(211)
+    for i_f in range(nf):
+        plt.plot(tau_point, -gf_legendre[:,i_f,i_f].real, color=color_list[i_f], marker='', label='flavor'+str(i_f), ls='--', markersize=0)
+
+    omega_point = np.array([(2*im+1)*np.pi/beta for im in xrange(gomega_l.shape[0])])
+    plt.subplot(212)
+    for i_f in range(nf):
+        plt.plot(omega_point, -gomega_l[:,i_f,i_f].imag, color=color_list[i_f], marker='', label='flavor'+str(i_f), ls='--', markersize=0)
+    plt.plot(omega_point, 1/omega_point, color='k', label=r'$1/\omega_n$', ls='-')
+
+
+plt.subplot(211)
+plt.legend(loc='best',shadow=True,frameon=False,prop={'size' : 12})
+
+plt.subplot(212)
+plt.legend(loc='best',shadow=True,frameon=False,prop={'size' : 12})
+
+plt.tight_layout()
+plt.savefig("GF.eps")
+plt.close(1)