How to use the mdtraj.Trajectory 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 msmbuilder / msmbuilder / test / test_msm.py View on Github external
i0 = trimmed_assignments[0][0]
    if i0 == 1:
        for m in trimmed_assignments:
            m *= -1
            m += 1
    
    pairs = msm.draw_samples(trimmed_assignments, 2000)

    samples = map_drawn_samples(pairs, data)
    mu = np.mean(samples, axis=1)
    eq(mu, np.array([[0., 0., 0.0], [25., 25., 25.]]), decimal=1)

    # We should make sure we can sample from Trajectory objects too...
    # Create a fake topology with 1 atom to match our input dataset
    top = md.Topology.from_dataframe(pd.DataFrame({"serial":[0], "name":["HN"], "element":["H"], "resSeq":[1], "resName":"RES", "chainID":[0]}), bonds=np.zeros(shape=(0, 2), dtype='int'))
    trajectories = [md.Trajectory(x[:, np.newaxis], top) for x in data]  # np.newaxis reshapes the data to have a 40000 frames, 1 atom, 3 xyz

    trj_samples = map_drawn_samples(pairs, trajectories)
    mu = np.array([t.xyz.mean(0)[0] for t in trj_samples])
    eq(mu, np.array([[0., 0., 0.0], [25., 25., 25.]]), decimal=1)
github mdtraj / mdtraj / tests / test_rmsd.py View on Github external
def test_superpose_refinds():
    # make one frame far from the origin
    normal = np.random.randn(1, 100, 3)
    normal_xyz = normal.copy()

    flipped = np.zeros_like(normal)
    flipped[:, :50, :] = normal[:, 50:, :]
    flipped[:, 50:, :] = normal[:, :50, :]
    flipped_xyz = flipped.copy()

    normal = md.Trajectory(xyz=normal, topology=None)
    flipped = md.Trajectory(xyz=flipped, topology=None)

    normal.superpose(flipped, atom_indices=np.arange(0, 50), ref_atom_indices=np.arange(50, 100))
    eq(normal.xyz, normal_xyz)

    flipped.superpose(normal, atom_indices=np.arange(50, 100), ref_atom_indices=np.arange(0, 50))
    eq(flipped.xyz, flipped_xyz)

    normal.superpose(flipped)
    assert not np.allclose(normal.xyz, normal_xyz)
github mdtraj / mdtraj / tests / test_distance.py View on Github external
def test_closest_contact():
    box_size = np.array([3.0, 4.0, 5.0])
    traj = md.Trajectory(xyz=xyz * box_size, topology=None)
    _verify_closest_contact(traj)
    traj.unitcell_lengths = np.array([box_size for i in range(N_FRAMES)])
    traj.unitcell_angles = np.array([[90.0, 90.0, 90.0] for i in range(N_FRAMES)])
    _verify_closest_contact(traj)
    traj.unitcell_angles = np.array([[80.0, 90.0, 100.0] for i in range(N_FRAMES)])
    _verify_closest_contact(traj)
