KO | EN

GitLite

Minimalist Open Source Repository Search Engine

gitlite — search

#dsh #dsh-plugin #ai-agents #developer-tools #codex #claude-code #mcp #rust #deepseek-harness #typescript #ai-agent #cli

pAPRika

★ 36 Open GitHub ↗

Advanced toolkit for binding free energy calculations

// repository documentation

Was this content helpful?

★ 0 (0 ratings)

# pAPRika
pAPRika is a toolkit for setting up, running, and analyzing free energy molecular dynamics simulations.

# Badges
[![Build Status](https://github.com/GilsonLabUCSD/pAPRika/workflows/ci/badge.svg)](https://github.com/GilsonLabUCSD/pAPRika/actions?query=branch%3Amaster+workflow%3Aci)
[![codecov](https://codecov.io/gh/GilsonLabUCSD/pAPRika/branch/master/graph/badge.svg?token=ZTjfSd3WK2)](https://codecov.io/gh/GilsonLabUCSD/pAPRika)
[![CodeQL](https://github.com/GilsonLabUCSD/pAPRika/actions/workflows/codeql.yml/badge.svg)](https://github.com/GilsonLabUCSD/pAPRika/actions/workflows/codeql.yml)

[![Last release tag](https://img.shields.io/github/release/GilsonLabUCSD/paprika.svg)](https://github.com/GilsonLabUCSD/paprika/releases)
[![Commits since release](https://img.shields.io/github/commits-since/GilsonLabUCSD/pAPRika/latest.svg)](https://github.com/GilsonLabUCSD/pAPRika/releases/tag/latest)
[![Releases](https://img.shields.io/badge/obtain-latest-green.svg)](https://github.com/GilsonLabUCSD/pAPRika/releases) [![Conda](https://img.shields.io/conda/v/conda-forge/paprika.svg)](https://anaconda.org/conda-forge/paprika) [![Last updated](https://anaconda.org/conda-forge/paprika/badges/latest_release_relative_date.svg)](https://anaconda.org/conda-forge/paprika) [![Anaconda Cloud downloads](https://anaconda.org/conda-forge/paprika/badges/downloads.svg)](https://anaconda.org/conda-forge/paprika)

[![Documentation Status](https://readthedocs.org/projects/paprika/badge/?version=latest)](https://paprika.readthedocs.io/en/latest/?badge=latest) [![platforms](https://img.shields.io/badge/Platforms-Linux%2C%20MacOS-orange.svg)](https://paprika.readthedocs.io/en/latest/install.html) [![python](https://img.shields.io/badge/python-3.8%2C%203.9%2C%203.10-blue.svg)](https://paprika.readthedocs.io/en/latest/install.html) [![License](https://img.shields.io/badge/License-BSD%203--Clause-green.svg)](https://opensource.org/licenses/BSD-3-Clause)

# Contributors
- Niel Henriksen (UCSD, Atomwise Inc.)
- David Slochower (UCSD, Vertex Pharmaceuticals)
- Simon Boothroyd (Sloan Kettering Institute, Open Force Field, XtalPi Inc.)
- Jeff Setiadi (UCSD)
- Willa Wang (UCSD)
- Jake Anderson (UCSD)

# Installation

We recommend installing pAPRika in a fresh `conda` environment if possible. There are three ways to install this package:

1. The latest release on `conda-forge`:
    1. `conda install -c conda-forge paprika`
    2. To use all features of pAPRika, you must either have [AmberTools](http://ambermd.org/AmberTools.php) in your `$PATH` or separately install AmberTools with `conda install -c conda-forge ambertools`.

2. The master branch on GitHub:
    1. Clone this `git` repository, then inside the `paprika` directory:
    2. Change the `name` field in `devtools/conda-envs/test_env.yaml` to be `paprika`.
    3. Create the environment: `conda env create -f devtools/conda-envs/test_env.yaml`.
    4. Activate the environment: `conda activate paprika`
    5. Install `paprika` in the environment: `pip install .`
    
3. The latest release on GitHub:
    1. Download [the latest release](https://github.com/GilsonLabUCSD/pAPRika/releases), extract it, and change to the `paprika` directory:
    2. Change the `name` field in `devtools/conda-envs/test_env.yaml` to be `paprika`.
    3. Create the environment: `conda env create -f devtools/conda-envs/test_env.yaml`.
    4. Activate the environment: `conda activate paprika`
    5. Install `paprika` in the environment: `pip install .`

# Example

In this example, we will configure and simulate butane (BUT) as a guest molecule for the host [cucurbit[6]uril](https://en.wikipedia.org/wiki/Cucurbituril) (CB6). CBs are rigid, symmetric, cyclic host molecules with oxygen atoms around the portal edge of the cavity. We will run the simulation in implicit solvent, using the [Generalized-Born](https://en.wikipedia.org/wiki/Implicit_solvation#Generalized_Born) model, for speed and simplicity, using AMBER. This tutorial assumes familiarity with basic MD procedures.

The `cb6-but` folder referred to below are located in the directory `paprika/data`.

pAPRika utilizes `Pint` units through the `openff-units` wrapper. Any float value that is passed to a function (mostly in the `build`, `restraints` and `analysis` modules) for masses, distances, angles, energy and force constants will be transformed to a `Pint` quantity based on the default units below (AMBER convention):
* mass: Dalton
* distance: angstrom
* angle: degrees
* energy: kcal/mol
* force constants: kcal/mol/angstrom^2 (distances) or kcal/mol/radians^2 (angles)

If a `Pint` quantity is specified then the units will be converted accordingly to match the list above.

```python
from openff.units import unit
```

## Initial Setup

The very first step in the calculation is creating a coordinate file for the bound host-guest complex (we usually use PDB format for this part because it works well with `tleap`). This does not have to perfectly match the bound state by any means (we will later minimize and equilibrate), but this should be a reasonable illustration of the bound complex. This file can be created by hand (in a program like Chimera, VMD, or PyMOL) or by docking the guest into the host cavity (with MOE, AutoDock, DOCK, ...).

In this example, the host-guest system coordinates (`cb6-but.pdb`) is located in `paprika/data`. Below is an image of the host-guest molecules.

![](https://raw.githubusercontent.com/GilsonLabUCSD/pAPRika/master/docs/tutorials/images/cb6-but.png)

In addition to the coordinate file, we need separate `mol2` files for the host and guest molecule that contain the partial atomic charges. For cyclic hosts like CB6, you can specify either a single residue for the entire molecule (this is what we do here) or you can provide coordinates and charges for a single monomer and have `tleap` build the structure (this is a little more tricky).

The `mol2` files that we use here (`cb6.mol2` and `but.mol2`) were created by running `antechamber -fi pdb -fo mol2 -i <pdb> -o <mol2> -c bcc` (I added `pl 10` for the host, which reduces the number of paths that `antechamber` takes to traverse the host; see `antechamber` help for more information).

### Create AMBER coordinate (`.rst7` or `.inpcrd`) and parameter files (`.prmtop` or `.topo`) for the host-guest complex

In this example, we will use GAFF parameters for both the host and guest. For the host, `parmchk2` has identified two parameters that are missing from GAFF and added the most similar ones into the supplementary `cb6.frcmod` file.

```python
from paprika.build.system import TLeap
```

```python
system = TLeap()
system.output_path = "cb6-but"
system.pbc_type = None
system.neutralize = False

system.template_lines = [
    "source leaprc.gaff",
    "loadamberparams cb6.frcmod",
    "CB6 = loadmol2 cb6.mol2",
    "BUT = loadmol2 but.mol2",
    "model = loadpdb cb6-but.pdb",
    "check model",
    "savepdb model vac.pdb",
    "saveamberparm model vac.prmtop vac.rst7"
]

system.build()
```

After running `tleap`, it is always a good idea to check `leap.log`. `pAPRika` does some automated checking of the ouptut, but sometimes things slip through. (By default, `tleap` will append to `leap.log`.)

## Prepare the complex for an APR calculation

Now we are ready to prepare the complex for the attach-pull-release calculation. This involves:

- Aligning the structure so the guest can be pulled along the *z*-axis, and
- Adding dummy atoms that are used to orient the host and guest.

To access the host-guest structure in Python, we use the ParmEd `Structure` class. So we start by loading the vacuum model that we just created. Then, we need to define two atoms on the guest that are placed along the *z*-axis. These should be heavy atoms on either end of the guest, the second atom leading the pulling.

These same atoms will be used later for the restraints, so I will name them `G1` and `G2`, using AMBER selection syntax.

### Align the pulling axis

```python
import parmed as pmd

```

```python
structure = pmd.load_file("cb6-but/vac.prmtop", "cb6-but/vac.rst7", structure=True)
```

```python

from paprika.build import align

```

```python
G1 = ":BUT@C"
G2 = ":BUT@C3"
```

```python
aligned_structure = align.zalign(structure, G1, G2)
aligned_structure.save("cb6-but/aligned.prmtop", overwrite=True)
aligned_structure.save("cb6-but/aligned.rst7", overwrite=True)
```

Here, the origin is shown as a grey sphere, with the *z*-axis drawn as a blue arrow. The coordinates used for this example were already aligned, so Python warns that the cross product is zero, but this won't be the case in general.

![](https://raw.githubusercontent.com/GilsonLabUCSD/pAPRika/master/docs/tutorials/images/aligned.png)

Next, we add the dummy atoms. The dummy atoms will be fixed in place during the simulation and are used to orient the host and guest in the lab frame. The dummy atoms are placed along the \ axis, behind the host. The dummy atoms are used in distance, angle, and torsion restraints and therefore, the exact positioning of these atoms affects the value of those restraints. For a typical host-guest system, like the one here, we generally place the first dummy atom 6 Angstroms behind the origin, the second dummy atom 9 Angstroms behind the origin, and the third dummy atom 11.2 Angstroms behind the origin and offset about 2.2 Angstroms along the *y*-axis. After we add restraints, the positioning of the dummy atoms should be more clear.

Note, these dummy atoms do not interact with the other atoms in the system, and therefore, there is no problem placing them near host atoms.

### Add dummy atoms

```python
from paprika.build import dummy

```

```python
structure = pmd.load_file("cb6-but/aligned.prmtop", "cb6-but/aligned.rst7", structure=True)
```

```python
structure = dummy.add_dummy(structure, residue_name="DM1", z=-6.0 * unit.angstrom)
structure = dummy.add_dummy(structure, residue_name="DM2", z=-9.0 * unit.angstrom)
structure = dummy.add_dummy(structure, residue_name="DM3", z=-11.2 * unit.angstrom, y=2.2 * unit.angstrom)
```

```python
structure.save("cb6-but/aligned_with_dummy.prmtop", overwrite=True)
structure.save("cb6-but/aligned_with_dummy.rst7", overwrite=True)
structure.save("cb6-but/aligned_with_dummy.pdb", overwrite=True)
```

![](https://raw.githubusercontent.com/GilsonLabUCSD/pAPRika/master/docs/tutorials/images/aligned_with_dummy.png)

When we solvate the system in `tleap`, we will need `frcmod` files for the dummy atoms (otherwise `tleap` will use GAFF parameters and the dummy atoms will *not* be non-interacting). There is a convenient method in `paprika` to write a `frcmod` file that only contains a `MASS` section. For convenience, I am also going to write `mol2` files for each of the dummy atoms. This makes it easy to build up the system, piece-by-piece, if we have a separate `mol2` file for each component of the system: host, guest, dummy atoms.

In principle, an APR calculate can be completed without dummy atoms, by simply lengthening the distance between the host and guest, but the addition of dummy atoms permits an easier way to think about dissociating the guest from the host. (Also, in the absence of dummy atoms, it is challenging to pull the guest straight out of the cavity without adding additional restraints.)

```python
dummy.write_dummy_frcmod(filepath="cb6-but/dummy.frcmod")
dummy.write_dummy_mol2(residue_name="DM1", filepath="cb6-but/dm1.mol2")
dummy.write_dummy_mol2(residue_name="DM2", filepath="cb6-but/dm2.mol2")
dummy.write_dummy_mol2(residue_name="DM3", filepath="cb6-but/dm3.mol2")
```

Now all the pieces are in place to build the system for an APR calculation.

### Put it all together

```python
system = TLeap()
system.output_path = "cb6-but"
system.pbc_type = None
system.neutralize = False

system.template_lines = [
    "source leaprc.gaff",
    "loadamberparams cb6.frcmod",
    "loadamberparams dummy.frcmod",
    "CB6 = loadmol2 cb6.mol2",
    "BUT = loadmol2 but.mol2",
    "DM1 = loadmol2 dm1.mol2",
    "DM2 = loadmol2 dm2.mol2",
    "DM3 = loadmol2 dm3.mol2",
    "model = loadpdb aligned_with_dummy.pdb",
    "check model",
    "savepdb model cb6-but-dum.pdb",
    "saveamberparm model cb6-but-dum.prmtop cb6-but-dum.rst7"
]

system.build()
```

Now we have AMBER coordinates and parameters for the `cb6-but` system with dummy atoms in the appropriate place and with the proper "dummy" parametesr.

![](https://raw.githubusercontent.com/GilsonLabUCSD/pAPRika/master/docs/tutorials/images/anchor-atoms-diagram2.png)

## Determine the number of windows

Before we add the restraints, it is helpful to set the λ fractions that control the strength of the force constants during attach and release, and to define the distances for the pulling phase.

The attach fractions go from 0 to 1 and we place more points at the bottom of the range to sample the curvature of ∂U/∂λ. Next, we generally apply a distance restraint until the guest is about 18 Å away from the host, in increments of 0.4 Å. This distance should be at least twice the Lennard-Jones cutoff in the system. These values have worked well for us, but this is one aspect that should be carefully checked for new systems.

```python
attach_string = "0.00 0.40 0.80 1.60 2.40 4.00 5.50 8.65 11.80 18.10 24.40 37.00 49.60 74.80 100.00"
attach_fractions = [float(i) / 100 for i in attach_string.split()]
```

```python
import numpy as np
initial_distance = 6.0
pull_distances = np.arange(0.0 + initial_distance, 18.0 + initial_distance, 1.0)
```

These values will be used to measure distance relative to the first dummy atom, hence the addition of `6.00`.

```python
release_fractions = []
```

Later, I will explain below why there are no release windows in this calculation.

```python
windows = [len(attach_fractions), len(pull_distances), len(release_fractions)]
print(f"There are {windows} windows in this attach-pull-release calculation.")
```

There are [15, 18, 0] windows in this attach-pull-release calculation.

Alternatively, we could specify the number of windows for each phase and the force constants and targets will be linearly interpolated. Other ways of specifying these values are documented in the code.

## Add restraints using pAPRika

`pAPRika` supports four different types of restraints:

- Static restraints: these six restraints keep the host and in the proper orientation during the simulation (necessary),
- Guest restraints: these restraints pull the guest away from the host along the *z*-axis (necessary),
- Conformational restraints: these restraints alter the conformational sampling of the host molecule (optional), and
- Wall restraints: these restraints help define the bound state of the guest (optional).

More information on these restraints can be found in:

Henriksen, N.M., Fenley, A.T., and Gilson, M.K. (2015). Computational Calorimetry: High-Precision Calculation of Host-Guest Binding Thermodynamics. J. Chem. Theory Comput. 11, 4377–4394. [DOI](http://doi.org/10.1021/acs.jctc.5b00405)

In this example, I will show how to setup the static restraints and the guest restraints.

We have already added the dummy atoms and we have already defined the guest anchor atoms. Now we need to define the host anchor atoms (H1, H2, and H3) in the above diagram. The host anchors should be heavy atoms distributed around the cavity (and around the pulling axis). One caveat is that the host anchors should be rigid relative to each other, so conformational restraints do not shift the alignment of the pulling axis relative to the solvation box. For CB6, I have chosen carbons around the central ridge.

```python
H1 = ":CB6@C"
H2 = ":CB6@C31"
H3 = ":CB6@C18"
```

I'll also make a shorthand for the dummy atoms.

```python
D1 = ":DM1"
D2 = ":DM2"
D3 = ":DM3"
```

### Static restraints

These harmonic restraints are constant throughout the entire simulation. These restraints are used to control the distances and angles between the host and guest relative to the dummy atom. We have created a special class for these restrains, `static_DAT_restraint`, that uses the initial value as the restraint target (this is why the starting structure should be a reasonable facsimile of the bound state).

Note that these restraints are not "attached" and they don't need to be "released" -- their force constants do not change in magnitude.

The first three static restraints affect the translational distance, angle, and torsion angle between the host and the dummy atoms. These control the position of the host, via the first anchor atom, from moving relative to the dummy atoms.

There is no *correct* value for the force constants. From experience, we know that a distance force constant of 5.0 kcal/mol/Angstrom^2 won't nail down the host and yet it also won't wander away. Likewise, we have had good results using 100.0 kcal/mol/radian^2 for the angle force constant.

![](https://raw.githubusercontent.com/GilsonLabUCSD/pAPRika/master/docs/tutorials/images/static-restraints-1.png)

```python
from paprika import restraints
```

```python
k_dist = 5.0 * unit.kcal / unit.mole / unit.angstrom ** 2
k_angle = 100.0 * unit.kcal / unit.mole / unit.radians ** 2
static_restraints = []
```

```python
structure = pmd.load_file("cb6-but/cb6-but-dum.prmtop",
                          "cb6-but/cb6-but-dum.rst7",
                          structure = True
                          )
```

```python
r = restraints.static_DAT_restraint(restraint_mask_list = [D1, H1],
                                    num_window_list = windows,
                                    ref_structure = structure,
                                    force_constant = k_dist,
                                    amber_index=True)