電位差一定のMD計算

基本的なアイデア

　拡張系による圧力制御、温度制御の方法に倣い、MD計算に仮想的なコンデンサを結合させることで電位差の制御を行う。

コンデンサ内の粒子の運動

　左図のようなコンデンサは、コンデンサ間に

$U^{Q-Q}={\frac {dQ^{2}}{2\varepsilon _{0}S}}={\frac {Q^{2}}{2C_{0}}}$

のポテンシャルが働く。またその中の粒子には

$U^{Q-MD}={\frac {Q}{\varepsilon _{0}S}}\sum _{i}q_{i}\left(z_{i}-{\frac {Z_{1}+Z_{2}}{2}}\right)$

のポテンシャルが働く（誘起双極子が存在する場合はその項も追加される。）。従ってこの系のハミルトニアンは、

$H(x,p,Q)=H(x,p)+{\frac {Q}{\varepsilon _{0}S}}\sum _{i}q_{i}\left(z_{i}-{\frac {Z_{1}+Z_{2}}{2}}\right)+{\frac {Q^{2}}{2C_{0}}}=H(x,p)$

(1)

と表される。

NVTVアンサンブル

　上記ハミルトニアンは NVTQ アンサンブルとして系を扱っている。我々が考えたいのは NVTV アンサンブルであるため、これの熱力学関数と分配関数の導出を行う。

完全な熱力学関数

　N,V,T,V を指定した際の完全な熱力学関数は、

$dE=TdS-pdV+VdQ$

のルジャンドル変換を考えることで、

$G=E-TS-QV$

(2)

となる（エネルギー E は電場 E の字体が異なるので注意）。

分布関数

　微視的状態 $i$ 　の出現確率を　 $f_{i}$ 　とすると、NVTV アンサンブルでは、

$\sum _{i}f_{i}=1,\qquad \sum _{i}E_{i}f_{i}=E,\qquad \sum _{i}Q_{i}f_{i}=Q$

が成り立つ。従って上記を拘束条件としてエントロピーを最大にすることを考えると、

$max\left[-k_{B}\sum _{i}f_{i}\ lnf_{i}+{\tilde {\alpha }}\left(\sum _{i}f_{i}-1\right)+{\tilde {\beta }}\left(\sum _{i}E_{i}f_{i}-E\right)+{\tilde {\gamma }}\left(\sum _{i}Q_{i}f_{i}-Q\right)\right]-k_{B}(lnf_{i}+1)+{\tilde {\alpha }}+{\tilde {\beta }}E_{i}+{\tilde {\gamma }}Q_{i}=0$

$f_{i}=exp\left(-1+{\frac {\tilde {\alpha }}{k_{B}}}+{\frac {\tilde {\beta }}{k_{B}}}E_{i}+{\frac {\tilde {\gamma }}{k_{B}}}Q_{i}\right)={\frac {1}{Z}}exp[-\beta (E_{i}-\gamma Q_{i})]$

が得られる。従って、

${\mathcal {S}}=-k_{B}\sum _{i}f_{i}\ lnf_{i}=-k_{B}\sum _{i}f_{i}[-\beta (E_{i}-\gamma Q_{i})-lnZ]={\frac {E-\gamma Q}{T}}+k_{B}lnZ$

$-k_{B}TlnZ=E-TS-\gamma Q$

(3)

となる。式 3 は式 2 に対応するので、

$\gamma =V$

であることがわかる。以上より分配関数は、

$Z(N,V,T,V)=\sum _{i}exp[-\beta (E_{i}-Q_{i}V)]$

(4)

となる。微視的なエネルギーは古典的には式1で表されるので式4の古典極限は、

$Z(N,V,T,V)=\sum _{i}exp[-\beta (E_{i}+Q_{i}V)]=\int dxdpdQ\ exp[-\beta {H(x,p,Q)-QV}]$

(5)

となる。

仮想系のハミルトニアンと運動方程式

　式1 に仮想的なコンデンサの運動エネルギーとポテンシャルを加える形で式 5 を実現する仮想系のハミルトニアンを作成する。今仮想系のハミルトニアンを

$H'(x,p,s,p_{s},Q,p_{Q})=H^{Nose}(x,p',s,p_{s})+{\frac {Q}{\varepsilon _{0}S}}\sum _{i}q_{i}\left(z_{i}-{\frac {Z_{1}+Z_{2}}{2}}\right)+{\frac {Q^{2}}{2C_{0}}}+{\frac {p_{Q}^{2}}{2m_{Q}}}+U_{Q}$

とする。この時分配関数は

