From a7fe6deb9dbd648a620aa3ddbf6c0457dce53a94 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 13:31:10 -0700 Subject: [PATCH 01/10] added alchemical splitting integrator --- blues/integrators.py | 243 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 243 insertions(+) create mode 100644 blues/integrators.py diff --git a/blues/integrators.py b/blues/integrators.py new file mode 100644 index 00000000..4ecc261f --- /dev/null +++ b/blues/integrators.py @@ -0,0 +1,243 @@ +from openmmtools.integrators import LangevinIntegrator +#TODO use Nonequilibrium baseclass directly +class AlchemicalLangevinSplittingIntegrator(LangevinIntegrator): + """Allows nonequilibrium switching based on force parameters specified in alchemical_functions. + Propagator is based on Langevin splitting, as described below. + One way to divide the Langevin system is into three parts which can each be solved "exactly:" + - R: Linear "drift" / Constrained "drift" + Deterministic update of *positions*, using current velocities + x <- x + v dt + - V: Linear "kick" / Constrained "kick" + Deterministic update of *velocities*, using current forces + v <- v + (f/m) dt + where f = force, m = mass + - O: Ornstein-Uhlenbeck + Stochastic update of velocities, simulating interaction with a heat bath + v <- av + b sqrt(kT/m) R + where + a = e^(-gamma dt) + b = sqrt(1 - e^(-2gamma dt)) + R is i.i.d. standard normal + We can then construct integrators by solving each part for a certain timestep in sequence. + (We can further split up the V step by force group, evaluating cheap but fast-fluctuating + forces more frequently than expensive but slow-fluctuating forces. Since forces are only + evaluated in the V step, we represent this by including in our "alphabet" V0, V1, ...) + When the system contains holonomic constraints, these steps are confined to the constraint + manifold. + Examples + -------- + - VVVR + splitting="O V R V O" + - BAOAB: + splitting="V R O R V" + - g-BAOAB, with K_r=3: + splitting="V R R R O R R R V" + - g-BAOAB with solvent-solute splitting, K_r=K_p=2: + splitting="V0 V1 R R O R R V1 R R O R R V1 V0" + Attributes + ---------- + _kinetic_energy : str + This is 0.5*m*v*v by default, and is the expression used for the kinetic energy + References + ---------- + [Nilmeier, et al. 2011] Nonequilibrium candidate Monte Carlo is an efficient tool for equilibrium simulation + [Leimkuhler and Matthews, 2015] Molecular dynamics: with deterministic and stochastic numerical methods, Chapter 7 + """ + + def __init__(self, + alchemical_functions, + splitting="V R O R V", + temperature=298.0 * simtk.unit.kelvin, + collision_rate=1.0 / simtk.unit.picoseconds, + timestep=1.0 * simtk.unit.femtoseconds, + constraint_tolerance=1e-8, + measure_shadow_work=False, + measure_heat=True, + direction="forward", + steps_per_propagation=1 + nsteps_neq=100): + """ + Parameters + ---------- + alchemical_functions : dict of strings + key: value pairs such as "global_parameter" : function_of_lambda where function_of_lambda is a Lepton-compatible + string that depends on the variable "lambda" + splitting : string, default: "V R O R V" + Sequence of R, V, O (and optionally V{i}), and { }substeps to be executed each timestep. + Forces are only used in V-step. Handle multiple force groups by appending the force group index + to V-steps, e.g. "V0" will only use forces from force group 0. "V" will perform a step using all forces. + ( will cause metropolization, and must be followed later by a ). + temperature : numpy.unit.Quantity compatible with kelvin, default: 298.0*simtk.unit.kelvin + Fictitious "bath" temperature + collision_rate : numpy.unit.Quantity compatible with 1/picoseconds, default: 91.0/simtk.unit.picoseconds + Collision rate + timestep : numpy.unit.Quantity compatible with femtoseconds, default: 1.0*simtk.unit.femtoseconds + Integration timestep + constraint_tolerance : float, default: 1.0e-8 + Tolerance for constraint solver + measure_shadow_work : boolean, default: False + Accumulate the shadow work performed by the symplectic substeps, in the global `shadow_work` + measure_heat : boolean, default: True + Accumulate the heat exchanged with the bath in each step, in the global `heat` + direction : str, default: "forward" + Whether to move the global lambda parameter from 0 to 1 (forward) or 1 to 0 (reverse). + nsteps_neq : int, default: 100 + Number of steps in nonequilibrium protocol. Default 100 + """ + + self._alchemical_functions = alchemical_functions + self._direction = direction + self._n_steps_neq = nsteps_neq + + # collect the system parameters. + self._system_parameters = {system_parameter for system_parameter in alchemical_functions.keys()} + + # call the base class constructor + super(AlchemicalLangevinSplittingIntegrator, self).__init__(splitting=splitting, temperature=temperature, + collision_rate=collision_rate, timestep=timestep, + constraint_tolerance=constraint_tolerance, + measure_shadow_work=measure_shadow_work, + measure_heat=measure_heat, + ) + + # add some global variables relevant to the integrator + self.add_global_variables(nsteps=nsteps_neq) + + def update_alchemical_parameters_step(self): + """ + Update Context parameters according to provided functions. + """ + for context_parameter in self._alchemical_functions: + if context_parameter in self._system_parameters: + self.addComputeGlobal(context_parameter, self._alchemical_functions[context_parameter]) + + def alchemical_perturbation_step(self): + """ + Add alchemical perturbation step, accumulating protocol work. + """ + # Store initial potential energy + self.addComputeGlobal("Eold", "Epert") + + # Use fractional state + if self._direction == 'forward': + self.addComputeGlobal('lambda', '(step+1)/nsteps') + elif self._direction == 'reverse': + self.addComputeGlobal('lambda', '(nsteps - step - 1)/nsteps') + + # Update all slaved alchemical parameters + self.update_alchemical_parameters_step() + + # Accumulate protocol work + self.addComputeGlobal("Enew", "energy") + self.addComputeGlobal("protocol_work", "protocol_work + (Enew-Eold)/kT") + + def sanity_check(self, splitting, allowed_characters="H{}RVO0123456789"): + super(AlchemicalLangevinSplittingIntegrator, self).sanity_check(splitting, allowed_characters=allowed_characters) + + def substep_function(self, step_string, measure_shadow_work, measure_heat, n_R, force_group_nV, mts): + """Take step string, and add the appropriate R, V, O, M, or H step with appropriate parameters. + The step string input here is a single character (or character + number, for MTS) + Parameters + ---------- + step_string : str + R, O, V, or Vn (where n is a nonnegative integer specifying force group) + measure_shadow_work : bool + Whether the steps should measure shadow work + measure_heat : bool + Whether the O step should measure heat + n_R : int + The number of R steps per integrator step + force_group_nV : dict + The number of V steps per integrator step per force group. {0: nV} if not mts + mts : bool + Whether the integrator is a multiple timestep integrator + """ + + if step_string == "O": + self.O_step(measure_heat) + elif step_string == "R": + self.R_step(measure_shadow_work, n_R) + elif step_string == "{": + self.begin_metropolize() + elif step_string == "}": + self.metropolize() + elif step_string[0] == "V": + # get the force group for this update--it's the number after the V + force_group = step_string[1:] + self.V_step(force_group, measure_shadow_work, force_group_nV, mts) + elif step_string == "H": + self.alchemical_perturbation_step() + + def add_integrator_steps(self, splitting, measure_shadow_work, measure_heat, ORV_counts, force_group_nV, mts): + """ + Override the base class to insert reset steps around the integrator. + """ + #if the step is zero, + self.beginIfBlock('step = 0') + self.addConstrainPositions() + self.addConstrainVelocities() + self.reset_work_step() + self.alchemical_reset_step() + self.addComputeGlobal("Epert", "energy") + self.endBlock() + + #call the superclass function to insert the appropriate steps, provided the step number is less than n_steps + self.beginIfBlock("step < nsteps") + super(AlchemicalLangevinSplittingIntegrator, self).add_integrator_steps(splitting, measure_shadow_work, + measure_heat, ORV_counts, + force_group_nV, mts) + + #increment the step number + self.addComputeGlobal("step", "step + 1") + + self.endBlock() + + + def add_global_variables(self, nsteps): + """Add the appropriate global parameters to the CustomIntegrator. nsteps refers to the number of + total steps in the protocol. + Parameters + ---------- + nsteps : int, greater than 0 + The number of steps in the switching protocol. + """ + self.addGlobalVariable('Eold', 0) #old energy value before perturbation + self.addGlobalVariable('Enew', 0) #new energy value after perturbation + self.addGlobalVariable('Epert', 0) #holder energy value after integrator step to keep track of non-alchemical work + self.addGlobalVariable('lambda', 0.0) # parameter switched from 0 <--> 1 during course of integrating internal 'nsteps' of dynamics + self.addGlobalVariable('kinetic', 0.0) # kinetic energy + self.addGlobalVariable('nsteps', nsteps) # total number of NCMC steps to perform + self.addGlobalVariable('step', 0) # current NCMC step number + self.addGlobalVariable('protocol_work', 0) # work performed by NCMC protocol + self.addGlobalVariable('pstep') + + def alchemical_reset_step(self): + """ + Reset the alchemical lambda to its starting value + This is 1 for reverse and 0 for forward + """ + if self._direction == "forward": + self.addComputeGlobal("lambda", "0") + if self._direction == "reverse": + self.addComputeGlobal("lambda", "1") + + self.addComputeGlobal("protocol_work", "0.0") + if self._measure_shadow_work: + self.addComputeGlobal("shadow_work", "0.0") + self.addComputeGlobal("step", "0.0") + #add all dependent parameters + self.update_alchemical_parameters_step() + + def reset_work_step(self): + """ + This step resets work statistics that have been accumulated. + """ + self.addComputeGlobal("protocol_work", "0.0") + self.addComputeGlobal("shadow_work", "0.0") + + def reset_integrator(self): + """ + Manually reset the work statistics and step + """ + self.setGlobalVariableByName("step", 0) + From 98f35b2e09dc639711bc35853651cddeb85da62f Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 13:53:18 -0700 Subject: [PATCH 02/10] fixed import statements --- blues/integrators.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/blues/integrators.py b/blues/integrators.py index 4ecc261f..72333e12 100644 --- a/blues/integrators.py +++ b/blues/integrators.py @@ -1,4 +1,15 @@ from openmmtools.integrators import LangevinIntegrator +import numpy +import logging +import re + +import simtk.unit + +import simtk.unit as units +import simtk.openmm as mm + +from openmmtools.constants import kB + #TODO use Nonequilibrium baseclass directly class AlchemicalLangevinSplittingIntegrator(LangevinIntegrator): """Allows nonequilibrium switching based on force parameters specified in alchemical_functions. @@ -54,7 +65,7 @@ def __init__(self, measure_shadow_work=False, measure_heat=True, direction="forward", - steps_per_propagation=1 + steps_per_propagation=1, nsteps_neq=100): """ Parameters From 40628a3193316b75b857bff9994394be39828e10 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 16:11:27 -0700 Subject: [PATCH 03/10] changed to inherit from base ncmcintegrator class --- blues/integrators.py | 237 ++++++++++++------------------------------- 1 file changed, 66 insertions(+), 171 deletions(-) diff --git a/blues/integrators.py b/blues/integrators.py index 72333e12..4dde06f2 100644 --- a/blues/integrators.py +++ b/blues/integrators.py @@ -1,18 +1,10 @@ -from openmmtools.integrators import LangevinIntegrator -import numpy -import logging -import re +from openmmtools.integrators import NonequilibriumLangevinIntegrator +import simtk -import simtk.unit - -import simtk.unit as units -import simtk.openmm as mm - -from openmmtools.constants import kB - -#TODO use Nonequilibrium baseclass directly -class AlchemicalLangevinSplittingIntegrator(LangevinIntegrator): +class NonequilibriumExternalLangevinIntegrator(NonequilibriumLangevinIntegrator): """Allows nonequilibrium switching based on force parameters specified in alchemical_functions. + A variable named lambda is switched from 0 to 1 linearly throughout the nsteps of the protocol. + The functions can use this to create more complex protocols for other global parameters. Propagator is based on Langevin splitting, as described below. One way to divide the Langevin system is into three parts which can each be solved "exactly:" - R: Linear "drift" / Constrained "drift" @@ -37,14 +29,14 @@ class AlchemicalLangevinSplittingIntegrator(LangevinIntegrator): manifold. Examples -------- + - g-BAOAB: + splitting="R V O H O V R" - VVVR - splitting="O V R V O" - - BAOAB: - splitting="V R O R V" - - g-BAOAB, with K_r=3: - splitting="V R R R O R R R V" - - g-BAOAB with solvent-solute splitting, K_r=K_p=2: - splitting="V0 V1 R R O R R V1 R R O R R V1 V0" + splitting="O V R H R V O" + - VV + splitting="V R H R V" + - An NCMC algorithm with Metropolized integrator: + splitting="O { V R H R V } O" Attributes ---------- _kinetic_energy : str @@ -57,24 +49,24 @@ class AlchemicalLangevinSplittingIntegrator(LangevinIntegrator): def __init__(self, alchemical_functions, - splitting="V R O R V", + splitting="R V O H O V R", temperature=298.0 * simtk.unit.kelvin, collision_rate=1.0 / simtk.unit.picoseconds, timestep=1.0 * simtk.unit.femtoseconds, constraint_tolerance=1e-8, measure_shadow_work=False, measure_heat=True, - direction="forward", - steps_per_propagation=1, - nsteps_neq=100): + nsteps_neq=100, + steps_per_propagation=1): """ Parameters ---------- alchemical_functions : dict of strings key: value pairs such as "global_parameter" : function_of_lambda where function_of_lambda is a Lepton-compatible string that depends on the variable "lambda" - splitting : string, default: "V R O R V" - Sequence of R, V, O (and optionally V{i}), and { }substeps to be executed each timestep. + splitting : string, default: "H V R O V R H" + Sequence of R, V, O (and optionally V{i}), and { }substeps to be executed each timestep. There is also an H option, + which increments the global parameter `lambda` by 1/nsteps_neq for each step. Forces are only used in V-step. Handle multiple force groups by appending the force group index to V-steps, e.g. "V0" will only use forces from force group 0. "V" will perform a step using all forces. ( will cause metropolization, and must be followed later by a ). @@ -90,165 +82,68 @@ def __init__(self, Accumulate the shadow work performed by the symplectic substeps, in the global `shadow_work` measure_heat : boolean, default: True Accumulate the heat exchanged with the bath in each step, in the global `heat` - direction : str, default: "forward" - Whether to move the global lambda parameter from 0 to 1 (forward) or 1 to 0 (reverse). nsteps_neq : int, default: 100 Number of steps in nonequilibrium protocol. Default 100 """ - self._alchemical_functions = alchemical_functions - self._direction = direction - self._n_steps_neq = nsteps_neq +# self._alchemical_functions = alchemical_functions +# self._n_steps_neq = nsteps_neq # collect the system parameters. - self._system_parameters = {system_parameter for system_parameter in alchemical_functions.keys()} +# self._system_parameters = {system_parameter for system_parameter in alchemical_functions.keys()} # call the base class constructor - super(AlchemicalLangevinSplittingIntegrator, self).__init__(splitting=splitting, temperature=temperature, - collision_rate=collision_rate, timestep=timestep, - constraint_tolerance=constraint_tolerance, - measure_shadow_work=measure_shadow_work, - measure_heat=measure_heat, - ) + super(EditIntegrator, self).__init__(alchemical_functions=alchemical_functions, + splitting=splitting, temperature=temperature, + collision_rate=collision_rate, timestep=timestep, + constraint_tolerance=constraint_tolerance, + measure_shadow_work=measure_shadow_work, + measure_heat=measure_heat, + nsteps_neq=nsteps_neq + ) # add some global variables relevant to the integrator - self.add_global_variables(nsteps=nsteps_neq) - - def update_alchemical_parameters_step(self): - """ - Update Context parameters according to provided functions. - """ - for context_parameter in self._alchemical_functions: - if context_parameter in self._system_parameters: - self.addComputeGlobal(context_parameter, self._alchemical_functions[context_parameter]) - - def alchemical_perturbation_step(self): - """ - Add alchemical perturbation step, accumulating protocol work. - """ - # Store initial potential energy - self.addComputeGlobal("Eold", "Epert") - - # Use fractional state - if self._direction == 'forward': - self.addComputeGlobal('lambda', '(step+1)/nsteps') - elif self._direction == 'reverse': - self.addComputeGlobal('lambda', '(nsteps - step - 1)/nsteps') - - # Update all slaved alchemical parameters - self.update_alchemical_parameters_step() - - # Accumulate protocol work - self.addComputeGlobal("Enew", "energy") - self.addComputeGlobal("protocol_work", "protocol_work + (Enew-Eold)/kT") - - def sanity_check(self, splitting, allowed_characters="H{}RVO0123456789"): - super(AlchemicalLangevinSplittingIntegrator, self).sanity_check(splitting, allowed_characters=allowed_characters) - - def substep_function(self, step_string, measure_shadow_work, measure_heat, n_R, force_group_nV, mts): - """Take step string, and add the appropriate R, V, O, M, or H step with appropriate parameters. - The step string input here is a single character (or character + number, for MTS) - Parameters - ---------- - step_string : str - R, O, V, or Vn (where n is a nonnegative integer specifying force group) - measure_shadow_work : bool - Whether the steps should measure shadow work - measure_heat : bool - Whether the O step should measure heat - n_R : int - The number of R steps per integrator step - force_group_nV : dict - The number of V steps per integrator step per force group. {0: nV} if not mts - mts : bool - Whether the integrator is a multiple timestep integrator - """ - - if step_string == "O": - self.O_step(measure_heat) - elif step_string == "R": - self.R_step(measure_shadow_work, n_R) - elif step_string == "{": - self.begin_metropolize() - elif step_string == "}": - self.metropolize() - elif step_string[0] == "V": - # get the force group for this update--it's the number after the V - force_group = step_string[1:] - self.V_step(force_group, measure_shadow_work, force_group_nV, mts) - elif step_string == "H": - self.alchemical_perturbation_step() + ##self.add_global_variables(nsteps=nsteps_neq) + kB = simtk.unit.BOLTZMANN_CONSTANT_kB * simtk.unit.AVOGADRO_CONSTANT_NA + kT = kB * temperature + self.kT = kT + self.addGlobalVariable("perturbed_pe", 0) + self.addGlobalVariable("unperturbed_pe", 0) + self.addGlobalVariable("first_step", 0) + self.addGlobalVariable("psteps", steps_per_propagation) + self.addGlobalVariable("pstep", 0) def add_integrator_steps(self, splitting, measure_shadow_work, measure_heat, ORV_counts, force_group_nV, mts): - """ - Override the base class to insert reset steps around the integrator. - """ - #if the step is zero, - self.beginIfBlock('step = 0') - self.addConstrainPositions() - self.addConstrainVelocities() - self.reset_work_step() - self.alchemical_reset_step() - self.addComputeGlobal("Epert", "energy") + self.addComputeGlobal("perturbed_pe", "energy") + # Assumes no perturbation is done before doing the initial MD step. + self.beginIfBlock("first_step < 1") + self.addComputeGlobal("first_step", "1") + self.addComputeGlobal("unperturbed_pe", "energy") self.endBlock() - - #call the superclass function to insert the appropriate steps, provided the step number is less than n_steps - self.beginIfBlock("step < nsteps") - super(AlchemicalLangevinSplittingIntegrator, self).add_integrator_steps(splitting, measure_shadow_work, - measure_heat, ORV_counts, - force_group_nV, mts) - - #increment the step number - self.addComputeGlobal("step", "step + 1") - + self.addComputeGlobal("perturbed_pe", "energy") + self.addComputeGlobal("protocol_work", "protocol_work + (perturbed_pe - unperturbed_pe)") + #repeat BAOAB integration psteps number of times per lambda + self.addComputeGlobal("pstep", "0") + self.beginWhileBlock('pstep < psteps') + super(EditIntegrator, self).add_integrator_steps(splitting, measure_shadow_work, measure_heat, ORV_counts, force_group_nV, mts) + self.beginIfBlock("pstep < psteps - 1") + self.addComputeGlobal("step", "step - 1") self.endBlock() + self.addComputeGlobal("pstep", "pstep + 1") + self.endBlock() + #update unperturbed_pe + self.addComputeGlobal("unperturbed_pe", "energy") + def getLogAcceptanceProbability(self, context): + protocol = self.getGlobalVariableByName("protocol_work") + shadow = self.getGlobalVariableByName("shadow_work") + logp_accept = -1.0*(protocol + shadow) + return logp_accept - def add_global_variables(self, nsteps): - """Add the appropriate global parameters to the CustomIntegrator. nsteps refers to the number of - total steps in the protocol. - Parameters - ---------- - nsteps : int, greater than 0 - The number of steps in the switching protocol. - """ - self.addGlobalVariable('Eold', 0) #old energy value before perturbation - self.addGlobalVariable('Enew', 0) #new energy value after perturbation - self.addGlobalVariable('Epert', 0) #holder energy value after integrator step to keep track of non-alchemical work - self.addGlobalVariable('lambda', 0.0) # parameter switched from 0 <--> 1 during course of integrating internal 'nsteps' of dynamics - self.addGlobalVariable('kinetic', 0.0) # kinetic energy - self.addGlobalVariable('nsteps', nsteps) # total number of NCMC steps to perform - self.addGlobalVariable('step', 0) # current NCMC step number - self.addGlobalVariable('protocol_work', 0) # work performed by NCMC protocol - self.addGlobalVariable('pstep') - - def alchemical_reset_step(self): - """ - Reset the alchemical lambda to its starting value - This is 1 for reverse and 0 for forward - """ - if self._direction == "forward": - self.addComputeGlobal("lambda", "0") - if self._direction == "reverse": - self.addComputeGlobal("lambda", "1") - - self.addComputeGlobal("protocol_work", "0.0") - if self._measure_shadow_work: - self.addComputeGlobal("shadow_work", "0.0") - self.addComputeGlobal("step", "0.0") - #add all dependent parameters - self.update_alchemical_parameters_step() - - def reset_work_step(self): - """ - This step resets work statistics that have been accumulated. - """ - self.addComputeGlobal("protocol_work", "0.0") - self.addComputeGlobal("shadow_work", "0.0") - - def reset_integrator(self): - """ - Manually reset the work statistics and step - """ + def reset(self): self.setGlobalVariableByName("step", 0) - + self.setGlobalVariableByName("lambda", 0.0) +# self.setGlobalVariableByName("total_work", 0.0) + self.setGlobalVariableByName("protocol_work", 0.0) + self.setGlobalVariableByName("shadow_work", 0.0) + self.setGlobalVariableByName("first_step", 0) From 1de9170bc3937aa6f738fc9dfb3b630fcff9bb28 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 17:21:58 -0700 Subject: [PATCH 04/10] fixed super self references --- blues/integrators.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/blues/integrators.py b/blues/integrators.py index 4dde06f2..b5c17ca8 100644 --- a/blues/integrators.py +++ b/blues/integrators.py @@ -93,7 +93,7 @@ def __init__(self, # self._system_parameters = {system_parameter for system_parameter in alchemical_functions.keys()} # call the base class constructor - super(EditIntegrator, self).__init__(alchemical_functions=alchemical_functions, + super(NonequilibriumExternalLangevinIntegrator, self).__init__(alchemical_functions=alchemical_functions, splitting=splitting, temperature=temperature, collision_rate=collision_rate, timestep=timestep, constraint_tolerance=constraint_tolerance, @@ -112,6 +112,10 @@ def __init__(self, self.addGlobalVariable("first_step", 0) self.addGlobalVariable("psteps", steps_per_propagation) self.addGlobalVariable("pstep", 0) + try: + self.getGlobalVariableByName("shadow_work") + except: + self.addGlobalVariable('shadow_work', 0) def add_integrator_steps(self, splitting, measure_shadow_work, measure_heat, ORV_counts, force_group_nV, mts): self.addComputeGlobal("perturbed_pe", "energy") @@ -125,7 +129,7 @@ def add_integrator_steps(self, splitting, measure_shadow_work, measure_heat, ORV #repeat BAOAB integration psteps number of times per lambda self.addComputeGlobal("pstep", "0") self.beginWhileBlock('pstep < psteps') - super(EditIntegrator, self).add_integrator_steps(splitting, measure_shadow_work, measure_heat, ORV_counts, force_group_nV, mts) + super(NonequilibriumExternalLangevinIntegrator, self).add_integrator_steps(splitting, measure_shadow_work, measure_heat, ORV_counts, force_group_nV, mts) self.beginIfBlock("pstep < psteps - 1") self.addComputeGlobal("step", "step - 1") self.endBlock() From 5417f31c008007f46934c2c553072eca68bbf39c Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 17:23:15 -0700 Subject: [PATCH 05/10] added test for new integrator --- blues/tests/test_integrator.py | 36 ++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 blues/tests/test_integrator.py diff --git a/blues/tests/test_integrator.py b/blues/tests/test_integrator.py new file mode 100644 index 00000000..4f6fff6a --- /dev/null +++ b/blues/tests/test_integrator.py @@ -0,0 +1,36 @@ +from openmmtools import testsystems, alchemy +from openmmtools.alchemy import AlchemicalFactory, AlchemicalRegion +from simtk.openmm.app import Simulation +import simtk +from integrators import NonequilibriumExternalLangevinIntegrator +testsystem = testsystems.AlanineDipeptideVacuum() +functions = { 'lambda_sterics' : '1', 'lambda_electrostatics' : '1'} +factory = AlchemicalFactory(consistent_exceptions=False) +alanine_vacuum = testsystem.system +alchemical_region = AlchemicalRegion(alchemical_atoms=range(22), + annihilate_electrostatics=True, annihilate_sterics=True) +alanine_alchemical_system = factory.create_alchemical_system(reference_system=alanine_vacuum, + alchemical_regions=alchemical_region) +nsteps_neq = 6 +integrator = NonequilibriumExternalLangevinIntegrator(alchemical_functions=functions, + timestep=0.05*simtk.unit.femtoseconds, + nsteps_neq=nsteps_neq, + measure_shadow_work=True) +for i in range(20): + print(integrator.getGlobalVariableName(i)) + +simulation = Simulation(testsystem.topology, alanine_alchemical_system, integrator) +simulation.context.setPositions(testsystem.positions) +for i in range(nsteps_neq): + simulation.step(1) + protocol_work = simulation.integrator.getGlobalVariableByName("protocol_work") + if i == 3: + state = simulation.context.getState(getPositions=True) + pos = state.getPositions(asNumpy=True) + pos[0,1] = pos[0,1] + 0.5*simtk.unit.nanometers + simulation.context.setPositions(pos) +print(protocol_work) +#protocol work is around 550.0 +assert 548.0 < protocol_work +assert 551.0 > protocol_work + From c6e6aa4fde18a7945493c078b2db2df79c797f56 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 17:25:48 -0700 Subject: [PATCH 06/10] added function to actually test --- blues/tests/test_integrator.py | 60 ++++++++++++++++++---------------- 1 file changed, 31 insertions(+), 29 deletions(-) diff --git a/blues/tests/test_integrator.py b/blues/tests/test_integrator.py index 4f6fff6a..43f576a1 100644 --- a/blues/tests/test_integrator.py +++ b/blues/tests/test_integrator.py @@ -3,34 +3,36 @@ from simtk.openmm.app import Simulation import simtk from integrators import NonequilibriumExternalLangevinIntegrator -testsystem = testsystems.AlanineDipeptideVacuum() -functions = { 'lambda_sterics' : '1', 'lambda_electrostatics' : '1'} -factory = AlchemicalFactory(consistent_exceptions=False) -alanine_vacuum = testsystem.system -alchemical_region = AlchemicalRegion(alchemical_atoms=range(22), - annihilate_electrostatics=True, annihilate_sterics=True) -alanine_alchemical_system = factory.create_alchemical_system(reference_system=alanine_vacuum, - alchemical_regions=alchemical_region) -nsteps_neq = 6 -integrator = NonequilibriumExternalLangevinIntegrator(alchemical_functions=functions, - timestep=0.05*simtk.unit.femtoseconds, - nsteps_neq=nsteps_neq, - measure_shadow_work=True) -for i in range(20): - print(integrator.getGlobalVariableName(i)) +def test_nonequilibrium_external_integrator(): + testsystem = testsystems.AlanineDipeptideVacuum() + functions = { 'lambda_sterics' : '1', 'lambda_electrostatics' : '1'} + factory = AlchemicalFactory(consistent_exceptions=False) + alanine_vacuum = testsystem.system + alchemical_region = AlchemicalRegion(alchemical_atoms=range(22), + annihilate_electrostatics=True, annihilate_sterics=True) + alanine_alchemical_system = factory.create_alchemical_system(reference_system=alanine_vacuum, + alchemical_regions=alchemical_region) + nsteps_neq = 6 + integrator = NonequilibriumExternalLangevinIntegrator(alchemical_functions=functions, + timestep=0.05*simtk.unit.femtoseconds, + nsteps_neq=nsteps_neq, + measure_shadow_work=True) + for i in range(20): + print(integrator.getGlobalVariableName(i)) + + simulation = Simulation(testsystem.topology, alanine_alchemical_system, integrator) + simulation.context.setPositions(testsystem.positions) + for i in range(nsteps_neq): + simulation.step(1) + protocol_work = simulation.integrator.getGlobalVariableByName("protocol_work") + if i == 3: + state = simulation.context.getState(getPositions=True) + pos = state.getPositions(asNumpy=True) + pos[0,1] = pos[0,1] + 0.5*simtk.unit.nanometers + simulation.context.setPositions(pos) + print(protocol_work) + #protocol work is around 550.0 + assert 548.0 < protocol_work + assert 551.0 > protocol_work -simulation = Simulation(testsystem.topology, alanine_alchemical_system, integrator) -simulation.context.setPositions(testsystem.positions) -for i in range(nsteps_neq): - simulation.step(1) - protocol_work = simulation.integrator.getGlobalVariableByName("protocol_work") - if i == 3: - state = simulation.context.getState(getPositions=True) - pos = state.getPositions(asNumpy=True) - pos[0,1] = pos[0,1] + 0.5*simtk.unit.nanometers - simulation.context.setPositions(pos) -print(protocol_work) -#protocol work is around 550.0 -assert 548.0 < protocol_work -assert 551.0 > protocol_work From 839c8911704ac933e4d3daf18856363322af7064 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 17:31:28 -0700 Subject: [PATCH 07/10] fixed import --- blues/tests/test_integrator.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/blues/tests/test_integrator.py b/blues/tests/test_integrator.py index 43f576a1..739a2368 100644 --- a/blues/tests/test_integrator.py +++ b/blues/tests/test_integrator.py @@ -2,7 +2,7 @@ from openmmtools.alchemy import AlchemicalFactory, AlchemicalRegion from simtk.openmm.app import Simulation import simtk -from integrators import NonequilibriumExternalLangevinIntegrator +from blues.integrators import NonequilibriumExternalLangevinIntegrator def test_nonequilibrium_external_integrator(): testsystem = testsystems.AlanineDipeptideVacuum() functions = { 'lambda_sterics' : '1', 'lambda_electrostatics' : '1'} From c91bc7c17f6bf1fc8babce757a4c7fbf5be04402 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 17:51:42 -0700 Subject: [PATCH 08/10] checks if steps_per_propogation is behaving --- blues/tests/test_integrator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/blues/tests/test_integrator.py b/blues/tests/test_integrator.py index 739a2368..395e67fb 100644 --- a/blues/tests/test_integrator.py +++ b/blues/tests/test_integrator.py @@ -16,7 +16,8 @@ def test_nonequilibrium_external_integrator(): integrator = NonequilibriumExternalLangevinIntegrator(alchemical_functions=functions, timestep=0.05*simtk.unit.femtoseconds, nsteps_neq=nsteps_neq, - measure_shadow_work=True) + measure_shadow_work=False, + steps_per_propagation=2) for i in range(20): print(integrator.getGlobalVariableName(i)) @@ -35,4 +36,3 @@ def test_nonequilibrium_external_integrator(): assert 548.0 < protocol_work assert 551.0 > protocol_work - From 5f3c673c6903e623b5a68903de50abe27bd7d385 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Mon, 24 Apr 2017 17:55:17 -0700 Subject: [PATCH 09/10] update doc string to reflect difference with orig --- blues/integrators.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/blues/integrators.py b/blues/integrators.py index b5c17ca8..ca938ce8 100644 --- a/blues/integrators.py +++ b/blues/integrators.py @@ -5,6 +5,8 @@ class NonequilibriumExternalLangevinIntegrator(NonequilibriumLangevinIntegrator) """Allows nonequilibrium switching based on force parameters specified in alchemical_functions. A variable named lambda is switched from 0 to 1 linearly throughout the nsteps of the protocol. The functions can use this to create more complex protocols for other global parameters. + This also takes into account work done outside the nonequilibrium switching between steps, + for example the work done if a molecule is rotated. Propagator is based on Langevin splitting, as described below. One way to divide the Langevin system is into three parts which can each be solved "exactly:" - R: Linear "drift" / Constrained "drift" From e1ffb8f9f210a8070c4076e268d4812c82324e45 Mon Sep 17 00:00:00 2001 From: sgill2 Date: Tue, 25 Apr 2017 11:02:08 -0700 Subject: [PATCH 10/10] made energy consistent (in kT) for external work --- blues/integrators.py | 2 +- blues/tests/test_integrator.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/blues/integrators.py b/blues/integrators.py index ca938ce8..293c828f 100644 --- a/blues/integrators.py +++ b/blues/integrators.py @@ -127,7 +127,7 @@ def add_integrator_steps(self, splitting, measure_shadow_work, measure_heat, ORV self.addComputeGlobal("unperturbed_pe", "energy") self.endBlock() self.addComputeGlobal("perturbed_pe", "energy") - self.addComputeGlobal("protocol_work", "protocol_work + (perturbed_pe - unperturbed_pe)") + self.addComputeGlobal("protocol_work", "protocol_work + (perturbed_pe - unperturbed_pe)/kT") #repeat BAOAB integration psteps number of times per lambda self.addComputeGlobal("pstep", "0") self.beginWhileBlock('pstep < psteps') diff --git a/blues/tests/test_integrator.py b/blues/tests/test_integrator.py index 395e67fb..73266ee1 100644 --- a/blues/tests/test_integrator.py +++ b/blues/tests/test_integrator.py @@ -18,8 +18,6 @@ def test_nonequilibrium_external_integrator(): nsteps_neq=nsteps_neq, measure_shadow_work=False, steps_per_propagation=2) - for i in range(20): - print(integrator.getGlobalVariableName(i)) simulation = Simulation(testsystem.topology, alanine_alchemical_system, integrator) simulation.context.setPositions(testsystem.positions) @@ -32,7 +30,9 @@ def test_nonequilibrium_external_integrator(): pos[0,1] = pos[0,1] + 0.5*simtk.unit.nanometers simulation.context.setPositions(pos) print(protocol_work) - #protocol work is around 550.0 - assert 548.0 < protocol_work - assert 551.0 > protocol_work + #protocol work is around 221.0 + assert 220.0 < protocol_work + assert 223.0 > protocol_work +if __name__ == '__main__': + test_nonequilibrium_external_integrator()