plot3d.py 技術ドキュメント

プログラムの動作

plot3d.py は、VASP (Vienna Ab initio Simulation Package) の出力ファイルである電荷密度（CHGCAR）や静電ポテンシャル（LOCPOT）などの三次元体積データを解析し、Matplotlib を用いて3Dグラフィックとして可視化するPythonスクリプトです。

主な機能は以下の通りです。

  1. VASPファイルの読み込み: CHGCARやLOCPOTといったVASPの体積データファイルを読み込み、格子情報や三次元グリッドデータを抽出します。

  2. 等値面 (Isosurfaces) の描画: 指定された値（レベル）において、電荷密度や電位が一定となる面を3Dで表示します。

  3. グリッド点 (Dots) の描画: 体積データが特定の閾値を超えるグリッド点を散布図として3D空間にプロットします。

  4. スライス面 (Slice Planes) の描画: 分率座標で指定された位置で体積データをスライスし、その断面図を色付きの平面として3D空間に表示します。

  5. ダウンサンプリング: 大規模なデータセットに対する処理速度向上のため、体積データを指定したステップで間引くことができます。

本プログラムは、VASP計算結果から得られる複雑な三次元的な物理量分布を直感的に理解しやすくするための強力な可視化ツールとして機能します。

原理

plot3d.py は、以下の数式やアルゴリズムを用いて三次元体積データを可視化します。

1. CHGCARファイルの解析

VASPのCHGCAR（または類似のフォーマット）ファイルは、結晶の格子情報、原子の位置、そして三次元グリッド上にマッピングされた体積データ（電荷密度など）を含んでいます。プログラムはこれらの情報を順次読み込み、特に格子ベクトル \(\mathbf{a}, \mathbf{b}, \mathbf{c}\) とグリッドデータ \(F(i, j, k)\) （\(i, j, k\) は各軸のグリッドインデックス）を抽出します。

2. 直交座標メッシュの生成

VASPのデータは分率座標系で定義されます。これを物理的な距離を示す直交座標系に変換します。グリッド点 \((i, j, k)\) の分率座標 \((f_x, f_y, f_z)\) は、グリッド数 \(N_x, N_y, N_z\) を用いて以下のように計算されます。 $\(f_x = \frac{i}{N_x}, \quad f_y = \frac{j}{N_y}, \quad f_z = \frac{k}{N_z}\)\( これらの分率座標は、格子ベクトル \)\mathbf{a}=(a_x, a_y, a_z)\(, \)\mathbf{b}=(b_x, b_y, b_z)\(, \)\mathbf{c}=(c_x, c_y, c_z)\( を用いて、直交座標 \)(X, Y, Z)\( に変換されます。 \)\(\mathbf{R} = X\mathbf{e}_x + Y\mathbf{e}_y + Z\mathbf{e}_z = f_x \mathbf{a} + f_y \mathbf{b} + f_z \mathbf{c}\)\( より具体的には、各成分について以下のように展開されます。 \)\(X = f_x a_x + f_y b_x + f_z c_x\)\( \)\(Y = f_x a_y + f_y b_y + f_z c_y\)\( \)\(Z = f_x a_z + f_y b_z + f_z c_z\)\( ここで \)\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z$ はデカルト座標系の単位ベクトルです。

3. 等値面 (Isosurfaces) の描画

等値面の描画には、scikit-image ライブラリの marching_cubes アルゴリズムが利用されます。このアルゴリズムは、三次元体積データ内の特定の閾値（level）を持つ点を結んでポリゴンメッシュ（頂点と面）を生成します。 生成されたメッシュの頂点座標は分率座標で得られるため、上記2.の変換式を用いて直交座標に変換され、matplotlib の Poly3DCollection を使って描画されます。

4. スライス面 (Slice Planes) の描画

スライス面は、特定の分率座標値で体積データを切り取った2D平面です。例えば、--slice-xy 0.5 は \(f_z = 0.5\) の位置でデータをスライスします。 スライスされた2Dグリッドの各セルについて、隣接するグリッド点の値の平均値を算出し、その平均値に応じて色を決定します。この色付けされた平面は、tklib.tkgraphic.tkPlot3d.plot_surface3d 関数を使用して描画されます。

5. グリッド点 (Dots) の描画

特定の閾値（--cutoff）またはパーセンタイル（--quantile）を超える体積データを持つグリッド点を抽出し、その直交座標を tklib.tkgraphic.tkPlot3d.plot_scatter3d 関数を用いて散布図として描画します。これにより、高密度または低密度の領域を視覚的に強調できます。

6. ダウンサンプリング

--ds DSX DSY DSZ オプションを使用すると、体積データは各軸で \(DSX, DSY, DSZ\) ステップごとに間引かれます。これにより、計算量が \(1/(DSX \cdot DSY \cdot DSZ)\) に削減され、特に大規模なデータセットの処理速度が向上します。

