Difference between revisions of "ALPS 2 Tutorials:TEBD-02 kink/ja"

From ALPS
Jump to: navigation, search
(Created page with "{{Languages|ALPS_2_Tutorials:TEBD-02_kink}} <!--= Evolution of a domain Wall =--> = 磁壁の展開 = <!--In this tutorial we will study the time evolution of a S=1/2 spin ch…")
 
 
Line 1: Line 1:
 
{{Languages|ALPS_2_Tutorials:TEBD-02_kink}}
 
{{Languages|ALPS_2_Tutorials:TEBD-02_kink}}
  
<!--=  Evolution of a domain Wall =-->
 
 
=  磁壁の展開 =
 
=  磁壁の展開 =
  
<!--In this tutorial we will study the time evolution of a S=1/2 spin chain prepared in a nonequilibrium state.  The particular state that we choose is that with all spins to the left of the chain center "down" and all of those to the right of the center "up," <math>| \downarrow \downarrow \dots \downarrow \uparrow \dots \uparrow \uparrow\rangle</math>.  This state can be chosen as the initial state by setting INITIAL_STATE to be 'kink'.  Some exact results are known regarding the evolution of this state under the 1D XX model, which allows for a detailed study of the errors present in TEBD.-->
 
 
このチュートリアルでは、非平衡状態で準備したS=1/2スピン鎖の時間発展について学習します。選択した特殊な状態は、次のように、真ん中で上下スピンをもち、左右でそれぞれ上スピン、下スピンを持つような構造です。<br>
 
このチュートリアルでは、非平衡状態で準備したS=1/2スピン鎖の時間発展について学習します。選択した特殊な状態は、次のように、真ん中で上下スピンをもち、左右でそれぞれ上スピン、下スピンを持つような構造です。<br>
 
<math>| \downarrow \downarrow \dots \downarrow \uparrow \dots \uparrow \uparrow\rangle</math><br>
 
<math>| \downarrow \downarrow \dots \downarrow \uparrow \dots \uparrow \uparrow\rangle</math><br>
Line 11: Line 9:
 
== XX modelの完全解 ==
 
== XX modelの完全解 ==
  
<!--The time evolution of the kink initial state under the XX model was solved exactly in [http://link.aps.org/doi/10.1103/PhysRevE.59.4912 Phys. Rev. E 59, 4912 (1999)] by a Jordan-Wigner transformation to free fermions.  It was found that the expectation value of the magnetization at any site as a function of time can be represented as a sum of Bessel functions, and the magnetization in the limit of long times and large distances from the initial domain wall approaches a scaling form in the variable <math>n/t </math>, where <math>n</math> is the distance from the center and <math>t</math> the time.  Explicitly, we have-->
+
XXモデルでのkink初期状態の時間発展はJordan-Wigner[http://link.aps.org/doi/10.1103/PhysRevE.59.4912 Phys. Rev. E 59, 4912 (1999)]の自由フェルミオンへの相転移によって厳密解が解かれています。時間関数として、どんなサイトでも磁化の期待値はBessel関数の合計として表現できることが明らかとなりました。そして、磁化と初期磁壁から大きな距離は、変数<math>n/t </math>でスケーリングします。<math>n</math>は中心からの距離、<math>t</math>は時間です。
XXモデルでのkink初期状態の時間発展はJordan-Wigner[http://link.aps.org/doi/10.1103/PhysRevE.59.4912 Phys. Rev. E 59, 4912 (1999)]の自由フェルミオンへの相転移によって厳密解が解かれています。
 
時間関数として、どんなサイトでも磁化の期待値はBessel関数の合計として表現できることが明らかとなりました。そして、磁化と初期磁壁から大きな距離は、変数<math>n/t </math>でスケーリングします。<math>n</math>は中心からの距離、<math>t</math>は時間です。
 
  
 
<math> M(n,t)=-\frac{1}{2}\sum_{i=1-n}^{n-1}j_i^2\left(t\right)</math>
 
<math> M(n,t)=-\frac{1}{2}\sum_{i=1-n}^{n-1}j_i^2\left(t\right)</math>
Line 19: Line 15:
 
<math>\lim_{n\to \infty} \lim_{t\to \infty} M(n,t)\to \phi\left(n/t\right)=-\frac{1}{\pi}\arcsin\left(n/t\right)</math>
 
<math>\lim_{n\to \infty} \lim_{t\to \infty} M(n,t)\to \phi\left(n/t\right)=-\frac{1}{\pi}\arcsin\left(n/t\right)</math>
  
<!--where <math> M(n,t)</math> is the magnetization a distance <math>n</math> from the center and <math> j_i(t)</math> is the Bessel function of order <math>i </math>.  In the first part of this tutorial we demonstrate these two results.-->
 
 
ここで、<math> M(n,t)</math>は磁化で、<math>n</math>は中心からの距離、<math> j_i(t)</math>はオーダ<math>i </math>のBessel関数を表します。このチュートリアルの最初のパートでは、これらの2つの結果を演習します。
 
ここで、<math> M(n,t)</math>は磁化で、<math>n</math>は中心からの距離、<math> j_i(t)</math>はオーダ<math>i </math>のBessel関数を表します。このチュートリアルの最初のパートでは、これらの2つの結果を演習します。
  
Line 54: Line 49:
 
         }]
 
         }]
  
