Monitor System

&monitor を FreeFlex 本体の入力ファイル中のネームリストに用いることで、 MD 計算で時々変化する変数（座標、速度、力、エネルギー等）をファイルに出力することができる。 出力ファイルを tools/analysis.exe に読み込ませることで、 グラフ表示や平均値、分布、時間相関関数等の計算を統一的に行うことができる。

Namelist

&monitor には以下のパラメータを指定できる。

Name	Type	Default	Description
interval	INTEGER	10	データを取得するMDステップ間隔。
file	STRING		出力するファイル名。 必須。
condition	STRING		出力変数の指定（詳細は後述）。必須。
separator	CHARACTER	“$”	出力条件文のセパレータ。 Default では ”$” であるが、他の文字を指定して、condition 中で使うことができる。
overwrite	LOGICAL	.false.	新しい出力を行う際にファイルの末尾にデータに追加するか (.false.) 、ファイルを上書きするか (.true.) を決めるパラメータ。
write_interval	integer	0(intervalと同値)	データを出力するMDステップ間隔。intervalより値が小さい、またintervalの値を約数に持たない場合、エラーとなる。特にoverwrite=trueでファイルへの出力が律速の場合有効である。

また、&monitor を FreeFlex 本体の入力ファイル中に複数置くことで、それぞれ別のファイルに別のデータを出力することができる。 設置できる &monitor ネームリストの最大数は、別途 &FreeFlex 中の max_monitor_number 変数で指定する（default: 1）。

条件指定

出力する変数は condition 変数に

condition = “ keyword1  option1  $ keyword2  option2  $ …”

の形で指定する。 指定できるキーワードとオプションは下記のとおり。 keyword によって出力が一行になるか複数行になるかが異なる。 一行で出力できるキーワードに限り separator による複数指定が可能であり、 同一ファイル中に複数のデータを出力できる。

一行出力

SPCOORD_Q id1 id2 … :

id のリストで指定されたユーザー定義座標。出力は

Q(id1) Q(id2) …

のように指定されたid の数だけ出力される。

SPCOORD_F id1 id2 … :

id のリストで指定されたユーザー定義座標方向に働く一般化力。出力は

F(id1) F(id2) …

のように指定されたid の数だけ出力される。 なお、一般化座標 ξ_i に対応する一般化力は、

$F_{\xi _{i}}\equiv -{\frac {\partial U}{\partial \xi _{i}}}=-\sum _{j}{\frac {\partial x_{j}}{\partial \xi _{i}}}{\frac {\partial U}{\partial x_{j}}}$

で定義される。 ここで$x_{j}$は系全体のデカルト座標で、${\frac {\partial x_{j}}{\partial \xi _{i}}}\equiv \mathbb {J} ^{+}$は疑似逆行列の計算

$\mathbb {J} ^{+}\equiv (\mathbb {JM} ^{-1}\mathbb {J} ^{T})^{-1}\mathbb {JM} ^{-1},\qquad J_{ij}\equiv {\frac {\partial \xi _{i}}{\partial x_{j}}},\qquad M_{ij}\equiv m_{i}\delta _{ij}$

により計算される。 これは指定されたユーザー定義座標同士以外が直交している場合に相当する。

REUSMK_qU id1 id2 … :

id のリストで指定されたユーザー定義座標およびバイアスポテンシャルの和。出力は

iRep Q(id1) Q(id2) U …

のように、レプリカにかかるバイアスの番号、指定されたid の数だけの座標、バイアスポテンシャルの和の順に出力される。

GROUP_F COM_id:

COM_id で指定した原子団の重心(COM)に働く力。重心の指定は、&addCOM を参照のこと。

Fx(COM_id) Fy(COM_id) Fz(COM_id) のようにx, y, zの3次元分が出力される。

SITE_POSITION site_id:

site_id で指定したサイトの座標。

COM_POSITION COM_id:

COM_id で指定した重心の座標。

SITE_VELOCITY site_id:

site_id で指定したサイトの速度。

SITE_FORCE site_id:

site_id で指定したサイトに働く力。 拘束力は考慮していない。外場をかけた計算では、重心に働く力をゼロにするよう差し引いていることに注意。

U

系全体のポテンシャルエネルギー。LJ、クーロン、バイアスのポテンシャルと、全ての和$U=U_{LJ}+U_{Coulomb}+U_{bias}$が出力される。

U U_LJ U_Coulomb U_bias

GROUP_M COM_id:

COM_id で指定した原子団の双極子モーメント[C・m]。Minimum image convention に基づいて計算された永久双極子と誘起双極子がそれぞれ出力される。

Mx_perm(COM_id) My_perm(COM_id) Mz_perm(COM_id) Mx_ind(COM_id) My_ind(COM_id) Mz_ind(COM_id)

