How to use the mdtraj.Topology.from_openmm function in mdtraj

old_mol_start_index, old_mol_len = self._proposal_engine._find_mol_start_index(old_complex.to_openmm())
        new_mol_start_index, new_mol_len = self._proposal_engine._find_mol_start_index(new_complex.to_openmm())

        old_pos = unit.Quantity(np.zeros([len(old_positions), 3]), unit=unit.nanometers)
        old_pos[:, :] = old_positions
        old_ligand_positions = old_pos[old_mol_start_index:(old_mol_start_index + old_mol_len), :]

        # subset the topologies:
        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_system(
            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)

def _remove_ligand_overlap(self):
        """
        Translate the ligand so that it is not overlapping with the receptor.

        Description
        -----------
        The bounding sphere of the ligand and receptor are computed, and the ligand translated along the x-direction to not overlap with the protein.

        TODO:
        * This is not guaranteed to work for periodic systems.

        """
        import mdtraj as md

        # Create an mdtraj Topology of the complex from OpenMM Topology object.
        mdtraj_complex_topology = md.Topology.from_openmm(self._topology)

        # Create an mdtraj instance of the complex.
        # TODO: Fix this when mdtraj can deal with OpenMM units.
        positions_in_mdtraj_format = np.array(self._positions / units.nanometers)
        mdtraj_complex = md.Trajectory(positions_in_mdtraj_format, mdtraj_complex_topology)

        # Compute centers of receptor and ligand.
        receptor_center = mdtraj_complex.xyz[0][self._receptor_atoms,:].mean(0)
        ligand_center = mdtraj_complex.xyz[0][self._ligand_atoms,:].mean(0)

        # Count number of receptor and ligand atoms.
        nreceptor_atoms = len(self._receptor_atoms)
        nligand_atoms = len(self._ligand_atoms)

        # Compute max radii of receptor and ligand.
        receptor_radius = (((mdtraj_complex.xyz[0][self._receptor_atoms,:] - np.tile(receptor_center, (nreceptor_atoms,1))) ** 2.).sum(1) ** 0.5).max()