<!--The math and scipy.special modules are required to generate the special functions needed to compare with the exact solution.  Note that we have chosen POWS to be zero, which corresponds to no quenching at all.  Thus, the values of GS, GIS, and GFS are arbitrary, and TAUS and NUMSTEPS give us the total simulation time and the number of time steps, respectively.  We write the input files, run the simulation, and get the output as usual:-->
 
 
mathとscipy.specialモジュールは、この厳密解を比較するために必要な特別な関数を生成するために使用されます。POWSを0としていることに注意してください。まったくクエンチしないことに対応しています。GS、GIS、GFSの値は任意です。TAUS、NUMSTEPSはそれぞれ、全計算時間とタイムステップを与えます。入力ファイルを準備し、計算を実行し、出力を得て、最後にプロットをおこないます。
 
mathとscipy.specialモジュールは、この厳密解を比較するために必要な特別な関数を生成するために使用されます。POWSを0としていることに注意してください。まったくクエンチしないことに対応しています。GS、GIS、GFSの値は任意です。TAUS、NUMSTEPSはそれぞれ、全計算時間とタイムステップを与えます。入力ファイルを準備し、計算を実行し、出力を得て、最後にプロットをおこないます。
  
Line 64: Line 58:
 
  Data=pyalps.load.loadTimeEvolution(pyalps.getResultFiles(prefix='tutorial_2a'), measurements=['Local Magnetization'])
 
  Data=pyalps.load.loadTimeEvolution(pyalps.getResultFiles(prefix='tutorial_2a'), measurements=['Local Magnetization'])
  
<!--We now must postprocess the raw output to compare with the exact solution.  To do this we first define empty arrays to hold the postprocessed data-->
 
 
厳密解と比較するために、計算結果の後処理が必要です。 まず最初に、後処理データを保持するための空の配列を定義します。
 
厳密解と比較するために、計算結果の後処理が必要です。 まず最初に、後処理データを保持するための空の配列を定義します。
  
Line 74: Line 67:
 
  scalingForm=[]
 
  scalingForm=[]
  
<!--we then calculate the exact result from the time data, and use the computed values of the magnetization at each site to compare with the exact solution.-->
 
 
時間データから厳密解を計算します。比較するために各サイトの磁化の計算された値を使用します。
 
時間データから厳密解を計算します。比較するために各サイトの磁化の計算された値を使用します。
  
Line 107: Line 99:
  
  
<!--Next, we calculate the exact scaling function, and then compute magnetization as a function of the scaling variable <math> n/t</math> to compare with the exact solution-->
 
 
次に、正確なスケーリング関数を計算します。スケーリング変数の関数として磁化を計算し、<math> n/t</math>厳密解と比較します。
 
次に、正確なスケーリング関数を計算します。スケーリング変数の関数として磁化を計算し、<math> n/t</math>厳密解と比較します。
  
  #compute the scaling form
+
#compute the scaling form
 
  # \phi(n/t)=-(1/pi)*arcsin(n/t) that M(n,t) approaches as n->infinity and t->infinity
 
  # \phi(n/t)=-(1/pi)*arcsin(n/t) that M(n,t) approaches as n->infinity and t->infinity
 
  # and compare it with the numerically computed values of M(n/t)
 
  # and compare it with the numerically computed values of M(n/t)
Line 137: Line 128:
  
  
<!--Finally, we plot the exact and numerical results for comparison-->
 
 
最後に、厳密解、数値解を比較のためプロットします。
 
最後に、厳密解、数値解を比較のためプロットします。
  
Line 169: Line 159:
 
  plt.show()
 
  plt.show()
  
<!--We see that the magnetization agrees very well to visual accuracy, and approaches the exact scaling form in the relevant limit.-->
 
 
視覚的には磁化が非常によく一致し、関連する制限されている正確なスケーリングフォームに近づくことがわかります。
 
視覚的には磁化が非常によく一致し、関連する制限されている正確なスケーリングフォームに近づくことがわかります。
  
Line 176: Line 165:
  
 
== TEBDの誤差解析:タイムステップ誤差 ==
 
== TEBDの誤差解析:タイムステップ誤差 ==
<!--We now use the exact solution to compute the error in a TEBD simulation as a function of time.  We first investigate the effects of changing the "infinitesimal" time step dt.-->
 
 
時間の関数としてTEBDシミュレーションの誤差を計算するために、厳密解を使用しています。最初に "無限小"時間ステップdtを変化の影響を調査します。
 
時間の関数としてTEBDシミュレーションの誤差を計算するために、厳密解を使用しています。最初に "無限小"時間ステップdtを変化の影響を調査します。
  
Line 214: Line 202:
 
                 })
 
                 })
  
<!--By changing the parameter NUMSTEPS we implicitly change the time step, since the total evolution time TAU is fixed.  We now write the input files, run the simulations, and collect data:-->
 
 
NUMSTEPSパラメータを代えることで、暗黙的にタイムステップを変化させます。それはTAUが固定値なので、結果、タイムステップを変化させることになります。入力ファイルを準備し、計算を実行させ、データを収集します。
 
NUMSTEPSパラメータを代えることで、暗黙的にタイムステップを変化させます。それはTAUが固定値なので、結果、タイムステップを変化させることになります。入力ファイルを準備し、計算を実行させ、データを収集します。
  
