How to use the mdtraj.compute_distances function in mdtraj

To help you get started, we’ve selected a few mdtraj examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

github mdtraj / mdtraj / tests / test_distance.py View on Github external
def _run_amber_traj(traj, ext_ref):
    # Test triclinic case where simple approach in Tuckerman text does not
    # always work
    distopt = md.compute_distances(traj, [[0, 9999]], opt=True)
    distslw = md.compute_distances(traj, [[0, 9999]], opt=False)
    dispopt = md.compute_displacements(traj, [[0, 9999]], opt=True)
    dispslw = md.compute_displacements(traj, [[0, 9999]], opt=False)

    eq(distopt, distslw, decimal=5)
    eq(dispopt, dispslw, decimal=5)

    assert_allclose(distopt.flatten(), ext_ref, atol=2e-5)

    # Make sure distances from displacements are the same
    eq(np.sqrt((dispopt.squeeze() ** 2).sum(axis=1)), distopt.squeeze())
    eq(np.sqrt((dispslw.squeeze() ** 2).sum(axis=1)), distslw.squeeze())
    eq(dispopt, dispslw, decimal=5)
github mdtraj / mdtraj / tests / test_distance.py View on Github external
def _run_amber_traj(traj, ext_ref):
    # Test triclinic case where simple approach in Tuckerman text does not
    # always work
    distopt = md.compute_distances(traj, [[0, 9999]], opt=True)
    distslw = md.compute_distances(traj, [[0, 9999]], opt=False)
    dispopt = md.compute_displacements(traj, [[0, 9999]], opt=True)
    dispslw = md.compute_displacements(traj, [[0, 9999]], opt=False)

    eq(distopt, distslw, decimal=5)
    eq(dispopt, dispslw, decimal=5)

    assert_allclose(distopt.flatten(), ext_ref, atol=2e-5)

    # Make sure distances from displacements are the same
    eq(np.sqrt((dispopt.squeeze() ** 2).sum(axis=1)), distopt.squeeze())
    eq(np.sqrt((dispslw.squeeze() ** 2).sum(axis=1)), distslw.squeeze())
    eq(dispopt, dispslw, decimal=5)
github ADicksonLab / wepy / src / wepy / boundary_conditions / receptor.py View on Github external
"""

        cell_lengths, cell_angles = box_vectors_to_lengths_angles(walker.state['box_vectors'])

        t2 = time.time()
        # make a traj out of it so we can calculate distances through
        # the periodic boundary conditions
        walker_traj = mdj.Trajectory(walker.state['positions'],
                                     topology=self._mdj_top,
                                     unitcell_lengths=cell_lengths,
                                     unitcell_angles=cell_angles)

        t3 = time.time()
        # calculate the distances through periodic boundary conditions
        # and get hte minimum distance
        min_distance = np.min(mdj.compute_distances(walker_traj,
                                                    it.product(self.ligand_idxs,
                                                               self.receptor_idxs),
                                                    periodic=self._periodic)
        )
        t4 = time.time()
        logging.info("Make a traj: {0}; Calc dists: {1}".format(t3-t2,t4-t3))

        return min_distance
github mdtraj / mdtraj / mdtraj / scripts / mdinspect.py View on Github external
self.section("Bond Check (without PBCs)")

        radii = []
        pairs = []
        for (a, b) in self.topology.bonds:
            try:
                radsum = COVALENT_RADII[a.element.symbol] + COVALENT_RADII[b.element.symbol]
            except KeyError:
                raise NotImplementedError("I don't have radii information for all of your atoms")
            radii.append(radsum)
            pairs.append((a.index, b.index))

        radii = np.array(radii)
        pairs = np.array(pairs)

        distances = md.compute_distances(self.t, pairs, periodic=False)
        low, high = self.bond_low * radii, self.bond_high * radii
        extreme = np.logical_or(distances < low, distances > high)

        if np.any(extreme):
            frames, bonds = np.nonzero(extreme)
            frame, bond = frames[0], bonds[0]
            a1 = self.topology.atom(pairs[bond][0])
            a2 = self.topology.atom(pairs[bond][0])

            self.log('error: atoms (%s) and (%s) are bonded according to the topology ' % (a1, a2))
            self.log('but they are a distance of %.3f nm apart in frame %d' % (distances[frame, bond], frame))
        else:
            self.log("All good.")
github bowman-lab / enspara / enspara / analysis / cdfs.py View on Github external
equivalent atoms, respectively.
    min_dist : boolean, optional, default: True
        Optionally reports the minimum distance between equivalent
        atoms. If False, will report the maximum distance.

    Returns
    ----------
    dists : array, shape [trj_length, n_pairs]
        Distances between specified pairs of atoms
    """
    top = trj.topology
    apairs = _convert_atom_names(top, apairs)
    dists = []
    for num in range(len(apairs)):
        dpairs = itertools.product(apairs[num][0], apairs[num][1])
        dists_tmp = md.compute_distances(trj, atom_pairs=dpairs)
        if min_dist:
            dists.append(dists_tmp.min(axis=1))
        else:
            dists.append(dists_tmp.max(axis=1))
    return np.array(dists).T