必要な非標準ライブラリとインストール方法

plot3d.py の実行には、以下の非標準Pythonライブラリが必要です。

  • numpy: 数値計算を効率的に行うためのライブラリです。

  • matplotlib: グラフ描画（特に3Dプロット）のためのライブラリです。

  • scikit-image: 等値面を生成する marching_cubes アルゴリズムを提供します。

  • tklib: tkPlot3d モジュールを含む、このプログラム作者が作成したカスタムライブラリです。3Dプロットのユーティリティ関数を提供します。

インストール方法

numpy, matplotlib, scikit-image は pip コマンドでインストールできます。

pip install numpy matplotlib scikit-image

tklib は一般的な pip リポジトリでは提供されていません。プログラムコード内の sys.path.append('d:/git/tkProg/tklib/python') から、このライブラリは開発者のローカル環境に配置されていると推測されます。

tklib のインストール（または配置）については、以下のいずれかの方法を取る必要があります。

  1. 既存の tklib モジュールを入手: もし tklib のソースコードが公開されている場合は、それをダウンロードし、上記の sys.path.append で指定されているディレクトリパス（例: d:/git/tkProg/tklib/python）またはそれに相当するパスに配置する必要があります。

  2. sys.path.append の変更: tklib モジュールを他のディレクトリに配置した場合は、plot3d.py ファイル内の sys.path.append('d:/git/tkProg/tklib/python') の行を、tklib モジュールが置かれているパスに修正してください。

必要な入力ファイル

plot3d.py は、VASPの体積データファイルを入力として期待します。

  • ファイル名: infile 引数で指定します。例: CHGCAR, LOCPOT。

  • ファイル形式: VASPの標準的な体積データファイル形式に準拠している必要があります。

典型的なファイル構造は以下の通りです。

  1. タイトル行: コメントやファイルの説明。

  2. スケール因子: 格子ベクトルのスケール因子。

  3. 格子ベクトル: 3行3列の行列で、各行が結晶の基本格子ベクトル \(\mathbf{a}, \mathbf{b}, \mathbf{c}\) を表します。

  4. 原子の種類: 各原子の元素記号（例: C O）。

  5. 原子の数: 各種類の原子の数（例: 10 20）。

  6. 座標タイプ: 座標が直交（Cartesian）か分率（Direct）か。

  7. 原子座標: 各原子の座標。

  8. 空行: 原子座標とグリッド情報の間には通常、空行が存在します。このプログラムは空行をスキップしてグリッド情報を探します。

  9. グリッドサイズ: \(N_x, N_y, N_z\) の3つの整数で、体積データの各軸のグリッド点数を表します。

  10. 体積データ: \(N_x \times N_y \times N_z\) 個の浮動小数点数値が複数行にわたって記述されます。これらが三次元グリッド上の物理量（電荷密度、電位など）の値を表します。

生成される出力ファイル

本プログラムは、標準ではMatplotlibのウィンドウに3Dプロットを表示します。 --save オプションが指定された場合、そのファイル名で3Dプロットが画像ファイルとして保存されます。

  • ファイル名: --save <filename> で指定されたファイル名（例: output.png, plot.pdf）。

  • ファイル形式: matplotlib がサポートする様々な画像形式（PNG, JPEG, PDF, SVGなど）に対応します。ファイル拡張子によって自動的に決定されます。

  • 内容: 入力されたVASP体積データの等値面、グリッド点、またはスライス面を可視化した3Dグラフィックが含まれます。図のタイトル、軸ラベル、凡例、カラーバーなどが表示されます。

コマンドラインでの使用例 (Usage)

基本的な実行コマンドと引数の説明は以下の通りです。

