QM9 MMFF94 Force-Field Failure Survey
Domain: physical_sciences — Chemistry & Materials Computation

Summary: Survey the QM9 quantum chemistry dataset to characterize systematic failures of the MMFF94 molecular mechanics force field relative to B3LYP/6-31G(2df,p) reference geometries, across a five-phase computational pipeline.

Task prompt:
You are surveying systematic failures of the MMFF94 molecular-mechanics force field relative to B3LYP/6-31G(2df,p) reference geometries across the QM9 dataset.

## Variant
`base`: Five-phase MMFF94 vs QM9 force-field failure survey on the full QM9 dataset.

## Input
- QM9 archive (bz2-compressed tar of 130,831 XYZ files): `base\input\dsgdb9nsd.xyz.tar.bz2`
- Stream the archive with `tarfile.open(..., "r:bz2")`. Do not extract it on disk.

## Software
- Python 3.12 is on PATH (`python`).
- RDKit, NumPy, pandas, and matplotlib are pre-installed in the system Python.
- No internet access is required or permitted.

## Output Directory
Write every result file directly inside:

```
base\output
```

## What You Must Produce

The following 10 files under the output directory:

1. `force_field_failures.csv` — Phase 1. For every QM9 molecule with exactly 2 heteroatoms (O/N/F): embed one ETKDG conformer with `randomSeed=42`, MMFF94-optimize, measure the heteroatom–heteroatom distance, and compare to the QM9 B3LYP distance. Emit rows with `|Discrepancy_A| >= 1.0` using columns `Molecule_ID, SMILES, QM9_Dist_A, RDKit_Dist_A, Discrepancy_A`. Flush continuously so a crash late in the scan does not lose earlier rows.
2. `phase2_classified.csv` — Phase 2. For every Phase 1 row, run a 200-conformer ETKDGv3 search (`randomSeed=42`), MMFF94-optimize each conformer, pick the lowest-energy conformer, and measure its heteroatom distance. Classify `genuine_ff_failure` when the residual discrepancy (|global-min dist − QM9 dist|) >= 0.5 Å, otherwise `sampling_artifact`. Columns: `Molecule_ID, SMILES, QM9_Dist_A, Naive_MMFF94_Dist_A, Naive_Discrepancy_A, Global_Min_Dist_A, Global_Min_Energy_kcal, Residual_Discrepancy_A, Classification`.
3. `phase3_scaffold_analysis.json` — Phase 3. Murcko-scaffold decomposition and SMARTS-based functional-group tagging of the genuine failures. Mean/max/min residual discrepancy per scaffold and per functional group, plus the top-5 worst molecules. Top-level keys: `survey_statistics`, `scaffold_analysis`, `top5_worst_molecules`, `all_genuine_failures`. **This JSON must not include QM9 molecule IDs** — identify molecules by SMILES and rank only.
4. `phase4_pes_results.json` — Phase 4. For each of the top-5 worst genuine failures ranked by residual discrepancy, run a constrained MMFF94 PES scan compressing the Phase 1 heteroatom pair from 8.0 Å down to 2.5 Å in 0.1 Å steps using `MMFFAddDistanceConstraint(relative=False, forceConstant=10000.0)`. Re-initialise the force field fresh at each step. Use the Phase 2 global-minimum geometry as the starting structure and RDKit indices after `AddHs`. Record the global minimum energy/distance, the quantum-forced energy at the exact QM9 distance, the conformational snap distance (arg-max |d²E/dd²|), the snap energy drop, and ΔE (quantum-forced − global-min). Top-level keys `rank_1`…`rank_5`, each with a `data` object containing `rank, molecule_id, canonical_isomeric_smiles, heteroatom_pair, qm9_heteroatom_distance, mmff94_global_min_distance, mmff94_global_min_energy, quantum_forced_energy, conformational_snap_distance, snap_energy_drop, delta_e`.
5. `pes_scan_rank1.png`, `pes_scan_rank2.png`, `pes_scan_rank3.png`, `pes_scan_rank4.png`, `pes_scan_rank5.png` — one PES plot per top-5 molecule.
6. `phase5_final_report.json` — Phase 5. Top-level keys `survey_statistics`, `scaffold_analysis`, `top5_characterization`, `top5_pes_summary`. The `top5_characterization` object must split the top-5 into `short_qm9_distance_molecules` (QM9 heteroatom distance < 3.5 Å) and `long_qm9_distance_molecules` (≥ 3.5 Å), plus a `failure_mode_distinction` text field explaining what physically distinguishes the two groups.

## Runtime Expectations
- Wall-clock budget: roughly 12–16 hours end to end. Phase 2 dominates (200-conformer search over ~635 candidates).
- Peak memory is observed in Phase 2 (~50–100 MB per molecule during conformer generation).
- Do not attempt to parallelise with OS-level process spawning inside the same Python process — sequential per-molecule processing is sufficient.


Required actions:
Phase 1: Scan all 130,831 QM9 molecules streamed from the bz2 archive; for those with exactly 2 heteroatoms (O/N/F), embed one ETKDG conformer (randomSeed=42), MMFF94-optimize, and save rows with `|Discrepancy_A| ≥ 1.0` to `force_field_failures.csv` with continuous flushing.
Phase 2: For every Phase 1 candidate run a 200-conformer ETKDGv3 search (randomSeed=42), MMFF94-optimize each conformer, pick the lowest-energy geometry, and classify each molecule as `genuine_ff_failure` (residual ≥ 0.5 Å) or `sampling_artifact` in `phase2_classified.csv`.
Phase 3: Perform Murcko scaffold decomposition and SMARTS-based functional-group analysis on the genuine failures, emitting `phase3_scaffold_analysis.json` with survey statistics, scaffold/group aggregates, and the top-5 worst molecules (SMILES + rank only, no QM9 IDs).
Phase 4: For the top-5 worst genuine failures, run a constrained MMFF94 PES scan from 8.0 Å to 2.5 Å in 0.1 Å steps (forceConstant 10000.0, relative=False, fresh force-field per step), save one PES plot per molecule (`pes_scan_rank1.png`…`pes_scan_rank5.png`), and record numerics in `phase4_pes_results.json`.
Phase 5: Consolidate everything into `phase5_final_report.json` with `survey_statistics`, `scaffold_analysis`, `top5_characterization` (split by <3.5 Å vs ≥3.5 Å QM9 heteroatom distance plus a `failure_mode_distinction` text field), and `top5_pes_summary`.

Software: Python
RDKit
NumPy
pandas
matplotlib

Input files:
name: dsgdb9nsd.xyz.tar.bz2
format: bzip2-compressed tar of XYZ files
description: QM9 dataset of 130,831 molecules with B3LYP/6-31G(2df,p) geometries; stream with tarfile, do not extract on disk.