Example 4: Nonequilibrium GCMC Moves

This is a simple example showing how grand can be used to enhance GCMC moves using NCMC. Rather than inserting or deleting a water and then immediately evaluating the energy change, here, waters are gradually inserted or deleted according to a nonequilibrium protocol, based on NCMC theory. The majority of the script below is taken from the previous BPTI example, with a few changes made to demonstrate the differences between conventional GCMC sampling and that enhanced by NCMC. These differences are summarised below:

• The StandardGCMCSphereSampler object is replaced with NonequilibriumGCMCSphereSampler (there is also an equivalent NonequilibriumGCMCSystemSampler object.

• The Sampler object now takes three additional arguments: the integrator argument is needed for the propagation steps, nPertSteps is the number of perturbation steps over which the insertion/deletion is carried out, and nPropStepsPerPert is the number of propagation/relaxation steps between perturbations (the total number of relaxation steps used over the move will be (nPertSteps + 1) * nPropStepsPerPert)

• The Simulation object is created using the gcmc_mover.compound_integrator object.

• Given the improved acceptance rate (and increased cost), fewer GCMC moves are required when using NCMC (note that the full benefit may not be seen for this BPTI system, as once the three waters are inserted, then the likelihood of further acceptances becomes very low).

The documentation for these functions can be found in the “grand package” section. The full script is included below.

"""
Description
-----------
Example script of how to run GCMC/MD in OpenMM for a BPTI system, showing
how the GCMC moves can be enhanced using nonequilibrium protocols (NCMC)

Note that this simulation is only an example, and is not long enough
to see equilibrated behaviour

Marley Samways
"""

from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
from sys import stdout

from openmmtools.integrators import BAOABIntegrator

import grand

# Load in PDB file
pdb = PDBFile('bpti-equil.pdb')

# Add ghost water molecules, which can be inserted
pdb.topology, pdb.positions, ghosts = grand.utils.add_ghosts(pdb.topology,
                                                             pdb.positions,
                                                             n=5,
                                                             pdb='bpti-ghosts.pdb')

# Create system
ff = ForceField('amber14-all.xml', 'amber14/tip3p.xml')
system = ff.createSystem(pdb.topology,
                         nonbondedMethod=PME,
                         nonbondedCutoff=12.0*angstroms,
                         switchDistance=10.0*angstroms,
                         constraints=HBonds)

# Define atoms around which the GCMC sphere is based
ref_atoms = [{'name': 'CA', 'resname': 'TYR', 'resid': '10'},
             {'name': 'CA', 'resname': 'ASN', 'resid': '43'}]

# BAOAB Langevin integrator
integrator = BAOABIntegrator(300*kelvin, 1.0/picosecond, 0.002*picoseconds)

# Define the NCMC Sampler
gcncmc_mover = grand.samplers.NonequilibriumGCMCSphereSampler(system=system,
                                                              topology=pdb.topology,
                                                              temperature=300*kelvin,
                                                              integrator=integrator,
                                                              # Make this a 10 ps protocol
                                                              nPertSteps=99, nPropStepsPerPert=50,
                                                              referenceAtoms=ref_atoms,
                                                              sphereRadius=4.2*angstroms,
                                                              log='bpti-gcmc.log',
                                                              dcd='bpti-raw.dcd',
                                                              rst='bpti-rst.rst7',
                                                              overwrite=False)

platform = Platform.getPlatformByName('CUDA')
platform.setPropertyDefaultValue('Precision', 'mixed')

simulation = Simulation(pdb.topology, system, gcncmc_mover.compound_integrator, platform)
simulation.context.setPositions(pdb.positions)
simulation.context.setVelocitiesToTemperature(300*kelvin)
simulation.context.setPeriodicBoxVectors(*pdb.topology.getPeriodicBoxVectors())

# Switch off ghost waters and those in sphere (to start fresh)
gcncmc_mover.initialise(simulation.context, ghosts)
gcncmc_mover.deleteWatersInGCMCSphere()

# Equilibrate water distribution - 10k moves over 5 ps
print("Equilibration...")
for i in range(50):
    # Carry out 2 moves every 100 fs
    gcncmc_mover.move(simulation.context, 1)
    simulation.step(50)
print("{}/{} equilibration GCMC moves accepted. N = {}".format(gcncmc_mover.n_accepted,
                                                               gcncmc_mover.n_moves,
                                                               gcncmc_mover.N))

# Add StateDataReporter for production
simulation.reporters.append(StateDataReporter(stdout,
                                              1000,
                                              step=True,
                                              potentialEnergy=True,
                                              temperature=True,
                                              volume=True))
# Reset GCMC statistics
gcncmc_mover.reset()

# Run simulation - 5k moves over 50 ps
print("\nProduction")
for i in range(50):
    # Carry out 5 GCMC moves per 1 ps of MD
    simulation.step(500)
    gcncmc_mover.move(simulation.context, 2)
    # Write data out
    gcncmc_mover.report(simulation)

#
# Need to process the trajectory for visualisation
#

# Shift ghost waters outside the simulation cell
trj = grand.utils.shift_ghost_waters(ghost_file='gcmc-ghost-wats.txt',
                                     topology='bpti-ghosts.pdb',
                                     trajectory='bpti-raw.dcd')

# Centre the trajectory on a particular residue
trj = grand.utils.recentre_traj(t=trj, resname='TYR', resid=10)

# Align the trajectory to the protein
grand.utils.align_traj(t=trj, output='bpti-gcmc.dcd')

# Write out a PDB trajectory of the GCMC sphere
grand.utils.write_sphere_traj(radius=4.2,
                              ref_atoms=ref_atoms,
                              topology='bpti-ghosts.pdb',
                              trajectory='bpti-gcmc.dcd',
                              output='gcmc_sphere.pdb',
                              initial_frame=True)