Line 224: Line 211:
 
  Magdata=pyalps.load.loadTimeEvolution( pyalps.getResultFiles(prefix='tutorial_2b'), measurements=['Local Magnetization'])
 
  Magdata=pyalps.load.loadTimeEvolution( pyalps.getResultFiles(prefix='tutorial_2b'), measurements=['Local Magnetization'])
  
<!--We now calculate the exact result from the time data, and then calculate the difference between the numerical and the exact result for the magnetization-->
 
 
時間データから正確な結果を計算します。そして、磁化の数値計算と厳密解を比較します。
 
時間データから正確な結果を計算します。そして、磁化の数値計算と厳密解を比較します。
  
Line 236: Line 222:
 
         q[0].y=[abs(q[0].y[syssize/2+1-1]-loc)]
 
         q[0].y=[abs(q[0].y[syssize/2+1-1]-loc)]
  
<!--Finally, we plot this magnetization error:-->
 
 
最後に磁化の誤差をプロットします。
 
最後に磁化の誤差をプロットします。
  
Line 254: Line 239:
 
  plt.show()
 
  plt.show()
  
 
<!--We see that, for short times, the errors are roughly proportional to dt^2, reflecting the contribution to the error from the trotter breakup of our exponential.  At long times, however, the simulations with the smallest dt have errors which become larger than those with larger dt, and eventually the errors blow up!  We will have more to say about this behavior in the next section.-->
 
 
次のことが言えます。短い時間だと、指数関数のトロッター分解からエラーへの寄与を反映し、、誤差はおおまかにdt^2に比例しています。しかし長い時間では、最小のdtでのシミュレーションでは、大きなdtと比較して大きな誤差が発生しており、最終的に誤差は膨大なものになってしまします。次のセクションでさらに追求します。
 
次のことが言えます。短い時間だと、指数関数のトロッター分解からエラーへの寄与を反映し、、誤差はおおまかにdt^2に比例しています。しかし長い時間では、最小のdtでのシミュレーションでは、大きなdtと比較して大きな誤差が発生しており、最終的に誤差は膨大なものになってしまします。次のセクションでさらに追求します。
  
Line 262: Line 245:
  
 
== TEBDの誤差解析:エンタングルメントカットオフ誤差 ==
 
== TEBDの誤差解析:エンタングルメントカットオフ誤差 ==
<!--We now investigate the effects of changing the entanglement cutoff parameter <math> \chi </math> on the errors in the magnetization.-->
 
 
続いて、エンタングルメントカットオフ-パラメータ<math> \chi </math>の誤差への影響を調査します。
 
続いて、エンタングルメントカットオフ-パラメータ<math> \chi </math>の誤差への影響を調査します。
  
Line 334: Line 316:
 
  plt.show()
 
  plt.show()
  
<!--We see that, for short times, the errors are roughly proportional to dt^2, again reflecting the contribution to the error from the trotter breakup of our exponential.  As time increases, however, a cascade of diverging errors ensues.  First the simulation with <math> \chi=10 </math> diverges around <math> t=5</math>, then the simulation with <math> \chi=20 </math> diverges around <math> t=9</math> and so on.  This breakdown is due to the fact that the protocol for finding the matrix product state which best approximates the time-evolved state is approximate when the state becomes highly entangled.  This approximation involves a renormalization of the wavefunction, and so the errors accumulate roughly exponentially in time.-->
 
 
短い時間だと、指数関数のトロッター分解からエラーへの寄与を反映し、、誤差はおおまかにdt^2に比例しています。時間が増加すると、多くの発散誤差が発生します。最初のシミュレーション<math> \chi=10 </math>では<math> t=5</math>あたりで、<math> \chi=20 </math>では<math> t=9</math>で発散します。この時、非常に状態が複雑になっている時の時間発展状態を近似する行列積を求めるための手法が近似式を使用していることに起因します。この近似は、波動関数の繰り込みを含み、誤差が時間内にほぼ指数関数的に蓄積されていってしまいます。
 
短い時間だと、指数関数のトロッター分解からエラーへの寄与を反映し、、誤差はおおまかにdt^2に比例しています。時間が増加すると、多くの発散誤差が発生します。最初のシミュレーション<math> \chi=10 </math>では<math> t=5</math>あたりで、<math> \chi=20 </math>では<math> t=9</math>で発散します。この時、非常に状態が複雑になっている時の時間発展状態を近似する行列積を求めるための手法が近似式を使用していることに起因します。この近似は、波動関数の繰り込みを含み、誤差が時間内にほぼ指数関数的に蓄積されていってしまいます。
  
<!--This exponential growth of errors also accounts for the failure of the simulations with smaller dt.  As dt becomes smaller we must apply the approximate propagation scheme more to reach the same fixed final time, and this means more accumulation of the exponentially growing truncation error.  Thus, we must strike a delicate balance between the error incurred by increasing the time step and the error incurred by taking more time steps.  All results should be carefully checked for convergence in both dt and <math> \chi </math>.-->
 
 
この誤差の指数関数的爆発は、小さなdtのときでも発生しています。dtがより小さい時でも、対応できるような近似伝搬法を採用しなくてはなりません。それは打ち切り誤差の蓄積による指数関数的な爆発を意味します。よって、タイムステップの増加によって発生する誤差と、さらなるタイムステップを取ることによって発生する誤差との間のデリケートなバランスを調整する必要があります。すべての結果は<math> \chi </math>とdtで収束のチェックを注意深くおこなうべきです。
 
