#!/usr/bin/env python import os, time from scipy import * from pylab import * def noninteracting_G0(beta, N): om = arange(1,2*N,2)*pi/beta return (om, 2j*om*(1-sqrt(1+1/om**2)) ) # Bethe lattice, half-filling def InverseFourier(Giom, om, tau): "Inverse Fourier computes G(tau) from G(iom)" beta = pi/om[0] Gtau = zeros(len(tau), dtype=float) for it,t in enumerate(tau): dsum=0 for iw,w in enumerate(om): dsum += cos(w*t)*Giom[iw].real + sin(w*t)*(Giom[iw].imag+1./w) Gtau[it] = 2*dsum/beta - 0.5 return Gtau def g0_2D(tau, G0): L = len(tau)-1 print 'L=', L g0 = zeros((L,L), dtype=float) for i in range(L): for j in range(L): if i>=j: g0[i,j] = -G0[i-j] else: g0[i,j] = G0[L+i-j] return g0 def CleanUpdate(g, g0, vn): L = len(vn) A = zeros((L,L), dtype=float) for i,s in enumerate([1,-1]): # a = e^V-1 a = [exp(vn[l]*s)-1 for l in range(L)] print 'a=', a # Matrix A = 1 + (1-g)(e^V-1) == (1+a)*I - g*a for l1 in range(L): for l2 in range(L): A[l1,l2] = -g0[l1,l2]*a[l2] A[l1,l1] += 1+a[l1] print 'A=', A g[i] = solve(A,g0) if __name__=='__main__': beta = 16. # inverse temperature N = 1500 # number of matsubara points L = 64 # number of time slices U = 2. # interaction U ####### nwarm0 = 100 # how many sweeps to warm up nmeas0 = 10 # how often to meassure (in sweeps) nsteps = 100000 # number of all MC steps ncout = 10000 # how often to print info nwarmup = int(nwarm0*L) measure = int(nmeas0*L) (om, G0start) = noninteracting_G0(beta, N) seed(1) S_ising = sign(rand(L)-0.5) G0iom = G0start tau = linspace(0,beta,L+1) G0tau = InverseFourier(G0iom, om, tau) xlam = arccosh(exp(0.5*(tau[1]-tau[0])*U)) # lambda vn = array(S_ising)*xlam # g0 in 2D representation g0 = g0_2D(tau, G0tau) g = zeros((2,L,L), dtype=float) CleanUpdate(g, g0, vn) print 'g=' print g #plot(tau[:-1], g0[:,0], label='1') #plot(tau[:-1], g0[0,:], label='2') #legend(loc='best') #show()