| 查看: 726 | 回复: 3 | ||
[求助]
python写的一个绝对对角化的一个程序跑出来能量是复数 已有1人参与
|
用python写了一个一维相互作用0 disorder的自旋为1/2的费米子系统，basis用Fock空间来表示，可是算出来的能量总是为复数，而且，博后要我画transition energy和相互作用U的曲线，在L个晶格，总数为L个电子的时候，不同晶格数的曲线是相互交叉的，不管怎么样我都得不到这样的答案，这个问题已经困扰我一个月了。大神们能不能帮我看一下 |
回复此楼» 猜你喜欢
284求调剂
已经有10人回复
-
一志愿山东大学药学学硕求调剂
已经有4人回复
-
07化学280分求调剂
已经有4人回复
-
298-一志愿中国农业大学-求调剂
已经有12人回复
-
求材料，环境专业调剂
已经有3人回复
-
335求调剂
已经有5人回复
-
求调剂
已经有7人回复
-
一志愿吉大化学322求调剂
已经有4人回复
-
环境学硕288求调剂
已经有8人回复
-
341求调剂(一志愿湖南大学070300)
已经有6人回复
2楼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 |
3楼2015-03-17 03:26:48
4楼2015-03-17 17:02:53
