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
-------------

.. code-block:: python

   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:

.. code-block:: python

   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:

.. code-block:: python

   {"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:

.. code-block:: python

   {"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:

.. code-block:: python

   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:

.. code-block:: python

   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.

.. code-block:: python

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