Diels–Alder Causality Notebooks

These four notebooks form a workflow for storing molecular-dynamics data, extracting reaction descriptors, plotting feature relationships, and estimating causal effects for three Diels–Alder systems.

The notebooks are documented and cleaned versions of:

• store_h5md.ipynb — build, update, inspect, and diagnose the HDF5 data store.
• analysis.ipynb — extract reactive trajectory windows and create molecular feature tables.
• plot.ipynb — generate distributions, feature-pair plots, dependence matrices, and conditional plots.
• dowhy.ipynb — define the causal graph and estimate/refute causal effects.

Recommended workflow

Run the notebooks in this order:

1. store_h5md.ipynb — optional. Run only when creating, repairing, or adding data to da.h5. Do not run every cell automatically.
2. analysis.ipynb. Read reactive data from da.h5 and produce feature CSV files.
3. plot.ipynb. Create descriptive and dependence plots from the feature CSV files.
4. dowhy.ipynb. Run causal estimation and robustness checks from the feature CSV files.

The current da.h5 already contains the reactive data required by analysis.ipynb, so the first notebook can normally be skipped.

WSL2 project location

The data file is located at:

/home/felix/causality/my_dowhy/da.h5

The notebooks use relative paths such as da.h5, ./data_01, and ./linear. Their kernel must therefore use /home/felix/causality/my_dowhy as its working directory.

If the notebooks remain on the Windows Desktop, run this once near the beginning of a Jupyter session:

%cd /home/felix/causality/my_dowhy

Alternatively, copy working versions of the notebooks into the WSL project directory before launching Jupyter.

Start JupyterLab from WSL2 with:

source ~/miniconda3/bin/activate
cd /home/felix/causality/my_dowhy
jupyter lab

The environment checked on 2026-07-16 contains Python 3.10.19, JupyterLab 4.4.4, and h5py 3.15.1.

Expected project layout

/home/felix/causality/my_dowhy/
├── da.h5
├── analysis.ipynb
├── dowhy.ipynb
├── plot.ipynb
├── store_h5md.ipynb
├── data_01/
├── data_02/
├── data_04/
├── linear/                 # Used by plot.ipynb's table-combination section
└── reactive/               # Raw CFDT files used by store_h5md.ipynb

Create output directories before running cells that save files:

mkdir -p data_01 data_02 data_04 linear

Python dependencies

Core analysis and plotting

• h5py
• numpy
• pandas
• ase
• matplotlib
• seaborn
• scikit-learn
• openpyxl for Excel output

Causal inference

• dowhy
• econml
• lightgbm
• networkx
• pygraphviz
• Graphviz system libraries required by pygraphviz

The standard-library modules json, glob, re, and os are also used.

Notebook details

1. store_h5md.ipynb

Purpose: populate and inspect the H5MD-style da.h5 hierarchy.

Main inputs:

• Whole-system and component ASE .traj files.
• Per-frame electrophilicity JSON files.
• CFDT text outputs under reactive/.

Main operations:

• Store whole-system positions, charges, species, energies, and steps.
• Store diene and dienophile charges, energies, and steps.
• Store per-atom omega and total molecular total_omega values.
• Rename legacy datasets from data to value.
• Print the HDF5 hierarchy and plot one selected observable.

Warning: h5md_modify_dataset deletes an existing dataset and recreates it. The legacy-renaming and group-moving cells also change the file. Make a backup and run only the intended cell.

Important configuration:

sys
trajs
center_steps
ranges
ranges_omega

The component trajectories currently contain 2,400 frames. Their stored step ranges begin at center_step - 2000; the endpoint is determined by the trajectory length.

2. analysis.ipynb

Purpose: convert reactive HDF5 data into one feature table per selected system.

Select the system in the configuration cell:

da = "da-01-bu"  # or da-02-bu / da-04-bu

Extracted features:

1. Interaction energy per atom.
2. Residue-normalized electrophilicity, $\omega_k$ R1.
3. Charge transfer.
4. First forming-bond distance.
5. Second forming-bond distance.
6. Torsion angle.
7. Pyramidalization descriptor.
8. $\Delta\phi$ descriptor.

Primary outputs:

data_01/da-01-BU_features.csv
data_02/da-02-BU_features.csv
data_04/da-04-BU_features.csv

The notebook also writes trajectory .xyz files and produces:

• Mutual-information heatmaps.
• Pearson-correlation heatmaps.
• Time-series grids.
• Pairwise scatter grids.
• Marginal histogram grids.

