お知らせ+活動記録+たわごと

HP と Twitter を補完するとともに、互いの密接な連携を図るため、本ブログを開設した。三位一体を目指す。情報提供、広報活動、教育・啓蒙活動の一環として、肩の力を抜き、冗長性を廃し、簡にして要を得た文章を書くよう心がける。
<< December 2017 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 >>
 
MOBILE
qrcode
PROFILE
無料ブログ作成サービス JUGEM
 
反射リストからデータを抽出して CIF を作成するシェルスクリプト
RIETAN-FP で NPRINT > 0 に設定してリートベルト解析を実行すると、"Summary of possible reflections (based on the refined parameters)" というタイトルの反射リストが標準出力ファイル hoge.lst の末尾近くに記録される。このリストから
  1. h
  2. k
  3. l
  4. 格子面間隔 d
  5. 観測積分強度 Io
  6. 結晶構造因子 |F|
  7. 半値全幅 H
を取り出し、各データをスペースで区切って Crystallographic Information File (CIF) の一部として出力するパートをシェルスクリプト lst2cif.command に追加してみた。当該部分は lst2cif.command が hoge.lst と hoge.dst から hoge.cif を作成する際、hoge.cif の末尾に反射リフトを付加する役割を受け持つ。

Cu Kα 特性X線で測定した粉末回折データの場合、Kα2反射のデータがリストに含まれる。また、部分的にプロファイルが欠落した高角反射(Iobs: H)の行も記録される。多相試料のリートベルト解析では、不純物相の反射も出力される。これらの不要な行は削除し、残った行から必要なデータをピックアップして標準出力に流せばよい。lst2cif.command 中の反射リスト出力部分(本エントリー用に若干変更)は次の通り:

# Labels for associated multiple data values
cat << EOS

_loop
   _pd_peak_id
   _refln_index_h
   _refln_index_k
   _refln_index_l
   _pd_peak_2theta_centroid
   _pd_peak_d_spacing
   _pd_peak_intensity
   _rietan_i100_meas
   _refln_F_cal
   _pd_peak_width_2theta
EOS

sed -n '/Iobs *Ical/,/^$/p' $1 | sed '1d;$d' | \ (1)
awk -v c9=$(grep "_pd_peak_intensity =" $1 | sed 's/^.* \* //') \ (2)
   '$2 == 1 && match($6, /[+-]2/) == 0 && $9 != "H" \ (3)
   { i++; if ($9 == "-") print i, $3, $4, $5, $7, $8, $9, $9, $11, $14; \ (4)
   else printf "%s %s %s %s %s %s %10.5f %6.2f %s %s\n", \ (5)
   i, $3, $4, $5, $7, $8, c9*$9, 0.001*$9, $11, $14 }' | \ (6)
sed -r 's/\s\s+/ /g' (7)

exit

_rietan_i100_meas は最強反射の観測積分強度を 100 としたときの相対強度の自家製定義である。本来はピーク強度を対象とすべきなのだろうが、hoge.lst を処理対象としているため、積分強度を対象として計算せざるを得なかった。

この bash スクリプトはややこしいリスト処理をわずか一行の呪文でやってのける。sed と awk の合わせ技である。その行はかなり長いので、便宜上 (1)(7) に 7 分割した。

(1) では、sed の連続技により hoge.lst から反射リストの部分を抜き取る。オプション -n は suppress automatic printing of pattern space を意味する。

(2) では、grep と sed をパイプでつなぎ hoge.lst から抽出した数値をオプション -v var=val で変数 c9 として awk に渡す。

(3) は awk のパターン部分であり、第 1 相以外の反射、Kα2反射、2θmax 付近の反射(9番目のフィールドが "H")を排除する。必要に応じて (3) を変更すれば、CIF に出力する行を変えられる。

(4)(6) が awk のアクション部分であり、フィールド(3〜5, 7〜9, 11, 14)を指示して出力する。まず (4) において Iobs が '-'(一部の観測データが欠落)の反射をそのまま出力する。(5)(6) がすべてのデータがそろっている反射の処理に相当する。最強反射で 100 000 となるように正規化した観測積分強度 Iobs にそれぞれ c9 と 0.001を乗じた値をフィールド No. 9 と 10 として出力する。

(7) では sed によりフィールド間を単一スペースに統一する。そこまで徹底する必要はないが、いざとなれば特定のフィールドを cut で抽出できるというメリットがある。

たとえば Fapatite.lst を処理すると、次の出力が得られる:

_loop
   _pd_peak_id
   _refln_index_h
   _refln_index_k
   _refln_index_l
   _pd_peak_2theta_centroid
   _pd_peak_d_spacing
   _pd_peak_intensity
   _refln_F_cal
   _pd_peak_width_2theta
1 1 0 1 16.877 5.24918 1.21318 5.35 11.0711 0.07402
2 2 0 0 21.891 4.05690 2.34131 10.33 29.9436 0.07424
3 1 1 1 22.945 3.87283 2.10185 9.27 21.1948 0.07430
4 2 0 1 25.464 3.49510 0.62246 2.75 12.8217 0.07449
5 0 0 2 25.864 3.44192 13.12474 57.88 147.1135 0.07452
6 1 0 2 28.139 3.16861 4.29895 18.96 37.4418 0.07472
7 2 1 0 29.095 3.06673 5.02640 22.17 58.8922 0.07482
8 1 2 0 29.095 3.06673 0.46282 2.04 17.8710 0.07482
9 1 2 1 31.921 2.80131 9.46481 41.74 66.0020 0.07514
・・・

エディターや表計算ソフトなどで GUI を通じて同様の表処理を繰り返すのは、あまりにも泥臭く面倒だ。場合分けを含む表形式テキストデータの複雑な処理は awk の独壇場といって過言でない。このスクリプトをひな形として書き換えれば、RIETAN-FP などが出力する種々のテキストファイルから、必要に応じて任意のデータを抽出してテキストファイルに出力できる。ぜひ活用していただきたい。
コメント
コメントする









 

(C) 2017 ブログ JUGEM Some Rights Reserved.