分子の座標から GIF animation を作る

わしの本職は高分子電解質のシミュレーション屋でして、右のような簡単なアニメーションを作ることができると、同業者にもアピールできますし、研究が捗ります。ここにいろいろなサンプルが展示してあります。んまぁ、分子シミュレーションの中でも低分子や非電解質、全原子モデルの計算をしていると、いろいろ描画ツールも揃っていたりしてありがたいのですが、わしらはもう少しマクロな系の計算をしているので、描画にもいろいろ苦労が多いです。でも、右のようなアニメーション程度でしたら Linux BOX に入っているオーソドックスなソフトを組み合わせれば簡単に作ることができます。

数値計算は自分でやろう

最初に、アニメーションの基礎になる分子の配置のデータを用意します。このアニメーションでは、
・DNA上の電荷(りん酸基):64塩基対のB型DNAなので64個の座標が二組
・Na⁺イオン:1152個の座標
・Cl_-イオン:1024個の座標
が、それぞれこんな感じに

  -0.260714E+01  -0.185021E+02   0.529618E+02
  -0.182008E+02   0.226066E+02   0.529646E+02
  -0.222052E+02   0.188900E+02   0.138810E+02
   0.439628E+02   0.119847E+02   0.167322E+02
   ............   ............   ............

入っているファイルをそれぞれ用意します。因みに単位は nm です。水和したイオンは剛体球として粗視化しているので、分子内振動などは考えていません。てことで一つの分子について (x, y, z) の 3つの座標があるということです。シミュレーションセルの大きさは半径 61 nm の球としています。

アニメーションでは、たとえば 20 シミュレーションステップの配置を使うと、上のようなファイルのセットが20組あることになります。

awk で分子座標のファイルを整形

上のような生座標のファイルは大型計算機で計算したりした結果ですが、その解析から先は Linux BOX でも充分可能です。たとえば、上のNa⁺イオンの座標のファイルの中から「原点を中心とする一辺の長さが 15 nm の立方体」に含まれるイオンの座標だけを取り出す場合は

awk '((NR != 1)&&($1*$1 < 225)&&($2*$2 < 225)&&($3*$3 < 225)){print $1,$2,$3}' < input > output

てな感じで「一行野郎」でおしまい。どうでもいいのですが、 one liner のことを「一行野郎」と最初に訳したのはどこの誰でしょう? ともかく、こういう処理のために表計算ソフトを使ったり、 FORTRAN や C でプログラム書くのは異常に面倒だと思うのは私だけでしょうか。 (NR != 1) の部分は、わしはデータファイルの一行目に簡単なコメントを書く習慣があるので、一行目は読まないという意味です。 input と output が沢山ある場合は pick.sh てな感じでシェルスクリプトにしてしまって、

% sh pick.sh

とすると、簡単簡単。この中身をエディタで一括変換しちまえば、別の処理にすることもできるし。

因みに、awk の参考書としては

プログラミング言語 AWK
A.V.エイホ,B.W.カーニハン,P.J.ワインバーガー 著
足立高徳 訳
トッパン(1989) ISBN4-8101-8008-5

が良いと思います。awk の作者様ですし。処理系に依存した部分は awk-jp.info がどこかにあったのを読めば良いんではないでしょうか。

さて、gnuplot でプロット

gnuplot も、ほとんどの Linux のディストリビューションに入っています。一見、使いにくそうなソフトですが、スクリプトファイルの書き方がわかれば、目的のグラフを大量生産できます。最初は

% gnuplot

        G N U P L O T
        Linux version 3.5 (pre 3.6)
        patchlevel beta 347 (+1.1.8 1998/07/04)
        last modified Mon Jun 22 13:22:33 BST 1998

        Copyright(C) 1986 - 1993, 1998
        Thomas Williams, Colin Kelley and many others

        Send comments and requests for help to info-gnuplot@dartmouth.edu
        Send bugs, suggestions and mods to bug-gnuplot@dartmouth.edu

Terminal type set to 'x11'
gnuplot>