$Z\ \propto \int \cdots \int dxdpdsdp_{s}dQdp_{Q}\delta [H'(x,p,s,p_{s},Q,p_{Q})-E']\ \propto \int \cdots \int dxdpdQexp[-\beta \{H(x,p,Q)+U_{Q}\}]$

となる。従ってこれが式5 に一致するには、

$U_{Q}=-QV$

であればよい。つまり、

$H'(x,p,s,p_{s},Q,p_{Q})=H^{Nose}(x,p',s,p_{s})+{\frac {p_{Q}^{2}}{2m_{Q}}}+U(x,Q)$

U(x,Q)={\frac {Q}{\varepsilon _{0}S}}\sum _{i}q_{i}\left(z_{i}-{\frac {Z_{1}+Z_{2}}{2}}\right)+{\frac {Q^{2}}{2C_{0}}}-QV

(6)

である。

　Q の時間発展は正準方程式により、

${\frac {dQ}{dt}}={\frac {p_{Q}}{m_{Q}}}$

${\frac {dp_{Q}}{dt}}=-{\frac {1}{\varepsilon _{0}S}}\sum _{i}q_{i}\left(z_{i}-{\frac {Z_{1}+Z_{2}}{2}}\right)-{\frac {Q}{C_{0}}}+V$

(7)

となる。今瞬間の誘電分極を

$P_{z}(t)={\frac {1}{V}}\sum _{i}q_{i}\left(z_{i}(t)-{\frac {Z_{1}+Z_{2}}{2}}\right)$

(8)

瞬間の比誘電率を、

$\varepsilon (t)={\frac {\varepsilon _{0}E_{z}(t)}{\varepsilon _{0}E_{z}(t)-P_{z}(t)}}={\frac {-\varepsilon _{0}Q(t)/\varepsilon _{0}S}{-\varepsilon _{0}Q(t)/\varepsilon _{0}S-P_{z}(t)}}={\frac {Q(t)}{Q(t)+P_{z}(t)S}}$

と置くと、式7 は、

${\frac {dp_{Q}}{dt}}=-{\frac {P_{z}(t)d}{\varepsilon _{0}}}-{\frac {Q(t)}{C_{0}}}+V=-{\frac {P_{z}(t)S+Q(t)}{C_{0}}}+V=-{\frac {Q(t)}{C_{0}\varepsilon (t)}}+V=V-{\frac {Q}{C(t)}}$

(9)

と表される。ここで、

$C(t)=C_{0}\varepsilon (t)$

とした。式9 は式6 で表される運動が、熱平衡における電位差とコンデンサに蓄えられる電荷の関係式

$Q=CV$

の周りで振動することを意味している。

一般化

本節では上記手法を任意の形状の電場 $E(z)$ が係る系に拡張する。

電場のエネルギー

真空中に $E(z)$ が存在するとき、そのエネルギーは

U^{\mathrm {E} }={\frac {\varepsilon _{0}S}{2}}\int E(z)\cdot E(z)dz

で表される。

周期系への適用

　上記手法を液液界面等の周期系に適用する場合、図のようにセル中の対称な形で荷電平板を設置する必要すればよい。

すると、各粒子には $Z_{2}$ 　に実対称な形で静電ポテンシャルが働き周期境界で静電ポテンシャルが連続となる。しかしながら、このままだとその境界面では電場が不連続になることから、各平板の電荷分布がデルタ関数ではなく、

$\rho _{i}(z)={\frac {Q_{i}}{S}}{\frac {1}{\sqrt {2\pi \sigma ^{2}}}}exp\left[-{\frac {(z-Z_{i})^{2}}{2\sigma ^{2}}}\right]$

のように広がりを持つと考え電場の連続性を確保する。各平板が作る電場、および静電ポテンシャルは、

${\begin{aligned}E_{i}(z,Q_{i})&={\frac {1}{\varepsilon _{0}}}\int dz\rho _{i}(z)={\frac {Q_{i}}{2\varepsilon _{0}S}}erf\left({\frac {z-Z_{i}}{\sqrt {2\sigma ^{2}}}}\right)\equiv {\frac {Q_{i}}{S}}T^{(1)}(z-Z_{i})\\\phi _{i}(z,Q_{i})&=-\int dzE_{i}(z)=-{\frac {Q_{i}}{\varepsilon _{0}S}}\left\{{\frac {1}{\sqrt {2\pi }}}\sigma exp\left[-{\frac {(z-Z_{i})^{2}}{2\sigma ^{2}}}\right]+{\frac {1}{2}}(z-Z_{i})erf\left({\frac {z-Z_{i}}{\sqrt {2\sigma ^{2}}}}\right)\right\}=-{\frac {Q_{i}}{S}}T^{(0)}(z-Z_{i})\end{aligned}}$

