qutip是个很有用的量子计算集成开发包，abbreviated from Quantum toolbox in Python.

但是确实使用时候遇到了很多困难和坑，正好目前所在的量子计算组和大创项目需要用到，就记录一下遇到的问题。

  1. 绘制Bloch Sphere用add_states()方法
1
2
3
sphere=qutip.Bloch()
sphere.add_states(result.states)
sphere.show()

报错

1
object has no attribute 'do_3d_projection'

原因是我用的matplotlib版本是3.5.1，的确没有了这个属性，需要降到3.3.2

  1. 在计算基态某个物理量的平均值时，用Qobj.groundstates()方法
1
print(qt.expect(qt.simgaz(),H.groundstate()))

报错

1
TypeError: Arguments must be quantum objects or eseries

改正方法，H.groundstate()返回本征值和本征矢

1
print(qt.expect(qt.simgaz(),H.groundstate()[1]))

尝试过将H.groundstate()强制转化成Qobj，也尝试过将qt.simgaz(),H.groundstate()转化成array用numpy恩算。

  1. Ising模型$\sigma^z_i \sigma^z_{i+1} $ for循环叠加
1
2
3
4
5
6
7
8
9
10
11
N=10

def sigma_x(index):
    sigma_x=qt.tensor(qt.qeye(2**index),qt.tensor(qt.tensor(qt.sigmax(),qt.qeye(2)), qt.qeye(2**(N-index-2))))
   
def hamiltionian_groundstate(g):
    H=0
    for i in range(N-1):
        H=H+g*sigma_x(i)-sigma_z(i)
    return H.groundstate()[1]

报错，thx to QuTiP不完全ガイド，日本人写的都比中文互联网全。

1
Incompatible quantum object dimensions

需要让每个sigma_x(i)都是相同dims的Qobj，解决方案

1
2
3
def sigma_x(index):
    sigma_x=qt.tensor(qt.qeye(2**index),qt.tensor(qt.tensor(qt.sigmax(),qt.qeye(2)), qt.qeye(2**(N-index-2))))
    return qt.Qobj(sigma_x, dims=[[2**N,2**N],[2**N,2**N]])
  1. 将数据保存成.csv and .xls等文件时若单元格空白则再用pandas读取时该数据为np.nan，应该用df.replace(np.nan,0),报错：
1
FutureWarning: elementwise comparison failed
  1. if 循环时，判断单元格是否为文本，否则加入列表中，结果最后列表都是空列表。noting that the judgement of type do not need the like isdigit() func isalpha(), just type ==type.
1
2
3
for i in range(len(df.columns)):
    if type(df[i][1])==type('a'):
        break

肯定不能用pass对吧，否则这条就没啥用了。用break的话只执行当前的i，后续的i都不执行力。那应该：

1
2
3
for i in range(len(df.columns)):
    if type(df[i][1])==type('a'):
        continue

这样跳过此i，执行下一个i。By the way it is worthing to pointing out that if the data is digit of string type, then transforming it to number type:

1
2
3
4
5
a='123'
a=int(''.join(a))
#or rather:
a=int(a)

join方法是字符串的操作函数，str.join(iterable)，对一个可迭代对象执行用str分割，最后返回一个str。这个函数很好玩，可以玩出来花来，比如把列表每个元素分行打印：

1
2
lst=['a','a','a']
'\n'.join(list)
1
2
3
4
print(*lst,sep='\n')

for i in lst:
    print(i,sep='\n')

matplotlib配色

  1. 高精度计算问题
    在计算

CB=(∂E∂T)B=−Nμ2Bˉ2ln⁡2kBT2csch(μBˉkBT)C_{B}=\left( \frac{\partial E}{\partial T}\right) _{B}=-\frac{N\mu^{2}\bar{B}^{2}\ln{2}}{k_{B}T^{2}}\text{csch}\left( \frac{\mu\bar{B}}{k_{B}T} \right)

时，因带着μ,kB\mu,k_{B}，其数量级为-e12，e-23,而cschx=2np.exp(x)/(np.exp(2x)-1) # 2/(np.exp(x)-np.exp(-x))上溢，为零。

1
2
RuntimeWarning: overflow encountered in exp
print(np.exp(1e30))

