概要
古典分子動力学計算のパッケージソフトウェアであるLAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)により得られたlog.lammps(デフォルト)をmatplotlibにより描画する方法を解説します。
→LAMMPSのリンク
対象ファイル
log.lammpsの注目部分は以下のようになっています。
またこのlogが記載される部分は、温度制御や圧力制御により分断されていることが往々にしてあります。
今回のシミュレーションは単純な系の酸化物ガラスにて実行したものです。
(300 Kで保持したのち、2000 Kまで昇温して保持)
Per MPI rank memory allocation (min/avg/max) = 13.6 | 13.6 | 13.6 Mbytes
Step CPU Temp PotEng KinEng TotEng Enthalpy Press Volume Lx Ly Lz Density
0 0 0 -4176.6949 0 -4176.6949 -4198.9215 -8520.0937 4179.6496 16.108186 16.108186 16.108186 2.3861959
100 0.450147 172.69601 -4181.7735 6.4512636 -4175.3223 -4173.1446 845.82115 4124.8893 16.037528 16.037528 16.037528 2.417874
200 0.896982 204.36575 -4180.706 7.6343242 -4173.0717 -4169.2759 1479.6702 4110.0746 16.018305 16.018305 16.018305 2.4265892
300 1.393903 235.47098 -4179.2291 8.7962967 -4170.4328 -4170.1291 117.43742 4143.0005 16.060966 16.060966 16.060966 2.4073042
Step CPU Temp PotEng KinEng TotEng Enthalpy Press Volume Lx Ly Lz Density
という文字列を発見したら（そうでない場合も動的に読み取ります）、それ以降の行に処理を行うという方針としました。
全体のコード
main.py
def read_lmp_log(filename):
with open(filename, "r") as f:
lines = [line.strip() for line in f.readlines()]
log_data = []
count_log_information = 0
for i, line in enumerate(lines):
if line.startswith("Step"):
log_information = line.split()
count_log_information = len(log_information)
continue
values = line.split()
if len(values) == count_log_information and values:
try:
# Create a list of valid numeric values
valid_row = [
float(value) if is_number(value) else None for value in values
]
if all(val is not None for val in valid_row):
log_data.append(valid_row)
except ValueError as e:
print(f"Skipping line due to conversion error: {line}")
log_dict = {key: [] for key in log_information}
for row in log_data:
for key, value in zip(log_information, row):
log_dict[key].append(value)
return log_dict
def is_number(value):
try:
float(value)
return True
except ValueError:
return False
def visualize_log_data(log_dict, target_property):
import matplotlib.pyplot as plt
if target_property not in log_dict:
print(f"Error: {target_property} not found in the log data.")
return
steps = log_dict["Step"]
target_data = log_dict[target_property]
# Create the plot
plt.figure(figsize=(7, 5))
plt.plot(steps, target_data, label=target_property)
plt.xlabel("Step")
plt.ylabel(target_property)
plt.title(f"{target_property} vs Step")
plt.legend()
plt.tight_layout()
plt.show()
def main():
filename_lmp_log = "log.lammps"
log_dict = read_lmp_log(filename_lmp_log)
target_property = "Density" # 必要により変更してください
visualize_log_data(log_dict, target_property)
if __name__ == "__main__":
main()
それぞれ解説
main関数
def main():
filename_lmp_log = "log.lammps"
log_dict = read_lmp_log(filename_lmp_log)
target_property = "Density" # 必要により変更してください
visualize_log_data(log_dict, target_property)
- read_lmp_log関数によりlog_dict(dict型)を取得
- vidualize_log_data関数により可視化（ステップに対して行います）
read_lmp_log関数
def read_lmp_log(filename):
with open(filename, "r") as f:
lines = [line.strip() for line in f.readlines()]
log_data = []
count_log_information = 0
for i, line in enumerate(lines):
if line.startswith("Step"):
log_information = line.split()
count_log_information = len(log_information)
continue
values = line.split()
if len(values) == count_log_information and values:
try:
# Create a list of valid numeric values
valid_row = [
float(value) if is_number(value) else None for value in values
]
if all(val is not None for val in valid_row):
log_data.append(valid_row)
except ValueError as e:
print(f"Skipping line due to conversion error: {line}")
log_dict = {key: [] for key in log_information}
for row in log_data:
for key, value in zip(log_information, row):
log_dict[key].append(value)
return log_dict
linesをfor文＋enumerate関数（今回indexは用いてないのでその必要はないですが）により処理する際にまずいったんlistに各行を保存していきます。
if len(values) == count_log_information and values:
try:
# Create a list of valid numeric values
valid_row = [
float(value) if is_number(value) else None for value in values
]
if all(val is not None for val in valid_row):
log_data.append(valid_row)
except ValueError as e:
print(f"Skipping line due to conversion error: {line}")
- Step・・・と続く行がcount_log_informationに入っており、その列数とそれに続く文字列の列数が一致したらその行の数値列をvalid_rowさらにはlog_dataに保存していきます。
try-exceptに関しては以下が参考になります。
→リンク
def is_number(value):
try:
float(value)
return True
except ValueError:
return False
- is_number関数を用いている理由は、注目している数値列以外にもlen(values)が一致する行があった際に例外としてはじくためです。
log_dict = {key: [] for key in log_information}
for row in log_data:
for key, value in zip(log_information, row):
log_dict[key].append(value)
- 可視化の際に便利なので、dict型に保存します。
可視化する(visualize_log_data関数)
def visualize_log_data(log_dict, target_property):
import matplotlib.pyplot as plt
if target_property not in log_dict:
print(f"Error: {target_property} not found in the log data.")
return
steps = log_dict["Step"]
target_data = log_dict[target_property]
# Create the plot
plt.figure(figsize=(7, 5))
plt.plot(steps, target_data, label=target_property)
plt.xlabel("Step")
plt.ylabel(target_property)
plt.title(f"{target_property} vs Step")
plt.legend()
plt.tight_layout()
plt.show()
横軸をStep数、縦軸をtarget_property(main関数にして指定)としてプロットします。
今回の例としてはStepに対してDensityをプロットしています。
まとめ
いかがでしたでしょうか。list内包表記や例外処理の練習になりました。
今後は他の第一原理計算パッケージの出力結果に対しても行っていきます。
誤っている部分がある場合は、コメントいただけると嬉しいです。