github markovmodel / molPX / projX / generate.py View on Github external
* What the dictionary actually contains

            * ``paths_dict[idxs][type_of_path]["proj"]`` : ndarray of shape (n_points, proj_dim) with the
            coordinates of the projection along the path

            * paths_dict[idxs][type_of_path]["proj"] : geometries along the path

    idata :
        list of ndarrays with the the data in  :obj:`projected_trajectories`
    """

    if not isinstance(MD_trajectories, list):
        MD_trajectories = [MD_trajectories]

    if isinstance(MD_trajectories[0], _md.Trajectory):
        src = MD_trajectories
    else:
        src = _source(MD_trajectories, top=MD_top)


    idata = _data_from_input(projected_trajectories)
    assert  len(MD_trajectories) == len(idata), "Mismatch between the number of " \
                                                "MD_trajectories and projected_trajectories: %u vs %u"%(len(MD_trajectories), len(idata))

    #TODO: assert total_n_frames (strided) coincies with the n_frames in data
    # What's the hightest dimensionlatiy that the input data allows?
    input_dim = idata[0].shape[1]
    if proj_idxs is None:
       proj_idxs = _np.arange(n_projs)
    else:
        if isinstance(proj_idxs, int):
github ADicksonLab / wepy / wepy / resampling / wexplore2_old.py View on Github external
def _maketraj(self, positions):
        Newxyz = np.zeros((1, self.ref.n_atoms, 3))
        
        for i in range(len(positions)):
            Newxyz[0,i,:] = ([positions[i]._value[0], positions[i]._value[1],
                                                        positions[i]._value[2]])
        return mdj.Trajectory(Newxyz, self.ref.topology)
github markovmodel / molPX / projX / bmutils.py View on Github external
atom_indices = _np.arange(now.n_atoms)
    else:
       atom_indices = selection

    # For the list of candidates, extract the closest one
    history = now
    for ii, cands in enumerate(path_of_candidates):
        closest_to_now = _np.argmin(_md.rmsd(cands, now, atom_indices=atom_indices))
        path_out.append(closest_to_now)
        #print("choose frame %u from %u cands"%(path_out[-1], len(cands)))
        now = cands[closest_to_now]
        history = history.join(now)
        if history_aware:
           history.superpose(history, atom_indices=atom_indices)
           xyz = history.xyz.mean(0)
           now = _md.Trajectory(xyz, history.top)
    return path_out
github markovmodel / PyEMMA / pyemma / msm / util / mapping.py View on Github external
states = np.unique(np.hstack(([np.unique(disctraj) for disctraj in disctrajs])))
    states = np.setdiff1d(states, [-1])  # exclude invalid states
    cluster = [None] * (np.max(states) + 1)
    for disctraj, traj, i in zip(disctrajs, trajs, xrange(len(trajs))):
        assert len(disctraj) == traj.n_frames, 'Length of disctraj[%d] doesn\'t match number of frames in traj[%d].' % (
            i, i)
        for s in states:
            match = (disctraj == s)
            if np.count_nonzero(match) > 0:
                if cluster[s] is None:
                    cluster[s] = traj.xyz[match, :, :]
                else:
                    cluster[s] = np.concatenate((cluster[s], traj.xyz[match, :, :]), axis=0)
    for i in xrange(len(cluster)):
        if not cluster[i] is None:
            cluster[i] = md.Trajectory(cluster[i], trajs[0].topology)
    return cluster
github deepchem / deepchem / deepchem / feat / atomic_coordinates.py View on Github external
def compute_neighbor_list(coords, neighbor_cutoff, max_num_neighbors,
                          periodic_box_size):
  """Computes a neighbor list from atom coordinates."""
  N = coords.shape[0]
  import mdtraj
  traj = mdtraj.Trajectory(coords.reshape((1, N, 3)), None)
  box_size = None
  if periodic_box_size is not None:
    box_size = np.array(periodic_box_size)
    traj.unitcell_vectors = np.array(
        [[[box_size[0], 0, 0], [0, box_size[1], 0], [0, 0, box_size[2]]]],
        dtype=np.float32)
  neighbors = mdtraj.geometry.compute_neighborlist(traj, neighbor_cutoff)
  neighbor_list = {}
  for i in range(N):
    if max_num_neighbors is not None and len(neighbors[i]) > max_num_neighbors:
      delta = coords[i] - coords.take(neighbors[i], axis=0)
      if box_size is not None:
        delta -= np.round(delta / box_size) * box_size
      dist = np.linalg.norm(delta, axis=1)
      sorted_neighbors = list(zip(dist, neighbors[i]))
      sorted_neighbors.sort()
github ADicksonLab / wepy / src / wepy / boundary_conditions / receptor.py View on Github external
Parameters
        ----------
        walker : object implementing the Walker interface

        Returns
        -------
        min_distance : float

        """

        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