BAMresearch/ZrV2O7-PDF-Refinement: 1.2

# PDF Structural Refinement Framework (v22)

**Author:** Tomasz Stawski  
**Contact:** tomasz.stawski@gmail.com | tomasz.stawski@bam.de  
**Version:** 1.2

---

## Overview

This repository contains a modular, object-oriented Python framework for automated, multi-stage structural refinement of crystalline materials using the **Pair Distribution Function (PDF)** technique. PDF analysis extracts real-space structural information from total X-ray scattering data, enabling the study of both long-range and local atomic structure.

The framework was developed in the context of temperature-dependent PDF studies of zirconium vanadate (ZrV₂O₇) and is designed to handle large, sequential datasets , e.g., measured during a single synchrotron experiment, in a fully automated manner. It supports progressive symmetry reduction, rigid-body geometric restraints, multiple optimisation strategies, structural simulation, and post-refinement forensic analysis.

The codebase is split into three files:

| File | Role |
|---|---|
| `sample_refinement_v22_classes.py` | Core library: all class definitions |
| `sample_refinement_v22_execution.py` | Main execution script: refinement workflow |
| `sample_refinement_v22_simulation_and_analysis.py` | Simulation and structural forensics workflow |

---

## Requirements

### Python Version

Due to `diffpy` dependencies, either Python 3.7 or 3.12 is required depending on the version.

### Dependencies

All dependencies are available via `pip` or `conda`. The following packages must be installed:

> **Note:**
> This program uses `pdfgetX3` (via the `diffpy.pdfgetx` package), which requires a special license. If `pdfgetX3` is not installed, the software can still be used, but no on-the-fly conversion of raw diffraction data to G(r) will be possible. In this case, the user must supply their own pre-calculated G(r) files matching the requested r-range.

| Package | Purpose |
|---|---|
| `diffpy.pdfgetx` | PDF generation from raw diffraction data (optional if providing G(r) files) |
| `diffpy.srfit` | Structural refinement framework (`FitRecipe`, `FitContribution`, `PDFGenerator`, etc.) |
| `diffpy.structure` | Crystal structure handling; CIF loading |
| `diffpy.srreal` | Real-space PDF calculator |
| `pyobjcryst` | Alternative crystal object loading |
| `numpy` | Numerical arrays and operations |
| `scipy` | Optimisation algorithms (`minimize`, `least_squares`, `basinhopping`) |
| `pandas` | DataFrame-based results export |
| `matplotlib` | Plotting |
| `seaborn` | Statistical histograms and KDE plots |
| `tqdm` | Progress bars |
| `dill` | Advanced object serialisation (for state checkpointing) |
| `psutil` | CPU utilisation monitoring for parallel setup |
| `multiprocessing` | Standard library; parallel PDF calculation |
| `logging` | Standard library; progress logging |

Installing via conda (recommended for `diffpy` packages):

It is highly recommended to install `diffpy` packages into a dedicated conda environment using the `conda-forge` channel:

```bash
# Add conda-forge to your channels
conda config --add channels conda-forge

# Create a new environment and install the main packages
conda create -n pdf_refinement_env diffpy.cmi diffpy.pdfgetx diffpy.srfit diffpy.structure pyobjcryst seaborn tqdm dill psutil
conda activate pdf_refinement_env
```

---

## Directory Structure

The framework expects the following directory layout relative to the working directory:

```
.
├── data/                          # Raw X-ray scattering data files (.dat, .xy, QA format)
├── CIFs/                          # Input CIF files for structural models
├── fits/                          # Output directory for all refinement results
│   └── <project_name>/
│       └── <timestamp>/           # Per-run timestamped output directory
│           ├── fitting_summary_<i>.txt
│           ├── fitting_<i>.png
│           ├── fitting_curve_<i>.csv
│           ├── Phase0_<i>.cif
│           ├── fit_state_<i>.txt
│           ├── rw_log_<i>.csv
│           └── <i>_Phase0_summary_plot.pdf
├── resultsSimulations/            # Output directory for simulation runs
├── sample_refinement_v22_classes.py
├── sample_refinement_v22_execution.py
└── sample_refinement_v22_simulation_and_analysis.py
```

---

