Full tutorial#

In this tutorial, we want to demonstrate the full workflow of shnitsel tools from import to final analysis and output. Due to size constraints in git, we will use an already compressed shnitsel-file input for Retinal as a basis. Please refer to the ./0_general_io tutorial notebook for information on how to read arbitrary input.

[1]:

import pandas as pd
import seaborn as sns

import shnitsel as st
import shnitsel.xarray

1) Import and preparation of data#

1.2) Optional: Converting the data#

For processing or analysis, you may prefer to work with a single xarray.Dataset. Shnitsel tree formats support conversion to either a stacked (i.e. transformed into a series of frames) or layered format (i.e. a new coordinate trajectory to index the various trajectories along) via two accessor properties:

[4]:

# Get a `stacked` representation of the data in the tree
ds_stacked = truncated_dt.as_stacked
# Get a `layered` representation of the data in the tree
ds_layered = truncated_dt.as_layered

/home/tpadmin/git/shnitsel-tools-official/shnitsel/data/traj_combiner_methods.py:491: FutureWarning: In a future version of xarray the default value for compat will change from compat='equals' to compat='override'. This change will result in the following ValueError: Cannot specify both coords='different' and compat='override'. The recommendation is to set compat explicitly for this case.
  frames = xr.concat(
/home/tpadmin/git/shnitsel-tools-official/shnitsel/data/traj_combiner_methods.py:735: FutureWarning: In a future version of xarray the default value for coords will change from coords='different' to coords='minimal'. This is likely to lead to different results when multiple datasets have matching variables with overlapping values. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set coords explicitly.
  layers = xr.concat(

2) Selection of relevant states and geometrical features#

2.1) Geometrical features#

We may not wish to perform the analysis of principal components or other statistics on the full structure. For example, we may wish to restrict our analysis to a subset. Shnitsel-tools provides a rich module for the selection of these geometric features in shnitsel.filtering. The StructureSelection class helps with the control of relevant structures.

First, you need to initialize the Structure selection:

[5]:

from shnitsel.filtering import StructureSelection

# Create a selection of all features
structure_all = StructureSelection.init_from_dataset(ds_stacked, ['atoms', 'bonds', 'angles', 'dihedrals', 'pyramids'])

# Create a selection of only bond lengths and dihedrals
structure_bd = StructureSelection.init_from_dataset(ds_stacked, ['bonds', 'dihedrals'])
structure_bd.draw()

[5]:

../_images/tutorial_symlinks_0_full_tutorial_9_0.svg

Generally, the initial selection includes all instances of the features in the array following the dataset. In the first case all features supported by the selection class. In the second case only bond pairs and dihedrals.

If you ever change your mind, you can add a specific feature back with the .select_all() method or (without a parameter) have all features activated on the selection:

[6]:

# Add `atoms`, i.e. positions back into the selection
structure_abd = structure_bd.select_all('atoms')
# Add all features back
structure_all2 = structure_bd.select_all()
structure_abd.draw()

[6]:

../_images/tutorial_symlinks_0_full_tutorial_11_0.svg

You can select specific features either with tuples representing them or using SMARTS strings. Additionally, you have the option to select all BATS (positions, bonds, angles, torsions and pyramidalization centers) at the same time with .select() (or .select_bats()) or modify the selection of individual features using .select_atoms(), .select_bonds(), .select_angles(), .select_dihedrals(), and .select_pyramids() respectively.

[7]:

# Select all bonds from a C or N atom to an H atom
structure_all.select_bonds('[#6,#7]~[#1]').draw()

[7]:

../_images/tutorial_symlinks_0_full_tutorial_13_0.svg

[8]:

# Pick one specific dihedral:
structure_all.select_dihedrals((6,5,7,8)).draw('dihedrals')

[8]:

../_images/tutorial_symlinks_0_full_tutorial_14_0.svg

You can combine multiple selections via simple set operations, i.e.

