# Restraints Guide

This guide covers how to define and apply distance restraints in PolyzyMD simulations.

## Overview

Restraints are commonly used to:

- Keep a substrate near the enzyme active site
- Maintain specific protein-ligand contacts
- Prevent dissociation during equilibration

PolyzyMD supports several restraint types with flexible atom selection.

## Restraint Types

### Flat-Bottom Restraint

The most common restraint type. No force is applied when atoms are within the threshold distance; a harmonic force kicks in beyond it.

```
U(r) = 0                      if r < r₀
U(r) = ½k(r - r₀)²           if r ≥ r₀
```

```yaml
restraints:
  - type: "flat_bottom"
    name: "substrate_active_site"
    atom1:
      selection: "resid 77 and name OG"
    atom2:
      selection: "resname LIG and name C1"
    distance: 3.3        # Angstroms - threshold
    force_constant: 10000.0  # kJ/mol/nm²
    enabled: true
```

### Harmonic Restraint

Standard harmonic potential - always applies force to maintain target distance.

```
U(r) = ½k(r - r₀)²
```

```yaml
restraints:
  - type: "harmonic"
    name: "distance_restraint"
    atom1:
      selection: "resid 50 and name CA"
    atom2:
      selection: "resid 100 and name CA"
    distance: 10.0       # Angstroms - target distance
    force_constant: 1000.0
    enabled: true
```

### Upper Wall

Prevents distance from exceeding threshold (same as flat-bottom).

```yaml
- type: "upper_wall"
```

### Lower Wall

Prevents distance from going below threshold.

```
U(r) = ½k(r₀ - r)²           if r < r₀
U(r) = 0                      if r ≥ r₀
```

```yaml
- type: "lower_wall"
```

---

## Atom Selection Syntax

PolyzyMD uses an MDAnalysis-inspired selection syntax.

### Available Keywords

| Keyword | Indexing | Description | Example |
|---------|----------|-------------|---------|
| `resid` | 1-indexed | Residue number (matches PDB/PyMOL) | `resid 77` |
| `resname` | - | Residue name | `resname SER` |
| `name` | - | Atom name | `name OG` |
| `index` | 0-indexed | OpenMM atom index | `index 2739` |
| `pdbindex` | 1-indexed | PDB atom serial (matches PyMOL) | `pdbindex 2740` |
| `chain` | - | Chain ID | `chain A` |

### Combining Selections

Use `and` / `or` to combine:

```yaml
# Both conditions must be true
selection: "resid 77 and name OG"

# Either condition
selection: "resname SER or resname THR"

# Complex combinations
selection: "(resid 77 and name OG) or (resid 110 and name NE2)"
```

### Index Conventions

```{important}
**Understanding indexing is critical for correct restraints!**
```

| What you see | Keyword to use | Value to use |
|--------------|----------------|--------------|
| PyMOL shows residue 77 | `resid` | 77 |
| PyMOL shows atom serial 2740 | `pdbindex` | 2740 |
| OpenMM internal index 2739 | `index` | 2739 |

**Rule of thumb:** Use `resid` and `pdbindex` - they match what you see in PyMOL.

---

## Understanding Atom Ordering in PolyzyMD

PolyzyMD builds systems in a **consistent, predictable order**:

1. **Protein** - atoms numbered as in original PDB
2. **Ligand/Substrate** - atoms numbered sequentially after protein  
3. **Polymers** (if enabled) - added after ligand
4. **Solvent and ions** - added last

This means if your protein PDB has atoms 1-2722, your ligand atoms will start at **2723**.

```{important}
The `solvated_system.pdb` file generated by PolyzyMD is your reference for atom indices.
Always use this file (not your input structures) when determining `pdbindex` values.
```

---

## Workflow: Finding Atom Indices with PyMOL

### Step 1: Run Your Simulation Build

First, build your system to generate `solvated_system.pdb`:

```bash
polyzymd build -c config.yaml
```

### Step 2: Open the Solvated PDB in PyMOL

```bash
pymol /path/to/scratch/dir/solvated_system.pdb
```

### Step 3: Select and Label Your Target Atom

1. **Click on the atom** you want to restrain
2. **Right-click** the atom in the viewer
3. Navigate to: **Label → atom identifiers → ID**

This displays the atom's **ID number** directly on the structure - this is the `pdbindex` value you need.

```{tip}
You can also see the atom ID in the PyMOL console when you click on an atom:
```
You clicked /solvated_system//RBY`1/C13x
 pk1: ... id 2739
```
```

### Step 4: Write Your Selection

Use the displayed ID with `pdbindex`:

```yaml
atom2:
  selection: "resname RBY and pdbindex 2739"
  description: "Substrate ester carbonyl carbon"
```

---

## Real Example: Restraining a Docked Ligand

Here's a real example from a LipA enzyme simulation with a Resorufin Butyrate substrate.

### The Solvated PDB Structure

After building the system, the `solvated_system.pdb` shows the ligand atoms following the protein:

```
TER    2722      ASN A 181
HETATM 2723  C1x RBY B   1      72.639  64.759  41.006  1.00  0.00           C  
HETATM 2724  C2x RBY B   1      73.757  65.093  40.238  1.00  0.00           C  
HETATM 2725  C3x RBY B   1      74.974  65.339  40.886  1.00  0.00           C  
HETATM 2726  C4x RBY B   1      75.074  65.248  42.284  1.00  0.00           C  
HETATM 2727  C5x RBY B   1      73.950  64.894  43.027  1.00  0.00           C  
HETATM 2728  C6x RBY B   1      72.739  64.634  42.394  1.00  0.00           C  
HETATM 2729  N1x RBY B   1      76.267  65.488  42.979  1.00  0.00           N  
HETATM 2730  C7x RBY B   1      77.297  65.796  42.258  1.00  0.00           C  
HETATM 2731  C8x RBY B   1      78.582  66.067  42.936  1.00  0.00           C  
HETATM 2732  C9x RBY B   1      79.679  66.376  42.240  1.00  0.00           C  
HETATM 2733 C10x RBY B   1      79.654  66.462  40.768  1.00  0.00           C  
HETATM 2734  O1x RBY B   1      80.663  66.736  40.132  1.00  0.00           O  
HETATM 2735 C11x RBY B   1      78.374  66.201  40.088  1.00  0.00           C  
HETATM 2736 C12x RBY B   1      77.264  65.893  40.770  1.00  0.00           C  
HETATM 2737  O2x RBY B   1      76.065  65.659  40.104  1.00  0.00           O  
HETATM 2738  O3x RBY B   1      71.447  64.422  40.347  1.00  0.00           O  
HETATM 2739 C13x RBY B   1      70.808  63.259  40.777  1.00  0.00           C   <- Ester carbonyl
HETATM 2740  O4x RBY B   1      69.599  63.071  40.828  1.00  0.00           O  
HETATM 2741 C14x RBY B   1      71.835  62.259  41.248  1.00  0.00           C  
...
TER    2757      RBY B   1
HETATM 2758  C1x SBM C   1      80.069  63.041  53.056  1.00  0.00           C   <- Polymers start
```

### The Configuration

To restrain atom **C13x** (the ester carbonyl carbon, pdbindex 2739) near the catalytic serine:

```yaml
restraints:
  - type: "flat_bottom"
    name: "substrate_active_site"
    atom1:
      selection: "resid 77 and name OG"
      description: "Catalytic serine oxygen"
    atom2:
      selection: "resname RBY and pdbindex 2739"
      description: "Substrate ester carbonyl carbon"
    distance: 10.0        # Angstroms - no force applied within this distance
    force_constant: 10000.0
    enabled: true
```

### What You'll See When Running

