GenIce

A swiss army knife to generate proton-disordered ice structures for GROMACS.

View the Project on GitHub vitroid/GenIce

[TOC]

GenIceについて

氷の中では、水分子はそれぞれが格子点の位置を占めており、隣接する水分子の間は水素結合がつないでいますが、その水素結合の向きに乱雑さがあります。氷の結晶はこの意味で単一の結晶構造ではなく、多数の異なる水素結合配列のアンサンブルとみなせます。氷の物性を評価する場合には、これらの構造を適切に生成し、アンサンブル平均をとるとともに、配置エントロピー(Pauling エントロピー)を考慮する必要があります。

GenIceは、氷の幾何学的・トポロジー的規則性(アイスルール)を満たした結晶構造のアンサンブルを生成するツールです。GenIceが生成する氷の構造は、全分極が0であることが保証されます(0でない構造を生成することも可能です)。

GenIceには、生成した構造が確かに十分ランダム化されているかどうかを検定するツールも含まれています。

GenIceは、クラスレートハイドレートの結晶構造を生成することもできます。この時、ケージの中に一定割合でゲスト分子(の混合物)を入れることができます。

さらに、氷の格子の一部をイオンで置きかえた、doped iceを生成することもできます。

GenIceは多様なデータ形式を出力できます。また、ユーザーがプラグインを追加して、新しいフォーマットに対応することも容易です。

氷の結晶構造の生成手順

7つのステージ

GenIceは、単位胞の構造情報から、7つのステージをへて氷の大きな結晶構造を生成します。各ステージでは、処理のあとにFormatプラグイン内で定義されたHook関数が呼びだされ、この関数が出力データを作りだします。

Hook関数が定義されていない場合には、不必要な処理は行われません。例えば、指定したFormatプラグインがHook 1 (水分子の重心位置だけに基いた処理)だけを定義している場合には、Stage 2以降は実行されません。

genice command

下の図にはおおまかな処理の流れが描かれています。実線がデータの流れ、破線は付帯するパラメータの流れやプログラムの流れを表しています。また、各オプションがどのプロセスに作用するかを示してあります。

genice data flow

analice command

下の図にはおおまかな処理の流れが描かれています。実線がデータの流れ、破線は付帯するパラメータの流れやプログラムの流れを表しています。analiceは複数のファイルを読みこんで処理を行い、複数のファイルを出力します。

analiceでは、7つのステージのうちいくつかは必要な場合のみ呼びだされます。読み込んだ座標ファイルが、個々の原子の座標まで含んでいる場合には、ステージ2、3、5は実行されませんが、重心位置のみ(あるいは酸素位置のみ)を含んでいる場合にはこれらのステージにより分子配向が決定されます。各ステージの処理が実行されるかどうかにかかわらず、Hookは実行されます。

genice data flow

分子のラベリング

GenIce内部では、水もゲストも分子として扱われ、0から順にラベルされます。analiceで座標ファイルを読みこんだ場合には、最初の水分子が0、その次の水分子が1、という風にラベル付けされます。

氷の生成

氷の結晶構造を作る場合には以下の情報が必要です。

  1. どの結晶構造を作るのか。
  2. どの水分子モデルを使うのか。(--waterオプション)
  3. どれぐらいの大きさの構造が必要か。(--repオプション)
  4. 密度はいくらか。(--densオプション)
  5. どんなファイル形式で出力するのか。(--formatオプション)

これらのうちいくつかは、省略することもできます。

さらに細かい設定もできます。

  1. 分極を0にするかどうか。(--nodepオプション)

--water モデル名

水分子モデルを指定します。デフォルトはtip3pです。

仮想相互作用サイトを持つ、tip4pモデルなども利用できます。

--rep x y z

単位胞をx,y,z軸方向に何回繰り返すかを指定します。デフォルトは1 1 1です。

--dens x

氷の密度を指定します。デフォルト値は結晶構造ごとにあらかじめ決まっています。

--format フォーマット名

出力ファイルの形式を指定します。デフォルトはgromacsです。

--nodep

全双極子モーメントを0にするプロセスをスキップします。イオンをドープした場合や、もともと4配位でない水分子を含む結晶構造の場合には、どのように水素結合ネットワークを再構成しても全双極子モーメントを0にできない場合がありますので、その場合にこのオプションを使用します。

Doped iceの生成

イオンをドープし、水分子を置換するかどうか。(--anion, --cationオプション)

今のところ、水素秩序氷にイオンをドープすることはできません。(dopeブランチでは可能)

ドープできるのは、単原子イオンのみ。

Doped hydrogen-ordered iceの生成

例えば、氷2にイオンをドープする場合には、まずgeniceで氷2の構造を作ります。

genice 2 > 2.gro

この構造に対して、analiceでドープします。

./analice.x 2.gro -c 0=NH4 -a 5=F > doped.gro

水分子0がカチオンNH4に、水分子5がアニオンFにそれぞれ置きかわります。(イオンはすべて単原子とみなされます)

途中、以下のようなメッセージが表示されます。

INFO   Anionize: {5: 'F'}.
INFO   Cationize: {0: 'NH4'}.
INFO   Number of inverted bonds by forced doping: 8

イオンを2個埋めこむために、8本の水素結合が反転されました。水素秩序氷にイオンをドープするためには、かならず一部の水素結合を反転して無秩序化する必要がありますが、その本数はイオンをドープする場所と個数に依存します。ここでは、反転する水素結合数が最小で済むような(かなり凝った)アルゴリズムを実装しています。 (水素無秩序氷にイオンをドープする場合はgeniceコマンドだけで作業することもできます。)

クラスレートハイドレートの生成

