From a6b5636e2cc92cd937dcb024222c9a8512e10544 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 14:09:12 -0500 Subject: [PATCH 01/57] replace rejection sampling with selection from precomputed values --- perses/rjmc/geometry.py | 28 +++++++++------------------- perses/tests/test_geometry_engine.py | 12 +++++++----- 2 files changed, 16 insertions(+), 24 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index e987e0b62..6c0f03626 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -10,6 +10,7 @@ import parmed import simtk.unit as units import logging +from numba import jit GeometryProposal = namedtuple('GeometryProposal',['new_positions','logp']) @@ -584,7 +585,7 @@ def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, in bond_atom_position = xyz if torsion.atom2 == atom else positions[torsion.atom2.idx] angle_atom_position = xyz if torsion.atom3 == atom else positions[torsion.atom3.idx] torsion_atom_position = xyz if torsion.atom4 == atom else positions[torsion.atom4.idx] - internal_coordinates, _ = self._autograd_ctoi(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position, calculate_jacobian=False) + internal_coordinates, _ = _ctoi_prim_numba(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position) phi = internal_coordinates[2]*units.radians ub_torsions += self._torsion_potential(torsion, phi, beta) return ub_angles+ub_torsions @@ -640,7 +641,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor #get the normalized probabilities for torsions p = q / Z - return p, Z + return p, Z, q, phis def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_position): """ @@ -680,25 +681,14 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, Propose a torsion angle, including energetic contributions from other torsions and angles """ #first, let's get the normalizing constant of this distribution - p, Z = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) - #now, use rejection sampling on the normalized distribution to randomly draw a phi - max_p = max(p) - logp = 0.0 - phi = 0 - accepted = False - while not accepted: - phi_samp = np.random.uniform(0.0, 2*np.pi)*units.radians - xyz, _ = self._autograd_itoc(bond_atom.idx, angle_atom.idx, torsion_atom.idx, r, theta, phi_samp, positions) - runif = np.random.uniform(0.0, max_p+1.0) - p_phi_samp = self._torsion_logp(atom, xyz, torsion, atoms_with_positions, positions, beta, Z) - if p_phi_samp > runif: - phi = phi_samp - logp = p_phi_samp - accepted = True - else: - continue + p, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + #choose from the set of possible torsion angles + phi_idx = np.random.choice(range(len(phis)), p=p) + logp = np.log(p[phi_idx]) + phi = phis[phi_idx] return phi, logp + class FFGeometryEngineOld(GeometryEngine): """ This is an implementation of the GeometryEngine class which proposes new dimensions based on the forcefield parameters diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index cb8ee846f..d4a116eac 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -89,12 +89,14 @@ 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) - context.setPositions(new_positions) - state2 = context.getState(getEnergy=True) - print("Energy after proposal is: %s" %str(state2.getPotentialEnergy())) + for i in range(10): + new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal) + context.setPositions(new_positions) + state2 = context.getState(getEnergy=True) + print("Energy after proposal is: %s" %str(state2.getPotentialEnergy())) if __name__=="__main__": - test_run_geometry_engine() \ No newline at end of file + for i in range(10): + test_run_geometry_engine() \ No newline at end of file From 8dbb9c14b4dca43d292e9787cc6bd1f57893c9ad Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 15:34:02 -0500 Subject: [PATCH 02/57] fixing probability calculations --- perses/rjmc/geometry.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 6c0f03626..c9369ab57 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -439,7 +439,7 @@ def _internal_to_cartesian(rthetaphi): #logging.debug("The spherical detJ is %f" % detj_spherical) return units.Quantity(atom_xyz, unit=units.nanometers), np.abs(jacobian_det) - def _bond_logp(self, r, bond, beta): + def _bond_logq(self, r, bond, beta): """ Calculate the log-probability of a given bond at a given inverse temperature @@ -456,11 +456,10 @@ def _bond_logp(self, r, bond, beta): """ k_eq = bond.type.k*units.kilocalories_per_mole/(units.angstrom**2) r0 = bond.type.req*units.nanometers - sigma = beta*2.0/np.sqrt(2.0*k_eq/k_eq.unit) - logp = stats.distributions.norm.logpdf(r/r.unit, r0/r0.unit, sigma) - return logp + logq = -beta*k_eq*(r-r0)**2 + return logq - def _angle_logp(self, theta, angle, beta): + def _angle_logq(self, theta, angle, beta): """ Calculate the log-probability of a given bond at a given inverse temperature @@ -475,11 +474,10 @@ def _angle_logp(self, theta, angle, beta): """ k_eq = angle.type.k*units.kilocalories_per_mole/(units.radians**2) theta0 = angle.type.theteq*units.radians - sigma = beta*2.0/np.sqrt(2.0*k_eq/k_eq.unit) - logp = stats.distributions.norm.logpdf(theta/theta.unit, theta0/theta0.unit, sigma) - return logp + logq = -beta*k_eq*(theta-theta0)**2 + return logq - def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta): + def _torsion_logq(self, atom, xyz, torsion, atoms_with_positions, positions, beta): """ Utility function for calculating the unnormalized probability of a torsion angle """ @@ -497,8 +495,8 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet gamma = torsion.type.phase V = torsion.type.phi_k n = torsion.type.per - q = np.exp(-beta*(V/2.0)*(1+np.cos(n*internal_coordinates[2]-gamma))) - return q + logq = -beta*(V/2.0)*(1+np.cos(n*internal_coordinates[2]-gamma)) + return logq def _choose_torsion(self, atoms_with_positions, atom_for_proposal): """ @@ -671,7 +669,7 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet involved_torsions = self._get_torsions(atoms_with_positions, atom) internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) if not Z: - p, Z = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + p, Z, _, _ = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) ub_torsion = self._torsion_and_angle_potential(xyz, atom, positions, involved_angles, involved_torsions, beta) p_torsion = np.exp(-ub_torsion) / Z return p_torsion From b11e9587302f785f550e03c87ce919afcc996622 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 15:36:00 -0500 Subject: [PATCH 03/57] moved old geometry engine to the attic --- attic/geometry_attic.py | 638 +++++++++++++++++++++++++++++++++++++++- perses/rjmc/geometry.py | 638 +--------------------------------------- 2 files changed, 638 insertions(+), 638 deletions(-) diff --git a/attic/geometry_attic.py b/attic/geometry_attic.py index c0af44ca1..9b016001c 100644 --- a/attic/geometry_attic.py +++ b/attic/geometry_attic.py @@ -204,4 +204,640 @@ def _cartesian_to_internal(self, atom, bond_atom, angle_atom, torsion_atom, posi detJ = sympy.lambdify("Xa, Ya, Za",jacobian_determinant, "numpy")(*atomic_position) - return np.array([r, theta, phi]), np.abs(detJ) \ No newline at end of file + return np.array([r, theta, phi]), np.abs(detJ) + + +class FFGeometryEngineOld(GeometryEngine): + """ + This is an implementation of the GeometryEngine class which proposes new dimensions based on the forcefield parameters + """ + + + def propose(self, top_proposal): + """ + Propose a new geometry for the appropriate atoms using forcefield parameters + + Arguments + ---------- + topology_proposal : TopologyProposal object + Object containing the relevant results of a topology proposal + + Returns + ------- + new_positions : [m, 3] np.array of floats + The positions of the m atoms in the new system (if not propose_positions, this is the old positions) + logp : float + The logp of the proposal, including the jacobian + """ + + #get the mapping between new and old atoms + n_atoms_new = top_proposal.n_atoms_new + n_atoms_old = top_proposal.n_atoms_old + + #get a list of the new atoms + new_atoms = top_proposal.unique_new_atoms + + #get a list of the old atoms + old_atoms = top_proposal.unique_old_atoms + + #create a new array with the appropriate size + new_positions = np.zeros([n_atoms_new, 3]) + + new_to_old_atom_map = top_proposal.new_to_old_atom_map + + #transfer known positions--be careful about units! + for atom in new_to_old_atom_map.keys(): + new_positions[atom] = top_proposal.old_positions[new_to_old_atom_map[atom]].in_units_of(units.nanometers) + + #restore units + new_positions = new_positions*units.nanometers + #get a new set of positions and a logp of that proposal + final_new_positions, logp_forward = self._propose_new_positions(new_atoms, top_proposal.new_system, new_positions, new_to_old_atom_map) + + #get the probability of a reverse proposal + #first get the atom map reversed + old_to_new_atom_map = top_proposal.old_to_new_atom_map + _, logp_reverse = self._propose_new_positions(old_atoms, top_proposal.old_system, top_proposal.old_positions, old_to_new_atom_map, propose_positions=False) + + #construct the return object + geometry_proposal = GeometryProposal(final_new_positions, logp_reverse - logp_forward) + + return geometry_proposal + + def _propose_new_positions(self, new_atoms, new_system, new_positions, new_to_old_atom_map, propose_positions=True): + """ + Propose new atomic positions, or get the probability of existing ones + (used to evaluate the reverse probability, where the old atoms become the new atoms, etc) + Arguments + --------- + new_atoms : list of int + Indices of atoms that need positions + topology_proposal : TopologyProposal namedtuple + The result of the topology proposal, containing the atom mapping and topologies + new_system : simtk.OpenMM.System object + The new system + new_positions : [m, 3] np.array of floats + The positions of the m atoms in the new system, with unpositioned atoms having [0,0,0] + + Returns + ------- + new_positions : [m, 3] np.array of floats + The positions of the m atoms in the new system (if not propose_positions, this is the old positions) + logp : float + The logp of the proposal, including the jacobian + """ + logp = 0.0 + atoms_with_positions = new_to_old_atom_map.keys() + forces = {new_system.getForce(index).__class__.__name__ : new_system.getForce(index) for index in range(new_system.getNumForces())} + torsion_force = forces['PeriodicTorsionForce'] + bond_force = forces['HarmonicBondForce'] + atomic_proposal_order = self._generate_proposal_order(new_atoms, new_system, new_to_old_atom_map) + for atomset in atomic_proposal_order: + for atom, possible_torsions in atomset.items(): + #randomly select a torsion from the available list, and calculate its log-probability + selected_torsion_index = np.random.randint(0, len(possible_torsions)) + selected_torsion = possible_torsions[selected_torsion_index] + logp_selected_torsion = - np.log(len(possible_torsions)) + torsion_parameters = torsion_force.getTorsionParameters(selected_torsion) + #check to see whether the atom is the first or last in the torsion + if atom == torsion_parameters[0]: + bonded_atom = torsion_parameters[1] + angle_atom = torsion_parameters[2] + torsion_atom = torsion_parameters[3] + else: + bonded_atom = torsion_parameters[2] + angle_atom = torsion_parameters[1] + torsion_atom = torsion_parameters[0] + + #get a new position for the atom, if specified + if propose_positions: + atomic_xyz, logp_atomic = self._propose_atomic_position(new_system, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters, new_positions) + new_positions[atom] = atomic_xyz + else: + logp_atomic = 0.0 + logp_atomic = self._calculate_logp_atomic_position(new_system, new_positions, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters) + logp += logp_atomic + logp_selected_torsion + + return new_positions, logp + + def _generate_proposal_order(self, new_atoms, new_system, new_to_old_atom_map): + """ + Generate the list of atoms (and corresponding torsions) to be proposed each round + + Arguments + --------- + new_atoms : list of int + Indices of atoms that need positions + topology_proposal : TopologyProposal namedtuple + The result of the topology proposal, containing the atom mapping and topologies + new_system : simtk.OpenMM.System object + The new system + new_positions : [m, 3] np.array of floats + The positions of the m atoms in the new system, with unpositioned atoms having [0,0,0] + + Returns + ------- + atom_proposal_list : list of dict + List of dict of form {atom_index : [torsion_list]}. each entry in the list is a + round of atoms eligible for proposal. + """ + atoms_with_positions = new_to_old_atom_map.keys() + forces = {new_system.getForce(index).__class__.__name__ : new_system.getForce(index) for index in range(new_system.getNumForces())} + torsion_force = forces['PeriodicTorsionForce'] + bond_force = forces['HarmonicBondForce'] + atom_proposal_list = [] + while len(new_atoms) > 0: + atom_torsion_proposals = self._atoms_eligible_for_proposal(torsion_force, bond_force, atoms_with_positions, new_atoms) + #now loop through the list of atoms found to be eligible for positions this round + atom_proposal_list.append(atom_torsion_proposals) + for atom in atom_torsion_proposals.keys(): + #remove atom from list of atoms needing positions + new_atoms.remove(atom) + atoms_with_positions.append(atom) + return atom_proposal_list + + + + def _propose_atomic_position(self, system, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters, positions): + """ + Method to propose the position of a single atom. + + Arguments + --------- + system : simtk.openmm.System + the system containing the atom of interest + atom : int + the index of the atom whose position should be proposed + bonded_atom : int + the index of the atom bonded to the atom of interest + angle_atom : int + the 1, 3 position atom index + torsion_parameters : dict + The parameters of the torsion that is being used for this proposal + + Returns + ------- + atomic_xyz : [1,3] np.array of floats + the cartesian coordinates of the atom + logp_atomic_position : float + the logp of this proposal, including jacobian correction + """ + logging.debug("Proposing position %d with bond atom %d, angle atom %d, torsion atom %d" % (atom, bonded_atom, angle_atom, torsion_atom)) + #get bond parameters and draw bond length + r0, k_eq = self._get_bond_parameters(system, atom, bonded_atom) + r, logp_bond = self._propose_bond_length(r0, k_eq) + + #get angle parameters and draw bond angle + theta0, k_eq = self._get_angle_parameters(system, atom, bonded_atom, angle_atom) + theta, logp_angle = self._propose_bond_angle(theta0, k_eq) + + #propose torsion + phi, logp_torsion = self._propose_torsion_angle(torsion_parameters[6], torsion_parameters[4], torsion_parameters[5]) + #convert spherical to cartesian coordinates + spherical_coordinates = np.asarray([r, theta, phi]) + atomic_xyz, detJ = self._autograd_itoc(bonded_atom, angle_atom, torsion_atom, r, theta, phi, positions) + #accumulate the forward logp with jacobian correction + logp_atomic_position = (logp_angle + logp_bond + logp_torsion + np.log(detJ)) + logging.debug("Proposed (r, theta, phi) of (%s, %s, %s)" % (str(r), str(theta), str(phi))) + return atomic_xyz, logp_atomic_position + + + def _calculate_logp_atomic_position(self, system, positions, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters): + """ + Calculate the log-probability of a given atom's position + + Arguments + --------- + system : simtk.openmm.System object + The system whose atomic position needs a logp + positions : [n, 3] np.ndarray of floats, Quantity nm + The positions of the particles in the system + atom : int + The index of the atom of interest + bonded_atom : int + The index of the bonded atom + angle_atom : int + The index of the angle atom + torsion_atom : int + The index of the torsion atom + torsion_parameters : array + The parameters of the torsion + + Returns + ------- + logp : float + logp of position, including det(J) correction + """ + #convert cartesian coordinates to spherical + relative_xyz = positions[atom] - positions[bonded_atom] + internal_coordinates, detJ = self._autograd_ctoi(atom, bonded_atom, angle_atom, torsion_atom, positions) + + #get bond parameters and calculate the probability of choosing the current positions + r0, k_eq = self._get_bond_parameters(system, atom, bonded_atom) + sigma_bond = 1.0/np.sqrt(2.0*k_eq/k_eq.unit) + logp_bond = self._ + + #calculate the probability of choosing this bond angle + theta0, k_eq = self._get_angle_parameters(system, atom, bonded_atom, angle_atom) + sigma_angle = 1.0/np.sqrt(2.0*k_eq/k_eq.unit) + logp_angle = stats.distributions.norm.logpdf(internal_coordinates[1], theta0/theta0.unit, sigma_angle) + + #calculate the probabilty of choosing this torsion angle + torsion_Z, _ = self._torsion_normalizer(torsion_parameters[6], torsion_parameters[4], torsion_parameters[5]) + logp_torsion = self._torsion_p(torsion_Z,torsion_parameters[6], torsion_parameters[4], torsion_parameters[5], internal_coordinates[2]) + + logp_atomic_position = (logp_angle + logp_bond + logp_torsion + np.log(detJ)) + + return logp_atomic_position + + def _atoms_eligible_for_proposal(self, torsion_force, bond_force, atoms_with_positions, new_atoms): + """ + Determine which atoms are eligible for position proposal this round. + The criteria are that the atom must not yet have a position, and must + be 1 or 4 in at least one torsion where the other atoms have positoins. + + Arguments + --------- + torsion_force : simtk.openmm.PeriodicTorsionForce + The torsion force from the system + bond_roce : simtk.openmm.HarmonicBondForce + The HarmonicBondForce from the system + atoms_with_positions : list of int + A list of the indices of atoms that have positions + + Returns + ------- + atom_torsion_proposals : dict + Dictionary containing {eligible_atom : [list_of_torsions]} + + """ + current_atom_proposals = [] + atom_torsion_proposals = {atom_index: [] for atom_index in new_atoms} + for torsion_index in range(torsion_force.getNumTorsions()): + atom1 = torsion_force.getTorsionParameters(torsion_index)[0] + atom2 = torsion_force.getTorsionParameters(torsion_index)[1] + atom3 = torsion_force.getTorsionParameters(torsion_index)[2] + atom4 = torsion_force.getTorsionParameters(torsion_index)[3] + #Only take torsions where the "new" statuses of atoms 1 and 4 are not equal + if (atom1 in atoms_with_positions) != (atom4 in atoms_with_positions): + #make sure torsion is not improper: + if not self._is_proper(torsion_force.getTorsionParameters(torsion_index), bond_force): + continue + #only take torsions where 2 and 3 have known positions + if (atom2 in atoms_with_positions) and (atom3 in atoms_with_positions): + #finally, append the torsion index to the list of possible atom sets for proposal + if atom1 in atoms_with_positions and (atom1 not in current_atom_proposals): + atom_torsion_proposals[atom4].append(torsion_index) + current_atom_proposals.append(atom4) + elif atom4 in atoms_with_positions and (atom4 not in current_atom_proposals): + atom_torsion_proposals[atom1].append(torsion_index) + current_atom_proposals.append(atom1) + else: + continue + return {atom: torsionlist for atom, torsionlist in atom_torsion_proposals.items() if len(torsionlist) > 0} + + + def _is_proper(self, torsion_parameters, bond_force): + """ + Utility function to determine if torsion is proper + """ + #get the atoms + torsion_atoms = torsion_parameters[:4] + is_proper = True + for i in range(3): + is_proper = is_proper and self._is_bond(torsion_atoms[i], torsion_atoms[i+1], bond_force) + return is_proper + + def _is_bond(self, atom1, atom2, bond_force): + """ + Utility function to determine if bond exists between two atoms + """ + for bond_index in range(bond_force.getNumBonds()): + parameters = bond_force.getBondParameters(bond_index) + if (parameters[0] == atom1 and parameters[1] == atom2) or (parameters[1] == atom1 and parameters[0] == atom2): + return True + return False + + + def _get_bond_parameters(self, system, atom1, atom2): + """ + Utility function to get the bonded parameters for a pair of atoms + + Arguments + --------- + system : simtk.openmm.System + the system containing the parameters of interest + atom1 : int + the index of the first atom in the bond + atom2 : int + the index of the second atom in the bond + + Returns + ------- + r0 : float + equilibrium bond length + k_eq : float + bond spring constant + """ + forces = {system.getForce(index).__class__.__name__ : system.getForce(index) for index in range(system.getNumForces())} + bonded_force = forces['HarmonicBondForce'] + for bond_index in range(bonded_force.getNumBonds()): + parameters = bonded_force.getBondParameters(bond_index) + if (parameters[0] == atom1 and parameters[1] == atom2) or (parameters[1] == atom1 and parameters[0] == atom2): + return parameters[2], parameters[3] + + def _get_angle_parameters(self, system, atom1, atom2, atom3): + """ + Utility function to retrieve harmonic angle parameters for + a given set of atoms in a system + + Arguments + --------- + system : simtk.openmm.System + the system with parameters + atom1 : int + the index of the first atom + atom2 : int + the index of the second atom + atom3 : int + the index of the third atom + + Returns + ------- + theta0 : float + equilibrium bond angle + k_eq : float + angle spring constant + """ + list_of_angle_atoms = [atom1, atom2, atom3] + #compare sorted lists of atoms + forces = {system.getForce(index).__class__.__name__ : system.getForce(index) for index in range(system.getNumForces())} + angle_force = forces['HarmonicAngleForce'] + for angle_index in range(angle_force.getNumAngles()): + parameters = angle_force.getAngleParameters(angle_index) + #the first three "parameters" are atom indices + atoms = parameters[:3] + if np.all(np.equal(list_of_angle_atoms, atoms)) or np.all(np.equal(list_of_angle_atoms[::-1], atoms)): + return parameters[3], parameters[4] + + def _get_unique_atoms(self, new_to_old_map, new_atom_list, old_atom_list): + """ + Get the set of atoms unique to both the new and old system. + + Arguments + --------- + new_to_old_map : dict + Dictionary of the form {new_atom : old_atom} + + Returns + ------- + unique_old_atoms : list of int + A list of the indices of atoms unique to the old system + unique_new_atoms : list of int + A list of the indices of atoms unique to the new system + """ + mapped_new_atoms = new_to_old_map.keys() + unique_new_atoms = [atom for atom in new_atom_list not in mapped_new_atoms] + mapped_old_atoms = new_to_old_map.values() + unique_old_atoms = [atom for atom in old_atom_list not in mapped_old_atoms] + return unique_old_atoms, unique_new_atoms + + def _propose_bond_length(self, r0, k_eq): + """ + Draw a bond length from the equilibrium bond length distribution + + Arguments + --------- + r0 : float + The equilibrium bond length + k_eq : float + The bond spring constant + + Returns + -------- + r : float + the proposed bond length + logp : float + the log-probability of the proposal + """ + sigma = 2.0/np.sqrt(2.0*k_eq/k_eq.unit) + r = sigma*np.random.random()*r0.unit + r0 + logp = stats.distributions.norm.logpdf(r/r.unit, r0/r0.unit, sigma) + return (r, logp) + + def _propose_bond_angle(self, theta0, k_eq): + """ + Draw a bond length from the equilibrium bond length distribution + + Arguments + --------- + theta0 : float + The equilibrium bond angle + k_eq : float + The angle spring constant + + Returns + -------- + r : float + the proposed bond angle + logp : float + the log-probability of the proposal + """ + sigma = 1.0/np.sqrt(2.0*k_eq/k_eq.unit) + theta = sigma*np.random.random()*theta0.unit + theta0 + logp = stats.distributions.norm.logpdf(theta / theta.unit, theta0/theta0.unit, sigma) + return (theta, logp) + + def _propose_torsion_angle(self, V, n, gamma): + """ + Draws a torsion angle from the distribution of one torsion. + Uses the functional form U(phi) = V/2[1+cos(n*phi-gamma)], + as described by Amber. + + Arguments + --------- + V : float + force constant + n : int + multiplicity of the torsion + gamma : float + phase angle + + Returns + ------- + phi : float + The proposed torsion angle + logp : float + The proposed log probability of the torsion angle + """ + (Z, max_p) = self._torsion_normalizer(V, n, gamma) + phi = 0.0 + logp = 0.0 + #sample from the distribution + accepted = False + #use rejection sampling + while not accepted: + phi_samp = np.random.uniform(0.0, 2*np.pi) + runif = np.random.uniform(0.0, max_p+1.0) + p_phi_samp = self._torsion_p(Z, V, n, gamma, phi_samp) + if p_phi_samp > runif: + phi = phi_samp + logp = np.log(p_phi_samp) + accepted = True + else: + continue + return (phi*units.radians, logp) + + def _torsion_p(self, Z, V, n, gamma, phi): + """ + Utility function for calculating the normalized probability of a torsion angle + """ + #must remove units + V = V/V.unit + gamma = gamma/gamma.unit + return (1.0/Z)*np.exp(-(V/2.0)*(1+np.cos(n*phi-gamma))) + + def _torsion_normalizer(self, V, n, gamma): + """ + Utility function to numerically normalize torsion angles. + Also return max_p to facilitate rejection sampling + """ + #generate a grid of 5000 points from 0 < phi < 2*pi + phis = np.linspace(0, 2.0*np.pi, 5000) + #evaluate the unnormalized probability at each of those points + #need to remove units--numexpr can't handle them otherwise + V = V/V.unit + gamma = gamma/gamma.unit + phi_q = ne.evaluate("exp(-(V/2.0)*(1+cos(n*phis-gamma)))") + #integrate the values + Z = np.trapz(phi_q, phis) + max_p = np.max(phi_q/Z) + return Z, max_p + + + + def _autograd_ctoi(self, atom, bond, angle, torsion, positions): + import autograd + import autograd.numpy as np + + positions = positions/positions.unit + atom_position = positions[atom] + + + def _cartesian_to_internal(xyz): + """ + Autograd-based jacobian of transformation from cartesian to internal. + Returns without units! + """ + a = positions[bond] - xyz + b = positions[angle] - positions[bond] + #3-4 bond + c = positions[angle] - positions[torsion] + a_u = a / np.linalg.norm(a) + b_u = b / np.linalg.norm(b) + c_u = c / np.linalg.norm(c) + + #bond length + r = np.linalg.norm(a) + + #bond angle + theta = np.arccos(np.dot(-a_u, b_u)) + + #torsion angle + plane1 = np.cross(a, b) + plane2 = np.cross(b, c) + + phi = np.arccos(np.dot(plane1, plane2)/(np.linalg.norm(plane1)*np.linalg.norm(plane2))) + + if np.dot(np.cross(plane1, plane2), b_u) < 0: + phi = -phi + + return np.array([r, theta, phi]) + + j = autograd.jacobian(_cartesian_to_internal) + internal_coords = _cartesian_to_internal(atom_position) + jacobian_det = np.linalg.det(j(atom_position)) + return internal_coords, np.abs(jacobian_det) + + def _autograd_itoc(self, bond, angle, torsion, r, theta, phi, positions): + """ + Autograd based coordinate conversion internal -> cartesian + + Arguments + --------- + bond : int + index of the bonded atom + angle : int + index of the angle atom + torsion : int + index of the torsion atom + r : float, Quantity nm + bond length + theta : float, Quantity rad + bond angle + phi : float, Quantity rad + torsion angle + positions : [n, 3] np.array Quantity nm + positions of the atoms in the molecule + + Returns + ------- + atom_xyz : [1, 3] np.array Quantity nm + The atomic positions in cartesian space + detJ : float + The absolute value of the determinant of the jacobian + """ + import autograd + import autograd.numpy as np + positions = positions/positions.unit + r = r/r.unit + theta = theta/theta.unit + phi = phi/phi.unit + rtp = np.array([r, theta, phi]) + + def _internal_to_cartesian(rthetaphi): + + a = positions[angle] - positions[bond] + b = positions[angle] - positions[torsion] + + a_u = a / np.linalg.norm(a) + b_u = b / np.linalg.norm(b) + + + d_r = rthetaphi[0]*a_u + + normal = np.cross(a, b) + + #construct the angle rotation matrix + axis_angle = normal / np.linalg.norm(normal) + a = np.cos(rthetaphi[1]/2) + b, c, d = -axis_angle*np.sin(rthetaphi[1]/2) + angle_rotation_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], + [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], + [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) + #apply it + d_ang = np.dot(angle_rotation_matrix, d_r) + + #construct the torsion rotation matrix and apply it + axis = a_u + a = np.cos(rthetaphi[2]/2) + b, c, d = -axis*np.sin(rthetaphi[2]/2) + torsion_rotation_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], + [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], + [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) + #apply it + d_torsion = np.dot(torsion_rotation_matrix, d_ang) + + #add the positions of the bond atom + xyz = positions[bond] + d_torsion + + return xyz + + j = autograd.jacobian(_internal_to_cartesian) + atom_xyz = _internal_to_cartesian(rtp) + jacobian_det = np.linalg.det(j(rtp)) + logging.debug("detJ is %f" %(jacobian_det)) + detj_spherical = r**2*np.sin(theta) + logging.debug("The spherical detJ is %f" % detj_spherical) + return units.Quantity(atom_xyz, unit=units.nanometers), np.abs(jacobian_det) + diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index c9369ab57..1fa917604 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -5,7 +5,6 @@ from collections import namedtuple import perses.rjmc.topology_proposal as topology_proposal import numpy as np -import scipy.stats as stats import numexpr as ne import parmed import simtk.unit as units @@ -583,7 +582,7 @@ def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, in bond_atom_position = xyz if torsion.atom2 == atom else positions[torsion.atom2.idx] angle_atom_position = xyz if torsion.atom3 == atom else positions[torsion.atom3.idx] torsion_atom_position = xyz if torsion.atom4 == atom else positions[torsion.atom4.idx] - internal_coordinates, _ = _ctoi_prim_numba(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position) + internal_coordinates, _ = self._autograd_ctoi(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position) phi = internal_coordinates[2]*units.radians ub_torsions += self._torsion_potential(torsion, phi, beta) return ub_angles+ub_torsions @@ -686,638 +685,3 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, phi = phis[phi_idx] return phi, logp - -class FFGeometryEngineOld(GeometryEngine): - """ - This is an implementation of the GeometryEngine class which proposes new dimensions based on the forcefield parameters - """ - - - def propose(self, top_proposal): - """ - Propose a new geometry for the appropriate atoms using forcefield parameters - - Arguments - ---------- - topology_proposal : TopologyProposal object - Object containing the relevant results of a topology proposal - - Returns - ------- - new_positions : [m, 3] np.array of floats - The positions of the m atoms in the new system (if not propose_positions, this is the old positions) - logp : float - The logp of the proposal, including the jacobian - """ - - #get the mapping between new and old atoms - n_atoms_new = top_proposal.n_atoms_new - n_atoms_old = top_proposal.n_atoms_old - - #get a list of the new atoms - new_atoms = top_proposal.unique_new_atoms - - #get a list of the old atoms - old_atoms = top_proposal.unique_old_atoms - - #create a new array with the appropriate size - new_positions = np.zeros([n_atoms_new, 3]) - - new_to_old_atom_map = top_proposal.new_to_old_atom_map - - #transfer known positions--be careful about units! - for atom in new_to_old_atom_map.keys(): - new_positions[atom] = top_proposal.old_positions[new_to_old_atom_map[atom]].in_units_of(units.nanometers) - - #restore units - new_positions = new_positions*units.nanometers - #get a new set of positions and a logp of that proposal - final_new_positions, logp_forward = self._propose_new_positions(new_atoms, top_proposal.new_system, new_positions, new_to_old_atom_map) - - #get the probability of a reverse proposal - #first get the atom map reversed - old_to_new_atom_map = top_proposal.old_to_new_atom_map - _, logp_reverse = self._propose_new_positions(old_atoms, top_proposal.old_system, top_proposal.old_positions, old_to_new_atom_map, propose_positions=False) - - #construct the return object - geometry_proposal = GeometryProposal(final_new_positions, logp_reverse - logp_forward) - - return geometry_proposal - - def _propose_new_positions(self, new_atoms, new_system, new_positions, new_to_old_atom_map, propose_positions=True): - """ - Propose new atomic positions, or get the probability of existing ones - (used to evaluate the reverse probability, where the old atoms become the new atoms, etc) - Arguments - --------- - new_atoms : list of int - Indices of atoms that need positions - topology_proposal : TopologyProposal namedtuple - The result of the topology proposal, containing the atom mapping and topologies - new_system : simtk.OpenMM.System object - The new system - new_positions : [m, 3] np.array of floats - The positions of the m atoms in the new system, with unpositioned atoms having [0,0,0] - - Returns - ------- - new_positions : [m, 3] np.array of floats - The positions of the m atoms in the new system (if not propose_positions, this is the old positions) - logp : float - The logp of the proposal, including the jacobian - """ - logp = 0.0 - atoms_with_positions = new_to_old_atom_map.keys() - forces = {new_system.getForce(index).__class__.__name__ : new_system.getForce(index) for index in range(new_system.getNumForces())} - torsion_force = forces['PeriodicTorsionForce'] - bond_force = forces['HarmonicBondForce'] - atomic_proposal_order = self._generate_proposal_order(new_atoms, new_system, new_to_old_atom_map) - for atomset in atomic_proposal_order: - for atom, possible_torsions in atomset.items(): - #randomly select a torsion from the available list, and calculate its log-probability - selected_torsion_index = np.random.randint(0, len(possible_torsions)) - selected_torsion = possible_torsions[selected_torsion_index] - logp_selected_torsion = - np.log(len(possible_torsions)) - torsion_parameters = torsion_force.getTorsionParameters(selected_torsion) - #check to see whether the atom is the first or last in the torsion - if atom == torsion_parameters[0]: - bonded_atom = torsion_parameters[1] - angle_atom = torsion_parameters[2] - torsion_atom = torsion_parameters[3] - else: - bonded_atom = torsion_parameters[2] - angle_atom = torsion_parameters[1] - torsion_atom = torsion_parameters[0] - - #get a new position for the atom, if specified - if propose_positions: - atomic_xyz, logp_atomic = self._propose_atomic_position(new_system, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters, new_positions) - new_positions[atom] = atomic_xyz - else: - logp_atomic = 0.0 - logp_atomic = self._calculate_logp_atomic_position(new_system, new_positions, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters) - logp += logp_atomic + logp_selected_torsion - - return new_positions, logp - - def _generate_proposal_order(self, new_atoms, new_system, new_to_old_atom_map): - """ - Generate the list of atoms (and corresponding torsions) to be proposed each round - - Arguments - --------- - new_atoms : list of int - Indices of atoms that need positions - topology_proposal : TopologyProposal namedtuple - The result of the topology proposal, containing the atom mapping and topologies - new_system : simtk.OpenMM.System object - The new system - new_positions : [m, 3] np.array of floats - The positions of the m atoms in the new system, with unpositioned atoms having [0,0,0] - - Returns - ------- - atom_proposal_list : list of dict - List of dict of form {atom_index : [torsion_list]}. each entry in the list is a - round of atoms eligible for proposal. - """ - atoms_with_positions = new_to_old_atom_map.keys() - forces = {new_system.getForce(index).__class__.__name__ : new_system.getForce(index) for index in range(new_system.getNumForces())} - torsion_force = forces['PeriodicTorsionForce'] - bond_force = forces['HarmonicBondForce'] - atom_proposal_list = [] - while len(new_atoms) > 0: - atom_torsion_proposals = self._atoms_eligible_for_proposal(torsion_force, bond_force, atoms_with_positions, new_atoms) - #now loop through the list of atoms found to be eligible for positions this round - atom_proposal_list.append(atom_torsion_proposals) - for atom in atom_torsion_proposals.keys(): - #remove atom from list of atoms needing positions - new_atoms.remove(atom) - atoms_with_positions.append(atom) - return atom_proposal_list - - - - def _propose_atomic_position(self, system, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters, positions): - """ - Method to propose the position of a single atom. - - Arguments - --------- - system : simtk.openmm.System - the system containing the atom of interest - atom : int - the index of the atom whose position should be proposed - bonded_atom : int - the index of the atom bonded to the atom of interest - angle_atom : int - the 1, 3 position atom index - torsion_parameters : dict - The parameters of the torsion that is being used for this proposal - - Returns - ------- - atomic_xyz : [1,3] np.array of floats - the cartesian coordinates of the atom - logp_atomic_position : float - the logp of this proposal, including jacobian correction - """ - logging.debug("Proposing position %d with bond atom %d, angle atom %d, torsion atom %d" % (atom, bonded_atom, angle_atom, torsion_atom)) - #get bond parameters and draw bond length - r0, k_eq = self._get_bond_parameters(system, atom, bonded_atom) - r, logp_bond = self._propose_bond_length(r0, k_eq) - - #get angle parameters and draw bond angle - theta0, k_eq = self._get_angle_parameters(system, atom, bonded_atom, angle_atom) - theta, logp_angle = self._propose_bond_angle(theta0, k_eq) - - #propose torsion - phi, logp_torsion = self._propose_torsion_angle(torsion_parameters[6], torsion_parameters[4], torsion_parameters[5]) - #convert spherical to cartesian coordinates - spherical_coordinates = np.asarray([r, theta, phi]) - atomic_xyz, detJ = self._autograd_itoc(bonded_atom, angle_atom, torsion_atom, r, theta, phi, positions) - #accumulate the forward logp with jacobian correction - logp_atomic_position = (logp_angle + logp_bond + logp_torsion + np.log(detJ)) - logging.debug("Proposed (r, theta, phi) of (%s, %s, %s)" % (str(r), str(theta), str(phi))) - return atomic_xyz, logp_atomic_position - - - def _calculate_logp_atomic_position(self, system, positions, atom, bonded_atom, angle_atom, torsion_atom, torsion_parameters): - """ - Calculate the log-probability of a given atom's position - - Arguments - --------- - system : simtk.openmm.System object - The system whose atomic position needs a logp - positions : [n, 3] np.ndarray of floats, Quantity nm - The positions of the particles in the system - atom : int - The index of the atom of interest - bonded_atom : int - The index of the bonded atom - angle_atom : int - The index of the angle atom - torsion_atom : int - The index of the torsion atom - torsion_parameters : array - The parameters of the torsion - - Returns - ------- - logp : float - logp of position, including det(J) correction - """ - #convert cartesian coordinates to spherical - relative_xyz = positions[atom] - positions[bonded_atom] - internal_coordinates, detJ = self._autograd_ctoi(atom, bonded_atom, angle_atom, torsion_atom, positions) - - #get bond parameters and calculate the probability of choosing the current positions - r0, k_eq = self._get_bond_parameters(system, atom, bonded_atom) - sigma_bond = 1.0/np.sqrt(2.0*k_eq/k_eq.unit) - logp_bond = stats.distributions.norm.logpdf(internal_coordinates[0], r0/r0.unit, sigma_bond) - - #calculate the probability of choosing this bond angle - theta0, k_eq = self._get_angle_parameters(system, atom, bonded_atom, angle_atom) - sigma_angle = 1.0/np.sqrt(2.0*k_eq/k_eq.unit) - logp_angle = stats.distributions.norm.logpdf(internal_coordinates[1], theta0/theta0.unit, sigma_angle) - - #calculate the probabilty of choosing this torsion angle - torsion_Z, _ = self._torsion_normalizer(torsion_parameters[6], torsion_parameters[4], torsion_parameters[5]) - logp_torsion = self._torsion_p(torsion_Z,torsion_parameters[6], torsion_parameters[4], torsion_parameters[5], internal_coordinates[2]) - - logp_atomic_position = (logp_angle + logp_bond + logp_torsion + np.log(detJ)) - - return logp_atomic_position - - def _atoms_eligible_for_proposal(self, torsion_force, bond_force, atoms_with_positions, new_atoms): - """ - Determine which atoms are eligible for position proposal this round. - The criteria are that the atom must not yet have a position, and must - be 1 or 4 in at least one torsion where the other atoms have positoins. - - Arguments - --------- - torsion_force : simtk.openmm.PeriodicTorsionForce - The torsion force from the system - bond_roce : simtk.openmm.HarmonicBondForce - The HarmonicBondForce from the system - atoms_with_positions : list of int - A list of the indices of atoms that have positions - - Returns - ------- - atom_torsion_proposals : dict - Dictionary containing {eligible_atom : [list_of_torsions]} - - """ - current_atom_proposals = [] - atom_torsion_proposals = {atom_index: [] for atom_index in new_atoms} - for torsion_index in range(torsion_force.getNumTorsions()): - atom1 = torsion_force.getTorsionParameters(torsion_index)[0] - atom2 = torsion_force.getTorsionParameters(torsion_index)[1] - atom3 = torsion_force.getTorsionParameters(torsion_index)[2] - atom4 = torsion_force.getTorsionParameters(torsion_index)[3] - #Only take torsions where the "new" statuses of atoms 1 and 4 are not equal - if (atom1 in atoms_with_positions) != (atom4 in atoms_with_positions): - #make sure torsion is not improper: - if not self._is_proper(torsion_force.getTorsionParameters(torsion_index), bond_force): - continue - #only take torsions where 2 and 3 have known positions - if (atom2 in atoms_with_positions) and (atom3 in atoms_with_positions): - #finally, append the torsion index to the list of possible atom sets for proposal - if atom1 in atoms_with_positions and (atom1 not in current_atom_proposals): - atom_torsion_proposals[atom4].append(torsion_index) - current_atom_proposals.append(atom4) - elif atom4 in atoms_with_positions and (atom4 not in current_atom_proposals): - atom_torsion_proposals[atom1].append(torsion_index) - current_atom_proposals.append(atom1) - else: - continue - return {atom: torsionlist for atom, torsionlist in atom_torsion_proposals.items() if len(torsionlist) > 0} - - - def _is_proper(self, torsion_parameters, bond_force): - """ - Utility function to determine if torsion is proper - """ - #get the atoms - torsion_atoms = torsion_parameters[:4] - is_proper = True - for i in range(3): - is_proper = is_proper and self._is_bond(torsion_atoms[i], torsion_atoms[i+1], bond_force) - return is_proper - - def _is_bond(self, atom1, atom2, bond_force): - """ - Utility function to determine if bond exists between two atoms - """ - for bond_index in range(bond_force.getNumBonds()): - parameters = bond_force.getBondParameters(bond_index) - if (parameters[0] == atom1 and parameters[1] == atom2) or (parameters[1] == atom1 and parameters[0] == atom2): - return True - return False - - - def _get_bond_parameters(self, system, atom1, atom2): - """ - Utility function to get the bonded parameters for a pair of atoms - - Arguments - --------- - system : simtk.openmm.System - the system containing the parameters of interest - atom1 : int - the index of the first atom in the bond - atom2 : int - the index of the second atom in the bond - - Returns - ------- - r0 : float - equilibrium bond length - k_eq : float - bond spring constant - """ - forces = {system.getForce(index).__class__.__name__ : system.getForce(index) for index in range(system.getNumForces())} - bonded_force = forces['HarmonicBondForce'] - for bond_index in range(bonded_force.getNumBonds()): - parameters = bonded_force.getBondParameters(bond_index) - if (parameters[0] == atom1 and parameters[1] == atom2) or (parameters[1] == atom1 and parameters[0] == atom2): - return parameters[2], parameters[3] - - def _get_angle_parameters(self, system, atom1, atom2, atom3): - """ - Utility function to retrieve harmonic angle parameters for - a given set of atoms in a system - - Arguments - --------- - system : simtk.openmm.System - the system with parameters - atom1 : int - the index of the first atom - atom2 : int - the index of the second atom - atom3 : int - the index of the third atom - - Returns - ------- - theta0 : float - equilibrium bond angle - k_eq : float - angle spring constant - """ - list_of_angle_atoms = [atom1, atom2, atom3] - #compare sorted lists of atoms - forces = {system.getForce(index).__class__.__name__ : system.getForce(index) for index in range(system.getNumForces())} - angle_force = forces['HarmonicAngleForce'] - for angle_index in range(angle_force.getNumAngles()): - parameters = angle_force.getAngleParameters(angle_index) - #the first three "parameters" are atom indices - atoms = parameters[:3] - if np.all(np.equal(list_of_angle_atoms, atoms)) or np.all(np.equal(list_of_angle_atoms[::-1], atoms)): - return parameters[3], parameters[4] - - def _get_unique_atoms(self, new_to_old_map, new_atom_list, old_atom_list): - """ - Get the set of atoms unique to both the new and old system. - - Arguments - --------- - new_to_old_map : dict - Dictionary of the form {new_atom : old_atom} - - Returns - ------- - unique_old_atoms : list of int - A list of the indices of atoms unique to the old system - unique_new_atoms : list of int - A list of the indices of atoms unique to the new system - """ - mapped_new_atoms = new_to_old_map.keys() - unique_new_atoms = [atom for atom in new_atom_list not in mapped_new_atoms] - mapped_old_atoms = new_to_old_map.values() - unique_old_atoms = [atom for atom in old_atom_list not in mapped_old_atoms] - return unique_old_atoms, unique_new_atoms - - def _propose_bond_length(self, r0, k_eq): - """ - Draw a bond length from the equilibrium bond length distribution - - Arguments - --------- - r0 : float - The equilibrium bond length - k_eq : float - The bond spring constant - - Returns - -------- - r : float - the proposed bond length - logp : float - the log-probability of the proposal - """ - sigma = 2.0/np.sqrt(2.0*k_eq/k_eq.unit) - r = sigma*np.random.random()*r0.unit + r0 - logp = stats.distributions.norm.logpdf(r/r.unit, r0/r0.unit, sigma) - return (r, logp) - - def _propose_bond_angle(self, theta0, k_eq): - """ - Draw a bond length from the equilibrium bond length distribution - - Arguments - --------- - theta0 : float - The equilibrium bond angle - k_eq : float - The angle spring constant - - Returns - -------- - r : float - the proposed bond angle - logp : float - the log-probability of the proposal - """ - sigma = 1.0/np.sqrt(2.0*k_eq/k_eq.unit) - theta = sigma*np.random.random()*theta0.unit + theta0 - logp = stats.distributions.norm.logpdf(theta / theta.unit, theta0/theta0.unit, sigma) - return (theta, logp) - - def _propose_torsion_angle(self, V, n, gamma): - """ - Draws a torsion angle from the distribution of one torsion. - Uses the functional form U(phi) = V/2[1+cos(n*phi-gamma)], - as described by Amber. - - Arguments - --------- - V : float - force constant - n : int - multiplicity of the torsion - gamma : float - phase angle - - Returns - ------- - phi : float - The proposed torsion angle - logp : float - The proposed log probability of the torsion angle - """ - (Z, max_p) = self._torsion_normalizer(V, n, gamma) - phi = 0.0 - logp = 0.0 - #sample from the distribution - accepted = False - #use rejection sampling - while not accepted: - phi_samp = np.random.uniform(0.0, 2*np.pi) - runif = np.random.uniform(0.0, max_p+1.0) - p_phi_samp = self._torsion_p(Z, V, n, gamma, phi_samp) - if p_phi_samp > runif: - phi = phi_samp - logp = np.log(p_phi_samp) - accepted = True - else: - continue - return (phi*units.radians, logp) - - def _torsion_p(self, Z, V, n, gamma, phi): - """ - Utility function for calculating the normalized probability of a torsion angle - """ - #must remove units - V = V/V.unit - gamma = gamma/gamma.unit - return (1.0/Z)*np.exp(-(V/2.0)*(1+np.cos(n*phi-gamma))) - - def _torsion_normalizer(self, V, n, gamma): - """ - Utility function to numerically normalize torsion angles. - Also return max_p to facilitate rejection sampling - """ - #generate a grid of 5000 points from 0 < phi < 2*pi - phis = np.linspace(0, 2.0*np.pi, 5000) - #evaluate the unnormalized probability at each of those points - #need to remove units--numexpr can't handle them otherwise - V = V/V.unit - gamma = gamma/gamma.unit - phi_q = ne.evaluate("exp(-(V/2.0)*(1+cos(n*phis-gamma)))") - #integrate the values - Z = np.trapz(phi_q, phis) - max_p = np.max(phi_q/Z) - return Z, max_p - - - - def _autograd_ctoi(self, atom, bond, angle, torsion, positions): - import autograd - import autograd.numpy as np - - positions = positions/positions.unit - atom_position = positions[atom] - - - def _cartesian_to_internal(xyz): - """ - Autograd-based jacobian of transformation from cartesian to internal. - Returns without units! - """ - a = positions[bond] - xyz - b = positions[angle] - positions[bond] - #3-4 bond - c = positions[angle] - positions[torsion] - a_u = a / np.linalg.norm(a) - b_u = b / np.linalg.norm(b) - c_u = c / np.linalg.norm(c) - - #bond length - r = np.linalg.norm(a) - - #bond angle - theta = np.arccos(np.dot(-a_u, b_u)) - - #torsion angle - plane1 = np.cross(a, b) - plane2 = np.cross(b, c) - - phi = np.arccos(np.dot(plane1, plane2)/(np.linalg.norm(plane1)*np.linalg.norm(plane2))) - - if np.dot(np.cross(plane1, plane2), b_u) < 0: - phi = -phi - - return np.array([r, theta, phi]) - - j = autograd.jacobian(_cartesian_to_internal) - internal_coords = _cartesian_to_internal(atom_position) - jacobian_det = np.linalg.det(j(atom_position)) - return internal_coords, np.abs(jacobian_det) - - def _autograd_itoc(self, bond, angle, torsion, r, theta, phi, positions): - """ - Autograd based coordinate conversion internal -> cartesian - - Arguments - --------- - bond : int - index of the bonded atom - angle : int - index of the angle atom - torsion : int - index of the torsion atom - r : float, Quantity nm - bond length - theta : float, Quantity rad - bond angle - phi : float, Quantity rad - torsion angle - positions : [n, 3] np.array Quantity nm - positions of the atoms in the molecule - - Returns - ------- - atom_xyz : [1, 3] np.array Quantity nm - The atomic positions in cartesian space - detJ : float - The absolute value of the determinant of the jacobian - """ - import autograd - import autograd.numpy as np - positions = positions/positions.unit - r = r/r.unit - theta = theta/theta.unit - phi = phi/phi.unit - rtp = np.array([r, theta, phi]) - - def _internal_to_cartesian(rthetaphi): - - a = positions[angle] - positions[bond] - b = positions[angle] - positions[torsion] - - a_u = a / np.linalg.norm(a) - b_u = b / np.linalg.norm(b) - - - d_r = rthetaphi[0]*a_u - - normal = np.cross(a, b) - - #construct the angle rotation matrix - axis_angle = normal / np.linalg.norm(normal) - a = np.cos(rthetaphi[1]/2) - b, c, d = -axis_angle*np.sin(rthetaphi[1]/2) - angle_rotation_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], - [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], - [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) - #apply it - d_ang = np.dot(angle_rotation_matrix, d_r) - - #construct the torsion rotation matrix and apply it - axis = a_u - a = np.cos(rthetaphi[2]/2) - b, c, d = -axis*np.sin(rthetaphi[2]/2) - torsion_rotation_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], - [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], - [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) - #apply it - d_torsion = np.dot(torsion_rotation_matrix, d_ang) - - #add the positions of the bond atom - xyz = positions[bond] + d_torsion - - return xyz - - j = autograd.jacobian(_internal_to_cartesian) - atom_xyz = _internal_to_cartesian(rtp) - jacobian_det = np.linalg.det(j(rtp)) - logging.debug("detJ is %f" %(jacobian_det)) - detj_spherical = r**2*np.sin(theta) - logging.debug("The spherical detJ is %f" % detj_spherical) - return units.Quantity(atom_xyz, unit=units.nanometers), np.abs(jacobian_det) - From a1b3e69274e0fcfd20bc1f2a9c76488696bc5411 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 17:42:42 -0500 Subject: [PATCH 04/57] add 1/2 prefactor --- perses/rjmc/geometry.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 1fa917604..654554d5d 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -193,14 +193,15 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) - logp_r = self._bond_logp(internal_coordinates[0], bond, top_proposal.beta) + bond_logz = np.sqrt(2.0*np.pi)* + logp_r = self._bond_logq(internal_coordinates[0], bond, top_proposal.beta) #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - logp_theta = self._angle_logp(internal_coordinates[1], angle, top_proposal.beta) + logp_theta = self._angle_logq(internal_coordinates[1], angle, top_proposal.beta) #calculate torsion probability - logp_phi = self._torsion_logp(atom, atom_coords, torsion, atoms_with_positions, reverse_proposal_coordinates, beta) + logp_phi = self._torsion_logq(atom, atom_coords, torsion, atoms_with_positions, reverse_proposal_coordinates, beta) logp = logp + logp_choice + logp_r + logp_theta + logp_phi + np.log(detJ) atoms_with_positions.append(atom) @@ -455,7 +456,7 @@ def _bond_logq(self, r, bond, beta): """ k_eq = bond.type.k*units.kilocalories_per_mole/(units.angstrom**2) r0 = bond.type.req*units.nanometers - logq = -beta*k_eq*(r-r0)**2 + logq = -beta*0.5*k_eq*(r-r0)**2 return logq def _angle_logq(self, theta, angle, beta): @@ -473,7 +474,7 @@ def _angle_logq(self, theta, angle, beta): """ k_eq = angle.type.k*units.kilocalories_per_mole/(units.radians**2) theta0 = angle.type.theteq*units.radians - logq = -beta*k_eq*(theta-theta0)**2 + logq = -beta*k_eq*0.5*(theta-theta0)**2 return logq def _torsion_logq(self, atom, xyz, torsion, atoms_with_positions, positions, beta): @@ -576,7 +577,7 @@ def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, in bond_atom_position = xyz if angle.atom2 == atom else positions[angle.atom2.idx] angle_atom_position = xyz if angle.atom3 == atom else positions[angle.atom3.idx] theta = self._calculate_angle(atom_position, bond_atom_position, angle_atom_position) - ub_angles += self._angle_logp(theta*units.radians, angle, beta) + ub_angles += self._angle_logq(theta*units.radians, angle, beta) for torsion in involved_torsions: atom_position = xyz if torsion.atom1 == atom else positions[torsion.atom1.idx] bond_atom_position = xyz if torsion.atom2 == atom else positions[torsion.atom2.idx] @@ -630,7 +631,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor ub[i]+=ub_i #exponentiate to get the unnormalized probability - q = np.exp(-ub) + q = np.exp(ub) #estimate the normalizing constant Z = np.trapz(q, phis) @@ -684,4 +685,3 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, logp = np.log(p[phi_idx]) phi = phis[phi_idx] return phi, logp - From 24ed133ff179ae479a2ab6937b1920c6ea9b8f69 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 18:58:28 -0500 Subject: [PATCH 05/57] calculate logZ separately --- perses/rjmc/geometry.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 654554d5d..9a94422c0 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -193,12 +193,15 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) - bond_logz = np.sqrt(2.0*np.pi)* - logp_r = self._bond_logq(internal_coordinates[0], bond, top_proposal.beta) + sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) + logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r)) + logp_r = self._bond_logq(internal_coordinates[0], bond, top_proposal.beta) - logZ_r #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - logp_theta = self._angle_logq(internal_coordinates[1], angle, top_proposal.beta) + sigma_theta = np.sqrt(1.0/(beta/beta.unit*angle.type.k)) + logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta)) + logp_theta = self._angle_logq(internal_coordinates[1], angle, top_proposal.beta) - logZ_theta #calculate torsion probability logp_phi = self._torsion_logq(atom, atom_coords, torsion, atoms_with_positions, reverse_proposal_coordinates, beta) From 1bb11e49509816151e07cb9f845c8a30d1a8016b Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 19:07:19 -0500 Subject: [PATCH 06/57] fixed bond and angle proposals --- perses/rjmc/geometry.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 9a94422c0..8c9f9d3b8 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -90,6 +90,7 @@ def propose(self, top_proposal): logp_proposal : float The log probability of the forward-only proposal """ + beta = top_proposa.beta logp_proposal = 0.0 structure = parmed.openmm.load_topology(top_proposal.new_topology, top_proposal.new_system) new_atoms = [structure.atoms[idx] for idx in top_proposal.unique_new_atoms] @@ -117,16 +118,18 @@ def propose(self, top_proposal): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) - r_proposed = self._propose_bond(bond, top_proposal.beta) - logp_r = self._bond_logp(r_proposed, bond, top_proposal.beta) + r_proposed = self._propose_bond(bond, beta) + sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) + logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r)) + logp_r = self._bond_logq(r_proposed, bond, beta) #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - theta_proposed = self._propose_angle(angle, top_proposal.beta) - logp_theta = self._angle_logp(theta_proposed, angle, top_proposal.beta) + theta_proposed = self._propose_angle(angle, beta) + logp_theta = self._angle_logp(theta_proposed, angle, beta) #propose a torsion angle and calcualate its probability - phi_proposed, logp_phi = self._propose_torsion(atom, r_proposed, theta_proposed, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, top_proposal.beta) + phi_proposed, logp_phi = self._propose_torsion(atom, r_proposed, theta_proposed, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, beta) #convert to cartesian xyz, detJ = self._autograd_itoc(bond_atom.idx, angle_atom.idx, torsion_atom.idx, r_proposed, theta_proposed, phi_proposed, new_positions) @@ -204,7 +207,7 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): logp_theta = self._angle_logq(internal_coordinates[1], angle, top_proposal.beta) - logZ_theta #calculate torsion probability - logp_phi = self._torsion_logq(atom, atom_coords, torsion, atoms_with_positions, reverse_proposal_coordinates, beta) + logp_phi = self._torsion_logp(atom, atom_coords, torsion, atoms_with_positions, reverse_proposal_coordinates, beta) logp = logp + logp_choice + logp_r + logp_theta + logp_phi + np.log(detJ) atoms_with_positions.append(atom) @@ -527,10 +530,9 @@ def _propose_bond(self, bond, beta): """ Bond length proposal """ - k_eq = bond.type.k*units.kilocalories_per_mole/(units.angstrom**2) r0 = bond.type.req*units.angstrom - sigma = beta*2.0/np.sqrt(2.0*k_eq/k_eq.unit) - r = sigma/sigma.unit*np.random.random()*r0.unit + r0 + sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) + r = sigma_r/sigma_r.unit*np.random.random()*r0.unit + r0 return r def _propose_angle(self, angle, beta): @@ -538,9 +540,8 @@ def _propose_angle(self, angle, beta): Bond angle proposal """ theta0 = angle.type.theteq*units.radians - k_eq = angle.type.k*units.kilocalories_per_mole/(units.radians**2) - sigma = beta*2.0/np.sqrt(2.0*k_eq/k_eq.unit) - theta = sigma/sigma.unit*np.random.random()*theta0.unit + theta0 + sigma_theta = np.sqrt(1.0/(beta/beta.unit*angle.type.k)) + theta = sigma_theta/sigma_theta.unit*np.random.random()*theta0.unit + theta0 return theta def _torsion_potential(self, torsion, phi, beta): @@ -556,6 +557,9 @@ def _torsion_potential(self, torsion, phi, beta): q = -beta*(V/2.0)*(1+np.cos(n*phi-gamma)) return q + def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta, Z=None): + pass + def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, positions, beta): pass From 48ffa76f450481886b758a4b6d79b23a2de8934f Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 19:09:45 -0500 Subject: [PATCH 07/57] make sure beta is in the right units --- perses/rjmc/geometry.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 8c9f9d3b8..7fa87ce9d 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -90,7 +90,7 @@ def propose(self, top_proposal): logp_proposal : float The log probability of the forward-only proposal """ - beta = top_proposa.beta + beta = top_proposal.beta.in_units_of(units.kilocalorie_per_mole**-1) logp_proposal = 0.0 structure = parmed.openmm.load_topology(top_proposal.new_topology, top_proposal.new_system) new_atoms = [structure.atoms[idx] for idx in top_proposal.unique_new_atoms] @@ -160,6 +160,7 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): logp : float The log probability of the proposal for the given transformation """ + beta = top_proposal.beta.in_units_of(units.kilocalorie_per_mole**-1) logp = 0.0 top_proposal = topology_proposal.SmallMoleculeTopologyProposal() structure = parmed.openmm.load_topology(top_proposal.old_topology, top_proposal.old_system) @@ -167,7 +168,6 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): new_atoms = [structure.atoms[idx] for idx in top_proposal.unique_old_atoms] #we'll need to copy the current positions of the core to the old system #In the case of C-A --> C/-A -> C/-B ---> C-B these are the same - beta = top_proposal.beta reverse_proposal_coordinates = units.Quantity(np.zeros([top_proposal.n_atoms_new, 3]), unit=units.nanometers) for atom in atoms_with_positions: new_index = top_proposal.old_to_new_atom_map[atom.idx] @@ -198,13 +198,13 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): bond = self._get_relevant_bond(atom, bond_atom) sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r)) - logp_r = self._bond_logq(internal_coordinates[0], bond, top_proposal.beta) - logZ_r + logp_r = self._bond_logq(internal_coordinates[0], bond, beta) - logZ_r #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) sigma_theta = np.sqrt(1.0/(beta/beta.unit*angle.type.k)) logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta)) - logp_theta = self._angle_logq(internal_coordinates[1], angle, top_proposal.beta) - logZ_theta + logp_theta = self._angle_logq(internal_coordinates[1], angle, beta) - logZ_theta #calculate torsion probability logp_phi = self._torsion_logp(atom, atom_coords, torsion, atoms_with_positions, reverse_proposal_coordinates, beta) From 38ac293cfc5754e04c7ff17f2d93bd74bcdbec46 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 20:54:33 -0500 Subject: [PATCH 08/57] fixed angle thing --- perses/rjmc/geometry.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 7fa87ce9d..9163784b2 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -119,14 +119,18 @@ def propose(self, top_proposal): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) r_proposed = self._propose_bond(bond, beta) - sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) - logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r)) - logp_r = self._bond_logq(r_proposed, bond, beta) + bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 + 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 #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) theta_proposed = self._propose_angle(angle, beta) - logp_theta = self._angle_logp(theta_proposed, angle, beta) + angle_k = angle.type.k*units.kilocalorie_per_mole/units.angstrom**2 + sigma_theta = units.sqrt(1/beta*angle_k) + logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) + logp_theta = self._angle_logq(theta_proposed, angle, beta) - logZ_theta #propose a torsion angle and calcualate its probability phi_proposed, logp_phi = self._propose_torsion(atom, r_proposed, theta_proposed, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, beta) @@ -196,14 +200,16 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) - sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) - logZ_r = np.log((np.sqrt(2*np.pi)*sigma_r)) + bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 + 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(internal_coordinates[0], bond, beta) - logZ_r #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - sigma_theta = np.sqrt(1.0/(beta/beta.unit*angle.type.k)) - logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta)) + angle_k = angle.type.k*units.kilocalorie_per_mole/(units.radians**2) + sigma_theta = units.sqrt(1/beta*angle_k) + logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) #need to eliminate unit to allow numpy log logp_theta = self._angle_logq(internal_coordinates[1], angle, beta) - logZ_theta #calculate torsion probability From 1344f99393aa0b04351660189ff0acf909d271fd Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:00:43 -0500 Subject: [PATCH 09/57] removed unused method --- perses/rjmc/geometry.py | 28 +++------------------------- 1 file changed, 3 insertions(+), 25 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 9163784b2..e88e183ce 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -489,27 +489,6 @@ def _angle_logq(self, theta, angle, beta): logq = -beta*k_eq*0.5*(theta-theta0)**2 return logq - def _torsion_logq(self, atom, xyz, torsion, atoms_with_positions, positions, beta): - """ - Utility function for calculating the unnormalized probability of a torsion angle - """ - if torsion.atom1 == atom: - bond_atom = torsion.atom2 - angle_atom = torsion.atom3 - torsion_atom = torsion.atom4 - else: - bond_atom = torsion.atom3 - angle_atom = torsion.atom2 - torsion_atom = torsion.atom1 - - internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) - beta = beta/beta.unit - gamma = torsion.type.phase - V = torsion.type.phi_k - n = torsion.type.per - logq = -beta*(V/2.0)*(1+np.cos(n*internal_coordinates[2]-gamma)) - return logq - def _choose_torsion(self, atoms_with_positions, atom_for_proposal): """ Pick an eligible torsion uniformly @@ -563,7 +542,7 @@ def _torsion_potential(self, torsion, phi, beta): q = -beta*(V/2.0)*(1+np.cos(n*phi-gamma)) return q - def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta, Z=None): + def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta): pass def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, positions, beta): @@ -665,7 +644,7 @@ def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_positio theta = np.arccos(np.dot(a_u, b_u)) return theta - def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta, Z=None): + def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta): """ Calculate the log-probability of a given torsion. This is calculated via a distribution that includes angle and other torsion potentials. @@ -681,8 +660,7 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet involved_angles = self._get_valid_angles(atoms_with_positions, atom) involved_torsions = self._get_torsions(atoms_with_positions, atom) internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) - if not Z: - p, Z, _, _ = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + p, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) ub_torsion = self._torsion_and_angle_potential(xyz, atom, positions, involved_angles, involved_torsions, beta) p_torsion = np.exp(-ub_torsion) / Z return p_torsion From d2adbb6748e6d3a351cdf415448e31761b9782f8 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:09:34 -0500 Subject: [PATCH 10/57] calculate reverse torsion_logp --- perses/rjmc/geometry.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index e88e183ce..c650665b6 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -661,9 +661,10 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet involved_torsions = self._get_torsions(atoms_with_positions, atom) internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) p, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) - ub_torsion = self._torsion_and_angle_potential(xyz, atom, positions, involved_angles, involved_torsions, beta) - p_torsion = np.exp(-ub_torsion) / Z - return p_torsion + #find the phi that's closest to the internal_coordinate phi: + phi_idx, phi = min(enumerate(phis), key=lambda x: abs(x[1]-internal_coordinates[2])) + p_phi = p[phi_idx] + return np.log(p_phi) def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, positions, beta): """ From bd08b30c9d434655eb5c58f2eb2e794788771c39 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:13:49 -0500 Subject: [PATCH 11/57] removed unnecessary imports --- perses/rjmc/geometry.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index c650665b6..aa066649b 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -5,11 +5,10 @@ from collections import namedtuple import perses.rjmc.topology_proposal as topology_proposal import numpy as np -import numexpr as ne import parmed import simtk.unit as units import logging -from numba import jit + GeometryProposal = namedtuple('GeometryProposal',['new_positions','logp']) From 25aedf77258a8ffc49a5ff8733c2bca679e8667d Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:18:46 -0500 Subject: [PATCH 12/57] fixed some unit stuff --- perses/rjmc/geometry.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index aa066649b..bbad5b987 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -515,8 +515,9 @@ def _propose_bond(self, bond, beta): Bond length proposal """ r0 = bond.type.req*units.angstrom - sigma_r = np.sqrt(1.0/(beta/beta.unit*bond.type.k)) - r = sigma_r/sigma_r.unit*np.random.random()*r0.unit + r0 + k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 + sigma_r = units.sqrt(1.0/(beta*k)) + r = sigma_r*np.random.random() + r0 return r def _propose_angle(self, angle, beta): @@ -524,8 +525,9 @@ def _propose_angle(self, angle, beta): Bond angle proposal """ theta0 = angle.type.theteq*units.radians - sigma_theta = np.sqrt(1.0/(beta/beta.unit*angle.type.k)) - theta = sigma_theta/sigma_theta.unit*np.random.random()*theta0.unit + theta0 + k = angle.type.k*units.kilocalorie_per_mole/units.radian**2 + sigma_theta = units.sqrt(1.0/(beta*k)) + theta = sigma_theta*np.random.random() + theta0 return theta def _torsion_potential(self, torsion, phi, beta): From 9346fff4dbe2196b59578fb426467c8f5013b8ec Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:27:06 -0500 Subject: [PATCH 13/57] remove loop --- perses/tests/test_geometry_engine.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index d4a116eac..98be99cb9 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -98,5 +98,4 @@ def test_run_geometry_engine(): if __name__=="__main__": - for i in range(10): - test_run_geometry_engine() \ No newline at end of file + test_run_geometry_engine() \ No newline at end of file From 52aa9d2a3cfe2785adf5e5a622fa7817819fc575 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:28:41 -0500 Subject: [PATCH 14/57] remove another loop --- perses/tests/test_geometry_engine.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 98be99cb9..cb8ee846f 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -89,11 +89,10 @@ def test_run_geometry_engine(): state = context.getState(getEnergy=True) print("Energy before proposal is: %s" % str(state.getPotentialEnergy())) - for i in range(10): - new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal) - context.setPositions(new_positions) - state2 = context.getState(getEnergy=True) - print("Energy after proposal is: %s" %str(state2.getPotentialEnergy())) + new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal) + context.setPositions(new_positions) + state2 = context.getState(getEnergy=True) + print("Energy after proposal is: %s" %str(state2.getPotentialEnergy())) From cf926735d2af1df6389626a9f67c7bda2a23cb28 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 21:35:16 -0500 Subject: [PATCH 15/57] make jacobian optional again --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index bbad5b987..1a655ff23 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -576,7 +576,7 @@ def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, in bond_atom_position = xyz if torsion.atom2 == atom else positions[torsion.atom2.idx] angle_atom_position = xyz if torsion.atom3 == atom else positions[torsion.atom3.idx] torsion_atom_position = xyz if torsion.atom4 == atom else positions[torsion.atom4.idx] - internal_coordinates, _ = self._autograd_ctoi(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position) + internal_coordinates, _ = self._autograd_ctoi(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position, calculate_jacobian=False) phi = internal_coordinates[2]*units.radians ub_torsions += self._torsion_potential(torsion, phi, beta) return ub_angles+ub_torsions From ec3e2187a1d60936d3605fe0b01c765317661814 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 1 Dec 2015 22:01:36 -0500 Subject: [PATCH 16/57] prevent integer division --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 1a655ff23..ebad3cca4 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -495,7 +495,7 @@ def _choose_torsion(self, atoms_with_positions, atom_for_proposal): eligible_torsions = self._get_torsions(atoms_with_positions, atom_for_proposal) torsion_idx = np.random.randint(0, len(eligible_torsions)) torsion_selected = eligible_torsions[torsion_idx] - return torsion_selected, np.log(1/len(eligible_torsions)) + return torsion_selected, np.log(1.0/len(eligible_torsions)) def _atoms_eligible_for_proposal(self, new_atoms, structure, atoms_with_positions): """ From d8ee10e1e7fb956260996a9823852687c3e35aa6 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 13:43:18 -0500 Subject: [PATCH 17/57] make some fixes --- perses/rjmc/geometry.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index ebad3cca4..50614d4d9 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -200,14 +200,14 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 - sigma_r = units.sqrt(1/beta*bond_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(internal_coordinates[0], bond, beta) - logZ_r #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - angle_k = angle.type.k*units.kilocalorie_per_mole/(units.radians**2) - sigma_theta = units.sqrt(1/beta*angle_k) + angle_k = angle.type.k*units.kilocalorie_per_mole/(units.degrees**2) + sigma_theta = units.sqrt(1/(beta*angle_k)) logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) #need to eliminate unit to allow numpy log logp_theta = self._angle_logq(internal_coordinates[1], angle, beta) - logZ_theta @@ -483,8 +483,8 @@ def _angle_logq(self, theta, angle, beta): beta : float 1/kT or inverse temperature """ - k_eq = angle.type.k*units.kilocalories_per_mole/(units.radians**2) - theta0 = angle.type.theteq*units.radians + k_eq = angle.type.k*units.kilocalories_per_mole/(units.degrees**2) + theta0 = angle.type.theteq*units.degrees logq = -beta*k_eq*0.5*(theta-theta0)**2 return logq @@ -524,7 +524,7 @@ def _propose_angle(self, angle, beta): """ Bond angle proposal """ - theta0 = angle.type.theteq*units.radians + theta0 = angle.type.theteq*units.degrees k = angle.type.k*units.kilocalorie_per_mole/units.radian**2 sigma_theta = units.sqrt(1.0/(beta*k)) theta = sigma_theta*np.random.random() + theta0 @@ -616,15 +616,15 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor xyzs[i], _ = self._autograd_itoc(bond_atom.idx, angle_atom.idx, torsion_atom.idx, r, theta, phi, positions, calculate_jacobian=False) #set up arrays for energies from angles and torsions - ub= np.zeros(n_divisions) + loq= np.zeros(n_divisions) for i, xyz in enumerate(xyzs): ub_i = self._torsion_and_angle_potential(xyz, atom, positions, involved_angles, involved_torsions, beta) if np.isnan(ub_i): ub_i = np.inf - ub[i]+=ub_i + loq[i]+=ub_i #exponentiate to get the unnormalized probability - q = np.exp(ub) + q = np.exp(loq) #estimate the normalizing constant Z = np.trapz(q, phis) From 9aa04c23b56af7361eddac87d8523b64f380292b Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 13:46:02 -0500 Subject: [PATCH 18/57] rename --- perses/rjmc/geometry.py | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 50614d4d9..27206ac76 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -530,7 +530,7 @@ def _propose_angle(self, angle, beta): theta = sigma_theta*np.random.random() + theta0 return theta - def _torsion_potential(self, torsion, phi, beta): + def _torsion_logq(self, torsion, phi, beta): """ Calculate the log-unnormalized probability of the torsion @@ -540,8 +540,8 @@ def _torsion_potential(self, torsion, phi, beta): gamma = torsion.type.phase V = torsion.type.phi_k n = torsion.type.per - q = -beta*(V/2.0)*(1+np.cos(n*phi-gamma)) - return q + logq = -beta*(V/2.0)*(1+np.cos(n*phi-gamma)) + return logq def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta): pass @@ -556,13 +556,13 @@ class FFAllAngleGeometryEngine(FFGeometryEngine): and torsion_p methods of the base. """ - def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, involved_torsions, beta): + def _torsion_and_angle_logq(self, xyz, atom, positions, involved_angles, involved_torsions, beta): """ Calculate the potential resulting from torsions and angles at a given cartesian coordinate """ - ub_angles = 0.0 - ub_torsions = 0.0 + logq_angles = 0.0 + logq_torsions= 0.0 if type(xyz) != units.Quantity: xyz = units.Quantity(xyz, units.nanometers) for angle in involved_angles: @@ -570,7 +570,7 @@ def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, in bond_atom_position = xyz if angle.atom2 == atom else positions[angle.atom2.idx] angle_atom_position = xyz if angle.atom3 == atom else positions[angle.atom3.idx] theta = self._calculate_angle(atom_position, bond_atom_position, angle_atom_position) - ub_angles += self._angle_logq(theta*units.radians, angle, beta) + logq_angles += self._angle_logq(theta*units.radians, angle, beta) for torsion in involved_torsions: atom_position = xyz if torsion.atom1 == atom else positions[torsion.atom1.idx] bond_atom_position = xyz if torsion.atom2 == atom else positions[torsion.atom2.idx] @@ -578,8 +578,8 @@ def _torsion_and_angle_potential(self, xyz, atom, positions, involved_angles, in torsion_atom_position = xyz if torsion.atom4 == atom else positions[torsion.atom4.idx] internal_coordinates, _ = self._autograd_ctoi(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position, calculate_jacobian=False) phi = internal_coordinates[2]*units.radians - ub_torsions += self._torsion_potential(torsion, phi, beta) - return ub_angles+ub_torsions + logq_torsions += self._torsion_logq(torsion, phi, beta) + return logq_angles+logq_torsions def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000): @@ -616,15 +616,15 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor xyzs[i], _ = self._autograd_itoc(bond_atom.idx, angle_atom.idx, torsion_atom.idx, r, theta, phi, positions, calculate_jacobian=False) #set up arrays for energies from angles and torsions - loq= np.zeros(n_divisions) + logq = np.zeros(n_divisions) for i, xyz in enumerate(xyzs): - ub_i = self._torsion_and_angle_potential(xyz, atom, positions, involved_angles, involved_torsions, beta) - if np.isnan(ub_i): - ub_i = np.inf - loq[i]+=ub_i + logq_i = self._torsion_and_angle_logq(xyz, atom, positions, involved_angles, involved_torsions, beta) + if np.isnan(logq_i): + logq_i = np.inf + logq[i]+=logq_i #exponentiate to get the unnormalized probability - q = np.exp(loq) + q = np.exp(logq) #estimate the normalizing constant Z = np.trapz(q, phis) @@ -632,7 +632,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor #get the normalized probabilities for torsions p = q / Z - return p, Z, q, phis + return p, Z, q, phis, logq def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_position): """ From 6542577371ba08b971765aa06755f140dc2072a3 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 13:50:12 -0500 Subject: [PATCH 19/57] cleaning up additional errors --- perses/rjmc/geometry.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 27206ac76..c51274a17 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -658,10 +658,8 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet bond_atom = torsion.atom3 angle_atom = torsion.atom2 torsion_atom = torsion.atom1 - involved_angles = self._get_valid_angles(atoms_with_positions, atom) - involved_torsions = self._get_torsions(atoms_with_positions, atom) internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) - p, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + p, Z, q, phis, logq = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #find the phi that's closest to the internal_coordinate phi: phi_idx, phi = min(enumerate(phis), key=lambda x: abs(x[1]-internal_coordinates[2])) p_phi = p[phi_idx] @@ -672,7 +670,7 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, Propose a torsion angle, including energetic contributions from other torsions and angles """ #first, let's get the normalizing constant of this distribution - p, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + p, Z, q, phis, logq = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #choose from the set of possible torsion angles phi_idx = np.random.choice(range(len(phis)), p=p) logp = np.log(p[phi_idx]) From a8072603425bcf223bfba4dc54599619d2a7d9d0 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 13:55:51 -0500 Subject: [PATCH 20/57] added some explicit units --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index c51274a17..44fa057df 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -517,7 +517,7 @@ def _propose_bond(self, bond, beta): r0 = bond.type.req*units.angstrom k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 sigma_r = units.sqrt(1.0/(beta*k)) - r = sigma_r*np.random.random() + r0 + r = sigma_r/sigma_r.unit*np.random.random()*units.angstrom + r0 return r def _propose_angle(self, angle, beta): @@ -527,7 +527,7 @@ def _propose_angle(self, angle, beta): theta0 = angle.type.theteq*units.degrees k = angle.type.k*units.kilocalorie_per_mole/units.radian**2 sigma_theta = units.sqrt(1.0/(beta*k)) - theta = sigma_theta*np.random.random() + theta0 + theta = sigma_theta/sigma_theta.unit*np.random.random()*units.radian + theta0 return theta def _torsion_logq(self, torsion, phi, beta): From 6a6fabcdf5584d41029a068f6516efbfa6ad78a4 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 14:00:17 -0500 Subject: [PATCH 21/57] subtracted max from logq --- perses/rjmc/geometry.py | 1 + 1 file changed, 1 insertion(+) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 44fa057df..095c865ce 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -622,6 +622,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor if np.isnan(logq_i): logq_i = np.inf logq[i]+=logq_i + logq - max(logq) #exponentiate to get the unnormalized probability q = np.exp(logq) From a741bf356171d53ec301ba3d0148b6049cd8a24f Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 14:32:42 -0500 Subject: [PATCH 22/57] fix denominator --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 095c865ce..d10c7092c 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -119,7 +119,7 @@ def propose(self, top_proposal): bond = self._get_relevant_bond(atom, bond_atom) r_proposed = self._propose_bond(bond, beta) bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 - sigma_r = units.sqrt(1/beta*bond_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 @@ -127,7 +127,7 @@ def propose(self, top_proposal): angle = self._get_relevant_angle(atom, bond_atom, angle_atom) theta_proposed = self._propose_angle(angle, beta) angle_k = angle.type.k*units.kilocalorie_per_mole/units.angstrom**2 - sigma_theta = units.sqrt(1/beta*angle_k) + sigma_theta = units.sqrt(1/(beta*angle_k)) logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) logp_theta = self._angle_logq(theta_proposed, angle, beta) - logZ_theta From 98e5673bf96e87b9a9c15cc162ae0a83236162a6 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 14:33:49 -0500 Subject: [PATCH 23/57] fix units --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index d10c7092c..ae09a8a57 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -126,7 +126,7 @@ def propose(self, top_proposal): #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) theta_proposed = self._propose_angle(angle, beta) - angle_k = angle.type.k*units.kilocalorie_per_mole/units.angstrom**2 + angle_k = angle.type.k*units.kilocalorie_per_mole/units.degree**2 sigma_theta = units.sqrt(1/(beta*angle_k)) logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) logp_theta = self._angle_logq(theta_proposed, angle, beta) - logZ_theta @@ -525,7 +525,7 @@ def _propose_angle(self, angle, beta): Bond angle proposal """ theta0 = angle.type.theteq*units.degrees - k = angle.type.k*units.kilocalorie_per_mole/units.radian**2 + k = angle.type.k*units.kilocalorie_per_mole/units.degree**2 sigma_theta = units.sqrt(1.0/(beta*k)) theta = sigma_theta/sigma_theta.unit*np.random.random()*units.radian + theta0 return theta From cc3257948d699af0a34a9e0e08e99c62f21ccd43 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 20:17:01 -0500 Subject: [PATCH 24/57] create functions to assign units properly --- perses/rjmc/geometry.py | 47 +++++++++++++++++++++++++++++++++++++---- 1 file changed, 43 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index ae09a8a57..fb72a4773 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -237,8 +237,10 @@ def _get_relevant_bond(self, atom1, atom2): """ bonds_1 = set(atom1.bonds) bonds_2 = set(atom2.bonds) - relevant_bond = bonds_1.intersection(bonds_2) - return relevant_bond.pop() + relevant_bond_set = bonds_1.intersection(bonds_2) + relevant_bond = relevant_bond_set.pop() + relevant_bond_with_units = self._add_bond_units(relevant_bond) + return relevant_bond_with_units def _get_relevant_angle(self, atom1, atom2, atom3): """ @@ -247,8 +249,45 @@ def _get_relevant_angle(self, atom1, atom2, atom3): atom1_angles = set(atom1.angles) atom2_angles = set(atom2.angles) atom3_angles = set(atom3.angles) - relevant_angle = atom1_angles.intersection(atom2_angles, atom3_angles) - return relevant_angle.pop() + relevant_angle_set = atom1_angles.intersection(atom2_angles, atom3_angles) + relevant_angle = relevant_angle_set.pop() + relevant_angle_with_units = self._add_angle_units(relevant_angle) + return relevant_angle_with_units + + def _add_bond_units(self, bond): + """ + Add the correct units to a harmonic bond + + Arguments + --------- + bond : parmed bond object + The bond to get units + + Returns + ------- + + """ + bond.type.req = units.Quantity(bond.type.req, unit=units.angstrom) + bond.type.k = units.Quantity(bond.type.k, unit=units.kilocalorie_per_mole/units.angstrom**2) + return bond + + def _add_angle_units(self, angle): + """ + Add the correct units to a harmonic angle + + Arguments + ---------- + angle : parmed angle object + the angle to get unit-ed + + Returns + ------- + angle_with_units : parmed angle + The angle, but with units on its parameters + """ + angle.type.theteq = units.Quantity(angle.type.theteq, unit=units.degree) + angle.type.k = units.Quantity(angle.type.k, unit=units.kilocalorie_per_mole/units.radian**2) + return angle def _get_torsions(self, atoms_with_positions, new_atom): """ From 90329101754ffe5da8d0208194561630e20bc9ac Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 22:16:30 -0500 Subject: [PATCH 25/57] continue fixing unit errors --- perses/rjmc/geometry.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index fb72a4773..a9cba74f5 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -118,7 +118,7 @@ def propose(self, top_proposal): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) r_proposed = self._propose_bond(bond, beta) - bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 + 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 @@ -126,7 +126,7 @@ def propose(self, top_proposal): #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) theta_proposed = self._propose_angle(angle, beta) - angle_k = angle.type.k*units.kilocalorie_per_mole/units.degree**2 + angle_k = angle.type.k sigma_theta = units.sqrt(1/(beta*angle_k)) logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) logp_theta = self._angle_logq(theta_proposed, angle, beta) - logZ_theta @@ -553,8 +553,8 @@ def _propose_bond(self, bond, beta): """ Bond length proposal """ - r0 = bond.type.req*units.angstrom - k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 + r0 = bond.type.req + k = bond.type.k sigma_r = units.sqrt(1.0/(beta*k)) r = sigma_r/sigma_r.unit*np.random.random()*units.angstrom + r0 return r @@ -563,8 +563,8 @@ def _propose_angle(self, angle, beta): """ Bond angle proposal """ - theta0 = angle.type.theteq*units.degrees - k = angle.type.k*units.kilocalorie_per_mole/units.degree**2 + theta0 = angle.type.theteq + k = angle.type.k sigma_theta = units.sqrt(1.0/(beta*k)) theta = sigma_theta/sigma_theta.unit*np.random.random()*units.radian + theta0 return theta From 63fae744620506a6cb52786c5264479735b7f7ee Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Wed, 2 Dec 2015 22:18:06 -0500 Subject: [PATCH 26/57] fix more units... --- perses/rjmc/geometry.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index a9cba74f5..dbd2f8b71 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -199,14 +199,14 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) - bond_k = bond.type.k*units.kilocalorie_per_mole/units.angstrom**2 + 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(internal_coordinates[0], bond, beta) - logZ_r #propose an angle and calculate its probability angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - angle_k = angle.type.k*units.kilocalorie_per_mole/(units.degrees**2) + angle_k = angle.type.k sigma_theta = units.sqrt(1/(beta*angle_k)) logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) #need to eliminate unit to allow numpy log logp_theta = self._angle_logq(internal_coordinates[1], angle, beta) - logZ_theta @@ -579,7 +579,7 @@ def _torsion_logq(self, torsion, phi, beta): gamma = torsion.type.phase V = torsion.type.phi_k n = torsion.type.per - logq = -beta*(V/2.0)*(1+np.cos(n*phi-gamma)) + logq = -beta*(V/2.0)*(1+units.cos(n*phi-gamma)) return logq def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta): From 045a1cbfaa41126a00f3a9620dbe67241d5c865a Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 11:27:24 -0500 Subject: [PATCH 27/57] removed superfluous units --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index dbd2f8b71..1c9701ad1 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -504,8 +504,8 @@ def _bond_logq(self, r, bond, beta): beta : float 1/kT or inverse temperature """ - k_eq = bond.type.k*units.kilocalories_per_mole/(units.angstrom**2) - r0 = bond.type.req*units.nanometers + k_eq = bond.type.k + r0 = bond.type.req logq = -beta*0.5*k_eq*(r-r0)**2 return logq From f49a85e90bf6e22d73df242c27611646fd5a30ca Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 11:28:58 -0500 Subject: [PATCH 28/57] still more units... --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 1c9701ad1..ae5360b98 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -522,8 +522,8 @@ def _angle_logq(self, theta, angle, beta): beta : float 1/kT or inverse temperature """ - k_eq = angle.type.k*units.kilocalories_per_mole/(units.degrees**2) - theta0 = angle.type.theteq*units.degrees + k_eq = angle.type.k + theta0 = angle.type.theteq logq = -beta*k_eq*0.5*(theta-theta0)**2 return logq From 625aedb088c8237e923e56b65ced1ace0c4a77b5 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 12:19:31 -0500 Subject: [PATCH 29/57] don't add another set of units... --- perses/rjmc/geometry.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index ae5360b98..3ab8f2ae6 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -353,7 +353,8 @@ def _get_valid_angles(self, atoms_with_positions, new_atom): else: if atoms_with_positions.issuperset(set([angle.atom1, angle.atom2])): eligible_angles.append(angle) - return eligible_angles + eligible_angles_with_units = [self._add_angle_units(angle) if type(angle.type.theteq) != units.Quantity else angle for angle in eligible_angles] + return eligible_angles_with_units def _autograd_ctoi(self, atom_position, bond_position, angle_position, torsion_position, calculate_jacobian=True): import autograd From 32bc181a29a52d84159dc1f2d6d1f5d01d2017f6 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 12:27:14 -0500 Subject: [PATCH 30/57] corrected sign --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 3ab8f2ae6..896ee72a7 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -660,7 +660,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor for i, xyz in enumerate(xyzs): logq_i = self._torsion_and_angle_logq(xyz, atom, positions, involved_angles, involved_torsions, beta) if np.isnan(logq_i): - logq_i = np.inf + logq_i = - np.inf logq[i]+=logq_i logq - max(logq) From 74b8a6e26b6061f4fd186d51d1b51d71d474d87e Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 13:18:52 -0500 Subject: [PATCH 31/57] changed sum --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 896ee72a7..8bfc04b3c 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -668,7 +668,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor q = np.exp(logq) #estimate the normalizing constant - Z = np.trapz(q, phis) + Z = np.sum(q) #get the normalized probabilities for torsions p = q / Z From 5cac44bfa84ccd832d931323abf96b6166b74e11 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 13:34:56 -0500 Subject: [PATCH 32/57] fix unit mistake --- perses/rjmc/geometry.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 8bfc04b3c..889333041 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -251,7 +251,8 @@ def _get_relevant_angle(self, atom1, atom2, atom3): atom3_angles = set(atom3.angles) relevant_angle_set = atom1_angles.intersection(atom2_angles, atom3_angles) relevant_angle = relevant_angle_set.pop() - relevant_angle_with_units = self._add_angle_units(relevant_angle) + if type(relevant_angle.type.k) != units.Quantity: + relevant_angle_with_units = self._add_angle_units(relevant_angle) return relevant_angle_with_units def _add_bond_units(self, bond): From e40e2a12671d58d2d1d7f87543fa07080c5a664a Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 13:37:04 -0500 Subject: [PATCH 33/57] another unit mistake. --- perses/rjmc/geometry.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 889333041..828fb3f8e 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -253,6 +253,8 @@ def _get_relevant_angle(self, atom1, atom2, atom3): relevant_angle = relevant_angle_set.pop() if type(relevant_angle.type.k) != units.Quantity: relevant_angle_with_units = self._add_angle_units(relevant_angle) + else: + relevant_angle_with_units = relevant_angle return relevant_angle_with_units def _add_bond_units(self, bond): From 8d3e5c3afc2303cdf42003f6cd95a6a6fdd5f80a Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 13:38:03 -0500 Subject: [PATCH 34/57] bond units --- perses/rjmc/geometry.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 828fb3f8e..c50427a1b 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -239,7 +239,10 @@ def _get_relevant_bond(self, atom1, atom2): bonds_2 = set(atom2.bonds) relevant_bond_set = bonds_1.intersection(bonds_2) relevant_bond = relevant_bond_set.pop() - relevant_bond_with_units = self._add_bond_units(relevant_bond) + if type(relevant_bond.type.k) != units.Quantity: + relevant_bond_with_units = self._add_bond_units(relevant_bond) + else: + relevant_bond_with_units = relevant_bond return relevant_bond_with_units def _get_relevant_angle(self, atom1, atom2, atom3): From 21a3182c663f649785b3d9b2fd5d7f50f1f9e5b3 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 14:40:06 -0500 Subject: [PATCH 35/57] fixed torsion units --- perses/rjmc/geometry.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index c50427a1b..d8ec45abc 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -581,10 +581,10 @@ def _torsion_logq(self, torsion, phi, beta): Calculate the log-unnormalized probability of the torsion """ - phi = phi/phi.unit - beta = beta/beta.unit - gamma = torsion.type.phase - V = torsion.type.phi_k + phi = phi + beta = beta + gamma = torsion.type.phase*units.degree + V = torsion.type.phi_k*units.kilocalorie_per_mole n = torsion.type.per logq = -beta*(V/2.0)*(1+units.cos(n*phi-gamma)) return logq From bc78b1cae99f853f32d6d6676e2b19c8a1645ef2 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 15:57:23 -0500 Subject: [PATCH 36/57] stop casting beta --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index d8ec45abc..915a8b902 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -89,7 +89,7 @@ def propose(self, top_proposal): logp_proposal : float The log probability of the forward-only proposal """ - beta = top_proposal.beta.in_units_of(units.kilocalorie_per_mole**-1) + beta = top_proposal.beta logp_proposal = 0.0 structure = parmed.openmm.load_topology(top_proposal.new_topology, top_proposal.new_system) new_atoms = [structure.atoms[idx] for idx in top_proposal.unique_new_atoms] @@ -163,7 +163,7 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): logp : float The log probability of the proposal for the given transformation """ - beta = top_proposal.beta.in_units_of(units.kilocalorie_per_mole**-1) + beta = top_proposal.beta logp = 0.0 top_proposal = topology_proposal.SmallMoleculeTopologyProposal() structure = parmed.openmm.load_topology(top_proposal.old_topology, top_proposal.old_system) From 8e2a4261dd04138116c11bcdcb6564e6cbb4e099 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 16:38:10 -0500 Subject: [PATCH 37/57] fix some unit stuff --- perses/rjmc/geometry.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 915a8b902..4c0c90c45 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -563,7 +563,7 @@ def _propose_bond(self, bond, beta): r0 = bond.type.req k = bond.type.k sigma_r = units.sqrt(1.0/(beta*k)) - r = sigma_r/sigma_r.unit*np.random.random()*units.angstrom + r0 + r = sigma_r*np.random.random() + r0 return r def _propose_angle(self, angle, beta): @@ -573,7 +573,7 @@ def _propose_angle(self, angle, beta): theta0 = angle.type.theteq k = angle.type.k sigma_theta = units.sqrt(1.0/(beta*k)) - theta = sigma_theta/sigma_theta.unit*np.random.random()*units.radian + theta0 + theta = sigma_theta*np.random.random() + theta0 return theta def _torsion_logq(self, torsion, phi, beta): @@ -668,7 +668,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor if np.isnan(logq_i): logq_i = - np.inf logq[i]+=logq_i - logq - max(logq) + logq -= max(logq) #exponentiate to get the unnormalized probability q = np.exp(logq) From 8f0db372579729f80191226af9449e94f9f7664a Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 16:41:45 -0500 Subject: [PATCH 38/57] fix force constants --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 4c0c90c45..4ecd25d41 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -274,7 +274,7 @@ def _add_bond_units(self, bond): """ bond.type.req = units.Quantity(bond.type.req, unit=units.angstrom) - bond.type.k = units.Quantity(bond.type.k, unit=units.kilocalorie_per_mole/units.angstrom**2) + bond.type.k = units.Quantity(2.0*bond.type.k, unit=units.kilocalorie_per_mole/units.angstrom**2) return bond def _add_angle_units(self, angle): @@ -292,7 +292,7 @@ def _add_angle_units(self, angle): The angle, but with units on its parameters """ angle.type.theteq = units.Quantity(angle.type.theteq, unit=units.degree) - angle.type.k = units.Quantity(angle.type.k, unit=units.kilocalorie_per_mole/units.radian**2) + angle.type.k = units.Quantity(2.0*angle.type.k, unit=units.kilocalorie_per_mole/units.radian**2) return angle def _get_torsions(self, atoms_with_positions, new_atom): From 30cce169f324500584fb4e885516995e8f0b525a Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 16:43:03 -0500 Subject: [PATCH 39/57] replace -inf with exception --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 4ecd25d41..0334c603b 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -666,7 +666,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor for i, xyz in enumerate(xyzs): logq_i = self._torsion_and_angle_logq(xyz, atom, positions, involved_angles, involved_torsions, beta) if np.isnan(logq_i): - logq_i = - np.inf + raise Exception("logq_i was NaN") logq[i]+=logq_i logq -= max(logq) From 4f691a2e645b381592bf66e26800a76d89ba7f68 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 16:53:08 -0500 Subject: [PATCH 40/57] change logp calculation --- perses/rjmc/geometry.py | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 0334c603b..0eafb356d 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -677,9 +677,9 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor Z = np.sum(q) #get the normalized probabilities for torsions - p = q / Z + logp = logq - np.log(Z) - return p, Z, q, phis, logq + return logp, Z, q, phis def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_position): """ @@ -706,20 +706,19 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet angle_atom = torsion.atom2 torsion_atom = torsion.atom1 internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) - p, Z, q, phis, logq = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + logp, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #find the phi that's closest to the internal_coordinate phi: phi_idx, phi = min(enumerate(phis), key=lambda x: abs(x[1]-internal_coordinates[2])) - p_phi = p[phi_idx] - return np.log(p_phi) + return logp[phi_idx] def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, positions, beta): """ Propose a torsion angle, including energetic contributions from other torsions and angles """ #first, let's get the normalizing constant of this distribution - p, Z, q, phis, logq = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + logp, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #choose from the set of possible torsion angles - phi_idx = np.random.choice(range(len(phis)), p=p) - logp = np.log(p[phi_idx]) + phi_idx = np.random.choice(range(len(phis)), p=np.exp(p)) + logp = logp[phi_idx] phi = phis[phi_idx] return phi, logp From e6bdbda2a39d638245f5a9835e9116b63b201bf9 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 17:12:30 -0500 Subject: [PATCH 41/57] fixed torsion units --- perses/rjmc/geometry.py | 27 ++++++++++++++++++++++----- 1 file changed, 22 insertions(+), 5 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 0eafb356d..32a53583f 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -295,6 +295,24 @@ def _add_angle_units(self, angle): angle.type.k = units.Quantity(2.0*angle.type.k, unit=units.kilocalorie_per_mole/units.radian**2) return angle + def _add_torsion_units(self, torsion): + """ + Add the correct units to a torsion + + Arguments + --------- + torsion : parmed.dihedral object + The torsion needing units + + Returns + ------- + torsion : parmed.dihedral object + Torsion but with units added + """ + torsion.type.phi_k = units.Quantity(torsion.type.phi_k, unit=units.kilocalorie_per_mole) + torsion.type.phase = units.Quantity(torsion.type.phase, unit=units.degree) + return torsion + def _get_torsions(self, atoms_with_positions, new_atom): """ Get the torsions that the new atom_index participates in, where all other @@ -325,7 +343,8 @@ def _get_torsions(self, atoms_with_positions, new_atom): continue if set(atoms_with_positions).issuperset(set(torsion_partners)): eligible_torsions.append(torsion) - return torsions + eligible_torsions_with_units = [self._add_torsion_units(torsion) if type(torsion.type.phi_k)!=units.Quantity else torsion for torsion in eligible_torsions] + return eligible_torsions_with_units def _get_valid_angles(self, atoms_with_positions, new_atom): """ @@ -581,10 +600,8 @@ def _torsion_logq(self, torsion, phi, beta): Calculate the log-unnormalized probability of the torsion """ - phi = phi - beta = beta - gamma = torsion.type.phase*units.degree - V = torsion.type.phi_k*units.kilocalorie_per_mole + gamma = torsion.type.phase + V = torsion.type.phi_k n = torsion.type.per logq = -beta*(V/2.0)*(1+units.cos(n*phi-gamma)) return logq From 2998ccbd17280c6432bbe06b7f52778496cf9830 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 17:19:14 -0500 Subject: [PATCH 42/57] fix torsion selection --- perses/rjmc/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 32a53583f..a9e07c3c0 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -336,9 +336,9 @@ def _get_torsions(self, atoms_with_positions, new_atom): if torsion.improper: continue if torsion.atom1 == new_atom: - torsion_partners = [torsion.atom2.idx, torsion.atom3.idx, torsion.atom4.idx] + torsion_partners = [torsion.atom2, torsion.atom3, torsion.atom4] elif torsion.atom4 == new_atom: - torsion_partners = [torsion.atom1.idx, torsion.atom2.idx, torsion.atom3.idx] + torsion_partners = [torsion.atom1, torsion.atom2, torsion.atom3] else: continue if set(atoms_with_positions).issuperset(set(torsion_partners)): From bc433e55ab4532c3389b506fae77c564c5590f1a Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 17:27:36 -0500 Subject: [PATCH 43/57] fix typo --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index a9e07c3c0..be86844b2 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -735,7 +735,7 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, #first, let's get the normalizing constant of this distribution logp, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #choose from the set of possible torsion angles - phi_idx = np.random.choice(range(len(phis)), p=np.exp(p)) + phi_idx = np.random.choice(range(len(phis)), p=np.exp(logp)) logp = logp[phi_idx] phi = phis[phi_idx] return phi, logp From c620e14ddede34ea51216330a45c0a12ab896aeb Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 17:32:05 -0500 Subject: [PATCH 44/57] make unit-ing idempotent --- perses/rjmc/geometry.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index be86844b2..2d4f06b52 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -239,10 +239,7 @@ def _get_relevant_bond(self, atom1, atom2): bonds_2 = set(atom2.bonds) relevant_bond_set = bonds_1.intersection(bonds_2) relevant_bond = relevant_bond_set.pop() - if type(relevant_bond.type.k) != units.Quantity: - relevant_bond_with_units = self._add_bond_units(relevant_bond) - else: - relevant_bond_with_units = relevant_bond + relevant_bond_with_units = self._add_bond_units(relevant_bond) return relevant_bond_with_units def _get_relevant_angle(self, atom1, atom2, atom3): @@ -273,6 +270,8 @@ def _add_bond_units(self, bond): ------- """ + if type(bond.type.k)==units.Quantity: + return bond bond.type.req = units.Quantity(bond.type.req, unit=units.angstrom) bond.type.k = units.Quantity(2.0*bond.type.k, unit=units.kilocalorie_per_mole/units.angstrom**2) return bond @@ -291,6 +290,8 @@ def _add_angle_units(self, angle): angle_with_units : parmed angle The angle, but with units on its parameters """ + if type(angle.type.k)==units.Quantity: + return angle angle.type.theteq = units.Quantity(angle.type.theteq, unit=units.degree) angle.type.k = units.Quantity(2.0*angle.type.k, unit=units.kilocalorie_per_mole/units.radian**2) return angle @@ -309,6 +310,8 @@ def _add_torsion_units(self, torsion): torsion : parmed.dihedral object Torsion but with units added """ + if type(torsion.type.phi_k)==units.Quantity: + return torsion torsion.type.phi_k = units.Quantity(torsion.type.phi_k, unit=units.kilocalorie_per_mole) torsion.type.phase = units.Quantity(torsion.type.phase, unit=units.degree) return torsion @@ -343,7 +346,7 @@ def _get_torsions(self, atoms_with_positions, new_atom): continue if set(atoms_with_positions).issuperset(set(torsion_partners)): eligible_torsions.append(torsion) - eligible_torsions_with_units = [self._add_torsion_units(torsion) if type(torsion.type.phi_k)!=units.Quantity else torsion for torsion in eligible_torsions] + eligible_torsions_with_units = [self._add_torsion_units(torsion) for torsion in eligible_torsions] return eligible_torsions_with_units def _get_valid_angles(self, atoms_with_positions, new_atom): @@ -378,7 +381,7 @@ def _get_valid_angles(self, atoms_with_positions, new_atom): else: if atoms_with_positions.issuperset(set([angle.atom1, angle.atom2])): eligible_angles.append(angle) - eligible_angles_with_units = [self._add_angle_units(angle) if type(angle.type.theteq) != units.Quantity else angle for angle in eligible_angles] + eligible_angles_with_units = [self._add_angle_units(angle) for angle in eligible_angles] return eligible_angles_with_units def _autograd_ctoi(self, atom_position, bond_position, angle_position, torsion_position, calculate_jacobian=True): From 2635d12efa35c1fca9edf58a2438fdcffd0171d8 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 21:33:21 -0500 Subject: [PATCH 45/57] replace linspace with arange --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 2d4f06b52..d4891b36b 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -674,7 +674,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor involved_torsions = self._get_torsions(atoms_with_positions, atom) #get an array of [0,2pi) - phis = units.Quantity(np.linspace(0, 2.0*np.pi, n_divisions), unit=units.radians) + phis = units.Quantity(np.arange(0, 2.0*np.pi, n_divisions), unit=units.radians) xyzs = np.zeros([n_divisions, 3]) #rotate atom about torsion angle, calculating an xyz for each From df0eb213d8d7eb32dc995daf3ae24d94d40f5f87 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Thu, 3 Dec 2015 21:55:11 -0500 Subject: [PATCH 46/57] fixed arange --- perses/rjmc/geometry.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index d4891b36b..e602e659c 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -674,8 +674,8 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor involved_torsions = self._get_torsions(atoms_with_positions, atom) #get an array of [0,2pi) - phis = units.Quantity(np.arange(0, 2.0*np.pi, n_divisions), unit=units.radians) - xyzs = np.zeros([n_divisions, 3]) + phis = units.Quantity(np.arange(0, 2.0*np.pi, (2.0*np.pi)/n_divisions), unit=units.radians) + xyzs = np.zeros([len(phis), 3]) #rotate atom about torsion angle, calculating an xyz for each for i, phi in enumerate(phis): @@ -726,7 +726,7 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet angle_atom = torsion.atom2 torsion_atom = torsion.atom1 internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) - logp, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + logp, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=500) #find the phi that's closest to the internal_coordinate phi: phi_idx, phi = min(enumerate(phis), key=lambda x: abs(x[1]-internal_coordinates[2])) return logp[phi_idx] From 1d1680699813e591c9bcd9813ec87926d4d4cfb4 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Fri, 4 Dec 2015 12:57:08 -0500 Subject: [PATCH 47/57] use more numerically stable cartesian to internal coordinate conversion --- perses/rjmc/geometry.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index e602e659c..ebbe0d176 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -417,10 +417,18 @@ def _cartesian_to_internal(xyz): plane1 = np.cross(a, b) plane2 = np.cross(b, c) - phi = np.arccos(np.dot(plane1, plane2)/(np.linalg.norm(plane1)*np.linalg.norm(plane2))) + plane1_u = plane1/np.linalg.norm(plane1) + plane2_u = plane2/np.linalg.norm(plane2) - if np.dot(np.cross(plane1, plane2), b_u) < 0: - phi = -phi + m = np.cross(plane1_u, b_u) + + x = np.dot(plane1_u, plane2_u) + y = np.dot(m, plane2_u) + + phi = np.arctan2(y, x) + + if np.isnan(phi): + raise Exception("phi is nan") return np.array([r, theta, phi]) internal_coords = _cartesian_to_internal(atom_position) @@ -726,7 +734,7 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet angle_atom = torsion.atom2 torsion_atom = torsion.atom1 internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) - logp, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=500) + logp, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #find the phi that's closest to the internal_coordinate phi: phi_idx, phi = min(enumerate(phis), key=lambda x: abs(x[1]-internal_coordinates[2])) return logp[phi_idx] From ca3f0ea48813d8b7a3f4f74569263548e82f5de1 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Sun, 6 Dec 2015 18:18:46 -0500 Subject: [PATCH 48/57] add coordinate conversion with units --- perses/rjmc/coordinate_tools.py | 108 ++++++++++++++++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 perses/rjmc/coordinate_tools.py diff --git a/perses/rjmc/coordinate_tools.py b/perses/rjmc/coordinate_tools.py new file mode 100644 index 000000000..805b36b09 --- /dev/null +++ b/perses/rjmc/coordinate_tools.py @@ -0,0 +1,108 @@ +import numpy as np +import simtk.unit as units + +def _rotation_matrix(axis, angle): + """ + This method produces a rotation matrix given an axis and an angle. + """ + axis = axis/units.norm(axis) + axis_squared = axis**2 + cos_angle = units.cos(angle) + sin_angle = units.sin(angle) + rot_matrix_row_one = np.array([cos_angle+axis_squared[0]*(1-cos_angle), + axis[0]*axis[1]*(1-cos_angle) - axis[2]*sin_angle, + axis[0]*axis[2]*(1-cos_angle)+axis[1]*sin_angle]) + + rot_matrix_row_two = np.array([axis[1]*axis[0]*(1-cos_angle)+axis[2]*sin_angle, + cos_angle+axis_squared[1]*(1-cos_angle), + axis[1]*axis[2]*(1-cos_angle) - axis[0]*sin_angle]) + + rot_matrix_row_three = np.array([axis[2]*axis[0]*(1-cos_angle)-axis[1]*sin_angle, + axis[2]*axis[1]*(1-cos_angle)+axis[0]*sin_angle, + cos_angle+axis_squared[2]*(1-cos_angle)]) + + rotation_matrix = np.array([rot_matrix_row_one, rot_matrix_row_two, rot_matrix_row_three]) + return rotation_matrix + +def _cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position): + """ + Cartesian to internal function (hah) + """ + a = bond_position - atom_position + b = angle_position - bond_position + #3-4 bond + c = angle_position - torsion_position + a_u = a / units.norm(a) + b_u = b / units.norm(b) + c_u = c / units.norm(c) + + #bond length + r = units.norm(a) + + #bond angle + theta = units.acos(units.dot(-a_u, b_u)) + + #torsion angle + plane1 = np.cross(a, b) + plane2 = np.cross(b, c) + + cos_phi = units.dot(plane1, plane2) / (units.norm(plane1)*units.norm(plane2)) + if cos_phi < -1.0: + cos_phi = -1.0 + elif cos_phi > 1.0: + cos_phi = 1.0 + + phi = units.acos(cos_phi) + + if units.dot(np.cross(plane1, plane2), b_u) > 0: + phi = -phi + + if np.isnan(phi/phi.unit): + raise Exception("phi is nan") + + return np.array([r, theta, phi]) + +def _internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi): + a = angle_position - bond_position + b = angle_position - torsion_position + + a_u = a / units.norm(a) + + + d_r = r*a_u + + normal = np.cross(a, b) + + #construct the angle rotation matrix + angle_axis = normal / units.norm(normal) + angle_rotation_matrix = _rotation_matrix(angle_axis, theta) + + #apply it + d_ang = units.dot(angle_rotation_matrix, d_r) + + #construct the torsion rotation matrix and apply it + torsion_axis = a_u + torsion_rotation_matrix = _rotation_matrix(torsion_axis, phi) + #apply it + d_torsion = units.dot(torsion_rotation_matrix, d_ang) + + #add the positions of the bond atom + xyz = bond_position + d_torsion + return xyz + +if __name__=="__main__": + example_coordinates = units.Quantity(np.random.normal(size=[100,3]), unit=units.nanometers) + #try to transform random coordinates to and from + for i in range(20): + indices = np.random.randint(100, size=4) + atom_position = example_coordinates[indices[0]] + bond_position = example_coordinates[indices[1]] + angle_position = example_coordinates[indices[2]] + torsion_position = example_coordinates[indices[3]] + rtp = _cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) + r = rtp[0] + theta = rtp[1] + phi = rtp[2] + xyz = _internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) + print(units.norm(atom_position)) + #print(units.norm(xyz-atom_position)) \ No newline at end of file From 73e137d758d9e8795695806494643a611455dbe8 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Sun, 6 Dec 2015 18:22:09 -0500 Subject: [PATCH 49/57] change parmtop input --- perses/tests/test_geometry_engine.py | 30 ++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index cb8ee846f..475c9e2f3 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -7,6 +7,7 @@ import openeye.oeomega as oeomega import simtk.openmm.app as app import simtk.unit as units +import numpy as np kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA beta = 1.0/ (300.0*units.kelvin*kB) @@ -36,7 +37,8 @@ def oemol_to_openmm_system(oemol, molecule_name): _ , tripos_mol2_filename = openmoltools.openeye.molecule_to_mol2(oemol, tripos_mol2_filename=molecule_name + '.tripos.mol2', conformer=0, residue_name='MOL') gaff_mol2, frcmod = openmoltools.amber.run_antechamber(molecule_name, tripos_mol2_filename) prmtop_file, inpcrd_file = openmoltools.amber.run_tleap(molecule_name, gaff_mol2, frcmod) - prmtop = app.AmberPrmtopFile(prmtop_file) + from parmed.amber import AmberParm + prmtop = AmberParm(prmtop_file) system = prmtop.createSystem(implicitSolvent=app.OBC1) crd = app.AmberInpcrdFile(inpcrd_file) return system, crd.getPositions(asNumpy=True), prmtop.topology @@ -82,6 +84,8 @@ def test_run_geometry_engine(): old_positions=pos1, logp_proposal=0.0, new_to_old_atom_map=new_to_old_atom_mapping, metadata={'test':0.0}) sm_top_proposal._beta = beta geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + test_pdb_file = open("npentanehexane.pdb", 'w') + integrator = openmm.VerletIntegrator(1*units.femtoseconds) context = openmm.Context(sys2, integrator) @@ -90,11 +94,33 @@ def test_run_geometry_engine(): print("Energy before proposal is: %s" % str(state.getPotentialEnergy())) new_positions, logp_proposal = geometry_engine.propose(sm_top_proposal) + app.PDBFile.writeFile(top2, new_positions, file=test_pdb_file) + test_pdb_file.close() context.setPositions(new_positions) state2 = context.getState(getEnergy=True) print("Energy after proposal is: %s" %str(state2.getPotentialEnergy())) +def test_coordinate_conversion(): + + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + example_coordinates = units.Quantity(np.random.normal(size=[100,3]), unit=units.nanometers) + #try to transform random coordinates to and from + for i in range(20): + indices = np.random.randint(100, size=4) + atom_position = example_coordinates[indices[0]] + bond_position = example_coordinates[indices[1]] + angle_position = example_coordinates[indices[2]] + torsion_position = example_coordinates[indices[3]] + rtp, detJ = geometry_engine._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) + r = rtp[0]*units.nanometers + theta = rtp[1]*units.radians + phi = rtp[2]*units.radians + xyz, _ = geometry_engine._internal_to_cartesian(indices[1], indices[2], indices[3], r, theta, phi, example_coordinates) + + if __name__=="__main__": - test_run_geometry_engine() \ No newline at end of file + test_coordinate_conversion() + test_run_geometry_engine() From 0ba012b775ec4bdea15916d9f33c135dbba1bd89 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Mon, 7 Dec 2015 16:37:01 -0500 Subject: [PATCH 50/57] fixes --- perses/rjmc/coordinate_tools.py | 70 +++---- perses/rjmc/geometry.py | 285 +++++++++++++-------------- perses/tests/test_geometry_engine.py | 204 +++++++++++++++++-- 3 files changed, 352 insertions(+), 207 deletions(-) diff --git a/perses/rjmc/coordinate_tools.py b/perses/rjmc/coordinate_tools.py index 805b36b09..909efd1f0 100644 --- a/perses/rjmc/coordinate_tools.py +++ b/perses/rjmc/coordinate_tools.py @@ -1,14 +1,13 @@ import numpy as np -import simtk.unit as units def _rotation_matrix(axis, angle): """ This method produces a rotation matrix given an axis and an angle. """ - axis = axis/units.norm(axis) + axis = axis/np.linalg.norm(axis) axis_squared = axis**2 - cos_angle = units.cos(angle) - sin_angle = units.sin(angle) + cos_angle = np.cos(angle) + sin_angle = np.sin(angle) rot_matrix_row_one = np.array([cos_angle+axis_squared[0]*(1-cos_angle), axis[0]*axis[1]*(1-cos_angle) - axis[2]*sin_angle, axis[0]*axis[2]*(1-cos_angle)+axis[1]*sin_angle]) @@ -26,39 +25,44 @@ def _rotation_matrix(axis, angle): def _cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position): """ - Cartesian to internal function (hah) + Cartesian to internal function """ - a = bond_position - atom_position + a = atom_position - bond_position b = angle_position - bond_position #3-4 bond c = angle_position - torsion_position - a_u = a / units.norm(a) - b_u = b / units.norm(b) - c_u = c / units.norm(c) + a_u = a / np.linalg.norm(a) + b_u = b / np.linalg.norm(b) + c_u = c / np.linalg.norm(c) #bond length - r = units.norm(a) + r = np.linalg.norm(a) #bond angle - theta = units.acos(units.dot(-a_u, b_u)) + cos_theta = np.dot(a_u, b_u) + if cos_theta > 1.0: + cos_theta = 1.0 + elif cos_theta < -1.0: + cos_theta = -1.0 + theta = np.arccos(cos_theta) #torsion angle - plane1 = np.cross(a, b) - plane2 = np.cross(b, c) + plane1 = np.cross(a_u, b_u) + plane2 = np.cross(b_u, c_u) - cos_phi = units.dot(plane1, plane2) / (units.norm(plane1)*units.norm(plane2)) + cos_phi = np.dot(plane1, plane2) / (np.linalg.norm(plane1)*np.linalg.norm(plane2)) if cos_phi < -1.0: cos_phi = -1.0 elif cos_phi > 1.0: cos_phi = 1.0 - phi = units.acos(cos_phi) + phi = np.arccos(cos_phi) - if units.dot(np.cross(plane1, plane2), b_u) > 0: + if np.dot(a, plane2) <= 0: phi = -phi - if np.isnan(phi/phi.unit): - raise Exception("phi is nan") + if np.isnan(phi): + raise Exception("phi is nan, cos_phi is %s" % str(cos_phi)) return np.array([r, theta, phi]) @@ -66,43 +70,27 @@ def _internal_to_cartesian(bond_position, angle_position, torsion_position, r, t a = angle_position - bond_position b = angle_position - torsion_position - a_u = a / units.norm(a) + a_u = a / np.linalg.norm(a) + b_u = b / np.linalg.norm(b) d_r = r*a_u - normal = np.cross(a, b) + normal = np.cross(a_u, b_u) #construct the angle rotation matrix - angle_axis = normal / units.norm(normal) + angle_axis = normal / np.linalg.norm(normal) angle_rotation_matrix = _rotation_matrix(angle_axis, theta) #apply it - d_ang = units.dot(angle_rotation_matrix, d_r) + d_ang = np.dot(angle_rotation_matrix, d_r) #construct the torsion rotation matrix and apply it torsion_axis = a_u - torsion_rotation_matrix = _rotation_matrix(torsion_axis, phi) + torsion_rotation_matrix = _rotation_matrix(torsion_axis, -(phi+np.pi)) #apply it - d_torsion = units.dot(torsion_rotation_matrix, d_ang) + d_torsion = np.dot(torsion_rotation_matrix, d_ang) #add the positions of the bond atom xyz = bond_position + d_torsion return xyz - -if __name__=="__main__": - example_coordinates = units.Quantity(np.random.normal(size=[100,3]), unit=units.nanometers) - #try to transform random coordinates to and from - for i in range(20): - indices = np.random.randint(100, size=4) - atom_position = example_coordinates[indices[0]] - bond_position = example_coordinates[indices[1]] - angle_position = example_coordinates[indices[2]] - torsion_position = example_coordinates[indices[3]] - rtp = _cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) - r = rtp[0] - theta = rtp[1] - phi = rtp[2] - xyz = _internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) - print(units.norm(atom_position)) - #print(units.norm(xyz-atom_position)) \ No newline at end of file diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index ebbe0d176..250735177 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -4,12 +4,12 @@ """ from collections import namedtuple import perses.rjmc.topology_proposal as topology_proposal -import numpy as np import parmed import simtk.unit as units import logging - - +import matplotlib.pyplot as pyplot +import numpy as np +import coordinate_tools GeometryProposal = namedtuple('GeometryProposal',['new_positions','logp']) @@ -124,27 +124,34 @@ def propose(self, top_proposal): logp_r = self._bond_logq(r_proposed, bond, beta) - logZ_r #propose an angle and calculate its probability - angle = self._get_relevant_angle(atom, bond_atom, angle_atom) - theta_proposed = self._propose_angle(angle, beta) - angle_k = angle.type.k - sigma_theta = units.sqrt(1/(beta*angle_k)) - logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) - logp_theta = self._angle_logq(theta_proposed, angle, beta) - logZ_theta - - #propose a torsion angle and calcualate its probability - phi_proposed, logp_phi = self._propose_torsion(atom, r_proposed, theta_proposed, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, beta) + propose_angle_separately = True + if propose_angle_separately: + angle = self._get_relevant_angle(atom, bond_atom, angle_atom) + theta_proposed = self._propose_angle(angle, beta) + angle_k = angle.type.k + sigma_theta = units.sqrt(1/(beta*angle_k)) + logZ_theta = np.log((np.sqrt(2*np.pi)*sigma_theta/sigma_theta.unit)) + logp_theta = self._angle_logq(theta_proposed, angle, beta) - logZ_theta + + #propose a torsion angle and calcualate its probability + phi_proposed, logp_phi = self._propose_torsion(atom, r_proposed, theta_proposed, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, beta) + else: + theta_proposed, phi_proposed, logp_theta_phi, xyz = self._joint_torsion_angle_proposal(atom, r_proposed, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, beta) #convert to cartesian - xyz, detJ = self._autograd_itoc(bond_atom.idx, angle_atom.idx, torsion_atom.idx, r_proposed, theta_proposed, phi_proposed, new_positions) + xyz, detJ = self._internal_to_cartesian(new_positions[bond_atom.idx], new_positions[angle_atom.idx], new_positions[torsion_atom.idx], r_proposed, theta_proposed, phi_proposed) + detJ=0.0 + rtp_test = self._cartesian_to_internal(xyz, new_positions[bond_atom.idx], new_positions[angle_atom.idx], new_positions[torsion_atom.idx]) #add new position to array of new positions new_positions[atom.idx] = xyz #accumulate logp - logp_proposal = logp_proposal + logp_choice + logp_r + logp_theta + logp_phi + np.log(detJ) + logp_proposal = logp_proposal + logp_choice + logp_r + logp_theta +logp_phi + np.log(detJ+1) atoms_with_positions.append(atom) new_atoms.remove(atom) return new_positions, logp_proposal + def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): """ Calculate the logp for the given geometry proposal @@ -195,7 +202,7 @@ def logp_reverse(self, top_proposal, new_coordinates, old_coordinates): bond_coords = reverse_proposal_coordinates[bond_atom.idx] angle_coords = reverse_proposal_coordinates[angle_atom.idx] torsion_coords = reverse_proposal_coordinates[torsion_atom.idx] - internal_coordinates, detJ = self._autograd_ctoi(atom_coords, bond_coords, angle_coords, torsion_coords) + internal_coordinates, detJ = self._cartesian_to_internal(atom_coords, bond_coords, angle_coords, torsion_coords) #propose a bond and calculate its probability bond = self._get_relevant_bond(atom, bond_atom) @@ -384,147 +391,59 @@ def _get_valid_angles(self, atoms_with_positions, new_atom): eligible_angles_with_units = [self._add_angle_units(angle) for angle in eligible_angles] return eligible_angles_with_units - def _autograd_ctoi(self, atom_position, bond_position, angle_position, torsion_position, calculate_jacobian=True): - import autograd - import autograd.numpy as np - - atom_position = atom_position/atom_position.unit - bond_position = bond_position/bond_position.unit - angle_position = angle_position/angle_position.unit - torsion_position = torsion_position/torsion_position.unit - - - def _cartesian_to_internal(xyz): - """ - Autograd-based jacobian of transformation from cartesian to internal. - Returns without units! - """ - a = bond_position - xyz - b = angle_position - bond_position - #3-4 bond - c = angle_position - torsion_position - a_u = a / np.linalg.norm(a) - b_u = b / np.linalg.norm(b) - c_u = c / np.linalg.norm(c) - - #bond length - r = np.linalg.norm(a) + def _rotation_matrix(self, axis, angle): + """ + This method produces a rotation matrix given an axis and an angle. + """ + axis = axis/np.linalg.norm(axis) + axis_squared = np.square(axis) + cos_angle = np.cos(angle) + sin_angle = np.sin(angle) + rot_matrix_row_one = np.array([cos_angle+axis_squared[0]*(1-cos_angle), + axis[0]*axis[1]*(1-cos_angle) - axis[2]*sin_angle, + axis[0]*axis[2]*(1-cos_angle)+axis[1]*sin_angle]) - #bond angle - theta = np.arccos(np.dot(-a_u, b_u)) + rot_matrix_row_two = np.array([axis[1]*axis[0]*(1-cos_angle)+axis[2]*sin_angle, + cos_angle+axis_squared[1]*(1-cos_angle), + axis[1]*axis[2]*(1-cos_angle) - axis[0]*sin_angle]) - #torsion angle - plane1 = np.cross(a, b) - plane2 = np.cross(b, c) + rot_matrix_row_three = np.array([axis[2]*axis[0]*(1-cos_angle)-axis[1]*sin_angle, + axis[2]*axis[1]*(1-cos_angle)+axis[0]*sin_angle, + cos_angle+axis_squared[2]*(1-cos_angle)]) - plane1_u = plane1/np.linalg.norm(plane1) - plane2_u = plane2/np.linalg.norm(plane2) + rotation_matrix = np.array([rot_matrix_row_one, rot_matrix_row_two, rot_matrix_row_three]) + return rotation_matrix - m = np.cross(plane1_u, b_u) + def _cartesian_to_internal(self, atom_position, bond_position, angle_position, torsion_position): + """ + Cartesian to internal function (hah) + """ + #ensure we have the correct units, then remove them + atom_position = atom_position.in_units_of(units.nanometers)/units.nanometers + bond_position = bond_position.in_units_of(units.nanometers)/units.nanometers + angle_position = angle_position.in_units_of(units.nanometers)/units.nanometers + torsion_position = torsion_position.in_units_of(units.nanometers)/units.nanometers - x = np.dot(plane1_u, plane2_u) - y = np.dot(m, plane2_u) + internal_coords = coordinate_tools._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) - phi = np.arctan2(y, x) - if np.isnan(phi): - raise Exception("phi is nan") + return internal_coords, 0 - return np.array([r, theta, phi]) - internal_coords = _cartesian_to_internal(atom_position) - if calculate_jacobian: - j = autograd.jacobian(_cartesian_to_internal) - jacobian_det = np.linalg.det(j(atom_position)) - else: - jacobian_det = 0.0 - return internal_coords, np.abs(jacobian_det) - def _autograd_itoc(self, bond, angle, torsion, r, theta, phi, positions, calculate_jacobian=True): + def _internal_to_cartesian(self, bond_position, angle_position, torsion_position, r, theta, phi): """ - Autograd based coordinate conversion internal -> cartesian - - Arguments - --------- - bond : int - index of the bonded atom - angle : int - index of the angle atom - torsion : int - index of the torsion atom - r : float, Quantity nm - bond length - theta : float, Quantity rad - bond angle - phi : float, Quantity rad - torsion angle - positions : [n, 3] np.array Quantity nm - positions of the atoms in the molecule - calculate_jacobian : boolean - Whether to calculate a jacobian--default True + Calculate the cartesian coordinates given the internal, as well as abs(detJ) + """ + r = r.in_units_of(units.nanometers)/units.nanometers + theta = theta.in_units_of(units.radians)/units.radians + phi = phi.in_units_of(units.radians)/units.radians + bond_position = bond_position.in_units_of(units.nanometers)/units.nanometers + angle_position = angle_position.in_units_of(units.nanometers)/units.nanometers + torsion_position = torsion_position.in_units_of(units.nanometers)/units.nanometers + xyz = coordinate_tools._internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) + xyz = units.Quantity(xyz, unit=units.nanometers) + return xyz, 0 - Returns - ------- - atom_xyz : [1, 3] np.array Quantity nm - The atomic positions in cartesian space - detJ : float - The absolute value of the determinant of the jacobian - """ - import autograd - import autograd.numpy as np - positions = positions/positions.unit - r = r/r.unit - theta = theta/theta.unit - phi = phi/phi.unit - rtp = np.array([r, theta, phi]) - - def _internal_to_cartesian(rthetaphi): - - a = positions[angle] - positions[bond] - b = positions[angle] - positions[torsion] - - a_u = a / np.linalg.norm(a) - b_u = b / np.linalg.norm(b) - - - d_r = rthetaphi[0]*a_u - - normal = np.cross(a, b) - - #construct the angle rotation matrix - axis_angle = normal / np.linalg.norm(normal) - a = np.cos(rthetaphi[1]/2) - b, c, d = -axis_angle*np.sin(rthetaphi[1]/2) - angle_rotation_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], - [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], - [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) - #apply it - d_ang = np.dot(angle_rotation_matrix, d_r) - - #construct the torsion rotation matrix and apply it - axis = a_u - a = np.cos(rthetaphi[2]/2) - b, c, d = -axis*np.sin(rthetaphi[2]/2) - torsion_rotation_matrix = np.array([[a*a+b*b-c*c-d*d, 2*(b*c-a*d), 2*(b*d+a*c)], - [2*(b*c+a*d), a*a+c*c-b*b-d*d, 2*(c*d-a*b)], - [2*(b*d-a*c), 2*(c*d+a*b), a*a+d*d-b*b-c*c]]) - #apply it - d_torsion = np.dot(torsion_rotation_matrix, d_ang) - - #add the positions of the bond atom - xyz = positions[bond] + d_torsion - - return xyz - - atom_xyz = _internal_to_cartesian(rtp) - if calculate_jacobian: - j = autograd.jacobian(_internal_to_cartesian) - jacobian_det = np.linalg.det(j(rtp)) - logging.debug("detJ is %f" %(jacobian_det)) - else: - jacobian_det = 0.0 - #detj_spherical = r**2*np.sin(theta) - #logging.debug("The spherical detJ is %f" % detj_spherical) - return units.Quantity(atom_xyz, unit=units.nanometers), np.abs(jacobian_det) def _bond_logq(self, r, bond, beta): """ @@ -623,6 +542,9 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, positions, beta): pass + def _joint_torsion_angle_proposal(self, atom, r, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, new_positions, beta): + return 0, 0, 0 + class FFAllAngleGeometryEngine(FFGeometryEngine): """ This is a forcefield-based geometry engine that takes all relevant angles @@ -644,15 +566,22 @@ def _torsion_and_angle_logq(self, xyz, atom, positions, involved_angles, involve bond_atom_position = xyz if angle.atom2 == atom else positions[angle.atom2.idx] angle_atom_position = xyz if angle.atom3 == atom else positions[angle.atom3.idx] theta = self._calculate_angle(atom_position, bond_atom_position, angle_atom_position) - logq_angles += self._angle_logq(theta*units.radians, angle, beta) + logq_i = self._angle_logq(theta*units.radians, angle, beta) + if np.isnan(logq_i): + raise Exception("Angle logq is nan") + logq_angles+=logq_i for torsion in involved_torsions: atom_position = xyz if torsion.atom1 == atom else positions[torsion.atom1.idx] bond_atom_position = xyz if torsion.atom2 == atom else positions[torsion.atom2.idx] angle_atom_position = xyz if torsion.atom3 == atom else positions[torsion.atom3.idx] torsion_atom_position = xyz if torsion.atom4 == atom else positions[torsion.atom4.idx] - internal_coordinates, _ = self._autograd_ctoi(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position, calculate_jacobian=False) + internal_coordinates, _ = self._cartesian_to_internal(atom_position, bond_atom_position, angle_atom_position, torsion_atom_position) phi = internal_coordinates[2]*units.radians - logq_torsions += self._torsion_logq(torsion, phi, beta) + logq_i = self._torsion_logq(torsion, phi, beta) + if np.isnan(logq_i): + raise Exception("Torsion logq is nan") + logq_torsions+=logq_i + return logq_angles+logq_torsions @@ -687,7 +616,7 @@ def _normalize_torsion_proposal(self, atom, r, theta, bond_atom, angle_atom, tor #rotate atom about torsion angle, calculating an xyz for each for i, phi in enumerate(phis): - xyzs[i], _ = self._autograd_itoc(bond_atom.idx, angle_atom.idx, torsion_atom.idx, r, theta, phi, positions, calculate_jacobian=False) + xyzs[i], _ = self._internal_to_cartesian(positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx], r, theta, phi) #set up arrays for energies from angles and torsions logq = np.zeros(n_divisions) @@ -717,7 +646,12 @@ def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_positio b = angle_atom_position - bond_atom_position a_u = a / np.linalg.norm(a) b_u = b / np.linalg.norm(b) - theta = np.arccos(np.dot(a_u, b_u)) + cos_theta = np.dot(a_u, b_u) + if cos_theta > 1.0: + cos_theta = 1.0 + elif cos_theta < -1.0: + cos_theta = -1.0 + theta = np.arccos(cos_theta) return theta def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, beta): @@ -733,7 +667,7 @@ def _torsion_logp(self, atom, xyz, torsion, atoms_with_positions, positions, bet bond_atom = torsion.atom3 angle_atom = torsion.atom2 torsion_atom = torsion.atom1 - internal_coordinates = self._autograd_ctoi(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) + internal_coordinates = self._cartesian_to_internal(xyz, positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx]) logp, Z, q, phis = self._normalize_torsion_proposal(atom, internal_coordinates[0], internal_coordinates[1], bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) #find the phi that's closest to the internal_coordinate phi: phi_idx, phi = min(enumerate(phis), key=lambda x: abs(x[1]-internal_coordinates[2])) @@ -744,9 +678,56 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, Propose a torsion angle, including energetic contributions from other torsions and angles """ #first, let's get the normalizing constant of this distribution - logp, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000) + logp, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=100) #choose from the set of possible torsion angles phi_idx = np.random.choice(range(len(phis)), p=np.exp(logp)) + pyplot.plot(phis, np.exp(logp)) logp = logp[phi_idx] phi = phis[phi_idx] return phi, logp + + def _normalize_torsion_and_angle(self, atom, r, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=5000): + """ + Construct a grid of points on the surface of a sphere bounded by the bond length + """ + involved_angles = self._get_valid_angles(atoms_with_positions, atom) + involved_torsions = self._get_torsions(atoms_with_positions, atom) + + xyz_grid = units.Quantity(np.zeros([n_divisions, n_divisions, 3]), unit=units.nanometers) + thetas = units.Quantity(np.arange(0, np.pi, np.pi/n_divisions), unit=units.radians) + phis = units.Quantity(np.arange(0, 2.0*np.pi, (2.0*np.pi)/n_divisions), unit=units.radians) + + #now, for every theta, phi combination, get the xyz coordinates: + for i in range(len(thetas)): + for j in range(len(phis)): + xyz_grid[i, j, :], _ = self._internal_to_cartesian(positions[bond_atom.idx], positions[angle_atom.idx], positions[torsion_atom.idx], r, thetas[i], phis[j]) + + logq = np.zeros([n_divisions, n_divisions]) + + #now, for each xyz point in the grid, calculate the log q(theta, phi) + for i in range(len(thetas)): + for j in range(len(phis)): + logq[i, j] = self._torsion_and_angle_logq(xyz_grid[i, j], atom, positions, involved_angles, involved_torsions, beta) + logq -= np.max(logq) + q = np.exp(logq) + Z = np.sum(q) + logp = logq - np.log(Z) + return logp, Z, q, thetas, phis, xyz_grid + + def _joint_torsion_angle_proposal(self, atom, r, bond_atom, angle_atom, torsion_atom, torsion, atoms_with_positions, positions, beta): + n_divisions = 50 + logp, Z, q, thetas, phis, xyz_grid = self._normalize_torsion_and_angle(atom, r, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=n_divisions) + logp_flat = np.ravel(logp) + #xyz_flat = np.ravel(xyz_grid) + #import matplotlib.pyplot as plt + #from mpl_toolkits.mplot3d import Axes3D + #fig = plt.figure() + #ax = fig.add_subplot(111, projection='3d') + #ax.scatter(xyz_grid[:,:,0], xyz_grid[:,:,1], xyz_grid[:,:,2]) + theta_phi_idx = np.random.choice(range(len(logp_flat)), p=np.exp(logp_flat)) + theta_idx, phi_idx = np.unravel_index(theta_phi_idx, [n_divisions, n_divisions]) + xyz = xyz_grid[theta_idx, phi_idx] + logp = logp_flat[theta_phi_idx] + theta = thetas[theta_idx] + phi = phis[phi_idx] + return theta, phi, logp, xyz \ No newline at end of file diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 475c9e2f3..563f258fb 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -8,6 +8,7 @@ import simtk.openmm.app as app import simtk.unit as units import numpy as np +import parmed kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA beta = 1.0/ (300.0*units.kelvin*kB) @@ -67,8 +68,8 @@ def test_run_geometry_engine(): Run the geometry engine a few times to make sure that it actually runs without exceptions. Convert n-pentane to 2-methylpentane """ - molecule_name_1 = 'n-pentane' - molecule_name_2 = '2-methylpentane' + molecule_name_1 = 'propane' + molecule_name_2 = 'butane' molecule1 = generate_initial_molecule(molecule_name_1) molecule2 = generate_initial_molecule(molecule_name_2) @@ -84,7 +85,7 @@ def test_run_geometry_engine(): old_positions=pos1, logp_proposal=0.0, new_to_old_atom_map=new_to_old_atom_mapping, metadata={'test':0.0}) sm_top_proposal._beta = beta geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) - test_pdb_file = open("npentanehexane.pdb", 'w') + test_pdb_file = open("pbutane_4.pdb", 'w') integrator = openmm.VerletIntegrator(1*units.femtoseconds) @@ -100,27 +101,202 @@ def test_run_geometry_engine(): state2 = context.getState(getEnergy=True) print("Energy after proposal is: %s" %str(state2.getPotentialEnergy())) -def test_coordinate_conversion(): +def test_existing_coordinates(): + """ + predUS + for each torsion, calculate position of atom1 + """ + molecule_name_2 = 'butane' + molecule2 = generate_initial_molecule(molecule_name_2) + sys, pos, top = oemol_to_openmm_system(molecule2, molecule_name_2) + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + structure = parmed.openmm.load_topology(top, sys) + torsions = [torsion for torsion in structure.dihedrals if not torsion.improper] + for torsion in torsions: + atom1_position = pos[torsion.atom1.idx] + atom2_position = pos[torsion.atom2.idx] + atom3_position = pos[torsion.atom3.idx] + atom4_position = pos[torsion.atom4.idx] + _internal_coordinates, _ = geometry_engine._cartesian_to_internal(atom1_position, atom2_position, atom3_position, atom4_position) + internal_coordinates = internal_in_units(_internal_coordinates) + recalculated_atom1_position, _ = geometry_engine._internal_to_cartesian(atom2_position, atom3_position, atom4_position, internal_coordinates[0], internal_coordinates[1], internal_coordinates[2]) + n = np.linalg.norm(atom1_position-recalculated_atom1_position) + print(n) + +def internal_in_units(internal_coords): + r = internal_coords[0]*units.nanometers if type(internal_coords[0]) != units.Quantity else internal_coords[0] + theta = internal_coords[1]*units.radians if type(internal_coords[1]) != units.Quantity else internal_coords[1] + phi = internal_coords[2]*units.radians if type(internal_coords[2]) != units.Quantity else internal_coords[2] + return [r, theta, phi] +def test_coordinate_conversion(): import perses.rjmc.geometry as geometry geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) example_coordinates = units.Quantity(np.random.normal(size=[100,3]), unit=units.nanometers) #try to transform random coordinates to and from - for i in range(20): + for i in range(200): indices = np.random.randint(100, size=4) - atom_position = example_coordinates[indices[0]] - bond_position = example_coordinates[indices[1]] - angle_position = example_coordinates[indices[2]] - torsion_position = example_coordinates[indices[3]] + atom_position = units.Quantity(np.array([ 0.80557722 ,-1.10424644 ,-1.08578826]), unit=units.nanometers) + bond_position = units.Quantity(np.array([ 0.0765, 0.1 , -0.4005]), unit=units.nanometers) + angle_position = units.Quantity(np.array([ 0.0829 , 0.0952 ,-0.2479]) ,unit=units.nanometers) + torsion_position = units.Quantity(np.array([-0.057 , 0.0951 ,-0.1863] ) ,unit=units.nanometers) rtp, detJ = geometry_engine._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) - r = rtp[0]*units.nanometers - theta = rtp[1]*units.radians - phi = rtp[2]*units.radians - xyz, _ = geometry_engine._internal_to_cartesian(indices[1], indices[2], indices[3], r, theta, phi, example_coordinates) + r = rtp[0] + theta = rtp[1] + phi = rtp[2] + xyz, _ = geometry_engine._internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) + assert np.linalg.norm(xyz-atom_position) < 1.0e-12 + +def test_dihedral_potential(): + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + molecule_name = 'ethane' + molecule2 = generate_initial_molecule(molecule_name) + sys, pos, top = oemol_to_openmm_system(molecule2, molecule_name) + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + structure = parmed.openmm.load_topology(top, sys) + +def test_openmm_dihedral(): + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + import simtk.openmm as openmm + integrator = openmm.VerletIntegrator(1.0*units.femtoseconds) + sys = openmm.System() + force = openmm.CustomTorsionForce("theta") + for i in range(4): + sys.addParticle(1.0*units.amu) + force.addTorsion(0,1,2,3) + sys.addForce(force) + atom_position = units.Quantity(np.array([ 0.10557722 ,-1.10424644 ,-1.08578826]), unit=units.nanometers) + bond_position = units.Quantity(np.array([ 0.0765, 0.1 , -0.4005]), unit=units.nanometers) + angle_position = units.Quantity(np.array([ 0.0829 , 0.0952 ,-0.2479]) ,unit=units.nanometers) + torsion_position = units.Quantity(np.array([-0.057 , 0.0951 ,-0.1863] ) ,unit=units.nanometers) + rtp, detJ = geometry_engine._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) + platform = openmm.Platform.getPlatformByName("Reference") + context = openmm.Context(sys, integrator, platform) + positions = [atom_position, bond_position, angle_position, torsion_position] + context.setPositions(positions) + state = context.getState(getEnergy=True) + potential = state.getPotentialEnergy() + + #rotate about the torsion: + n_divisions = 100 + phis = units.Quantity(np.arange(0, 2.0*np.pi, (2.0*np.pi)/n_divisions), unit=units.radians) + omm_phis = np.zeros(n_divisions) + for i, phi in enumerate(phis): + xyz_atom1, _ = geometry_engine._internal_to_cartesian(bond_position, angle_position, torsion_position, rtp[0]*units.nanometers, rtp[1]*units.radians, phi) + context.setPositions([xyz_atom1, bond_position, angle_position, torsion_position]) + state = context.getState(getEnergy=True) + omm_phis[i] = state.getPotentialEnergy()/units.kilojoule_per_mole + + return 0 + +def _try_random_itoc(): + + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + import simtk.openmm as openmm + integrator = openmm.VerletIntegrator(1.0*units.femtoseconds) + sys = openmm.System() + force = openmm.CustomTorsionForce("theta") + for i in range(4): + sys.addParticle(1.0*units.amu) + force.addTorsion(0,1,2,3) + sys.addForce(force) + atom_position = units.Quantity(np.array([ 0.10557722 ,-1.10424644 ,-1.08578826]), unit=units.nanometers) + bond_position = units.Quantity(np.array([ 0.0765, 0.1 , -0.4005]), unit=units.nanometers) + angle_position = units.Quantity(np.array([ 0.0829 , 0.0952 ,-0.2479]) ,unit=units.nanometers) + torsion_position = units.Quantity(np.array([-0.057 , 0.0951 ,-0.1863] ) ,unit=units.nanometers) + for i in range(10): + atom_position += units.Quantity(np.random.normal(size=3), unit=units.nanometers) + r, theta, phi = _get_internal_from_omm(atom_position, bond_position, angle_position, torsion_position) + r = (r/r.unit)*units.nanometers + theta = (theta/theta.unit)*units.radians + phi = (phi/phi.unit)*units.radians + recomputed_xyz, _ = geometry_engine._internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) + new_r, new_theta, new_phi = _get_internal_from_omm(recomputed_xyz,bond_position, angle_position, torsion_position) + crtp = geometry_engine._cartesian_to_internal(recomputed_xyz,bond_position, angle_position, torsion_position) + print(atom_position-recomputed_xyz) + + +def _get_internal_from_omm(atom_coords, bond_coords, angle_coords, torsion_coords): + import copy + #master system, will be used for all three + sys = openmm.System() + platform = openmm.Platform.getPlatformByName("Reference") + for i in range(4): + sys.addParticle(1.0*units.amu) + + #first, the bond length: + bond_sys = openmm.System() + bond_sys.addParticle(1.0*units.amu) + bond_sys.addParticle(1.0*units.amu) + bond_force = openmm.CustomBondForce("r") + bond_force.addBond(0, 1) + bond_sys.addForce(bond_force) + bond_integrator = openmm.VerletIntegrator(1*units.femtoseconds) + bond_context = openmm.Context(bond_sys, bond_integrator, platform) + bond_context.setPositions([atom_coords, bond_coords]) + bond_state = bond_context.getState(getEnergy=True) + r = bond_state.getPotentialEnergy() + del bond_sys, bond_context, bond_integrator + + #now, the angle: + angle_sys = copy.deepcopy(sys) + angle_force = openmm.CustomAngleForce("theta") + angle_force.addAngle(0,1,2) + angle_sys.addForce(angle_force) + angle_integrator = openmm.VerletIntegrator(1*units.femtoseconds) + angle_context = openmm.Context(angle_sys, angle_integrator, platform) + angle_context.setPositions([atom_coords, bond_coords, angle_coords, torsion_coords]) + angle_state = angle_context.getState(getEnergy=True) + theta = angle_state.getPotentialEnergy() + del angle_sys, angle_context, angle_integrator + + #finally, the torsion: + torsion_sys = copy.deepcopy(sys) + torsion_force = openmm.CustomTorsionForce("theta") + torsion_force.addTorsion(0,1,2,3) + torsion_sys.addForce(torsion_force) + torsion_integrator = openmm.VerletIntegrator(1*units.femtoseconds) + torsion_context = openmm.Context(torsion_sys, torsion_integrator, platform) + torsion_context.setPositions([atom_coords, bond_coords, angle_coords, torsion_coords]) + torsion_state = torsion_context.getState(getEnergy=True) + phi = torsion_state.getPotentialEnergy() + del torsion_sys, torsion_context, torsion_integrator + + return r, theta, phi + + +def test_rotation_matrix(): + """ + Test the rotation matrix code to ensure that it does what we expect + """ + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + #get 3 points: + points = units.Quantity(np.random.normal(size=[3,3]), unit=units.nanometers) + #calculate the angle between them: + a = points[1] - points[0] + b = points[2] - points[1] + a_u = a / units.norm(a) + b_u = b / units.norm(b) + theta_initial = units.acos(units.dot(a_u, b_u)) + #now, rotate point 0 about an axis: + theta_rotate = np.pi/2.0*units.radians + normal = np.cross(a, b) + angle_axis = normal / units.norm(normal) + angle_rotation_matrix = geometry_engine._rotation_matrix(angle_axis, theta_rotate) + d_ang = units.dot(angle_rotation_matrix, a) if __name__=="__main__": - test_coordinate_conversion() + #test_coordinate_conversion() test_run_geometry_engine() + #test_existing_coordinates() + #test_openmm_dihedral() + #_try_random_itoc() \ No newline at end of file From 806adc0e28a60ff095acf43d697e2b2dc379841c Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 13:44:38 -0500 Subject: [PATCH 51/57] fix angles --- perses/rjmc/geometry.py | 2 +- perses/tests/test_geometry_engine.py | 21 +++++++++++++++++---- 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 250735177..a44052d21 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -646,7 +646,7 @@ def _calculate_angle(self, atom_position, bond_atom_position, angle_atom_positio b = angle_atom_position - bond_atom_position a_u = a / np.linalg.norm(a) b_u = b / np.linalg.norm(b) - cos_theta = np.dot(a_u, b_u) + cos_theta = np.dot(-a_u, b_u) if cos_theta > 1.0: cos_theta = 1.0 elif cos_theta < -1.0: diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 563f258fb..03bbc9f9c 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -12,6 +12,19 @@ kB = units.BOLTZMANN_CONSTANT_kB * units.AVOGADRO_CONSTANT_NA beta = 1.0/ (300.0*units.kelvin*kB) + +def generate_molecule_from_smiles(smiles): + """ + Generate oemol with geometry from smiles + """ + mol = oechem.OEMol() + oechem.OESmilesToMol(mol, smiles) + oechem.OEAddExplicitHydrogens(mol) + omega = oeomega.OEOmega() + omega.SetMaxConfs(1) + omega(mol) + return mol + def generate_initial_molecule(iupac_name): """ Generate an oemol with a geometry @@ -52,7 +65,7 @@ def align_molecules(mol1, mol2): atomexpr = oechem.OEExprOpts_AtomicNumber bondexpr = 0 mcs.Init(mol1, atomexpr, bondexpr) - mcs.SetMCSFunc(oechem.OEMCSMaxBondsCompleteCycles()) + mcs.SetMCSFunc(oechem.OEMCSMaxAtomsCompleteCycles()) unique = True match = [m for m in mcs.Match(mol2, unique)][0] new_to_old_atom_mapping = {} @@ -68,8 +81,8 @@ def test_run_geometry_engine(): Run the geometry engine a few times to make sure that it actually runs without exceptions. Convert n-pentane to 2-methylpentane """ - molecule_name_1 = 'propane' - molecule_name_2 = 'butane' + molecule_name_1 = 'pentane' + molecule_name_2 = '2-methylpentane' molecule1 = generate_initial_molecule(molecule_name_1) molecule2 = generate_initial_molecule(molecule_name_2) @@ -85,7 +98,7 @@ def test_run_geometry_engine(): old_positions=pos1, logp_proposal=0.0, new_to_old_atom_map=new_to_old_atom_mapping, metadata={'test':0.0}) sm_top_proposal._beta = beta geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) - test_pdb_file = open("pbutane_4.pdb", 'w') + test_pdb_file = open("ethene12diol_3.pdb", 'w') integrator = openmm.VerletIntegrator(1*units.femtoseconds) From d956d885216b37cc9b068870987587b84fe66267 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 14:12:28 -0500 Subject: [PATCH 52/57] test angle computation --- perses/tests/test_geometry_engine.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 03bbc9f9c..21e078c63 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -282,6 +282,17 @@ def _get_internal_from_omm(atom_coords, bond_coords, angle_coords, torsion_coord return r, theta, phi +def test_angle(): + """ + Test the _calculate_angle function in the geometry engine to make sure it gets the same number as openmm + """ + import perses.rjmc.geometry as geometry + geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) + example_coordinates = units.Quantity(np.random.normal(size=[4,3])) + r, theta, phi = _get_internal_from_omm(example_coordinates[0], example_coordinates[1], example_coordinates[2], example_coordinates[3]) + theta_g = geometry_engine._calculate_angle(example_coordinates[0], example_coordinates[1], example_coordinates[2]) + assert abs(theta / theta.unit - theta_g) < 1.0e-12 + def test_rotation_matrix(): """ @@ -309,7 +320,8 @@ def test_rotation_matrix(): if __name__=="__main__": #test_coordinate_conversion() - test_run_geometry_engine() + #test_run_geometry_engine() #test_existing_coordinates() #test_openmm_dihedral() - #_try_random_itoc() \ No newline at end of file + #_try_random_itoc() + test_angle() \ No newline at end of file From 5ca5c9f6901b04279e1368f33c49244eaa695e75 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 14:13:53 -0500 Subject: [PATCH 53/57] remove plotting --- perses/rjmc/geometry.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index a44052d21..954cf8003 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -7,7 +7,6 @@ import parmed import simtk.unit as units import logging -import matplotlib.pyplot as pyplot import numpy as np import coordinate_tools @@ -681,7 +680,6 @@ def _propose_torsion(self, atom, r, theta, bond_atom, angle_atom, torsion_atom, logp, Z, q, phis = self._normalize_torsion_proposal(atom, r, theta, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=100) #choose from the set of possible torsion angles phi_idx = np.random.choice(range(len(phis)), p=np.exp(logp)) - pyplot.plot(phis, np.exp(logp)) logp = logp[phi_idx] phi = phis[phi_idx] return phi, logp @@ -718,12 +716,6 @@ def _joint_torsion_angle_proposal(self, atom, r, bond_atom, angle_atom, torsion_ n_divisions = 50 logp, Z, q, thetas, phis, xyz_grid = self._normalize_torsion_and_angle(atom, r, bond_atom, angle_atom, torsion_atom, atoms_with_positions, positions, beta, n_divisions=n_divisions) logp_flat = np.ravel(logp) - #xyz_flat = np.ravel(xyz_grid) - #import matplotlib.pyplot as plt - #from mpl_toolkits.mplot3d import Axes3D - #fig = plt.figure() - #ax = fig.add_subplot(111, projection='3d') - #ax.scatter(xyz_grid[:,:,0], xyz_grid[:,:,1], xyz_grid[:,:,2]) theta_phi_idx = np.random.choice(range(len(logp_flat)), p=np.exp(logp_flat)) theta_idx, phi_idx = np.unravel_index(theta_phi_idx, [n_divisions, n_divisions]) xyz = xyz_grid[theta_idx, phi_idx] From cba5b45e1309fb163301768c212943a6290e5d9d Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 14:35:07 -0500 Subject: [PATCH 54/57] temprorary unit fix --- perses/tests/test_geometry_engine.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 21e078c63..9faeaa2ea 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -155,9 +155,9 @@ def test_coordinate_conversion(): angle_position = units.Quantity(np.array([ 0.0829 , 0.0952 ,-0.2479]) ,unit=units.nanometers) torsion_position = units.Quantity(np.array([-0.057 , 0.0951 ,-0.1863] ) ,unit=units.nanometers) rtp, detJ = geometry_engine._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) - r = rtp[0] - theta = rtp[1] - phi = rtp[2] + r = rtp[0]*units.nanometers + theta = rtp[1]*units.radians + phi = rtp[2]*units.radians xyz, _ = geometry_engine._internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) assert np.linalg.norm(xyz-atom_position) < 1.0e-12 @@ -301,13 +301,10 @@ def test_rotation_matrix(): import perses.rjmc.geometry as geometry geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) #get 3 points: - points = units.Quantity(np.random.normal(size=[3,3]), unit=units.nanometers) + points = units.Quantity(np.random.normal(size=[4,3]), unit=units.nanometers) + + #have openmm calculate the angle between them: - #calculate the angle between them: - a = points[1] - points[0] - b = points[2] - points[1] - a_u = a / units.norm(a) - b_u = b / units.norm(b) theta_initial = units.acos(units.dot(a_u, b_u)) #now, rotate point 0 about an axis: @@ -318,10 +315,13 @@ def test_rotation_matrix(): d_ang = units.dot(angle_rotation_matrix, a) + + if __name__=="__main__": - #test_coordinate_conversion() + test_coordinate_conversion() #test_run_geometry_engine() #test_existing_coordinates() #test_openmm_dihedral() #_try_random_itoc() - test_angle() \ No newline at end of file + #test_angle() + #test_rotation_matrix() \ No newline at end of file From 390efd13d0906fa38c949a512f6d67d74b58b826 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 15:39:55 -0500 Subject: [PATCH 55/57] change name of test to make run --- perses/tests/test_geometry_engine.py | 35 +++++----------------------- 1 file changed, 6 insertions(+), 29 deletions(-) diff --git a/perses/tests/test_geometry_engine.py b/perses/tests/test_geometry_engine.py index 9faeaa2ea..3a95b9f0c 100644 --- a/perses/tests/test_geometry_engine.py +++ b/perses/tests/test_geometry_engine.py @@ -206,7 +206,7 @@ def test_openmm_dihedral(): return 0 -def _try_random_itoc(): +def test_try_random_itoc(): import perses.rjmc.geometry as geometry geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) @@ -294,34 +294,11 @@ def test_angle(): assert abs(theta / theta.unit - theta_g) < 1.0e-12 -def test_rotation_matrix(): - """ - Test the rotation matrix code to ensure that it does what we expect - """ - import perses.rjmc.geometry as geometry - geometry_engine = geometry.FFAllAngleGeometryEngine({'test': 'true'}) - #get 3 points: - points = units.Quantity(np.random.normal(size=[4,3]), unit=units.nanometers) - - #have openmm calculate the angle between them: - - theta_initial = units.acos(units.dot(a_u, b_u)) - - #now, rotate point 0 about an axis: - theta_rotate = np.pi/2.0*units.radians - normal = np.cross(a, b) - angle_axis = normal / units.norm(normal) - angle_rotation_matrix = geometry_engine._rotation_matrix(angle_axis, theta_rotate) - d_ang = units.dot(angle_rotation_matrix, a) - - - if __name__=="__main__": test_coordinate_conversion() - #test_run_geometry_engine() - #test_existing_coordinates() - #test_openmm_dihedral() - #_try_random_itoc() - #test_angle() - #test_rotation_matrix() \ No newline at end of file + test_run_geometry_engine() + test_existing_coordinates() + test_openmm_dihedral() + test_try_random_itoc() + test_angle() From d465db1bf7a818e786aac77efa8010481afd2fa5 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 15:42:01 -0500 Subject: [PATCH 56/57] calculate determinant of jacobian --- perses/rjmc/geometry.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index 954cf8003..c6304476b 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -415,7 +415,7 @@ def _rotation_matrix(self, axis, angle): def _cartesian_to_internal(self, atom_position, bond_position, angle_position, torsion_position): """ - Cartesian to internal function (hah) + Cartesian to internal function """ #ensure we have the correct units, then remove them atom_position = atom_position.in_units_of(units.nanometers)/units.nanometers @@ -426,7 +426,7 @@ def _cartesian_to_internal(self, atom_position, bond_position, angle_position, t internal_coords = coordinate_tools._cartesian_to_internal(atom_position, bond_position, angle_position, torsion_position) - return internal_coords, 0 + return internal_coords, internal_coords[0]**2*np.sin(internal_coords) def _internal_to_cartesian(self, bond_position, angle_position, torsion_position, r, theta, phi): @@ -441,7 +441,7 @@ def _internal_to_cartesian(self, bond_position, angle_position, torsion_position torsion_position = torsion_position.in_units_of(units.nanometers)/units.nanometers xyz = coordinate_tools._internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) xyz = units.Quantity(xyz, unit=units.nanometers) - return xyz, 0 + return xyz, (r/r.unit)**2*np.sin(theta/theta.unit) def _bond_logq(self, r, bond, beta): From 6d5788ebec85644d9c3020067f060cc7fbe8baf1 Mon Sep 17 00:00:00 2001 From: pgrinaway Date: Tue, 8 Dec 2015 15:47:19 -0500 Subject: [PATCH 57/57] fix unit --- perses/rjmc/geometry.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/perses/rjmc/geometry.py b/perses/rjmc/geometry.py index c6304476b..7f6df50a6 100644 --- a/perses/rjmc/geometry.py +++ b/perses/rjmc/geometry.py @@ -441,7 +441,7 @@ def _internal_to_cartesian(self, bond_position, angle_position, torsion_position torsion_position = torsion_position.in_units_of(units.nanometers)/units.nanometers xyz = coordinate_tools._internal_to_cartesian(bond_position, angle_position, torsion_position, r, theta, phi) xyz = units.Quantity(xyz, unit=units.nanometers) - return xyz, (r/r.unit)**2*np.sin(theta/theta.unit) + return xyz, r**2*np.sin(theta) def _bond_logq(self, r, bond, beta):