When the simulation runs, you'll see confirmation that the restraint was applied:

```
Applying 1 restraint(s)...
  - substrate_active_site: atom 1192 <-> atom 2738 (type=flat_bottom, d=10.0 Å, k=10000.0 kJ/mol/nm²)
Successfully applied 1 restraint(s)
Running energy minimization...
```

```{note}
The logged atom indices are **0-indexed** (OpenMM internal). They will be one less than the `pdbindex` values.
```

---

## Common Restraint Patterns

### Substrate in Active Site

Keep substrate within bonding distance of catalytic residue:

```yaml
restraints:
  - type: "flat_bottom"
    name: "substrate_catalytic"
    atom1:
      selection: "resid 77 and name OG"
      description: "Catalytic serine"
    atom2:
      selection: "resname LIG and name C1"
      description: "Substrate carbonyl carbon"
    distance: 3.5
    force_constant: 10000.0
    enabled: true
```

### Multiple Active Site Contacts

```yaml
restraints:
  - type: "flat_bottom"
    name: "contact_1"
    atom1:
      selection: "resid 77 and name OG"
    atom2:
      selection: "resname LIG and name C1"
    distance: 3.5
    force_constant: 10000.0
    enabled: true
    
  - type: "flat_bottom"
    name: "contact_2"
    atom1:
      selection: "resid 110 and name NE2"
    atom2:
      selection: "resname LIG and name O2"
    distance: 3.5
    force_constant: 5000.0
    enabled: true
    
  - type: "flat_bottom"
    name: "contact_3"
    atom1:
      selection: "resid 12 and name N"
    atom2:
      selection: "resname LIG and name O1"
    distance: 4.0
    force_constant: 5000.0
    enabled: false  # Disabled for now
```

### Protein-Protein Distance

Monitor or restrain distance between protein regions:

```yaml
restraints:
  - type: "harmonic"
    name: "domain_distance"
    atom1:
      selection: "resid 50 and name CA"
    atom2:
      selection: "resid 150 and name CA"
    distance: 25.0
    force_constant: 100.0  # Weak restraint
    enabled: true
```

---

## Force Constant Guidelines

| Use Case | Force Constant (kJ/mol/nm²) |
|----------|---------------------------|
| Strong restraint (prevent dissociation) | 10000 - 50000 |
| Moderate restraint | 1000 - 5000 |
| Weak/guiding restraint | 100 - 500 |
| Very weak (biasing) | 10 - 50 |

```{note}
Start with stronger restraints during equilibration, then reduce or disable for production if desired.
```

---

## Troubleshooting

### "No atoms match selection"

```
ValueError: No atoms match selection: 'resid 77 and name OG'
```

**Causes:**
1. Residue numbering doesn't match (check PDB)
2. Atom name is different (check PDB atom names)
3. Typo in selection string

**Solution:** Open PDB in PyMOL and verify exact residue numbers and atom names.

### "Requires exactly one atom per selection"

```
ValueError: Restraint 'my_restraint' requires exactly one atom per selection.
Got 5 for atom1, 1 for atom2
```

**Cause:** Selection matches multiple atoms.

**Solution:** Make selection more specific:

```yaml
# Too broad - matches all atoms in residue
selection: "resid 77"

# Better - matches specific atom
selection: "resid 77 and name OG"

# Most specific - use atom index
selection: "resid 77 and pdbindex 1193"
```

### Restraint Not Taking Effect

1. Check `enabled: true` is set
2. Verify distance units (should be Angstroms)
3. Check force constant is reasonable (not too small)
4. Verify atom selections resolve correctly

---

## Disabling Restraints

### In Configuration

```yaml
restraints:
  - type: "flat_bottom"
    name: "my_restraint"
    # ... other settings ...
    enabled: false  # Disabled
```

### No Restraints at All

```yaml
restraints: []
```

Or simply omit the `restraints` section entirely.