## Architecture

The library (`sample_refinement_v22_classes.py`) defines seven classes, each with a distinct responsibility:

```
RefinementConfig
    └── Validates configuration; creates output directories

StructureAnalyzer
    └── Geometric calculations: bonds, angles, dihedrals
    └── Generates algebraic constraint expressions

ResultsManager
    └── Logging, plotting, CIF export, results saving

PDFManager
    └── PDF generation from raw data (via PDFGetter)
    └── PDFContribution construction from CIF files

RefinementHelper
    └── Static utilities: space group parameter mapping

PDFRefinement
    └── Core class: builds and executes a single FitRecipe
    └── Manages symmetry changes, rigid-body restraints

PDFWorkflowManager  (inherits from PDFRefinement)
    └── High-level orchestration: sequential multi-dataset workflow
    └── Simulation workflow
    └── Structural forensics
```

The two execution scripts instantiate these classes, configure them via a dictionary, and call high-level workflow methods.

---

## Configuration

All refinement parameters are defined in a single Python dictionary (`project_config`) at the top of the execution script. No separate configuration files are required. The complete set of required keys is validated at runtime by `RefinementConfig`; missing keys raise a `KeyError` with a descriptive message.

### Required Configuration Keys

| Key | Type | Description |
|---|---|---|
| `project_name` | `str` | Project label; used as the output subdirectory name under `fit_directory` |
| `xrd_directory` | `str` | Path to the directory containing raw scattering data |
| `cif_directory` | `str` | Path to the directory containing input CIF files |
| `fit_directory` | `str` | Root output directory for refinement results |
| `dataset_list` | `list[str]` | Ordered list of data filenames to refine sequentially |
| `ciffile` | `dict` | Maps each CIF filename to `[space_group, periodic, supercell]` |
| `composition` | `str` | Chemical formula for PDF generation (e.g., `'O7 V2 Zr1'`) |
| `detailed_composition` | `dict` | Per-element metadata (see below) |
| `qdamp` | `float` | Instrumental resolution damping parameter (Å⁻¹) |
| `qbroad` | `float` | Instrumental peak broadening parameter (Å⁻¹) |
| `qmax` | `float` | Maximum Q value used in the Fourier transform (Å⁻¹) |
| `adp_mode` | `str` | Atomic displacement parameter mode: `'isotropic'`, `'anisotropic'`, or `'fixed_shape'` |
| `sgoffset` | `list[float]` | Space group origin offset `[dx, dy, dz]` for `constrainAsSpaceGroup` |
| `myrange` | `tuple[float, float]` | Full r-range for PDF calculation `(rmin, rmax)` in Å |
| `myrstep` | `float` | Step size Δr for the PDF grid (Å) |
| `convergence_options` | `dict` | Options passed to the SciPy optimiser (e.g., `{'disp': True}`) |
| `pdfgetx_config` | `dict` | PDFGetter configuration (see below) |
| `log_file` | `str` | Path to the progress log file |
| `refinement_plan` | `dict` | Stage-by-stage refinement strategy (see below) |
| `refine_qdamp` | `bool` | Whether to include `qdamp` as a free variable |
| `refine_qbroad` | `bool` | Whether to include `qbroad` as a free variable |
| `start_each_dataset_fresh` | `bool` | If `True`, each dataset begins from the original CIF; if `False`, from the previous dataset's refined structure |
| `optimizer` | `str` | Optimiser choice: `'minimize'`, `'least_squares'`, or `'basinhopping'` |
| `optimizer_method` | `str` | Method string passed to the chosen optimiser |

### `detailed_composition` Dictionary

This dictionary provides per-element structural metadata for bond detection and ADP initialisation:

```python
'detailed_composition': {
    'Zr': {
        'symbol': 'Zr',
        'Uiso': 0.0065,           # Initial isotropic ADP (Ų)
        'polyhedron_center': True, # Acts as a polyhedral centre (bond source)
        'polyhedron_vertex': False,
        'cutoff': (1.8, 2.2)      # Bond length range (Å) for centre–vertex detection
    },
    'V': {
        'symbol': 'V',
        'Uiso': 0.0100,
        'polyhedron_center': True,
        'polyhedron_vertex': False,
        'cutoff': (1.5, 2.4)
    },
    'O': {
        'symbol': 'O',
        'Uiso': 0.0250,
        'polyhedron_center': False,
        'polyhedron_vertex': True  # Acts as a polyhedral vertex (bond target)
    },
}
```