所以就成了计算exp(-e30)问题，python精度最大为2^52，直接溢出，lg(exp(1e30))大概等于e29这么多位，，，，能算出来就有鬼了，解决办法：

1
2
3
4
5
6
7
8
1. 每一步计算都转换科学记数法分别算前面和10的次数部分，只计算科学计数法前半部分，后半部分就带着

2. 字符串处理，和科学计数法差不多，有点类似于c语言字符串处理高精,double11位指数位可以支持2^-2047,调gmp或者mpfr的计算后端

3. c=ab d=lgc=lga+lgb c=exp(d) 
    lg能降低次数，还保证单调递增,把要存的实数转化为浮点数最精确的幅度内进行计算,别有特别大和特别小的数相加减,否则损失精度，如果用lg的话也可能会遇到比如这个数字有1e100那么多位，这也爆了

4. 转化单位制，就像自然单位制一样，c=hbar=G=1
  1. 量子多体数值就别tm用qutip了
    因为之前做的量子优化控制偏量子光学，所以现在来做量子多体习惯性的继续用qutip，但是遇到到大失败！以下是用qutip and scipy.sparse来做相关计算的用时记录。对于橫场伊辛模型来说，发现qutip生成多体Hamiltionian的速度极快，比scipy.sparse还要快，但是求解矩阵本征态有着多项式级别的差异，同时矩阵点乘的速度也相差不多，而且qutip还有内存的问题。qutip老老实实算量子光学得了。

L=8

1
2
3
4
5
6
t0=datetime.datetime.now()
for i in range(6,22):
    E,states=sl.eigs(test_H(),k=i，which='SA') #the parameter(often translated as formal parameter), while argument often refers to actual parameter. 
    t2=datetime.datetime.now()
    print(E)
    print(t2-t0)

求解一个本征矢需要多花0:00:00.05

1
2
3
E,states=sl.eigs(test_H(),k=2**8-2,which='SR')
E1,states1=sl.eigsh(test_H(),k=2**8-2,which='SA') #最大只能到N-1，不太理解为什么不能生成N个本征矢。t1=0:00:00.278994
t2=0:00:00.199048

所以可能还是eigsh faster than eigs.

L=12 10times
qutip H.eigenstats() 0:00:59.843034
sl.eigsh 0:00:01.549344
sl.eigs 0:00:01.491331
E1,states1=scipy.linalg.eig(test_H().toarray()) 0:00:45.689474 #one time

L=14
eigsh 0:00:00.806247
eigs 0:00:00.908200

L=16 10 times
eigsh 0:00:27.331040
eigs 0:00:24.588102
此时qutip会爆内存，因为32GB的电脑最大存储16个格点的系统, 2^162^162^3=32GB，一个float类型变量占8B.

dot
np.dot(H,v[:,i]) 0:00:00.090018 #取本征矢时用 v[:,i]
H*states 0:00:00.096038

conclusion: Viva la scipy, go bitch qutip

In fact, the eigsh func is a packaging of Arpack package of C and Fortran, using the lanczos algorithm.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def lanczos(hm,length):
    dimension = 2**length
    diagElements = []
    subDiagElements = []
    
    f0 = np.random.rand(dimension)
    f0 = f0 / norm(f0)
    
    f1 = hm @ f0
    lastEnergy = 0
    groundEnergy = 0
    for truncationIndex in range(0,dimension):
        hdiag = f0 @ f1
        f1 = f1 - hdiag * f0
        
        hsub = norm(f1)
        f1 = f1 / hsub 
        
        [f0 , f1] = [f1 , hm @ f1 -hsub*f0]
        diagElements.append(hdiag)

        # 常规方法对角化
        # [va,ve] = scipy.linalg.eigh(self.fillElements(diagElements,subDiagElements,len(diagElements)))
        # groundEnergy = va[0]

        # if abs(groundEnergy - lastEnergy) < 1e-10:
        #     break
        
        # 三对角矩阵方法
        [va,ve] = eigh_tridiagonal(np.array(diagElements),np.array(subDiagElements))
        groundEnergy = va[0]
        if(abs(groundEnergy - lastEnergy) < 1e-10):
            break
        
        lastEnergy = groundEnergy
        subDiagElements.append(hsub)

    return groundEnergy

Reference

User Guide for QuTip