pAPRika tutorial 3 - K-Cl dissociation¶

In this example, we will setup and simulate KCl dissociation using a single distance restraints. This tutorial will demonstrate the use of pAPRika besides host-guest systems.

[17]:

import os
import json

import numpy as np
import parmed as pmd

from paprika.analysis import fe_calc
from paprika.build.system import TLeap
from paprika.io import NumpyEncoder
from paprika.restraints.restraints import DAT_restraint
from paprika.restraints.amber import amber_restraint_linefrom paprika.restraints.utils import create_window_listfrom paprika.simulate import AMBER

Specify directory for data¶

In this case, we just need some initial coordinates. In other cases, we might need mol2 or frcmod files.

[13]:

from pathlib import Path
cwd = Path().resolve()
k_cl_pdb = os.path.abspath(os.path.join(cwd, "../../paprika/data/k-cl/k-cl.pdb"))

Setup the calculation¶

Build the vacuum `prmtop` and `inpcrd` files¶

[18]:

# Build the model in vacuum

system = TLeap()
system.template_lines = [
    "source leaprc.water.tip3p",
    "loadamberparams frcmod.ionsjc_tip3p",
    f"model = loadpdb {k_cl_pdb}",
]
system.output_path = "tmp"
system.output_prefix = "k-cl"
system.pbc_type = None
system.target_waters = None
system.neutralize = False
system.build()

Specify the number of windows for the umbrella sampling¶

These are overkill; I have been testing how much data we need to converge this calculation and how quickly we can run a stripped-down version on Travis.

[19]:

attach_fractions = np.linspace(0, 1.0, 25)
initial_distance = 2.65
pull_distances = np.linspace(0 + initial_distance, 16.0 + initial_distance, 40)

Add a single distance restraint between K+ and Cl-¶

[20]:

restraint = DAT_restraint()
restraint.continuous_apr = True
restraint.amber_index = True
restraint.topology = k_cl_pdb
restraint.mask1 = "@K+"
restraint.mask2 = "@Cl-"

restraint.attach["target"] = initial_distance
restraint.attach["fraction_list"] = attach_fractions
restraint.attach["fc_final"] = 10.0
restraint.pull["fc"] = restraint.attach["fc_final"]
restraint.pull["target_list"] = pull_distances
restraint.initialize()

Optionally, add a “wall restraint” to define the bound state and speed convergence¶

This will apply a harmonic potential that keeps the “guest” Cl- within 3.5 Angstroms.

[21]:

wall = DAT_restraint()
wall.auto_apr = False
wall.amber_index = True
wall.topology = k_cl_pdb
wall.mask1 = "@K+"
wall.mask2 = "@Cl-"

wall.attach["fc_initial"] = 1.0
wall.attach["fc_final"] = 1.0

wall.custom_restraint_values["rk2"] = 1.0
wall.custom_restraint_values["rk3"] = 1.0
wall.custom_restraint_values["r1"] = 0.0
wall.custom_restraint_values["r2"] = 3.5
wall.custom_restraint_values["r3"] = 3.5
wall.custom_restraint_values["r4"] = 999

wall.attach["target"] = 3.5
wall.attach["num_windows"] = len(attach_fractions)

wall.initialize()

Create the directories for each window and write the AMBER-style restraint input file¶

This makes it easy to run each window in parallel as a separate simulation.

[22]:

# Create folder for the working directory
window_list = create_window_list([restraint])
for window in window_list:
    os.makedirs(f"tmp/windows/{window}", exist_ok=True)

# Write the AMBER restraint files for each window
for window in window_list:
    with open(f"tmp/windows/{window}/disang.rest", "a") as file:
        if window[0] == "a":
            for r in [restraint, wall]:
                string = amber_restraint_line(r, window)
                if string is not None:
                    file.write(string)
        else:
            string = amber_restraint_line(restraint, window)
            file.write(string)

