はじめに

分子動力学法における数値積分法の時間反転対称性とシンプレクティック性についてのまとめ。

用語の定義など

演算子の指数分解とシンプレクティック性

簡単のため、一次元系$(p,q)$を考える。運動方程式$\{ \dot{p}, \dot{q}\}$と初期条件$\{p(0), q(0)\}$が与えられた時、時刻$t$後の座標が知りたい。微分方程式から、この系のリュービル演算子は

i \mathcal{L} = \dot{p} \frac{\partial }{\partial p} + \dot{q} \frac{\partial }{\partial q}

と書くことができる¹。リュービル演算子を使って、時間を$h$だけ進める時間発展演算子は

U(h) = \exp\left(ih \mathcal{L} \right)

と定義される。しかし、一般に上式を厳密に評価するのは難しいので、

i \mathcal{L} = i \mathcal{L}_K + \mathcal{L}_V

のように分離し、

U_V(h) = \exp\left(ih \mathcal{L}_V \right) = 1 + ih \mathcal{L}_V

などと、指数関数の肩に演算子を載せたものが厳密に評価できる部分に分割したものを利用して数値積分を構築するのがシンプレクティック積分であった。

いちいちmathcalとか書くのが面倒なので、以下、

\begin{align}
X &= i \mathcal{L}_K \\
Y &= i \mathcal{L}_V 
\end{align}

とする。こう書いた時のLie-Trotter公式(トロッター分解)

\exp(h (X+Y)) = \left( \exp(hX/n) \exp(hY/n)  \right)^n + O(h^2/n)

については以前の記事を参照のこと。

さて、トロッター分解の最低次は$n=1$の場合、すなわち、

\exp(h (X+Y)) \sim \exp(hX) \exp(hY)

である。左の指数関数をバラして、$X$と$Y$の係数を比較すれば、これが指数関数のテイラー展開としては1次まで一致していることがわかる。トロッター分解は、$n$を大きくすると元の演算子の近似としては$1/n$で向上していくが、指数関数のテイラー展開の次数という意味では1次までしか一致しない。これを、係数を工夫して高次まで一致させることを考える。

例えば、

\exp(h (X+Y)) \sim \exp(h_1 X) \exp(h_2 Y) \exp(h_3 X) \exp(h_4 Y)

とし、左辺と右辺のテイラー展開の係数が二次まで正しい公式を作る。これは非線形連立方程式であり、かつ解も一意ではないため解くのは容易ではないが、天下り式に解の一つを与えると、

(h_1, h_2, h_3, h_4) = \left(\frac{h}{2}, h, \frac{h}{2}, 0\right)

となる(参考文献参照)。これを代入すると、

\exp(h (X+Y)) \sim \exp(h X/2) \exp(h Y) \exp(h X/2)

これは二次の対称分解と呼ばれる。$X,Y$として運動項、ポテンシャル項のどちらを選ぶかにより、それぞれVelocity Verlet, Position Verlet法と呼ばれる数値積分公式に対応する。

数値積分がシンプレクティックであるとは、(近似された)時間発展演算子$U(h)$によって$(p,q)$が$(P,Q)$に変換された時、シンプレクティック二形式が保存すること、すなわち変換のヤコビアンが1であるということであった。

先程のように、演算子の指数分解を用いて時間発展演算子を構築した場合、そのビルディングブロックである$\exp(t A)$などがシンプレクティックであることから、その組み合わせも必ずシンプレクティックとなる²。従って、指数分解公式で作られた数値積分は(そのビルディングブロックがシンプレクティックである限り)必ずシンプレクティックとなる。

時間反転対称性

数値積分においては、度々時間反転対称性が問題となる³。時間反転対称性とは、ある程度時間発展させた状態で、時間を逆転させた場合、ちゃんと初期条件に戻る、という条件である。これは、時間を$h$だけ発展させる時間発展演算子$U(h)$を用いれば、

U^{-1}(h) = U(-h)

と定義される。ハミルトンの運動方程式は時間反転対称性を満たしているため、厳密な時間発展演算子を作れば、もちろん時間反転対称性は満たされる。しかし、離散化により近似された時間発展演算子が時間反転対称性を満たしているかどうかは自明ではない。

指数分解公式のビルディングブロックを以下のように定義しよう。

\begin{align}
\exp(h (X+Y)) &\equiv U(h) \\
\exp(h X) &\equiv U_X(h) \\
\exp(h Y) &\equiv U_Y(h) 
\end{align}

ビルディングブロックはユニタリ演算子なので、それぞれ$U^{-1}_X(h) = U_X(-h), U^{-1}_Y(h) = U_Y(-h)$を満たしている。

さて、指数分解の一次の公式を考える。

U(h) \sim U_X(h) U_Y(h) \equiv U_1(h)

右辺のinverseをとると

\begin{align}
U_1(h) &= (U_X(h) U_Y(h))^{-1} \\
 &= U_Y^{-1}(h) U_X^{-1}(h) \\
