[WIP] Modifications to acceptance ratio from Moves

blues/engine.py

@@ -74,4 +74,8 @@ def runEngine(self, context):
             print(e)
             raise SystemExit
 
-        return new_context
+        return new_context
+
+    def resetAcceptanceRatios(self):
+        for move in self.moves:
+            move.acceptance_ratio = 1.0

blues/moves.py

@@ -6,18 +6,16 @@
 a combination of other pre-defined moves such as via instances of Move.
 
 Authors: Samuel C. Gill
-Contributors: Nathan M. Lim, Kalistyn Burley, David L. Mobley 
+Contributors: Nathan M. Lim, Kalistyn Burley, David L. Mobley
 """
 
 import parmed
 from simtk import unit
 import mdtraj
 import numpy as np
-import sys, traceback
 import math
 import copy
 import random
-import os
 from openeye.oechem import *
 
 
@@ -36,8 +34,9 @@ class Move(object):
 
     def __init__(self):
         """Initialize the Move object
-        Currently empy.
         """
+        self.acceptance_ratio = 1.0
+
 
     def initializeSystem(self, system, integrator):
         """If the system or integrator needs to be modified to perform the move
@@ -79,7 +78,7 @@ def beforeMove(self, context):
 
         """
         return context
-        
+
     def afterMove(self, context):
         """This method is called at the end of the NCMC portion if the
         context needs to be checked or modified before performing the move
@@ -98,7 +97,7 @@ def afterMove(self, context):
         """
 
         return context
-        
+
     def _error(self, context):
         """This method is called if running during NCMC portion results
         in an error. This allows portions of the context, such as the
@@ -147,6 +146,8 @@ class RandomLigandRotationMove(Move):
     """
 
     def __init__(self, structure, resname='LIG'):
+        super(RandomLigandRotationMove, self).__init__()
+
         """Initialize the model.
         Parameters
         ----------
