pAPRika tutorial 5 - APR/Amber with Plumed restraints¶
In this tutorial, we will perform APR calculations for the butane (BUT)–cucurbit[6]uril (CB6) host-guest system. This is a repeat of Tutorial 1 using Plumed-based restraints and the AMBER MD engine. Plumed is a plugin for MD codes that can analyze trajectories and perform free-energy calculations on collective variables. It is a versatile plugin that can interface with a number of MD engines. Here, we will go through the process of converting APR
restraints constructed with pAPRika to a Plumed file and run a short calculation with sander.
Initialize¶
Before you start¶
We will run the simulations in this tutorial using sander (Ambertools) and Plumed. Both of these should be installed in your conda environment if you installed pAPRika though the conda route. However, for Plumed to work with sander 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.
[1]:
import os
'PLUMED_KERNEL' in os.environ.keys()
[1]:
False
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 Mac 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.
[2]:
os.environ['PLUMED_KERNEL'] = f"{os.environ['CONDA_PREFIX']}/lib/libplumedKernel.so"
For running Plumed with Amber outside of this notebook it might be better to export the PLUMED_KERNEL variable into your .bashrc file.
Note: we can run Plumed with AMBER versions 18 and 20, but version 18 requires you to patch the source code first and recompile. Older versions of Amber are not supported. See the Plumed documentation for more details https://www.plumed.org/doc-v2.6/user-doc/html/index.html.
🔵 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.
Define names¶
We will store the files created in this tutorial in a folder called plumed so we don’t mix files with the previous tutorial.
[3]:
base_name = "cb6-but-dum"
work_dir = "plumed"
complex_dir = "complex"
Configure APR Restraints¶
Define anchor atoms¶
See tutorial 1 for the choice of selection
[4]:
# 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.
[5]:
import numpy as np
[6]:
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¶
[7]:
import parmed as pmd
Load complex structure
[8]:
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
[9]:
import paprika.restraints as restraints
[10]:
static_restraints = []
[11]:
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)
[12]:
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)
[13]:
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)
[14]:
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)
[15]:
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)
[16]:
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
[17]:
guest_restraints = []
[18]:
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)
[19]:
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)
[20]:
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.
[21]:
from paprika.restraints.restraints import create_window_list
[22]:
window_list = create_window_list(guest_restraints)
[23]:
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) but we will not do that here. Instead we
will use the built-in position restraints feature in Amber (see Simulation section below).
Note: be careful when specifiying 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.
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).
[24]:
from paprika.restraints.plumed import Plumed
[25]:
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.
[26]:
import shutil
[27]:
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.
Simulation¶
Since we are going to run an implicit solvent simulation, we have everything ready to go. pAPRika has an AMBER 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_gb_min() to setup reasonable default simulation parameters for a minimization in the Generalized-Born ensemble. After that, we directly modify the simulation cntrl section to apply the
positional restraints on the dummy atoms.
We will run the simulations with sander but it is also possible and faster to run this with pmemd or pmemd.cuda if you have them installed.
Note: The difference here compared to Tutorial 1 is that instead of specifying a simulation.restraint_file we will specify simulation.plumed_file.
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.
[28]:
from paprika.simulate import AMBER
Initialize logger
[29]:
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¶
Run a quick minimization in every window. Note that we need to specify simulation.cntrl["ntr"] = 1 to enable the positional restraints on the dummy atoms.
[37]:
for window in window_list:
simulation = AMBER()
simulation.executable = "sander"
simulation.path = f"{work_dir}/{window}/"
simulation.prefix = "minimize"
simulation.topology = "cb6-but-dum.prmtop"
simulation.coordinates = "cb6-but-dum.rst7"
simulation.ref = "cb6-but-dum.rst7"
simulation.plumed_file = "plumed.dat"
simulation.config_gb_min()
simulation.cntrl["ntr"] = 1
simulation.cntrl["restraint_wt"] = 50.0
simulation.cntrl["restraintmask"] = "'@DUM'"
logger.info(f"Running minimization in window {window}...")
simulation.run(overwrite=True)
2020-10-01 11:18:20 AM Running minimization in window a000...
2020-10-01 11:18:31 AM Running minimization in window a001...
2020-10-01 11:18:42 AM Running minimization in window a002...
2020-10-01 11:18:54 AM Running minimization in window a003...
2020-10-01 11:19:05 AM Running minimization in window a004...
2020-10-01 11:19:17 AM Running minimization in window a005...
2020-10-01 11:19:28 AM Running minimization in window a006...
2020-10-01 11:19:42 AM Running minimization in window a007...
2020-10-01 11:19:54 AM Running minimization in window a008...
2020-10-01 11:20:06 AM Running minimization in window a009...
2020-10-01 11:20:18 AM Running minimization in window a010...
2020-10-01 11:20:30 AM Running minimization in window a011...
2020-10-01 11:20:42 AM Running minimization in window a012...
2020-10-01 11:20:53 AM Running minimization in window a013...
2020-10-01 11:21:04 AM Running minimization in window p000...
2020-10-01 11:21:14 AM Running minimization in window p001...
2020-10-01 11:21:26 AM Running minimization in window p002...
2020-10-01 11:21:38 AM Running minimization in window p003...
2020-10-01 11:21:50 AM Running minimization in window p004...
2020-10-01 11:22:01 AM Running minimization in window p005...
2020-10-01 11:22:15 AM Running minimization in window p006...
2020-10-01 11:22:28 AM Running minimization in window p007...
2020-10-01 11:22:40 AM Running minimization in window p008...
2020-10-01 11:22:54 AM Running minimization in window p009...
2020-10-01 11:23:06 AM Running minimization in window p010...
2020-10-01 11:23:19 AM Running minimization in window p011...
2020-10-01 11:23:31 AM Running minimization in window p012...
2020-10-01 11:23:44 AM Running minimization in window p013...
2020-10-01 11:23:55 AM Running minimization in window p014...
2020-10-01 11:24:07 AM Running minimization in window p015...
2020-10-01 11:24:19 AM Running minimization in window p016...
2020-10-01 11:24:30 AM Running minimization in window p017...
Production Run¶
Here we will skip the equilibration step and go straight to production!
[31]:
for window in window_list:
simulation = AMBER()
simulation.executable = "sander"
simulation.path = f"{work_dir}/{window}/"
simulation.prefix = "production"
simulation.topology = "cb6-but-dum.prmtop"
simulation.coordinates = "minimize.rst7"
simulation.ref = "cb6-but-dum.rst7"
simulation.plumed_file = "plumed.dat"
simulation.config_gb_md()
simulation.cntrl["ntr"] = 1
simulation.cntrl["restraint_wt"] = 50.0
simulation.cntrl["restraintmask"] = "'@DUM'"
logger.info(f"Running production in window {window}...")
simulation.run(overwrite=True)
2020-10-01 11:12:15 AM Running production in window a000...
2020-10-01 11:12:25 AM Running production in window a001...
2020-10-01 11:12:34 AM Running production in window a002...
2020-10-01 11:12:44 AM Running production in window a003...
2020-10-01 11:12:55 AM Running production in window a004...
2020-10-01 11:13:06 AM Running production in window a005...
2020-10-01 11:13:19 AM Running production in window a006...
2020-10-01 11:13:29 AM Running production in window a007...
2020-10-01 11:13:39 AM Running production in window a008...
2020-10-01 11:13:49 AM Running production in window a009...
2020-10-01 11:13:59 AM Running production in window a010...
2020-10-01 11:14:10 AM Running production in window a011...
2020-10-01 11:14:19 AM Running production in window a012...
2020-10-01 11:14:29 AM Running production in window a013...
2020-10-01 11:14:39 AM Running production in window p000...
2020-10-01 11:14:49 AM Running production in window p001...
2020-10-01 11:14:59 AM Running production in window p002...
2020-10-01 11:15:10 AM Running production in window p003...
2020-10-01 11:15:21 AM Running production in window p004...
2020-10-01 11:15:33 AM Running production in window p005...
2020-10-01 11:15:43 AM Running production in window p006...
2020-10-01 11:15:55 AM Running production in window p007...
2020-10-01 11:16:05 AM Running production in window p008...
2020-10-01 11:16:18 AM Running production in window p009...
2020-10-01 11:16:31 AM Running production in window p010...
2020-10-01 11:16:45 AM Running production in window p011...
2020-10-01 11:16:56 AM Running production in window p012...
2020-10-01 11:17:05 AM Running production in window p013...
2020-10-01 11:17:16 AM Running production in window p014...
2020-10-01 11:17:26 AM Running production in window p015...
2020-10-01 11:17:37 AM Running production in window p016...
2020-10-01 11:17:49 AM 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.
[32]:
import paprika.analysis as analysis
[33]:
free_energy = analysis.fe_calc()
free_energy.topology = "cb6-but-dum.prmtop"
free_energy.trajectory = 'production*.nc'
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.
[34]:
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
[35]:
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
)
[36]:
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.24 +/- 6.52 kcal/mol