pAPRika tutorial 6 - APR/Gromacs with Plumed restraints

In this tutorial, we will perform APR calculations for the butane (BUT)–cucurbit[6]uril (CB6) host-guest system. This tutorial is similar to tutorial 5 where Plumed will be used for the restraints but we will be using the GROMACS MD engine instead. One other difference in this tutorial is that we will be performing the simulation with periodic boundary condition (PBC) and explicit water molecules. This is because GROMACS does not support implicit solvent calculations.

Initialize

Before you start

We will run the simulations in this tutorial using GROMACS and Plumed. GROMACS will need to be installed before you can run the simulation part of this tutorial. Plumed, however, should be installed in your conda environment if you installed pAPRika though the conda route. To use Plumed we first need to make sure the PLUMED_KERNEL environment variable is loaded (the library is called libplumedKernel.so). pAPRika should load the Plumed kernel automatically but let’s make sure it is loaded and run the cell below.

import os
'PLUMED_KERNEL' in os.environ.keys()

True

If it does not exists we will load the environment in this Jupyter Notebook. Since Plumed is installed through conda the kernel will be located in your conda environment library folder. If you are running this on a MacOS, replace the kernel library in the cell below to libplumedKernel.dylib. If you compiled Plumed yourself and then you will need to change the path below.

os.environ['PLUMED_KERNEL'] = f"{os.environ['CONDA_PREFIX']}/lib/libplumedKernel.so"

For running Plumed outside of this notebook it might be better to export the PLUMED_KERNEL variable into your .bashrc file.

Note: we can run Plumed with GROMACS versions 2020.2 but the user is required to patch the source code. Unlike AMBER and LAMMPS, GROMACS does not come with a patched code out-of-the-box and the conda-forge repository also do not contain a patched version of GROMACS. See the Plumed documentation for more details on how to patch and recompile the GROMACS code https://www.plumed.org/doc-v2.6/user-doc/html/index.html.

Note: We will be running the CB6-BUT host-guest system in explicit solvent because GROMACS does not support implicit solvent calculation. Thus, we recommend running the simulation on a GPU.

🔵 Since we have a prepared the host-guest-dummy setup from the first tutorial, we will skip the initial tleap steps and go right into initializing the restraints. If you haven’t completed tutorial 1 please go back and run the notebook to generate the starting structure.

Define names

We will store the files created in this tutorial in a folder called gromacs-plumed so we don’t mix files with the previous tutorial.

base_name = "cb6-but-dum"
work_dir = "gromacs-plumed"
complex_dir = "complex"
data = "../../../paprika/data/cb6-but"

Configure APR Restraints

Define anchor atoms

See tutorial 1 for the choice of selection

# Guest atoms
G1 = ":BUT@C"
G2 = ":BUT@C3"

# Host atoms
H1 = ":CB6@C"
H2 = ":CB6@C31"
H3 = ":CB6@C18"

# Dummy atoms
D1 = ":DM1"
D2 = ":DM2"
D3 = ":DM3"

Determine the number of windows

Before we add the restraints, it is helpful to set the \(\lambda\) 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 \(dU/d \lambda\). Next, we generally apply a distance restraint until the guest is ~18 Angstroms away from the host, in increments of 0.4 Angstroms. 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.

import numpy as np

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()]

initial_distance = 6.0
pull_distances = np.arange(0.0 + initial_distance, 18.0 + initial_distance, 1.0)

release_fractions = []

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.

Load structure

import parmed as pmd

• Load complex structure

structure = pmd.load_file(
    os.path.join(complex_dir, f"{base_name}.prmtop"),
    os.path.join(complex_dir, f"{base_name}.rst7"),
    structure = True,
)

Host Static Restraints

See tutorial 1 for an explanation of the static restraints.

import paprika.restraints as restraints

static_restraints = []

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

static_restraints.append(r)

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

static_restraints.append(r)

r = restraints.static_DAT_restraint(restraint_mask_list = [D3, D2, D1, H1],
                                    num_window_list = windows,
                                    ref_structure = structure,
                                    force_constant = 100.0,
                                    amber_index=True)

static_restraints.append(r)

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

static_restraints.append(r)

r = restraints.static_DAT_restraint(restraint_mask_list = [D2, D1, H1, H2],
                                    num_window_list = windows,
                                    ref_structure = structure,
                                    force_constant = 100.0,
                                    amber_index=True)

static_restraints.append(r)

