Analysis workflow

The Sample class coordinates post-processing of trajectories and bond files. You can configure multiple samplers and execute them in serial or parallel.

The easiest way to analyse a simulation is to run the ana.py script generated by the Simulate class. All you need to do is to set the atoms and bonds you want to sample. This creates .obj files that storing the sampled data.

Available sampler APIs

  • add_molecule_structure_sampling

  • add_charge_sampling

  • add_angle_sampling

  • add_density_sampling

  • add_bond_density_sampling

  • add_bond_length_sampling

  • add_rdf_sampling

  • add_reaction_sampling

Example setup

from porereax import Sample
atom_lib = {'H': 1, 'O': 2, 'Si': 3}
masses = {'Si': 28.086, 'O': 15.9994, 'H': 2.016}
sampler = Sample(
    atom_lib=atom_lib,
    masses=masses,
    trajectory_file="run_0.lammpstrj",
    bond_file="run_0.bonds",
)
sampler.add_molecule_structure_sampling(
    name_out="molecule_structure"
)
sampler.add_charge_sampling(
    name_out="charge",
    atoms=[{"atom": "O", "bonds": ["H", "H"]}]
)
sampler.add_rdf_sampling(
    name_out="rdf",
    pairs=[({"atom": "O", "bonds": ["Si", "H"]}, {"atom": "O", "bonds": ["H", "H"]})],
)
sampler.sample(is_parallel=True)

This generates the following output files:

  • molecule_structure.obj

  • charge.obj

  • rdf.obj

Execution notes

  • Parallel execution uses multiprocessing. If your system only supports fork it is not possible to run with the ovito module loaded (or other modules that use multithreaded engine)

  • Depending on environment and trajectory size, serial mode may be preferable:

sampler.sample(is_parallel=False)

Setting up atoms and bonds

Atoms and bonds are defined as dictionaries. For atoms, the dictionary must contain a key atom with an atom identifier. Optionally, a key bonds can be added to specify the bonding environment of the atom as a list of atom identifiers. If no bonding environment is specified, all atoms of the given type will be sampled. For example, to sample the charge of oxygen atoms bonded to two hydrogen atoms, the following dictionary can be used:

{"atom": "O", "bonds": ["H", "H"]}

Bonds must contain a bond key with two atom identifiers separated by a - character. Optionally, a bonds_A and bonds_B key can be added to specify the additional bonding environment of the first or second atom in the bond. For example, to sample the bond length of O-H bonds in a water molecule, the following dictionary can be used:

{"bond": "O-H", "bonds_A": ["H"], "bonds_B": []}

Visualise the results

All results are stored in .obj files. To get the raw data, you can load the files with the load_object function returning a dictionary. For example, to load the charge data:

from porereax import load_object
charge_data = load_object("charge.obj")

The sampled atoms and bonds are maped to unique identifiers as keys in the dictionary. All identifiers can be retrieved as a list with the get_identifiers function. To get only the data for a specific atom or bond, you can use the get_data function. For example, to get the charge data for oxygen atoms bonded to two hydrogen atoms:

from porereax import get_data
charge_O_H2 = get_data("charge.obj", "O(H+H)")

The easiest way to visualize the results, stored in .obj files, is to use the plot(...) function. It takes the path of the output file and returns a matplotlib figure and axes object. The figure can be displayed or saved to a file.

from porereax import plot
import matplotlib.pyplot as plt
fig, ax = plot("charge.obj")
plt.show()