--guest

セミクラスレートハイドレートの生成

プラグイン

GenIceでは、さまざまな機能をプラグインで実現しています。プラグインにはformat, lattice, molecule, loaderの4種類があり、それぞれが適切なタイミングで呼び出されます。

ユーザーが独自に作成したプラグインを利用することも容易です。

Latticeプラグイン

結晶構造(単位胞の大きさ、形、分子の位置、結合)を定義するプラグインです。例えば、氷IIの構造を生成する場合には、次のようなコマンドを実行します。

% genice 2 > 2.gro

この時、GenIce内部では、lattices/2.pyプラグインが呼び出されて、水分子の配置を定義しています。GenIceには数百種類の氷の構造があらかじめlatticeプラグインとして準備されています。また、ユーザーが独自の氷の構造を定義することもできます。

Moleculeプラグイン

水分子のモデルを定義するプラグインです。例えば、水分子モデルとしてTIP4P (4サイトモデル)を採用したい場合には、以下のように-wオプションで分子モデルを指定します。

% genice CS2 -w tip4p > CS2.gro

moleculeプラグインは、クラスレートハイドレートのゲスト分子を指定する場合にも呼び出されます。また、ユーザーが独自の分子を定義することもできます。

Formatプラグイン

出力形式を定義するプラグインです。例えば、水分子の原子座標ではなく、重心位置とオイラー角を得たい場合には、次のように出力形式を指定します。

% genice 3 -f euler > 3.euler

また、分子の座標を出力する代わりに、解析結果を出力するプラグインもあります。_RDFプラグインは、すべての原子種の間の動径分布関数を出力します。

% genice 5 -f _RDF > 5.rdf

ユーザーが独自の出力フォーマットや解析手法をプラグインとして実装することも可能です。

formatプラグインは、結晶構造を生成する7つのステージの各段階で呼び出され、段階的に出力データを組みたてていきます。例えば、第2ステージまでで出力データが構成しおえた場合には、第3ステージ以降の処理は行われません。それぞれのプラグインが、そのステージでどんな処理を行うかを規定しています。

Loaderプラグイン

geniceがlatticeプラグインを使って結晶構造を新規に生成するのに対し、analiceは既存の座標ファイルを読みこんで、それに対して処理を行います。analiceがさまざまなフォーマットの座標ファイルを読み込むために利用するのがloaderプラグインです。

通常、analiceは読み込むファイルの拡張子を見て、自動的にプラグインの種類を選びます。例えば、拡張子が.groのファイルが指定された場合には、loaders/gro.pyプラグインが呼び出されます。-sオプションで、プラグインの種類を明示することも可能です。また、ユーザーが独自のloaderプラグインを実装することもできます。

拡張子 プラグイン名 ファイル形式
gro gro.py Gromacs
mdv mdv.py MDView (Angstrom)
mdva mdva.py MDView (Atomic unit)

ユーザー定義プラグイン

genice/analiceコマンドを実行するディレクトリ内に、formats/molecules/といったサブディレクトリを作り、その中に自作のformatプラグインやmoleculeプラグインを入れておくと、それらもGenIceから利用できます。

自作プラグインのプログラム方法は、プラグインの種類によってかなり違います。

検定

geniceはさまざまな出力フォーマットに対応するformatプラグインをそなえています。そのうち、以下のプラグインは、生成した氷の構造が正しくランダムになっているかどうかを検定する目的に使用できます。

_ringstatプラグイン

このプラグインは、geniceで生成したネットワーク構造に含まれる環を探し、その結合の向きの統計をとり、それが理想的な分布からどれぐらいずれているかを評価します。

水素結合には向きがあります。例えば、水素結合でできた四員環の辺にそって一周すすむと、その水素結合の向きは、順方向(F)か逆方向(B)の2通りで以下のような16種類の可能性があります。

これらには、始点をずらしたり環を逆に回ったりすると同一視できるものが含まれています。重複を省くと、実質的には以下の4種類の向きがありえます。

これらの出現頻度は、次のような近似で。モデル化できます。

  1. すべての節点は4結合である。
  2. すべての節点で入結合2、出結合2である。(ice rule)
  3. 注目する四員環以外には環がない。(ベーテ近似)

すると、これらの出現確率は、それぞれ

となります。

この近似と、geniceで生成したネットワーク内での環の向きの統計とを比較することで、ネットワーク構造が十分ランダム化されているかどうかを検定します。

以下の例では、氷VIの四員環の出現頻度を検定しています。

$ genice 6 -f _ringstat -r 4 4 4  |less

4 0 0000 16/41 0.39024 265/640 0.41406 
4 1 0001 16/41 0.39024 224/640 0.35000 
4 3 0011 8/41 0.19512 143/640 0.22344 
4 5 0101 1/41 0.02439 8/640 0.01250 
0.01205979926917721 4 dKL[4-ring]
...

第1カラムが環の大きさ、第3カラムが結合の向き(二進表現)、第2カラムはそれの10進表現、第4、5カラムがモデル分布(分数と小数)、第6、7カラムが生成した構造での分布となっています。また最後の行は、モデル分布と生成分布の間のKullback-Leibler Divergenceで、これが小さいほど2つの分布は近いと言えます。

_KGプラグイン

Kirkwood Gは次のように定義され、双極子の遠距離相関を計量します。

TIPS

Installation

Execution

実行時にUnicodeEncodeErrorが発生してファイルに書きこめない。

GenIceを実行しているTerminalのエンコーディング設定がUTF-8でない場合に起こります。環境変数PYTHONIOENCODINGUTF-8に設定して下さい。

% export PYTHONIOENCODING=UTF-8