r = restraints.static_DAT_restraint(restraint_mask_list = [D1, H1, H2, H3],
                                    num_window_list = windows,
                                    ref_structure = structure,
                                    force_constant = 100.0,
                                    amber_index=True)

static_restraints.append(r)

Guest translational and rotational restraints

See tutorial 1 for an explanation of the guest restraints

guest_restraints = []

r = restraints.DAT_restraint()
r.mask1 = D1
r.mask2 = G1
r.topology = structure
r.auto_apr = True
r.continuous_apr = True
r.amber_index = True

r.attach["target"] = pull_distances[0]          # Angstroms
r.attach["fraction_list"] = attach_fractions
r.attach["fc_final"] = 5.0                      # kcal/mol/Angstroms**2

r.pull["target_final"] = 24.0                   # Angstroms
r.pull["num_windows"] = windows[1]

r.initialize()
guest_restraints.append(r)

r = restraints.DAT_restraint()
r.mask1 = D2
r.mask2 = D1
r.mask3 = G1
r.topology = structure
r.auto_apr = True
r.continuous_apr = True
r.amber_index = True

r.attach["target"] = 180.0                      # Degrees
r.attach["fraction_list"] = attach_fractions
r.attach["fc_final"] = 100.0                    # kcal/mol/radian**2

r.pull["target_final"] = 180.0                  # Degrees
r.pull["num_windows"] = windows[1]

r.initialize()
guest_restraints.append(r)

r = restraints.DAT_restraint()
r.mask1 = D1
r.mask2 = G1
r.mask3 = G2
r.topology = structure
r.auto_apr = True
r.continuous_apr = True
r.amber_index = True

r.attach["target"] = 180.0                      # Degrees
r.attach["fraction_list"] = attach_fractions
r.attach["fc_final"] = 100.0                    # kcal/mol/radian**2

r.pull["target_final"] = 180.0                  # Degrees
r.pull["num_windows"] = windows[1]

r.initialize()
guest_restraints.append(r)

Create APR windows

We use the guest restraints to create a list of windows with the appropriate names and then create the directories.

from paprika.restraints.utils import create_window_list

window_list = create_window_list(guest_restraints)

if not os.path.isdir(work_dir):
    os.makedirs(work_dir)

for window in window_list:
    folder = os.path.join(work_dir, window)
    if not os.path.isdir(folder):
        os.makedirs(os.path.join(work_dir, window))

Write APR restraints to Plumed format

In this section we create an instance of Plumed() from paprika.restraints.plumed, which is a class to generate the Plumed restraint files. We need to specify the list of restraints used throughout the APR calculations and the corresponding windows list. In this tutorial we will print the host static restraints and the guest restraints. The Plumed class includes a method to add restraints to dummy atoms (add_dummy_atoms_to_file). We will add the positional restraints on the dummy atoms after solvation of the host-guest complex.

Note: be careful when specifying the force constants in DAT_restraints. We follow the AMBER (and CHARMM) convention where the force constant is already multiplied by a factor of 1/2 but Plumed requires the user to specify the force constant without this factor, i.e.

\[\begin{split}\begin{align} U_{amber} &= K_{amber} (r-r_{0})^2 \\ U_{plumed} &= \frac{1}{2} k_{plumed} (r - r_{0})^2 \end{align}\end{split}\]

thus \(k_{plumed} = 2 \times K_{amber}\). If AMBER force constants was used in generating the DAT_restraints (the case in this tutorial) we need to set the variable uses_legacy_k to True (this is on by default).

from paprika.restraints.plumed import Plumed

restraints_list = (static_restraints + guest_restraints)

plumed = Plumed()
plumed.file_name = 'plumed.dat'
plumed.path = work_dir
plumed.window_list = window_list
plumed.restraint_list = restraints_list
plumed.uses_legacy_k = True

plumed.dump_to_file()

Prepare host-guest system

Translate guest molecule

For the attach windows, we will use the initial, bound coordinates for the host-guest complex. Only the force constants change during this phase, so a single set of coordinates is sufficient. For the pull windows, we will translate the guest to the target value of the restraint before solvation, and for the release windows, we will use the coordinates from the final pull window.

import shutil

