diff --git a/blues/engine.py b/blues/engine.py index 62c3a59f..62dfc641 100644 --- a/blues/engine.py +++ b/blues/engine.py @@ -74,4 +74,8 @@ def runEngine(self, context): print(e) raise SystemExit - return new_context \ No newline at end of file + return new_context + + def resetAcceptanceRatios(self): + for move in self.moves: + move.acceptance_ratio = 1.0 \ No newline at end of file diff --git a/blues/moves.py b/blues/moves.py index 5cb3df3a..460cdf00 100644 --- a/blues/moves.py +++ b/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 diff --git a/blues/simulation.py b/blues/simulation.py index 5631017b..ad95de74 100644 --- a/blues/simulation.py +++ b/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): diff --git a/blues/version.py b/blues/version.py index 291a313a..b8db839f 100644 --- a/blues/version.py +++ b/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