この誤差の指数関数的爆発は、小さなdtのときでも発生しています。dtがより小さい時でも、対応できるような近似伝搬法を採用しなくてはなりません。それは打ち切り誤差の蓄積による指数関数的な爆発を意味します。よって、タイムステップの増加によって発生する誤差と、さらなるタイムステップを取ることによって発生する誤差との間のデリケートなバランスを調整する必要があります。すべての結果は<math> \chi </math>とdtで収束のチェックを注意深くおこなうべきです。
  
Line 347: Line 327:
 
== XXZモデルの解法 ==
 
== XXZモデルの解法 ==
  
<!--We saw from the exact solution that the magnetization profile had a well defined front which expanded ballistically with velocity <math> v=1 </math>.  The XX model has many special properties and so it is natural to ask if this same magnetization behavior holds under more general conditions.  In this part of the tutorial we investigate the effects of adding a <math> J_z S_i^z S_{i+1}^z </math> term to the Hamiltonian, corresponding to the XXZ model.  In the limit as this term dominates the spins become frozen in a parallel configuration, and so the initial state becomes an exact eigenstate of the Hamiltonian.  The XX terms in the Hamiltonian try to flip the spins, and are responsible for the propagating magnetization wavefront we saw in the pure XX model.  As a quantitative measure of the ability of the system to transport spin, we consider the integrated flow of magnetization through the center defined in [http://pre.aps.org/abstract/PRE/v71/i3/e036102 Phys. Rev. E 71 036102(2005)] as
+
磁化は、速度<math> v=1 </math>で弾道敵に拡大していくことが厳密解からわかりました。XXモデルでは多くの特別な性質を持っているので、このような磁化の振る舞いが一般的なことかどうか興味が持たれます。このチュートリアルのパートでは、ハミルトニアンに<math> J_z S_i^z S_{i+1}^z </math>を加え、その影響を調査します。これはXXZモデルに相当します。この追加によって、平行構造のスピンは凍結スピンとなり、初期状態はハミルトニアンの厳密な固有状態となります。ハミルトニアンのXX項はスピンを反転しようとし、ピュアなXXモデルでみられたように磁化の波面を伝搬する役割をはたします。スピン輸送する系の特性の定量的な測定として、[http://pre.aps.org/abstract/PRE/v71/i3/e036102 Phys. Rev. E 71 036102(2005)]
  
<math> \Delta M(t)=\sum_{n>L/2}^{L} (\langle S_n^z(t)\rangle+1/2) </math>-->
+
<math> \Delta M(t)=\sum_{n>L/2}^{L} (\langle S_n^z(t)\rangle+1/2) </math>
  
磁化は、速度<math> v=1 </math>で弾道敵に拡大していくことが厳密解からわかりました。XXモデルでは多くの特別な性質を持っているので、このような磁化の振る舞いが一般的なことかどうか興味が持たれます。このチュートリアルのパートでは、ハミルトニアンに<math> J_z S_i^z S_{i+1}^z </math>を加え、その影響を調査します。これはXXZモデルに相当します。この追加によって、平行構造のスピンは凍結スピンとなり、初期状態はハミルトニアンの厳密な固有状態となります。ハミルトニアンのXX項はスピンを反転しようとし、ピュアなXXモデルでみられたように磁化の波面を伝搬する役割をはたします。スピン輸送する系の特性の定量的な測定として、[http://pre.aps.org/abstract/PRE/v71/i3/e036102 Phys. Rev. E 71 036102(2005)]の<math> \Delta M(t)=\sum_{n>L/2}^{L} (\langle S_n^z(t)\rangle+1/2) </math>-->で定義された手法を元に磁化の総合的なフローを考慮します。
+
で定義された手法を元に磁化の総合的なフローを考慮します。
  
  
Line 357: Line 337:
  
 
Pythonスクリプトを用いた実行方法です。tutorial2d.pyを参照してください。
 
Pythonスクリプトを用いた実行方法です。tutorial2d.pyを参照してください。
 
  
 
  import pyalps
 
  import pyalps
Line 391: Line 370:
 
          })
 
          })
  
Note that we are simulating a range of Jz-couplings.  We then write the input files, run the simulation, and get the output as usual:
+
<!--Note that we are simulating a range of Jz-couplings.  We then write the input files, run the simulation, and get the output as usual:-->
 +
Jz-結合のシミュレーションをおこないます。入力ファイルを準備し、計算を実行し、出力結果を得ます。
  
 
  baseName='tutorial_2d'
 
  baseName='tutorial_2d'
Line 411: Line 391:
 
  q[0].y=[0.5*(syssize/2)+sum(q[0].y[syssize/2:syssize])]
 
  q[0].y=[0.5*(syssize/2)+sum(q[0].y[syssize/2:syssize])]
 
   
 
   
 
 
最後に、Jzカップリング幅の総合磁化をプロットします。
 
最後に、Jzカップリング幅の総合磁化をプロットします。
  
Line 425: Line 404:
 
  plt.show()
 
  plt.show()
  