for window in window_list:
    if window[0] == "a":
        shutil.copy(os.path.join(complex_dir, f"{base_name}.prmtop"),
                    os.path.join(work_dir, window, f"{base_name}.prmtop"))
        shutil.copy(os.path.join(complex_dir, f"{base_name}.rst7"),
                    os.path.join(work_dir, window, f"{base_name}.rst7"))

    elif window[0] == "p":
        structure = pmd.load_file(
            os.path.join(complex_dir, f"{base_name}.prmtop"),
            os.path.join(complex_dir, f"{base_name}.rst7"),
            structure = True
        )
        target_difference = guest_restraints[0].phase['pull']['targets'][int(window[1:])] -\
                            guest_restraints[0].pull['target_initial']
        print(f"In window {window} we will translate the guest {target_difference.magnitude:0.1f}.")

        for atom in structure.atoms:
            if atom.residue.name == "BUT":
                atom.xz += target_difference.magnitude

        structure.save(os.path.join(work_dir, window, f"{base_name}.prmtop"), overwrite=True)
        structure.save(os.path.join(work_dir, window, f"{base_name}.rst7"), overwrite=True)

In window p000 we will translate the guest 0.0 Angstroms.
In window p001 we will translate the guest 1.1 Angstroms.
In window p002 we will translate the guest 2.1 Angstroms.
In window p003 we will translate the guest 3.2 Angstroms.
In window p004 we will translate the guest 4.2 Angstroms.
In window p005 we will translate the guest 5.3 Angstroms.
In window p006 we will translate the guest 6.4 Angstroms.
In window p007 we will translate the guest 7.4 Angstroms.
In window p008 we will translate the guest 8.5 Angstroms.
In window p009 we will translate the guest 9.5 Angstroms.
In window p010 we will translate the guest 10.6 Angstroms.
In window p011 we will translate the guest 11.6 Angstroms.
In window p012 we will translate the guest 12.7 Angstroms.
In window p013 we will translate the guest 13.8 Angstroms.
In window p014 we will translate the guest 14.8 Angstroms.
In window p015 we will translate the guest 15.9 Angstroms.
In window p016 we will translate the guest 16.9 Angstroms.
In window p017 we will translate the guest 18.0 Angstroms.

Solvate the host-guest system

Here we will solvate the host-guest system in each window with 2000 water molecules. Since Gromacs cannot read in Amber topology and coordinates files we will have to convert them to the appropriate formats. TLeap has a built-in function called TLeap.convert_to_gromacs(), which will convert the Amber files to Gromacs. We also need to add positional restraints on dummy atoms to the plumed.dat file and this is done by invoking the function plumed.add_dummy_atoms_to_file

from paprika.build.system import TLeap

for window in window_list:
    print(f"Solvating system in {window}.")

    structure = pmd.load_file(
        os.path.join(work_dir, window, f"{base_name}.prmtop"),
        os.path.join(work_dir, window, f"{base_name}.rst7"),
        structure=True,
    )

    if not os.path.exists(os.path.join(work_dir, window, f"{base_name}.pdb")):
        structure.save(os.path.join(work_dir, window, f"{base_name}.pdb"))

    system = TLeap()
    system.output_path = os.path.join(work_dir, window)
    system.output_prefix = f"{base_name}-sol"

    system.target_waters = 2000
    system.neutralize = True
    system.add_ions = ["Na+", 6, "Cl-", 6]
    system.pbc_type = "rectangular"
    system.template_lines = [
        "source leaprc.gaff",
        "source leaprc.water.tip3p",
        f"loadamberparams {data}/cb6.frcmod",
        "loadamberparams ../../complex/dummy.frcmod",
        f"CB6 = loadmol2 {data}/cb6.mol2",
        f"BUT = loadmol2 {data}/but.mol2",
        "DM1 = loadmol2 ../../complex/dm1.mol2",
        "DM2 = loadmol2 ../../complex/dm2.mol2",
        "DM3 = loadmol2 ../../complex/dm3.mol2",
        "model = loadpdb cb6-but-dum.pdb",
    ]
    system.build()

    # Convert Amber to Gromacs
    system.convert_to_gromacs(overwrite=True)

    # Add dummy atom restraints - load in the Gromacs files instead
    # of the amber files because the coordinates a shifted.
    structure = pmd.load_file(
        os.path.join(work_dir, window, f"{base_name}-sol.top"),
        xyz=os.path.join(work_dir, window, f"{base_name}-sol.gro"),
        structure=True
    )
    plumed.add_dummy_atoms_to_file(structure,window=window)