+ or | to add two selections together
- to remove the second selection from the first
& to obtain the intersection between two selections
~ to obtain an inverted selection.

For example, if we want to select the two central dihedrals in retinal to calculate HOOP (hydrogen out of plane):

[9]:

# Only keep the selection of dihedrals
dih_6578 = structure_all.select_dihedrals((6,5,7,8)).only('dihedrals')
dih_CCCC = structure_all.select_dihedrals("[#6][#6]=[#6][#6]").only('dihedrals')

combined = dih_6578 | dih_CCCC
combined.draw('dihedrals')

[9]:

../_images/tutorial_symlinks_0_full_tutorial_16_0.svg

We will now also present a quite visually pleasing way of combining our selection capabilities with the tools to calculate certain derived observables like distances and use that to verify that our filtering was successful.

E.g, our sanity_check() by default enforces a C-H bondlength of at most 2 angstroms. Let us get an overview of C-H bond lengths in the system:

[10]:

import rdkit
from shnitsel.geo.geocalc import get_distances
from rdkit.Chem import Mol
from shnitsel.bridges import default_mol
from shnitsel.units.conversion import convert_length

# Get all lengths of C-H bonds by creating an on-the-go structural selection using SMARTS strings directly
c_h_distances = get_distances(ds_stacked, structure_selection="[#6][#1]")

# Get a mol object for our dataset:
draw_mol = default_mol(ds_stacked)

# Calculate proportions per bond and annotate the molecule.
for bond_descriptor in c_h_distances.feature_indices.values:
    bond_index = structure_all.bond_descriptor_to_mol_index(
        bond_descriptor
    )
    dist = convert_length(c_h_distances.sel(feature_indices=bond_descriptor), "angstrom")

    # Assign the label to the bond
    draw_mol.GetBondWithIdx(bond_index).SetProp(
        'bondNote', f"<l>={dist.min().item():.2f}A"
    )
draw_mol

[10]:

../_images/tutorial_symlinks_0_full_tutorial_18_0.png

So on average our C-H bonds are around 0.9A long. But let us check that all of them are below 2A:

[11]:

import rdkit
from shnitsel.geo.geocalc import get_distances
from rdkit.Chem import Mol
from shnitsel.bridges import default_mol
from shnitsel.units.conversion import convert_length

# Get all lengths of C-H bonds by creating an on-the-go structural selection using SMARTS strings directly
c_h_distances = get_distances(ds_stacked, structure_selection="[#6][#1]")

# Get a mol object for our dataset:
draw_mol = default_mol(ds_stacked)

# Calculate proportions per bond and annotate the molecule.
for bond_descriptor in c_h_distances.feature_indices.values:
    bond_index = structure_all.bond_descriptor_to_mol_index(
        bond_descriptor
    )
    dist = convert_length(c_h_distances.sel(feature_indices=bond_descriptor), "angstrom")

    # Get ratio of bonds below 2 A
    proportion =  ((dist < 2).sum() / dist.sizes['frame']).item()

    # Assign the label to the bond
    draw_mol.GetBondWithIdx(bond_index).SetProp(
        'bondNote', f"{proportion * 100:.0f}%"
    )
draw_mol

[11]:

../_images/tutorial_symlinks_0_full_tutorial_20_0.png

And indeed, all of our bonds are within the constraints, proving that our filtration has worked.

3) Structural analysis#

3.1) General Principal Component Analysis (PCA)#

Shnitsel tools has extensive support for PCA analysis combined with the strong selection/filtration features presented above. Let us demonstrate this on the example of retinal and try and identify the most significant contributing features of the retinal energy landscape.

First, we will perform a PCA to identify key features of the Retinal and investigate clustering behavior. The following box demonstrates multiple options for doing so, using the biplot_kde() function that plots the PCA results, highlights clusters based on a KDE density estimation and plots the most relevant features of the PCAs as visualizations on the molecular structure:

[12]:

from shnitsel.analyze.pca import PCA
from shnitsel.vis.plot.kde import biplot_kde

