cif形式の結晶構造をpythonで描画してみたよ

pythonでcifの描画

化学、物理で結晶の研究をしているとcif形式の結晶構造データを扱うことが多いです。今回はこれをpythonに読み込んで描画してみようって試みです。

インプットとアウトプット

・インプット
-cifファイル
有機結晶の例としてアスピリンの結晶構造¹を使います。Crystallography Open Database(COD)²から無料でダウンロードして使います。ファイル名は1515581.cifでした。

・アウトプット
こんな感じのpng画像が得られます。結晶軸も一緒に描画しています。赤線がa軸、緑線がb軸、青線がc軸です。

準備

pythonプログラムはAnaconda³のspyderで作成してます。プログラムの実行にあたって次のパッケージのインストールが必要です。

1)pymatgen⁴
2)plotly⁵

pymatgenは物理・化学系の解析に特化したpythonライブラリです。cifの読み込みに使います。
plotlyは描画用のライブラリです。pythonの描画といえばmatplotlib一択ですが、結晶構造の単位格子に含まれる全原子を三次元空間にプロットするような場合、結晶がちょっと大きくなるとmatplotlibでは動作が重くなってしまうきらいがあります。三次元空間で描画した結晶構造をぬるぬる回したいので描画をplotlyで行います。

実行

from pymatgen.core.structure import Structure
import numpy as np
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default='browser'

##1) inport cif file
aspirin = Structure.from_file("./1515581.cif")
                        
aspirin_lattice = aspirin.lattice.matrix
aspirin_atom_num = np.size(aspirin.atomic_numbers)
aspirin_atom_posi = aspirin.cart_coords
aspirin_atom_kind = aspirin.atomic_numbers

a_axis = [[0.0, aspirin_lattice[0][0]], [0.0, aspirin_lattice[0][1]], [0.0, aspirin_lattice[0][2]]]
b_axis = [[0.0, aspirin_lattice[1][0]], [0.0, aspirin_lattice[1][1]], [0.0, aspirin_lattice[1][2]]]
c_axis = [[0.0, aspirin_lattice[2][0]], [0.0, aspirin_lattice[2][1]], [0.0, aspirin_lattice[2][2]]]

##2) make bond list

#covalent radius[AA]
#reference... https://www.hulinks.co.jp/support/c-maker/qa_005.html
covalent_radius_list = [0.00, 0.37, 0.32,
                        1.34, 0.90, 0.82, 0.77, 0.75, 0.73, 0.71, 0.69]

dstance = []
bond_list = []
bond_length = []

for j in range(aspirin_atom_num):
    j_cov_rad = covalent_radius_list[aspirin_atom_kind[j]]
    
    for k in range(aspirin_atom_num):
        if j < k:
            distance = np.linalg.norm(aspirin_atom_posi[j] - aspirin_atom_posi[k])
            
            k_cov_rad = covalent_radius_list[aspirin_atom_kind[k]]
                
            if distance <= j_cov_rad + k_cov_rad:
                bond_list.append([j, k])
                bond_length.append([distance, j_cov_rad + k_cov_rad])
            else:
                pass
            
        else:
            pass

##3) 3D plot by plotly

colorscale = ["powderblue", "black", "black", "black", "green", "lightgrey", "blue", "red"]

trace = go.Scatter3d(
   x = aspirin_atom_posi[:, 0],
   y = aspirin_atom_posi[:, 1],
   z = aspirin_atom_posi[:, 2],
   name = "atoms",
   mode = 'markers', 
   marker = dict(
       size = 7,
       color = aspirin_atom_kind,
       colorscale = ["powderblue", "black", "black", "black", "green", "lightgrey", "blue", "red"]
       )
   )

a_axisp = go.Scatter3d(
   x = a_axis[0], y = a_axis[1], z = a_axis[2],
   name = 'a_axis',
   mode = 'lines',
    marker=dict(size=0, color='red', opacity=0.8)
    )