<!--We see that for Jz<1 the integrated magnetization increases roughly linearly in time, and so the magnetization transport is ballistic as in the XX case.  For Jz around 1, we see a change in the qualitative behavior to one in which the integrated magnetization eventually saturates.-->
 
 
Jz>1では、総合磁化はおおよそリニアに増加していきます。その磁化の相転移は、XXの場合のように弾道的(ballistic)です。1のJzでは、統合磁化は最終的に飽和するといった定性的な振る舞いがみられます。
 
Jz>1では、総合磁化はおおよそリニアに増加していきます。その磁化の相転移は、XXの場合のように弾道的(ballistic)です。1のJzでは、統合磁化は最終的に飽和するといった定性的な振る舞いがみられます。
  
Line 432: Line 410:
  
 
=== 質疑  ===
 
=== 質疑  ===
<!--* The point Jz=1 where the behavior of the integrated magnetization undergoes a distinct qualitative change is the point at which the XXZ model transitions from a critical phase to the Antiferromagnetic phase.  However, this phase transition is a priori a low-energy phenomenon, affecting the ground state.  Can you deduce how this low energy change affects the dynamical properties of our high-energy initial state?-->
 
 
* Jz=1での総合磁化が異なる性質に変化する振る舞いをするポイントは、XXZモデルは臨界点から反強磁性相に転移しているということです。しかし、この相転移は、定常状態に影響与える低エネルギーでの先験的現象です。この低エネルギーでの変化がどのように、今回おこなった高エネルギーの初期状態のダイナミックスに影響を与えるのか、推測できますか?
 
* Jz=1での総合磁化が異なる性質に変化する振る舞いをするポイントは、XXZモデルは臨界点から反強磁性相に転移しているということです。しかし、この相転移は、定常状態に影響与える低エネルギーでの先験的現象です。この低エネルギーでの変化がどのように、今回おこなった高エネルギーの初期状態のダイナミックスに影響を与えるのか、推測できますか?

Latest revision as of 04:20, 19 March 2012

磁壁の展開

このチュートリアルでは、非平衡状態で準備したS=1/2スピン鎖の時間発展について学習します。選択した特殊な状態は、次のように、真ん中で上下スピンをもち、左右でそれぞれ上スピン、下スピンを持つような構造です。
| \downarrow \downarrow \dots \downarrow \uparrow \dots \uparrow \uparrow\rangle
この状態は初期状態としてINITIAL_STATEに'kink'と設定されています。いくつかの厳密解は一次元XXモデルに関して知られています。

XX modelの完全解

XXモデルでのkink初期状態の時間発展はJordan-WignerPhys. Rev. E 59, 4912 (1999)の自由フェルミオンへの相転移によって厳密解が解かれています。時間関数として、どんなサイトでも磁化の期待値はBessel関数の合計として表現できることが明らかとなりました。そして、磁化と初期磁壁から大きな距離は、変数n/t でスケーリングします。nは中心からの距離、tは時間です。

 M(n,t)=-\frac{1}{2}\sum_{i=1-n}^{n-1}j_i^2\left(t\right)

\lim_{n\to \infty} \lim_{t\to \infty} M(n,t)\to \phi\left(n/t\right)=-\frac{1}{\pi}\arcsin\left(n/t\right)

ここで、 M(n,t)は磁化で、nは中心からの距離、 j_i(t)はオーダi のBessel関数を表します。このチュートリアルの最初のパートでは、これらの2つの結果を演習します。

Pythonを用いた実行法

Pythonスクリプトを用いた実行方法です。tutorial2a.pyを参照してください。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import copy
import math
import scipy.special
#prepare the input parameters
parms = [{ 
         'L'                         : 50,
         'MODEL'                     : 'spin',
         'local_S'                   : 0.5,
         'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
         'Jxy'                         : 1,
         'INITIAL_STATE' : 'kink',
         'CHI_LIMIT' : 40,
         'TRUNC_LIMIT' : 1E-12,
         'NUM_THREADS' : 1,
         'TAUS' : [20.0],
         'POWS' : [0.0],
         'GS' : ['H'],
         'GIS' : [0.0],
         'GFS' : [0.0],
         'NUMSTEPS' : [500],
         'STEPSFORSTORE' : [2]
       }]

mathとscipy.specialモジュールは、この厳密解を比較するために必要な特別な関数を生成するために使用されます。POWSを0としていることに注意してください。まったくクエンチしないことに対応しています。GS、GIS、GFSの値は任意です。TAUS、NUMSTEPSはそれぞれ、全計算時間とタイムステップを与えます。入力ファイルを準備し、計算を実行し、出力を得て、最後にプロットをおこないます。

baseName='tutorial_2a'
nmlname=pyalps.writeTEBDfiles(parms, baseName)
res=pyalps.runTEBD(nmlname)
#Get the results of the simulation
Data=pyalps.load.loadTimeEvolution(pyalps.getResultFiles(prefix='tutorial_2a'), measurements=['Local Magnetization'])

厳密解と比較するために、計算結果の後処理が必要です。 まず最初に、後処理データを保持するための空の配列を定義します。

#define a dataset numericalSolution to contain the numerical result
numericalResult=[]
#define a dataset exactSolution to contain the exact solution
exactResult=[]
#define a dataset scalingForm to contain the scaling form
scalingForm=[]

時間データから厳密解を計算します。比較するために各サイトの磁化の計算された値を使用します。

 #Compute the exact result M(n,t)=<S_n^z>=-(1/2)*sum_{i=1-n}^{n-1} j_i(t)^2, where
