Molecular Vibration¶

Molecular Vibration on Prism Hexamer example Image

Pipeline description¶

This example computes the Bond Occurrence Rate from an ensemble of electron density fields modeling a molecular system (the so-called Prism water hexamer) under molecular vibration. Specifically, this example partly reproduces the Figure 7 of the paper BondMatcher: H-Bond Stability Analysis in Molecular Systems.

The pipeline includes the following steps.

First, the appropriate datasets are read from disk using the filters CinemaReader, CinemaQuery, and CinemaProductReader. Specifically, the dataset GeomID=10 (top right) corresponds to the equilibrium state of the molecule, while the datasets GeomID=0 (bottom left) and GeomID=20 (top left) correspond to its most extreme vibrations.

Second, for convenience, the opposite of the density is considered and normalized with the ScalarFieldNormalizer to facilitate the subsequent processing steps.

Third, TopologicalSimplificationByPersistence is used to remove spurious extrema from the data, such that only 18 minima remains for each dataset (for the 18 atoms of the water hexamers).

Next, for each scalar field, the MorseSmaleComplex is computed. This complex is post-processed to cancel the least salient saddle-saddle pairs (see the checkbox Return Saddle Connectors). At this stage, the 1-dimensional separatrices of the complex are filtered, to only keep the saddle-minimum separatrices. This yields an Extremum Graph, which includes the minima of the scalar field and which connects them according to the saddle-minimum separatrices (representing chemical interactions, e.g., an H-bond or a covalent bond).

Finally, the Bond Occurrence Rate is computed from this extremum graph with the filter SeparatrixStability and it is visualized for the equilibrium state (bottom right). For each chemical interactionn, this quantity documents the rate of its occurrence in the ensemble (see the publication BondMatcher: H-Bond Stability Analysis in Molecular Systems for the formalization of this measure). Specifically, each occurrence rate is saved as the CellDataArray named Occurrence, itself associated to the output extremum graphs.

ParaView¶

To reproduce the above screenshot, go to your ttk-data directory and enter the following command:

paraview states/molecularVibration.pvsm

Python code¶

#!/usr/bin/env python
from paraview.simple import *

# create a new 'TTK CinemaReader'
tTKCinemaReader1 = TTKCinemaReader(DatabasePath="molecularVibration.cdb")

# create a new 'TTK CinemaQuery'
tTKCinemaQuery1 = TTKCinemaQuery(InputTable=tTKCinemaReader1)
tTKCinemaQuery1.SQLStatement = """SELECT * FROM InputTable0
WHERE ModeId=0 and GeomID in (10, 0, 20)"""

# create a new 'TTK CinemaProductReader'
tTKCinemaProductReader1 = TTKCinemaProductReader(Input=tTKCinemaQuery1)

# create a new 'Calculator'
calculator1 = Calculator(Input=tTKCinemaProductReader1)
calculator1.ResultArrayName = "opposite"
calculator1.Function = "-rho"

# create a new 'TTK ScalarFieldNormalizer'
tTKScalarFieldNormalizer1 = TTKScalarFieldNormalizer(Input=calculator1)
tTKScalarFieldNormalizer1.ScalarField = ["POINTS", "opposite"]

# create a new 'TTK TopologicalSimplificationByPersistence'
tTKTopologicalSimplificationByPersistence1 = TTKTopologicalSimplificationByPersistence(
    Input=tTKScalarFieldNormalizer1
)
tTKTopologicalSimplificationByPersistence1.InputArray = ["POINTS", "opposite"]
tTKTopologicalSimplificationByPersistence1.PersistenceThreshold = 5e-05

# create a new 'TTK MorseSmaleComplex'
tTKMorseSmaleComplex1 = TTKMorseSmaleComplex(
    Input=tTKTopologicalSimplificationByPersistence1
)
tTKMorseSmaleComplex1.ScalarField = ["POINTS", "opposite"]
tTKMorseSmaleComplex1.DiscreteGradientBackend = "Stochastic algorithm (IEEE TVCG 2012)"
tTKMorseSmaleComplex1.ReturnSaddleConnectors = 1
tTKMorseSmaleComplex1.SaddleConnectorsPersistenceThreshold = 2.5e-05

# create a new 'Threshold'
threshold1 = Threshold(Input=OutputPort(tTKMorseSmaleComplex1, 1))
threshold1.Scalars = ["CELLS", "SeparatrixType"]
threshold1.LowerThreshold = 0
threshold1.UpperThreshold = 0

# create a new 'TTK SeparatrixStability'
tTKSeparatrixStability1 = TTKSeparatrixStability(Input=threshold1)

# Extract the entry #1, corresponding to the equilibrium state
extractBlock11 = ExtractBlock(Input=tTKSeparatrixStability1)
extractBlock11.Assembly = "Hierarchy"
extractBlock11.Selectors = ["/Root/Block1"]

SaveData("bondOccurrenceRate.vtm", extractBlock11)