LAMMPS計算結果の温度やエネルギーなどが書かれたログをMatplotlibでグラフにします。
非GUI環境でMatplotlibを使う方法も解説しています。
ログをNumpyArrayにする
以下のようなコードを書きます。
import numpy as np
with open("log.lammps", "r") as f:
text = f.read()
lines = text.split("\n")
LAMMPS_stat = False
Label = ""
data = []
for line in lines:
if "Loop time of " in line:
LAMMPS_stat = False
if LAMMPS_stat:
data.append(line.split())
if "Step " in line:
LAMMPS_stat = True
Label = line.split()
if "timestep " in line:
line1 = line.split()
dt = float(line1[1])
data = np.array(data, dtype=np.float64)
print(Label)
print(data[:5])仕組みとしてはログの始まりを検出したらLAMMPS_statをTrueにして、終わりを検出したらFalseにして、それがTrueの間のテキストを読んでsplitしてdata入れる、という感じです。dataの中身が空っぽになるなど上手く動かない場合、18行目の”Step ”の文字列が一致していないか別の箇所で反応している可能性があります。前後のスペースを追加したり”Step Temp”のように次のキーワードも一緒にするなどして対応してください。
dataは浮動小数点型のNumpyArrayに変換して、Labelと一緒に5行だけprintしています。
また、グラフの横軸で必要になるtimestepの値も取得しています。
Matplotlibでグラフ化
thermo_style custom step temp press etotal density pe ke で出力されるdensityをプロットしています。
コードはこんな感じです。
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
with open("log.lammps", "r") as f:
text = f.read()
lines = text.split("\n")
LAMMPS_stat = False
Label = ""
data = []
for line in lines:
if "Loop time of " in line:
LAMMPS_stat = False
if LAMMPS_stat:
data.append(line.split())
if "Step " in line:
LAMMPS_stat = True
Label = line.split()
if "timestep " in line:
line1 = line.split()
dt = float(line1[1])
data = np.array(data, dtype=np.float64)
x = data[:,0] * dt * 0.001
y = data[:,4]
i1 = int(len(x)/2)
i2 = len(x)
print("Average ", np.mean(y[i1:i2]))
plt.plot(x, y)
plt.xlabel("Time [ps]")
plt.ylabel("Density [g/cm3]")
plt.savefig("Density.png")最初の matplotlib.use("Agg") によって、非GUI環境でもグラフが描けるようにしています。
サーバにSSHでつないだりするときは必ずこうしています。
密度の平均値を計算の後半部分で出力していますが、不要であれば削除してください。
このコードでは以下のようなグラフが描けます。
こちらは水の密度を計算したものです。
詳細はこちら
まとめ
LAMMPSのログ出力をPythonのNumpyに読み込んでMatplotlibでグラフ化しました。
また非GUI環境でMatplotlibを使う方法も解説しました。
コメント