Skip to content
Masakazu Matsumoto edited this page Dec 14, 2020 · 5 revisions

Cage information is missing in the crystal structure

The lattice plug-in, which produces the crystal structure of a typical classlet hydrate, such as CS1 (sI) or CS2 (sII), also contains information about the cage position, making it easy to position the guests, as Clathrate hydrates also notes.

In other cases, such as Zeolitic ices, guests cannot be placed because the cage is not located.

Use the cage plug-in to estimate the cage position based on the network structure.

pip install genice-cage

In this example, the genice-cif plug-in is used to generate ice of the zeolite SOD structure with genice, and the cage plug-in is used to locate the cage.

pip install genice-cif genice-cage
genice zeolite[SOD] -r 2 2 2 -f cage[python] > cages.py

We add `-r 2 2 2 'because if the unit cell is too small, one cage will span the unit cell and the cage will not be detected normally.

You will get the output.

cages="""
14         0.5000 0.5000 0.5000
14         0.7500 0.2500 0.7500
14         0.5000 0.0000 0.4286
14         0.7500 0.7500 0.2500
14         0.0000 0.0000 0.0000
14         0.7500 0.6786 0.7500
14         0.2500 0.7500 0.2500
14         0.0000 0.5000 0.9286
14         0.2500 0.6786 0.7500
14         0.5000 0.5000 0.9286
14         0.7500 0.2500 0.1786
14         0.0000 0.0000 0.4286
14         0.5000 0.0000 0.0000
14         0.0000 0.5000 0.5000
14         0.2500 0.2500 0.1786
14         0.2500 0.2500 0.7500
"""

The cell-relative coordinates of the center of 16 Kelvin 14-hedral cages are listed.

In addition, the unit cell of the 2x2x2 SOD structure is made into a plugin.

genice zeolite[SOD] -r 2 2 2 -f python > SOD.py

Combine these two to create a plug-in that generates a 2x2x2 SOD structure containing cage information. (Warning: this does not work with GenIce2.)

mkdir lattices
cat SOD.py cages.py > lattices/SOD222.py

The new plug-in is used to generate SOD structures containing methane in all 14-hedral cages.

genice SOD222 -g 14=me > SOD222+me.gro

In GenIce2, the --assess_cages option has been added for cases where cage information is not included in the crystal lattice data. In the example below, one of the dodecahedral cages of the DOH (sH) hydrate (there are two types) is filled with methane.

genice.x DOH --assess_cages -r 2 2 2 -g A12a=me > output.gro

Look up cages in Gromacs outputs

When Gromacs simulates the spontaneous formation of a clathrate hydrate from an aqueous solution, the trajectory may have a cage-like structure. You can also use `genice-cage 'to find it out.

A typical clathrate hydrate cage is a 10- to 20-hedral cage consisting of only 4-, 5-, or 6-membered rings. To detect this, use a command such as the following:

analice output.gro -H H -O O -w tip4p -f cage[10-20:ring=4-6]

-H and -O specifies the labels for hydrogen and oxygen in the Gromacs file.

Then you will get the output like following.

Cage 0: (0.3210629804603331, 0.5472966563158128, 0.296569238873769) 10 hedron
  Ring 3112: (0.26907623255152646, 0.5402898613866303, 0.28480659081662124) 6 gon
    Nodes: [3446, 6827, 6821, 6760, 2754, 1726]
  Ring 3113: (0.30123819073913627, 0.5320125746829685, 0.259455021335314) 6 gon
    Nodes: [3446, 6827, 6821, 999, 4561, 1681]
  Ring 3115: (0.2952168818936886, 0.5672722935117082, 0.33551997756963636) 6 gon
    Nodes: [2754, 1726, 5151, 6150, 5701, 7444]
  Ring 3116: (0.31356774656125147, 0.5816953265156833, 0.2869234915193971) 6 gon
    Nodes: [3446, 1726, 5151, 1410, 7420, 1681]
  Ring 3448: (0.2906858945665161, 0.511557985043145, 0.3025484439403703) 6 gon
    Nodes: [2754, 6760, 6821, 999, 969, 7444]
...

It shows the center of mass position (cell relative) and number of faces for the cage, the center of mass position and number of vertices (water molecule) for each face, and labels for the vertices that make up the face. If you find it a bit verbose for readability, choose the JSON format. The JSON format is easy to read in a variety of programming languages.

