diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 43b816e9b..045bef31e 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -18,7 +18,8 @@ requirements: - numexpr - autograd - pymbar - - openmm >=7.2.0 + - cuda92 + - openmm ==7.3.0 - parmed - openmoltools - alchemy >=1.2.3 @@ -43,7 +44,8 @@ requirements: - numexpr - autograd - pymbar - - openmm >=7.2.0 + - cuda92 + - openmm ==7.3.0 - parmed - openmoltools - alchemy >=1.2.3 diff --git a/examples/cdk2-example/cdk2_setup.yaml b/examples/cdk2-example/cdk2_setup.yaml index 79c53149c..3691f5ed4 100644 --- a/examples/cdk2-example/cdk2_setup.yaml +++ b/examples/cdk2-example/cdk2_setup.yaml @@ -13,8 +13,8 @@ new_ligand_index: 15 #be provided with a full path forcefield_files: - gaff.xml - - amber99sbildn.xml - - tip3p.xml + - amber14/protein.ff14SB.xml + - amber14/tip3p.xml #the temperature and pressure of the simulation, as well as how much solvent paddding to add #units: @@ -31,7 +31,9 @@ save_setup_pickle_as: fesetup_hbonds.pkl #the type of free energy calculation. #currently, this could be either nonequilibrium or sams -fe_type: nonequilibrium +fe_type: sams + +n_states: 100 #the forward switching functions. The reverse ones will be computed from this. Only valid for nonequilibrium switching forward_functions: @@ -42,17 +44,17 @@ forward_functions: lambda_torsions: lambda #the number of equilibration iterations: -n_equilibration_iterations: 100 +n_equilibration_iterations: 10 #The number of equilibrium steps to take between nonequilibrium switching events (only valid for nonequiibrium switching) -n_equilibrium_steps_per_iteration: 100 +n_equilibrium_steps_per_iteration: 1000 #The length of the ncmc protocol (only valid for nonequilibrium switching) n_steps_ncmc_protocol: 50 #The number of NCMC steps per move application. This controls the output frequency #1 step/move application means writing out the work at every step. (only valid for nonequilibrium switching) -n_steps_per_move_application: 10 +n_steps_per_move_application: 2500 #where to put the trajectories trajectory_directory: cdk2_neq_hbonds @@ -66,6 +68,7 @@ atom_selection: not water #which phases do we want to run (solvent, complex, or both solvent and complex in the list are acceptable): phases: - solvent + - complex #timestep in fs timestep: 1.0 @@ -83,5 +86,5 @@ measure_shadow_work: True scheduler_address: null #how many iterations to run (n_cycles*n_iterations_per_cycle) (only valid for nonequilibrium switching) -n_cycles: 5 +n_cycles: 10000 n_iterations_per_cycle: 1 \ No newline at end of file diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index c7954ac94..0c6562ff3 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -45,7 +45,7 @@ class HybridTopologyFactory(object): _known_forces = {'HarmonicBondForce', 'HarmonicAngleForce', 'PeriodicTorsionForce', 'NonbondedForce', 'MonteCarloBarostat'} _known_softcore_methods = ['default', 'amber', 'classic'] - def __init__(self, topology_proposal, current_positions, new_positions, use_dispersion_correction=False, functions=None, softcore_method='default', softcore_alpha=None, softcore_beta=None): + def __init__(self, topology_proposal, current_positions, new_positions, use_dispersion_correction=False, functions=None, softcore_method='amber', softcore_alpha=None, softcore_beta=None): """ Initialize the Hybrid topology factory. @@ -73,6 +73,7 @@ def __init__(self, topology_proposal, current_positions, new_positions, use_disp "alpha" parameter of softcore sterics. If None is provided, value will be set to 0.5 softcore_beta: unit, default None "beta" parameter of softcore electrostatics. If None is provided, value will be set to 12*unit.angstrom**2 + """ self._topology_proposal = topology_proposal self._old_system = copy.deepcopy(topology_proposal.old_system) @@ -84,17 +85,17 @@ def __init__(self, topology_proposal, current_positions, new_positions, use_disp self._new_positions = new_positions self._use_dispersion_correction = use_dispersion_correction - + if softcore_alpha is None: self.softcore_alpha = 0.5 else: self.softcore_alpha = softcore_alpha - + if softcore_beta is None: self.softcore_beta = 12*unit.angstrom**2 else: self.softcore_beta = softcore_beta - + if softcore_method not in self._known_softcore_methods: raise ValueError("Softcore method {} is not a valid method. Acceptable options are default, amber, and classic".format(softcore_method)) @@ -700,32 +701,6 @@ def _add_nonbonded_force_terms(self): custom_nonbonded_method = self._translate_nonbonded_method_to_custom(self._nonbonded_method) - # Create CustomNonbondedForce to handle interactions between alchemically-modified atoms and rest of system. - total_electrostatics_energy = "U_electrostatics;" + electrostatics_energy_expression + electrostatics_mixing_rules - if self._has_functions: - try: - total_electrostatics_energy += 'lambda_electrostatics = ' + self._functions['lambda_electrostatics'] - except KeyError as e: - print("Functions were provided, but there is no entry for electrostatics") - raise e - - electrostatics_custom_nonbonded_force = openmm.CustomNonbondedForce(total_electrostatics_energy) - electrostatics_custom_nonbonded_force.addGlobalParameter("softcore_beta", self.softcore_beta) - electrostatics_custom_nonbonded_force.addPerParticleParameter("chargeA") # partial charge initial - electrostatics_custom_nonbonded_force.addPerParticleParameter("chargeB") # partial charge final - - if self._has_functions: - electrostatics_custom_nonbonded_force.addGlobalParameter("lambda", 0.0) - electrostatics_custom_nonbonded_force.addEnergyParameterDerivative('lambda') - else: - electrostatics_custom_nonbonded_force.addGlobalParameter("lambda_electrostatics", 0.0) - - - electrostatics_custom_nonbonded_force.setNonbondedMethod(custom_nonbonded_method) - - self._hybrid_system.addForce(electrostatics_custom_nonbonded_force) - self._hybrid_system_forces['core_electrostatics_force'] = electrostatics_custom_nonbonded_force - total_sterics_energy = "U_sterics;" + sterics_energy_expression + sterics_mixing_rules if self._has_functions: try: @@ -767,18 +742,10 @@ def _add_nonbonded_force_terms(self): standard_nonbonded_force.setSwitchingDistance(switching_distance) sterics_custom_nonbonded_force.setUseSwitchingFunction(True) sterics_custom_nonbonded_force.setSwitchingDistance(switching_distance) - electrostatics_custom_nonbonded_force.setUseSwitchingFunction(True) - electrostatics_custom_nonbonded_force.setSwitchingDistance(switching_distance) else: standard_nonbonded_force.setUseSwitchingFunction(False) - electrostatics_custom_nonbonded_force.setUseSwitchingFunction(False) sterics_custom_nonbonded_force.setUseSwitchingFunction(False) - #Add a CustomBondForce for exceptions: - custom_nonbonded_bond_force = self._nonbonded_custom_bond_force(sterics_energy_expression, electrostatics_energy_expression) - self._hybrid_system.addForce(custom_nonbonded_bond_force) - self._hybrid_system_forces['core_nonbonded_bond_force'] = custom_nonbonded_bond_force - def _nonbonded_custom_sterics_common(self): """ Get a custom sterics expression that is common to all nonbonded methods @@ -937,48 +904,6 @@ def _nonbonded_custom_mixing_rules(self): electrostatics_mixing_rules += "chargeprodB = chargeB1*chargeB2;" # mixing rule for charges return sterics_mixing_rules, electrostatics_mixing_rules - def _nonbonded_custom_bond_force(self, sterics_energy_expression, electrostatics_energy_expression): - """ - Add a CustomBondForce to represent the exceptions in the NonbondedForce - - Parameters - ---------- - sterics_energy_expression : str - The complete energy expression being used for sterics - electrostatics_energy_expression : str - The complete energy expression being used for electrostatics - - Returns - ------- - custom_bond_force : openmm.CustomBondForce - The custom bond force for the nonbonded exceptions - """ - #Create the force and add its relevant parameters. - #we don't need to check that the keys exist, since by the time this is called, these are already checked. - if self._has_functions: - sterics_energy_expression += 'lambda_sterics = ' + self._functions['lambda_sterics'] - electrostatics_energy_expression += 'lambda_electrostatics = ' + self._functions['lambda_electrostatics'] - custom_bond_force = openmm.CustomBondForce("U_sterics + U_electrostatics;" + sterics_energy_expression + electrostatics_energy_expression) - #custom_bond_force.addGlobalParameter("lambda_electrostatics", 0.0) - #custom_bond_force.addGlobalParameter("lambda_sterics", 0.0) - custom_bond_force.addGlobalParameter("softcore_alpha", self.softcore_alpha) - custom_bond_force.addGlobalParameter("softcore_beta", self.softcore_beta) - custom_bond_force.addPerBondParameter("chargeprodA") - custom_bond_force.addPerBondParameter("sigmaA") - custom_bond_force.addPerBondParameter("epsilonA") - custom_bond_force.addPerBondParameter("chargeprodB") - custom_bond_force.addPerBondParameter("sigmaB") - custom_bond_force.addPerBondParameter("epsilonB") - - if self._has_functions: - custom_bond_force.addGlobalParameter('lambda', 0.0) - custom_bond_force.addEnergyParameterDerivative('lambda') - else: - custom_bond_force.addGlobalParameter("lambda_electrostatics", 0.0) - custom_bond_force.addGlobalParameter("lambda_sterics", 0.0) - - return custom_bond_force - def _find_bond_parameters(self, bond_force, index1, index2): """ This is a convenience function to find bond parameters in another system given the two indices. @@ -1259,11 +1184,11 @@ def handle_nonbonded(self): #add the particle to the hybrid custom sterics and electrostatics. self._hybrid_system_forces['core_sterics_force'].addParticle([sigma, epsilon, sigma, 0.0]) - self._hybrid_system_forces['core_electrostatics_force'].addParticle([charge, 0.0]) - #Add the particle to the regular nonbonded force as required, but zero out interaction - #it will be handled by an exception - self._hybrid_system_forces['standard_nonbonded_force'].addParticle(0.0, 1.0, 0.0) + # Add particle to the regular nonbonded force, but Lennard-Jones will be handled by CustomNonbondedForce + particle_index = self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, 0.0) + # Charge will be turned on at lambda_electrostatics = 0, off at lambda_electrostatics = 1 + self._hybrid_system_forces['standard_nonbonded_force'].addParticleParameterOffset('lambda_electrostatics', particle_index, -charge, 0, 0) elif particle_index in self._atom_classes['unique_new_atoms']: #get the parameters in the new system @@ -1272,12 +1197,11 @@ def handle_nonbonded(self): #add the particle to the hybrid custom sterics and electrostatics self._hybrid_system_forces['core_sterics_force'].addParticle([sigma, 0.0, sigma, epsilon]) - self._hybrid_system_forces['core_electrostatics_force'].addParticle([0.0, charge]) - - #Add the particle to the regular nonbonded force as required, but zero out interaction - #it will be handled by an exception - self._hybrid_system_forces['standard_nonbonded_force'].addParticle(0.0, 1.0, 0.0) + # Add particle to the regular nonbonded force, but Lennard-Jones will be handled by CustomNonbondedForce + particle_index = self._hybrid_system_forces['standard_nonbonded_force'].addParticle(0.0, sigma, 0.0) + # Charge will be turned off at lambda_electrostatics = 0, on at lambda_electrostatics = 1 + self._hybrid_system_forces['standard_nonbonded_force'].addParticleParameterOffset('lambda_electrostatics', particle_index, +charge, 0, 0) elif particle_index in self._atom_classes['core_atoms']: #get the parameters in the new and old systems: @@ -1288,10 +1212,12 @@ def handle_nonbonded(self): #add the particle to the custom forces, interpolating between the two parameters self._hybrid_system_forces['core_sterics_force'].addParticle([sigma_old, epsilon_old, sigma_new, epsilon_new]) - self._hybrid_system_forces['core_electrostatics_force'].addParticle([charge_old, charge_new]) #still add the particle to the regular nonbonded force, but with zeroed out parameters. - self._hybrid_system_forces['standard_nonbonded_force'].addParticle(0.0, 1.0, 0.0) + particle_index = self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge_old, 0.5*(sigma_old+sigma_new), 0.0) + # Charge is charge_old at lambda_electrostatics = 0, charge_new at lambda_electrostatics = 1 + # TODO: We could also interpolate the Lennard-Jones here instead of core_sterics force so that core_sterics_force could just be softcore + self._hybrid_system_forces['standard_nonbonded_force'].addParticleParameterOffset('lambda_electrostatics', particle_index, (charge_new - charge_old), 0, 0) #otherwise, the particle is in the environment else: @@ -1301,7 +1227,6 @@ def handle_nonbonded(self): #add the particle to the hybrid custom sterics and electrostatics, but they dont change self._hybrid_system_forces['core_sterics_force'].addParticle([sigma, epsilon, sigma, epsilon]) - self._hybrid_system_forces['core_electrostatics_force'].addParticle([charge, charge]) #add the environment atoms to the regular nonbonded force as well: self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, epsilon) @@ -1350,7 +1275,6 @@ def _handle_interaction_groups(self): Must be called after particles are added to the Nonbonded forces """ #get the force objects for convenience: - electrostatics_custom_force = self._hybrid_system_forces['core_electrostatics_force'] sterics_custom_force = self._hybrid_system_forces['core_sterics_force'] #also prepare the atom classes @@ -1360,32 +1284,24 @@ def _handle_interaction_groups(self): environment_atoms = self._atom_classes['environment_atoms'] - electrostatics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) sterics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) - electrostatics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) sterics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) - electrostatics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) sterics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) - electrostatics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) sterics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) - electrostatics_custom_force.addInteractionGroup(core_atoms, environment_atoms) sterics_custom_force.addInteractionGroup(core_atoms, environment_atoms) - electrostatics_custom_force.addInteractionGroup(core_atoms, core_atoms) sterics_custom_force.addInteractionGroup(core_atoms, core_atoms) def _handle_hybrid_exceptions(self): """ Instead of excluding interactions that shouldn't occur, we provide exceptions for interactions that were zeroed out but should occur. - Returns ------- - """ print("handling exceptions") @@ -1397,8 +1313,6 @@ def _handle_hybrid_exceptions(self): unique_old_atoms = self._atom_classes['unique_old_atoms'] unique_new_atoms = self._atom_classes['unique_new_atoms'] - nonbonded_force = self._hybrid_system_forces['standard_nonbonded_force'] - #get the list of interaction pairs for which we need to set exceptions: unique_old_pairs = list(itertools.combinations(unique_old_atoms, 2)) unique_new_pairs = list(itertools.combinations(unique_new_atoms, 2)) @@ -1411,21 +1325,27 @@ def _handle_hybrid_exceptions(self): #now we check if the pair is in the exception dictionary if old_index_atom_pair in self._old_system_exceptions: [chargeProd, sigma, epsilon] = self._old_system_exceptions[old_index_atom_pair] - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + chargeProd *= 0.0 + self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['core_sterics_force'].addExclusion(atom_pair[0], atom_pair[1]) # add exclusion to ensure exceptions are consistent #check if the pair is in the reverse order and use that if so elif old_index_atom_pair[::-1] in self._old_system_exceptions: [chargeProd, sigma, epsilon] = self._old_system_exceptions[old_index_atom_pair[::-1]] - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + chargeProd *= 0.0 + self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['core_sterics_force'].addExclusion(atom_pair[0], atom_pair[1]) # add exclusion to ensure exceptions are consistent #If it's not handled by an exception in the original system, we just add the regular parameters as an exception else: [charge0, sigma0, epsilon0] = self._old_system_forces['NonbondedForce'].getParticleParameters(old_index_atom_pair[0]) [charge1, sigma1, epsilon1] = self._old_system_forces['NonbondedForce'].getParticleParameters(old_index_atom_pair[1]) chargeProd = charge0*charge1 + chargeProd *= 0.0 epsilon = unit.sqrt(epsilon0*epsilon1) sigma = 0.5*(sigma0+sigma1) - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['core_sterics_force'].addExclusion(atom_pair[0], atom_pair[1]) # add exclusion to ensure exceptions are consistent #add back the interactions of the new unique atoms, unless there are exceptions for atom_pair in unique_new_pairs: @@ -1435,21 +1355,28 @@ def _handle_hybrid_exceptions(self): #now we check if the pair is in the exception dictionary if new_index_atom_pair in self._new_system_exceptions: [chargeProd, sigma, epsilon] = self._new_system_exceptions[new_index_atom_pair] - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + chargeProd *= 0.0 + self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['core_sterics_force'].addExclusion(atom_pair[0], atom_pair[1]) # add exclusion to ensure exceptions are consistent #check if the pair is present in the reverse order and use that if so elif new_index_atom_pair[::-1] in self._new_system_exceptions: [chargeProd, sigma, epsilon] = self._new_system_exceptions[new_index_atom_pair[::-1]] - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + chargeProd *= 0.0 + self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['core_sterics_force'].addExclusion(atom_pair[0], atom_pair[1]) # add exclusion to ensure exceptions are consistent #If it's not handled by an exception in the original system, we just add the regular parameters as an exception else: [charge0, sigma0, epsilon0] = self._new_system_forces['NonbondedForce'].getParticleParameters(new_index_atom_pair[0]) [charge1, sigma1, epsilon1] = self._new_system_forces['NonbondedForce'].getParticleParameters(new_index_atom_pair[1]) chargeProd = charge0*charge1 + chargeProd *= 0.0 epsilon = unit.sqrt(epsilon0*epsilon1) sigma = 0.5*(sigma0+sigma1) - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) + self._hybrid_system_forces['core_sterics_force'].addExclusion(atom_pair[0], atom_pair[1]) # add exclusion to ensure exceptions are consistent + print("done handling exceptions") def _handle_original_exceptions(self): @@ -1478,20 +1405,16 @@ def _handle_original_exceptions(self): #in the unique-old case, it is handled elsewhere due to internal peculiarities regarding exceptions if index_set.issubset(self._atom_classes['environment_atoms']): self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_old, sigma_old, epsilon_old) + self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #otherwise, check if one of the atoms in the set is in the unique_old_group: - elif len(index_set.intersection(self._atom_classes['unique_old_atoms'])) > 0: - #if it is, we should add it to the CustomBondForce for the nonbonded exceptions, and have it remain on - #by having the two endpoints with the same parameters. - #Currently, we keep sigma at the same value - self._hybrid_system_forces['core_nonbonded_bond_force'].addBond(index1_hybrid, index2_hybrid, - [chargeProd_old, sigma_old, - epsilon_old, chargeProd_old, - sigma_old, epsilon_old]) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting + #we have already handled unique old - unique old exceptions + elif len(index_set.intersection(self._atom_classes['unique_old_atoms'])) == 2: + continue + + #otherwise, check if one of the atoms in the set is in the unique_old_group and the other is not: + elif len(index_set.intersection(self._atom_classes['unique_old_atoms'])) == 1: + self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_old, sigma_old, epsilon_old) self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) #If the exception particles are neither solely old unique, solely environment, nor contain any unique old atoms, they are either core/environment or core/core #In this case, we need to get the parameters from the exception in the other (new) system, and interpolate between the two @@ -1504,15 +1427,11 @@ def _handle_original_exceptions(self): [index1_new, index2_new, chargeProd_new, sigma_new, epsilon_new] = self._find_exception( new_system_nonbonded_force, index1_new, index2_new) - #Now add a term to the CustomBondForce to interpolate between the new and old systems: - self._hybrid_system_forces['core_nonbonded_bond_force'].addBond(index1_hybrid, index2_hybrid, - [chargeProd_old, sigma_old, - epsilon_old, chargeProd_new, - sigma_new, epsilon_new]) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting + #interpolate between old and new + exception_index = self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_old, sigma_old, epsilon_old) + self._hybrid_system_forces['standard_nonbonded_force'].addExceptionParameterOffset('lambda_electrostatics', exception_index, (chargeProd_new - chargeProd_old), 0, 0) + self._hybrid_system_forces['standard_nonbonded_force'].addExceptionParameterOffset('lambda_sterics', exception_index, 0, (sigma_new - sigma_old), (epsilon_new - epsilon_old)) self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) #now, loop through the new system to collect remaining interactions. The only that remain here are #uniquenew-uniquenew, uniquenew-core, and uniquenew-environment. @@ -1533,14 +1452,8 @@ def _handle_original_exceptions(self): #look for the final class- interactions between uniquenew-core and uniquenew-environment. They are treated #similarly: they are simply on and constant the entire time (as a valence term) elif len(index_set.intersection(self._atom_classes['unique_new_atoms'])) > 0: - self._hybrid_system_forces['core_nonbonded_bond_force'].addBond(index1_hybrid, index2_hybrid, - [chargeProd_new, sigma_new, - epsilon_new, chargeProd_new, - sigma_new, epsilon_new]) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting + self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_new, sigma_new, epsilon_new) self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) def _find_exception(self, force, index1, index2): """ diff --git a/perses/dispersed/relative_setup.py b/perses/dispersed/relative_setup.py index 773aa926a..5a9eb8f36 100644 --- a/perses/dispersed/relative_setup.py +++ b/perses/dispersed/relative_setup.py @@ -187,7 +187,8 @@ def __init__(self, ligand_file, old_ligand_index, new_ligand_index, forcefield_f barostat = None self._system_generator = SystemGenerator(forcefield_files, barostat=barostat, forcefield_kwargs={'nonbondedMethod': self._nonbonded_method, - 'constraints': app.HBonds}) + 'constraints': app.HBonds, + 'hydrogenMass': 4 * unit.amus}) else: self._system_generator = SystemGenerator(forcefield_files, forcefield_kwargs={'constraints': app.HBonds}) @@ -878,11 +879,11 @@ def __init__(self, *args, hybrid_factory=None, **kwargs): super(HybridSAMSSampler, self).__init__(*args, hybrid_factory=hybrid_factory, **kwargs) self._factory = hybrid_factory - def setup(self, n_states, temperature, storage_file, checkpoint_interval): + def setup(self, n_states, temperature, storage_file): hybrid_system = self._factory.hybrid_system initial_hybrid_positions = self._factory.hybrid_positions lambda_zero_alchemical_state = alchemy.AlchemicalState.from_system(hybrid_system) - lambda_zero_alchemical_state.set_alchemical_parameters(1.0) + #lambda_zero_alchemical_state.set_alchemical_parameters(1.0) thermostate = states.ThermodynamicState(hybrid_system, temperature=temperature) compound_thermodynamic_state = states.CompoundThermodynamicState(thermostate, composable_states=[lambda_zero_alchemical_state]) @@ -1050,7 +1051,7 @@ def run_setup(setup_options): ne_fep[phase] = NonequilibriumSwitchingFEP(top_prop['%s_topology_proposal' % phase], top_prop['%s_old_positions' % phase], top_prop['%s_new_positions' % phase], - forward_functions=forward_functions, + forward_functions=forward_functions, n_equil_steps=n_equilibrium_steps_per_iteration, ncmc_nsteps=n_steps_ncmc_protocol, nsteps_per_iteration=n_steps_per_move_application, @@ -1074,7 +1075,7 @@ def run_setup(setup_options): for phase in phases: htf[phase] = HybridTopologyFactory(top_prop['%s_topology_proposal' % phase], top_prop['%s_old_positions' % phase], - top_prop['%s_new_positions' % phase]) + top_prop['%s_new_positions' % phase], softcore_method='amber') if atom_selection: selection_indices = htf[phase].hybrid_topology.select(atom_selection) @@ -1082,16 +1083,17 @@ def run_setup(setup_options): selection_indices = None storage_name = "-".join([trajectory_prefix, '%s.nc' % phase]) - reporter = MultiStateReporter(storage_name, analysis_particle_indices=selection_indices) - - hss[phase] = HybridSAMSSampler(mcmc_moves=mcmc.LangevinDynamicsMove(timestep=2.0 * unit.femtosecond, - collision_rate=5.0 / unit.picosecond, - n_steps=n_steps_per_move_application, - reassign_velocities=False, - n_restart_attempts=6), - hybrid_factory=htf[phase]) - hss[phase].setup(n_states=n_states, temperature=300.0 * unit.kelvin, - storage_file=reporter, - checkpoint_interval=10) + reporter = MultiStateReporter(storage_name, analysis_particle_indices=selection_indices, + checkpoint_interval=10) + + hss[phase] = HybridSAMSSampler(mcmc_moves=mcmc.LangevinSplittingDynamicsMove(timestep=4.0 * unit.femtosecond, + collision_rate=5.0 / unit.picosecond, + n_steps=n_steps_per_move_application, + reassign_velocities=False, + n_restart_attempts=6, + splitting="V R R R O R R R V"), + hybrid_factory=htf[phase], online_analysis_interval=10, + online_analysis_target_error=0.2, online_analysis_minimum_iterations=10) + hss[phase].setup(n_states=n_states, temperature=300.0 * unit.kelvin, storage_file=reporter) return {'topology_proposals': top_prop, 'hybrid_topology_factories': htf, 'hybrid_sams_samplers': hss} diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 05a891624..7bc8ebba2 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -373,6 +373,7 @@ def _oemol_from_residue(res, verbose=True): oemol : openeye.oechem.OEMol an oemol representation of the residue with topology indices """ + # TODO: This seems to be broken. Can we fix it? from openmoltools.forcefield_generators import generateOEMolFromTopologyResidue external_bonds = list(res.external_bonds()) for bond in external_bonds: diff --git a/perses/rjmc/topology_proposal.py b/perses/rjmc/topology_proposal.py index 931661b31..f73236d31 100644 --- a/perses/rjmc/topology_proposal.py +++ b/perses/rjmc/topology_proposal.py @@ -1324,7 +1324,7 @@ def __init__(self, list_of_smiles, system_generator, residue_name='MOL', super(SmallMoleculeSetProposalEngine, self).__init__(system_generator, proposal_metadata=proposal_metadata, always_change=always_change) - def propose(self, current_system, current_topology, current_metadata=None): + def propose(self, current_system, current_topology, current_mol=None, proposed_mol=None, current_metadata=None): """ Propose the next state, given the current state @@ -1336,6 +1336,10 @@ def propose(self, current_system, current_topology, current_metadata=None): the topology of the current state current_metadata : dict dict containing current smiles as a key + current_mol : OEMol, optional, default=None + If specified, use this OEMol instead of converting from topology + proposed_mol : OEMol, optional, default=None + If specified, use this OEMol instead of converting from topology Returns ------- @@ -1343,7 +1347,11 @@ def propose(self, current_system, current_topology, current_metadata=None): topology proposal object """ # Determine SMILES string for current small molecule - current_mol_smiles, current_mol = self._topology_to_smiles(current_topology) + if current_mol is None: + current_mol_smiles, current_mol = self._topology_to_smiles(current_topology) + else: + # TODO: Make sure we're using canonical mol to smiles conversion + current_mol_smiles = oechem.OEMolToSmiles(current_mol) # Remove the small molecule from the current Topology object current_receptor_topology = self._remove_small_molecule(current_topology) @@ -1355,7 +1363,12 @@ def propose(self, current_system, current_topology, current_metadata=None): old_alchemical_atoms = range(old_mol_start_index, len_old_mol) # Select the next molecule SMILES given proposal probabilities - proposed_mol_smiles, proposed_mol, logp_proposal = self._propose_molecule(current_system, current_topology, current_mol_smiles) + if proposed_mol is None: + proposed_mol_smiles, proposed_mol, logp_proposal = self._propose_molecule(current_system, current_topology, current_mol_smiles) + else: + # TODO: Make sure we're using canonical mol to smiles conversion + proposed_mol_smiles = oechem.OEMolToSmiles(current_mol) + logp_proposal = 0.0 # Build the new Topology object, including the proposed molecule new_topology = self._build_new_topology(current_receptor_topology, proposed_mol) diff --git a/perses/tests/test_hybrid_builder.py b/perses/tests/test_hybrid_builder.py index 1ce3b17a8..1292d8ccb 100644 --- a/perses/tests/test_hybrid_builder.py +++ b/perses/tests/test_hybrid_builder.py @@ -2,6 +2,12 @@ from simtk import unit, openmm import numpy as np import os + +try: + from StringIO import StringIO +except ImportError: + from io import StringIO + from perses.annihilation.new_relative import HybridTopologyFactory from perses.rjmc.geometry import FFAllAngleGeometryEngine from perses.rjmc.topology_proposal import SmallMoleculeSetProposalEngine, SystemGenerator, TopologyProposal @@ -123,18 +129,18 @@ forcefield = app.ForceField('amber99sbildn.xml') -def generate_vacuum_topology_proposal(mol_name="naphthalene", ref_mol_name="benzene"): +def generate_vacuum_topology_proposal(current_mol_name="benzene", proposed_mol_name="toluene"): """ Generate a test vacuum topology proposal, current positions, and new positions triplet from two IUPAC molecule names. Parameters ---------- - mol_name : str, optional + current_mol_name : str, optional name of the first molecule - ref_mol_name : str, optional + proposed_mol_name : str, optional name of the second molecule - + Returns ------- topology_proposal : perses.rjmc.topology_proposal @@ -148,11 +154,11 @@ def generate_vacuum_topology_proposal(mol_name="naphthalene", ref_mol_name="benz from perses.tests.utils import createOEMolFromIUPAC, createSystemFromIUPAC, get_data_filename - m, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(mol_name) - refmol = createOEMolFromIUPAC(ref_mol_name) + current_mol, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(current_mol_name) + proposed_mol = createOEMolFromIUPAC(proposed_mol_name) - initial_smiles = oechem.OEMolToSmiles(m) - final_smiles = oechem.OEMolToSmiles(refmol) + initial_smiles = oechem.OEMolToSmiles(current_mol) + final_smiles = oechem.OEMolToSmiles(proposed_mol) gaff_xml_filename = get_data_filename("data/gaff.xml") forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml') @@ -161,31 +167,31 @@ def generate_vacuum_topology_proposal(mol_name="naphthalene", ref_mol_name="benz solvated_system = forcefield.createSystem(top_old, removeCMMotion=False) gaff_filename = get_data_filename('data/gaff.xml') - system_generator = SystemGenerator([gaff_filename, 'amber99sbildn.xml', 'tip3p.xml']) + system_generator = SystemGenerator([gaff_filename, 'amber99sbildn.xml', 'tip3p.xml'], forcefield_kwargs={'removeCMMotion': False, 'nonbondedMethod': app.NoCutoff}) geometry_engine = FFAllAngleGeometryEngine() proposal_engine = SmallMoleculeSetProposalEngine( - [initial_smiles, final_smiles], system_generator, residue_name=mol_name) + [initial_smiles, final_smiles], system_generator, residue_name=current_mol_name) #generate topology proposal - topology_proposal = proposal_engine.propose(solvated_system, top_old) + topology_proposal = proposal_engine.propose(solvated_system, top_old, current_mol=current_mol, proposed_mol=proposed_mol) #generate new positions with geometry engine new_positions, _ = geometry_engine.propose(topology_proposal, pos_old, beta) return topology_proposal, pos_old, new_positions -def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name="benzene"): +def generate_solvated_hybrid_test_topology(current_mol_name="naphthalene", proposed_mol_name="benzene"): """ Generate a test solvated topology proposal, current positions, and new positions triplet from two IUPAC molecule names. Parameters ---------- - mol_name : str, optional + current_mol_name : str, optional name of the first molecule - ref_mol_name : str, optional + proposed_mol_name : str, optional name of the second molecule - + Returns ------- topology_proposal : perses.rjmc.topology_proposal @@ -200,11 +206,11 @@ def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name= from perses.tests.utils import createOEMolFromIUPAC, createSystemFromIUPAC, get_data_filename - m, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(mol_name) - refmol = createOEMolFromIUPAC(ref_mol_name) + current_mol, unsolv_old_system, pos_old, top_old = createSystemFromIUPAC(current_mol_name) + proposed_mol = createOEMolFromIUPAC(proposed_mol_name) - initial_smiles = oechem.OEMolToSmiles(m) - final_smiles = oechem.OEMolToSmiles(refmol) + initial_smiles = oechem.OEMolToSmiles(current_mol) + final_smiles = oechem.OEMolToSmiles(proposed_mol) gaff_xml_filename = get_data_filename("data/gaff.xml") forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml') @@ -219,13 +225,12 @@ def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name= solvated_system.addForce(barostat) - gaff_filename = get_data_filename('data/gaff.xml') system_generator = SystemGenerator([gaff_filename, 'amber99sbildn.xml', 'tip3p.xml'], barostat=barostat, forcefield_kwargs={'removeCMMotion': False, 'nonbondedMethod': app.PME}) geometry_engine = FFAllAngleGeometryEngine() proposal_engine = SmallMoleculeSetProposalEngine( - [initial_smiles, final_smiles], system_generator, residue_name=mol_name) + [initial_smiles, final_smiles], system_generator, residue_name=current_mol_name) #generate topology proposal topology_proposal = proposal_engine.propose(solvated_system, solvated_topology) @@ -306,7 +311,7 @@ def check_result(results, threshold=3.0, neffmin=10): def test_simple_overlap(): """Test that the variance of the endpoint->nonalchemical perturbation is sufficiently small for pentane->butane in vacuum""" - topology_proposal, current_positions, new_positions = generate_vacuum_topology_proposal() + topology_proposal, current_positions, new_positions = generate_vacuum_topology_proposal(current_mol_name='imatinib', proposed_mol_name='nilotinib') results = run_hybrid_endpoint_overlap(topology_proposal, current_positions, new_positions) for idx, lambda_result in enumerate(results): @@ -320,7 +325,23 @@ def test_simple_overlap(): @skipIf(istravis, "Skip expensive test on travis") def test_difficult_overlap(): """Test that the variance of the endpoint->nonalchemical perturbation is sufficiently small for imatinib->nilotinib in solvent""" - topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(mol_name='imatinib', ref_mol_name='nilotinib') + name1 = 'imatinib' + name2 = 'nilotinib' + + print(name1, name2) + topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(current_mol_name=name1, proposed_mol_name=name2) + results = run_hybrid_endpoint_overlap(topology_proposal, solvated_positions, new_positions) + + for idx, lambda_result in enumerate(results): + try: + check_result(lambda_result) + except Exception as e: + message = "solvated imatinib->nilotinib failed at lambda %d \n" % idx + message += str(e) + raise Exception(message) + + print(name2, name1) + topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(current_mol_name=name2, proposed_mol_name=name1) results = run_hybrid_endpoint_overlap(topology_proposal, solvated_positions, new_positions) for idx, lambda_result in enumerate(results): diff --git a/scripts/setup_relative_calculation.py b/scripts/setup_relative_calculation.py index 9c31f88f2..117e2248b 100644 --- a/scripts/setup_relative_calculation.py +++ b/scripts/setup_relative_calculation.py @@ -32,8 +32,9 @@ setup_dict = relative_setup.run_setup(setup_options) print("setup complete") + print('setup_dict keys: {}'.format(setup_dict.keys())) - n_equilibration_iterations = setup_dict['n_equilibration_iterations'] + n_equilibration_iterations = setup_options['n_equilibration_iterations'] trajectory_prefix = setup_options['trajectory_prefix'] #write out topology proposals @@ -97,13 +98,17 @@ setup_dict['hybrid_topology_factories']) hss = setup_dict['hybrid_sams_samplers'] + logZ = dict() free_energies = dict() for phase in phases: hss_run = hss[phase] hss_run.minimize() hss_run.equilibrate(n_equilibration_iterations) - hss_run.extend(1000) - free_energies[phase] = hss_run._logZ[-1] - hss_run._logZ[0] - print("Finished phase %s with dG estimated as %.4f kT" % (phase, free_energies[phase])) + hss_run.extend(setup_options['n_cycles']) + logZ[phase] = hss_run._logZ[-1] - hss_run._logZ[0] + free_energies[phase] = hss_run._last_mbar_f_k[-1] - hss_run._last_mbar_f_k[0] + print("Finished phase %s with dG estimated as %.4f +/- %.4f kT" % ( + phase, free_energies[phase], hss_run._last_err_free_energy)) + print("Finished phase %s with logZ dG estimated as %.4f kT" % (phase, logZ[phase])) print("Total ddG is estimated as %.4f kT" % (free_energies['complex'] - free_energies['solvent'])) \ No newline at end of file