1

QuTiP を使用して、量子力学行列微分方程式 (リンドブラッド方程式) を解こうとしています。コードは次のとおりです。

from qutip import *
from matplotlib import *
import numpy as np

hamiltonian = np.array([[215, -104.1, 5.1, -4.3  ,4.7,-15.1 ,-7.8 ],
[-104.1,  220.0, 32.6 ,7.1, 5.4, 8.3, 0.8],
      [ 5.1, 32.6, 0.0, -46.8, 1.0 , -8.1, 5.1],
     [-4.3, 7.1, -46.8, 125.0, -70.7, -14.7, -61.5],
       [ 4.7, 5.4, 1.0, -70.7, 450.0, 89.7, -2.5],
    [-15.1, 8.3, -8.1, -14.7, 89.7, 330.0, 32.7],
     [-7.8, 0.8, 5.1, -61.5,  -2.5, 32.7,  280.0]])
H=Qobj(hamiltonian)
ground=Qobj(np.array([[ 0.0863685 ],
  [ 0.17141713],
  [-0.91780802],
  [-0.33999268],
  [-0.04835763],
  [-0.01859027],
  [-0.05006013]]))

rho0 = ground*ground.dag()
from scipy.constants import *
ktuple=physical_constants['Boltzmann constant in eV/K']
k = ktuple[0]* 8065.6
htuple = physical_constants['Planck constant in eV s'] 
hbar = (htuple[0]* 8065.6)/(2*pi)
gamma=(2*pi)*((k*300)/hbar)*(35/(150*hbar))

L1 = Qobj(np.array([[1,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L2 = Qobj(np.array([[0,0,0,0,0,0,0],[0,1,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L3 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,1,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L4 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,1,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L5 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,1,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]]))
L6 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,1,0],[0,0,0,0,0,0,0]]))
L7 = Qobj(np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,1]]))
#Since our gamma variable cannot be directly applied onto the Lindblad operator, we must multiply it with the collapse operators:

L1=gamma*L1
L2=gamma*L2
L3=gamma*L3
L4=gamma*L4
L5=gamma*L5
L6=gamma*L6
L7=gamma*L7

options = Options(nsteps=100000)

results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
print results

このコードは、次の方程式を解くことになっています。

リンドブラッド方程式

ここで、L_i は行列 (リスト: [L1,L2,L3,L4,L5,L6,L7])、H はハミルトニアン、別の行列、$\rho$密度行列、T は温度に$\ガンマ$等しい定数です。 $2\pi kT/\hbar*E_{R}/(\hbar\omega_{c})$、k はボルツマン定数、および$\hbar$ = $h/2\pi$、h はプランク定数です。

コードを実行するたびに、次のエラーが表示されます。

ZVODE--  At T (=R1) and step size H (=R2), the
       corrector convergence failed repeatedly
       or with abs(H) = HMIN
      In above,  R1 =  0.0000000000000D+00   R2 =  0.1202322246215D-36
/usr/local/lib/python2.7/dist-packages/scipy/integrate/_ode.py:853: UserWarning: zvode: Repeated convergence failures. (Perhaps bad Jacobian supplied or wrong choice of MF o
r tolerances.)
  'Unexpected istate=%s' % istate))
Traceback (most recent call last):
  File "lindbladqutip.py", line 48, in <module>
    results = mesolve(H, rho0, np.linspace(0.0, 1000, 20),[L1,L2,L3,L4,L5,L6,L7],[],options=options)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 264, in mesolve
    progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 692, in _mesolve_const
    return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
  File "/projects/d6138712-e5f4-4d85-9d4d-77ce0a7b4a61/.local/lib/python2.7/site-packages/qutip/mesolve.py", line 866, in _generic_ode_solve
    raise Exception("ODE integration error: Try to increase "
Exception: ODE integration error: Try to increase the allowed number of substeps by increasing the nsteps parameter in the Options class.

デバッグ分析を行った後、最初または 2 番目の統合が失敗したようです。このエラーは、私が試した nsteps パラメータを増やすように指示しています。それでも失敗します。時間のリストを変更しても (np.linspace 関数が時間のリストを作成します)、効果はありません。

このエラーを修正するために何ができるか知りたいです。詳細が必要な場合はコメントしてください。

ご助力いただきありがとうございます!

4

1 に答える 1

1

プログラミングの観点からは、問題はあなたの値にgammaあり、したがって崩壊演算子のサイズにあるようです。プリントアウトgamma- これは正しい順序です10**25- これがソルバーの収束を妨げているようです。

テストのために(私はエンジニアであり、量子物理学者ではありません...)、より小さな値gamma(たとえば0.1)を入力すると、ソルバーは機能しているようで、明らかに妥当な出力が得られますresults.states

私はあなたのことをよく理解していませんgamma-あなたが設定したように、単位はcm -1 s -2のようです。たぶん、1回だけ割りたいのだろうかhbar。私が言うように、私は量子物理学者ではないので、ここでは、プログラミングがうまくいかない理由と少しの次元分析を組み合わせたものに基づいて推測しているだけです。

編集

OPはコメントで、大きさ/単位の間違った順序がgammaプログラミングの問題(つまり、数値計算の収束の妨げ)であるように思われることを示していますが、ガンマの計算方法について完全には明確ではありません。この段階で、それについてhttp://physics.stackexchange.comまたはhttp://math.stackexchange.comに質問を投稿する価値があるかもしれません- 必要に応じて文脈のためにこれを参照してください。

編集2

物理サイトでこの関連する質問をしたことに注意してください。これにより、ガンマの式がどこから来たのかが明確になり、それによって単純に提示された定数項30150この質問で実際に単位があることが明確になります (それぞれエネルギーと周波数)。これにより、次元分析が変更されます。ガンマの単位は s -1または、適切に変換すると cm -1になります。

また、コメントで言及した値 - 300 cm -1も表示されます。

于 2015-06-03T01:00:23.283 に答える