【Python】LAMMPSの座標ファイル(dump)から座標を取得

Last updated at 2025-02-25Posted at 2025-02-23

概要

古典分子動力学計算のパッケージソフトウェアであるLAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)により得られた計算結果を取得し、構造解析を行うということがあります。その際にPythonによりファイル内容を読み取る方法を解説します。
LAMMPSのリンク

ファイルの概要

対象ファイルは以下のdumpコマンドで作成されたテキストファイルであることを想定しています。

dump 1 all custom 100000 lammps_dump_file.txt id type x y z

これは100000ステップごとにid type x y zをlammps_dump_file.txtに出力する例です。

中身は例えば以下のようになっています。

ITEM: TIMESTEP
1300000
ITEM: NUMBER OF ATOMS
290
ITEM: BOX BOUNDS pp pp pp
1.2765723376990704e-01 1.6235842766230988e+01
1.2765723376990704e-01 1.6235842766230988e+01
1.2765723376990704e-01 1.6235842766230988e+01
ITEM: ATOMS id type x y z
1 1 5.10984 5.01571 6.37821
2 1 1.71332 8.87056 5.12845
3 1 12.6497 0.996619 5.99052
4 1 2.29157 0.38054 9.74531
5 1 11.821 10.846 1.15356
6 1 8.33773 5.39669 3.05033
7 1 9.96088 6.39124 0.928865

1300000はそのフレームのステップ
290は粒子数
セルは、立方体セルを想定しています。
またITEM: ATOMS id type x y z以降に各原子座標がプロットされています。

全体のコード

main.py

def read_lmp_dump(filename):

    with open(filename, "r") as f:
        lines = [line.strip() for line in f.readlines()]

    step_everyframes = []
    box_length_everyframes = []
    coordination_everyframes = []

    for i, line in enumerate(lines):

        if line == "ITEM: TIMESTEP":
            next_line = lines[i + 1]
            step_frame = int(next_line)
            step_everyframes.append(step_frame)

        if line == "ITEM: NUMBER OF ATOMS":
            next_line = lines[i + 1]
            N_tot = int(next_line)

        if line == "ITEM: BOX BOUNDS pp pp pp":
            next_line = lines[i + 1]
            box_min, box_max = next_line.split()
            box_length = float(box_max) - float(box_min)
            box_length_everyframes.append(box_length)

        if line == "ITEM: ATOMS id type x y z":

            coordinate_frame = []

            for atom_id in range(N_tot):
                next_line = lines[i + atom_id + 1]
                coordinate_line = next_line.split()
                coordinate = [float(coord) for coord in coordinate_line[-3:]]
                coordinate_frame.append(coordinate)

            coordination_everyframes.append(coordinate_frame)

    return N_tot, step_everyframes, box_length_everyframes, coordination_everyframes


def main():
    filename_lmp_dump = "lammps_dump_file.txt"
    N_tot, step_everyframes, box_length_everyframes, coordination_everyframes = (
        read_lmp_dump(filename_lmp_dump)
    )

    print(f"Number of atoms:{N_tot}")
    print(f"Step:{step_everyframes}")
    print(f"box_length:{box_length_everyframes}")
    print(f"coordinations:{coordination_everyframes}")


if __name__ == "__main__":
    main()

それぞれ解説

メイン関数

def main():
    filename_lmp_dump = "lammps_dump_file.txt"
    N_tot, step_everyframes, box_length_everyframes, coordination_everyframes = (
        read_lmp_dump(filename_lmp_dump)
    )

    print(f"Number of atoms:{N_tot}")
    print(f"Step:{step_everyframes}")
    print(f"box_length:{box_length_everyframes}")
    print(f"coordinations:{coordination_everyframes}")

先ほどのdumpコマンドにより出力されたファイル名をread_lmp_dump関数に渡し、その結果をf文字列により出力しています。

read_lmp_dump関数

def read_lmp_dump(filename):

    with open(filename, "r") as f:
        lines = [line.strip() for line in f.readlines()]

    step_everyframes = []
    box_length_everyframes = []
    coordination_everyframes = []

    for i, line in enumerate(lines):

        if line == "ITEM: TIMESTEP":
            next_line = lines[i + 1]
            step_frame = int(next_line)
            step_everyframes.append(step_frame)

        if line == "ITEM: NUMBER OF ATOMS":
            next_line = lines[i + 1]
            N_tot = int(next_line)

        if line == "ITEM: BOX BOUNDS pp pp pp":
            next_line = lines[i + 1]
            box_min, box_max = next_line.split()
            box_length = float(box_max) - float(box_min)
            box_length_everyframes.append(box_length)

        if line == "ITEM: ATOMS id type x y z":

            coordinate_frame = []

            for atom_id in range(N_tot):
                next_line = lines[i + atom_id + 1]
                coordinate_line = next_line.split()
                coordinate = [float(coord) for coord in coordinate_line[-3:]]
                coordinate_frame.append(coordinate)

            coordination_everyframes.append(coordinate_frame)

    return N_tot, step_everyframes, box_length_everyframes, coordination_everyframes