# We can create the plot with a feature selection and make it calculate the PCA directly:
# biplot_kde(ds_stacked, 6,5,7,8, structure_selection=structure_all.only(['bonds', 'dihedrals', 'angles']), scatter_color_property='geo')

# Or we can calculate the PCA ourselves and provide it to the biplot function:
pca_res = PCA(ds_stacked, structure_selection=structure_all.only(['bonds', 'dihedrals', 'angles']))

# We can plot the PCA analysis with coloring using the `time` coordinate to see how trajectories evolve over time in PCA space:
# biplot_kde(ds_stacked, pca_data=pca_res, scatter_color_property='time')

# We can highlight the dihedral (6,5,7,8) as the coloring property to differentiate the different isomers
# biplot_kde(ds_stacked, 6,5,7,8, pca_data=pca_res, scatter_color_property='geo')

# Or we can use the distance between H atoms 4 and 10 as a coloring parameter to indicate isomerization
# but then we may have to provide our own range for the clustering in `geo_kde_ranges` because default
# clustering ranges are for bond-lengths and 4-10 is not a true bond.
plot = biplot_kde(ds_stacked, 4, 10, pca_data=pca_res, scatter_color_property='geo', geo_kde_ranges=[(0,6), (8,15)])

../_images/tutorial_symlinks_0_full_tutorial_23_0.png

We can see, that the molecule isomerizes based on the distance between the H-atoms with indices 4 and 10. However, the PCA result provides us with more in-depth results. From the Biplot figures a-d we can already see that the central dihedrals around C-C bond (5,7) are most relevant for the isomerization classification.

If we want to gain deeper insights into the most relevant features, we can retrieve a description of the PCA result in textual form:

[13]:

# The parameters for `.explain_loadings(n1, n2)` are the number of most-weighted features per PCA component (n1) and overall across all components (n2)
print(pca_res.explain_loadings(5,5))

Maximum contributing features overall:
 dih(0,3,5,7) (weight: 0.4002896845340729) (Idxs: (0, 3, 5, 7))
 dih(3,5,7,9) (weight: 0.45366981625556946) (Idxs: (3, 5, 7, 9))
 dih(6,5,7,9) (weight: 0.4903942048549652) (Idxs: (6, 5, 7, 9))
 dih(3,5,7,8) (weight: 0.501832127571106) (Idxs: (3, 5, 7, 8))
 dih(6,5,7,8) (weight: 0.5060983896255493) (Idxs: (6, 5, 7, 8))


Maximum contributing features to component 0 :
 dih(8,7,9,10)  (weight: -0.16361398994922638) (Idxs: (8, 7, 9, 10))
 dih(3,5,7,9)  (weight: -0.41540634632110596) (Idxs: (3, 5, 7, 9))
 dih(3,5,7,8)  (weight: 0.44011834263801575) (Idxs: (3, 5, 7, 8))
 dih(6,5,7,8)  (weight: -0.44119080901145935) (Idxs: (6, 5, 7, 8))
 dih(6,5,7,9)  (weight: 0.4422222375869751) (Idxs: (6, 5, 7, 9))

Maximum contributing features to component 1 :
 dih(5,7,9,10)  (weight: -0.2583121359348297) (Idxs: (5, 7, 9, 10))
 dih(4,3,5,6)  (weight: 0.30467310547828674) (Idxs: (4, 3, 5, 6))
 dih(0,3,5,6)  (weight: -0.31193774938583374) (Idxs: (0, 3, 5, 6))
 dih(4,3,5,7)  (weight: -0.3276280462741852) (Idxs: (4, 3, 5, 7))
 dih(0,3,5,7)  (weight: 0.3666488528251648) (Idxs: (0, 3, 5, 7))

This confirms our suspicion that, indeed, the first principal component has its highest contributions from all dihedrals surrounding the bond 5,7. Similarly, the second component has its most significant contributions from dihedrals around the bond 3,5. This already tells us, that the isomerization around 5,7 is key to the structural and energetic landscape of the molecule. Similarly, the torsions around 3,5 provide degrees of freedom without a full isomerization of the molecule.

