# Cyclized Type-2/3/4 Particle Cluster Rg Analysis

This script analyzes HOOMD GSD trajectories and measures the size and radius of gyration of **clusters of type-2/3/4 particles in cyclized molecules**.
Cyclization is decided at the block level, and cluster statistics are then extracted at the particle level over the last up to 20 frames.

## What this script does

For each experiment folder, the script performs the following steps:

1. Reads `init.txt` and reconstructs the particle-type layout for one molecule.
2. Builds contiguous folding-analysis blocks from particle types 2, 3, and 4.
3. Opens the last up to 20 frames of `sample.gsd`.
4. In each frame, classifies all molecules as `cyclized`, `misfolded`, or `unfolded` based on block-block contacts.
5. For cyclized molecules only, collects all type-2/3/4 particles and clusters them using a particle-particle minimum-image distance cutoff.
6. Computes cluster radius of gyration (Rg) and cluster size for each cluster.
7. Pools cluster statistics over the selected frames and writes a tab-delimited summary file.

The reported observables are, for type-2/3/4 particle clusters in cyclized molecules:

- mean cluster Rg ± SEM  
- mean cluster size (number of particles)  
- total number of clusters

## Cyclization criterion (block level)

Cyclization is determined **before** particle-level clustering and is based on contacts between contiguous folding blocks of particle types 2, 3, and 4.

For each molecule, the script:

- builds contiguous blocks of types 2/3/4 from `init.txt`,
- computes block centers of mass (COM) and equal-weight radii of gyration \(R_g\),
- evaluates minimum-image distances between all block pairs.

For each block pair \((i, j)\), it defines:

- \(\mathrm{COM}_i, \mathrm{COM}_j\): block COMs,
- \(R_{g,i}, R_{g,j}\): block radii of gyration,
- \(d_{ij}\): minimum-image distance between block centers,
- \(d_{\mathrm{cut},ij} = c (R_{g,i} + R_{g,j})\), where `c = RG_C`.

A block pair is considered in contact when

\[
d_{ij} < c (R_{g,i} + R_{g,j})
\]

Using these contacts, a molecule is labeled:

- `cyclized`  
  at least one (3,4) or (4,3) block pair is in contact, and no other block-block contact exists.
- `misfolded`  
  at least one contact exists, but not all contacts are 3–4 pairs.
- `unfolded`  
  no block-block contacts exist.

Only molecules labeled `cyclized` proceed to the particle-level clustering stage.

## Particle-level clustering and cluster Rg

After cyclized molecules have been identified in a frame, the script:

1. Collects all particles of type 2, 3, or 4 in each cyclized molecule.
2. Builds an undirected adjacency matrix where two particles are connected if their minimum-image distance is less than `CLUSTER_CUTOFF`.
3. Finds connected components of this graph; each connected component is treated as one cluster.
4. For each cluster:
   - computes cluster size (number of particles),
   - unwraps cluster coordinates relative to the first particle using the minimum-image convention,
   - computes cluster radius of gyration \(R_g\) from the unwrapped coordinates.

Cluster Rg is defined as the root-mean-square distance of cluster particles from the cluster center of mass in the unwrapped coordinate frame, using equal particle weights.

## Frames and averaging

- The script analyzes the **last up to 20 frames** of each trajectory:
  - If a trajectory has \(N\) frames, it uses frames \(\max(0, N-20), \dots, N-1\).
- In each analyzed frame, cyclized molecules are identified and clusters are built for those molecules.
- Cluster quantities from all frames are pooled before computing a single mean and SEM.

Important interpretation detail:

- `n_cyclized` in the output is a **frame-counted quantity**:
  - If the same molecule is cyclized in multiple frames, it contributes multiple times.
- `n_clusters_total` is the total number of clusters across all analyzed frames and all cyclized molecules.

## Input layout

The script expects a directory structure like:

```text
.
├── B2/
│   ├── run1/
│   │   ├── sample.gsd
│   │   └── init.txt
│   └── run2/
│       ├── sample.gsd
│       └── init.txt
├── B3/
│   └── ...
├── B4/
├── B5/
└── B6/
```

Each top-level folder listed in `TOP_B_FOLDERS` is scanned, and each subfolder is treated as one experiment if both `sample.gsd` and `init.txt` are present.

## Required files

### `sample.gsd`

- HOOMD trajectory file readable via `gsd.hoomd.open(...)`.
- The script uses random indexing into the trajectory to access frames.

### `init.txt`

- Whitespace-delimited file defining the one-molecule template.
- The first column is interpreted as integer particle type (`type_id`).
- The number of rows is interpreted as the number of particles per molecule.
- The script assumes:
  - all molecules in the trajectory share this layout,
  - particles are stored contiguously molecule-by-molecule in the GSD file.

## Assumptions

This script makes the following assumptions:

- The simulation box is **orthorhombic**, so a simple minimum-image wrapping formula is valid.
- Particles are stored contiguously by molecule in `sample.gsd`.
- All molecules have the same number of particles, equal to the row count of `init.txt`.
- Folding-analysis blocks are contiguous runs of particle types in `ALLOWED_TYPES = {2, 3, 4}`.
- Particle clusters are built only from particles of types 2, 3, and 4.
- Cluster Rg is computed after a local unwrap relative to the first particle in each cluster, which is appropriate for **compact clusters**.

## Output

For each top-level folder such as `B2`, the script writes:

```text
B2_summary_cyclized_particle_cluster_rg_234.txt
```

The output table has the following columns:

| Column | Meaning |
| --- | --- |
| `top_folder` | Top-level directory being processed. |
| `folder` | Experiment subfolder path. |
| `n_molecules` | Number of molecules inferred from the trajectory and `init.txt`. |
| `n_cyclized` | Frame-counted number of cyclized molecules over the analyzed frames. |
| `cluster_cutoff` | Particle-particle distance cutoff used for clustering. |
| `n_clusters_total` | Total number of clusters pooled across all frames and cyclized molecules. |
| `mean_cluster_rg` | Mean cluster radius of gyration over all clusters. |
| `sem_cluster_rg` | Standard error of the mean for the cluster Rg values. |
| `mean_cluster_size` | Mean cluster size (number of particles per cluster). |

## Dependencies

Install the Python packages used by the script:

```bash
pip install numpy pandas gsd
```

The trajectory reader relies on the `gsd.hoomd` Python API for HOOMD-native GSD files.

## Usage

Run the script from the directory that contains the `B2`, `B3`, `B4`, `B5`, and `B6` folders:

```bash
python your_cluster_script_name.py
```

If needed, edit the user-configurable section near the top of the script to change:

- `TOP_B_FOLDERS`
- `TARGET_GSD_FILENAME`
- `TARGET_INIT_FILENAME`
- `RG_C`
- `CLUSTER_CUTOFF`
- `ALLOWED_TYPES`

## Interpretation tips

- `RG_C` controls how strictly cyclization is defined at the **block level**. Increasing `RG_C` makes the contact criterion looser and generally increases the number of molecules labeled as cyclized.
- `CLUSTER_CUTOFF` controls cluster connectivity at the **particle level**:
  - Smaller values favor more, smaller clusters.
  - Larger values tend to merge clusters into bigger aggregates.
- Because `n_cyclized` and `n_clusters_total` are frame-counted, comparisons between simulations with different trajectory lengths or sampling frequencies should be made with care.

The script is intended for post-processing HOOMD trajectories to quantify the internal clustering of folding-related particles (types 2/3/4) inside molecules that satisfy a geometry-based cyclization rule.