Molecule and Bond Counts from Reactive MD

Molecular Dynamics Trajectory analysis ReaxFF Technical

Run a short reactive molecular dynamics simulation with ReaxFF and tabulate the time evolution of molecule and bond-type counts from the generated trajectory.

Downloads: Notebook | Script ?

Requires: AMS2027 or later

Related examples
Related documentation
../_images/bond_counts_rkf_7_0_1b59eb01.png

Molecule and Bond Counts from Reactive MD

This example performs molecule and bond counting from a reactive MD trajectory. It first runs a short ReaxFF simulation, then plots molecule populations versus time and bond counts per frame from the generated ams.rkf.

The MD run requires a ReaxFF license. If you already have a reactive MD trajectory with molecule data, you can skip the MD cell and set ams_rkf_path to your own file.

Initial Imports

from collections import Counter
from pathlib import Path
from typing import Dict, List, Literal, Optional, Sequence, Set, Tuple

import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from scm.plams import AMSJob, Molecule, Settings, Trajectory, from_smiles, packmol, plot_molecule_counts, view

Run a short reactive MD simulation

This produces an ams.rkf trajectory with molecule and bond information for analysis.

o2 = from_smiles("O=O")
h2 = from_smiles("[HH]")
mixture = packmol(molecules=[o2, h2], n_molecules=[4, 4], density=1.0)
view(mixture, direction="tilt_z", width=300, height=300)

md_settings = Settings()
md_settings.input.reaxff.forcefield = "CHO.ff"
md_settings.input.ams.task = "MolecularDynamics"
md_settings.input.ams.MolecularDynamics.NSteps = 3000
md_settings.input.ams.MolecularDynamics.TimeStep = 0.5
md_settings.input.ams.MolecularDynamics.InitialVelocities.Temperature = 3500
md_settings.input.ams.MolecularDynamics.Trajectory.WriteMolecules = "Yes"
md_settings.input.ams.MolecularDynamics.Trajectory.WriteBonds = "Yes"
md_settings.runscript.nproc = 1
md_settings.runscript.preamble_lines = ["export OMP_NUM_THREADS=1"]

md_job = AMSJob(name="reactive_md_toy", settings=md_settings, molecule=mixture)
md_results = md_job.run();

[10:35:05] WARNING: not removing hydrogen atom without neighbors
[10:35:05] WARNING: not removing hydrogen atom without neighbors


[13.08|10:35:06] JOB reactive_md_toy STARTED
[13.08|10:35:06] JOB reactive_md_toy RUNNING
[13.08|10:35:09] JOB reactive_md_toy FINISHED
[13.08|10:35:09] JOB reactive_md_toy SUCCESSFUL

Molecule populations versus time

plot_molecule_counts reads the molecule analysis from ams.rkf and plots the population of every detected molecular formula against time. Pass species=["H2", "O2"] to plot only selected formulas, or set x_axis="frame" to use trajectory frames instead.

ax = plot_molecule_counts(md_job, time_unit="ps")
ax.figure.tight_layout()
ax;
image generated from notebook

Extract the molecule populations

get_molecule_count_history returns the same data as a dictionary containing frame indices, times, and one array per requested formula. Here it is converted to a table for further analysis or export.

import pandas as pd

molecule_population_history = md_job.results.get_molecule_count_history(
    species=["H2", "O2", "H2O"], time_unit="ps"
)
molecule_populations = pd.DataFrame(molecule_population_history).rename(columns={"time": "time_ps"})
molecule_populations

frame

time_ps

H2

O2

H2O

0

1

0.00

4

4

0

1

2

0.05

3

2

0

2

3

0.10

2

2

1

3

4

0.15

2

1

0

4

5

0.20

1

1

1

5

6

0.25

1

1

0

6

7

0.30

0

1

2

7

8

0.35

0

1

2

8

9

0.40

0

1

1

9

10

0.45

0

1

3

10

11

0.50

0

2

1

11

12

0.55

0

1

3

12

13

0.60

0

1

2

13

14

0.65

0

0

2

14

15

0.70

0

1

3

15

16

0.75

0

1

3

16

17

0.80

0

1

2

17

18

0.85

0

2

3

18

19

0.90

0

1

2

19

20

0.95

0

1

2

20

21

1.00

0

1

3

21

22

1.05

0

1

2

22

23

1.10

0

1

3

23

24

1.15

0

1

3

24

25

1.20

0

1

2

25

26

1.25

0

1

0

26

27

1.30

0

1

2

27

28

1.35

0

1

2

28

29

1.40

0

0

0

29

30

1.45