# Generate the topology and coordinates file for each window
for window in window_list:
    if window[0] == "a":
        structure = pmd.load_file("tmp/k-cl.prmtop", "tmp/k-cl.rst7", structure=True)
        for atom in structure.atoms:
            if atom.name == "Cl-":
                atom.xz = 2.65
        structure.save(f"tmp/windows/{window}/k-cl.prmtop", overwrite=True)
        structure.save(f"tmp/windows/{window}/k-cl.rst7", overwrite=True)

    elif window[0] == "p":
        structure = pmd.load_file("tmp/k-cl.prmtop", "tmp/k-cl.rst7", structure=True)
        target_difference = (
            restraint.phase["pull"]["targets"][int(window[1:])]
            - restraint.phase["pull"]["targets"][0]
        )
        print(
            f"In window {window} we will translate the guest {target_difference.magnitude:0.1f}."
        )
        for atom in structure.atoms:
            if atom.name == "Cl-":
                atom.xz += target_difference.magnitude
        structure.save(f"tmp/windows/{window}/k-cl.prmtop", overwrite=True)
        structure.save(f"tmp/windows/{window}/k-cl.rst7", overwrite=True)

In window p000 we will translate the guest 0.0 Angstroms.
In window p001 we will translate the guest 0.4 Angstroms.
In window p002 we will translate the guest 0.8 Angstroms.
In window p003 we will translate the guest 1.2 Angstroms.
In window p004 we will translate the guest 1.6 Angstroms.
In window p005 we will translate the guest 2.1 Angstroms.
In window p006 we will translate the guest 2.5 Angstroms.
In window p007 we will translate the guest 2.9 Angstroms.
In window p008 we will translate the guest 3.3 Angstroms.
In window p009 we will translate the guest 3.7 Angstroms.
In window p010 we will translate the guest 4.1 Angstroms.
In window p011 we will translate the guest 4.5 Angstroms.
In window p012 we will translate the guest 4.9 Angstroms.
In window p013 we will translate the guest 5.3 Angstroms.
In window p014 we will translate the guest 5.7 Angstroms.
In window p015 we will translate the guest 6.2 Angstroms.
In window p016 we will translate the guest 6.6 Angstroms.
In window p017 we will translate the guest 7.0 Angstroms.
In window p018 we will translate the guest 7.4 Angstroms.
In window p019 we will translate the guest 7.8 Angstroms.
In window p020 we will translate the guest 8.2 Angstroms.
In window p021 we will translate the guest 8.6 Angstroms.
In window p022 we will translate the guest 9.0 Angstroms.
In window p023 we will translate the guest 9.4 Angstroms.
In window p024 we will translate the guest 9.8 Angstroms.
In window p025 we will translate the guest 10.3 Angstroms.
In window p026 we will translate the guest 10.7 Angstroms.
In window p027 we will translate the guest 11.1 Angstroms.
In window p028 we will translate the guest 11.5 Angstroms.
In window p029 we will translate the guest 11.9 Angstroms.
In window p030 we will translate the guest 12.3 Angstroms.
In window p031 we will translate the guest 12.7 Angstroms.
In window p032 we will translate the guest 13.1 Angstroms.
In window p033 we will translate the guest 13.5 Angstroms.
In window p034 we will translate the guest 13.9 Angstroms.
In window p035 we will translate the guest 14.4 Angstroms.
In window p036 we will translate the guest 14.8 Angstroms.
In window p037 we will translate the guest 15.2 Angstroms.
In window p038 we will translate the guest 15.6 Angstroms.
In window p039 we will translate the guest 16.0 Angstroms.

Optionally, tweak some parameters, like changing the charge of K+ to 1.3 and Cl- to -1.3¶

[23]:

for window in window_list:
    structure = pmd.load_file(
        f"tmp/windows/{window}/k-cl.prmtop",
        f"tmp/windows/{window}/k-cl.rst7",
        structure=True,
    )
    for atom in structure.atoms:
        if atom.name == "Cl-":
            atom.charge = -1.3
        elif atom.name == "K+":
            atom.charge = 1.3
    structure.save(f"tmp/windows/{window}/k-cl.prmtop", overwrite=True)
    structure.save(f"tmp/windows/{window}/k-cl.rst7", overwrite=True)

Solvate the structure in each window to the same number of waters¶

[24]:

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

    if os.path.exists(f"tmp/windows/{window}/k-cl-sol.prmtop"):
        print("Skipping...")
        continue

    structure = pmd.load_file(
        f"tmp/windows/{window}/k-cl.prmtop", f"tmp/windows/{window}/k-cl.rst7"
    )

    if not os.path.exists(f"tmp/windows/{window}/k-cl.pdb"):
        structure.save(f"tmp/windows/{window}/k-cl.pdb")

    system = TLeap()
    system.output_path = os.path.join("tmp", "windows", window)
    system.output_prefix = "k-cl-sol"

    system.target_waters = 2000
    system.neutralize = False
    system.template_lines = [
        "source leaprc.water.tip3p",
        "model = loadpdb k-cl.pdb"
    ]
    system.build()

