【Python】LAMMPSのlog.lammpsから、logデータを抽出して可視化

# 概要

古典分子動力学計算のパッケージソフトウェアであるLAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)により得られたlog.lammps(デフォルト)をmatplotlibにより描画する方法を解説します。
→[LAMMPSのリンク](https://docs.lammps.org/Manual.html)

# 対象ファイル
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関数
```python
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関数
```python
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に各行を保存していきます。

```python
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に関しては以下が参考になります。
→[リンク](https://qiita.com/suipy/items/9d02e197df813f3c2fab)

```python
def is_number(value):
    try:
        float(value)
        return True
    except ValueError:
        return False
```
- is_number関数を用いている理由は、注目している数値列以外にもlen(values)が一致する行があった際に例外としてはじくためです。

```python
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関数)

```python
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()
```

![forqiita.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3956079/1495fb6c-9ed3-4c74-bced-5d4cdba1fd6a.png)


横軸をStep数、縦軸をtarget_property(main関数にして指定)としてプロットします。
今回の例としてはStepに対してDensityをプロットしています。


# まとめ

いかがでしたでしょうか。list内包表記や例外処理の練習になりました。
今後は他の第一原理計算パッケージの出力結果に対しても行っていきます。
誤っている部分がある場合は、コメントいただけると嬉しいです。