TOTAL_M

系全体の双極子モーメント[C・m]。 全電荷が非ゼロのグループは含まれない。 Mx_perm My_perm Mz_perm Mx_ind My_ind Mz_ind

Conducting boundary condition を用いた Ewald 計算では、以下の式から系の（光学誘電率を含まない）誘電率を計算することができる。

光学誘電率の計算には系の分極率を計算する必要がある（未実装）。

$\epsilon =1+{\frac {4\pi }{3}}{\frac {\langle M^{2}\rangle }{Vk_{\mathrm {B} }T}}$

• Neumann, M. Dipole moment fluctuation formulas in computer simulations of polar systems. Molecular Physics, 1983, 50, 841-858.

wbond :

input:
condition = "wbond COM1 COM2 wcorr"
wcorr: 単位はÅ。COM1～COM2までの分子を対象にして分子間距離を計算して
3Å以下を結合していると判定する。
COM1をイオンとして扱いイオン-水間の距離はwcorrで補正する。
水-水間だけを対象とする場合はwcorrを0を与えればよい。

output:
edge, eff
水和構造を解析する。水和クラスターの結合の本数をedge、頂点の個数あたりのedgeをeffで出力する。

PRESS, PRESSCON, PRESSBIAS:

タイムステップごとの圧力テンソルを出力する。

このうちPRESSは運動エネルギー由来の項、ビリアル項(coulomb+LJ)を出力する。

また、PRESSCONは拘束力由来の項を、PRESSBIASはバイアス由来（分子内ポテンシャルを含む）の項を出力する。

EQ_OTHER_WN:

SPCOORD の keyword に "WATFING_N" を指定したときにのみ使える。 出力の形は

$w_{n}$ と $w_{n+1}$ が等しいかの真偽値 $c_{1}$ $c_{2}$ ... $c_{n}$ $c_{n+1}$

となっている。ここで $n$ は SPCOORD の keyword に "WATFING_N" を指定したときの option の active_N で指定した整数である。 $w_{n}$, $c_{n}$ の詳しい定義は この資料を参照。

SITE_RELATED_WN:

SPCOORD の keyword に "WATFING_N" を指定したときにのみ使える。 出力は $w$, $w_{n}$, $w~{'}$, $c_{1}$, $c_{2}$ ... $c_{n}$ のそれぞれの長さ、 構成するサイトのサイト番号 $i$, $j$ を左から順に出力する。

ここで $n$ は SPCOORD の keyword に "WATFING_N" を指定したときの option の active_N で指定した整数である。 $w_{n}$, $w~{'}$, $c_{n}$ の詳しい定義は この資料を参照。

以下未作成

SITE_SITE_DISTANCE site_id1 site_id2: サイト-サイト間距離。

複数行出力

phi_z nline COM_id0 site_id1 site_id2:

z軸方向の電位（V）。点電荷による電位（phi_q）および点双極子による電位（phi_m）が出力される。 conducting boundary conditionを用いる場合には系全体の電荷が中性でない場合を許容する（background chargeを考慮した静電ポテンシャルを計算する）。

z phi_q(z) phi_m(z)

conducting boundary conditionに基づいたEwald法において、以下の式により静電ポテンシャルが計算される。

$\phi (z)=\phi (0)-{\frac {4\pi }{V}}{\Biggl (}M_{z}-{\frac {QL_{z}}{2}}{\Biggr )}z-4\pi \int _{0}^{z}dz^{\prime }(L_{z}-z^{\prime })\rho (z^{\prime })+{\frac {2\pi Qz^{2}}{V}}$

ここで系の全電荷および全双極子は以下の式で定義され、全双極子は周期境界条件に基づいて計算される。

構文解析に失敗 (構文エラー): {\displaystyle Q = \int_0^{L_z} dz \rho(z) \\ M_z = \int_0^{L_z} dz \rho(z)z }

open boundary condition（surface termを用いた計算）では、上式に次式の項が加わる。

$-4\pi {V}M_{z}^{surf}z$

ここで$M_{z}^{surf}$は周期境界条件を適用せずに計算された全双極子である。

binの幅dzはz軸方向のセル長を(nline)分割するように決まる。

分割する際に中心に取るための重心を引数COM_id0として指定する必要がある。 site_id1, site_id2によりヒストグラムに含めるサイトを指定することができる。 すなわち site_id1=1, site_id2=100 とした場合には 1～100番目のサイトのみを考慮する。 重心番号とサイト番号を混同しないように注意。

出力ファイルを小さくするため、各ステップまでの平均値を出力するようにしており、 overwrite=.true.とともに用いることを推奨する。 従って、平衡化のステップを含めないように注意する必要がある。

phi_multi nline COM_id0 COM_id1 COM_id2:

テスト用であり、通常はphi_zを用いること。

分子の多重極（四重極）により計算されたz軸方向の電位（V）。

conducting boundary conditionに基づいたEwald法での計算にのみ対応。

系全体が中性でない場合には更に拡張が必要。

binの幅dzはz軸方向のセル長を(nline)分割するように決まる。

分割する際に中心に取るための重心を引数COM_id0として指定する必要がある。

COM_id1, COM_id2によりヒストグラムに含める分子を指定することができる。

すなわち COM_id1=1, COM_id2=100 とした場合には 1～100番目の分子のみを考慮する。

出力ファイルを小さくするため、各ステップまでの平均値を出力するようにしており、 overwrite=.true.とともに用いることを推奨する。

TOTAL_zrho COM_id0 COM_id1 COM_id2:

z-dz/2～z+dz/2に含まれる分子の密度 (kg/m^3)。

binの幅dzはz軸方向のセル長を100分割するように決まる。

分割する際に中心に取るための重心を引数COM_id0として指定する必要がある。

COM_id1, COM_id2によりヒストグラムに含める分子を指定することができる。

すなわち COM_id1=1, COM_id2=100 とした場合には 1～100番目の分子のみを考慮する。

TOTAL_zM nline COM_id0 COM_id1 COM_id2:

z-dz/2～z+dz/2に含まれる分子から計算された平均双極子 (C・m)。

binの幅dzはz軸方向にセルサイズをnline分割したものになる。

分割する際に中心に取るための重心を引数COM_id0として指定する必要がある。

ここで、COM_id0 = 0とすると、中心を絶対座標でとることができる。

（外部電場の値が絶対座標に依存している場合などに用いる）

COM_id1, COM_id2によりヒストグラムに含める分子を指定することができる。

すなわち COM_id1=1, COM_id2=100 とした場合には 1～100番目の分子のみを考慮する。

&addCOMを用いると重心をユーザー変数として分子番号の最後に追加できる。（詳しくはtutorial7を参照）

例）condition="TOTAL_zM 100 877 1 523" (全分子数＋1 水の分子番号の範囲（1~523))　

この例ではtutorial1のような系での水を中心とした水の平均双極子が計算できる。

TOTAL_qrho nbin COM_id0 COM_id1 COM_id2:

z 方向の電荷分布および誘起分極分布をヒストグラムとして出力する。

COM_id0 で指定された重心の z 座標を基準として重心番号 COM_id1 から COM_id2 に属するサイトの最近接のクローンを対象とする。

nbin でヒストグラムのビンの個数を設定できる。

出力ファイルの最初に表示されるラベル Mq Mmu Qq はそれぞれの列が

Qq: 電荷分布 $\rho (z)$

Mq: 電荷分布の積分値 $-\int _{-\infty }^{z}dz'\rho (z')$

Mmu: 誘起分極分布 $P(z)$

であることを意味している。ここで電荷分布の積分値に 負符号 - がついているのは