Elements designated as `polyhedron_center` are paired with elements designated as `polyhedron_vertex` to identify bonds. The `cutoff` tuple `(min, max)` specifies the Cartesian distance range (Å) used to determine valid coordination bonds.

### `pdfgetx_config` Dictionary

Passed directly to `diffpy.pdfgetx.PDFConfig` during PDF generation:

```python
'pdfgetx_config': {
    'mode': 'xray',     # Scattering type: 'xray' or 'neutron'
    'dataformat': 'QA', # Input data format (Q in Å⁻¹)
    'rpoly': 1.3,       # Low-r polynomial correction cutoff (Å)
    'qmin': 0.0         # Minimum Q value (Å⁻¹)
}
```

### `refinement_plan` Dictionary

Defines the multi-stage refinement strategy. Each integer key corresponds to one refinement stage:

```python
refinement_plan = {
    0: {
        'description': 'Initial fit with Pa-3 symmetry',
        'space_group': ['Pa-3'],          # Space group(s), one per phase
        'enforce_other_lattice': True,     # Enforce high-symmetry lattice during symmetry reduction
        'constraints': {
            'constrain_bonds':     (True, 0.001),  # (enable, sigma)
            'constrain_angles':    (True, 0.001),
            'constrain_dihedrals': (False, 0.001),
            'adaptive': True               # Recalculate bond vectors at each stage
        },
        'fitting_range': [1.5, 27],        # r-range for fitting (Å)
        'fitting_order': ['lat', 'scale', 'psize', 'delta2', 'adp', 'xyz', 'all']
    },
    1: {
        'description': 'Symmetry reduction to P23',
        'space_group': ['P23'],
        ...
    },
    2: {
        'description': 'Lowest symmetry (P1)',
        'space_group': ['P1'],
        ...
    }
}
```

The `fitting_order` list controls the sequential freeing of parameter groups. Valid tags include:

| Tag | Parameters freed |
|---|---|
| `'lat'` | Lattice parameters |
| `'scale'` | Scale factors |
| `'psize'` | Nanoparticle size parameter |
| `'delta2'` | Correlated motion broadening |
| `'adp'` | Atomic displacement parameters |
| `'xyz'` | Atomic coordinates |
| `'all'` | All free parameters simultaneously |

---

## Usage

### 1. Running a Sequential Refinement

The primary use case is automated sequential refinement of multiple datasets:

```python
# sample_refinement_v22_execution.py (simplified)

project_config = { ... }   # Define all required parameters
refinement_plan = { ... }  # Define the multi-stage strategy
project_config['refinement_plan'] = refinement_plan

import multiprocessing, numpy as np, psutil
from multiprocessing import Pool
from sample_refinement_v22_classes import (
    RefinementConfig, StructureAnalyzer, ResultsManager,
    PDFManager, RefinementHelper, PDFWorkflowManager
)

if __name__ == '__main__':
    config = RefinementConfig(project_config)

    syst_cores = multiprocessing.cpu_count()
    cpu_percent = psutil.cpu_percent()
    avail_cores = np.floor((100 - cpu_percent) / (100.0 / syst_cores))
    ncpu = int(np.max([1, avail_cores]))
    pool = Pool(processes=ncpu)

    analyzer        = StructureAnalyzer(config.detailed_composition)
    results_manager = ResultsManager(config, analyzer)
    pdf_manager     = PDFManager(config, ncpu, pool)
    helper          = RefinementHelper()

    workflow = PDFWorkflowManager(
        config, pdf_manager, results_manager, helper, analyzer, ncpu, pool
    )
    workflow.run_sequential_workflow()
```

The `run_sequential_workflow()` method iterates through all datasets in `dataset_list`, applies the full `refinement_plan` to each, and saves all results.

### 2. Choosing an Optimiser

Three optimisation strategies are available, selected via the `optimizer` configuration key:

| `optimizer` | `optimizer_method` options | Recommended use |
|---|---|---|
| `'minimize'` | `'L-BFGS-B'`, `'SLSQP'`, `'trust-constr'` | General-purpose; fast local optimisation |
| `'least_squares'` | `'trf'`, `'dogbox'`, `'lm'` | Dedicated least-squares; well-suited to residual minimisation |
| `'basinhopping'` | `'L-BFGS-B'` | Global optimisation; use when the fit is likely trapped in a local minimum |

For `'basinhopping'`, additional options are passed via:

```python
'basinhopping_options': {'stepsize': 50, 'niter': 100, 'T': 200.0}
```

### 3. Atomic Displacement Parameter Modes

The `adp_mode` key controls how ADPs are handled:

| Mode | Description |
|---|---|
| `'isotropic'` | All atoms described by a single isotropic Uiso per element; anisotropy tensor is ignored |
| `'anisotropic'` | Full U-tensor refinement (U11, U22, U33, U12, U13, U23); most parameters |
| `'fixed_shape'` | Anisotropic shape preserved from the CIF; Uiso serves as a scale factor for the tensor |

### 4. Applying Rigid-Body Geometric Restraints

Geometric restraints are soft penalties added to the cost function to maintain chemically realistic geometry. They are specified per stage in `refinement_plan` and applied automatically during the workflow. They can also be applied manually:

```python
workflow.apply_rigid_body_constraints(
    constrain_bonds=(True, 0.001),      # (enable, sigma)
    constrain_angles=(True, 0.001),
    constrain_dihedrals=(False, 0.001),
    adaptive=True  # Recalculate bond vectors from current structure
)
```

- **Bond restraints** enforce that centre–vertex bond lengths remain within ±5% of their initial values. Bonds near cell edges (`EDGE`), shared-vertex bonds, and bonds outside the defined cutoff range (`PROBLEMATIC`) are automatically assigned a tighter sigma of 10⁻⁷.
- **Angle restraints** enforce bond angles within ±1° of their initial values (in radians).
- **Dihedral restraints** enforce dihedral angles within ±2° of their initial values.
- The `adaptive` flag, when `True`, recalculates all bond vectors from the current refined structure before applying restraints; when `False`, the bond vectors calculated at the start of the workflow are reused.

### 5. Controlling Sequential Dataset Initialisation

Two modes exist for how each dataset initialises its structure:

- **`start_each_dataset_fresh: True`**: Every dataset begins from the original CIF (or from `special_structure` if defined).
- **`start_each_dataset_fresh: False`**: Each dataset, after the first, is initialised from the refined structure of the previous dataset. This is appropriate for temperature-series data where parameters evolve smoothly.

An optional `special_structure` entry allows a specific pre-refined CIF to be used as the starting point for the first dataset:

```python
'special_structure': {
    'file_path': 'fits/previous_run/Phase0_2.cif',
    'phase_index_to_update': 0
}
```

### 6. Multi-Phase Modelling

Multiple structural phases are specified by adding entries to the `ciffile` dictionary:

```python
'ciffile': {
    'phase_A.cif': ['Pa-3', True, (1, 1, 1)],  # CIF, periodic, supercell
    'phase_B.cif': ['Pa-3', True, (1, 1, 1)]
}
```

The PDF model equation is automatically constructed as:

```
G(r) = s_Phase0 · Phase0(r) · sphere_Phase0(r) + s_Phase1 · Phase1(r) · sphere_Phase1(r) + ...
```

Each phase has independent scale (`s_PhaseN`), nanoparticle size (`psize_PhaseN`), and correlated motion broadening (`delta2_PhaseN`) parameters. At finalisation, the framework saves partial PDFs for each individual phase.

### 7. Supercell Expansion

A supercell can be specified as a three-element tuple in the `ciffile` entry:

```python
'ciffile': {'model.cif': ['Pa-3', True, (3, 3, 3)]}
```

This triples the unit cell in all three directions before refinement. Both the standard (`PDFGenerator`) and Debye (`DebyePDFGenerator`) approaches are available.

---

## Simulation Workflow

For forward simulation — i.e., calculating a theoretical PDF from a known structure without refinement — use `sample_refinement_v22_simulation_and_analysis.py`. This is useful for model validation and for generating baseline patterns.