txtファイルの文字列を読み取るコードは人によってさまざまあるかと思うのですが、私は基本的にwith openでファイルを開き、読み取ったlines(list)の各要素、lineに対して操作を行うという方法がわかりやすいのでそちらを用いています。

余談ですが、C++に比べPythonはファイル操作がかなりしやすいと感じています。

 with open(filename, "r") as f:
        lines = [line.strip() for line in f.readlines()]

関数の引数で受け取ったfilenameをopenにより開き、リスト内包表記により全行の改行文字\nや余計な空白を消去しています。

if line == "ITEM: TIMESTEP":
            next_line = lines[i + 1]
            step_frame = int(next_line)
            step_everyframes.append(step_frame)

ステップを読み取り、step_everyframesにappendしています。

if line == "ITEM: NUMBER OF ATOMS":
            next_line = lines[i + 1]
            N_tot = int(next_line)

ステップ同様に粒子数も拾っておきます。

if line == "ITEM: BOX BOUNDS pp pp pp":
            next_line = lines[i + 1]
            box_min, box_max = next_line.split()
            box_length = float(box_max) - float(box_min)
            box_length_everyframes.append(box_length)

split()メソッドを用いて、空白により分割しそれぞれ変数に格納しています。立方体セルでない場合は、このブロックは変更する必要があります。

if line == "ITEM: ATOMS id type x y z":

            coordinate_frame = []

            for atom_id in range(N_tot):
                next_line = lines[i + atom_id + 1]
                coordinate_line = next_line.split()
                coordinate = [float(coord) for coord in coordinate_line[-3:]]
                coordinate_frame.append(coordinate)

            coordination_everyframes.append(coordinate_frame)

座標ブロックでは後半3列がそれぞれx, y, zに対応しているので、それをsplit()で分割し後半3列をスライスにより取得しています。

出力結果

Number of atoms:290
Step:[1300000, 1400000, 1500000, 1600000, 1700000, 1800000, 1900000, 2000000, 2100000, 2200000, 2300000, 2400000, 2500000, 2600000, 2700000, 2800000, 2900000, 3000000, 3100000, 3200000, 3300000, 3400000, 3500000, 3600000, 3700000, 3800000, 3900000, 4000000, 4100000, 4200000]
box_length:[16.108185532461082, 16.156499923147596, 16.180855058162663, 16.197134857049537, 16.212271293043283, 16.162278978628663, 16.184759405573097, 16.17435708910438, 16.15618711119981, 16.190179930302918, 16.10362307788567, 16.128279514667735, 16.16309571670343, 16.183055162790104, 16.107036232437366, 16.130639434668485, 16.195124837809733, 16.158793448537683, 16.11934551893416, 16.113227470838982, 16.1998050002098, 16.09434822934957, 16.20957272411758, 16.137734029980265, 16.128412539852846, 16.161546251437773, 16.203477174478948, 16.13786947214681, 16.16112588415096, 16.103949856772253]
coordinations:[[[5.10984, 5.01571, 6.37821], [1.71332, 8.87056, 5.12845], [12.6497, 0.996619, 5.99052], [2.29157, 0.38054, 9.74531], [11.821, 10.846, 1.15356], [8.33773, 5.39669, 3.05033], [9.96088, 6.39124, 0.928865], [13.7562, 13.3305, 1.98335], [0.729708, 4.34297, 1.08828], [11.5645, 1.93099, 11.869], [10.1722, 7.06451, 5.15459], [0.764709, 13.7475, 10.4021], [2.00048, 4.67229, 14.2882], [1.47787, 7.70778, 14.3221], [2.73055, 11.4671, 9.41419], [8.4903, 11.9397, 5.50156], ......

各フレームの情報を抽出できています。

まとめ

LAMMPSの座標ファイルから各フレームの情報を取得しました。こちらの情報を用いて、構造解析を行うことができます。元素種の情報を取得したい場合は、座標データを取得する部分に追加で書き加えてください。
また今後、自分の勉強のために計算化学やデータサイエンスに関する情報を発信していきます。
よろしくお願いいたします。

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up