Solvating window a000...
Solvating window a001...
Solvating window a002...
Solvating window a003...
Solvating window a004...
Solvating window a005...
Solvating window a006...
Solvating window a007...
Solvating window a008...
Solvating window a009...
Solvating window a010...
Solvating window a011...
Solvating window a012...
Solvating window a013...
Solvating window a014...
Solvating window a015...
Solvating window a016...
Solvating window a017...
Solvating window a018...
Solvating window a019...
Solvating window a020...
Solvating window a021...
Solvating window a022...
Solvating window a023...
Solvating window p000...
Solvating window p001...
Solvating window p002...
Solvating window p003...
Solvating window p004...
Solvating window p005...
Solvating window p006...
Solvating window p007...
Solvating window p008...
Solvating window p009...
Solvating window p010...
Solvating window p011...
Solvating window p012...
Solvating window p013...
Solvating window p014...
Solvating window p015...
Solvating window p016...
Solvating window p017...
Solvating window p018...
Solvating window p019...
Solvating window p020...
Solvating window p021...
Solvating window p022...
Solvating window p023...
Solvating window p024...
Solvating window p025...
Solvating window p026...
Solvating window p027...
Solvating window p028...
Solvating window p029...
Solvating window p030...
Solvating window p031...
Solvating window p032...
Solvating window p033...
Solvating window p034...
Solvating window p035...
Solvating window p036...
Solvating window p037...
Solvating window p038...
Solvating window p039...

Run the calculation¶

We have a few helper functions – like config_pbc_min() and config_pbc_md() – that help setup some smart defaults for AMBER. (I’ll make a note to work on adding this for the OpenMM side of things.) The simulations can either be run directly, as indicated below, with simulation.run() or the input file can be written using _amber_write_input_file() and then wrapped using a cluster script (like PBS or whatever).

Energy Minimization¶

[ ]:

for window in window_list:
    simulation = AMBER()
    simulation.executable = "pmemd.cuda"

    simulation.path = f"tmp/windows/{window}/"
    simulation.prefix = "minimize"

    simulation.topology = "k-cl-sol.prmtop"
    simulation.coordinates = "k-cl-sol.rst7"
    simulation.ref = "k-cl-sol.rst7"
    simulation.restraint_file = "disang.rest"

    simulation.config_pbc_min()
    simulation.cntrl["ntr"] = 1
    simulation.cntrl["restraint_wt"] = 50.0
    simulation.cntrl["restraintmask"] = "':1-2'"
    print(f"Running minimization in window {window}...")
    simulation.run()

Production Run¶

[ ]:

# Simulate
for window in window_list:
    simulation = AMBER()
    simulation.executable = "pmemd.cuda"

    simulation.path = f"tmp/windows/{window}/"
    simulation.prefix = "production"

    simulation.topology = "k-cl-sol.prmtop"
    simulation.coordinate = "minimize.rst7"
    simulation.ref = "k-cl-sol.rst7"
    simulation.restraint_file = "disang.rest"

    simulation.config_pbc_md()
    simulation.cntrl["nstlim"] = 50000

    print(f"Running production in window {window}...")
    simulation.run()

Analyze the simulation¶

Setup the analysis¶

The analysis needs to know about:

The parameter file that was used for the molecules,
The simulation path,
The trajectories, and
The method to do the analysis (e.g., TI for the free energy with blocking analysis for the SEM)

[ ]:

free_energy = fe_calc()
free_energy.topology = "k-cl-sol.prmtop"
free_energy.trajectory = "production*.nc"
free_energy.path = "tmp/windows"
free_energy.restraint_list = [restraint]
free_energy.collect_data()
free_energy.methods = ["ti-block", "mbar-block"]
free_energy.ti_matrix = "full"
free_energy.boot_cycles = 1000

Run the analysis and save the results¶

[ ]:

free_energy.compute_free_energy()
free_energy.compute_ref_state_work([restraint, None, None, None, None, None])

[ ]:

with open("./tmp/results.json", "w") as f:
    dumped = json.dumps(free_energy.results, cls=NumpyEncoder)
    f.write(dumped)

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 for K+ (+1.3) and Cl- (-1.3) = {binding_affinity.magnitude:0.2f} +/- {sem.magnitude:0.2f} kcal/mol"
)