であることから、系に働く電場は、

$E(z,Q)=E_{1}(z,Q)+E_{2}(z,-2Q)+E_{3}(z,Q)$

静電ポテンシャルは、

$\phi (z,Q)=\phi _{1}(z,Q)+\phi _{2}(z,-2Q)+\phi _{3}(z,Q)$

となる。以上よりコンデンサ間に働くポテンシャルは、

$U^{Q-Q}=\sum _{i}\sum _{j\neq i}Q_{i}\phi _{j}(z,Q_{j})$

系の中の粒子に働くポテンシャルは、

$U^{Q-MD}=\sum _{i}q_{i}\phi (z_{i},Q)-\sum _{i}\mu _{z,i}E(z_{i},Q)$

により計算される。また、一般化座標としての Q に働くポテンシャルはコンデンサが２つ分と考え、

$U_{Q}=-2QV$

とする。従って拡張系のハミルトニアンは

$H'(x,p,s,p_{s},Q,p_{Q})=H^{Nose}(x,p',s,p_{s})+{\frac {p_{Q}^{2}}{2m_{Q}}}+U(x,Q)$

$U(x,Q)=\sum _{i}q_{i}\phi (z_{i},Q)-\sum _{i}\mu _{z,i}E(z_{i},Q)+{\frac {1}{2}}\sum _{i}\sum _{j\neq i}Q_{i}\phi _{j}(Z_{i},Q_{j})-2QV$

となる。

時間発展

　上記拡張系ハミルトニアンに基づく運動方程式はフーバー形式 $(\zeta \equiv {\dot {s}}/s)$ 　で表すと、

${\begin{aligned}{\frac {dx}{dt}}&={\frac {p}{m}}\\{\frac {dp}{dt}}&=F(x)+qE(z,Q)-\mu _{z}{\frac {\partial E(z,Q)}{\partial x}}-\zeta p\\{\frac {d\zeta }{dt}}&={\frac {p_{\zeta }}{m_{\zeta }}}\\{\frac {dp_{\zeta }}{dt}}&=gk_{B}[T(p)-T_{0}]\\{\frac {dQ}{dt}}&={\frac {p_{Q}}{m_{Q}}}\\{\frac {dp_{Q}}{dt}}&=F_{Q}(x,Q)\\&=\sum _{i}q_{i}{\frac {T^{(0)}(z_{i}-Z_{1})-2T^{(0)}(z_{i}-Z_{2})+T^{(0)}(z_{i}-Z_{3})}{S}}\\&+\sum _{i}\mu _{z,i}{\frac {T^{(1)}(z-Z_{1})-2T^{(1)}(z-Z_{2})+T^{(1)}(z_{i}-Z_{3})}{S}}-\sum _{i}\sum _{j\neq i}{\frac {dQ_{i}}{dQ}}\phi _{j}(Z_{i},Q_{j})+2V\end{aligned}}$

で与えられる。最後の式において

${\begin{aligned}\sum _{i}q_{i}&{\frac {T^{(0)}(z_{i}-Z_{1})-2T^{(0)}(z_{i}-Z_{2})+T^{(0)}(z_{i}-Z_{3})}{S}}\\&+\sum _{i}\mu _{z,i}{\frac {T^{(1)}(z-Z_{1})-2T^{(1)}(z-Z_{2})+T^{(1)}(z_{i}-Z_{3})}{S}}=-{\frac {2P_{z}(x)d}{\varepsilon _{0}}}\end{aligned}}$

に相当する。以上より（疑似的な）リウビル演算子は

$i{\mathcal {L}}={\frac {p}{m}}{\frac {\partial }{\partial x}}\ +\ \left[F(x)+qE(z,Q)+\mu _{z}{\frac {\partial E(z,Q)}{\partial x}}-\zeta p\right]{\frac {\partial }{\partial p}}\ +\ {\frac {p_{\zeta }}{m_{\zeta }}}{\frac {\partial }{\partial \zeta }}\ +\ gk_{B}[T(p)-T_{0}]{\frac {\partial }{\partial p_{\zeta }}}\ +\ {\frac {p_{Q}}{m_{Q}}}{\frac {\partial }{\partial Q}}\ +\ F_{Q}(x,Q){\frac {\partial }{\partial p_{Q}}}$

となる。これを非可換性と計算の利便性によりいくつかの項に分離し