analice output.gro -H H -O O -w tip4p -f cage[12-20:ring=4-6:JSON]

If the input file contains multiple frames, you should specify the template of the output file name with -o option.

JSON format contains the following information.

  • rings: A list of lists of water labels that make up each ring.
  • cages: A list of lists of labels for the rings that make up each cage.
  • ringpos, cagepos: A list of fractional positions of the centers of mass of rings and cages.

Here is a sample python code to show the contents.

import numpy as np
import json
with open("cages00000.json") as f:
    cages = json.load(f)

# the dimension of the cell; should be read from the gro file.
cell = np.array([10.01000, 2.51000, 2.51000])

# convert the vector to a diagonal matrix/
cellmat = np.diag(cell)

# fractional to absolute coordinate
cagepos = np.array(cages["cagepos"]) @ cellmat
print("Centers of the cages")
print(cagepos)

# rings is a list of the labels of molecules.
# Only the water molecules are labelled and the number starts from zero.
rings = cages["rings"]

for cage in cages["cages"]:
    # count the number of cycles (faces)
    print(f"{len(cage)}-hedra")
    for ringid in cage:
        vertices = rings[ringid]
        print(f"{len(vertices)}-gon", vertices)

The "polyhedron" that the cage plug-in can detect satisfies the following conditions and, unlike geometric polyhedra, allows curved faces. It does not detect cavities whose Euler characteristic is not 2, such as channels or pores.[MBO2007]

  1. A vertex is shared by two or three faces.
  2. Two faces share an edge.
  3. The Euler characteristic F-E+V is 2, where F, E, and V are the number of faces, edges, and vertices constituting the polyhedron.

References

[MBO2007] M. Matsumoto, A. Baba, and I. Ohmine, Topological building blocks of hydrogen bond network in water, J. Chem. Phys. 127, 134504 (2007). http://doi.org/10.1063/1.2772627


結晶構造にケージの情報が含まれていない場合

CS1 (sI)やCS2 (sII)のような典型的なクラスレートハイドレートの結晶構造を生成するlatticeプラグインには、ケージの位置に関する情報も含まれていますので、Clathrate hydratesにも書いたように、ゲストを簡単に配置することができます。

一方、Zeolite DBから入手した氷などの場合には、ケージの位置が特定されていないので、ゲストを配置できません。

ネットワークの構造をもとに、ケージの位置を推定するには、cageプラグインを使用します。

pip install genice-cage

ここでは、genice-cifプラグインを使ってゼオライトSODの構造の氷をgeniceで生成し、cageプラグインでそのケージの位置を特定します。

pip install genice-cif genice-cage
genice zeolite[SOD] -r 2 2 2 -f cage[python] > cages.py

ここで-r 2 2 2を加えたのは、単位胞が小さすぎると1つのケージが単位胞をまたいでしまい、正常にケージを検出できないためです。

次のような出力が得られます。

cages="""
14         0.5000 0.5000 0.5000
14         0.7500 0.2500 0.7500
14         0.5000 0.0000 0.4286
14         0.7500 0.7500 0.2500
14         0.0000 0.0000 0.0000
14         0.7500 0.6786 0.7500
14         0.2500 0.7500 0.2500
14         0.0000 0.5000 0.9286
14         0.2500 0.6786 0.7500
14         0.5000 0.5000 0.9286
14         0.7500 0.2500 0.1786
14         0.0000 0.0000 0.4286
14         0.5000 0.0000 0.0000
14         0.0000 0.5000 0.5000
14         0.2500 0.2500 0.1786
14         0.2500 0.2500 0.7500
"""

Kelvin 14面体ケージ16個の中心のセル内相対座標が列挙されました。

また、SOD構造の2x2x2倍構造の単位胞をplugin化します。

genice zeolite[SOD] -r 2 2 2 -f python > SOD.py

これら2つを連結して、ケージ情報を含む、SOD構造の2x2x2サイズの氷を生成するプラグインを作ります。(この方法はGenIce2では使えません)

mkdir lattices
cat SOD.py cages.py > lattices/SOD222.py

新しいプラグインを使い、すべての14面体ケージにメタンが入ったSOD構造を生成します。

genice SOD222 -g 14=me > SOD222+me.gro