0

2

2

30

31

1.50

0

1

2

To focus on selected species, pass their molecular formulas with the species argument.

ax = plot_molecule_counts(md_job, species=["H2", "O2", "H2O"], time_unit="ps")
ax.figure.tight_layout()
ax;
image generated from notebook

Bond counts per frame

Count bond types per frame, using bond order rounded to either integer or half-integer values.

import pandas as pd

trajectory = Trajectory(md_job.results.rkfpath())
n_frames = len(trajectory)


def bo_to_symbol(bo):
    s = "---"
    if bo == 0.5:
        s = "--"
    if bo == 1.0:
        s = "-"
    if bo == 1.5:
        s = "=="
    if bo == 2.0:
        s = "="
    if bo == 2.5:
        s = "≡≡"
    if bo == 3.0:
        s = "≡"
    return s


all_bonds = {}
for frame_index, mol in enumerate(trajectory):
    frame_counts: Counter = Counter()
    seen_pairs: Set[Tuple[int, int]] = set()

    for bond in mol.bonds:
        pair = tuple(sorted((mol.index(bond.atom1), mol.index(bond.atom2))))
        if pair in seen_pairs:
            continue
        seen_pairs.add(pair)

        symbols = sorted((bond.atom1.symbol, bond.atom2.symbol))
        bo = float(round(bond.order))
        label = f"{symbols[0]}{bo_to_symbol(bo)}{symbols[1]} (BO {bo:.1f})"
        frame_counts[label] += 1

    for label, count in frame_counts.items():
        if label not in all_bonds:
            all_bonds[label] = [0] * n_frames
        all_bonds[label][frame_index] = int(count)

bond_names = sorted(all_bonds.keys())
bond_counts = []
for frame_index in range(n_frames):
    counts = {"Frame": frame_index + 1}
    for name in bond_names:
        counts[name] = all_bonds[name][frame_index]
    bond_counts.append(counts)

bond_counts_df = pd.DataFrame(bond_counts)
bond_counts_df

Frame

H—H (BO 0.0)

H—O (BO 0.0)

H-H (BO 1.0)

H-O (BO 1.0)

O-O (BO 1.0)

O=O (BO 2.0)

0

1

0

0

4

0

0

4

1

2

0

0

3

2

3

1

2

3

0

0

2

4

2

1

3

4

0

0

2

5

1

1

4

5

0

0

1

6

1

1

5

6

0

1

1

6

1

1

6

7

0

0

0

8

1

1

7

8

0

0

1

7

1

1

8

9

0

0

0

8

1

1

9

10

0

0

0

8

2

0

10

11

1

0

0

7

0

2

11

12

0

0

0

8

2

0

12

13

1

0

0

7

1

1

13

14

0

0

0

8

1

1

14

15

0

1

0

8

2

0

15

16

0

2

0

8

1

1

16

17

0

0

0

10

1

1

17

18

0

1

0

7

0

2

18

19

0

0

0

8

1

0

19

20

0

2

0

7

2

0

20

21

0

1

0

8

0

1

21

22

0

1

0

9

0

1

22

23

0

0

0

9

0

1

23

24

0

1

0

8

1

0

24

25

0

2

0

9

0

1

25

26

0

0

0

11

0

1

26

27

0

0

0

8

2

0

27

28

0

0

0

8

1

1

28

29

0

0

0

10

3

0

29

30

0

0

0

8

0

2

30

31

0

0

0

7

2

0

fig, ax = plt.subplots(figsize=(8, 4))

x = bond_counts_df["Frame"]

for i, name in enumerate(bond_names):
    y = bond_counts_df[name]
    ax.plot(x, y, label=name)

ax.set_xlabel("Frame")
ax.set_ylabel("Bond Count")
ax.legend(loc="best", fontsize=8)
fig.tight_layout()
ax;
image generated from notebook

See also

Python Script

#!/usr/bin/env python
# coding: utf-8

# ## Molecule and Bond Counts from Reactive MD
# 
# This example performs molecule and bond counting from a reactive MD trajectory. It first runs a short ReaxFF simulation, then plots molecule populations versus time and bond counts per frame from the generated `ams.rkf`.

# The MD run requires a ReaxFF license. If you already have a reactive MD trajectory with molecule data, you can skip the MD cell and set `ams_rkf_path` to your own file.
# 

# ### Initial Imports

from collections import Counter
from pathlib import Path
from typing import Dict, List, Literal, Optional, Sequence, Set, Tuple

import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from scm.plams import AMSJob, Molecule, Settings, Trajectory, from_smiles, packmol, plot_molecule_counts, view