A separate `simulation_data` dictionary configures the simulation:

```python
simulation_data = {
    'cif_directory': '/path/to/refined/structure/',
    'ciffile': {'refined_model.cif': ['P1', True, (1, 1, 1)]},
    'powder_data_file': 'PDF_sample_25C.dat',  # Provides the r-grid
    'output_path': 'resultsSimulations/25C_run',
    'optimized_params': {
        'Phase0': {
            's': 4.92e-01,        # Scale factor
            'psize': 266.7,       # Nanoparticle size (Å)
            'delta2': 2.54        # Correlated motion parameter
        }
    },
    'default_Uiso': {
        'Zr': 5.79e-04,           # Isotropic ADP (Ų)
        'V':  3.22e-03,
        'O':  7.22e-03
    },
    'fitting_range': [1.5, 27],   # r-range for output (Å)
    'csv_filename': 'sim_vs_obs.csv'
}
```

Execute:

```python
workflow.simulate_pdf_workflow(main_config=config, sim_config=simulation_data)
```

This produces:
- A CSV file with observed, simulated, and difference G(r).
- A summary plot (PDF comparison, bond-length histograms, bond-angle histograms).

---

## Structural Forensics

Post-refinement structural forensics analyses bond lengths in any CIF file (e.g., a checkpoint from a refinement), identifies outliers, and categorises bonds by connectivity type (terminal, bridging, linking).

### Text Report

```python
workflow.run_structural_forensics(
    cif_path='fits/run_01/Phase0_2.cif',
    target_bond='V-O',           # Bond type to analyse
    outlier_threshold=1.8        # Upper threshold for flagging outliers (Å)
)
```

Output includes:
- Global statistics (count, mean, standard deviation, range).
- Detailed tables of threshold outliers and statistical deviations (outside mean ± 1σ).
- Per-connectivity-type outlier percentages.

### Graphical Report

```python
workflow.visualize_structural_forensics(
    cif_path='fits/run_01/Phase0_2.cif',
    target_bond='V-O',
    outlier_threshold=1.8,
    output_dir='forensics_results/Dataset_105C'
)
```

Saves plots to `output_dir`.

---

## Output Files

For each dataset and each refinement stage `i`, the following files are generated in the timestamped output directory:

| File | Description |
|---|---|
| `fitting_summary_<i>.txt` | Full DiffPy `FitResults` text report; includes refined values, uncertainties, and Rw |
| `fitting_<i>.png` | Plot of G_obs, G_calc, and G_diff |
| `fitting_curve_<i>.csv` | G_obs, G_calc, G_diff as a function of r |
| `Phase0_<i>.cif` | Refined crystal structure exported as a CIF |
| `fit_state_<i>.txt` | Snapshot of all variable names and values from `fit.show()` |
| `rw_log_<i>.csv` | Convergence log: Rw value at each optimiser iteration |
| `<i>_Phase0_summary_plot.pdf` | Three-panel summary: PDF fit, bond-length histogram, bond-angle histogram |
| `_final_extrapolated_fit.png` | Fit extrapolated over the full r-range (`myrange`) |
| `_final_extrapolated_fit.csv` | Data for the extrapolated fit |
| `<phase>_partial_fit.png` | Partial PDF contribution for each phase (multi-phase models only) |
| `<phase>_partial.csv` | Data for the partial PDF contribution |

---

## Key Classes: Reference

### `RefinementConfig`

Validates a configuration dictionary and exposes all parameters as instance attributes. The method `new_output_directory(subdir_name=None)` creates a new timestamped output directory during a run.

### `StructureAnalyzer`

Performs all crystallographic geometry calculations:

| Method | Description |
|---|---|
| `get_polyhedral_bond_vectors(phase)` | Returns all centre–vertex and vertex–vertex bond vectors and lengths |
| `find_bond_pairs(bond_vectors, added_params)` | Filters bonds to those involving refined atoms only |
| `find_angle_triplets(bond_pairs)` | Identifies unique X–Y–Z angle triplets |
| `find_dihedral_quadruplets(bond_pairs)` | Identifies unique A–B–C–D dihedral quadruplets |
| `create_constrain_expressions_for_bonds(...)` | Generates algebraic string expressions for bond-length restraints |
| `create_constrain_expressions_for_angles(...)` | Generates algebraic string expressions for angle restraints |
| `create_constrain_expressions_for_dihedrals(...)` | Generates algebraic string expressions for dihedral restraints |
| `detect_edge_bonds(bond_vectors)` | Identifies bonds near unit cell boundaries |

