目的
結晶構造可視化ソフトウェアであるVESTA をPython から外部操作することにより、良い感じのGif 動画を作成したいと思います。苦労して解析した構造や、第一原理計算等で得られた美しいVolumetric Dataを動画として提示して、プレゼンテーションでドヤ顔しましょう。
動作環境
Windows10, Anaconda, Python3.6.10 + JupyterLab0.35.4
VESTA v3.5.0
VESTA の概要と外部プログラムからの操作
公式サイト：https://jp-minerals.org/vesta/jp/
VESTA は結晶構造の可視化ができるソフトウェアで、単なる構造だけでなく、電子密度・核密度分布等の三次元データや結晶外形も可視化できます。かなりの多機能ですが、なんとフリーソフトです。Windows, Mac, Linux で軽快に動作します。
~~最近、外部のプログラムから操作するためのインターフェースが備わったようで、開発者の一人である泉先生がある日Twitter につぶやいた共有Evernote で解説されています。~~Evernote の公開はやめたようですね。何か不都合があったのかもしれません。公式ドキュメントには記載がないので、まだ試験的な機能ということなのかもしれません。
外部操作は、コマンドラインで以下のオプションを利用することで実現できます。Evernoteのノート「VESTA によるリアルタイム構造・等値曲面表示」https://t.co/4qPkUZHLPP
— 泉 富士夫 (@Izumi_Fujio) February 26, 2020
では、最高クラスに役立つ機能を書き漏らしていたことが判明した。東北大講習会の無期延期で意気消沈しているので、今週はそれで遊ぶ。
- -open hoge ：ファイルhogeを開く。
- -reopen hoge ：ファイルhogeを開く。
- -close ：現在表示されているファイルを閉じる。
- -rotate_x phi：x軸回りにphi°回転させて構造を表示する。 (phi は角度、単位は°)
- -rotate_y phi：y軸回りにphi°回転させて構造を表示する。(phi は角度、単位は°)
- -rotate_z phi：z軸回りにphi°回転させて構造を表示する。 (phi は角度、単位は°)
-rotate_x,y,z は単独では動かず、-openと組み合わせて使います。また、-reopenとの組み合わせでも動作すべきな気がしますが、なぜか動きません。（Windows10+VESTA v3.5.0の場合で、他の環境では試していません。）
これらの機能をPython から呼び出して使うことが、今回の記事の目的の一つです。(pymatgen 等と連携しても便利だと思いませんか？)
構造の作成
なにはともあれ、まずは表示するための構造ファイルが必要です。ここでは、Si の結晶構造を作ってみます。
VESTA の[File]→[New Structure]を開き、以下のように設定して構造を作成します。
また、必須ではないですが、[Edit]→[Bonds]で以下のようにして結合を描いておきましょう。見栄えがよくなります。
ここまでうまく出来ていれば、以下のような構造ができるはずです。
最後に、[File]→[Save]からファイル名Si.vesta として保存しておきます。
Python からVESTA を操作
Python には外部のプログラムを動作させるための標準モジュールsubprocess が用意されていますので、これを利用します。
さて、subprocess を使うには、VESTA の場所を指定しなければなりません。Windows の場合は環境変数PathにVESTA のインストールディレクトリを設定すれば良いのですが、VESTA 自体は特別なインストールを必要としないですし、会社のPC等で管理者権限を持っていない場合には面倒なので、コード内にベタ書きすることにします。
VESTA = r"C:\Program Files\VESTA-win64\VESTA.exe"
raw 文字列部分には各自のPC におけるVESTA の実行パス(置いてある場所+VESTA.exe)を指定してください。
subprocess で外部プログラムを実行する方法はいくつかありますが、色々試した限り非同期で実行できるPopen を使うと良さそうです。Popen では、第一引数にコマンドを表すリストを指定します。このリストには入力ファイルやオプション等を要素で分けて入れていきますが、実行時には文字列として結合されるようです（細かいところは調べていないので、厳密には違うかもしれませんが)。単にVESTA を起動するだけならば、VESTA の実行パスをリストに入れるだけでOK です。
import subprocess
subprocess.Popen([VESTA])
VESTA が開きましたでしょうか。次に、先程作成した構造ファイルSi.vesta を開いてみましょう。
Pythonの実行ディレクトリにSi.vesta ファイルをおいて、以下のようなコードを実行します。
subprocess.Popen([VESTA, '-open', 'Si.vesta'])
先程起動していたVESTA のウィンドウにSi.vesta の構造が描かれるはずです。次に、この構造を一旦閉じてみます。
subprocess.Popen([VESTA, '-close'])
Si.vesta が描かれていたタブが消えます。このように、VESTA のウィンドウが多重起動することなく、Pythonからの指示で構造を描いたり消したりすることができるわけです。
また、前述のように、-rotate_x, y, z を使うとSi.vesta を開く際に構造を回転して表示できます。ここでは、x, y, z にそれぞれ10°ずつ傾けた状態で表示してみましょう。subprocess.Popen で複数のオプションを設定したい場合には、単にリストの要素を増やしていけば良いのですが、スペースで区切って数値を入力するようなオプションがある場合は、オプションと数値でリストの要素を分けておく必要があるという点に注意が必要です。たとえば、今回の場合には以下のようにします。
subprocess.Popen([VESTA, '-open', 'Si.vesta', '-rotate_x', '10', '-rotate_y', '10', '-rotate_z', '10'])
このコードは、コマンドプロンプトで以下のコマンドが実行されることに対応します。
C:\Program Files\VESTA-win64\VESTA.exe -open Si.vesta -rotate_x 10 -rotate_y 10 -rotate_z 10
問題なく動作していれば、以下のような回転した画像が表示されるはずです。
動画の作成
さて、本題の動画の作成に移りましょう。現状のVESTA の外部操作インターフェースには画像の保存機能がないため、これ単独では動画の作成はできません。そこで、ここではpyautogui(pip, conda等でインストールして下さい) のスクリーンショット機能を活用することで、Python 側で画像を取得することで動画を作成したいと思います。
大まかな流れは以下のとおりです。
- VESTA を起動しておく。（普通に起動しても良いし、Python から起動しても良い）
- VESTA の適当なボタンの画像を用意し、pyautogui でウィンドウの画面上の座標を取得する。
- スクリーンショットを取り、結果を見ながら範囲を調整する。
- 範囲が決まったら、動画にしたい構造を順に表示しながら、スクリーンショットを取得する。
- スクリーンショットをまとめてPIL の機能でGif 動画として出力する。
ここでは、スクリーンショットの確認等で画像 をinline で表示できると便利なので、jupyter notebook/lab の利用を推奨します。
それでは、Siの構造を回転させていく動画を作ってみます。まず、必要なものをimport します。
import subprocess
import time
import pyautogui
PIL はpyautogui が内部的に使っているだけなので、明示的にimport する必要はありません。次に、なんでも良いのですがVESTAのウインドウからボタン等の部分だけのスクリーンショットをWindows の機能等で取ります。(最近だと[Win]+[Shift]+[S]のショートカットが便利です)今回は、以下のような画像にしました。（結晶構造を見る方向の選択ボタンです。）
ここから先は、VESTA のウインドウは動かさないようにしましょう。また、他のウインドウがかぶっているとややこしいことになるので、それも避けましょう。
pyautogui を用いてVESTA のウィンドウ位置、というか先に取得したボタンのスクリーンショット画像が、実際の画面上のどの座標にあるかを判別します。
loc = pyautogui.locateOnScreen('VESTA_button.png')
loc
Box(left=836, top=207, width=185, height=32)
画像のleft(左), top(上), width(幅), height(高さ)が取得できました。ここで、試しにこの座標を基準に、pyautogui を用いてスクリーンショットを取ってみましょう。region に所望の(left, top, width, height)を指定します。
left = loc.left
top = loc.top
size = (700, 700)
sc = pyautogui.screenshot(region=(left+200, top, size[0], size[1]))
sc
pyautogui.screenshot の実行結果はPIL.Image クラスのインスタンスになっているため、Jupyter Notebook/Lab では上記のように評価するだけで画像が表示されます。(もちろん、matplotlib のpyplot.imshow を使っても良いです)
ここでは以下のような画像が得られましたが、当然VESTA のウィンドウサイズにもよりますので、適宜調整してください。後々PIL の機能でcrop することもできるので、大きめに取っておいても良いかもしれません。
さて、ここまで来たら、あとは構造を次々に表示して、スクリーンショットを取得していくだけです。x軸周りに0～180° まで1° 刻みで回転させてみます。空のリストimages を作っておいて、スクリーンショットで得られた結果をappend していきます。また、time.sleep （単位は秒で指定します。）を使って適当な待ち時間を入れています。待ち時間の値は実行するPC の性能によっても変わってくると思うので、試しながら実行してみてください。また、待ち時間を試す際にはloopの数を少なめ(2～3)にして試すことをおすすめします。
# 結構時間がかかるので注意！
images = []
for i in range(180):
subprocess.Popen([VESTA, "-open", "Si.vesta", "-rotate_x", str(i)])
time.sleep(0.5)
sc = pyautogui.screenshot(region=(left+200, top, size[0], size[1]))
images.append(sc)
subprocess.Popen([VESTA, "-close"])
time.sleep(0.5)
大量に構造が開いてしまうとメモリを食うので、スクリーンショットを撮ったあとは一々構造を閉じています。
さて、これで動画の元になる画像群がimages リストとして出来上がりました。そのままGif動画にしても良いのですが、見栄えを良くするため結晶構造の表示部分のみをcrop して、またファイルサイズが大きくなるので画像サイズを縮小してから動画にしてみます。
new_images = []
for i in images:
tmp_image = i.crop((60, 60, 680, 600))
size1, size2 = tmp_image.size
new_images.append(tmp_image.resize((int(size1/2), int(size2/2))))
new_images[0].save("Si_rotate_x.gif", save_all=True, append_images=new_images[1:], duration=20, loop=0)
以下のようなGif動画が得られるはずです。
回転以外の操作について
VESTA の外部プログラムとのインターフェースは未だ発展途上なためか、回転以外のオプションは存在しません。
しかし、表示している.vesta ファイル自体を編集してしまうことで、だいたい何でも対応が可能です。.vesta ファイルはテキストファイルなので、編集するのに苦労はありません。
さて、しかし、.vesta ファイルのフォーマットに関する情報は私の知る限り、どこにもまとまっていません。なので、作りたい動画によって試行錯誤が必要になるでしょう。ここでは、試しに作ってみた動画の例を3つ紹介します。
ズームイン・ズームアウト
最も簡単な例として、ズームイン/ズームアウトを行う動画を作成します。構造は先程のSi.vesta を引き続き使っていきます。
まずは、ズームイン/ズームアウトに対応する.vesta ファイルの書式を特定しましょう。これは実は簡単で、ズームイン・ズームアウトの条件だけを変えた.vesta ファイルを作り、それらのdiff を取れば良いのです。各自試して頂くと良いかと思いますが、ここで結果から書きますと.vesta ファイル内のSCENEというブロックの最後の行がズームに対応します。
SCENE
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
0.000000 0.000000 0.000000 1.000000
0.000 0.000
0.000
1.000
これを書き換えたファイルを逐次作成しながら表示すれば、ズームを変えた動画が作れるわけです。先程の回転の例と組み合わせて「回転させながらズームを変える」動画を作ってみました。以下のコードで作成できます。
images = []
rotate_angle = np.arange(0, 180, 1)
zoom = np.linspace(0.5, 1.5, len(rotate_angle))
for a, zm in zip(rotate_angle, zoom):
with open("Si.vesta") as f:
lines = f.readlines()
idx_SCENE_kw = [i for (i, l) in enumerate(lines) if "SCENE" in l][0]
lines[idx_SCENE_kw+7] = " {zoom}".format(zoom=str(zm))
with open("tmp.vesta", "w")as f:
f.writelines(lines)
subprocess.Popen([VESTA, "-open", "tmp.vesta", "-rotate_x", str(a), "-rotate_y", str(a), "-rotate_z", str(a)])
time.sleep(0.5)
sc = pyautogui.screenshot(region=(left, top, size[0], size[1]))
images.append(sc)
subprocess.Popen([VESTA, "-close"])
time.sleep(0.5)
new_images = []
for i in images:
tmp_image = i.crop((60, 120, 660, 660)) # この辺は適宜変えてください。
size1, size2 = tmp_image.size
new_images.append(tmp_image.resize((int(size1/2), int(size2/2))))
new_images = new_images + new_images[::-1]
new_images[0].save("Si_rotate_x_zoom.gif", save_all=True, append_images=new_images[1::2], duration=20, loop=0)
3Dデータの等値面(Isosurface)
これは泉先生のEvernote にもありますが、Isosurface の値を連続的に変えた動画が作れます。先程と同様の方法で確認すると、ISURFというブロックにある値を変えれば良いようです。
ISURF
1 1 0.04 255 255 0 127 255
0 0 0 0
この例では0.04 という値になっている部分を書き換えればOK です。
以下のコードで動画を生成できます。ただし、当然ですがIsosurface の元になるデータは別途用意してください。ここでは、PWscf で適当に生成した電子密度分布（価電子のみ）を用いています。また、.vesta ファイルには3D データは含まれない点に注意してください。なので、3D データを同じディレクトリに置いておく必要がありますこれを忘れると、Isosurface が描かれませんのでご注意ください。
3Dデータを用意し、そのデータを用いてIsosurface を描いた状態を.vesta ファイルで保存(Si_edens.vesta)した上で、以下のコードを実行してください。
images = []
rotate_angle = np.arange(0, 180, 4)
isurf_ary = np.linspace(0.01, 0.08, len(rotate_angle))
for a, s in zip(rotate_angle, isurf_ary):
with open("Si_edens.vesta") as f:
lines = f.readlines()
idx_ISURF_kw = [i for (i, l) in enumerate(lines) if "ISURF" in l][0]
#以下は、元にする.vesta ファイルによって若干変わる可能性があります。
lines[idx_ISURF_kw+1] = " 1 1 {s} 255 255 0 127 255\n".format(s=str(s))
with open("tmp.vesta", "w")as f:
f.writelines(lines)
subprocess.Popen([VESTA, "-open", "tmp.vesta", "-rotate_x", str(a), "-rotate_y", str(a), "rotate_z", str(a)])
time.sleep(0.5)
sc = pyautogui.screenshot(region=(left, top, size[0], size[1]))
images.append(sc)
subprocess.Popen([VESTA, "-close"])
time.sleep(0.5)
new_images = []
for i in images:
tmp_image = i.crop((60, 120, 660, 660))# この辺は適宜変えてください。
size1, size2 = tmp_image.size
new_images.append(tmp_image.resize((int(size1/1), int(size2/1))))
new_images = new_images + new_images[::-1]
new_images[0].save("Si_eledens.gif", save_all=True, append_images=new_images[1::], duration=60, loop=0)
電子密度のデータは適当に計算しただけで単位も定かではない(PWScf は普段使わないので...)ので、値はIsosurface の見栄えが良くなるように適当に設定しています。まずはVESTA 上で値の範囲を確認してから、動画の作成をする方が良いと思います。
格子面の連続的な描画
最後に、3D データを読み込んだ状態で格子面を描き、その位置を連続的に変えて断面像を可視化してみます。
格子面の情報はSPLANというブロックに書かれているようです。また、格子面の位置は絶対値指定(単位：Å)で書かれています。以下の例では、1.000000E+000 0.000000E+000 0.000000E+000 が面指数(100)、2.7155 が格子面の位置に相当します。
SPLAN
1 1.000000E+000 0.000000E+000 0.000000E+000 2.7155 255 0 255 192
0 0 0 0
さきほどのIsosurface の例と同様に、あらかじめ3Dデータを読み込み、Isosurface が表示された場合は消し、さらに(100)面を描いておきます。その状態で.vesta ファイルを保存(Si_edense2.vesta)した上で、以下のコードを実行します。
d_ary = np.arange(0.0, 5.43, 0.1)
images = []
for d in d_ary:
with open("Si_edens2.vesta") as f:
lines = f.readlines()
idx_ISURF_kw = [i for (i, l) in enumerate(lines) if "SPLAN" in l][0]
lines[idx_ISURF_kw+1] = " 1 1.000000E+000 0.000000E+000 0.000000E+000 {d} 255 0 255 192\n".format(d=str(d))
with open("test2.vesta", "w")as f:
f.writelines(lines)
subprocess.Popen([VESTA, "-open", "test2.vesta"])
time.sleep(0.5)
sc = pyautogui.screenshot(region=(left, top, size[0], size[1]))
images.append(sc)
subprocess.Popen([VESTA, "-close"])
time.sleep(0.5)
new_images = []
for i in images:
tmp_image = i.crop((60, 60, 680, 600))
size1, size2 = tmp_image.size
new_images.append(tmp_image.resize((int(size1/1), int(size2/1))))
new_images = new_images + new_images[::-1]
new_images[0].save("Si_eledens_plane.gif", save_all=True, append_images=new_images[1:], duration=80, loop=0)
まとめ
VESTA の外部プログラム用のインターフェースをPython から操作することで、結晶構造や3Dデータを含む動画の作成を行いました。なかなか良い感じの動画が出来たんじゃないかと思います。
pyautogui を用いるとスクリーンショットの取得からGif 動画の作成までPython だけで完結できるため、すっきりとした構成になります。今回はGif動画にしましたが、たとえばOpenCV を用いればMP4フォーマットの動画にすることもできると思います。そちらの方がプレゼンテーション用途だと良いかもしれません。(OpenCV はなにかAnaconda と相性が悪いのかインストールがうまくいかなかったトラウマがあるので試していませんが...)
課題というか、できるかと思ってできなかった点として、subprocess.Popen でVESTA を実行する際にPython のファイルオブジェクトをそのまま渡せるかと思って色々試したのですが、無理でした。これができれば一々テンポラリファイルを作成しなくても良いのに...と思いましたが、冷静に考えれば実行時間のボトルネックにはなってないので、出来たとしてもあまり良いことは無いかもしれません。
あとは、.vesta ファイルのフォーマットを一々調べるのが面倒なので、情報を誰かまとめてくれると助かります。（それとも、どこかに公開されているんでしょうか？）専門家諸氏、どうぞよろしくお願いいたします。