# ### Run a short reactive MD simulation
# 
# This produces an `ams.rkf` trajectory with molecule and bond information for analysis.
# 

o2 = from_smiles("O=O")
h2 = from_smiles("[HH]")
mixture = packmol(molecules=[o2, h2], n_molecules=[4, 4], density=1.0)
view(mixture, direction="tilt_z", width=300, height=300)

md_settings = Settings()
md_settings.input.reaxff.forcefield = "CHO.ff"
md_settings.input.ams.task = "MolecularDynamics"
md_settings.input.ams.MolecularDynamics.NSteps = 3000
md_settings.input.ams.MolecularDynamics.TimeStep = 0.5
md_settings.input.ams.MolecularDynamics.InitialVelocities.Temperature = 3500
md_settings.input.ams.MolecularDynamics.Trajectory.WriteMolecules = "Yes"
md_settings.input.ams.MolecularDynamics.Trajectory.WriteBonds = "Yes"
md_settings.runscript.nproc = 1
md_settings.runscript.preamble_lines = ["export OMP_NUM_THREADS=1"]

md_job = AMSJob(name="reactive_md_toy", settings=md_settings, molecule=mixture)
md_results = md_job.run();


# ### Molecule populations versus time
# 
# `plot_molecule_counts` reads the molecule analysis from `ams.rkf` and plots the population of every detected molecular formula against time. Pass `species=["H2", "O2"]` to plot only selected formulas, or set `x_axis="frame"` to use trajectory frames instead.
# 

ax = plot_molecule_counts(md_job, time_unit="ps")
ax.figure.tight_layout()
ax;
ax.figure.savefig("picture1.png")


# ### Extract the molecule populations
# 
# `get_molecule_count_history` returns the same data as a dictionary containing frame indices, times, and one array per requested formula. Here it is converted to a table for further analysis or export.
# 

import pandas as pd

molecule_population_history = md_job.results.get_molecule_count_history(
    species=["H2", "O2", "H2O"], time_unit="ps"
)
molecule_populations = pd.DataFrame(molecule_population_history).rename(columns={"time": "time_ps"})
molecule_populations


# To focus on selected species, pass their molecular formulas with the `species` argument.
# 

ax = plot_molecule_counts(md_job, species=["H2", "O2", "H2O"], time_unit="ps")
ax.figure.tight_layout()
ax;
ax.figure.savefig("picture2.png")


# ## Bond counts per frame
# 
# Count bond types per frame, using bond order rounded to either integer or half-integer values.
# 

import pandas as pd

trajectory = Trajectory(md_job.results.rkfpath())
n_frames = len(trajectory)


def bo_to_symbol(bo):
    s = "---"
    if bo == 0.5:
        s = "--"
    if bo == 1.0:
        s = "-"
    if bo == 1.5:
        s = "=="
    if bo == 2.0:
        s = "="
    if bo == 2.5:
        s = "≡≡"
    if bo == 3.0:
        s = "≡"
    return s


all_bonds = {}
for frame_index, mol in enumerate(trajectory):
    frame_counts: Counter = Counter()
    seen_pairs: Set[Tuple[int, int]] = set()

    for bond in mol.bonds:
        pair = tuple(sorted((mol.index(bond.atom1), mol.index(bond.atom2))))
        if pair in seen_pairs:
            continue
        seen_pairs.add(pair)

        symbols = sorted((bond.atom1.symbol, bond.atom2.symbol))
        bo = float(round(bond.order))
        label = f"{symbols[0]}{bo_to_symbol(bo)}{symbols[1]} (BO {bo:.1f})"
        frame_counts[label] += 1

    for label, count in frame_counts.items():
        if label not in all_bonds:
            all_bonds[label] = [0] * n_frames
        all_bonds[label][frame_index] = int(count)

bond_names = sorted(all_bonds.keys())
bond_counts = []
for frame_index in range(n_frames):
    counts = {"Frame": frame_index + 1}
    for name in bond_names:
        counts[name] = all_bonds[name][frame_index]
    bond_counts.append(counts)

bond_counts_df = pd.DataFrame(bond_counts)
print(bond_counts_df)


fig, ax = plt.subplots(figsize=(8, 4))

x = bond_counts_df["Frame"]

for i, name in enumerate(bond_names):
    y = bond_counts_df[name]
    ax.plot(x, y, label=name)

ax.set_xlabel("Frame")
ax.set_ylabel("Bond Count")
ax.legend(loc="best", fontsize=8)
fig.tight_layout()
ax;
ax.figure.savefig("picture3.png")