# Polygon-based Patchy Particle Assembly (ALJ Model)

This script runs a HOOMD-blue simulation of **rigid polygon particles with patchy interactions** built from external geometry files.  
Each rigid body consists of a polygonal core and two patch types, and the system is compressed to a target packing fraction and then driven through a thermal cycle to explore assembly behavior.

## Overview

The script:

- reads a polygon vertex file (`polygon_verts.txt`) and a rigid-body patch file (`polygon_rigid.txt`),
- builds a rigid-body model with:
  - a polygon core (type `A`),
  - two patch species (types `B` and `C`),
- defines an **anisotropic ALJ potential** for polygon cores and core–patch excluded volume,
- uses isotropic LJ (`gauss`) interactions for patch–patch attractions,
- uses Yukawa interactions for patch–patch repulsions,
- compresses the system to a **target packing fraction** in 2D,
- turns on patch attractions/repulsions,
- applies a temperature cycle to probe different assembly states,
- writes trajectories and shape information to `sample.gsd`, along with a simple status table.

## Inputs

The script expects two text files in the working directory:

### `polygon_verts.txt`

Each line corresponds to one polygon vertex and contains:

```text
x  y  z  sigma_in  sigma_out  v_poly  Lbox
```

Typical interpretation:

- `x, y, z` – vertex position relative to the polygon center,
- `sigma_in` – inner length scale (e.g. in-radius),
- `sigma_out` – outer length scale (e.g. circumradius),
- `v_poly` – scalar derived from the polygon geometry (e.g. polygon area),
- `Lbox` – reference box length, carried along for convenience.

These vertices are used to define the anisotropic shape of the core via HOOMD’s `md.pair.aniso.ALJ` interaction.

### `polygon_rigid.txt`

Each line has the format:

```text
type  x  y  z  patch_size
```

where:

- `type` – integer index mapping into the type list `['A','B','C']`:
  - `0` → `A` (core),
  - `1` → `B` (patch type 1),
  - `2` → `C` (patch type 2),
- `x, y, z` – patch position relative to the core center,
- `patch_size` – size parameter for the patch (used as a length scale).

These entries are used to attach patch particles to the rigid body.

## Particle types and roles

The script defines three particle types:

- `A` – polygon core:
  - anisotropic ALJ shape using the vertex list from `polygon_verts.txt`,
  - carries the rigid body and rotational degrees of freedom.

- `B` – patch type 1:
  - spherical “cap” modeled as an ALJ shape with a finite rounding radius,
  - receives **attractive LJ interactions** with itself (B–B) after the compression stage.

- `C` – patch type 2:
  - also modeled as a spherical cap,
  - receives **repulsive Yukawa interactions** with itself (C–C) after the compression stage.

All patches belonging to a given rigid body share the same body ID as the core and move rigidly with it.

## Interactions

The model uses three main interaction terms:

1. **Anisotropic ALJ (ALJ)**  
   - `A–A`: core–core anisotropic interactions based on the polygon shape.
   - `A–B`, `A–C`: core–patch excluded volume via spherical caps.
   - Controlled by:
     - `epsilon_alj` – strength,
     - `sigma_core`, `sigma_patch` – size parameters.

2. **Isotropic LJ (`gauss`)**  
   - Initially zero for all pairs; later updated to:
   - `B–B`: attractive LJ with `epsilon_att` and `sigma = sigma_patch`.
   - Used as a simple isotropic attraction between patches to promote clustering or face-face contacts.

3. **Yukawa**  
   - Initially zero for all pairs; later updated to:
   - `C–C`: repulsive Yukawa with `epsilon_yukawa` and `kappa_yukawa`.
   - Used to introduce longer-range, soft repulsion for one patch species.

Neighbor lists are built with `hoomd.md.nlist.Cell` excluding interactions within the same rigid body.

## Simulation protocol

The simulation proceeds in several stages:

### 1. System initialization

- Load one polygon + patch rigid body at the origin.
- Replicate it in a dense 2D grid:
  - `nrep` times along x,
  - `nrep` times along y.
- The initial box is taken from `Lbox` and then replicated using `sim.state.replicate`.

The total polygon “area” is estimated as:

\[
V_{\text{total}} = n_{\text{rep}}^2 \, v_{\text{poly}}
\]

and is used to define a target packing fraction.

### 2. Equilibration at initial density

- Short equilibration run at temperature `kBT = 1.0`.
- Integrator: `hoomd.md.methods.ConstantVolume` with MTTK thermostat.
- Only cores (`A`) are thermostatted directly; patches follow rigid body motion.

### 3. Compression

- The simulation box is compressed in 2D using `hoomd.update.BoxResize` with a `Ramp` variant.
- The final volume is set such that:

\[
\phi_{\text{target}} = \frac{V_{\text{total}}}{V_{\text{box}}}
\]

where `phi_target` is a user-defined packing fraction (e.g. 0.4).

### 4. Turn on patch interactions

After compression:

- `B–B` LJ attraction is enabled:

  - `epsilon = epsilon_att`,
  - `sigma = sigma_patch`,
  - cutoff set to `3 * rcut_patch_att_lj`.

- `C–C` Yukawa repulsion is enabled:

  - `epsilon = epsilon_yukawa`,
  - `kappa = kappa_yukawa`,
  - cutoff set to `3 * sigma_patch`.

A short production run is performed under these conditions.

### 5. Thermal cycle

- The NVT method is reset and a **temperature cycle** is defined using `hoomd.variant.Cycle`:

  - lower temperature `kBT_low = kBT`,
  - higher temperature `kBT_high = 1.5`,
  - hold times and ramp times set by `time_scale`.

- A long run is performed over many cycles to explore different assembly states and potential ordering transitions.

## Output

The script writes:

- `sample.gsd`  
  A GSD-format trajectory containing particle positions, types, rigid body information, and ALJ type shapes.  
  The writer attaches a `Logger` that includes `alj.type_shapes`, so shape metadata is stored in the GSD.

It also prints to the terminal:

- current timestep, TPS (timesteps per second),
- estimated time remaining (ETR),

using a `hoomd.write.Table` writer and a small `Status` helper class.

## Dependencies

The script requires:

- `hoomd` (HOOMD-blue, v3+),
- `gsd`,
- `numpy`,
- standard Python 3 modules (`datetime`, `random`).

Install, for example, with:

```bash
pip install hoomd gsd numpy
```

(Ensure you follow HOOMD’s installation instructions for your platform.)

## Usage

1. Prepare or generate the geometry files:
   - `polygon_verts.txt`
   - `polygon_rigid.txt`

2. Place the run script (e.g. `run_polygon_alj.py`) and these geometry files in the same directory.

3. Run:

```bash
python run_polygon_alj.py
```

4. After completion, analyze `sample.gsd` with your preferred tools (e.g. `freud`, custom scripts) to:

   - measure local polygon packing,
   - identify clusters,
   - compute order parameters for assembled structures.

## Notes and context

- The box is effectively 2D (z-extent set to 0 in the initial snapshot). In HOOMD this still corresponds to a 3D box with thin z, so make sure your analysis scripts reflect that convention.
- This script is intended as a **protocol** for polygon-based patchy particle assembly, complementing:
  - MATLAB scripts that generate perfect polygon lattices,
  - init generators for linear block copolymer chains,
  - analysis scripts that characterize polygon folding and membrane-like assembly.

You can adapt the parameters (`epsilon_*`, `kappa_yukawa`, `phi_target`, `nrep`, temperature cycle) to explore different assembly regimes while keeping the overall workflow and file formats unchanged.