python输入一个复数_python写的一个绝对对角化的一个程序跑出来能量是复数

引用回帖:

KalaShayminS at 2015-03-17 00:21:45

这没有程序想看也看不了啊

def dim(L, Np, Nd):

# number of basis states,Np denotes the number of spin up fermions and

# Nd denote the spin-down fermions

return ((math.factorial(L) * math.factorial(L)) /

(math.factorial(Nd) * math.factorial(Np) *

math.factorial(L - Np) * math.factorial(L - Nd)))

def dim_compact(L, N):

return (math.factorial(L)) / (math.factorial(N) * math.factorial(L - N))

def matglue(a, b):

r = []

for i in a:

for j in b:

r.append(concatenate([i, j]))

u = np.matrix(r)

return u

def NinL(L, N):

# Ns=dim(L,N)

Ns = dim_compact(L, N)

# array that stores all basis states in terms of lattice occupation

# numbers (list of occupation number states):

nl = zeros((Ns, L), int)

# first basis state with all particles ordered to the left:

for l in range(N):

nl[0, l] = 1

for a in range(Ns-1):

if nl[a, L-1] != 1:

for j in reversed(range(L-1)):

if nl[a, j] == 1 and nl[a, j + 1] == 0:

for i in range(j):

nl[a + 1, i] = nl[a, i]

nl[a+1, j] = 0

nl[a+1, j+1] = 1

break

else:

for j in reversed(range(L-1)):

if nl[a, j] == 1 and nl[a, j+1] == 0:

for i in range(j):

nl[a+1, i] = nl[a, i]

nl[a+1, j] = 0

nl[a+1, j+1] = 1

c = 0

for i in range(j+2, L):

c += nl[a, i]

for i in range(j+2, j+2+c):

nl[a+1, i] = 1

break

return nl

def basis(L, Nd, Np):

Mp = NinL(L, Np)

Md = NinL(L, Nd)

nl = matglue(Mp, Md)

return nl

# single-particle Hamiltonian definition:

def H1(V, b, phi, j):

return V * cos(2 * pi * b * j + phi)

# doublon number of a configuration comprised of spinless-half fermions

def doub(nl,L):

conf = 0

nl = nl.tolist()[0]

# check wether a site is occupied by a doublon

for i in range(L):

if nl== 1 and nl[i+L] == 1:

conf += 1

return conf

# adding hopping element in the Hamiltonian

def new_hopmat(L, Nd, Np):

nl = basis(L, Nd, Np)

Ns = dim(L, Np, Nd)

H = zeros((Ns, Ns))

L2 = 2 * L

#N = 0

lookup = {}

for a1 in range(Ns):

m = nl[a1]

l = m.tolist()[0]

#print l

k = ''.join(str(i) for i in l)

lookup.setdefault(k,  a1)

#print lookup

for a1 in range(Ns):

for i in range(L2 - 1):

if nl[a1, i] != nl[a1, i + 1]:

tmp = [] + nl[a1].tolist()[0]

tmp= 1 - nl[a1, i]

tmp[i + 1] = 1 - nl[a1, i + 1]

k = ''.join(str(j) for j in tmp)

if lookup.has_key(k):

a2 = lookup[k]

if nl[a1, i] == 1:

H[a1, a2] += -1.

else:

H[a1,a2] += 1.

#N += -1.

#print N

#H[a2, a1] += -1.

return H, Ns, nl

# many-body Hamiltonian definition:

def Hmany(L, Nd, Np, V, b, phi, U):

#H, Ns, nl = Hopmat(L, Nd, Np)

H, Ns, nl = new_hopmat(L, Nd, Np)

for a in range(Ns):

H[a, a] += U * doub(nl[a],L)

#print(H[a,a], a)

# sum over all particle and spin of single particle potential

#for j in range(2 * L):

#H[a, a] += H1(V, b, phi, j) * nl[a, j]

return H

，