# j_i(t) is the Bessel function of order i and compare to the numerically obtained result
for q in Data:
        syssize=q[0].props['L']
        #Assign a label 'Distance' denoting the distance from the center n (only do the first two sites
        #to avoid cluttering the plot)
        for n in range(1,3):
                #Create copies of the data for postprocessing
                numericalCopy=copy.deepcopy(q)
                exactCopy=copy.deepcopy(q)
                
                numericalCopy[0].props['Distance']=n
                numericalCopy[0].props['SIMID']='Numerical at n='+str(n)
                exactCopy[0].props['Distance']=n
                exactCopy[0].props['SIMID']='Exact at n='+str(n)

                #compute the exact result of the manetization n sites from the center
                loc=0.0
                for i in range(1-n,n):
                        loc-=0.5*scipy.special.jn(i,q[0].props['Time'])*scipy.special.jn(i,q[0].props['Time'])                        
                exactCopy[0].y=[loc]
                #add to the the exact dataset
                exactResult.extend(exactCopy)

                #get the numerical result of the magnetization n sites from the center
                numericalCopy[0].y=[q[0].y[syssize/2+n-1]]
                #add to the the numerical dataset
                numericalResult.extend(numericalCopy)


次に、正確なスケーリング関数を計算します。スケーリング変数の関数として磁化を計算し、 n/t厳密解と比較します。

#compute the scaling form
# \phi(n/t)=-(1/pi)*arcsin(n/t) that M(n,t) approaches as n->infinity and t->infinity
# and compare it with the numerically computed values of M(n/t)
for q in Data:
        syssize=q[0].props['L']
        #Assign a label 'Distance' denoting the distance from the center n (only do the first few sites
        #to avoid cluttering the plot)
        for n in range(0,5):
                #Create a copy of the data for postprocessing
                scalingCopy=copy.deepcopy(q)
                scalingCopy[0].props['Distance']=n

                #The first distance contains the exact scaling form \phi(n/t)=-(1/pi)*arcsin(n/t)
                if n==0:
                        scalingCopy[0].props['Time']=1.0/scalingCopy[0].props['Time']
                        scalingCopy[0].y=[-(1.0/3.1415926)*math.asin(min(scalingCopy[0].props['Time'],1.0))]
	                scalingCopy[0].props['SIMID']='Exact'

                #The other distances contain the numerical data as a function of the scaling variable M(n/t)
                else:
                        scalingCopy[0].props['Time']=n/scalingCopy[0].props['Time']
                        scalingCopy[0].y=[scalingCopy[0].y[syssize/2+n-1] ]
	                scalingCopy[0].props['SIMID']='Numerical at n='+str(n)
                #add to the scaling dataset
                scalingForm.extend(scalingCopy)


最後に、厳密解、数値解を比較のためプロットします。

#Plot the numerical and exact magnetization for comparison
exactMag=pyalps.collectXY(exactResult, x='Time', y='Local Magnetization',foreach=['SIMID'])
for q in exactMag:
	q.props['label']=q.props['SIMID']
numericalMag=pyalps.collectXY(numericalResult, x='Time', y='Local Magnetization',foreach=['SIMID'])
for q in numericalMag:
	q.props['label']=q.props['SIMID']

plt.figure()
pyalps.plot.plot([exactMag, numericalMag])
plt.xlabel('Time $t$')
plt.ylabel('Magnetization')
plt.legend(loc='lower right')
plt.title('Magnetization vs. time')
#Plot the scaling form with the numerical data for comparison
Scal=pyalps.collectXY(scalingForm, x='Time', y='Local Magnetization', foreach=['SIMID'])
for q in Scal:
	q.props['label']=q.props['SIMID']
plt.figure()
pyalps.plot.plot(Scal)
plt.xlabel('Scaling variable $n/t$')
plt.ylabel('Magnetization$(n,t)$')
plt.legend()
plt.xlim(0,1.5)
plt.title('Magnetization scaling function; numerical and exact results')
plt.show()

視覚的には磁化が非常によく一致し、関連する制限されている正確なスケーリングフォームに近づくことがわかります。

Vistrailsを用いた実行法

Vistrailsを用いた実行は、tutorial1a.vtを開き、ワークフロー"tutorial2b"を選択してください。

TEBDの誤差解析:タイムステップ誤差

時間の関数としてTEBDシミュレーションの誤差を計算するために、厳密解を使用しています。最初に "無限小"時間ステップdtを変化の影響を調査します。

Pythonを用いた実行法

Pythonスクリプトを用いた実行方法です。tutorial2b.pyを参照してください。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import numpy as np
import math
import scipy.special
#prepare the input parameters
parms=[]
count=0
for nsteps in [100, 250, 500, 750, 1000]:
       count+=1
       parms.append({ 
                 'L'                         : 50,
                 'MODEL'                     : 'spin',
                 'local_S'                   : 0.5,
                 'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
                 'Jxy'                         : 1,
                 'INITIAL_STATE' : 'kink',
                 'CHI_LIMIT' : 20,
                 'TRUNC_LIMIT' : 1E-12,
                 'NUM_THREADS' : 1,
                 'TAUS' : [20.0],
                 'POWS' : [0.0],
                 'GS' : ['H'],
                 'GIS' : [0.0],
                 'GFS' : [0.0],
                 'NUMSTEPS' : [nsteps],
                 'STEPSFORSTORE' : [int(math.floor(nsteps/100))],
                 'SIMID': count
               })