Solvating system in a000.
Solvating system in a001.
Solvating system in a002.
Solvating system in a003.
Solvating system in a004.
Solvating system in a005.
Solvating system in a006.
Solvating system in a007.
Solvating system in a008.
Solvating system in a009.
Solvating system in a010.
Solvating system in a011.
Solvating system in a012.
Solvating system in a013.
Solvating system in p000.
Solvating system in p001.
Solvating system in p002.
Solvating system in p003.
Solvating system in p004.
Solvating system in p005.
Solvating system in p006.
Solvating system in p007.
Solvating system in p008.
Solvating system in p009.
Solvating system in p010.
Solvating system in p011.
Solvating system in p012.
Solvating system in p013.
Solvating system in p014.
Solvating system in p015.
Solvating system in p016.
Solvating system in p017.

Simulation

Now that we have a solvated structure and a plumed.dat file with dummy atoms position restraints, we have everything ready to go. pAPRika has a GROMACS wrapper module that can help setting default parameters for the simulation. There are some high level options that we set directly, like simulation.path, and then we call the function config_pbc_min() to setup reasonable default simulation parameters for a minimization and production run. See the pAPRika documentation for more details.

GROMACS, unlike AMBER and OpenMM, have a higher dependency on CPUs, so more than one processor for every GPU is preferred to deliver better speed. Setting the number of processors requires using different keyword depending on how GROMACS was compiled. The cells below were run using GROMACS that is compiled with MPI and hence the suffix, i.e., gmx_mpi (for a non-MPI version, the executable is usually called gmx). We set the number of processors with simulation.n_threads = 8 and the GPU device with simulation.gpu_devices = 0.

Note: For host-guest systems I find that the non-MPI version is considerably much faster. Also, 12 CPUs with a single GPU seems to be a good combination with a non-MPI version.

By default, pAPRika will use -ntomp to specify the number of CPUs when executing mdrun if the executable is gmx_mpi and -nt for gmx. If these default settings do not work with the GROMACS currently installed then we can use simulation.custom_mdrun_command to override the default, e.g.,

simulation.custom_mdrun_command = "-ntmpi 8 -gpu_id 1"

or

simulation.custom_mdrun_command = "mpirun -np 4 gmx_mpi mdrun -ntomp 8 -gpu_id 4"

Note: We do not need to run grompp beforehand because the simulation wrapper in pAPRika takes care of this automatically when we call the simulation.run() method. We also set simulation.grompp_maxwarn = 5 to ignore warnings about the total charge of the system (the total charge is slightly off from neutral).

Note: as explained at the start of this tutorial, make sure that the PLUMED_KERNEL environment variable is set otherwise the simulation will fail to run.

from paprika.simulate import GROMACS

Initialize logger

import logging
from importlib import reload
reload(logging)

logger = logging.getLogger()
logging.basicConfig(
    format='%(asctime)s %(message)s',
    datefmt='%Y-%m-%d %I:%M:%S %p',
    level=logging.INFO
)

Energy Minimization

for window in window_list:
    simulation = GROMACS()
    simulation.executable = "gmx_mpi"
    simulation.n_threads = 8
    simulation.gpu_devices = 0
#     simulation.custom_mdrun_command = "-ntomp 8 -gpu_id 6"
    simulation.path = os.path.join(work_dir, window)
    simulation.maxwarn = 5

    simulation.prefix = "minimize"
    simulation.topology = f"{base_name}-sol.top"
    simulation.coordinates = f"{base_name}-sol.gro"
    simulation.plumed_file = "plumed.dat"

    simulation.config_pbc_min()

    logging.info(f"Running minimization in window {window}...")
    simulation.run(overwrite=True)