@@ -284,6 +285,7 @@ class SideChainMove(Move):
         and applies a rotation matrix to the target atoms to update their coordinates"""
 
     def __init__(self, structure, residue_list):
+        super(SideChainMove, self).__init__()
         self.structure = structure
         self.molecule = self._pmdStructureToOEMol()
         self.residue_list = residue_list

blues/simulation.py

@@ -8,9 +8,7 @@
 from simtk import unit, openmm
 from simtk.openmm import app
 import parmed, math
-import mdtraj
-import sys, time
-from datetime import datetime
+import time
 from openmmtools import alchemy
 from blues.integrators import AlchemicalExternalLangevinIntegrator
 import logging
@@ -257,12 +255,12 @@ def __init__(self, simulations, move_engine, **opt):
 
         self.accept = 0
         self.reject = 0
-        self.accept_ratio = 0
+        self.accept_probability = 0
 
         #if nstepsNC not specified, set it to 0
         #will be caught if NCMC simulation is run
         if (self.opt['nstepsNC'] % 2) != 0:
-            raise ValueError('nstepsNC needs to be even to ensure the protocol is symmetric (currently %i)' % (nstepsNC))
+            raise ValueError('nstepsNC needs to be even to ensure the protocol is symmetric (currently %i)' % (self.opt['nstepsNC']))
         else:
             self.movestep = int(self.opt['nstepsNC']) / 2
 
@@ -455,19 +453,34 @@ def acceptRejectNCMC(self, temperature=300, write_move=False, **opt):
             correction_factor = (nc_state0['potential_energy'] - md_state0['potential_energy'] + alch_state1['potential_energy'] - nc_state1['potential_energy']) * (-1.0/self.nc_integrator.kT)
             log_ncmc = log_ncmc + correction_factor
 
-        if log_ncmc > randnum:
-            self.accept += 1
-            self.log.info('NCMC MOVE ACCEPTED: log_ncmc {} > randnum {}'.format(log_ncmc, randnum) )
-            self.md_sim.context.setPositions(nc_state1['positions'])
-            if write_move:
-            	self.writeFrame(self.md_sim, '{}acc-it{}.pdb'.format(self.opt['outfname'],self.current_iter))
-
-        else:
+        # If the move has an acceptance ratio factor, apply that to the overall acceptance
+        # but first check if the move acceptance ratio is 0 since that will always be rejected
+        if self.move_engine.moves[self.move_engine.selected_move].acceptance_ratio == 0:
             self.reject += 1
-            self.log.info('NCMC MOVE REJECTED: log_ncmc {} < {}'.format(log_ncmc, randnum) )
+            self.log.info('NCMC MOVE REJECTED: move acceptance ratio is 0')
             self.nc_context.setPositions(md_state0['positions'])
+        else:
+            # Since we're working with the log of the energy difference we can add the log of the move acceptance ratio
+            # ln(xy) = ln(x) + ln(y), where y is acceptance_ratio and x is the protocol work
+            log_ncmc = log_ncmc + math.log(self.move_engine.moves[self.move_engine.selected_move].acceptance_ratio)
+
+            if log_ncmc > randnum:
+                self.accept += 1
+                self.log.info('NCMC MOVE ACCEPTED: log_ncmc {} > randnum {}'.format(log_ncmc, randnum) )
+                self.md_sim.context.setPositions(nc_state1['positions'])
+                if write_move:
+                    self.writeFrame(self.md_sim, '{}acc-it{}.pdb'.format(self.opt['outfname'],self.current_iter))
+
+            else:
+                self.reject += 1
+                self.log.info('NCMC MOVE REJECTED: log_ncmc {} < {}'.format(log_ncmc, randnum) )
+                self.nc_context.setPositions(md_state0['positions'])
+
+
+
 
         self.nc_integrator.reset()
+        self.move_engine.resetAcceptanceRatios()
         self.md_sim.context.setVelocitiesToTemperature(temperature)
 
     def simulateNCMC(self, nstepsNC=5000, ncmc_traj=None,
@@ -486,12 +499,18 @@ def simulateNCMC(self, nstepsNC=5000, ncmc_traj=None,
                 #Attempt anything related to the move before protocol is performed
                 if nc_step == 0:
                     self.nc_context = self.move_engine.moves[self.move_engine.selected_move].beforeMove(self.nc_context)
+                    #if the acceptance ratio is 0, stop running, since the move would be rejected anyway
+                    if self.move_engine.moves[self.move_engine.selected_move].acceptance_ratio == 0:
+                        break
                 # Attempt selected MoveEngine Move at the halfway point
                 #to ensure protocol is symmetric
                 if self.movestep == nc_step:
                     #Do move
                     self.log.info('Performing %s...' % move_name)
                     self.nc_context = self.move_engine.runEngine(self.nc_context)
+                    #if the acceptance ratio is 0, stop running, since the move would be rejected anyway
+                    if self.move_engine.moves[self.move_engine.selected_move].acceptance_ratio == 0:
+                        break
 
                 # Do 1 NCMC step with the integrator
                 self.nc_integrator.step(1)
@@ -571,8 +590,8 @@ def run(self, nIter=100):
             self.simulateMD(**self.opt)
 
         # END OF NITER
-        self.accept_ratio = self.accept/float(nIter)
-        self.log.info('Acceptance Ratio: %s' % self.accept_ratio)
+        self.accept_probability = self.accept/float(nIter)
+        self.log.info('Acceptance Probability: %s' % self.accept_probability)
         self.log.info('nIter: %s ' % nIter)
 
     def simulateMC(self):

blues/version.py

@@ -2,6 +2,6 @@
 # This file is automatically generated by setup.py
 short_version = '0.1.3'
 version = '0.1.3.dev1'
-full_version = '0.1.3.dev1-973bbc9'
-git_revision = '973bbc975a8a8223d5b8a49edc2f89207a8390eb'
+full_version = '0.1.3.dev1-47cb523'
+git_revision = '47cb5230a3748a7a0019af2b3e1ce393f51ba0ec'
 release = False