GenIce2では、結晶格子のデータにケージの情報が含まれていない場合のために、--assess_cagesオプションが追加されました。下の例では、DOH (sH) ハイドレートの12面体ケージ(2種類ある)の一方をメタンで埋めます。

genice.x DOH --assess_cages -r 2 2 2 -g A12a=me > output.gro

Gromacsの出力の中にケージを探す

クラスレートハイドレートが水溶液から自発的に生成するような条件をGromacsでシミュレーションすると、トラジェクトリの中に、ケージのような構造ができているかもしれません。これを見付けだすのにも、genice-cageが利用できます。

典型的なクラスレートハイドレートのケージは10〜20面体で、4、5、6員環のみからできています。これを検出するには、以下のようなコマンドを使います。

analice output.gro -H H -O O -w tip4p -f cage[10-20:ring=4-6]

-H-OはGromacsファイル内での水素と酸素のラベルです。

すると、以下のような出力が得られます。

Cage 0: (0.3210629804603331, 0.5472966563158128, 0.296569238873769) 10 hedron
  Ring 3112: (0.26907623255152646, 0.5402898613866303, 0.28480659081662124) 6 gon
    Nodes: [3446, 6827, 6821, 6760, 2754, 1726]
  Ring 3113: (0.30123819073913627, 0.5320125746829685, 0.259455021335314) 6 gon
    Nodes: [3446, 6827, 6821, 999, 4561, 1681]
  Ring 3115: (0.2952168818936886, 0.5672722935117082, 0.33551997756963636) 6 gon
    Nodes: [2754, 1726, 5151, 6150, 5701, 7444]
  Ring 3116: (0.31356774656125147, 0.5816953265156833, 0.2869234915193971) 6 gon
    Nodes: [3446, 1726, 5151, 1410, 7420, 1681]
  Ring 3448: (0.2906858945665161, 0.511557985043145, 0.3025484439403703) 6 gon
    Nodes: [2754, 6760, 6821, 999, 969, 7444]
...

ケージの重心位置(セル内相対)と面の数、それぞれの面の重心位置と頂点(水分子)の数、面を構成する頂点のラベルが表示されています。可読性のために少し冗長だと感じる場合は、JSON形式を選びましょう。JSON形式なら、さまざまなプログラム言語で簡単に読みこめます。

analice output.gro -H H -O O -w tip4p -f cage[12-20:ring=4-6:JSON]

入力ファイルに複数のフレームが含まれている場合、出力ファイル名のテンプレートを -o オプションで指定する必要があります。

ここでは、内容を表示するためのpythonのサンプルコードを示します。

JSON形式の内容は次の通りです。

  • rings: それぞれの環を構成する水分子のラベルのリストのリスト
  • cages: それぞれのケージを構成する環のラベルのリストのリスト
  • ringpos, cagepos: 環とケージの重心位置(セル相対)のリスト
import numpy as np
import json
with open("cages00000.json") as f:
    cages = json.load(f)

# the dimension of the cell; should be read from the gro file.
cell = np.array([10.01000, 2.51000, 2.51000])

# convert the vector to a diagonal matrix/
cellmat = np.diag(cell)

# fractional to absolute coordinate
cagepos = np.array(cages["cagepos"]) @ cellmat
print("Centers of the cages")
print(cagepos)

# rings is a list of the labels of molecules.
# Only the water molecules are labelled and the number starts from zero.
rings = cages["rings"]

for cage in cages["cages"]:
    # count the number of cycles (faces)
    print(f"{len(cage)}-hedra")
    for ringid in cage:
        vertices = rings[ringid]
        print(f"{len(vertices)}-gon", vertices)

なお、cageプラグインが検出できる「多面体」とは以下の条件をみたすもので、幾何学的な多面体とは異なり、曲がった面を許容します。また、チャネルや細孔などの、Euler標数が2でない形状の空洞は検出しません。[MBO2007]

  1. 1つの頂点を2つまたは3つの面が共有する。
  2. 1つの辺を2つの面が共有する。
  3. Euler標数F-E+Vが2である。ただし、F、E、Vは多面体を構成する面、辺、頂点の数である。

References

[MBO2007] M. Matsumoto, A. Baba, and I. Ohmine, Topological building blocks of hydrogen bond network in water, J. Chem. Phys. 127, 134504 (2007). http://doi.org/10.1063/1.2772627