python plot3d.py INFILE [OPTIONS]

  • INFILE: VASPの体積データファイル（例: CHGCAR, LOCPOT）。

  • [OPTIONS]: 以下の引数が利用可能です。

    • --mode {iso,dots,both}: 表示モード。iso（等値面）、dots（グリッド点）、both（等値面とスライス面）から選択。デフォルトは both。

    • --levels N [N ...]: 等値面の絶対値レベルをスペース区切りで指定（例: 0.05 0.1）。

    • --nlevels N: 最小値から最大値の間にN個の等値面を自動生成。

    • --alpha FLOAT: 等値面またはグリッド点の透明度（0.0～1.0）。デフォルトは 0.30。

    • --edge STR: 等値面のエッジの色。none で非表示。

    • --lw FLOAT: 等値面のエッジの線幅。デフォルトは 0.0。

    • --no-legend: 凡例を非表示にする。

    • --cutoff FLOAT: dots モードで、この値以上のグリッド点のみをプロット。

    • --quantile FLOAT: dots モードで、このパーセンタイル値以上のグリッド点のみをプロット（例: 0.9 は上位10%）。

    • --subsample INT: dots モードで、選択された点をさらに間引くステップ数。

    • --max-points INT: dots モードで、プロットする最大点の数。

    • --size FLOAT: dots モードで、プロットする点のサイズ。デフォルトは 0.4。

    • --cmap STR: dots モードで使用するカラーマップ名。デフォルトは viridis。

    • --slice-xy N [N ...]: 分率\(z\)座標を固定したXYスライス面の表示（例: 0.25 0.5）。

    • --slice-yz N [N ...]: 分率\(x\)座標を固定したYZスライス面の表示。

    • --slice-zx N [N ...]: 分率\(y\)座標を固定したZXスライス面の表示。

    • --slice-cmap STR: スライス面で使用するカラーマップ名。デフォルトは viridis。

    • --slice-alpha FLOAT: スライス面の透明度。デフォルトは 0.6。

    • --float32: 体積データの処理に float32 を使用してメモリを節約し高速化。

    • --ds DSX DSY DSZ: 各軸のダウンサンプリングステップ数（例: 2 2 2）。

    • --ortho: 正射影モードを有効にする。

    • --pad FLOAT: 軸の表示範囲にパディングを追加。

    • --save FILENAME: 生成された画像をファイルに保存。

    • --title STR: プロットのタイトルを指定。

コマンドラインでの具体的な使用例

ここでは、架空の CHGCAR および LOCPOT ファイルを例として、具体的な使用例とその結果について説明します。

例1: 等値面を複数レベルでプロットし、PNGファイルに保存

電荷密度 CHGCAR の0.05と0.10という2つの等値面を描画し、結果を chgcar_iso.png として保存します。

python plot3d.py CHGCAR --levels 0.05 0.1 --save chgcar_iso.png

実行結果: chgcar_iso.png ファイルが生成されます。この画像には、0.05レベルの等値面（通常はより大きく広がる）と0.10レベルの等値面（通常は原子核に近い領域に限定される）が異なる色で描画され、それぞれの凡例が表示されます。等値面は透明度0.3で半透明に表示され、3D空間における電荷密度の分布が視覚化されます。

例2: 静電ポテンシャルの高値領域を散布図でプロット

静電ポテンシャル LOCPOT の値が上位10%にあるグリッド点を散布図としてプロットし、hot カラーマップを使用してプロットサイズを小さくします。結果を locpot_dots.png として保存します。

python plot3d.py LOCPOT --mode dots --quantile 0.9 --cmap hot --size 0.2 --save locpot_dots.png

実行結果: locpot_dots.png ファイルが生成されます。画像には、高い静電ポテンシャルを持つ空間領域が小さな点（サイズ0.2）の集まりとして表現されます。点の色はhotカラーマップに基づき、値が高いほど明るい色で表示され、カラーバーによって値と色の対応が示されます。これにより、ポテンシャルが特に高い、または低い領域の空間分布を把握できます。

例3: 複数スライス面を描画し、アルファ値を調整

電荷密度 CHGCAR の \(z=0.25, 0.75\) (分率座標)、\(x=0.5\) (分率座標)、\(y=0.5\) (分率座標) の位置でスライス面を描画し、スライス面の透明度を0.2に設定します。結果を chgcar_slices.png として保存します。

python plot3d.py CHGCAR --slice-xy 0.25 0.75 --slice-yz 0.5 --slice-zx 0.5 --slice-alpha 0.2 --save chgcar_slices.png

実行結果: chgcar_slices.png ファイルが生成されます。この画像には、結晶構造の異なる断面が3枚の半透明な平面として表示されます。各スライス面は電荷密度の分布に応じて色付けされ、その透明度により内部の構造が透けて見えるようになります。これにより、特定の平面上での電荷密度の変化を詳細に観察できます。

例4: 等値面とスライス面を同時にプロットし、データをダウンサンプリング

電荷密度 CHGCAR の等値面3つと \(z=0.5\) のスライス面を同時にプロットします。描画負荷軽減のため、データを各軸で2ステップずつダウンサンプリング（グリッド数を1/8に削減）します。結果を chgcar_both_ds.png として保存します。

python plot3d.py CHGCAR --mode both --nlevels 3 --slice-xy 0.5 --ds 2 2 2 --save chgcar_both_ds.png

実行結果: chgcar_both_ds.png ファイルが生成されます。画像には、ダウンサンプリングされたデータから生成された3つの等値面と、\(z=0.5\) の位置にあるスライス面が同時に描画されます。等値面とスライス面が組み合わされることで、体積全体の概観と特定の断面の詳細を一度に確認でき、データの間引きにより大規模なファイルでも比較的スムーズに可視化されます。