#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
VASP体積データ可視化ツール (plot3d.py)
概要:
VASPのCHGCARファイルから電荷密度やELFなどの体積データを読み込み、
3D空間に等値面、高密度点群、またはスライス面として可視化するスクリプト。
詳細説明:
本スクリプトは、VASP計算で得られた電荷密度やELFなどの体積データファイル（CHGCAR形式）を解析し、
そのデータをmatplotlibとscikit-image (skimage) のmarching_cubesアルゴリズム、
およびカスタムのtkPlot3dユーティリティを使用して3次元で描画します。
ユーザーは、等値面レベルの指定、データ点群のフィルタリング、特定の分率座標でのスライス表示、
ダウンサンプリングによるパフォーマンス向上などのオプションを通じて、
データの詳細な分析と美しい図の生成が可能です。
matplotlibの3Dプロット機能に依存しており、等値面生成にはscikit-imageが必要です。
関連リンク:
:doc:`plot3d_usage` （本スクリプトの使用方法に関する詳細ドキュメント）
"""
import sys
import argparse
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
try:
from skimage.measure import marching_cubes
except:
print(f"\nError: Can not import skimage")
print(f" Install: pip install scikit-image")
input("\nPress ENTER to terminate>>\n")
exit()
sys.path.append('d:/git/tkProg/tklib/python')
import tklib.tkgraphic.tkPlot3d as tk3d
# --- CHGCAR 読み込み（空行スキップ対応 + dtype 選択） ---
[ドキュメント]
def read_chgcar(filename, dtype=np.float32):
"""
VASPのCHGCARファイルを読み込み、ヘッダー情報と体積データを抽出する。
概要:
指定されたVASP CHGCARファイルから結晶情報と3次元の体積データを読み込む。
詳細説明:
CHGCARファイルのフォーマットに従い、タイトル、スケール因子、格子ベクトル、
原子の種類と数、座標タイプ、原子座標、グリッドサイズ、そして電荷密度などの
体積データをパースする。空行のスキップに対応し、データの読み込み時にNumPyの
データ型 (dtype) を指定できる。
:param filename: str, 読み込むCHGCARファイルのパス。
:param dtype: numpy.dtype, 体積データの格納に使用するNumPyデータ型 (例: np.float32, np.float64)。
:returns: tuple,
- info (dict): CHGCARヘッダーから抽出された情報を格納した辞書。
'title', 'scale', 'lattice', 'atoms', 'numbers',
'coord_type', 'coords', 'grid' のキーを含む。
- P (list of lists of lists): グリッドデータ `F` を [[[F[i,j,k] for k] for j] for i] 形式にしたもの。
（既存コードの形式を維持）
- F (numpy.ndarray): 読み込まれた体積データを格納した (nx, ny, nz) 形状の3次元NumPy配列。
:rtype: tuple[dict, list, numpy.ndarray]
"""
info = {}
with open(filename, 'r') as f:
info['title'] = f.readline().strip()
info['scale'] = float(f.readline().strip())
lattice = []
for _ in range(3):
lattice.append([float(x) for x in f.readline().split()])
info['lattice'] = np.array(lattice, dtype=dtype)
info['atoms'] = f.readline().split()
info['numbers'] = [int(x) for x in f.readline().split()]
info['coord_type'] = f.readline().strip()
natoms = sum(info['numbers'])
coords = []
for _ in range(natoms):
l = f.readline()
coords.append([float(x) for x in l.split()[:3]])
info['coords'] = np.array(coords, dtype=dtype)
# 空行スキップしてグリッド
while True:
l = f.readline()
if not l:
raise ValueError("Unexpected EOF before grid size.")
if l.strip():
break
nx, ny, nz = [int(x) for x in l.split()]
info['grid'] = (nx, ny, nz)
# 体積データ
ngrid = nx * ny * nz
data = []
while len(data) < ngrid:
line = f.readline()
if not line:
break
data.extend([float(x) for x in line.split()])
if len(data) != ngrid:
raise ValueError("CHGCAR data size mismatch")
F = np.array(data, dtype=dtype).reshape((nx, ny, nz), order='F')
P = [[[F[i, j, k] for k in range(nz)] for j in range(ny)] for i in range(nx)]
return info, P, F
# --- 直交座標メッシュ生成（分率→直交） ---
[ドキュメント]
def build_cartesian_grid(lattice, nx, ny, nz, dtype=np.float32):
"""
分率座標から直交座標への変換に基づき、3次元直交座標グリッドを構築する。
概要:
指定された格子ベクトルとグリッド分割数を用いて、空間の各点に対応する
直交座標のメッシュデータを生成する。
詳細説明:
まず、各軸に沿って正規化された分率座標 (0から1まで) の配列を生成し、
これらを `np.meshgrid` で組み合わせて3次元分率座標グリッドを作成する。
その後、この分率座標を `lattice` (格子ベクトル) を使って直交座標に変換する。
これにより、体積データが定義されている空間の物理的なX, Y, Z座標が計算される。
:param lattice: numpy.ndarray, (3, 3) 形状の格子ベクトル行列。各行が基底ベクトル。
:param nx: int, X軸方向のグリッド点数。
:param ny: int, Y軸方向のグリッド点数。
:param nz: int, Z軸方向のグリッド点数。
:param dtype: numpy.dtype, 座標値の格納に使用するNumPyデータ型。
:returns: tuple,
- X (numpy.ndarray): X座標の3次元グリッド (nx, ny, nz)。
- Y (numpy.ndarray): Y座標の3次元グリッド (nx, ny, nz)。
- Z (numpy.ndarray): Z座標の3次元グリッド (nx, ny, nz)。
:rtype: tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]
"""
fx = np.arange(nx, dtype=dtype) / dtype(nx)
fy = np.arange(ny, dtype=dtype) / dtype(ny)
fz = np.arange(nz, dtype=dtype) / dtype(nz)
FX, FY, FZ = np.meshgrid(fx, fy, fz, indexing='ij') # (nx,ny,nz)
a, b, c = lattice[0], lattice[1], lattice[2]
X = FX * a[0] + FY * b[0] + FZ * c[0]
Y = FX * a[1] + FY * b[1] + FZ * c[1]
Z = FX * a[2] + FY * b[2] + FZ * c[2]
return X, Y, Z
# --- ダウンサンプリング（各軸のステップ） ---
[ドキュメント]
def downsample_volume(F, dsx, dsy, dsz):
"""
3次元体積データを指定されたステップでダウンサンプリングする。
概要:
NumPy配列として与えられた体積データ `F` を、各軸 (X, Y, Z) に沿って
指定された間隔で間引くことでダウンサンプリングする。
詳細説明:
データの各軸に対してスライシング `[::step]` を適用し、
より少ないデータ点数で元のボリュームを表現する。
これにより、可視化処理のパフォーマンスを向上させることができる。
:param F: numpy.ndarray, (nx, ny, nz) 形状の元の3次元体積データ。
:param dsx: int, X軸方向のダウンサンプリングステップサイズ。1はダウンサンプリングなし。
:param dsy: int, Y軸方向のダウンサンプリングステップサイズ。1はダウンサンプリングなし。
:param dsz: int, Z軸方向のダウンサンプリングステップサイズ。1はダウンサンプリングなし。
:returns: numpy.ndarray, ダウンサンプリングされた3次元体積データ。
:rtype: numpy.ndarray
"""
return F[::dsx, ::dsy, ::dsz]
[ドキュメント]
def fractional_spacing_after_downsampling(nx, ny, nz, dsx, dsy, dsz):
"""
ダウンサンプリング後のデータにおける分率座標でのグリッド間隔を計算する。
概要:
元のグリッドサイズとダウンサンプリングステップから、
ダウンサンプリング後の論理的な分率間隔を算出する。
詳細説明:
`skimage.measure.marching_cubes` 関数に渡す `spacing` パラメータは、
実際のデータ点間の相対的な間隔を表す。ダウンサンプリングを行った場合、
この間隔は元のグリッド間隔にダウンサンプリングステップを掛けたものになる。
例えば、nx=100でdsx=2の場合、元の分率間隔は1/100だが、ダウンサンプリング後は
2/100となる。
:param nx: int, 元のX軸方向のグリッド点数。
:param ny: int, 元のY軸方向のグリッド点数。
:param nz: int, 元のZ軸方向のグリッド点数。
:param dsx: int, X軸方向のダウンサンプリングステップサイズ。
:param dsy: int, Y軸方向のダウンサンプリングステップサイズ。
:param dsz: int, Z軸方向のダウンサンプリングステップサイズ。
:returns: tuple, 各軸の分率間隔 (spacing_x, spacing_y, spacing_z)。
:rtype: tuple[float, float, float]
"""
return (dsx / nx, dsy / ny, dsz / nz)
# --- 等値面（分率で抽出→直交に変換） ---
[ドキュメント]
def plot_isosurfaces_fractional(ax, F, lattice, levels, alpha=0.3, edgecolor='none', linewidth=0.0,
colors=None, spacing_frac=(1,1,1), legend=True, antialias=False):
"""
3次元体積データから等値面を抽出し、直交座標系で描画する。
概要:
`skimage.measure.marching_cubes` を使用して指定された等値レベルの面を抽出し、
それを直交座標に変換してmatplotlibの3D軸上に描画する。
詳細説明:
`F` は分率座標での体積データと見なされ、`marching_cubes` 関数に
`spacing_frac` を指定して等値面を抽出する。得られた頂点データは
分率座標であるため、`lattice` (格子ベクトル) を用いて直交座標に変換される。
その後、`tk3d.Poly3DCollection` を使ってメッシュとして描画される。
複数の等値レベルに対応し、それぞれ異なる色で描画することが可能。
凡例の表示もサポートする。
:param ax: matplotlib.axes.Axes, 3D描画を行うmatplotlibの軸オブジェクト。
:param F: numpy.ndarray, (nx, ny, nz) 形状の3次元体積データ。
:param lattice: numpy.ndarray, (3, 3) 形状の格子ベクトル行列。各行が基底ベクトル。
:param levels: list of float, 描画する等値レベルのリスト。
:param alpha: float, 等値面の透明度 (0.0:完全に透明, 1.0:完全に不透明)。
:param edgecolor: str or tuple, 等値面の辺の色。'none' で辺を描画しない。
:param linewidth: float, 等値面の辺の太さ。
:param colors: list of str or tuple, 各等値レベルに対応する色のリスト。Noneの場合、デフォルトのカラーマップが使用される。
:param spacing_frac: tuple of float, `marching_cubes` に渡す各軸の分率グリッド間隔 (dx, dy, dz)。
:param legend: bool, 等値レベルの凡例を表示するかどうか。
:param antialias: bool, 等値面をアンチエイリアス処理するかどうか。
:returns: None
"""
patches = []
if colors is None:
cmap = plt.get_cmap('tab10')
colors = [cmap(i % 10) for i in range(len(levels))]
for il, level in enumerate(levels):
try:
verts_frac, faces, _, _ = marching_cubes(F, level=level, spacing=spacing_frac)
except Exception:
continue
verts_cart = verts_frac @ lattice # (N,3) @ (3,3)
mesh = tk3d.Poly3DCollection(verts_cart[faces], alpha=alpha,
facecolor=colors[il], edgecolor=edgecolor,
linewidth=linewidth, antialiased=antialias)
ax.add_collection3d(mesh)
if legend:
patches.append(Patch(facecolor=colors[il], edgecolor='k', label=f"Level={level:.4g}"))
if legend and patches:
ax.legend(handles=patches, loc="upper right", fontsize=9)
# --- スライス面を正しく直交座標に配置して描画 ---
def _facecolors_from_vals(vals, cmap):
"""
スカラー値グリッドからmatplotlibの`plot_surface3d`用のFaceColorsを生成する。
概要:
2次元のスカラー値配列 `vals` を受け取り、各セルの中心での値を計算し、
指定されたカラーマップを適用してRGBA形式のFaceColors配列を生成する。
詳細説明:
`vals` は `(M, N)` 形状のグリッドデータで、これはグリッド点での値を示す。
`plot_surface3d` の `facecolors` 引数は、各面の色を `(M-1, N-1, 4)` 形状の
RGBA配列として期待する。この関数では、隣接する4つの点の値の平均をセルの中心値とし、
その値を正規化してカラーマップで色に変換する。
:param vals: numpy.ndarray, (M, N) 形状の2次元スカラー値グリッド。
:param cmap: str or matplotlib.colors.Colormap, 使用するカラーマップの名前またはオブジェクト。
:returns: numpy.ndarray, (M-1, N-1, 4) 形状のRGBA値のNumPy配列。
:rtype: numpy.ndarray
"""
import numpy as np
import matplotlib.pyplot as plt
norm = plt.Normalize(vmin=float(vals.min()), vmax=float(vals.max()))
vc = 0.25 * (vals[:-1, :-1] + vals[1:, :-1] + vals[:-1, 1:] + vals[1:, 1:])
return plt.get_cmap(cmap)(norm(vc))
[ドキュメント]
def render_slice_plane(ax, lattice, F, which, frac, cmap='viridis', alpha=0.6):
"""
指定された分率座標に沿ってスライス面を抽出し、直交座標系で描画する。
概要:
3D体積データ `F` から、特定の分率座標 (`frac`) で切り取られた2Dスライス面を抽出し、
それを直交座標に変換して3Dプロットに描画する。
詳細説明:
`which` パラメータ ('xy', 'yz', 'zx') に応じて、対応する軸の分率座標 (`frac`) を固定し、
他の2軸に沿ったスライスデータを抽出する。このスライスデータから `_facecolors_from_vals`
を用いて各面のRGBA色を生成し、格子ベクトル `lattice` を使って直交座標を計算する。
最終的に `tk3d.plot_surface3d` を用いて、色付けされたスライス面を3D軸上に配置する。
:param ax: matplotlib.axes.Axes, 3D描画を行うmatplotlibの軸オブジェクト。
:param lattice: numpy.ndarray, (3, 3) 形状の格子ベクトル行列。各行が基底ベクトル。
:param F: numpy.ndarray, (nx, ny, nz) 形状の3次元体積データ。
:param which: str, スライスする平面の方向 ('xy', 'yz', 'zx')。
- 'xy': fz (Z軸の分率座標) が `frac` に固定される。
- 'yz': fx (X軸の分率座標) が `frac` に固定される。
- 'zx': fy (Y軸の分率座標) が `frac` に固定される。
:param frac: float, スライスする分率座標 (0.0から1.0の範囲)。
:param cmap: str or matplotlib.colors.Colormap, スライス面の描画に使用するカラーマップ。
:param alpha: float, スライス面の透明度 (0.0:完全に透明, 1.0:完全に不透明)。
:returns: None
:raises ValueError: `which` が 'xy', 'yz', 'zx' のいずれでもない場合。
"""
nx, ny, nz = F.shape
a, b, c = lattice[0], lattice[1], lattice[2]
if which == 'xy':
# fz = frac
fx = np.arange(nx) / nx
fy = np.arange(ny) / ny
FX, FY = np.meshgrid(fx, fy, indexing='ij') # (nx, ny)
fz = np.full_like(FX, frac, dtype=FX.dtype)
X = FX * a[0] + FY * b[0] + fz * c[0]
Y = FX * a[1] + FY * b[1] + fz * c[1]
Z = FX * a[2] + FY * b[2] + fz * c[2]
iz = int(round(frac * (nz - 1)))
vals = F[:, :, iz] # (nx, ny)
elif which == 'yz':
# fx = frac
fy = np.arange(ny) / ny
fz = np.arange(nz) / nz
FY, FZ = np.meshgrid(fy, fz, indexing='ij') # 形状: (ny, nz)
fx = np.full_like(FY, frac, dtype=FY.dtype)
X = fx * a[0] + FY * b[0] + FZ * c[0]
Y = fx * a[1] + FY * b[1] + FZ * c[1]
Z = fx * a[2] + FY * b[2] + FZ * c[2]
ix = int(round(frac * (nx - 1)))
# ★転置は不要です。FY/FZ は (ny, nz) なので、(ny, nz) のままにする
vals = F[ix, :, :] # ← ここを .T しない！
elif which == 'zx':
# fy = frac
fx = np.arange(nx) / nx
fz = np.arange(nz) / nz
FX, FZ = np.meshgrid(fx, fz, indexing='ij') # (nx, nz)
fy = np.full_like(FX, frac, dtype=FX.dtype)
X = FX * a[0] + fy * b[0] + FZ * c[0]
Y = FX * a[1] + fy * b[1] + FZ * c[1]
Z = FX * a[2] + fy * b[2] + FZ * c[2]
iy = int(round(frac * (ny - 1)))
vals = F[:, iy, :] # (nx, nz)
else:
raise ValueError("which must be one of 'xy', 'yz', 'zx'")
# セル中心の facecolors を作る（(M-1,N-1,4)）
facecolors = _facecolors_from_vals(vals, cmap)
# 描画（X,Y,Z は (M,N)、facecolors は (M-1,N-1,4) でOK）
tk3d.plot_surface3d(ax, X, Y, Z, facecolors=facecolors, edgecolor='none', alpha=alpha, shade=False)
# ユーティリティ：小数配列パース（例: "0.25 0.5"）
[ドキュメント]
def parse_fracs(str_list):
"""
文字列のリストから浮動小数点数のリストをパースするユーティリティ関数。
概要:
コマンドライン引数などで渡される文字列のリストを浮動小数点数のリストに変換する。
詳細説明:
例えば、`["0.25", "0.5"]` のような文字列のリストを受け取り、
`[0.25, 0.5]` のような浮動小数点数のリストを返す。
入力が `None` の場合は空のリストを返す。
:param str_list: list of str or None, パースする文字列のリスト。
:returns: list of float, パースされた浮動小数点数のリスト。
:rtype: list[float]
"""
if str_list is None:
return []
out = []
for s in str_list:
out.append(float(s))
return out
[ドキュメント]
def main():
"""
VASP体積データの可視化スクリプトのメインエントリポイント。
概要:
コマンドライン引数を解析し、CHGCARファイルの読み込み、データ処理、
および3Dプロットによる可視化（等値面、高密度点群、スライス面）を実行する。
詳細説明:
`argparse` を使用して、入力ファイル、表示モード、等値レベル、ダウンサンプリング、
スライス平面の位置、カラーマップ、透明度、出力ファイル名などの様々なオプションを
受け付ける。読み込まれたCHGCARデータは、必要に応じてダウンサンプリングされ、
指定されたモードでmatplotlibの3Dプロットとして描画される。
結果は画面に表示されるか、画像ファイルとして保存される。
:param: None, コマンドライン引数からパラメータを受け取るため、関数としての引数はない。
:returns: None
"""
ap = argparse.ArgumentParser(description="Visualize VASP volumetric data with isosurfaces and slice planes.")
ap.add_argument("infile")
# 表示モード
ap.add_argument("--mode", choices=["iso", "dots", "both"], default="both",
help="Isosurfaces, dots, or both (isosurface + slices).")
# 等値面
ap.add_argument("--levels", type=float, nargs="+", help="Iso-surface levels (absolute values)")
ap.add_argument("--nlevels", type=int, default=None, help="Auto-generate N iso levels between min/max")
ap.add_argument("--alpha", type=float, default=0.30, help="Alpha for isosurfaces/dots")
ap.add_argument("--edge", default="none")
ap.add_argument("--lw", type=float, default=0.0)
ap.add_argument("--no-legend", action="store_true")
# dots
ap.add_argument("--cutoff", type=float, default=None)
ap.add_argument("--quantile", type=float, default=None)
ap.add_argument("--subsample", type=int, default=1)
ap.add_argument("--max-points", type=int, default=None)
ap.add_argument("--size", type=float, default=0.4)
ap.add_argument("--cmap", default="viridis")
# スライス（分率座標で指定：0～1）
ap.add_argument("--slice-xy", nargs="+", help="Fractions for XY slices (fix fractional z), e.g., 0.25 0.5")
ap.add_argument("--slice-yz", nargs="+", help="Fractions for YZ slices (fix fractional x)")
ap.add_argument("--slice-zx", nargs="+", help="Fractions for ZX slices (fix fractional y)")
ap.add_argument("--slice-cmap", default="viridis")
ap.add_argument("--slice-alpha", type=float, default=0.6)
# 高速化
ap.add_argument("--float32", action="store_true")
ap.add_argument("--ds", type=int, nargs=3, metavar=("DSX","DSY","DSZ"), default=[1,1,1],
help="Downsample steps per axis (e.g., --ds 2 2 2)")
# 表示
ap.add_argument("--ortho", action="store_true")
ap.add_argument("--pad", type=float, default=0.0)
ap.add_argument("--save", default=None)
ap.add_argument("--title", default=None)
args = ap.parse_args()
dtype = np.float32 if args.float32 else np.float64
info, _, F = read_chgcar(args.infile, dtype=dtype)
lattice = (info['lattice'] * info['scale']).astype(dtype, copy=False)
nx, ny, nz = info['grid']
# ダウンサンプリング
dsx, dsy, dsz = args.ds
if dsx > 1 or dsy > 1 or dsz > 1:
F = downsample_volume(F, dsx, dsy, dsz)
spacing_frac = fractional_spacing_after_downsampling(nx, ny, nz, dsx, dsy, dsz)
nx2, ny2, nz2 = F.shape
else:
spacing_frac = (1.0/nx, 1.0/ny, 1.0/nz)
nx2, ny2, nz2 = nx, ny, nz
# 直交グリッド（dots 用と軸範囲決定用）
X, Y, Z = build_cartesian_grid(lattice, nx2, ny2, nz2, dtype=dtype)
# 図
fig = plt.figure(figsize=(9, 7))
ax = fig.add_subplot(111, projection="3d")
if args.ortho:
try:
ax.set_proj_type('ortho')
except Exception:
pass
title = args.title if args.title else info['title']
ax.set_title(title)
# 軸範囲（pad 余白）
minx, maxx = float(X.min()), float(X.max())
miny, maxy = float(Y.min()), float(Y.max())
minz, maxz = float(Z.min()), float(Z.max())
if args.pad != 0.0:
minx -= args.pad; maxx += args.pad
miny -= args.pad; maxy += args.pad
minz -= args.pad; maxz += args.pad
ax.set_xlim([minx, maxx]); ax.set_ylim([miny, maxy]); ax.set_zlim([minz, maxz])
ax.set_box_aspect([1, 1, 1])
ax.set_xlabel("X (Å)"); ax.set_ylabel("Y (Å)"); ax.set_zlabel("Z (Å)")
# 等値面
if args.mode in ("iso", "both"):
vmin, vmax = float(F.min()), float(F.max())
if args.levels:
levels = args.levels
elif args.nlevels:
levels = np.linspace(vmin, vmax, args.nlevels + 2, dtype=float)[1:-1]
else:
levels = np.linspace(vmin, vmax, 5, dtype=float)[1:-1]
cmap = plt.get_cmap('tab10')
colors = [cmap(i % 10) for i in range(len(levels))]
plot_isosurfaces_fractional(
ax, F, lattice, levels,
alpha=args.alpha,
edgecolor=(None if args.edge == 'none' else args.edge),
linewidth=args.lw,
colors=colors,
spacing_frac=spacing_frac,
legend=not args.no_legend if hasattr(args, "no_legend") else True,
antialias=False
)
# dots
if args.mode in ("dots",):
# 値選別
vals = F.ravel()
xs = X.ravel(); ys = Y.ravel(); zs = Z.ravel()
if args.quantile is not None:
thr = np.quantile(vals, args.quantile)
mask = vals >= thr
elif args.cutoff is not None:
mask = vals >= float(args.cutoff)
else:
thr = np.median(vals)
mask = vals >= thr
xs = xs[mask]; ys = ys[mask]; zs = zs[mask]; cs = vals[mask]
if args.subsample > 1:
xs = xs[::args.subsample]; ys = ys[::args.subsample]; zs = zs[::args.subsample]; cs = cs[::args.subsample]
if args.max_points is not None and xs.size > args.max_points:
idx = np.random.default_rng(0).choice(xs.size, size=args.max_points, replace=False)
xs = xs[idx]; ys = ys[idx]; zs = zs[idx]; cs = cs[idx]
norm = plt.Normalize(vmin=float(cs.min()), vmax=float(cs.max()))
tk3d.plot_scatter3d(ax, xs, ys, zs,
minx, maxx, miny, maxy, minz, maxz,
cmap=args.cmap, c=cs, norm=norm, size=args.size, alpha=args.alpha)
if not args.no_legend:
tk3d.make_colorbar(ax, cs, plt.get_cmap(args.cmap), vmin=float(cs.min()), vmax=float(cs.max()), label='Value')
# スライス
fr_xy = parse_fracs(args.slice_xy)
fr_yz = parse_fracs(args.slice_yz)
fr_zx = parse_fracs(args.slice_zx)
for fz in fr_xy:
render_slice_plane(ax, lattice, F, which='xy', frac=fz, cmap=args.slice_cmap, alpha=args.slice_alpha)
for fx in fr_yz:
render_slice_plane(ax, lattice, F, which='yz', frac=fx, cmap=args.slice_cmap, alpha=args.slice_alpha)
for fy in fr_zx:
render_slice_plane(ax, lattice, F, which='zx', frac=fy, cmap=args.slice_cmap, alpha=args.slice_alpha)
if args.save:
plt.tight_layout()
plt.savefig(args.save, dpi=220)
plt.show()
if __name__ == "__main__":
main()