NUMSTEPSパラメータを代えることで、暗黙的にタイムステップを変化させます。それはTAUが固定値なので、結果、タイムステップを変化させることになります。入力ファイルを準備し、計算を実行させ、データを収集します。

baseName='tutorial_2b_'
nmlnameList=pyalps.writeTEBDfiles(parms, baseName)
res=pyalps.runTEBD(nmlnameList)
#Get magnetization data
Magdata=pyalps.load.loadTimeEvolution( pyalps.getResultFiles(prefix='tutorial_2b'), measurements=['Local Magnetization'])

時間データから正確な結果を計算します。そして、磁化の数値計算と厳密解を比較します。

#Postprocessing-get the exact result for comparison
for q in Magdata:
        syssize=q[0].props['L']
        #Get the exact result of M(1,t)=-(1/2)*(j_0(t)^2), where j_0(t) is the 0^{th} order
        # bessel function and M(1,t) is the magnetization one site to the right of the chain center
        loc=-0.5*scipy.special.jn(0,q[0].props['Time'])*scipy.special.jn(0,q[0].props['Time'])
        #Get the difference between the computed and exact results
        q[0].y=[abs(q[0].y[syssize/2+1-1]-loc)]

最後に磁化の誤差をプロットします。

#Plot the Error in the magnetization one site to the right of the chain center
Mag=pyalps.collectXY(Magdata, x='Time', y='Local Magnetization', foreach=['SIMID'])
for q in Mag:
	dt=round(q.props['TAUS']/q.props['NUMSTEPS'],3)
	q.props['label']='dt='+str(dt)
plt.figure()
pyalps.plot.plot(Mag)
plt.xlabel('Time $t$')
plt.yscale('log')
plt.ylabel('Magnetization Error')
plt.title('Error in the magnetization vs. time')
plt.legend(loc='lower left')
plt.show()

次のことが言えます。短い時間だと、指数関数のトロッター分解からエラーへの寄与を反映し、、誤差はおおまかにdt^2に比例しています。しかし長い時間では、最小のdtでのシミュレーションでは、大きなdtと比較して大きな誤差が発生しており、最終的に誤差は膨大なものになってしまします。次のセクションでさらに追求します。

Vistrailsを用いた実行法

Vistrailsを用いた実行は、tutorial1a.vtを開き、ワークフロー"tutorial2a"を選択してください。

TEBDの誤差解析:エンタングルメントカットオフ誤差

続いて、エンタングルメントカットオフ-パラメータ \chi の誤差への影響を調査します。

Pythonを用いた実行法

Pythonスクリプトを用いた実行方法です。tutorial2c.pyを参照してください。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import math
import scipy.special
#prepare the input parameters
parms=[]
count=0
for chi in [10, 20, 30, 40]:
       count+=1
       parms.append({ 
                 'L'                         : 50,
                 'MODEL'                     : 'spin',
                 'local_S'                   : 0.5,
                 'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
                 'Jxy'                         : 1,
                 'INITIAL_STATE' : 'kink',
                 'CHI_LIMIT' : chi,
                 'TRUNC_LIMIT' : 1E-12,
                 'NUM_THREADS' : 1,
                 'TAUS' : [20.0],
                 'POWS' : [0.0],
                 'GS' : ['H'],
                 'GIS' : [0.0],
                 'GFS' : [0.0],
                 'NUMSTEPS' : [500],
                 'STEPSFORSTORE' : [5],
                 'SIMID': count
               })


入力ファイルを準備し、計算を実行し、出力を得て、最後にプロットをおこないます。

baseName='tutorial_2c_'
nmlnameList=pyalps.writeTEBDfiles(parms, baseName)
res=pyalps.runTEBD(nmlnameList)

#Get magnetization data
Magdata=pyalps.load.loadTimeEvolution( pyalps.getResultFiles(prefix='tutorial_2c'), measurements=['Local Magnetization'])

#Postprocessing-get the exact result for comparison
for q in Magdata:
        syssize=q[0].props['L']
        #Get the exact result of M(1,t)=-(1/2)*(j_0(t)^2), where j_0(t) is the 0^{th} order
        # bessel function and M(1,t) is the magnetization one site to the right of the chain center
        loc=-0.5*scipy.special.jn(0,q[0].props['Time'])*scipy.special.jn(0,q[0].props['Time'])
        #Get the difference between the computed and exact results
        q[0].y=[abs(q[0].y[syssize/2+1-1]-loc)]


最後に、磁化誤差のプロットをおこないます。

#Plot the Error in the magnetization one site to the right of the chain center
Mag=pyalps.collectXY(Magdata, x='Time', y='Local Magnetization', foreach=['SIMID'])
for q in Mag:
	q.props['label']='$\chi$='+str(q.props['CHI_LIMIT'])
plt.figure()
pyalps.plot.plot(Mag)
plt.xlabel('Time $t$')
plt.yscale('log')
plt.ylabel('Magnetization Error')
plt.title('Error in the magnetization vs. time')
plt.legend(loc='lower left')
plt.show()