てな感じで、おそるおそるグラフを描いてみるのが良かろうと思いますが、気に入った設定になったら、

gnuplot> save "filename"

で、スクリプトファイルをセーブすることを覚えれば、あとは吐き出したスクリプトファイルをいじっていく間に使い方を覚えます。

今回のアニメーションでは、こんなスクリプトで、 pbm 形式のグラフを出力します。これを後で gif 形式に変換します。スクリプトの最後の方を見ていただけるとわかりますが、上のようなグラフを一気に20枚書き出しています。気に入るような設定になるように、最初の一枚を描くためのスクリプトを作ったら、あとは一気にできるという次第です。

因みに、gnuplot の参考書としては

使いこなす GNUPLOT
大竹敢 著
テクノプレス(1996)
ISBN4-924998-11-7

などが良いかも。Web 上にリファレンスの日本語訳があります。

GNUPLOT 日本語リファレンス・ガイド

pbm のグラフを gif に変換

先程の、 pbm 形式のグラフを見ていただければわかるように、そろそろ完成が見えてきてるなぁという感じです。で、次の目標としてはこんな感じに gif 形式のグラフにしたいです。
やりたいことは、
・配色を変える
・大きさを変える
・フォーマットを変える
ということですが、 Image Magic に含まれている convert でなんとかならないかと思ったのですが、なんだか処理が重いし出力がモノクロになっちまうし、。。
% convert -crop 790x740+76+152 @xyz01.pbm gif:@xyz01.gif
といった処理を期待しているのですが。。あ、convert では
% convert | grep PBM
PBM Portable bitmap format (black and white).
ということで、白黒になってしまうのでした。ということで、最初から gnuplot で吐き出すファイルを postscript にしてみたりしました。が、 convert の書き出す gif はファイルサイズが非常に大きいです。うー、これは gif のライセンスの問題だったという話をどこかで読んだような気がするが全くの嘘かもしれません。

ということで gimpで処理しました。
・ファイル→開く
・画像→インデックス画像
・画像→サイズの変更、で 790x740+76+152 に。
・ダイアログ→インデックスパレット、で右クリックで一つ一つの色(全6色)を変えて
・ファイル→別名で保存
わしの gimp はすっかり日本語化されているので、説明が windows のソフトの使い方みたいですみません。 gimp のスクリプトの書き方を学ばないと。。でもまぁ、とりあえず出来ればいいんだ、いいんだ。。

でも、やっぱり大量生産するためにはスクリプトの書き方を学ぼうということで、 gimp でバッチ処理をしてみました。これは、このページの続編です。

20 枚の gif を gifmerge でアニメーションに

gimp でレイヤを駆使してアニメーションを作っても良いのですが (実際、最初はわしはそうした)、コマ数が多くなると gimp が喘ぐ。。ということで、 gifmerge なんですが、siag を日本語化された metalさんのページにある plamo 用 siag のバイナリの中に入っていたものを流用させていただいています。使い方は簡単で
% gifmerge -5 -l0 @xyz*.gif > animation.gif
で、あっというまに出来ます。最初の -5 は 5/100 秒ごとにコマを送る、次の -l0 はループは無限回、@xyz*.gif は @xyz01.gif から @xyz20.gif まで 20枚の元の画像があるということです。ファイル名の頭に @ をつけている理由は、こうした処理の中間ファイル全てを @ から始まるファイル名にしておくと、cron などで要らんファイルを消去するときの目印になるからです。これ、 @xyz1.gif,@xyz2.gif ,...,@xyz20.gif というファイル名にしてしまうと @xyz1.gif,@xyz10.gif,@xyz11.gif,...,@xyz19.gif,@xyz2.gif,@xyz3.gif,... という順番になってしまいます(わしはそうしてしまった。。)。
まぁ、そういう間抜けな数々の苦労をした挙げ句、めでたく完成ということで。一度、大きなアニメが出来てしまえば上のような小さな画像に切り取りたい場合は gimp で普通に crop すれば出来上がりです。