$\Delta \phi =-{\frac {1}{\varepsilon _{0}}}\int _{-\infty }^{\infty }{\biggl (}\int _{-\infty }^{z}\rho (z')dz'-P(z){\biggr )}dz$

に従って電位を計算する際にそのまま使用できるからと考えられる。

gt:

系全体のgauche/trans比を計算する。
1,2-dichloroethaneバルクを調べる用。
input: 
condition = "gt natom q1 q2 q3 q4"
output:
g, t

gt_z:

z-dz/2～z+dz/2に分子の重心があるときの各binでのgauche/trans比を計算する。
input:
condition = "gt_z natom q1 q2 q3 q4 id0 nline id1 id2 id3 id4"
natom: 分子を構成する原子数
q1～q4: 二面角を構成する原子の名前
id0: 基準にとる重心番号
id1: 開始点 (重心番号)
id2: 終了点 (重心番号)
id3: 開始点 (サイト番号)
id4: 終了点 (サイト番号)

output:
z, g, t
g: gaucheの比
t: transの比
n: binに含まれる分子数

(例)水/1,2-dichloroethane界面でのgauche/transの分布を調べる。

water 523個/ DCE 283個
condition = "gt 8 Cl1 C1 C2 Cl2 807 50 524 806 1570 3833"
binは50くらい少なめにとった方が収束しやすい。

PRESS_z, PRESSCON_z, PRESSBIAS_z:

各ビンにおける圧力テンソルを出力する。ここで出力される行数はFreeFlexのnbin_pressで決まる。

こちらはoverwrite = .true.とすることで、毎ステップ平均を取って更新していく。

また、PRESS_zは運動エネルギー由来の項、ビリアル項(coulomb+LJ)を出力する。

PRESSCONは拘束力由来の項を、PRESSBIASはバイアス由来（分子内ポテンシャルを含む）の項を出力する。

PRESSEXT:

外部電場由来の圧力テンソルを出力する。ここで出力される行数はFreeFlexのnbin_pressで決まる。

z,P_ext

overwrite = .true. 推奨。ここではz方向に電場をかけることを想定しているので、zz方向のみ出力される。

dielectric,dielectric_elect:

波数ベクトル依存の電気感受率を計算するプログラム。計算時に波数 $k=2\pi n/L_{z}$ のnを指定することで計算する波数の範囲を指定することができる。

電気感受率は以下の計算式で表すことができ、dielectricで第一項の配向分極成分、第二項で電子分極成分を計算できる。

なお、dielectric内ではバルク系とスラブ系で計算が異なっており、バルク系では電気感受率の対角成分のみを計算し、スラブ系では非対角項を含めた平衡からのずれを計算しているため 用途に応じてプログラムを切り替える必要がある。

$\chi ^{0}(k_{1},k_{2})={\frac {\langle M(k_{1,z})M(-k_{2,z})\rangle }{Vk_{\mathrm {B} }T}}+{\frac {\langle A_{zz}(k_{1,z},k_{2,z})\rangle }{V}}$

k1,k2,chi

overwrite = .true. 推奨。詳細は メディア:誘電体論new.pdf　に書かれている。dielectric_electはかなり計算が重く、メモリを大量に要求するためコンパイル方法を変更する必要がある。

total_Q,total_zQ:

四重極計算に関するプログラム。total_Qは一分子あたりの四重極のトレースを出力し、total_zQは指定分子の四重極密度を出力する。

total_zQは仮想セル分割数(nline)、重心番号(id0)、特定分子のmsファイル上の分子番号(id1,id2)を指定し、重心をz=0として四重極密度を計算する。

N_HYDRATION:

spcoord が watfing もしくは watfing_n の時、自動に検出しそのとき、spcoor_core 内のwatfing.F もしくは watfing_n.F の中で計算されている nwsh を出力する。 nwsh : 輸送分子の水和数

出力の形式

monitor システムで出力されるファイルには、冒頭にネームリスト &monitor_output および指定した出力変数のタグ行が出力される。 ネームリスト &monitor_output は、このファイルを後続の解析プログラムが読み込む際に引き渡す情報が出力され、また続くタグ行では出力内容に対応するタグが空白区切りで出力される。

ネームリスト &monitor_output の内容

 version: 出力内容のバージョン(*1)
 timestep: 出力の時間刻み（単位：s）
 condition: 条件指定の文字列(*2)
 separator: セパレータ(*2)

(*1) 本 monitor システムと解析用の analysis システムの今後のバージョンアップに備えて、整合性をとるため。

(*2) 上の &monitor ネームリストの内容と同じ。

複数種のデータの出力

FreeFlex 中で複数個の変数をモニターすることができる。複数個の指定には、

(i) ネームリスト &monitor が複数個ある。

(ii) ネームリスト内の condition に、セパレータで区切ったキーワードが複数個ある。

の 2 通りがある。

(i) 複数個の &monitor ネームリストを使うときには別々のファイルに出力する。 &monitor ネームリストは、最大 &FreeFlex 中の max_monitor_number の数まで置くことができる。 なお、file に同じファイルにしないようにすること。 同じファイル名を指定した場合、それぞれのデータが別の行に出力され、順序もバラバラになってしまう。

(ii)&monitor 内でセパレータで区切って複数個を指定すると、同じ出力ファイル内の1行で複数のデータが出力される。

プログラム構造と出力の追加

プログラムの構造は spcoord や constraint2 と類似しており、

のような構造をしている。 モニター変数全体は monitors_struct 構造体で管理される。 その中で monitor_struct 構造体は個々の &monitor ネームリストを管理する。 セパレータごとの変数は、monitor_data_struct 構造体に対応する。

出力の追加は xxx.F に monitor_core_struct の拡張構造体を作成し、monitor_data_struct 中のポインタ指定を設定すればよい。 拡張構造体では、

subroutine init(condition)

condition の内容の処理

character(len=256) function get_tag() result(tag)

出力に応じたタグの指定

character(len=256) function get_data() result(data(:))

データに対応する文字列を返す。 配列のサイズは出力する行数に対応する。

の 3 つのプロシージャを再定義する必要がある。 なお、データは FreeFlex_module を読み込むことで取得する。

複数行の出力を行いたい場合には上記の 3 つのプロシージャの作成の他に拡張構造体中で、 monitor_core_struct%nline に出力行数を指定する必要がある。

具体的な実装の流れは、

①monitor_core/内にプログラムを用意する。

②monitor_data.Fに構造体を追加する。

③Makefileにプログラム名を追加する。

となっている。