From db8534aa1d2205df850485ed9b129ad11274c821 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 00:50:01 -0400 Subject: [PATCH 01/18] Modifications to use new 'offset' OpenMM feature --- perses/annihilation/new_relative.py | 139 +++++++++++++++------------- 1 file changed, 76 insertions(+), 63 deletions(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index c7954ac94..cac356078 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -84,17 +84,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,31 +700,32 @@ def _add_nonbonded_force_terms(self): custom_nonbonded_method = self._translate_nonbonded_method_to_custom(self._nonbonded_method) + # TODO: Delete # 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 + #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 + #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) + #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) + #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 + #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: @@ -767,17 +768,18 @@ 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) + #electrostatics_custom_nonbonded_force.setUseSwitchingFunction(True) # TODO: Delete + #electrostatics_custom_nonbonded_force.setSwitchingDistance(switching_distance) # TODO: Delete else: standard_nonbonded_force.setUseSwitchingFunction(False) - electrostatics_custom_nonbonded_force.setUseSwitchingFunction(False) + #electrostatics_custom_nonbonded_force.setUseSwitchingFunction(False) # TODO: Delete sterics_custom_nonbonded_force.setUseSwitchingFunction(False) + # TODO: Delete #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 + #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): """ @@ -1259,11 +1261,12 @@ 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]) + #self._hybrid_system_forces['core_electrostatics_force'].addParticle([charge, 0.0]) # TODO: Remove once core_electrostatics_force is eliminated - #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 +1275,12 @@ 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) + #self._hybrid_system_forces['core_electrostatics_force'].addParticle([0.0, charge]) # TODO: Remove once core_electrostatics_force is eliminated + # 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 +1291,13 @@ 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]) + #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,13 +1307,13 @@ 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]) + #self._hybrid_system_forces['core_electrostatics_force'].addParticle([charge, charge]) # TODO: Remove once core_electrostatics_force is eliminated #add the environment atoms to the regular nonbonded force as well: self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, epsilon) self._handle_interaction_groups() - self._handle_hybrid_exceptions() + #self._handle_hybrid_exceptions() # TODO: Delete? self._handle_original_exceptions() def _generate_dict_from_exceptions(self, force): @@ -1350,7 +1356,7 @@ 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'] + #electrostatics_custom_force = self._hybrid_system_forces['core_electrostatics_force'] # TODO: Delete sterics_custom_force = self._hybrid_system_forces['core_sterics_force'] #also prepare the atom classes @@ -1360,22 +1366,22 @@ def _handle_interaction_groups(self): environment_atoms = self._atom_classes['environment_atoms'] - electrostatics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) + #electrostatics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) - electrostatics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) + #electrostatics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) - electrostatics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) + #electrostatics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) - electrostatics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) + #electrostatics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) - electrostatics_custom_force.addInteractionGroup(core_atoms, environment_atoms) + #electrostatics_custom_force.addInteractionGroup(core_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(core_atoms, environment_atoms) - electrostatics_custom_force.addInteractionGroup(core_atoms, core_atoms) + #electrostatics_custom_force.addInteractionGroup(core_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(core_atoms, core_atoms) def _handle_hybrid_exceptions(self): @@ -1484,14 +1490,16 @@ def _handle_original_exceptions(self): #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]) + #self._hybrid_system_forces['core_nonbonded_bond_force'].addBond(index1_hybrid, index2_hybrid, # TODO: Delete this when we remove core_nonbonded_bond_force + # [chargeProd_old, sigma_old, + # epsilon_old, chargeProd_old, + # sigma_old, epsilon_old]) + + self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_old, sigma_old, epsilon_old) #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) + #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated #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 @@ -1505,14 +1513,17 @@ def _handle_original_exceptions(self): 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]) + #self._hybrid_system_forces['core_nonbonded_bond_force'].addBond(index1_hybrid, index2_hybrid, # TODO: Delete this when we remove core_nonbonded_bond_force + # [chargeProd_old, sigma_old, + # epsilon_old, chargeProd_new, + # sigma_new, epsilon_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)) #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) + #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated #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 +1544,16 @@ 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]) + #self._hybrid_system_forces['core_nonbonded_bond_force'].addBond(index1_hybrid, index2_hybrid, # TODO: Delete this when we remove core_nonbonded_bond_force + # [chargeProd_new, sigma_new, + # epsilon_new, chargeProd_new, + # sigma_new, epsilon_new]) + self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, 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['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) + #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated def _find_exception(self, force, index1, index2): """ From 1e1cae9b97916aaadead3f3e0d3fbfcf7b43a1e0 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 00:54:09 -0400 Subject: [PATCH 02/18] Increment openmm min version to 7.3.0 --- devtools/conda-recipe/meta.yaml | 4 ++-- perses/rjmc/geometry.py | 1 + perses/rjmc/topology_proposal.py | 19 ++++++++++++++--- perses/tests/test_hybrid_builder.py | 33 ++++++++++++++++++----------- 4 files changed, 40 insertions(+), 17 deletions(-) diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 43b816e9b..9e9d86d7f 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -18,7 +18,7 @@ requirements: - numexpr - autograd - pymbar - - openmm >=7.2.0 + - openmm >=7.3.0 - parmed - openmoltools - alchemy >=1.2.3 @@ -43,7 +43,7 @@ requirements: - numexpr - autograd - pymbar - - openmm >=7.2.0 + - openmm >=7.3.0 - parmed - openmoltools - alchemy >=1.2.3 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..1e6f1301e 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,26 +154,29 @@ 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') - forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) + ffxml = forcefield_generators.generateForceFieldFromMolecules([current_mol, proposed_mol]) + forcefield.loadFile(StringIO(ffxml)) + #forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) 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._forcefield.loadFile(StringIO(ffxml)) 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) @@ -185,7 +194,7 @@ def generate_solvated_hybrid_test_topology(mol_name="naphthalene", ref_mol_name= name of the first molecule ref_mol_name : str, optional name of the second molecule - + Returns ------- topology_proposal : perses.rjmc.topology_proposal From 0c9c22449e9bd4f8f476db4355884b9d7fe8b286 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 21:09:57 -0400 Subject: [PATCH 03/18] Refine exceptions --- perses/annihilation/new_relative.py | 86 +---------------------------- 1 file changed, 1 insertion(+), 85 deletions(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index cac356078..b1d020ba0 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -1313,7 +1313,6 @@ def handle_nonbonded(self): self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, epsilon) self._handle_interaction_groups() - #self._handle_hybrid_exceptions() # TODO: Delete? self._handle_original_exceptions() def _generate_dict_from_exceptions(self, force): @@ -1384,80 +1383,6 @@ def _handle_interaction_groups(self): #electrostatics_custom_force.addInteractionGroup(core_atoms, core_atoms) # TODO: Delete 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") - - old_system_nonbonded_force = self._old_system_forces['NonbondedForce'] - new_system_nonbonded_force = self._new_system_forces['NonbondedForce'] - - import itertools - #prepare the atom classes - 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)) - - #add back the interactions of the old unique atoms, unless there are exceptions - for atom_pair in unique_old_pairs: - #since the pairs are indexed in the dictionary by the old system indices, we need to convert - old_index_atom_pair = (self._hybrid_to_old_map[atom_pair[0]], self._hybrid_to_old_map[atom_pair[1]]) - - #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) - - #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) - - #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 - epsilon = unit.sqrt(epsilon0*epsilon1) - sigma = 0.5*(sigma0+sigma1) - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) - - #add back the interactions of the new unique atoms, unless there are exceptions - for atom_pair in unique_new_pairs: - #since the pairs are indexed in the dictionary by the new system indices, we need to convert - new_index_atom_pair = (self._hybrid_to_new_map[atom_pair[0]], self._hybrid_to_new_map[atom_pair[1]]) - - #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) - - #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) - - #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 - epsilon = unit.sqrt(epsilon0*epsilon1) - sigma = 0.5*(sigma0+sigma1) - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) - print("done handling exceptions") - def _handle_original_exceptions(self): """ This method ensures that exceptions present in the original systems are present in the hybrid appropriately. @@ -1484,6 +1409,7 @@ 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: @@ -1496,10 +1422,7 @@ def _handle_original_exceptions(self): # sigma_old, epsilon_old]) self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_old, sigma_old, epsilon_old) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated #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 @@ -1520,10 +1443,7 @@ def _handle_original_exceptions(self): 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)) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated #now, loop through the new system to collect remaining interactions. The only that remain here are #uniquenew-uniquenew, uniquenew-core, and uniquenew-environment. @@ -1549,11 +1469,7 @@ def _handle_original_exceptions(self): # epsilon_new, chargeProd_new, # sigma_new, epsilon_new]) self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, 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['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated def _find_exception(self, force, index1, index2): """ From 94445269ce96c7930cd8a6e06e3eb0e715fe38b2 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 21:17:29 -0400 Subject: [PATCH 04/18] Cleanup --- perses/annihilation/new_relative.py | 67 +---------------------------- 1 file changed, 1 insertion(+), 66 deletions(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index b1d020ba0..b33e6baf7 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -700,33 +700,6 @@ def _add_nonbonded_force_terms(self): custom_nonbonded_method = self._translate_nonbonded_method_to_custom(self._nonbonded_method) - # TODO: Delete - # 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: @@ -768,19 +741,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) # TODO: Delete - #electrostatics_custom_nonbonded_force.setSwitchingDistance(switching_distance) # TODO: Delete else: standard_nonbonded_force.setUseSwitchingFunction(False) - #electrostatics_custom_nonbonded_force.setUseSwitchingFunction(False) # TODO: Delete sterics_custom_nonbonded_force.setUseSwitchingFunction(False) - # TODO: Delete - #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 @@ -961,8 +925,6 @@ def _nonbonded_custom_bond_force(self, sterics_energy_expression, electrostatics 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") @@ -1261,7 +1223,6 @@ 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]) # TODO: Remove once core_electrostatics_force is eliminated # 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) @@ -1275,7 +1236,6 @@ 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]) # TODO: Remove once core_electrostatics_force is eliminated # 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) @@ -1291,7 +1251,6 @@ 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. particle_index = self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge_old, 0.5*(sigma_old+sigma_new), 0.0) @@ -1307,7 +1266,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]) # TODO: Remove once core_electrostatics_force is eliminated #add the environment atoms to the regular nonbonded force as well: self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, epsilon) @@ -1355,7 +1313,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'] # TODO: Delete sterics_custom_force = self._hybrid_system_forces['core_sterics_force'] #also prepare the atom classes @@ -1365,22 +1322,16 @@ def _handle_interaction_groups(self): environment_atoms = self._atom_classes['environment_atoms'] - #electrostatics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) - #electrostatics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) - #electrostatics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) - #electrostatics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) - #electrostatics_custom_force.addInteractionGroup(core_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(core_atoms, environment_atoms) - #electrostatics_custom_force.addInteractionGroup(core_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(core_atoms, core_atoms) def _handle_original_exceptions(self): @@ -1413,14 +1364,6 @@ def _handle_original_exceptions(self): #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, # TODO: Delete this when we remove core_nonbonded_bond_force - # [chargeProd_old, sigma_old, - # epsilon_old, chargeProd_old, - # sigma_old, epsilon_old]) - 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) @@ -1435,11 +1378,7 @@ 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, # TODO: Delete this when we remove core_nonbonded_bond_force - # [chargeProd_old, sigma_old, - # epsilon_old, chargeProd_new, - # sigma_new, epsilon_new]) + #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)) @@ -1464,10 +1403,6 @@ 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, # TODO: Delete this when we remove core_nonbonded_bond_force - # [chargeProd_new, sigma_new, - # epsilon_new, chargeProd_new, - # sigma_new, epsilon_new]) 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) From c7b3deeaa790907e74b4b578d742f50b57065802 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 21:56:40 -0400 Subject: [PATCH 05/18] Clean up vacuum and solvent tests --- perses/tests/test_hybrid_builder.py | 32 +++++++++++++++-------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/perses/tests/test_hybrid_builder.py b/perses/tests/test_hybrid_builder.py index 1e6f1301e..785ff06bb 100644 --- a/perses/tests/test_hybrid_builder.py +++ b/perses/tests/test_hybrid_builder.py @@ -162,14 +162,14 @@ def generate_vacuum_topology_proposal(current_mol_name="benzene", proposed_mol_n gaff_xml_filename = get_data_filename("data/gaff.xml") forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml') - ffxml = forcefield_generators.generateForceFieldFromMolecules([current_mol, proposed_mol]) - forcefield.loadFile(StringIO(ffxml)) - #forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) + #ffxml = forcefield_generators.generateForceFieldFromMolecules([current_mol, proposed_mol]) + #forcefield.loadFile(StringIO(ffxml)) + forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) 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}) system_generator._forcefield.loadFile(StringIO(ffxml)) geometry_engine = FFAllAngleGeometryEngine() proposal_engine = SmallMoleculeSetProposalEngine( @@ -183,16 +183,16 @@ def generate_vacuum_topology_proposal(current_mol_name="benzene", proposed_mol_n 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 @@ -209,14 +209,16 @@ 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') + #ffxml = forcefield_generators.generateForceFieldFromMolecules([current_mol, proposed_mol]) + #forcefield.loadFile(StringIO(ffxml)) forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) modeller = app.Modeller(top_old, pos_old) @@ -228,13 +230,13 @@ 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}) + #system_generator._forcefield.loadFile(StringIO(ffxml)) 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) @@ -315,7 +317,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): @@ -329,7 +331,7 @@ 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') + topology_proposal, solvated_positions, new_positions = generate_solvated_hybrid_test_topology(current_mol_name='imatinib', proposed_mol_name='nilotinib') results = run_hybrid_endpoint_overlap(topology_proposal, solvated_positions, new_positions) for idx, lambda_result in enumerate(results): From 1e4626c2ea42c6d295ecb3646ee24e7630d23835 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 22:18:35 -0400 Subject: [PATCH 06/18] Pin openmm to 7.3.0 for testing --- devtools/conda-recipe/meta.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 9e9d86d7f..9dc2a3c5d 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -18,7 +18,7 @@ requirements: - numexpr - autograd - pymbar - - openmm >=7.3.0 + - openmm ==7.3.0 - parmed - openmoltools - alchemy >=1.2.3 @@ -43,7 +43,7 @@ requirements: - numexpr - autograd - pymbar - - openmm >=7.3.0 + - openmm ==7.3.0 - parmed - openmoltools - alchemy >=1.2.3 From 29cebab613c483147fa2c10cf024c05269f64fe4 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 22:18:53 -0400 Subject: [PATCH 07/18] Fix bug in setup_relative_calculation script where setup_dict was erroneously used instead of setup_options --- scripts/setup_relative_calculation.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/scripts/setup_relative_calculation.py b/scripts/setup_relative_calculation.py index 9c31f88f2..4d527b852 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 @@ -106,4 +107,4 @@ 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])) - print("Total ddG is estimated as %.4f kT" % (free_energies['complex'] - free_energies['solvent'])) \ No newline at end of file + print("Total ddG is estimated as %.4f kT" % (free_energies['complex'] - free_energies['solvent'])) From 3c42fdc7d80a519026f72bd13df51302ab92de68 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Fri, 29 Jun 2018 23:04:47 -0400 Subject: [PATCH 08/18] Fix some test failures --- perses/dispersed/relative_setup.py | 2 +- perses/tests/test_hybrid_builder.py | 6 ------ 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/perses/dispersed/relative_setup.py b/perses/dispersed/relative_setup.py index 773aa926a..20da104fd 100644 --- a/perses/dispersed/relative_setup.py +++ b/perses/dispersed/relative_setup.py @@ -1050,7 +1050,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, diff --git a/perses/tests/test_hybrid_builder.py b/perses/tests/test_hybrid_builder.py index 785ff06bb..18e844cb9 100644 --- a/perses/tests/test_hybrid_builder.py +++ b/perses/tests/test_hybrid_builder.py @@ -162,15 +162,12 @@ def generate_vacuum_topology_proposal(current_mol_name="benzene", proposed_mol_n gaff_xml_filename = get_data_filename("data/gaff.xml") forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml') - #ffxml = forcefield_generators.generateForceFieldFromMolecules([current_mol, proposed_mol]) - #forcefield.loadFile(StringIO(ffxml)) forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) 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'], forcefield_kwargs={'removeCMMotion': False, 'nonbondedMethod': app.NoCutoff}) - system_generator._forcefield.loadFile(StringIO(ffxml)) geometry_engine = FFAllAngleGeometryEngine() proposal_engine = SmallMoleculeSetProposalEngine( [initial_smiles, final_smiles], system_generator, residue_name=current_mol_name) @@ -217,8 +214,6 @@ def generate_solvated_hybrid_test_topology(current_mol_name="naphthalene", propo gaff_xml_filename = get_data_filename("data/gaff.xml") forcefield = app.ForceField(gaff_xml_filename, 'tip3p.xml') - #ffxml = forcefield_generators.generateForceFieldFromMolecules([current_mol, proposed_mol]) - #forcefield.loadFile(StringIO(ffxml)) forcefield.registerTemplateGenerator(forcefield_generators.gaffTemplateGenerator) modeller = app.Modeller(top_old, pos_old) @@ -233,7 +228,6 @@ def generate_solvated_hybrid_test_topology(current_mol_name="naphthalene", propo 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}) - #system_generator._forcefield.loadFile(StringIO(ffxml)) geometry_engine = FFAllAngleGeometryEngine() proposal_engine = SmallMoleculeSetProposalEngine( [initial_smiles, final_smiles], system_generator, residue_name=current_mol_name) From 84ada3e3224317d63e4e50c53a351ac1bc3c8ddb Mon Sep 17 00:00:00 2001 From: brycestx Date: Sat, 30 Jun 2018 14:47:30 -0400 Subject: [PATCH 09/18] Fixes bug in relative_setup where first state lambda is set to 0. Also adds more support for sams as default method and online analysis. --- examples/cdk2-example/cdk2_setup.yaml | 17 +-- perses/annihilation/new_relative.py | 153 +------------------------- perses/dispersed/relative_setup.py | 32 +++--- scripts/setup_relative_calculation.py | 20 ++-- 4 files changed, 42 insertions(+), 180 deletions(-) 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 cac356078..b33e6baf7 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -700,33 +700,6 @@ def _add_nonbonded_force_terms(self): custom_nonbonded_method = self._translate_nonbonded_method_to_custom(self._nonbonded_method) - # TODO: Delete - # 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: @@ -768,19 +741,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) # TODO: Delete - #electrostatics_custom_nonbonded_force.setSwitchingDistance(switching_distance) # TODO: Delete else: standard_nonbonded_force.setUseSwitchingFunction(False) - #electrostatics_custom_nonbonded_force.setUseSwitchingFunction(False) # TODO: Delete sterics_custom_nonbonded_force.setUseSwitchingFunction(False) - # TODO: Delete - #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 @@ -961,8 +925,6 @@ def _nonbonded_custom_bond_force(self, sterics_energy_expression, electrostatics 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") @@ -1261,7 +1223,6 @@ 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]) # TODO: Remove once core_electrostatics_force is eliminated # 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) @@ -1275,7 +1236,6 @@ 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]) # TODO: Remove once core_electrostatics_force is eliminated # 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) @@ -1291,7 +1251,6 @@ 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. particle_index = self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge_old, 0.5*(sigma_old+sigma_new), 0.0) @@ -1307,13 +1266,11 @@ 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]) # TODO: Remove once core_electrostatics_force is eliminated #add the environment atoms to the regular nonbonded force as well: self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, epsilon) self._handle_interaction_groups() - #self._handle_hybrid_exceptions() # TODO: Delete? self._handle_original_exceptions() def _generate_dict_from_exceptions(self, force): @@ -1356,7 +1313,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'] # TODO: Delete sterics_custom_force = self._hybrid_system_forces['core_sterics_force'] #also prepare the atom classes @@ -1366,98 +1322,18 @@ def _handle_interaction_groups(self): environment_atoms = self._atom_classes['environment_atoms'] - #electrostatics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_old_atoms, core_atoms) - #electrostatics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_old_atoms, environment_atoms) - #electrostatics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_new_atoms, core_atoms) - #electrostatics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(unique_new_atoms, environment_atoms) - #electrostatics_custom_force.addInteractionGroup(core_atoms, environment_atoms) # TODO: Delete sterics_custom_force.addInteractionGroup(core_atoms, environment_atoms) - #electrostatics_custom_force.addInteractionGroup(core_atoms, core_atoms) # TODO: Delete 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") - - old_system_nonbonded_force = self._old_system_forces['NonbondedForce'] - new_system_nonbonded_force = self._new_system_forces['NonbondedForce'] - - import itertools - #prepare the atom classes - 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)) - - #add back the interactions of the old unique atoms, unless there are exceptions - for atom_pair in unique_old_pairs: - #since the pairs are indexed in the dictionary by the old system indices, we need to convert - old_index_atom_pair = (self._hybrid_to_old_map[atom_pair[0]], self._hybrid_to_old_map[atom_pair[1]]) - - #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) - - #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) - - #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 - epsilon = unit.sqrt(epsilon0*epsilon1) - sigma = 0.5*(sigma0+sigma1) - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) - - #add back the interactions of the new unique atoms, unless there are exceptions - for atom_pair in unique_new_pairs: - #since the pairs are indexed in the dictionary by the new system indices, we need to convert - new_index_atom_pair = (self._hybrid_to_new_map[atom_pair[0]], self._hybrid_to_new_map[atom_pair[1]]) - - #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) - - #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) - - #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 - epsilon = unit.sqrt(epsilon0*epsilon1) - sigma = 0.5*(sigma0+sigma1) - nonbonded_force.addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) - print("done handling exceptions") - def _handle_original_exceptions(self): """ This method ensures that exceptions present in the original systems are present in the hybrid appropriately. @@ -1484,22 +1360,12 @@ 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, # TODO: Delete this when we remove core_nonbonded_bond_force - # [chargeProd_old, sigma_old, - # epsilon_old, chargeProd_old, - # sigma_old, epsilon_old]) - self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, chargeProd_old, sigma_old, epsilon_old) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated #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 @@ -1512,18 +1378,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, # TODO: Delete this when we remove core_nonbonded_bond_force - # [chargeProd_old, sigma_old, - # epsilon_old, chargeProd_new, - # sigma_new, epsilon_new]) + #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)) - - #We also need to exclude this interaction from the custom nonbonded forces, otherwise we'll be double counting self._hybrid_system_forces['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated #now, loop through the new system to collect remaining interactions. The only that remain here are #uniquenew-uniquenew, uniquenew-core, and uniquenew-environment. @@ -1544,16 +1403,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, # TODO: Delete this when we remove core_nonbonded_bond_force - # [chargeProd_new, sigma_new, - # epsilon_new, chargeProd_new, - # sigma_new, epsilon_new]) self._hybrid_system_forces['standard_nonbonded_force'].addException(index1_hybrid, index2_hybrid, 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['core_sterics_force'].addExclusion(index1_hybrid, index2_hybrid) - #self._hybrid_system_forces['core_electrostatics_force'].addExclusion(index1_hybrid, index2_hybrid) # TODO: Remove once core_electrostatics_force is eliminated def _find_exception(self, force, index1, index2): """ diff --git a/perses/dispersed/relative_setup.py b/perses/dispersed/relative_setup.py index 773aa926a..70367766e 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]) @@ -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/scripts/setup_relative_calculation.py b/scripts/setup_relative_calculation.py index 9c31f88f2..74418d5b6 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,18 @@ setup_dict['hybrid_topology_factories']) hss = setup_dict['hybrid_sams_samplers'] + logZ = dict() free_energies = dict() - for phase in phases: + for phase in ['complex', 'solvent']: 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])) - - print("Total ddG is estimated as %.4f kT" % (free_energies['complex'] - free_energies['solvent'])) \ No newline at end of file + 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'])) + print("Total logZ ddG is estimated as %.4f kT" % (logZ['complex'] - logZ['solvent'])) From 09e853c5bdf11629b37678a8f1d21c076b429f7e Mon Sep 17 00:00:00 2001 From: John Chodera Date: Sun, 1 Jul 2018 01:23:03 -0400 Subject: [PATCH 10/18] Make sure to install correct version of OpenMM 7.3.0 --- devtools/conda-recipe/meta.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index 9dc2a3c5d..045bef31e 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -18,6 +18,7 @@ requirements: - numexpr - autograd - pymbar + - cuda92 - openmm ==7.3.0 - parmed - openmoltools @@ -43,6 +44,7 @@ requirements: - numexpr - autograd - pymbar + - cuda92 - openmm ==7.3.0 - parmed - openmoltools From 29f36f0a79a9f2373f5d7118020a4428d63f79ba Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 16:23:24 -0400 Subject: [PATCH 11/18] Restore erroneously deleted exceptions function --- perses/annihilation/new_relative.py | 87 ++++++++++++++++++++++++++++- 1 file changed, 86 insertions(+), 1 deletion(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index b33e6baf7..e2d407638 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='default', softcore_alpha=None, softcore_beta=None, annihilate_intraunique_charges=False): """ Initialize the Hybrid topology factory. @@ -73,6 +73,9 @@ 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 + annihilate_intraunique_charges: bool, optional, default False + If True, unique old/new atoms will be discharged as part of the switching process + """ self._topology_proposal = topology_proposal self._old_system = copy.deepcopy(topology_proposal.old_system) @@ -1271,6 +1274,7 @@ def handle_nonbonded(self): self._hybrid_system_forces['standard_nonbonded_force'].addParticle(charge, sigma, epsilon) self._handle_interaction_groups() + self._handle_hybrid_exceptions() self._handle_original_exceptions() def _generate_dict_from_exceptions(self, force): @@ -1334,6 +1338,87 @@ def _handle_interaction_groups(self): 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") + + old_system_nonbonded_force = self._old_system_forces['NonbondedForce'] + new_system_nonbonded_force = self._new_system_forces['NonbondedForce'] + + import itertools + #prepare the atom classes + unique_old_atoms = self._atom_classes['unique_old_atoms'] + unique_new_atoms = self._atom_classes['unique_new_atoms'] + + #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)) + + #add back the interactions of the old unique atoms, unless there are exceptions + for atom_pair in unique_old_pairs: + #since the pairs are indexed in the dictionary by the old system indices, we need to convert + old_index_atom_pair = (self._hybrid_to_old_map[atom_pair[0]], self._hybrid_to_old_map[atom_pair[1]]) + + #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] + 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]] + 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 + if self.annihilate_intraunique_charges: + chargeProd *= 0.0 + epsilon = unit.sqrt(epsilon0*epsilon1) + sigma = 0.5*(sigma0+sigma1) + 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: + #since the pairs are indexed in the dictionary by the new system indices, we need to convert + new_index_atom_pair = (self._hybrid_to_new_map[atom_pair[0]], self._hybrid_to_new_map[atom_pair[1]]) + + #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] + 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]] + 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 + if self.annihilate_intraunique_charges: + chargeProd *= 0.0 + epsilon = unit.sqrt(epsilon0*epsilon1) + sigma = 0.5*(sigma0+sigma1) + 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): """ This method ensures that exceptions present in the original systems are present in the hybrid appropriately. From c982764f6cabc944e6baab5a1e49dfeebadcfb77 Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 16:29:09 -0400 Subject: [PATCH 12/18] Fix bug in annihilate_intraunique_charges option --- perses/annihilation/new_relative.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index e2d407638..b4e19fa63 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -87,6 +87,7 @@ def __init__(self, topology_proposal, current_positions, new_positions, use_disp self._new_positions = new_positions self._use_dispersion_correction = use_dispersion_correction + self._annihilate_intraunique_charges = annihilate_intraunique_charges if softcore_alpha is None: self.softcore_alpha = 0.5 @@ -1367,12 +1368,16 @@ 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] + if self._annihilate_intraunique_charges: + 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]] + if self._annihilate_intraunique_charges: + 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 @@ -1381,7 +1386,7 @@ def _handle_hybrid_exceptions(self): [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 - if self.annihilate_intraunique_charges: + if self._annihilate_intraunique_charges: chargeProd *= 0.0 epsilon = unit.sqrt(epsilon0*epsilon1) sigma = 0.5*(sigma0+sigma1) @@ -1396,12 +1401,16 @@ 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] + if self._annihilate_intraunique_charges: + 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]] + if self._annihilate_intraunique_charges: + 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 @@ -1410,7 +1419,7 @@ def _handle_hybrid_exceptions(self): [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 - if self.annihilate_intraunique_charges: + if self._annihilate_intraunique_charges: chargeProd *= 0.0 epsilon = unit.sqrt(epsilon0*epsilon1) sigma = 0.5*(sigma0+sigma1) From 73306d0012c6da975eb1379fccf4bad75a02b339 Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 16:41:24 -0400 Subject: [PATCH 13/18] Try to fix issue with duplicate exceptions --- perses/annihilation/new_relative.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index b4e19fa63..19a599816 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -1457,7 +1457,7 @@ def _handle_original_exceptions(self): 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: + 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) From 017f1c9bd96268150d1dfc396cc732d49c90c808 Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 16:44:12 -0400 Subject: [PATCH 14/18] Fix bug in handling of unique old - unique old exceptions --- perses/annihilation/new_relative.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index 19a599816..b9d48e328 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -1456,7 +1456,11 @@ def _handle_original_exceptions(self): 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: + #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) From 47bf368002ef27b17a9bfa981491ea4a14e9d1e4 Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 16:48:37 -0400 Subject: [PATCH 15/18] Test eliminating charges from intra-unique sets --- perses/annihilation/new_relative.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index b9d48e328..be71c1c37 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, annihilate_intraunique_charges=False): + def __init__(self, topology_proposal, current_positions, new_positions, use_dispersion_correction=False, functions=None, softcore_method='default', softcore_alpha=None, softcore_beta=None, annihilate_intraunique_charges=True): """ Initialize the Hybrid topology factory. From 2069ad7584525816b11580db091e04e6d06f3c0a Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 16:58:29 -0400 Subject: [PATCH 16/18] Debug test: Comment out electrostatics from CustomBondForce --- perses/annihilation/new_relative.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index be71c1c37..afb70a3b0 100644 --- a/perses/annihilation/new_relative.py +++ b/perses/annihilation/new_relative.py @@ -928,7 +928,8 @@ def _nonbonded_custom_bond_force(self, sterics_energy_expression, electrostatics 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 = openmm.CustomBondForce("U_sterics + U_electrostatics;" + sterics_energy_expression + electrostatics_energy_expression) + custom_bond_force = openmm.CustomBondForce("U_sterics;" + sterics_energy_expression + electrostatics_energy_expression) # DEBUG custom_bond_force.addGlobalParameter("softcore_alpha", self.softcore_alpha) custom_bond_force.addGlobalParameter("softcore_beta", self.softcore_beta) custom_bond_force.addPerBondParameter("chargeprodA") From b15fbfad0955048f5c7d5b021fd4488e649a7045 Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 17:08:33 -0400 Subject: [PATCH 17/18] Update default HybridTopologyFactory treatment to "amber" for testing --- perses/annihilation/new_relative.py | 43 +---------------------------- 1 file changed, 1 insertion(+), 42 deletions(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index afb70a3b0..c4f369bb1 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, annihilate_intraunique_charges=True): + def __init__(self, topology_proposal, current_positions, new_positions, use_dispersion_correction=False, functions=None, softcore_method='amber', softcore_alpha=None, softcore_beta=None, annihilate_intraunique_charges=True): """ Initialize the Hybrid topology factory. @@ -907,47 +907,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 = openmm.CustomBondForce("U_sterics;" + sterics_energy_expression + electrostatics_energy_expression) # DEBUG - 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. From aa3bbac42ffcb8f085c6099f2800cf445f4438ab Mon Sep 17 00:00:00 2001 From: jchodera Date: Tue, 3 Jul 2018 17:20:45 -0400 Subject: [PATCH 18/18] Test both directions with test_difficult_overlap --- perses/annihilation/new_relative.py | 23 +++++++---------------- perses/tests/test_hybrid_builder.py | 18 +++++++++++++++++- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/perses/annihilation/new_relative.py b/perses/annihilation/new_relative.py index c4f369bb1..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='amber', softcore_alpha=None, softcore_beta=None, annihilate_intraunique_charges=True): + 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,8 +73,6 @@ 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 - annihilate_intraunique_charges: bool, optional, default False - If True, unique old/new atoms will be discharged as part of the switching process """ self._topology_proposal = topology_proposal @@ -87,7 +85,6 @@ def __init__(self, topology_proposal, current_positions, new_positions, use_disp self._new_positions = new_positions self._use_dispersion_correction = use_dispersion_correction - self._annihilate_intraunique_charges = annihilate_intraunique_charges if softcore_alpha is None: self.softcore_alpha = 0.5 @@ -1328,16 +1325,14 @@ 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] - if self._annihilate_intraunique_charges: - chargeProd *= 0.0 + 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]] - if self._annihilate_intraunique_charges: - chargeProd *= 0.0 + 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 @@ -1346,8 +1341,7 @@ def _handle_hybrid_exceptions(self): [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 - if self._annihilate_intraunique_charges: - chargeProd *= 0.0 + chargeProd *= 0.0 epsilon = unit.sqrt(epsilon0*epsilon1) sigma = 0.5*(sigma0+sigma1) self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) @@ -1361,16 +1355,14 @@ 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] - if self._annihilate_intraunique_charges: - chargeProd *= 0.0 + 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]] - if self._annihilate_intraunique_charges: - chargeProd *= 0.0 + 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 @@ -1379,8 +1371,7 @@ def _handle_hybrid_exceptions(self): [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 - if self._annihilate_intraunique_charges: - chargeProd *= 0.0 + chargeProd *= 0.0 epsilon = unit.sqrt(epsilon0*epsilon1) sigma = 0.5*(sigma0+sigma1) self._hybrid_system_forces['standard_nonbonded_force'].addException(atom_pair[0], atom_pair[1], chargeProd, sigma, epsilon) diff --git a/perses/tests/test_hybrid_builder.py b/perses/tests/test_hybrid_builder.py index 18e844cb9..1292d8ccb 100644 --- a/perses/tests/test_hybrid_builder.py +++ b/perses/tests/test_hybrid_builder.py @@ -325,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(current_mol_name='imatinib', proposed_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):