${\begin{aligned}i{\mathcal {L}}&=i{\mathcal {L}}_{1}+i{\mathcal {L}}_{2}+i{\mathcal {L}}_{3}+i{\mathcal {L}}_{4}+i{\mathcal {L}}_{5}+i{\mathcal {L}}_{6}\\i{\mathcal {L}}_{1}&=-\zeta p{\frac {\partial }{\partial p}}\\i{\mathcal {L}}_{2}&=\left[F(x)+qE(z,Q)+\sum _{i}\mu _{z}{\frac {\partial E(z,Q)}{\partial x}}\right]{\frac {\partial }{\partial p}}\\i{\mathcal {L}}_{3}&=F_{Q}(x,Q){\frac {\partial }{\partial p_{Q}}}\\i{\mathcal {L}}_{4}&={\frac {p}{m}}{\frac {\partial }{\partial x}}+{\frac {p_{Q}}{m_{Q}}}{\frac {\partial }{\partial Q}}\\i{\mathcal {L}}_{5}&=gk_{B}[T(p)-T_{0}]{\frac {\partial }{\partial p_{\zeta }}}\\i{\mathcal {L}}_{6}&={\frac {p_{\zeta }}{m_{\zeta }}}{\frac {\partial }{\partial \zeta }}\\\end{aligned}}$

とする。するとトロッタ分解により時間発展は

${\begin{aligned}exp(i{\mathcal {L}}\Delta t)=exp\left(i{\mathcal {L}}_{1}{\frac {\Delta t}{2}}\right)exp&\left(i{\mathcal {L}}_{2}{\frac {\Delta t}{2}}\right)exp\left(i{\mathcal {L}}_{3}{\frac {\Delta t}{2}}\right)exp(i{\mathcal {L}}_{4}\Delta t)exp\left(i{\mathcal {L}}_{3}{\frac {\Delta t}{2}}\right)\\\\&exp\left(i{\mathcal {L}}_{5}{\frac {\Delta t}{2}}\right)exp(i{\mathcal {L}}_{6}\Delta t)exp\left(i{\mathcal {L}}_{5}{\frac {\Delta t}{2}}\right)exp\left(i{\mathcal {L}}_{2}{\frac {\Delta t}{2}}\right)exp\left(i{\mathcal {L}}_{1}{\frac {\Delta t}{2}}\right)\end{aligned}}$

のように表される（ ${\mathcal {L}}_{3}$ 　及び ${\mathcal {L}}_{4}$ 　と ${\mathcal {L}}_{5},{\mathcal {L}}_{6}$ 　が可換なことを利用している。）。 ${\mathcal {L}}_{1},{\mathcal {L}}_{2},{\mathcal {L}}_{5},{\mathcal {L}}_{6}$ 　の部分は力に外部電場の項が加わったことを除き、通常のNVTの計算と同様となる。残りの部分については、

${\begin{aligned}p_{Q}&\gets p_{Q}+{\frac {\Delta t}{2}}F_{Q}(x,Q)\\x&\gets x+{\frac {p}{m}}\Delta t\\Q&\gets Q+{\frac {p_{Q}}{m_{Q}}}\Delta t\\&---\\p_{Q}&\gets p_{Q}+{\frac {\Delta t}{2}}F_{Q}(x,Q)\\\end{aligned}}$

のように時間発展させればよい。なお、 $Q$ 　の時間発展後、 $p_{Q}$ 　の時間発展の前に

相互作用テンソル $T^{(0)}(z_{i}-Z_{j}),T^{(1)}(z_{i}-Z_{j}),T^{(2)}(z_{i}-Z_{j})$ 　を計算する。
コンデンサが作る静電ポテンシャル $\phi (z_{i},Q)$ 　、及び、電場 $E(z_{i},Q)$ 　を計算する。
$\phi (z_{i},Q),E(z_{i},Q),dEdz(z_{i},Q)$ 　を外場としてクーロン相互作用を計算（誘起双極子計算を含む）。
得られた誘起双極子を用いて $F_{Q}(x,Q)$ 　を計算する。

の処理を行う必要がある。また、 $p_{Q}$ 　の時間発展後、 $U(x,Q),{\frac {p_{Q}^{2}}{2m_{Q}}}$ 　等の出力変数の計算も必要になる。

実装

　本計算をFreeFlexに実装する場合、move_struct内に実装するべきである。しかしこの場合問題点として、

coulomb_struct 内にある電荷や誘起双極子が move_struct 内で必要になる。
コンデンサによる外部電場が move_struct 内で計算されるため、計算した電場を coulomb_struct に渡す必要がある。

