Merge pull request #116 from issp-center-dev/cTPQ

.github/workflows/main.yml

@@ -36,7 +36,9 @@ jobs:
           brew install openmpi scalapack
 
       - name: pip
-        run: python -m pip install numpy
+        run: |
+          python -m pip install numpy
+          python3 -m pip install numpy
 
       - name: make workspace
         run: cmake -E make_directory ${{runner.workspace}}/build

doc/en/X.sh

@@ -0,0 +1,2 @@
+make latexpdf
+open build/latex/userguide_HPhi_en.pdf 

doc/en/source/algorithm/TPQ_en.rst

@@ -5,43 +5,80 @@ Finite-temperature calculations by the TPQ method
 
 Sugiura and Shimizu showed that
 it is possible to calculate the finite-temperature properties
-from a few wavefunctions (in the thermodynamic limit, only one wave function is necessary) [#]_ .
+from a few wavefunctions (in the thermodynamic limit, only one wave function is necessary) [1]_, [2]_ .
 The wavefunction is called the thermal pure quantum (TPQ) state.
 Because the TPQ state can be generated by operating the Hamiltonian 
 to the random initial wavefunction,
 we directly use the routine Lanczos method for the TPQ calculations.
-Here, we explain how to construct the micro canonical TPQ (mTPQ) state,
+The TPQ state is given by
+
+.. math::
+
+   |\Phi(\beta)\rangle\equiv\exp(-\frac{\beta\hat{\mathcal H}}{2})|\Phi_{\rm rand}\rangle.
+
+Here,  :math:`\beta` represents inverse temperature and :math:`|\Phi_{\rm rand}\rangle`
+is a random initial vector.
+It is shown that it is possible to
+calculate physical quantities at inverse temperature  :math:`\beta` 
+as expectation values with respect to :math:`|\Phi(\beta)\rangle`.
+
+Details of implementation
+-------------------------
+
+**Construction of the micro canonical TPQ (mTPQ) state**
+^^^^^^^^^^^^^^^^^^
+Here, we explain how to construct the micro canonical TPQ (mTPQ) state [1]_,
 which offers the simplest method for calculating finite-temperature properties.
 
-Let :math:`|\psi_{0}\rangle` be a random initial vector.
+Let :math:`|\Phi_{\rm rand}\rangle` be a random initial vector.
 By operating :math:`(l-\hat{\mathcal H }/N_{s})^{k}`\(:math:`l` is constant and :math:`N_{s}` represents the number of sites) 
-to :math:`|\psi_{0}\rangle`,
+to :math:`|\Phi_{\rm rand}\rangle`,
 we obtain the :math:`k`\th TPQ states as
 
 .. math::
 
-   |\psi_{k}\rangle \equiv \frac{(l-\hat{\mathcal H }/N_{s})|\psi_{k-1}\rangle}{|(l-\hat{\mathcal H }/N_{s})|\psi_{k-1}\rangle|}.
+   |\Phi_{k}\rangle \equiv \frac{(l-\hat{\mathcal H }/N_{s})|\Phi_{k-1}\rangle}{|(l-\hat{\mathcal H }/N_{s})|\Phi_{k-1}\rangle}.
 
-From :math:`|\psi_{k}\rangle`, we estimate the corresponding inverse temperature :math:`\beta_{k}` as
+From :math:`|\Phi_{k}\rangle`, we estimate the corresponding inverse temperature :math:`\beta_{k}` as
 
 .. math::
 
-   \beta_{k}\sim \frac{2k/N_{s}}{l-u_{k}}, u_{k} = \langle \psi_{k}|\hat{\mathcal H }|\psi_{k}\rangle/N_{s},
+   \beta_{k}\sim \frac{2k/N_{s}}{l-u_{k}}, u_{k} = \langle \Phi_{k}|\hat{\mathcal H }|\Phi_{k}\rangle/N_{s},
 
 where :math:`u_{k}` is the internal energy.
 The arbitrary local physical properties at :math:`\beta_{k}` are also estimated as
 
 .. math::
 
-   \langle \hat{A}\rangle_{\beta_{k}} =  \langle \psi_{k}|\hat{A}|\psi_{k}\rangle/N_{s}.
+   \langle \hat{A}\rangle_{\beta_{k}} =  \langle \Phi_{k}|\hat{A}|\Phi_{k}\rangle.
 
 In a finite-size system,
 error is caused by the choice of the initial random vector.
 To estimate the average value and error of the physical properties,
-we perform some independent calculations by changing :math:`|\psi_{0}\rangle`.
+we perform some independent calculations by changing :math:`|\Phi_{\rm rand}\rangle`.
+Since the temperature
+
+**Construction of the canonical TPQ (cTPQ) state**
+^^^^^^^^^^^^^^^^^^
+Here, we explain how to construct the canonical TPQ (cTPQ) state [2]_,
+which is another way to construct the TPQ state.
+In the cTPQ method, :math:`\exp[-\beta\hat{\mathcal H}/2]` is 
+approximated as follows:
+
+.. math::
+
+  &\exp(-\frac{\beta\hat{\mathcal H}}{2})|\Phi_{\rm rand}\rangle\sim|\Phi_{k}\rangle = U(\Delta\tau)^{k}|\Phi_{\rm rand}\rangle\\
+  &U(\Delta\tau)=\sum_{n=0}^{n_{\rm max}}\frac{1}{n!}(-\frac{\Delta\tau}{2}\hat{\mathcal H})^{n}\\
+  &\beta_{k}=k\Delta \tau
+   
+As in the mTPQ method,
+physical quantities are calculated as 
+
+.. math::
+
+   \langle \hat{A}\rangle_{\beta_{k}} =  \frac{\langle\Phi_{k}|\hat{A}|\Phi_{k}\rangle}{\langle\Phi_{k}|\Phi_{k}\rangle}.
+
 
-Details of implementation
--------------------------
 
 **Initial vector**
 ^^^^^^^^^^^^^^^^^^
@@ -50,12 +87,13 @@ For the TPQ method, the initial vector is given by using a random generator, i.e
 
 .. math::
 
-   123432+(n_{\rm run}+1)\times  |r_s|+k_{\rm Thread}+N_{\rm Thread} \times k_{\rm Process},
+   123432+(n_{\rm run}+1)\times|r_s|+k_{\rm Thread}+N_{\rm Thread} \times k_{\rm Process},
 
-where :math:`r_s` is the number given by an input file and :math:`n_{\rm run}` is the number of runs. :math:`r_s` and the maximum value of :math:`n_{\rm run}` are defined by ``initial_iv`` and ``NumAve`` in an input file for Standard mode or a ModPara file for Expert mode, respectively. Random numbers are generated by using SIMD-oriented Fast Mersenne Twister (dSFMT) [#]_ . We can select the type of initial vector as a real number or complex number by using ``InitialVecType`` in a ModPara file.
+where :math:`r_s` is the number given by an input file and :math:`n_{\rm run}` is the number of runs. :math:`r_s` and the maximum value of :math:`n_{\rm run}` are defined by ``initial_iv`` and ``NumAve`` in an input file for Standard mode or a ModPara file for Expert mode, respectively. Random numbers are generated by using SIMD-oriented Fast Mersenne Twister (dSFMT) [3]_ . We can select the type of initial vector as a real number or complex number by using ``InitialVecType`` in a ModPara file.
 :math:`k_{\rm Thread}, N_{\rm Thread}, and k_{\rm Process}` indicate 
 the thread ID, number of threads, process ID, respectively;
 the initial vector depends both on ``initial_iv`` and the number of parallelisms.
 
-.. [#] \S. Sugiura, A. Shimizu, Phys. Rev. Lett. **108**, 240401 (2012).
-.. [#] \http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT.
+.. [1] \S. Sugiura and A. Shimizu, Phys. Rev. Lett. **108**, 240401 (2012).
+.. [2] \S. Sugiura and A. Shimizu, Phys. Rev. Lett. **111**, 010401 (2013).
+.. [3] \http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT.

doc/en/source/filespecification/expertmode_en/CalcMod_file_en.rst

@@ -68,9 +68,11 @@ The parameters correlated with the keywords are as follows.
    | **Description :** Select the method for calculation from the
      following list:
    | 0: Lanczos method
-   | 1: Analysis of the physical properties by using TPQ
-   | 2: Full diagonalization method.
-   | 3: LOBCG for the ground state.
+   | 1: mTPQ method
+   | 2: Full diagonalization method
+   | 3: LOBCG for the ground state
+   | 4: Time-evolution
+   | 5: cTPQ method
 
 *  ``CalcModel``
 
@@ -108,9 +110,12 @@ The parameters correlated with the keywords are as follows.
 
    **Type :** Int (default value: 0)
 
-   | **Description :** Select the type of an initial vector:
-   | 0: Complex type
-   | 1: Real type.
+   | **Description :** Select the type of an initial vector (:math:`v0`):
+   | -1: Real part (:math:`{\rm Re}[v0]]`) and imaginary part  (:math:`{\rm Re}[v0]]`) of the initial 
+    vector are give as the normally distributed random numbers. Thus, the normalized initial vectors are uniformly distributed
+    on the :math:`N_{\rm H}` dimensional super sphere (:math:`N_{\rm H}` is the dimension of the Hilbert space). 
+   | 0: Complex type (:math:`{\rm Re}[v0]\in[-1:1]`, :math:`{\rm Im}[v0]\in[-1:1]` ).
+   | 1: Real type (:math:`{\rm Re}[v0]\in[-1:1]`, :math:`{\rm Im}[v0]=0`).
 
 *  ``OutputEigenVec``
 

doc/en/source/filespecification/expertmode_en/ModPara_file_en.rst

@@ -231,7 +231,7 @@ CG method
 
 
 
-TPQ method
+TPQ (mTPQ/cTPQ) method
 ~~~~~~~~~~
 
 *  ``Lanczos_max``
@@ -247,7 +247,9 @@ TPQ method
    **Type :** Double
 
    **Description :** An integer giving :math:`l` of 
-   :math:`l-{\mathcal H}/N_{s}` used in the TPQ method.
+   :math:`l-{\mathcal H}/N_{s}` used in the mTPQ method.
+   In the cTPQ method,  :math:`l` is used as the width of the imaginary-time evolution, i.e.,
+   :math:`\Delta\tau=1/l`. See :ref:`Ch:algorithm` for the details of the cTPQ method.
 
 *  ``NumAve``
 
@@ -264,6 +266,18 @@ TPQ method
      calculating the correlation functions in the TPQ method.
    | **Note:** A small interval increases the time cost of calculations.
 
+*  ``ExpandCoef``
+
+   **Type :** Int (positive integer)
+
+   **Description :** An integer giving the expansion order :math:`n_{\rm max}` for
+   cTPQ method;
+
+   .. math:: U(\Delta\tau) = \sum_{n=0}^{n_{\rm max}}\frac{1}{n!}\left(-\frac{\Delta\tau}{2}\mathcal{H}\right)^n .
+
+   See :ref:`Ch:algorithm` for the details of the cTPQ method.
+
+
 
 
 Calculating dynamical Green’s functions
@@ -328,9 +342,11 @@ Real time evolution method
    **Type :** Int (positive integer)
 
    **Description :** An integer giving the expansion order :math:`n` for
-   real time evolution method;
+   real-time evolution method;
+
+   .. math:: \exp\left(-i \mathcal{H} \Delta t \right) = \sum_{i=0}^{N}\frac{1}{n!}\left(-i \mathcal{H} \Delta t \right)^n.
+
 
-   .. math:: \exp\left(-i \mathcal{H} \Delta t \right) = \sum_{i=0}^{N}\frac{1}{n!}\left(-i \mathcal{H} \Delta t \right)^n .
 
 *  ``ExpecInterval``
 

doc/en/source/filespecification/standardmode_en/Parameters_for_the_numerical_condition.rst

@@ -93,11 +93,13 @@ Parameters for the numerical condition
 
    **Type :** Double (the default value is provided below)
 
-   **Description :** (Only for TPQ) :math:`l` as :math:`l-{\mathcal H}/N_{s}`
-   is used in the TPQ calculation. Usually, the largest eigenvalue of
+   **Description :** (Only for TPQ) :math:`l` as :math:`(l-{\mathcal H}/N_{s})`
+   is used in the mTPQ calculation. Usually, the largest eigenvalue of
    the Hamiltonian is used as :math:`l`. Thus, the default value of
    :math:`l` is taken as the summation of the absolute values of each
    coefficient in the Hamiltonian divided by the number of sites.
+   In the cTPQ calculation, ``LargeValue`` is used as :math:`\Delta\tau=1/LargeValue`.
+   The definition of :math:`\Delta\tau` is explained in :ref:`Ch:algorithm`.
 
 *  ``NumAve``
 

doc/ja/X.sh

@@ -0,0 +1,2 @@
+make latexpdf
+open build/latex/userguide_HPhi_ja.pdf 

doc/ja/source/algorithm/TPQ_ja.rst

@@ -4,44 +4,76 @@
 ----------------------------------
 
 杉浦・清水によって、 少数個（サイズが大きい場合はほぼ一つ）の
-波動関数から有限温度の物理量を計算する方法が提案されました [1]_ 。
+波動関数から有限温度の物理量を計算する方法が提案されました [1]_ [2]_ 。
 その状態は熱的純粋量子状態(TPQ)と呼ばれています。
 TPQはハミルトニアンを波動関数に順次作用させて得られるので、
 Lanczos法の技術がそのまま使うことができます。
+TPQ状態は次のように与えられます。
+
+.. math::
+
+   |\Phi(\beta)\rangle\equiv\exp(-\frac{\beta\hat{\mathcal H}}{2})|\Phi_{\rm rand}\rangle.
+
+ここで :math:`\beta` は逆温度、:math:`|\Phi_{\rm rand}\rangle`  はランダムな初期ベクトルです。
+:math:`|\Phi(\beta)\rangle` の期待値として、有限温度の物理量が計算できることが
+示されています.
+
+実際の実装について
+~~~~~~~~~~~~~~~~~~
+
+mTPQ状態の構成について
+^^^^^^^^^^^^^^^^^^^^^^^^^^
 ここでは、とくに計算が簡単な, micro canonical TPQ(mTPQ)の
-概要を述べます。
+概要を述べます [1]_。
 
-:math:`|\psi_{0}\rangle`\ をあるランダムベクトルとします。
+:math:`|\Phi_{\rm rand}\rangle`\ をあるランダムベクトルとします。
 これに\ :math:`(l-{\mathcal H}/N_{s})`\ (:math:`l`\ はある定数、\ :math:`N_{s}`\ はサイト数)を\ :math:`k`\ 回作用させた
 （規格化された）ベクトルは次のように与えられます。
 
 .. math::
 
    \begin{aligned}
-   |\psi_{k}\rangle \equiv \frac{(l-{\mathcal H}/N_{s})|\psi_{k-1}\rangle}{|(l-{\mathcal H}/N_{s})|\psi_{k-1}\rangle|}.\end{aligned}
+   |\Phi_{k}\rangle \equiv \frac{(l-{\mathcal H}/N_{s})|\Phi_{k-1}\rangle}{|(l-{\mathcal H}/N_{s})|\Phi_{k-1}\rangle|}.\end{aligned}
 
-この\ :math:`|\psi_{k}\rangle`\ がmTPQ状態で、このmTPQ状態に対応する逆温度\ :math:`\beta_{k}`\ は
+この\ :math:`|\Phi_{k}\rangle`\ がmTPQ状態で、このmTPQ状態に対応する逆温度\ :math:`\beta_{k}`\ は
 以下のように内部エネルギー\ :math:`u_{k}`\ から求めることができます。
 
 .. math::
 
    \begin{aligned}
    \beta_{k}\sim \frac{2k/N_{s}}{l-u_{k}},~~
-   u_{k} = \langle \psi_{k}|{\mathcal H}|\psi_{k}\rangle/N_{s}.\end{aligned}
+   u_{k} = \langle \Phi_{k}|{\mathcal H}|\Phi_{k}\rangle/N_{s}.\end{aligned}
 
-そして、任意 [2]_ の物理量\ :math:`\hat{A}`\ の\ :math:`\beta_{k}`\ での平均値は
+そして、任意 [3]_ の物理量\ :math:`\hat{A}`\ の\ :math:`\beta_{k}`\ での平均値は
 
 .. math::
 
    \begin{aligned}
-   \langle \hat{A}\rangle_{\beta_{k}} =  \langle \psi_{k}|\hat{A}|\psi_{k}\rangle/N_{s}\end{aligned}
+   \langle \hat{A}\rangle_{\beta_{k}} =  \langle \Phi_{k}|\hat{A}|\Phi_{k}\rangle/N_{s}\end{aligned}
 
 となります。 有限系では最初の乱数ベクトルによる誤差がありますので、
 いくつか独立な計算を行って、\ :math:`|\psi_{0}\rangle`
 に関する平均値および標準偏差を見積もっています。
 
-実際の実装について
-~~~~~~~~~~~~~~~~~~
+cTPQ状態の構成について
+^^^^^^^^^^^^^^^^^^^^^^^^^^
+カノニカルTPQ(cTPQ)状態の構成方法について述べます [2]_,
+cTPQ法では :math:`\exp[-\beta\hat{\mathcal H}/2]` は
+次にように近似されます.
+
+.. math::
+
+  &\exp(-\frac{\beta\hat{\mathcal H}}{2})|\Phi_{\rm rand}\rangle\sim|\Phi_{k}\rangle = U(\Delta\tau)^{k}|\Phi_{\rm rand}\rangle\\
+  &U(\Delta\tau)=\sum_{n=0}^{n_{\rm max}}\frac{1}{n!}(-\frac{\Delta\tau}{2}\hat{\mathcal H})^{n}\\
+  &\beta_{k}=k\Delta \tau
+   
+mTPQと同じように物理量はcTPQ状態の期待値として計算できます。
+
+.. math::
+
+   \langle \hat{A}\rangle_{\beta_{k}} =  \frac{\langle\Phi_{k}|\hat{A}|\Phi_{k}\rangle}{\langle\Phi_{k}|\Phi_{k}\rangle}.
+
+
 
 初期ベクトルの設定について
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -60,11 +92,12 @@ Lanczos法の技術がそのまま使うことができます。
 で与えられます。ここで、\ :math:`n_{\rm run}`\ はrunの回数であり、runの最大回数はスタンダードモード用入力ファイル、
 もしくはModParaで指定される入力ファイルの\ ``NumAve``\ で指定します。
 ``initial_iv``\ はスタンダードモード用の入力ファイル、もしくはエキスパートモードではModParaで指定される入力ファイルで指定します。乱数はSIMD-oriented
-Fast Mersenne Twister(dSFMT)を用い発生させています [3]_ 。
+Fast Mersenne Twister(dSFMT)を用い発生させています [4]_ 。
 また、\ :math:`k_{\rm Thread}, N_{\rm Thread}, k_{\rm Process}`\ 
 はそれぞれスレッド番号、スレッド数、プロセス番号を表します。
 したがって同じ\ ``initial_iv``\ を用いても、並列数が異なる場合には別の初期波動関数が生成されます。
 
 .. [1] \S. Sugiura, A. Shimizu, Phys. Rev. Lett. **108**, 240401 (2012).
-.. [2] 局所的にという条件がつきます。
-.. [3] \http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT.
+.. [2] \S. Sugiura and A. Shimizu, Phys. Rev. Lett. **111**, 010401 (2013).
+.. [3] 局所的にという条件がつきます。
+.. [4] \http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/SFMT.

doc/ja/source/filespecification/expertmode_ja/CalcMod_file_ja.rst

@@ -64,10 +64,11 @@ CalcModファイル
 
    | **説明 :** 計算手法の指定を行います。
    | 0: Lanczos法
-   | 1: TPQを利用した有限温度解析
+   | 1: mTPQ計算
    | 2: 全対角化
    | 3: LOBCG法
    | 4: 実時間発展
+   | 5: cTPQ計算
    | を選択することができます。
 
 -  ``CalcModel``
@@ -100,9 +101,10 @@ CalcModファイル
 
    **形式 :** int型 (デフォルト値 0)
 
-   | **説明 :** 固有ベクトルの初期値の型の指定を行います。
-   | 0: 複素数
-   | 1: 実数
+   | **説明 :** 固有ベクトル(:math:`v0`)の初期値の型の指定を行います。
+   | -1:実部、虚部ともに正規分布乱数で与えられます。したがって,規格化した初期ベクトルは :math:`N_{\rm H}` 次元の超球面上に一様分布したものとなります。
+   | 0: 複素数(:math:`{\rm Re}[v0]\in[-1:1]`, :math:`{\rm Im}[v0]\in[-1:1]` ).
+   | 1: 実数(:math:`{\rm Re}[v0]\in[-1:1]`, :math:`{\rm Im}[v0]=0`).
    | から選択することが出来ます。
 
 -  ``OutputEigenVec``

doc/ja/source/filespecification/expertmode_ja/ModPara_file_ja.rst

@@ -220,7 +220,9 @@ TPQ法で使用するパラメータ
    **形式 :** double型 (実数)
 
    **説明 :**
-   TPQで使用する\ :math:`l-\mathcal H/N_{s}`\ の\ :math:`l`\ を指定する整数。
+   TPQで使用する\ :math:`l-\mathcal H/N_{s}`\ の\ :math:`l`\ を指定する実数。
+   cTPQ計算では, :math:`l` は虚時間発展の刻み幅として使われます(:math:`\Delta\tau=1/l`). 
+   cTPQ計算の詳細は :ref:`Ch:algorithm` を参照してください.
 
 -  ``NumAve``
 
@@ -236,6 +238,16 @@ TPQ法で使用するパラメータ
    相関関数の計算を何回のTPQステップおきに行うかの指定する整数。
    頻度を上げると計算コストが増大するので注意してください。
 
+*  ``ExpandCoef``
+
+   **Type :** Int (自然数)
+
+   **Description :** cTPQ計算における展開次数を指定する整数です(:math:`n_{\rm max}`);
+
+   .. math:: U(\Delta\tau) = \sum_{n=0}^{n_{\rm max}}\frac{1}{n!}\left(-\frac{\Delta\tau}{2}\mathcal{H}\right)^n .
+
+   詳しくは, :ref:`Ch:algorithm`を参照してください.
+
 
 
 動的グリーン関数の計算で使用するパラメータ

doc/ja/source/filespecification/standardmode_ja/Parameters_for_the_numerical_condition.rst

@@ -87,6 +87,8 @@
    **説明 :**
    (TPQ法のみで使用) :math:`l-\mathcal{H}/N_{s}`\ の\ :math:`l`\ 。
    ハミルトニアンの各項の係数の絶対値の総和をサイト数で割ったものがデフォルト値になります。
+   cTPQ計算では, ``LargeValue`` は虚時間発展の刻み幅として使われます( :math:`\Delta\tau=1/LargeValue`).
+   :math:`\Delta\tau` の定義はアルゴリズムの章を見てください(:ref:`Ch:algorithm`) 。
 
 *  ``NumAve``