The ASE view(molecules[0]) call opens an interactive molecular viewer. Comment it out when running headlessly.

3. plot.ipynb

Purpose: create publication-oriented summaries from the feature CSV files.

Available sections:

• Marginal histograms with total and per-trajectory populations.
• Pairwise scatter plots with 95% Gaussian covariance ellipses.
• Pairwise mutual-information matrices.
• Pearson-correlation matrices.
• Combined Excel/CSV causal-result tables.
• Gaussian conditional-distribution summaries and histograms.

Run the shared configuration cell first. Later plotting cells can then be run independently.

The table-combination section expects files named like:

linear/da-01-bu-2000-200_linear.csv

This differs from the current batch filename produced by dowhy.ipynb, which ends in _testforpaper.csv. Align the directory and filename templates before using the table-combination section.

4. dowhy.ipynb

Purpose: estimate the effect of each molecular descriptor on interaction energy using a domain-informed directed acyclic graph.

Supported datasets:

das = ["da-01-bu", "da-02-bu", "da-04-bu", "da-total-bu"]

The workflow:

1. Harmonizes feature names for graph tools.
2. Constructs the assumed causal graph.
3. Identifies an adjustment estimand with DoWhy.
4. Fits LinearDML with LightGBM nuisance models.
5. Provides an optional CausalForestDML implementation.
6. Runs placebo-treatment and random-common-cause refuters.
7. Displays or exports estimated and refuted effects.

The causal graph is a scientific assumption, not a graph learned from the data. Interpret causal estimates only after checking the graph, overlap, measurement assumptions, and robustness results.

The batch cell skips an analysis when its target output file already exists.

Systems and reactive trajectories

System	Primary analysis trajectories	Whole-system atoms
da_01_BU	2091, 3456	17
da_02_BU	1493, 2033	18
da_04_BU	0877, 1069	18

da.h5 also contains da_04_BU reactive trajectories 2314 and 3641. They are configured in store_h5md.ipynb but are not currently included in analysis.ipynb.

The analysis range is reaction-centered and inclusive:

-2000 through +200 = 2,201 frames per trajectory

Current da.h5 audit notes

The WSL2 data file was inspected read-only on 2026-07-16.

• HDF5 opens successfully and contains 370 groups and 430 datasets.
• All value/step lengths match, and step arrays are strictly increasing.
• All reactive windows required by analysis.ipynb are present.
• A full floating-point scan found no NaN or infinity values.
• The shorter da_01_BU/2091 omega series still covers the complete analysis window.
• There are 48 empty component placeholder groups for data not used by the four-notebook workflow.
• The layout is H5MD-style but lacks the /h5md metadata group required for complete H5MD compliance.

Known cleanup candidates:

• /observables/da_04_BU/non-reactive/0517/whole contains three exact duplicate groups whose names end in :.
• total_nucleophilicity: is the Hartree-scale version of the canonical eV-scale total_nucleophilic data; the fitted conversion factor is approximately 27.2113. Units should be recorded explicitly before removing or renaming either dataset.
• The file is about 933 MiB, but only about 257 MiB is active raw data. Approximately 709 MB is unaccounted HDF5 space, likely left by delete/recreate operations.
• None of the 430 datasets is compressed or chunked.

If space must be reclaimed, create a new file rather than repacking over the original:

h5repack da.h5 da.repacked.h5

Validate the new file fully before replacing or deleting the original.

Troubleshooting

FileNotFoundError: da.h5

Check the kernel's directory:

%pwd

Then switch to the WSL project:

%cd /home/felix/causality/my_dowhy

Output directory does not exist

Create it before saving:

mkdir -p data_01 data_02 data_04 linear

ModuleNotFoundError

Confirm that Jupyter is using the Miniconda environment at /home/felix/miniconda3 and install the missing package in that same environment.

pygraphviz or Graphviz layout errors

Both the Python pygraphviz package and Graphviz system libraries must be available. The manual NetworkX DAG plot can still run without pygraphviz_layout.

Unexpected HDF5 growth

Repeated delete/recreate operations can leave unused space inside an HDF5 file. Keep an original backup and periodically repack to a new file after validating all datasets.

Safety checklist

Before changing da.h5:

• Close other kernels or programs using the file.
• Make a timestamped backup.
• Confirm the selected system, trajectory, molecular part, and dataset path.
• Prefer creating a new dataset or repacked file over overwriting the only copy.
• Re-run the structural and numerical audit after modification.