### `PDFManager`

| Method | Description |
|---|---|
| `generatePDF(...)` | Computes G(r) from a raw data file using `diffpy.pdfgetx.PDFGetter` |
| `contribution(name, qmin, qmax, ciffile, periodic, super_cell)` | Creates a `PDFGenerator` from a CIF |
| `DebyeContribution(...)` | Creates a `DebyePDFGenerator` from a CIF |
| `build_contribution(r0, g0, cfg, cif_info, fitRange)` | Builds a multi-phase `PDFContribution` |

### `PDFRefinement`

| Method | Description |
|---|---|
| `build_initial_recipe()` | Constructs the initial `FitRecipe` with all parameters |
| `run_refinement_step(i, fitting_range, ...)` | Executes one full refinement stage |
| `modify_recipe_spacegroup(spacegroup_list, enforce_other_lattice)` | Switches the active space group and regenerates the parameter set |
| `apply_rigid_body_constraints(...)` | Adds geometric restraints to the recipe |
| `update_recipe_from_initial(fit_old, fit_new, cpdf_new)` | Transfers refined values from one recipe to the next |
| `compare_fits(fitA, fitB)` | Compares two `FitRecipe` objects and reports differences |
| `rebuild_contribution(old_cpdf, r_obs, g_obs)` | Rebuilds a `PDFContribution` with new data, preserving structure |

### `PDFWorkflowManager` (inherits `PDFRefinement`)

| Method | Description |
|---|---|
| `run_sequential_workflow()` | Main entry point; iterates over all datasets and all refinement stages |
| `simulate_pdf_workflow(main_config, sim_config)` | Standalone simulation from a fixed structure |
| `run_structural_forensics(cif_path, target_bond, outlier_threshold)` | Text-based bond analysis |
| `visualize_structural_forensics(cif_path, target_bond, ...)` | Graphical bond analysis |

---

## Parallelisation

Parallel PDF calculation is enabled automatically via Python's `multiprocessing.Pool`. The number of processes is determined at runtime from the number of available CPU cores and current CPU utilisation:

```python
syst_cores = multiprocessing.cpu_count()
cpu_percent = psutil.cpu_percent()
avail_cores = np.floor((100 - cpu_percent) / (100.0 / syst_cores))
ncpu = int(np.max([1, avail_cores]))
pool = Pool(processes=ncpu)
```

The pool is passed to `PDFManager` and activated per phase generator via `PDFGenerator.parallel(ncpu=ncpu, mapfunc=pool.map)`.

**Important:** All code that creates a `Pool` and invokes the workflow must be enclosed in a `if __name__ == '__main__':` guard to prevent recursive subprocess spawning on systems using the `'spawn'` start method (Windows, macOS).

---

## Goodness-of-Fit Metric

The weighted R-factor (R_w) is used as the primary goodness-of-fit metric:

$$R_w = \sqrt{\frac{\sum_i (G_{\mathrm{obs}}(r_i) - G_{\mathrm{calc}}(r_i))^2}{\sum_i G_{\mathrm{obs}}(r_i)^2}}$$

R_w is logged at every optimiser iteration to `rw_log_<i>.csv` and is reported in the title of the `summary_plot.pdf` for each stage.

---

## Licence

This software is distributed under the terms specified in `LICENSE`. When used in or to support a scientific publication, appropriate citation is requested.

---

## Notes on Reproducibility

- All refined structures are exported as CIF files at every stage, enabling reproducibility and re-initialisation from intermediate points.
- The `special_structure` entry in the configuration allows any prior CIF to be injected as the starting model.
- The `start_each_dataset_fresh: False` mode propagates refined structures through a temperature series, which is appropriate for studies of continuous structural evolution.
- Log files (`refinement_log.txt` and `rw_log_<i>.csv`) provide a complete record of the refinement history.