b_axisp = go.Scatter3d(
   x = b_axis[0], y = b_axis[1], z = b_axis[2],
   name = 'b_axis',
   mode = 'lines',
    marker=dict(size=0, color='green', opacity=0.8)
   )

c_axisp = go.Scatter3d(
   x = c_axis[0], y = c_axis[1], z = c_axis[2],
   name = 'c_axis',
   mode = 'lines',
    marker=dict(size=0, color='blue', opacity=0.8)
  )

axis = [a_axisp] + [b_axisp] + [c_axisp]
bond = []

for l in range(len(bond_list)):
    bond_x = [aspirin_atom_posi[bond_list[l][0]][0], aspirin_atom_posi[bond_list[l][1]][0]]
    bond_y = [aspirin_atom_posi[bond_list[l][0]][1], aspirin_atom_posi[bond_list[l][1]][1]]
    bond_z = [aspirin_atom_posi[bond_list[l][0]][2], aspirin_atom_posi[bond_list[l][1]][2]]
    
    bond_i = go.Scatter3d(
        x = bond_x, y = bond_y, z = bond_z,
        name = 'bond',
        mode = 'lines',
        marker=dict(size=5,line=dict(width=15,), color='gray', opacity=0.8)
        )
    
    bond = bond + [bond_i]
    
fig = go.Figure(data = bond + [trace] + axis)

noaxis=dict(
    showbackground=False,
    showgrid=False,
    showline=False,
    showticklabels=False,
    ticks='',
    title='',        
    zeroline=False)

fig.update_layout( 
    width=1500, height=500,
    margin=dict(l=0, r=0, b=0, t=0),
    scene_aspectmode='data',
    scene = dict(
        xaxis = noaxis,
        yaxis = noaxis,
        zaxis = noaxis),
    plot_bgcolor="gray",
    )

fig.show()

少し解説します。

1)inport cif file
pymatgenを用いて1515581.cifを読み込みます。結晶軸も描画するのでこれらの情報もpymatgenを用いて読み込みます

2)make bond list
共有結合の描画のためにinport cif fileで読み込んだ原子位置から共有結合のリスト(bond_list)を作ります。
全部の原子間距離を計算して、原子間距離が二つの原子の共有結合半径の和より小さければ結合しているとみなします。各原子の共有結合半径はcovalent_radius_listに手動で打ち込んでいます。

3)3D plot with plotly
plotlyで描画します。spyderからだとplotlyがデフォルトでは動作しません。そこで、次のように指定することでplotがブラウザ上で動作させます。

pio.renderers.default='browser'

結果として次のような画像がブラウザ上で表示されます。

ブラウザ上では拡大縮小、ぐりぐり回す等の動作がスムーズに動きます。

まとめ

やってみてPythonの勉強には非常になりました。しかし現状、cifを無料で表示できるソフトがVESTA⁶、Mercury⁷、CrystalExplorer⁸など多種ある中でそれらより綺麗に表示するのは難しいなあと思った次第です。結晶構造から自力でプログラム作って解析したいよって人には役に立つかもです。

1. COD ID：1515581
   アスピリンの結晶構造の文献。
   Varughese, Sunil, et al. Chemical Science 2.11 (2011): 2236-2242. ↩

2. Crystallography Open Database(COD)
   http://www.crystallography.net/cod/ ↩

3. Anaconda公式HP、ダウンロードもここから
   https://www.anaconda.com/products/distribution ↩

4. pymatgenのcondaによるインストール
   https://anaconda.org/conda-forge/pymatgen ↩

5. plotlyはターミナルからcondaでインストールしようとしてもできないみないです。インストールはこのブログを参照。
   http://python.geo.jp/archives/51 ↩

6. VESTA公式HP
   https://jp-minerals.org/vesta/jp/ ↩

7. Mercury公式HP
   https://www.ccdc.cam.ac.uk/solutions/csd-core/components/mercury/ ↩

8. crystalexplorer公式HP
   https://crystalexplorer.net/ ↩