How to use the mdtraj.Topology 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 choderalab / openmmtools / openmmtools / testsystems.py View on Github external
import mdtraj as md
        import pandas as pd
    except ImportError as e:
        print("Error: generate_dummy_trajectory() requires mdtraj and pandas!")
        raise(e)

    n_atoms = len(xyz)
    data = []

    for i in range(n_atoms):
        data.append(dict(serial=i, name="H", element="H", resSeq=i + 1, resName="UNK", chainID=0))

    data = pd.DataFrame(data)
    unitcell_lengths = box * np.ones((1, 3))
    unitcell_angles = 90 * np.ones((1, 3))
    top = md.Topology.from_dataframe(data, np.zeros((0, 2), dtype='int'))
    traj = md.Trajectory(xyz, top, unitcell_lengths=unitcell_lengths, unitcell_angles=unitcell_angles)

    return traj
github ADicksonLab / wepy / examples / Lennard_Jones_Pair / gpu_run.py View on Github external
get_state_kwargs = dict(GET_STATE_KWARG_DEFAULTS)
    init_sim_state = context.getState(**get_state_kwargs)

    thermostat = omm.AndersenThermostat(300.0 * unit.kelvin, 1/unit.picosecond)
    barostat = omm.MonteCarloBarostat(1.0*unit.atmosphere, 300.0*unit.kelvin, 50)

    runner = OpenMMRunner(test_sys.system, test_sys.topology, integrator, platform='CUDA')

    num_walkers = 10
    init_weight = 1.0 / num_walkers

    init_walkers = [Walker(OpenMMState(init_sim_state), init_weight) for i in range(num_walkers)]


    # the mdtraj here is needed for the distance function
    mdtraj_topology = mdj.Topology.from_openmm(test_sys.topology)

    # make a distance object which can be used to compute the distance
    # between two walkers, for our scorer class
    distance = PairDistance()

    # make a WExplore2 resampler with default parameters and our
    # distance metric
    resampler = WExplore1Resampler(max_region_sizes=(0.5, 0.5, 0.5, 0.5),
                                   max_n_regions=(10, 10, 10, 10),
                                   distance=distance,
                                   pmin=1e-12, pmax=0.5)

    ubc = UnbindingBC(cutoff_distance=2.0,
                      initial_state=init_walkers[0].state,
                      topology=mdtraj_topology,
                      ligand_idxs=np.array(test_sys.ligand_indices),
github choderalab / perses / perses / app / relative_point_mutation_setup.py View on Github external
protein_pdbfile = open(protein_filename, 'r')
        protein_pdb = app.PDBFile(protein_pdbfile)
        protein_pdbfile.close()
        protein_positions, protein_topology, protein_md_topology = protein_pdb.positions, protein_pdb.topology, md.Topology.from_openmm(protein_pdb.topology)
        protein_topology = protein_md_topology.to_openmm()
        protein_n_atoms = protein_md_topology.n_atoms

        # Load the ligand, if present
        molecules = []
        if ligand_file:
            if ligand_file.endswith('.sdf'): # small molecule
                ligand_mol = createOEMolFromSDF(ligand_file, index=ligand_index)
                ligand_mol = generate_unique_atom_names(ligand_mol)
                molecules.append(Molecule.from_openeye(ligand_mol, allow_undefined_stereo=False))
                ligand_positions, ligand_topology = extractPositionsFromOEMol(ligand_mol),  forcefield_generators.generateTopologyFromOEMol(ligand_mol)
                ligand_md_topology = md.Topology.from_openmm(ligand_topology)
                ligand_n_atoms = ligand_md_topology.n_atoms

            elif ligand_file.endswith('pdb'): # protein
                ligand_pdbfile = open(ligand_file, 'r')
                ligand_pdb = app.PDBFile(ligand_pdbfile)
                ligand_pdbfile.close()
                ligand_positions, ligand_topology, ligand_md_topology = ligand_pdb.positions, ligand_pdb.topology, md.Topology.from_openmm(
                    ligand_pdb.topology)
                ligand_n_atoms = ligand_md_topology.n_atoms

            else:
                _logger.warning(f'ligand filetype not recognised. Please provide a path to a .pdb or .sdf file')
                return

            # Now create a complex
            complex_md_topology = protein_md_topology.join(ligand_md_topology)
github ADicksonLab / wepy / mocks / init_run_groups.py View on Github external
hdf5_filename = 'tmp.wepy.h5'
wepy_h5 = WepyHDF5(hdf5_filename, mode='w', topology=top_json)


# make the test system from openmmtools
test_sys = LennardJonesPair()

# integrator
TEMPERATURE= 300.0*unit.kelvin
FRICTION_COEFFICIENT = 1/unit.picosecond
STEP_SIZE = 0.002*unit.picoseconds
integrator = omm.LangevinIntegrator(TEMPERATURE, FRICTION_COEFFICIENT, STEP_SIZE)

# the mdtraj here is needed for unbininding BC
mdtraj_topology = mdj.Topology.from_openmm(test_sys.topology)

# the initial state for the BC
context = omm.Context(test_sys.system, copy(integrator))
context.setPositions(test_sys.positions)
get_state_kwargs = dict(GET_STATE_KWARG_DEFAULTS)
init_sim_state = context.getState(**get_state_kwargs)
init_state = OpenMMState(init_sim_state)

# initialize the unbinding boundary conditions
ubc = UnbindingBC(cutoff_distance=0.5,
                  initial_state=init_state,
                  topology=mdtraj_topology,
                  ligand_idxs=np.array(test_sys.ligand_indices),
                  binding_site_idxs=np.array(test_sys.receptor_indices))
github choderalab / openmmtools / openmmtools / forcefactories.py View on Github external
def restrain_atoms_by_dsl(thermodynamic_state, sampler_state, topology, atoms_dsl, **kwargs):
    # Make sure the topology is an MDTraj topology.
    if isinstance(topology, mdtraj.Topology):
        mdtraj_topology = topology
    else:
        mdtraj_topology = mdtraj.Topology.from_openmm(topology)

    # Determine indices of the atoms to restrain.
    restrained_atoms = mdtraj_topology.select(atoms_dsl).tolist()
    restrain_atoms(thermodynamic_state, sampler_state, restrained_atoms, **kwargs)
github ADicksonLab / wepy / examples / Lennard_Jones_Pair / wexplore2.py View on Github external
get_state_kwargs = dict(GET_STATE_KWARG_DEFAULTS)
    init_sim_state = context.getState(**get_state_kwargs)

    thermostat = omm.AndersenThermostat(300.0 * unit.kelvin, 1/unit.picosecond)
    barostat = omm.MonteCarloBarostat(1.0*unit.atmosphere, 300.0*unit.kelvin, 50)

    runner = OpenMMRunner(test_sys.system, test_sys.topology, integrator, platform='Reference')

    num_walkers = 10
    init_weight = 1.0 / num_walkers

    init_walkers = [Walker(OpenMMState(init_sim_state), init_weight) for i in range(num_walkers)]

    # the mdtraj here is needed for the distance function
    mdtraj_topology = mdj.Topology.from_openmm(test_sys.topology)

    # make a distance object which can be used to compute the distance
    # between two walkers, for our scorer class
    distance = PairDistance()

    # we need a scorer class to perform the all-to-all distances
    # using our distance class
    scorer = AllToAllScorer(distance=distance)

    # make a WExplore2 resampler with default parameters and our
    # distance metric
    resampler = WExplore2Resampler(scorer=scorer,
                                   pmax=0.5)

    ubc = UnbindingBC(cutoff_distance=1.5,
                      initial_state=init_walkers[0].state,
github mdtraj / mdtraj / mdtraj / formats / lh5.py View on Github external
def _topology_from_arrays(AtomID, AtomNames, ChainID, ResidueID, ResidueNames):
    """Build topology object from the arrays stored in the lh5 file"""
    # Delayed import due to wacky recursive imports in compatibilty
    from mdtraj import Topology
    topology = Topology()

    # assert that the ChainID is just an array of empty strings, which appears
    # to be the case in our test systems for this legacy format
    if not np.all(chainid == '' for chainid in ChainID):
        raise NotImplementedError('Im not prepared to parse multiple chains')
    chain0 = topology.add_chain()

    # register the residues
    registered_residues = {}
    for i in np.argsort(ResidueID):
        residue_name = ResidueNames[i]
        if not isinstance(residue_name, basestring):
            residue_name = residue_name.decode()
        if ResidueID[i] not in registered_residues:
            res = topology.add_residue(residue_name, chain0)
            registered_residues[ResidueID[i]] = res
github markovmodel / PyEMMA / pyemma / coordinates / data / featurization / featurizer.py View on Github external
def topologyfile(self, topfile):
        self._topologyfile = topfile
        if isinstance(topfile, str):
            self.topology = load_topology_cached(topfile)
            self._topologyfile = topfile
        elif isinstance(topfile, mdtraj.Topology):
            self.topology = topfile
        elif isinstance(topfile, mdtraj.Trajectory):
            self.topology = topfile.topology
        else:
            raise ValueError("no valid topfile arg: type was %s, "
                             "but only string or mdtraj.Topology allowed." % type(topfile))
github choderalab / perses / perses / app / experiments.py View on Github external
old_ligand_topology = old_complex.subset(old_complex.select("resname == 'MOL' "))
        new_ligand_topology = new_complex.subset(new_complex.select("resname == 'MOL' "))

        # solvate the old ligand topology:
        old_solvated_topology, old_solvated_positions, old_solvated_system = self._solvate(old_ligand_topology.to_openmm(), old_ligand_positions)

        old_solvated_md_topology = md.Topology.from_openmm(old_solvated_topology)

        # now remove the old ligand, leaving only the solvent
        solvent_only_topology = old_solvated_md_topology.subset(old_solvated_md_topology.select("not resname MOL"))
        # append the solvent to the new ligand-only topology:
        new_solvated_ligand_md_topology = new_ligand_topology.join(solvent_only_topology)
        nsl, b = new_solvated_ligand_md_topology.to_dataframe()

        # dirty hack because new_solvated_ligand_md_topology.to_openmm() was throwing bond topology error
        new_solvated_ligand_md_topology = md.Topology.from_dataframe(nsl, b)

        new_solvated_ligand_omm_topology = new_solvated_ligand_md_topology.to_openmm()
        new_solvated_ligand_omm_topology.setPeriodicBoxVectors(old_solvated_topology.getPeriodicBoxVectors())

        # create the new ligand system:
        new_solvated_system = self.system_generator.create_system(new_solvated_ligand_omm_topology)

        new_to_old_atom_map = {atom_map[x] - new_mol_start_index: x - old_mol_start_index for x in
                               old_complex.select("resname == 'MOL' ") if x in atom_map.keys()}

        # adjust the atom map to account for the presence of solvent degrees of freedom:
        # By design, all atoms after the ligands are water, and should be mapped.
        n_water_atoms = solvent_only_topology.to_openmm().getNumAtoms()
        for i in range(n_water_atoms):
            new_to_old_atom_map[new_mol_len + i] = old_mol_len + i