github dwhswenson / contact_map / contact_map / concurrence.py View on Github external
def __init__(self, trajectory, atom_contacts, cutoff=0.45):
        atom_contacts = _regularize_contact_input(atom_contacts, "atom")
        atom_pairs = [[contact[0][0].index, contact[0][1].index]
                      for contact in atom_contacts]
        labels = [str(contact[0]) for contact in atom_contacts]
        distances = md.compute_distances(trajectory, atom_pairs=atom_pairs)
        vector_f = np.vectorize(lambda d: d < cutoff)
        # transpose because distances is ndarray shape (n_frames,
        # n_contacts); values should be list shape (n_contacts, n_frames)
        values = vector_f(distances).T.tolist()
        super(AtomContactConcurrence, self).__init__(values=values,
                                                     labels=labels)
github getcontacts / getcontacts / v1 / contact_calc / vanderwaal.py View on Github external
"""
	filtered_candidate_pairs = {} # outer dictionary key = time, value = list of tuples (ca_1, ca_2)
	ca_keys, ca_indices = [], []
	for ca_atom in chain_top.atoms_by_name('CA'):
		ca_keys.append(ca_atom)
		ca_indices.append(ca_atom.index)
	ca_pair_keys, ca_pair_indices = [], []
	n = len(ca_indices)
	if(n == 0): return 
	for i1 in range(n):
		for i2 in range(i1 + 1, n):
			ca_pair_keys.append((ca_keys[i1], ca_keys[i2]))
			ca_pair_indices.append([ca_indices[i1], ca_indices[i2]])
	ca_pair_indices = np.array(ca_pair_indices)
	if(len(ca_pair_indices) == 0): return 
	pairDistances = md.compute_distances(traj, ca_pair_indices)
	for time in range(len(pairDistances)):
		t_distances = pairDistances[time]
		potential_indices = [i for i in range(len(t_distances)) if t_distances[i] <= ALPHA_CARBON_DIST_CUTOFF]
		filtered_candidate_pairs[time] = [ca_pair_keys[i] for i in potential_indices]
	filtered_candidate_pairs = fillTimeGaps(range(len(pairDistances)), filtered_candidate_pairs)
	return filtered_candidate_pairs
github getcontacts / getcontacts / v1 / contact_calc / aromatic.py View on Github external
Perform soft distance cutoff
	"""
	filtered_candidate_pairs = {}
	pairKeys = []
	atomPairs = []
	tempDict = {}
	for k1 in cand_aromatic_dict.keys():
		for k2 in cand_aromatic_dict.keys():
			if(k1 != k2 and (k1, k2) not in tempDict and (k2, k1) not in tempDict):
				tempDict[(k1, k2)] = 1
				aromatic1 = cand_aromatic_dict[k1]
				aromatic2 = cand_aromatic_dict[k2]
				generatePairKeys(pairKeys, atomPairs, k1, aromatic1, k2, aromatic2)
	# calculate atom pair distances
	atomPairs = np.array(atomPairs)
	pairDistances = md.compute_distances(traj, atomPairs)
	for time in range(len(pairDistances)):
		t_distances = pairDistances[time]
		for i in range(0, len(t_distances), 9):
			arom_arom_distances = t_distances[i:i+9]
			if(min(arom_arom_distances) <= soft_cutoff_dist):
				arom1_key = pairKeys[i][0][0]
				arom2_key = pairKeys[i][1][0]
				if(time not in filtered_candidate_pairs.keys()):
					filtered_candidate_pairs[time] = {arom1_key:[arom2_key]}
				else:
					if(arom1_key not in filtered_candidate_pairs[time].keys()):
						filtered_candidate_pairs[time][arom1_key] = [arom2_key]
					else:
						filtered_candidate_pairs[time][arom1_key].append(arom2_key)
	filtered_candidate_pairs = fillTimeGaps(range(len(pairDistances)), filtered_candidate_pairs)
	return filtered_candidate_pairs