&= U_Y(-h) U_X(-h) \\
&\neq U_X(-h)U_Y(-h)
\end{align}

$U_X(h)$と$U_Y(h)$は非可換なので$U_1^{-1}(h) \neq U_1(-h)$である。すなわち、一次のシンプレクティック積分公式は、シンプレクティックではあるが時間反転対称ではない。

同様に二次の対称分解を考える。

U(h) \sim U_X(h/2) U_Y(h) U_X(h/2) \equiv U_2(h)

これが時間反転対称性$U_2^{-1}(h) = U_2(-h)$を満たすのは自明であろう。

以上を適当な系でちゃんと確認してみよう、というのが本稿の趣旨である(前置き超長かった)。

シンプレクティック積分については以前の記事も参照のこと。

確認

調和振動子

一次元調和振動子にシンプレクティック積分を適用すると、線形シンプレクティック変換になるので教育上都合がよろしい。というわけで、まずは一次のシンプレクティック積分$(p,q) \rightarrow (P,Q)$を作ってみる。

In[1]:= P := p - q h
In[2]:= Q := q + P h
In[3]:= P
Out[3]= p - h q
In[4]:= Q
Out[4]= q + h (p - h q)

念のためシンプレクティック性の確認。

In[5]:= Det[D[{P, Q}, {{p, q}}]]
Out[5]= 1

シンプレクティックですね。

次に、時間反転対称性の確認。逆変換を考えるということは$(P,Q)$を$(p,q)$について解けば良い。

r := Solve[{R == P, S == Q}, {p, q}]

適当な記号がなかったので$(R,S)$で代用。答えがrにリストとして入ってるので、それぞれ比較してみる。

In[11]:= r[[1]][[1]]
Out[11]= p -> R - h^2 R + h S
In[12]:= P
Out[12]= p - h q
In[13]:= r[[1]][[2]]
Out[13]= q -> -h R + S
In[14]:= Q
Out[14]= q + h (p - h q)

明らかに$h \rightarrow -h$という変換で一致しない。すなわち、この時間積分はシンプレクティックだが時間反転対称ではない。

次に二次の対称分解を試す。

In[1]:= tp := p - q h /2
In[2]:= Q := q + tp h
In[3]:= P := tp - Q h /2
In[4]:= Det[D[{P, Q}, {{p, q}}]]
Out[4]= 1

$p(t + h/2)$をtpと表記し、シンプレクティック条件まで確認した。これを逆に解いてみる。

In[14]:= Collect[r[[1]][[1]], {R, S}]
Out[14]= p -> 1/4 (4 - 2 h^2) R + 1/4 (4 h - h^3) S
In[12]:= Collect[P, {p, q}]
Out[12]= (1 - h^2/2) p + (-h + h^3/4) q
In[15]:= Collect[r[[1]][[2]], {R, S}]
Out[15]= q -> -h R + 1/2 (2 - h^2) S
In[16]:= Collect[Q, {p, q}]
Out[16]= h p + (1 - h^2/2) q

それぞれちゃんと書き下して並べてみると

\begin{align}
P &= \left( 1 - \frac{h^2}{2} \right) p + \left(-h + \frac{h^3}{4} \right)q \\
p &= \left( 1 - \frac{h^2}{2} \right) P + \left(h - \frac{h^3}{4} \right)Q \\
Q &= -h p + \left( 1- \frac{h^2}{2} \right) q \\
q &=  h P + \left( 1- \frac{h^2}{2} \right) q
\end{align}

$h \rightarrow -h$という変換で、一致することが確認できる。

まとめ

数値積分のシンプレクティック性と時間反転対称性について確認した。離散化された時間発展演算子において、一般にはシンプレクティック性と時間反転対称性は等価ではない。例えば一次のシンプレクティック積分はシンプレクティックだが時間反転対称ではない。

ここでは調和振動子を考えたが、これだと変換が線形になってしまうので簡単すぎで、いちいちMathematicaのSolveとか使わなくても、変換行列の逆行列を考えれば十分である。本当は変換が線形にならない例、例えば$H = p^2/2 + q^4/4$みたいな系で似たようなことをやろうと思っていたのだがここで力尽きた。

ムネンアトヲタノム

参考文献

Symplectic数値積分法
あとは「吉田　シンプレクティック」とか「Yoshida Symplectic」で検索かければいろいろ見つかるんじゃないかな。

本当はちゃんと使い分けるべきですが、気分でリュービル演算子に負符号つけたりつけなかったりしています。すみません。 ↩
行列で考えればわかりやすい。演算子がシンプレクティックであるとは、行列でいえば行列式が1であるということであり、行列式が1となるような行列の積として得られる行列の行列式もやはり1となる。 ↩
実際の分子動力学法においては離散化した時間発展における厳密な時間反転対称性を気にする人はほとんどいないが、まぁ気にする人はいる。 ↩