2020-09-25 05:22:17 PM Running minimization in window a000...
2020-09-25 05:22:36 PM Running minimization in window a001...
2020-09-25 05:22:56 PM Running minimization in window a002...
2020-09-25 05:23:15 PM Running minimization in window a003...
2020-09-25 05:23:37 PM Running minimization in window a004...
2020-09-25 05:23:56 PM Running minimization in window a005...
2020-09-25 05:24:16 PM Running minimization in window a006...
2020-09-25 05:24:33 PM Running minimization in window a007...
2020-09-25 05:24:52 PM Running minimization in window a008...
2020-09-25 05:25:12 PM Running minimization in window a009...
2020-09-25 05:25:32 PM Running minimization in window a010...
2020-09-25 05:25:50 PM Running minimization in window a011...
2020-09-25 05:26:10 PM Running minimization in window a012...
2020-09-25 05:26:29 PM Running minimization in window a013...
2020-09-25 05:26:49 PM Running minimization in window p000...
2020-09-25 05:27:08 PM Running minimization in window p001...
2020-09-25 05:27:28 PM Running minimization in window p002...
2020-09-25 05:27:46 PM Running minimization in window p003...
2020-09-25 05:28:05 PM Running minimization in window p004...
2020-09-25 05:28:24 PM Running minimization in window p005...
2020-09-25 05:28:46 PM Running minimization in window p006...
2020-09-25 05:29:04 PM Running minimization in window p007...
2020-09-25 05:29:23 PM Running minimization in window p008...
2020-09-25 05:29:41 PM Running minimization in window p009...
2020-09-25 05:30:01 PM Running minimization in window p010...
2020-09-25 05:30:20 PM Running minimization in window p011...
2020-09-25 05:30:39 PM Running minimization in window p012...
2020-09-25 05:30:58 PM Running minimization in window p013...
2020-09-25 05:31:18 PM Running minimization in window p014...
2020-09-25 05:31:37 PM Running minimization in window p015...
2020-09-25 05:31:56 PM Running minimization in window p016...
2020-09-25 05:32:15 PM Running minimization in window p017...

NVT - Equilibration

Here, we will perform 50,000 steps of equilibration with the NVT ensemble.

for window in window_list:
    simulation = GROMACS()
    simulation.executable = "gmx_mpi"
    simulation.n_threads = 8
    simulation.gpu_devices = 0
#     simulation.custom_mdrun_command = "-ntomp 8 -gpu_id 6"
    simulation.path = os.path.join(work_dir, window)
    simulation.maxwarn = 5

    simulation.prefix = "equilibration"
    simulation.topology = f"{base_name}-sol.top"
    simulation.coordinates = "minimize.gro"
    simulation.plumed_file = "plumed.dat"

    simulation.config_pbc_md(ensemble='nvt')
    simulation.control["nsteps"] = 50000

    logging.info(f"Running equilibration in window {window}...")
    simulation.run(overwrite=True)

2020-09-25 05:41:42 PM Running equilibration in window a000...
2020-09-25 05:42:26 PM Running equilibration in window a001...
2020-09-25 05:43:07 PM Running equilibration in window a002...
2020-09-25 05:43:49 PM Running equilibration in window a003...
2020-09-25 05:44:33 PM Running equilibration in window a004...
2020-09-25 05:45:17 PM Running equilibration in window a005...
2020-09-25 05:46:01 PM Running equilibration in window a006...
2020-09-25 05:46:36 PM Running equilibration in window a007...
2020-09-25 05:47:13 PM Running equilibration in window a008...
2020-09-25 05:47:54 PM Running equilibration in window a009...
2020-09-25 05:48:36 PM Running equilibration in window a010...
2020-09-25 05:49:12 PM Running equilibration in window a011...
2020-09-25 05:49:56 PM Running equilibration in window a012...
2020-09-25 05:50:31 PM Running equilibration in window a013...
2020-09-25 05:51:11 PM Running equilibration in window p000...
2020-09-25 05:51:53 PM Running equilibration in window p001...
2020-09-25 05:52:36 PM Running equilibration in window p002...
2020-09-25 05:53:20 PM Running equilibration in window p003...
2020-09-25 05:54:02 PM Running equilibration in window p004...
2020-09-25 05:54:44 PM Running equilibration in window p005...
2020-09-25 05:55:26 PM Running equilibration in window p006...
2020-09-25 05:56:09 PM Running equilibration in window p007...
2020-09-25 05:56:44 PM Running equilibration in window p008...
2020-09-25 05:57:27 PM Running equilibration in window p009...
2020-09-25 05:58:05 PM Running equilibration in window p010...
2020-09-25 05:58:48 PM Running equilibration in window p011...
2020-09-25 05:59:31 PM Running equilibration in window p012...
2020-09-25 06:00:14 PM Running equilibration in window p013...
2020-09-25 06:00:59 PM Running equilibration in window p014...
2020-09-25 06:01:41 PM Running equilibration in window p015...
2020-09-25 06:02:20 PM Running equilibration in window p016...
2020-09-25 06:03:02 PM Running equilibration in window p017...

Production Run

We will run the production phase for 50,000 steps. Note: a proper production run would require simulation in the tens of nanoseconds per window.

for window in window_list:
    simulation = GROMACS()
    simulation.executable = "gmx_mpi"
    simulation.n_threads = 8
    simulation.gpu_devices = 0