-------        
        topology : mdtraj.Topology
            the topology

        Notes
        -----
        This is taken from the configuration of the first frame. Otherwise there is still un ugly fall-back to look
        for an openmm.Simulation object in Trajectory.simulator. and construct an mdtraj.Topology from this.
        """        

        if len(self) > 0 and self[0].topology is not None:
            # if no topology is defined
            topology = self[0].topology
        else:
            # TODO: kind of ugly fall-back, but helps for now
            topology = md.Topology.from_openmm(Trajectory.simulator.simulation.topology)
        
        if self.atom_indices is not None:
            topology = topology.subset(self.atom_indices)       
        
        return topology

"""
        Create an mdtraj topology corresponding to the hybrid system.
        This is purely for writing out trajectories--it is not expected to be parameterized.

        Returns
        -------
        hybrid_topology : mdtraj.Topology
        """
        #first, make an md.Topology of the old system:
        old_topology = md.Topology.from_openmm(self._topology_proposal.old_topology)

        #now make a copy for the hybrid:
        hybrid_topology = copy.deepcopy(old_topology)

        #next, make a topology of the new system:
        new_topology = md.Topology.from_openmm(self._topology_proposal.new_topology)

        added_atoms = dict()

        #get the core atoms in the new index system (as opposed to the hybrid index system). We will need this later
        core_atoms_new_indices = {self._hybrid_to_new_map[core_atom] for core_atom in self._atom_classes['core_atoms']}

        #now, add each unique new atom to the topology (this is the same order as the system)
        for particle_idx in self._topology_proposal.unique_new_atoms:
            new_system_atom = new_topology.atom(particle_idx)

            #first, we get the residue in the new system associated with this atom
            new_system_residue = new_system_atom.residue

            #next, we have to enumerate the other atoms in that residue to find mapped atoms
            new_system_atom_set = {atom.index for atom in new_system_residue.atoms}

if protein_pdb_filename:
            self._protein_pdb_filename = protein_pdb_filename
            protein_pdbfile = open(self._protein_pdb_filename, 'r')
            pdb_file = app.PDBFile(protein_pdbfile)
            protein_pdbfile.close()
            self._receptor_positions_old = pdb_file.positions
            self._receptor_topology_old = pdb_file.topology
            self._receptor_md_topology_old = md.Topology.from_openmm(self._receptor_topology_old)

        elif receptor_mol2_filename:
            self._receptor_mol2_filename = receptor_mol2_filename
            self._receptor_mol = createOEMolFromSDF(self._receptor_mol2_filename)
            mol_list.append(self._receptor_mol)
            self._receptor_positions_old = extractPositionsFromOEMol(self._receptor_mol)
            self._receptor_topology_old = forcefield_generators.generateTopologyFromOEMol(self._receptor_mol)
            self._receptor_md_topology_old = md.Topology.from_openmm(self._receptor_topology_old)
        else:
            raise ValueError("You need to provide either a protein pdb or a receptor mol2 to run a complex simulation.")

        self._complex_md_topology_old = self._receptor_md_topology_old.join(self._ligand_md_topology_old)

        n_atoms_spectators = 0
        if self._spectator_filenames:
            for i, spectator_topology in enumerate(self._spectator_md_topologies,1):
                _logger.debug(f'Appending spectator number {i} to complex topology')
                self._complex_md_topology_old = self._complex_md_topology_old.join(spectator_topology)
                n_atoms_spectators += spectator_topology.n_atoms
        self._complex_topology_old = self._complex_md_topology_old.to_openmm()

        n_atoms_total_old = self._complex_topology_old.getNumAtoms()
        n_atoms_protein_old = self._receptor_topology_old.getNumAtoms()
        n_atoms_ligand_old = n_atoms_total_old - n_atoms_protein_old - n_atoms_spectators

old_mol_start_index, old_mol_len = self.proposal_engine._find_mol_start_index(old_complex.to_openmm())
        new_mol_start_index, new_mol_len = self.proposal_engine._find_mol_start_index(new_complex.to_openmm())

        old_pos = unit.Quantity(np.zeros([len(old_positions), 3]), unit=unit.nanometers)
        old_pos[:, :] = old_positions
        old_ligand_positions = old_pos[old_mol_start_index:(old_mol_start_index + old_mol_len), :]

        # subset the topologies:
        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)

If complex is specified, the topology proposal is recycled for all phases; else, the proposal is conducted
        in solvent and recycled for vacuum.  If the only phase is vacuum, then we generate a single proposal without recycling.

        This method requires some refactoring...perhaps use classes to generate the phases...

        Arguments
        ---------

        """
        proposals = {}
        if 'complex' in self.proposal_arguments['phases']:
            _logger.debug(f"\t\tcomplex:")
            assert self.nonbonded_method == app.PME, f"Complex phase is specified, but the nonbonded method is not {app.PME} (is currently {self.nonbonded_method})."
            complex_md_topology, complex_topology, complex_positions = self._setup_complex_phase(ligand_oemol = current_oemol,
                                                                                                 ligand_positions = current_positions,
                                                                                                 ligand_topology = md.Topology.from_openmm(current_topology))

            solvated_complex_topology, solvated_complex_positions, solvated_complex_system = self._solvate(complex_topology,
                                                                                                           complex_positions,
                                                                                                           model = self.proposal_arguments['water_model'],
                                                                                                           vacuum = False)
            solvated_complex_md_topology = md.Topology.from_openmm(solvated_complex_topology)
            complex_topology_proposal = self.proposal_engine.propose(current_system = solvated_complex_system,
                                                                     current_topology = solvated_complex_topology,
                                                                     current_mol = current_oemol,
                                                                     proposed_mol = proposed_oemol)

            proposed_solvated_complex_positions, complex_logp_proposal = self.geometry_engine.propose(complex_topology_proposal,
                                                                                                      solvated_complex_positions,
                                                                                                      self.beta)
            complex_logp_reverse = self.geometry_engine.logp_reverse(complex_topology_proposal,
                                                             proposed_solvated_complex_positions,

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()

    # 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),