Conductor struct

概要

conductor_struct は導体分子の設定とMD中の分極計算を管理する構造体。

変数

conductor_struct の主な構成変数を以下の表に示す。

項目	変数の型	Description
nconductor	INTEGER	系に含まれる導体の総数。ネームリスト &conductor の個数に等しい。
xi(1:coord_core%nsite)	DOUBLE_PRECISION	導体モデルにおけるサイト電荷の広がり $\xi$（Å）。導体以外のサイトも点電荷（$\xi$ = 0.0）として含む。
tensor(1:nconductor)	tensor_struct (後述)	導体1つ分の相互作用テンソルを管理する構造体。

上の変数で用いられる構造型 tensor_struct は主に以下の構成変数を持つ。

項目	変数の型	Description
charge	DOUBLE PRECISION	導体の全電荷。
xi	DOUBLE_PRECISION	導体モデルにおけるサイト電荷の広がり $\xi$（Å）。
nsite	INTEGER	各導体のサイト数
index(1:nsite)	INTEGER	導体のサイト番号から系全体のサイト番号への変換インデックス。
Tij(1:nsite+1,1:nsite+1), Tijinv(:,:)	DOUBLE PRECISION	導体内部の相互作用行列とその逆行列。

Tij は次節で示す行列方程式の係数行列であり、Tijinv はその逆行列である。 MD の最初にこれらを一度求めておき、各ステップでは Tijinv のみを用いて等電位条件を満たす導体サイトの分極を求める。

計算の詳細

導体の各サイト上の電位が等しくなるようにサイト電荷が分極する導体モデル 。 原論文では導体サイト電荷を変数に加えた拡張ラグランジアンの方法が採用されているが、後述のアルゴリズムでは各ステップで厳密に等電位条件が満たされる。

導体の $i$ 番目のサイトに働く電位が以下の式で与えられているとする。

$V_{i}=V_{i}^{0}+\sum _{j(\neq i)}^{n}{\frac {f(r_{ij})Q_{j}}{r_{ij}}}+\alpha Q_{i}\quad (i=1,2,\cdots ,n)$

ここで $V_{i}^{0}$ は注目する導体以外の分子によって作られる電位、$r_{ij}$ は導体内部のサイト $i$, $j$ 間の距離を表す。 右辺第3項は自己相互作用項であり、係数 $\alpha$ はガウス型の電荷分布の広がりパラメータ $\xi$ を用いて以下の式で与えられる。

$\alpha ={\frac {1}{{\sqrt {\pi }}\xi }}$

また、短距離のダンピング係数は以下の式で表される。

$f(r_{ij})=\mathrm {erf} \left({\frac {r_{ij}}{2\xi }}\right)$

原論文 のEq. (8) 右辺第2項に対応するが、係数が異なる。

導体サイトの電荷は等電位条件

$V_{1}=V_{2}=\cdots =V_{n}=-\lambda$

を満たすように決定される。ただし全電荷の一定条件

$\sum _{i}^{n}Q_{i}=Q$

が拘束条件として課される。以上の条件を満たすサイト電荷は以下の行列方程式を解くことによって得られる。

${\begin{aligned}{\begin{pmatrix}\alpha &{\frac {f(r_{12})}{r_{12}}}&\cdots &{\frac {f(r_{1n})}{r_{1n}}}&1\\{\frac {f(r_{21})}{r_{21}}}&\alpha &\cdots &{\frac {f(r_{2n})}{r_{2n}}}&1\\\vdots &\vdots &\ddots &\vdots &\vdots \\{\frac {f(r_{n1})}{r_{n1}}}&{\frac {f(r_{n2})}{r_{n2}}}&\cdots &\alpha &1\\1&1&\cdots &1&0\end{pmatrix}}{\begin{pmatrix}Q_{1}\\Q_{2}\\\vdots \\Q_{n}\\\lambda \end{pmatrix}}={\begin{pmatrix}-V_{1}^{0}\\-V_{2}^{0}\\\vdots \\-V_{n}^{0}\\Q\end{pmatrix}}\end{aligned}}$

このとき、係数行列は導体のサイト座標のみによって決まるため、導体が固定されている場合は MD の中で係数行列が不変となる。 FreeFlex では導体の座標を常に固定し、係数行列とその逆行列を MD のはじめに一度だけ求めておく。

系に分極する溶媒や2個以上の導体が含まれる場合には、全ての分極電荷/双極子がSCF条件を満たすように繰り返し計算によって分極が決定される。 複数個の導体を含む場合は係数行列を拡張することでSCF計算によらずサイト電荷を決めることができるが、未対応。 係数行列は PME 法を用いて求められ、計算オーダーは導体のサイト数 $n$, 逆格子のグリッド数 $N_{\mathrm {grid} }$ に対して $O(nN_{\mathrm {grid} }\mathrm {log} N_{\mathrm {grid} })$。

出力

conductor_structはconductor.Fにまとめられており、以下のサブルーチンから構成される。

 subroutine read_conductor(@)
 !!!
 !!! ネームリスト &conductor を読み込み、初期設定を行う。
 !!!
 
 subroutine calc_matrix(@)
 !!!
 !!! 導体分子内部の相互作用行列を計算する。
 !!! 導体の座標が固定されていることを仮定し、MDの最初に一度だけ求める。
 !!!
 
 subroutine calc_matrix_real(@)
 !!!
 !!! 導体分子内のサイトによるEwald実空間項と自己相互作用項の計算。
 !!!
 
 subroutine calc_matrix_recip(@)
 !!!
 !!! 導体分子内のサイトによるEwald逆空間項とsurface termの計算。
 !!!
 
 subroutine inverse(@)
 !!!
 !!! calc_matrix で求めた相互作用行列の逆行列を求める。
 !!!
 
 subroutine solv(@)
 !!!
 !!! MD の各ステップで等電位条件から導体の電荷を設定する。
 !!!
 
 subroutine set_charge(@)
 !!!
 !!! 導体の総電荷を &conductor で指定した値に補正する。
 !!!
 
 subroutine correct_phi(@)
 !!!
 !!! 導体サイトの自己相互作用を電位に加える。
 !!!
 
 subroutine correct_tensor(@)
 !!!
 !!! 通常の Ewald 構造体で点電荷同士の相互作用として計算された相互作用テンソルを
 !!! ガウス分布の相互作用に補正するためのサブルーチン。
 !!!