という2点を解決しなければならない。後者については外部電場を指定するサブルーチンとして set_external_electric_field サブルーチンが coulomb_struct 内に用意されているのでFreeFlex.Fでこれを毎ステップ呼び出せばよい。一方前者についても set_electric_information サブルーチンを用意し、FreeFlex.F 上で呼び出せばとりあえずは解決できる。 ^[1]

↑ （本質的には）電荷や誘起双極子のような一般的なパラメータが coulomb_struct という深い位置の構造体に記憶されていることに問題があると思われる。そのため、これらを coulomb_struct から coord_struct に移すことも考えた。しかしながらこの場合、FEP 計算で電荷を変化させる際に、複製した coulomb_struct を使用している状況も変更しなくてはならない。勿論、coord_struct を複製する形に変更することも考えられるが、この場合は時間発展後の座標のコピーを毎ステップ行わなくてはいけなくなる。この変更の煩雑さを避けるため、上記のように set_electric_infomation を使用する形で実装することにした。

　以上のことから、move_struct 内に次に示す電位差一定のMD計算の構造体 move_NVTV_struct を用意する。

type move_NVTV_struct

   REAL8 :: Q    !!! 一般化座標としてのコンデンサの電荷
   REAL8 :: vQ    !!! 一般化座標の速度
   REAL8 :: mQ    !!! 一般化座標の質量
   REAL8 :: FQ    !!! 一般化座標にかかる力

   REAL8 :: V    !!! 電位差
   REAL8 :: U    !!! コンデンサが関わるポテンシャル。 テキストの U(x,Q) のこと。
   REAL8 :: K    !!! 一般化座標の運動エネルギー。

   REAL8 :: q(1:nsite)    !!! 各サイトの電荷
   REAL8 :: mu(1:nsite)    !!! 各サイトの誘起双極子のz方向成分
   REAL8 :: phi(1:nsite)
   !!! コンデンサが各サイトに作る静電ポテンシャル。テキストの ϕ(z_i,Q) のこと。
   REAL8 :: E(1:nsite)
   !!! コンデンサが各サイトに作る電場。 テキストの E(z_i,Q) のこと。
   REAL8 :: dEdz(1:nsite)
   !!! コンデンサが各サイトに作る電場の座標微分。 テキストの ∂E(z_i,Q)⁄(∂z_i ) のこと
   REAL8 :: T0(1:nsite), T1(1:nsite), T2(1:nsite)
   !!! 相互作用テンソル。

   REAL8 :: phiQ(1:3)
   !!! 各コンデンサ位置の静電ポテンシャル。 テキストの ∑_(j≠i)▒〖ϕ_j (Z_i,Q_j ) 〗 のこと。

   REAL8 :: sigma    !!! 各コンデンサの電荷の広がり
   integer :: iCOM    !!! コンデンサを置く位置の重心番号。テキストの Z_2 に対応。
   REAL8 :: S    !!! コンデンサの面積

contains

   subroutine init(…)
   !!! 初期設定
   endsubroutne

   subroutine set_electric_info(q,mu)
   !!! 電荷と誘起双極子の値を受け取るサブルーチン
     REAL8 :: q(1:nsite)
     REAL8 :: mu(1:3,1:nsite)
   endsubroutine

   subroutine vQ_update1(coord,dt)
   !!! vQ の時間発展を行うサブルーチン
     type(coord_struct) :: coord
     REAL8 :: dt
   endsubroutine

   subroutine Q_update(coord,dt)
   !!! Q の時間発展を行うサブルーチン。
   !!! 同時に相互作用テンソル、phi、E、dEdz、phiQの計算も行う。
     type(coord_struct) :: coord
     REAL8 :: dt
   endsubroutine

   subroutine vQ_update2(coord,dt)
   !!! vQ の時間発展を行うサブルーチン。
   !!! 同時に FQ、及び出力変数の計算も行う。
     type(coord_struct) :: coord
     REAL8 :: dt
   endsubroutine

endtype

</source>

[1] （本質的には）電荷や誘起双極子のような一般的なパラメータが coulomb_struct という深い位置の構造体に記憶されていることに問題があると思われる。そのため、これらを coulomb_struct から coord_struct に移すことも考えた。しかしながらこの場合、FEP 計算で電荷を変化させる際に、複製した coulomb_struct を使用している状況も変更しなくてはならない。勿論、coord_struct を複製する形に変更することも考えられるが、この場合は時間発展後の座標のコピーを毎ステップ行わなくてはいけなくなる。この変更の煩雑さを避けるため、上記のように set_electric_infomation を使用する形で実装することにした。

[1]