短い時間だと、指数関数のトロッター分解からエラーへの寄与を反映し、、誤差はおおまかにdt^2に比例しています。時間が増加すると、多くの発散誤差が発生します。最初のシミュレーション \chi=10 では t=5あたりで、 \chi=20 では t=9で発散します。この時、非常に状態が複雑になっている時の時間発展状態を近似する行列積を求めるための手法が近似式を使用していることに起因します。この近似は、波動関数の繰り込みを含み、誤差が時間内にほぼ指数関数的に蓄積されていってしまいます。

この誤差の指数関数的爆発は、小さなdtのときでも発生しています。dtがより小さい時でも、対応できるような近似伝搬法を採用しなくてはなりません。それは打ち切り誤差の蓄積による指数関数的な爆発を意味します。よって、タイムステップの増加によって発生する誤差と、さらなるタイムステップを取ることによって発生する誤差との間のデリケートなバランスを調整する必要があります。すべての結果は \chi とdtで収束のチェックを注意深くおこなうべきです。

Vistrailsを用いた実行法

Vistrailsを用いた実行は、tutorial1a.vtを開き、ワークフロー"tutorial2c"を選択してください。


XXZモデルの解法

磁化は、速度 v=1 で弾道敵に拡大していくことが厳密解からわかりました。XXモデルでは多くの特別な性質を持っているので、このような磁化の振る舞いが一般的なことかどうか興味が持たれます。このチュートリアルのパートでは、ハミルトニアンに J_z S_i^z S_{i+1}^z を加え、その影響を調査します。これはXXZモデルに相当します。この追加によって、平行構造のスピンは凍結スピンとなり、初期状態はハミルトニアンの厳密な固有状態となります。ハミルトニアンのXX項はスピンを反転しようとし、ピュアなXXモデルでみられたように磁化の波面を伝搬する役割をはたします。スピン輸送する系の特性の定量的な測定として、Phys. Rev. E 71 036102(2005)

 \Delta M(t)=\sum_{n>L/2}^{L} (\langle S_n^z(t)\rangle+1/2)

で定義された手法を元に磁化の総合的なフローを考慮します。


Pythonを用いた実行法

Pythonスクリプトを用いた実行方法です。tutorial2d.pyを参照してください。

import pyalps
import matplotlib.pyplot as plt
import pyalps.plot
import math
import scipy.special
#prepare the input parameters
parms=[]
count=0
for z in [0.0, 0.3, 0.9, 1.0, 1.1, 1.5]:
	count+=1
	parms.append({ 
	          'L'                         : 50,
	          'MODEL'                     : 'spin',
	          'local_S'                   : 0.5,
	          'CONSERVED_QUANTUMNUMBERS'  : 'Sz',
	          'Jxy'                         : 1,
	          'Jz'                         : z,
		  'INITIAL_STATE' : 'kink',
		  'CHI_LIMIT' : 40,
		  'TRUNC_LIMIT' : 1E-12,
		  'NUM_THREADS' : 1,
		  'TAUS' : [20.0],
		  'POWS' : [0.0],
		  'GS' : ['H'],
		  'GIS' : [0.0],
		  'GFS' : [0.0],
		  'NUMSTEPS' : [500],
		  'STEPSFORSTORE' : [5],
		  'SIMID': count
	        })

Jz-結合のシミュレーションをおこないます。入力ファイルを準備し、計算を実行し、出力結果を得ます。

baseName='tutorial_2d'
nmlnameList=pyalps.writeTEBDfiles(parms, baseName)
res=pyalps.runTEBD(nmlnameList)
#Get magnetization data
Magdata=pyalps.load.loadTimeEvolution( pyalps.getResultFiles(prefix='tutorial_2d'), measurements=['Local Magnetization'])

計算された磁化データから、上記で定義したように総合的な磁化を計算します。

#Compute the integrated magnetization across the center
for q in Magdata:
	syssize=q[0].props['L']
	#Compute the integrated flow of magnetization through the center \Delta M=\sum_{n>L/2}^{L} (<S_n^z(t)>+1/2)
	#\Delta M= L/4
	loc=0.5*(syssize/2)
	#\Delta M-=<S_n^z(t)> from n=L/2 to L
	q[0].y=[0.5*(syssize/2)+sum(q[0].y[syssize/2:syssize])]

最後に、Jzカップリング幅の総合磁化をプロットします。

#Plot the integrated magnetization
Mag=pyalps.collectXY(Magdata, x='Time', y='Local Magnetization', foreach=['Jz'])
plt.figure()
pyalps.plot.plot(Mag)
plt.xlabel('Time $t$')
plt.ylabel('Integrated Magnetization $\Delta M(t)$')
plt.title('Integrated Magnetization vs. Time')
plt.legend(loc='upper left')
plt.show()

Jz>1では、総合磁化はおおよそリニアに増加していきます。その磁化の相転移は、XXの場合のように弾道的(ballistic)です。1のJzでは、統合磁化は最終的に飽和するといった定性的な振る舞いがみられます。

Vistrailsを用いた実行法

Vistrailsを用いた実行は、tutorial1a.vtを開き、ワークフロー"tutorial2d"を選択してください。

質疑

  • Jz=1での総合磁化が異なる性質に変化する振る舞いをするポイントは、XXZモデルは臨界点から反強磁性相に転移しているということです。しかし、この相転移は、定常状態に影響与える低エネルギーでの先験的現象です。この低エネルギーでの変化がどのように、今回おこなった高エネルギーの初期状態のダイナミックスに影響を与えるのか、推測できますか?