#     simulation.custom_mdrun_command = "-ntomp 8 -gpu_id 6"
    simulation.path =  os.path.join(work_dir, window)
    simulation.maxwarn = 5

    simulation.prefix = "production"
    simulation.topology = f"{base_name}-sol.top"
    simulation.coordinates = "equilibration.gro"
    simulation.checkpoint = "equilibration.cpt"
    simulation.plumed_file = "plumed.dat"

    simulation.config_pbc_md(ensemble='npt')
    simulation.control["nsteps"] = 50000

    logging.info(f"Running production in window {window}...")
    simulation.run(overwrite=True)

2020-09-25 06:04:07 PM Running production in window a000...
2020-09-25 06:04:49 PM Running production in window a001...
2020-09-25 06:05:33 PM Running production in window a002...
2020-09-25 06:06:17 PM Running production in window a003...
2020-09-25 06:07:02 PM Running production in window a004...
2020-09-25 06:07:45 PM Running production in window a005...
2020-09-25 06:08:30 PM Running production in window a006...
2020-09-25 06:09:14 PM Running production in window a007...
2020-09-25 06:09:56 PM Running production in window a008...
2020-09-25 06:10:37 PM Running production in window a009...
2020-09-25 06:11:22 PM Running production in window a010...
2020-09-25 06:12:04 PM Running production in window a011...
2020-09-25 06:12:49 PM Running production in window a012...
2020-09-25 06:13:39 PM Running production in window a013...
2020-09-25 06:14:23 PM Running production in window p000...
2020-09-25 06:15:07 PM Running production in window p001...
2020-09-25 06:15:50 PM Running production in window p002...
2020-09-25 06:16:35 PM Running production in window p003...
2020-09-25 06:17:20 PM Running production in window p004...
2020-09-25 06:18:03 PM Running production in window p005...
2020-09-25 06:18:55 PM Running production in window p006...
2020-09-25 06:19:40 PM Running production in window p007...
2020-09-25 06:20:24 PM Running production in window p008...
2020-09-25 06:21:13 PM Running production in window p009...
2020-09-25 06:21:57 PM Running production in window p010...
2020-09-25 06:22:43 PM Running production in window p011...
2020-09-25 06:23:28 PM Running production in window p012...
2020-09-25 06:24:16 PM Running production in window p013...
2020-09-25 06:24:58 PM Running production in window p014...
2020-09-25 06:25:43 PM Running production in window p015...
2020-09-25 06:26:24 PM Running production in window p016...
2020-09-25 06:27:06 PM Running production in window p017...

Analysis

Once the simulation is completed, we can using the Analysis module to determine the binding free energy. We supply the location of the parameter information, a string or list for the file names (wildcards supported), the location of the windows, and the restraints on the guest.

In this example, we use the method ti-block which determines the free energy using thermodynamic iintegration and then estimates the standard error of the mean at each data point using blocking analysis. Bootstrapping it used to determine the uncertainty of the full thermodynamic integral for each phase.

After running compute_free_energy(), a dictionary called results will be populated, that contains the free energy and SEM for each phase of the simulation.

import paprika.analysis as analysis

free_energy = analysis.fe_calc()
free_energy.topology = "cb6-but-dum-sol.pdb"
free_energy.trajectory = 'production*.trr'
free_energy.path = work_dir
free_energy.restraint_list = guest_restraints
free_energy.collect_data()
free_energy.methods = ['ti-block']
free_energy.ti_matrix = "full"
free_energy.boot_cycles = 1000
free_energy.compute_free_energy()

We also need to calculate the free-energy cost of releasing the restraints on the guest molecule.

free_energy.compute_ref_state_work([
    guest_restraints[0], guest_restraints[1], None, None,
    guest_restraints[2], None
])

Then we add the free-energies together and combine the uncertainties to get the binding-free energy

binding_affinity = -1 * (
    free_energy.results["attach"]["ti-block"]["fe"] + \
    free_energy.results["pull"]["ti-block"]["fe"] + \
    free_energy.results["ref_state_work"]
)

sem = np.sqrt(
    free_energy.results["attach"]["ti-block"]["sem"]**2 + \
    free_energy.results["pull"]["ti-block"]["sem"]**2
)

print(f"The binding affinity of butane to cucurbit[6]uril = {binding_affinity.magnitude:0.2f} +/- {sem.magnitude:0.2f} kcal/mol")

The binding affinity of butane to cucurbit[6]uril = -7.55 +/- 1.58 kcal/mol