Internal coordinate darting

.travis.yml

@@ -64,6 +64,7 @@ before_install:
 
   # Add conda-forge channel to top priority
   - conda config --add channels conda-forge
+  - conda config --add channels sgill2
 
 install:
   # Create test environment for package

README.md

@@ -30,12 +30,17 @@ BLUES is compatible with MacOSX/Linux with Python>=3.5 (blues<=1.1 still works w
 Install [miniconda](http://conda.pydata.org/miniconda.html) according to your system.
 
 ## Requirements
-Starting from v1.2, you will need the OpenEye toolkits and related tools:
+Starting from v1.2, you will need the OpenEye toolkits and related tools to use the `SideChainMove` class:
 ```bash
 conda install -c openeye/label/Orion -c omnia oeommtools packmol
 
 # Requires OpenEye License
 conda install -c openeye openeye-toolkits
+
+# Requirements for the MolDartMove class
+conda install --yes -c sgill2 chemcoord
+conda install -c omnia mdtraj parmed yank openmmtools=0.14.0
+conda install --yes -c conda-forge future
 ```
 
 ## Installation

blues/__init__.py

@@ -4,6 +4,8 @@
 """
 # Add imports here
 from blues import integrators, moves, reporters, settings, simulation, utils
+#from blues import integrators, moves, reporters, settings, simulation, utils, moldart
+
 
 # Handle versioneer
 from ._version import get_versions

blues/examples/example.py

@@ -0,0 +1,109 @@
+"""
+example.py: Provides an example script to run BLUES and
+benchmark the run on a given platform
+
+Authors: Samuel C. Gill
+Contributors: Nathan M. Lim, David L. Mobley
+
+* Benchmarking related code adapted from:
+https://github.com/pandegroup/openmm/blob/master/examples/benchmark.py
+(Authors: Peter Eastman)
+"""
+
+from __future__ import print_function
+from blues.moves import RandomLigandRotationMove
+from blues.engine import MoveEngine
+from blues import utils
+from blues.simulation import Simulation, SimulationFactory
+import parmed
+from simtk import openmm
+from optparse import OptionParser
+import sys
+import logging
+from blues.reporters import init_logger, BLUESHDF5Reporter, BLUESStateDataReporter
+
+def runNCMC(platform_name, nstepsNC, nprop, outfname):
+
+    #Generate the ParmEd Structure
+    prmtop = utils.get_data_filename('blues', 'tests/data/eqToluene.prmtop')
+    inpcrd = utils.get_data_filename('blues', 'tests/data/eqToluene.inpcrd')
+    struct = parmed.load_file(prmtop, xyz=inpcrd)
+
+    #Define some options
+    opt = { 'temperature' : 300.0, 'friction' : 1, 'dt' : 0.002,
+            'nIter' : 100, 'nstepsNC' : 10000, 'nstepsMD' : 10000, 'nprop' : 1,
+            'nonbondedMethod' : 'PME', 'nonbondedCutoff': 10,
+            'constraints': 'HBonds', 'freeze_distance' : 5.0,
+            'trajectory_interval' : 2000, 'reporter_interval' : 1000,
+            'write_move' : False,
+            'platform' : platform_name,
+            'outfname' : 't4-toluene'}
+
+
+    logger = init_logger(logging.getLogger(), level=logging.INFO, outfname=opt['outfname'])
+    opt['Logger'] = logger
+
+    #Define the 'model' object we are perturbing here.
+    # Calculate particle masses of object to be moved
+    ligand = RandomLigandRotationMove(struct, 'LIG')
+
+    # Initialize object that proposes moves.
+    ligand_mover = MoveEngine(ligand)
+
+    # Generate the MD, NCMC, ALCHEMICAL Simulation objects
+    simulations = SimulationFactory(struct, ligand_mover, **opt)
+    simulations.createSimulationSet()
+
+    # Add reporters to MD simulation.
+    traj_reporter = openmm.app.DCDReporter(opt['outfname']+'-nc{}.dcd'.format(nstepsNC), opt['trajectory_interval'])
+    md_progress_reporter = BLUESStateDataReporter(logger, separator="\t", title='md',
+                                reportInterval=opt['reporter_interval'],
+                                step=True, totalSteps=opt['nIter']*opt['nstepsMD'],
+                                time=False, speed=True, progress=True, remainingTime=True)
+    simulations.md.reporters.append(traj_reporter)
+    simulations.md.reporters.append(md_progress_reporter)
+
+    # Add reporters to NCMC simulation.
+    ncmc_reporter = BLUESHDF5Reporter(file=opt['outfname']+'-pmoves.h5',
+                                    reportInterval=1,
+                                    coordinates=True, frame_indices=[1,opt['nstepsNC']],
+                                    time=False, cell=True, temperature=False,
+                                    potentialEnergy=False, kineticEnergy=False,
+                                    velocities=False, atomSubset=None,
+                                    protocolWork=True, alchemicalLambda=True,
+                                    parameters=opt, environment=True)
+    ncmc_progress_reporter = BLUESStateDataReporter(logger, separator="\t", title='ncmc',
+                                reportInterval=opt['reporter_interval'],
+                                step=True, totalSteps=opt['nstepsNC'],
+                                time=False, speed=True, progress=True, remainingTime=True)
+    simulations.nc.reporters.append(ncmc_reporter)
+    simulations.nc.reporters.append(ncmc_progress_reporter)
+
+    # Run BLUES Simulation
+    blues = Simulation(simulations, ligand_mover, **opt)
+    blues.run(opt['nIter'])
+
+parser = OptionParser()
+parser.add_option('-f', '--force', action='store_true', default=False,
+                  help='run BLUES example without GPU platform')
+parser.add_option('-n','--ncmc', dest='nstepsNC', type='int', default=100,
+                  help='number of NCMC steps')
+parser.add_option('-p','--nprop', dest='nprop', type='int', default=1,
+                  help='number of propgation steps')
+parser.add_option('-o','--output', dest='outfname', type='str', default="blues",
+                  help='Filename for output DCD')
+(options, args) = parser.parse_args()
+
+
+
+platformNames = [openmm.Platform.getPlatform(i).getName() for i in range(openmm.Platform.getNumPlatforms())]
+if 'CUDA' in platformNames:
+    runNCMC('CUDA', options.nstepsNC, options.nprop, options.outfname)
+elif 'OpenCL' in platformNames:
+    runNCMC('OpenCL',options.nstepsNC, options.nprop, options.outfname)
+else:
+    if options.force:
+        runNCMC('CPU', options.nstepsNC, options.outfname)
+    else:
+        print('WARNING: Could not find a valid CUDA/OpenCL platform. BLUES is not recommended on CPUs.')
+        print("To run on CPU: 'python blues/example.py -f'")

blues/examples/example_moldart.py

@@ -0,0 +1,109 @@
+"""
+example.py: Provides an example script to run BLUES and
+benchmark the run on a given platform
+
+Authors: Samuel C. Gill
+Contributors: Nathan M. Lim, David L. Mobley
+
+* Benchmarking related code adapted from:
+https://github.com/pandegroup/openmm/blob/master/examples/benchmark.py
+(Authors: Peter Eastman)
+"""
+
+from __future__ import print_function
+from blues.moves import MolDartMove
+from blues.engine import MoveEngine
+from blues import utils
+from blues.simulation import Simulation, SimulationFactory
+import parmed
+from simtk import openmm
+from optparse import OptionParser
+import mdtraj as md
+import logging, sys
+from blues.reporters import init_logger, BLUESHDF5Reporter, BLUESStateDataReporter
+
+
+def runNCMC(platform_name):
+    #Define some options
+    opt = { 'temperature' : 300.0, 'friction' : 1, 'dt' : 0.002,
+            'nIter' : 10, 'nstepsNC' : 2000, 'nstepsMD' : 5000,
+            'nonbondedMethod' : 'PME', 'nonbondedCutoff': 10, 'constraints': 'HBonds',
+            'trajectory_interval' : 1000, 'reporter_interval' : 250,
+            'platform' : 'CUDA',
+            'outfname' : 't4-tol',
+            'write_move' : False,
+            'nprop':5,
+            'freeze_distance' : 5.0,
+            'verbose' : True
+ }
+
+    logger = init_logger(logging.getLogger(), level=logging.INFO, outfname=opt['outfname'])
+    opt['Logger'] = logger
+
+    #Generate the ParmEd Structure
+    prmtop = utils.get_data_filename('blues', 'tests/data/eqToluene.prmtop')#
+    inpcrd = utils.get_data_filename('blues', 'tests/data/eqToluene.inpcrd')
+    struct = parmed.load_file(prmtop, xyz=inpcrd)
+
+    #Define the 'model' object we are perturbing here.
+    # Calculate particle masses of object to be moved
+    posA = utils.get_data_filename('blues', 'tests/data/posA.pdb')
+    posB = utils.get_data_filename('blues', 'tests/data/posB.pdb')
+    traj = md.load(inpcrd, top=prmtop)
+    fit_atoms = traj.top.select("protein")
+
+    ligand = MolDartMove(structure=struct, resname='LIG',
+                                      pdb_files=[posA, posB],
+                                      fit_atoms=fit_atoms,
+                                      restrained_receptor_atoms=[1605, 1735, 1837],
+                                      rigid_ring=True
+                                      )
+
+    # Initialize object that proposes moves.
+    ligand_mover = MoveEngine(ligand)
+
+    # Generate the MD, NCMC, ALCHEMICAL Simulation objects
+    simulations = SimulationFactory(struct, ligand_mover, **opt)
+    simulations.createSimulationSet()
+
+    # Add reporters to MD simulation.
+    #traj_reporter = openmm.app.DCDReporter(opt['outfname']+'-nc{}.dcd'.format(opt['nstepsNC']), 1)
+    from blues.reporters import NetCDF4Reporter
+    traj_reporter = NetCDF4Reporter(opt['outfname']+'-nc{}.nc'.format(opt['nstepsNC']), 1000)
+    md_progress_reporter = openmm.app.StateDataReporter(sys.stdout, separator="\t",
+                                reportInterval=opt['reporter_interval'],
+                                step=True, totalSteps=opt['nIter']*opt['nstepsMD'],
+                                time=True, speed=True, progress=True,
+                                elapsedTime=True, remainingTime=True)
+    simulations.md.reporters.append(traj_reporter)
+    simulations.md.reporters.append(md_progress_reporter)
+
+    # Add reporters to NCMC simulation.
+    ncmc_progress_reporter = openmm.app.StateDataReporter(sys.stdout, separator="\t",
+                                reportInterval=opt['reporter_interval'],
+                                step=True, totalSteps=opt['nstepsNC'],
+                                time=True, speed=True, progress=True,
+                                elapsedTime=True, remainingTime=True)
+    #simulations.nc.reporters.append(traj_reporter)
+    simulations.nc.reporters.append(ncmc_progress_reporter)
+
+    blues = Simulation(simulations, ligand_mover, **opt)
+    blues.run(opt['nIter'])
+
+parser = OptionParser()
+parser.add_option('-f', '--force', action='store_true', default=False,
+                  help='run BLUES example without GPU platform')
+(options, args) = parser.parse_args()
+
+platformNames = [openmm.Platform.getPlatform(i).getName() for i in range(openmm.Platform.getNumPlatforms())]
+
+if 'OpenCL' in platformNames:
+    runNCMC('OpenCL')
+elif 'CUDA' in platformNames:
+    runNCMC('CUDA')
+else:
+    if options.force:
+        runNCMC('CPU')
+    else:
+        print('WARNING: Could not find a valid CUDA/OpenCL platform. BLUES is not recommended on CPUs.')
+        print("To run on CPU: 'python blues/example.py -f'")

blues/examples/example_sidechain.py

@@ -0,0 +1,35 @@
+from blues.moves import SideChainMove
+from blues.engine import MoveEngine
+from blues.simulation import *
+from blues.settings import *
+
+# Parse a YAML configuration, return as Dict
+cfg = Settings('sidechain_cuda.yaml').asDict()
+structure = cfg['Structure']
+
+#Select move type
+sidechain = SideChainMove(structure, [1])
+#Iniitialize object that selects movestep
+sidechain_mover = MoveEngine(sidechain)
+
+#Generate the openmm.Systems outside SimulationFactory to allow modifications
+systems = SystemFactory(structure, sidechain.atom_indices, cfg['system'])
+
+#Generate the OpenMM Simulations
+simulations = SimulationFactory(systems, sidechain_mover, cfg['simulation'], cfg['md_reporters'], cfg['ncmc_reporters'])
+
+# Run BLUES Simulation
+blues = BLUESSimulation(simulations, cfg['simulation'])
+blues.run()
+
+#Analysis
+import mdtraj as md
+import numpy as np
+
+traj = md.load_netcdf('vacDivaline-test/vacDivaline.nc', top='tests/data/vacDivaline.prmtop')
+indicies = np.array([[0, 4, 6, 8]])
+dihedraldata = md.compute_dihedrals(traj, indicies)
+with open("vacDivaline-test/dihedrals.txt", 'w') as output:
+    for value in dihedraldata:
+        output.write("%s\n" % str(value)[1:-1])
+

blues/integrators.py

@@ -1,6 +1,6 @@
 import simtk
 from openmmtools.integrators import AlchemicalNonequilibriumLangevinIntegrator
-
+import inspect
 # Energy unit used by OpenMM unit system
 _OPENMM_ENERGY_UNIT = simtk.unit.kilojoules_per_mole
 
@@ -122,6 +122,10 @@ def __init__(self,
             nsteps_neq=nsteps_neq)
 
         self._prop_lambda = self._get_prop_lambda(prop_lambda)
+        frame = inspect.currentframe()
+        args, _, _, values = inspect.getargvalues(frame)
+        inputs = dict([(i, values[i]) for i in args if i is not 'self'])
+        self.int_kwargs = inputs
 
         # add some global variables relevant to the integrator
         kB = simtk.unit.BOLTZMANN_CONSTANT_kB * simtk.unit.AVOGADRO_CONSTANT_NA
@@ -246,4 +250,4 @@ def reset(self):
         self.setGlobalVariableByName("perturbed_pe", 0.0)
         self.setGlobalVariableByName("unperturbed_pe", 0.0)
         self.setGlobalVariableByName("prop", 1)
-        super(AlchemicalExternalLangevinIntegrator, self).reset()
+        super(AlchemicalExternalLangevinIntegrator, self).reset()

blues/moldart/__init__.py

@@ -0,0 +1 @@
+from blues.moldart import chemcoord, darts, lin_math, move