From 6e66933ec14062e50b5921808d76014cc661f5aa Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:13:29 -0500 Subject: [PATCH 1/7] add constraint stuff --- perses/rjmc/geometry.py | 82 ++++++++++++++++++++++++++++++++++------- 1 file changed, 68 insertions(+), 14 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index e897a9115..106637c35 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -8,6 +8,8 @@ import numpy as np import copy from perses.rjmc import coordinate_tools +import simtk.openmm as openmm + @@ -119,12 +121,17 @@ def propose(self, top_proposal, current_positions, beta): torsion_atom = torsion.atom1 #propose a bond and calculate its probability + #if it's not a bond, it's a constraint bond = self._get_relevant_bond(atom, bond_atom) - r_proposed = self._propose_bond(bond, beta) - bond_k = bond.type.k - sigma_r = units.sqrt(1/(beta*bond_k)) - logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r/sigma_r.unit)) - logp_r = self._bond_logq(r_proposed, bond, beta) - logZ_r + if type(bond)==parmed.Atom: + r_proposed = self._propose_bond(bond, beta) + bond_k = bond.type.k + sigma_r = units.sqrt(1/(beta*bond_k)) + logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r/sigma_r.unit)) + logp_r = self._bond_logq(r_proposed, bond, beta) - logZ_r + else: + r_proposed = bond #set bond length to exactly constraint + logp_r = 0.0 #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) @@ -208,13 +215,19 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates, beta): theta = internal_coordinates[1]*units.radian phi = internal_coordinates[2]*units.radian - #propose a bond and calculate its probability + #propose a bond and calculate its probability, if it's a bond. + #Otherwise, it's a constraint. bond = self._get_relevant_bond(atom, bond_atom) - bond_k = bond.type.k - sigma_r = units.sqrt(1/(beta*bond_k)) - logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r/sigma_r.unit)) #need to eliminate unit to allow numpy log - logp_r = self._bond_logq(r, bond, beta) - logZ_r - + if type(bond)==parmed.Atom: + bond_k = bond.type.k + sigma_r = units.sqrt(1/(beta*bond_k)) + logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r/sigma_r.unit)) #need to eliminate unit to allow numpy log + logp_r = self._bond_logq(r, bond, beta) - logZ_r + else: + if np.abs(r-bond) > 0: + logp_r = -np.inf + else: + logp_r = 0.0 #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) angle_k = angle.type.k @@ -247,9 +260,11 @@ def _get_proposal_order(self, topology_for_proposal, reference_topology, new_to_ """ - def _get_relevant_bond(self, atom1, atom2): + def _get_relevant_bond(self, atom1, atom2, system): """ - utility function to get the bond connecting atoms 1 and 2 + utility function to get the bond connecting atoms 1 and 2. + Returns either a bond object or a constraint distance + (since there is no constraint class) Arguments --------- @@ -257,19 +272,58 @@ def _get_relevant_bond(self, atom1, atom2): One of the atoms in the bond atom2 : parmed.atom object The other atom in the bond + system : openmm.System object + The system (containing constraints) Returns ------- bond : bond object - Bond connecting the two atoms + Bond connecting the two atoms, if there is one + constraint_parameter : unit.Quantity, nm + Constraint parameter, if there is a constraint """ bonds_1 = set(atom1.bonds) bonds_2 = set(atom2.bonds) relevant_bond_set = bonds_1.intersection(bonds_2) + if len(relevant_bond_set)==0: #this means there is a bond constraint + constraint_parameter = self._get_bond_constraint(atom1, atom2, system) + return constraint_parameter relevant_bond = relevant_bond_set.pop() relevant_bond_with_units = self._add_bond_units(relevant_bond) return relevant_bond_with_units + def _get_bond_constraint(self, atom1, atom2, system): + """ + Get the constraint parameters corresponding to the bond + between the given atoms + + Parameters + ---------- + atom1 : parmed.Atom object + the first atom of the constrained bond + atom2 : parmed.Atom object + the second atom of the constrained bond + system : openmm.System object + The system containing the constraint + + Returns + ------- + constraint : float, quantity nm + the parameters of the bond constraint + """ + atom_indices = {atom1.idx, atom2.idx} + n_constraints = system.getNumConstraints() + constraint = None + for i in range(n_constraints): + constraint_parameters = system.getConstraintParameters(i) + constraint_atoms = set(constraint_parameters[:2]) + if len(constraint_atoms.intersection(atom_indices))==2: + constraint = constraint_parameters[2] + return constraint + + + + def _get_relevant_angle(self, atom1, atom2, atom3): """ Get the angle containing the 3 given atoms From cbf2d8e11363a1cb926be6583137935896e2376c Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:18:13 -0500 Subject: [PATCH 2/7] made api better --- perses/rjmc/geometry.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 106637c35..9260ded97 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -123,14 +123,15 @@ def propose(self, top_proposal, current_positions, beta): #propose a bond and calculate its probability #if it's not a bond, it's a constraint bond = self._get_relevant_bond(atom, bond_atom) - if type(bond)==parmed.Atom: + if bond is not None: r_proposed = self._propose_bond(bond, beta) bond_k = bond.type.k sigma_r = units.sqrt(1/(beta*bond_k)) logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r/sigma_r.unit)) logp_r = self._bond_logq(r_proposed, bond, beta) - logZ_r else: - r_proposed = bond #set bond length to exactly constraint + constraint = self._get_bond_constraint(atom, bond_atom, top_proposal.new_system) + r_proposed = constraint #set bond length to exactly constraint logp_r = 0.0 #propose an angle and calculate its probability @@ -218,13 +219,14 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates, beta): #propose a bond and calculate its probability, if it's a bond. #Otherwise, it's a constraint. bond = self._get_relevant_bond(atom, bond_atom) - if type(bond)==parmed.Atom: + if bond is not None: bond_k = bond.type.k sigma_r = units.sqrt(1/(beta*bond_k)) logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r/sigma_r.unit)) #need to eliminate unit to allow numpy log logp_r = self._bond_logq(r, bond, beta) - logZ_r else: - if np.abs(r-bond) > 0: + constraint = self._get_bond_constraint(atom, bond_atom, top_proposal.old_system) + if np.abs(r-constraint) > 0: logp_r = -np.inf else: logp_r = 0.0 @@ -263,7 +265,7 @@ def _get_proposal_order(self, topology_for_proposal, reference_topology, new_to_ def _get_relevant_bond(self, atom1, atom2, system): """ utility function to get the bond connecting atoms 1 and 2. - Returns either a bond object or a constraint distance + Returns either a bond object or None (since there is no constraint class) Arguments @@ -278,16 +280,14 @@ def _get_relevant_bond(self, atom1, atom2, system): Returns ------- bond : bond object - Bond connecting the two atoms, if there is one - constraint_parameter : unit.Quantity, nm - Constraint parameter, if there is a constraint + Bond connecting the two atoms, if there is one. None if constrained or + no bond. """ bonds_1 = set(atom1.bonds) bonds_2 = set(atom2.bonds) relevant_bond_set = bonds_1.intersection(bonds_2) if len(relevant_bond_set)==0: #this means there is a bond constraint - constraint_parameter = self._get_bond_constraint(atom1, atom2, system) - return constraint_parameter + return None relevant_bond = relevant_bond_set.pop() relevant_bond_with_units = self._add_bond_units(relevant_bond) return relevant_bond_with_units From 51c2ca0ab320a36cc20c53eb48c109240e916b9b Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:18:46 -0500 Subject: [PATCH 3/7] remove incomplete function --- perses/rjmc/geometry.py | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 9260ded97..ef568e793 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -245,23 +245,6 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates, beta): old_unique_atoms.remove(atom) return logp - def _get_proposal_order(self, topology_for_proposal, reference_topology, new_to_old_atom_map): - """ - Get the order of proposal for the atoms needing proposal. Also get a list of torsions - that could be used in each atom's proposal. - - Parameters - ---------- - topology_for_proposal - reference_topology - new_to_old_atom_map - - Returns - ------- - atoms_for_proposal : dict of {atom : - - """ - def _get_relevant_bond(self, atom1, atom2, system): """ utility function to get the bond connecting atoms 1 and 2. From b5712c5de65acec96d6c581494f4cbd4dff872e9 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:19:05 -0500 Subject: [PATCH 4/7] remove extraneous parameter --- perses/rjmc/geometry.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index ef568e793..aebcf9923 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -245,7 +245,7 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates, beta): old_unique_atoms.remove(atom) return logp - def _get_relevant_bond(self, atom1, atom2, system): + def _get_relevant_bond(self, atom1, atom2): """ utility function to get the bond connecting atoms 1 and 2. Returns either a bond object or None @@ -257,8 +257,6 @@ def _get_relevant_bond(self, atom1, atom2, system): One of the atoms in the bond atom2 : parmed.atom object The other atom in the bond - system : openmm.System object - The system (containing constraints) Returns ------- From 000075c3830a052e4237a7a24cbe496ee9ac1e45 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:30:35 -0500 Subject: [PATCH 5/7] add bond constraints --- perses/rjmc/geometry.py | 4 ++-- perses/tests/test_geometry_engine.py | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index aebcf9923..54cc3ead0 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -267,9 +267,9 @@ def _get_relevant_bond(self, atom1, atom2): bonds_1 = set(atom1.bonds) bonds_2 = set(atom2.bonds) relevant_bond_set = bonds_1.intersection(bonds_2) - if len(relevant_bond_set)==0: #this means there is a bond constraint - return None relevant_bond = relevant_bond_set.pop() + if relevant_bond.type is None: + return None relevant_bond_with_units = self._add_bond_units(relevant_bond) return relevant_bond_with_units diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 52097bb28..900026542 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -67,7 +67,7 @@ def oemol_to_openmm_system(oemol, molecule_name): from perses.rjmc import topology_proposal from openmoltools import forcefield_generators gaff_xml_filename = get_data_filename('data/gaff.xml') - system_generator = topology_proposal.SystemGenerator([gaff_xml_filename]) + system_generator = topology_proposal.SystemGenerator([gaff_xml_filename], forcefield_kwargs={'constraints' : app.HBonds}) topology = forcefield_generators.generateTopologyFromOEMol(oemol) system = system_generator.build_system(topology) positions = extractPositionsFromOEMOL(oemol) @@ -133,8 +133,8 @@ def test_run_geometry_engine(): molecule2 = generate_initial_molecule(molecule_name_2) new_to_old_atom_mapping = align_molecules(molecule1, molecule2) - sys1, pos1, top1 = oemol_to_openmm_system_amber(molecule1, molecule_name_1) - sys2, pos2, top2 = oemol_to_openmm_system_amber(molecule2, molecule_name_2) + sys1, pos1, top1 = oemol_to_openmm_system(molecule1, molecule_name_1) + sys2, pos2, top2 = oemol_to_openmm_system(molecule2, molecule_name_2) import perses.rjmc.geometry as geometry import perses.rjmc.topology_proposal as topology_proposal @@ -152,7 +152,7 @@ def test_run_geometry_engine(): state = context.getState(getEnergy=True) print("Energy before proposal is: %s" % str(state.getPotentialEnergy())) - new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal) + new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal, pos1, beta) app.PDBFile.writeFile(top2, new_positions, file=test_pdb_file) test_pdb_file.close() context.setPositions(new_positions) @@ -375,9 +375,9 @@ def test_angle(): if __name__=="__main__": #test_coordinate_conversion() - #test_run_geometry_engine() + test_run_geometry_engine() #test_existing_coordinates() #test_openmm_dihedral() #test_try_random_itoc() #test_angle() - test_logp_reverse() \ No newline at end of file + #test_logp_reverse() \ No newline at end of file From 7d81ec2deaae4309bea3311074967ddb70633933 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:39:07 -0500 Subject: [PATCH 6/7] fix constraints --- perses/rjmc/geometry.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 54cc3ead0..fd46868bc 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -226,10 +226,9 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates, beta): logp_r = self._bond_logq(r, bond, beta) - logZ_r else: constraint = self._get_bond_constraint(atom, bond_atom, top_proposal.old_system) - if np.abs(r-constraint) > 0: - logp_r = -np.inf - else: - logp_r = 0.0 + if constraint is None: + raise Exception("No Bond or constraint between atom %d and %d" % (atom.idx, bond_atom.idx)) + logp_r = 0.0 #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) angle_k = angle.type.k From 4a77185401188de7c37079cabfc258f9e1bdd238 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 12 Feb 2016 17:41:30 -0500 Subject: [PATCH 7/7] fix test --- perses/tests/test_geometry_engine.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 900026542..eefa96b41 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -298,7 +298,7 @@ def test_logp_reverse(): import perses.rjmc.topology_proposal as topology_proposal sm_top_proposal = topology_proposal.TopologyProposal(new_topology=top2, new_system=sys2, old_topology=top1, old_system=sys1, - logp_proposal=0.0, new_to_old_atom_map=new_to_old_atom_mapping, metadata={'test':0.0}) + logp_proposal=0.0, new_to_old_atom_map=new_to_old_atom_mapping, new_chemical_state_key="CCC", old_chemical_state_key="CC", metadata={'test':0.0}) geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal, pos1, beta) @@ -306,7 +306,7 @@ def test_logp_reverse(): #reverse the atom map: old_to_new_atom_mapping = {value : key for key, value in new_to_old_atom_mapping.items()} sm_reverse_proposal = topology_proposal.TopologyProposal(new_topology=top1, new_system=sys1, old_topology=top2, old_system=sys2, - logp_proposal=0.0, new_to_old_atom_map=old_to_new_atom_mapping, metadata={'test':0.0}) + logp_proposal=0.0, new_to_old_atom_map=old_to_new_atom_mapping, new_chemical_state_key="CC", old_chemical_state_key="CCCC", metadata={'test':0.0}) logp_reverse = geometry_engine.logp_reverse(sm_reverse_proposal, pos1, new_positions, beta) print(logp_proposal) print(logp_reverse) @@ -380,4 +380,4 @@ def test_angle(): #test_openmm_dihedral() #test_try_random_itoc() #test_angle() - #test_logp_reverse() \ No newline at end of file + test_logp_reverse() \ No newline at end of file