4) Getting a general overview with the Datasheet#

If there is no need for detailed analysis but instead, we want to get a quick overview of data within a dataset, shnitsel-tools offer a default visualization class, the Datasheet that provides a certain set of insights. For example, it can provide us with an overview of singlet-state statistics as well as triplet-state statistics, PCA-insights and general trajectory-level metadata.

As a last point, we will here generate a datasheet for our filtered Retinal dataset:

[20]:

from shnitsel.vis.datasheet import Datasheet

# Initialize the datasheet structures
sheet = Datasheet(truncated_dt,
                  structure_selection=structure_bd,
                  )

# To actually plot the datasheet, we can use the `plot()` command with which we can enable certain pages of the datasheet.
# With the `path` parameter, the output can be written to a pdf.
# Note, that this may take some time due to a wide range of analyses being run in parallel (can reach a few minutes depending on the dataset size)
sheet.plot(include_coupling_page=True, include_meta_page=True, path='./datasheet.pdf')

/home/tpadmin/git/shnitsel-tools-official/shnitsel/data/traj_combiner_methods.py:491: FutureWarning: In a future version of xarray the default value for compat will change from compat='equals' to compat='override'. This change will result in the following ValueError: Cannot specify both coords='different' and compat='override'. The recommendation is to set compat explicitly for this case.
  frames = xr.concat(
100%|██████████| 3/3 [00:00<00:00, 359.51it/s]
Written: 100%|██████████| 1/1 [00:01<00:00,  2.00s/page]

[20]:

{'/I02/agg/': [<Figure size 827x974.167 with 13 Axes>,
  <Figure size 150x150 with 9 Axes>,
  <Figure size 827x1169 with 7 Axes>]}

../_images/tutorial_symlinks_0_full_tutorial_40_2.png

../_images/tutorial_symlinks_0_full_tutorial_40_3.png

../_images/tutorial_symlinks_0_full_tutorial_40_4.png

Of course we also have the option to pick and choose which states and combinations to include in the analysis, if we so choose:

[21]:

from shnitsel.vis.datasheet import Datasheet

# Initialize the datasheet structures
sheet = Datasheet(truncated_dt,
                  structure_selection=structure_abd, # Use positions, bond lengths and dihedrals where necessary
                  state_selection="S, 1<>2, 2<>3" # Use all singlet states, transitions between states 1 and 2 as well as between states 2 and 3 but not between 1 and 3.
                  )

# We can also leave out certain pages of the datasheet in the plot outpu, like here, for example, we discard the PCA page and the meta information page bu
# due to the presence of an explicit selection, we will get a `selection` page before singlets.
sheet.plot(include_coupling_page=True, include_pca_page=False, include_meta_page=False, path='./datasheet_states_selected.pdf')

/home/tpadmin/git/shnitsel-tools-official/shnitsel/data/traj_combiner_methods.py:491: FutureWarning: In a future version of xarray the default value for compat will change from compat='equals' to compat='override'. This change will result in the following ValueError: Cannot specify both coords='different' and compat='override'. The recommendation is to set compat explicitly for this case.
  frames = xr.concat(
100%|██████████| 2/2 [00:00<00:00, 265.41it/s]
Written: 100%|██████████| 1/1 [00:02<00:00,  2.36s/page]

[21]:

{'/I02/agg/': [<Figure size 827x974.167 with 13 Axes>,
  <Figure size 827x974.167 with 13 Axes>,
  <Figure size 150x150 with 9 Axes>]}

../_images/tutorial_symlinks_0_full_tutorial_42_2.png

../_images/tutorial_symlinks_0_full_tutorial_42_3.png

../_images/tutorial_symlinks_0_full_tutorial_42_4.png

Full tutorial#

1) Import and preparation of data#

1.1) Loading data#

1.2) Cleaning data#