def __init__(self, filename, mode='r', force_overwrite=True):
        self._closed = True   # is the file currently closed?
        self._mode = mode      # what mode were we opened in
        if StrictVersion(import_('scipy.version').short_version) < StrictVersion('0.12.0'):
            raise ImportError('MDTraj NetCDF support requires scipy>=0.12.0. '
                              'You have %s' % import_('scipy.version').short_version)
        netcdf = import_('scipy.io').netcdf_file

        if mode not in ['r', 'w']:
            raise ValueError("mode must be one of ['r', 'w']")

        if mode == 'w' and not force_overwrite and os.path.exists(filename):
            raise IOError('"%s" already exists' % filename)

        # AMBER uses the NetCDF3 format, with 64 bit encodings, which
        # for scipy.io.netcdf_file is "version=2"
        self._handle = netcdf(filename, mode=mode, version=2)
        self._closed = False

        # self._frame_index is the current frame that we're at in the
        #     file
        # self._needs_initialization indicates whether we need to set the
        #     global properties of the file. This is required before the first
        #     write operation on a new file

        if mode == 'w':
            self._frame_index = 0
            self._needs_initialization = True
        elif mode == 'r':
            self._frame_index = 0
            self._needs_initialization = False
        else:
            raise RuntimeError()

    def __init__(self, filename, mode='r', force_overwrite=True):
        self._open = False
        self.filename = filename
        self.mode = mode
        if mode == 'w' and not force_overwrite and os.path.exists(filename):
            raise IOError('"%s" already exists' % filename)
        # import tables
        self.tables = import_('tables')

        if mode == 'w':
            print("Warning: The LH5 trajectory format is deprecated.", file=sys.stderr)
            # what frame are we currently reading or writing at?
            self._frame_index = 0
            # do we need to write the header information?
            self._needs_initialization = True
            if not filename.endswith('.lh5'):
                warnings.warn('The .lh5 extension is recommended.')
        elif mode == 'r':
            self._frame_index = 0
            self._needs_initialization = False
        else:
            raise ValueError("mode must be one of ['r', 'w']")

        # Compression style of legacy MSMBuilder2 lh5 trajectory format
        compression = self.tables.Filters(
            complib='blosc', shuffle=True, complevel=1)
        self._handle = self._open_file(
            filename, mode=mode, filters=compression)
        self._open = True

def in_units_of(quantity, units_out, units_in=None):
    """Convert a quantity between unit systems

    Parameters
    ----------
    quantity : number, np.ndarray, or simtk.unit.Quantity
        quantity can either be a unitted quantity -- i.e. instance of
        simtk.unit.Quantity, or just a bare number or numpy array
    units_out : str
        A string description of the units you want out. This should look
        like "nanometers/picosecondsecond" or "nanometers**3" or whatever
    units_in : str
        If you supply a quantity that's not a simtk.unit.Quantity, you should
        tell me what units it is in. If you don't, i'm just going to echo you
        back your quantity without doing any unit checking.

    Examples
    --------
    >>> in_units_of(1*units.meter**2/units.second, 'nanometers**2/picosecond')  # doctest: +SKIP
    1000000.0
    """
    units = import_('simtk.unit')

    if quantity is None:
        return quantity

    if isinstance(quantity, units.Quantity):
        return quantity.value_in_unit(_str_to_unit(units_out))
    else:
        if units_in is None:
            return quantity
        united_quantity = units.Quantity(quantity, _str_to_unit(units_in))
        return united_quantity.value_in_unit(_str_to_unit(units_out))

def _find_chains(bond_list):
    """Given a set of bonds, find unique molecules, with the assumption that
    there are no bonds between separate chains (i.e., only INTRAmolecular
    bonds), which also implies that each atom can be in exactly one chain.
    
    Parameters
    ----------
    bond_list : list of (int, int)
        The list of bonds

    Returns
    _______
    chains : list of list of int
        List of atoms in each chain

    Notes
    -----
    This function requires the NetworkX python package.
    """
    nx = import_('networkx')
    chains = []
    bond_list = np.asarray(bond_list)
    molecules = nx.Graph()
    molecules.add_nodes_from(set(bond_list.flatten()))
    molecules.add_edges_from(bond_list)
    return list(nx.connected_components(molecules))

    def to_openmm(self):
        """Convert this topology into OpenMM topology

        Returns
        -------
        topology : simtk.openmm.app.Topology
           This topology, as an OpenMM topology
        """
        app = import_('simtk.openmm.app')

        out = app.Topology()
        atom_mapping = {}

        for chain in self.chains:
            c = out.addChain()
            for residue in chain.residues:
                r = out.addResidue(residue.name, c)
                for atom in residue.atoms:
                    a = out.addAtom(atom.name, app.Element.getBySymbol(atom.element.symbol), r)
                    atom_mapping[atom] = a

        for a1, a2 in self.bonds:
            out.addBond(atom_mapping[a1], atom_mapping[a2])

        return out

    def from_openmm(cls, value):
        """Create a mdtraj topology from an OpenMM topology

        Parameters
        ----------
        value : simtk.openmm.app.Topology
            An OpenMM topology that you wish to convert to a
            mdtraj topology.
        """
        app = import_('simtk.openmm.app')

        if not isinstance(value, app.Topology):
            raise TypeError('value must be an OpenMM Topology. '
                            'You supplied a %s' % type(value))

        out = cls()
        atom_mapping = {}

        for chain in value.chains():
            c = out.add_chain()
            for residue in chain.residues():
                r = out.add_residue(residue.name, c)
                for atom in residue.atoms():
                    if atom.element is None:
                        element = elem.virtual
                    else:
                        element = elem.get_by_symbol(atom.element.symbol)
                    a = out.add_atom(atom.name, element, r)
                    atom_mapping[atom] = a

        for a1, a2 in value.bonds():
            out.add_bond(atom_mapping[a1], atom_mapping[a2])

        return out

def entry_point():
    subparsers = parser.add_subparsers(dest="subparser_name")
    scriptfiles = {}
    argv = sys.argv[:]
    if len(argv) == 1:
        argv.append('-h')

    for scriptname in scripts.__all__:
        # get the name and first sentence of the description from each of the
        # msmbuilder commands
        with warnings.catch_warnings():
            warnings.filterwarnings("ignore")
            script = import_('msmbuilder.scripts.%s' % scriptname)
            scriptparser = getattr(script, 'parser', None)
        scriptfiles[scriptname] = script.__file__

        try:
            description = scriptparser.description
        except:
            description = scriptparser.parser.description

        # http://stackoverflow.com/a/17124446/1079728
        first_sentence = ' '.join(' '.join(re.split(r'(?<=[.:;])\s', description)[:1]).split())
        subparsers.add_parser(scriptname, help=first_sentence)

    args = parser.parse_args(argv[1:2])        
    sys.argv = argv[1:]
    getattr(scripts, args.subparser_name).entry_point()

    def to_dataframe(self):
        """Convert this topology into a pandas dataframe

        Returns
        -------
        atoms : pandas.DataFrame
            The atoms in the topology, represented as a data frame.
        bonds : np.ndarray
            The bonds in this topology, represented as an n_bonds x 2 array
            of the indices of the atoms involved in each bond.
        """
        pd = import_('pandas')
        data = []
        for atom in self.atoms:
            if atom.element is None:
                element_symbol = ""
            else:
                element_symbol = atom.element.symbol
            data.append((atom.serial, atom.name, element_symbol,
                         atom.residue.resSeq, atom.residue.name,
                         atom.residue.chain.index))

        atoms = pd.DataFrame(data, columns=["serial", "name", "element",
                                            "resSeq", "resName", "chainID"])

        bonds = np.array([(a.index, b.index) for (a, b) in self.bonds])
        return atoms, bonds

def get_dihedral_connectivity(ibonds):
    """Given the bonds, get the indices of the atoms defining all the dihedral
    angles

    Parameters
    ----------
    ibonds : np.ndarray, shape=[n_bonds, 2], dtype=int
        n_bonds x 2 array of indices, where each row is the index of two
        atom who participate in a bond.

    Returns
    -------
    idihedrals : np.ndarray, shape[n_dihedrals, 4], dtype=int
        All sets of 4 atoms A,B,C,D such that A is bonded to B, B is bonded
        to C, and C is bonded to D
    """
    nx = import_('networkx')
    graph = nx.from_edgelist(ibonds)
    n_atoms = graph.number_of_nodes()
    idihedrals = []

    # TODO: CHECK FOR DIHEDRAL ANGLES THAT ARE 180 and recover
    # conf : msmbuilder.Trajectory
    #    An msmbuilder trajectory, only the first frame will be used. This
    #    is used purely to make the check for angle(ABC) != 180.

    for a in xrange(n_atoms):
        for b in graph.neighbors(a):
            for c in filter(lambda c: c not in [a, b], graph.neighbors(b)):
                for d in filter(lambda d: d not in [a, b, c], graph.neighbors(c)):
                    idihedrals.append((a, b, c, d))

    return np.array(idihedrals)

def get_angle_connectivity(ibonds):
    """Given the bonds, get the indices of the atoms defining all the bond
    angles

    Parameters
    ----------
    ibonds : np.ndarray, shape=[n_bonds, 2], dtype=int
        n_bonds x 2 array of indices, where each row is the index of two
        atom who participate in a bond.

    Returns
    -------
    iangles : np.ndarray, shape[n_angles, 3], dtype=int
        n_angles x 3 array of indices, where each row is the index of three
        atoms m,n,o such that n is bonded to both m and o.
    """
    nx = import_('networkx')
    graph = nx.from_edgelist(ibonds)
    n_atoms = graph.number_of_nodes()
    iangles = []

    for i in xrange(n_atoms):
        for (m, n) in combinations(graph.neighbors(i), 2):
            # so now the there is a bond angle m-i-n
            iangles.append((m, i, n))

    return np.array(iangles)

def _str_to_unit(unit_string):
    """eval() based transformer that extracts a simtk.unit object
    from a string description.

    Parameters
    ----------
    unit_string : str
        string description of a unit. this may contain expressions with
        multiplication, division, powers, etc.

    Examples
    --------
    >>> type(_str_to_unit('nanometers**2/meters*gigajoules'))
    <class 'simtk.unit.unit.Unit'>
    >>> str(_str_to_unit('nanometers**2/meters*gigajoules'))
    'nanometer**2*gigajoule/meter'

    """
    units = import_('simtk.unit')
    # parse the string with the ast, and then run out unit context
    # visitor on it, which will basically change bare names like
    # "nanometers" into "unit.nanometers" and simulataniously check that
    # there's no nefarious stuff in the expression.


    node = _unit_context.visit(ast.parse(unit_string, mode='eval'))
    fixed_node = ast.fix_missing_locations(node)
    output = eval(compile(fixed_node, '<string>', mode='eval'))

    return output

    def to_openmm(self, traj=None):
        """Convert this topology into OpenMM topology

        Parameters
        ----------
        traj : MDTraj.Trajectory, optional, default=None
            If specified, use the first frame from this trajectory to
            set the unitcell information in the openmm topology.

        Returns
        -------
        topology : simtk.openmm.app.Topology
           This topology, as an OpenMM topology
        """
        app = import_('simtk.openmm.app')
        mm = import_('simtk.openmm')
        u = import_('simtk.unit')

        out = app.Topology()
        atom_mapping = {}

        for chain in self.chains:
            c = out.addChain()
            for residue in chain.residues:
                r = out.addResidue(residue.name, c)
                for atom in residue.atoms:
                    if atom.element is elem.virtual:
                        element = None
                    else:
                        element = app.Element.getBySymbol(atom.element.symbol)
                    a = out.addAtom(atom.name, element, r)
                    atom_mapping[atom] = a

        for a1, a2 in self.bonds:
            out.addBond(atom_mapping[a1], atom_mapping[a2])

        if traj is not None:
            angles = traj.unitcell_angles[0]

            if np.linalg.norm(angles - 90.0) > 1E-4:
                raise(ValueError("Unitcell angles must be 90.0 to use "
                                 "in OpenMM topology."))

            box_vectors = mm.Vec3(*traj.unitcell_lengths[0]) * u.nanometer
            out.setUnitCellDimensions(box_vectors)

        return out

def chemical_shifts_shiftx2(trj, pH=5.0, temperature=298.00):
    """Predict chemical shifts of a trajectory using ShiftX2.

    Parameters
    ----------
    trj : Trajectory
        Trajectory to predict shifts for.
    pH : float, optional, default=5.0
        pH value which gets passed to the ShiftX2 predictor.
    temperature : float, optional, default=298.00
        Temperature which gets passed to the ShiftX2 predictor.

    Returns
    -------
    results : pandas DataFrame
        Dataframe containing results, with index consisting of
        (resSeq, atom_name) pairs and columns for each frame in trj.

    Notes
    -----
    You must have ShiftX2 available on your path; see (http://www.shiftx2.ca/).

    Chemical shift prediction is for PROTEIN atoms; trajectory objects
    with ligands, solvent, ions, or other non-protein components may give
    UNKNOWN RESULTS.

    Please cite the appropriate reference below.

    References
    ----------
    .. [1] Beomsoo Han, Yifeng Liu, Simon Ginzinger, and David Wishart.
       "SHIFTX2: significantly improved protein chemical shift
       prediction." J. Biomol. NMR, 50, 1 43-57 (2011)
    """
    pd = import_('pandas')
    binary = find_executable(SHIFTX2)
    if binary is None:
        raise OSError('External command not found. Looked for %s in PATH. `chemical_shifts_shiftx2` requires the external program SHIFTX2, available at http://www.shiftx2.ca/' % ', '.join(SHIFTX2))

    results = []
    with enter_temp_directory():
        for i in range(trj.n_frames):
            trj[i].save("./trj%d.pdb" % i)
        cmd = "%s -b 'trj*.pdb' -p %.1f -t %.2f" % (binary, pH, temperature)

        return_flag = os.system(cmd)

        if return_flag != 0:
            raise(IOError("Could not successfully execute command '%s', check your ShiftX2 installation or your input trajectory." % cmd))

        for i in range(trj.n_frames):
            d = pd.read_csv("./trj%d.pdb.cs" % i)
            d.rename(columns={"NUM": "resSeq", "RES": "resName", "ATOMNAME": "name"}, inplace=True)
            d["frame"] = i
            results.append(d)

    results = pd.concat(results)
    results = results.pivot_table(rows=["resSeq", "name"], cols="frame", values="SHIFT")
    return results

def chemical_shifts_ppm(trj):
    """Predict chemical shifts of a trajectory using ppm.

    Parameters
    ----------
    trj : Trajectory
        Trajectory to predict shifts for.

    Returns
    -------
    results : pandas.DataFrame
        Dataframe containing results, with index consisting of
        (resSeq, atom_name) pairs and columns for each frame in trj.

    Notes
    -----
    You must have ppm available on your path; see
    (http://spin.ccic.ohio-state.edu/index.php/download/index).

    Chemical shift prediction is for PROTEIN atoms; trajectory objects
    with ligands, solvent, ions, or other non-protein components may give
    UNKNOWN RESULTS.

    Please cite the appropriate reference below.

    References
    ----------
    .. [1] Li, DW, and Bruschweiler, R. "PPM: a side-chain and backbone chemical
       shift predictor for the assessment of protein conformational ensembles."
       J Biomol NMR. 2012 Nov;54(3):257-65.
    """
    pd = import_('pandas')
    binary = find_executable(PPM)

    first_resSeq = trj.top.residue(0).resSeq

    if binary is None:
        raise OSError('External command not found. Looked for %s in PATH. `chemical_shifts_ppm` requires the external program PPM, available at http://spin.ccic.ohio-state.edu/index.php/download/index' % ', '.join(PPM))

    with enter_temp_directory():
        trj.save("./trj.pdb")
        cmd = "%s -pdb trj.pdb -mode detail" % binary

        return_flag = os.system(cmd)

        if return_flag != 0:
            raise(IOError("Could not successfully execute command '%s', check your PPM installation or your input trajectory." % cmd))

        d = pd.read_table("./bb_details.dat", index_col=False, header=None, sep="\s*").drop([3], axis=1)

        d = d.rename(columns={0: "resSeq", 1: "resName", 2: "name"})
        d["resSeq"] += first_resSeq - 1  # Fix bug in PPM that reindexes to 1
        d = d.drop("resName", axis=1)
        d = d.set_index(["resSeq", "name"])
        d.columns = np.arange(trj.n_frames)
        d.columns.name = "frame"

    return d

def mol2_to_dataframes(filename):
    """Convert a GAFF (or sybyl) mol2 file to a pair of pandas dataframes.

    Parameters
    ----------
    filename : str
        Name of mol2 filename

    Returns
    -------
    atoms_frame : pd.DataFrame
        DataFrame containing atom information
    bonds_frame : pd.DataFrame
        DataFrame containing bond information

    Notes
    -----
    These dataframes may contain force field information as well as the
    information necessary for constructing the coordinates and molecular
    topology.  This function has been tested for GAFF and sybyl-style
    mol2 files but has been primarily tested on GAFF mol2 files.
    This function does NOT accept multi-structure MOL2 files!!!

    See Also
    --------
    If you just need the coordinates and bonds, use load_mol2(filename)
    to get a Trajectory object.
    """
    pd = import_('pandas')
    with open(filename) as f:
        data = dict((key, list(grp)) for key, grp in itertools.groupby(f, _parse_mol2_sections))

    # Mol2 can have "status bits" at the end of the bond lines. We don't care
    # about these, but they interfere with using pd_read_table because it looks
    # like one line has too many columns. So we just regex out the offending
    # text.
    status_bit_regex = "BACKBONE|DICT|INTERRES|\|"
    data["@<TRIPOS>BOND\n"] = [re.sub(status_bit_regex, lambda _: "", s)
                               for s in data["@<TRIPOS>BOND\n"]]

    if len(data["@<TRIPOS>BOND\n"]) > 1:
        csv = StringIO()
        csv.writelines(data["@<TRIPOS>BOND\n"][1:])
        csv.seek(0)
        bonds_frame = pd.read_table(csv, names=["bond_id", "id0", "id1", "bond_type"],
            index_col=0, header=None, sep="\s*", engine='python')
    else:
        bonds_frame = None

    csv = StringIO()
    csv.writelines(data["@<TRIPOS>ATOM\n"][1:])
    csv.seek(0)
    atoms_frame = pd.read_csv(csv, sep="\s*", engine='python',  header=None)
    ncols = atoms_frame.shape[1]
    names=["serial", "name", "x", "y", "z", "atype", "code", "resName", "charge", "status"]
    atoms_frame.columns = names[:ncols]
    
    return atoms_frame, bonds_frame

    def visit_Name(self, node):
        # we want to prefix all names to look like unit.nanometers instead
        # of just "nanometers", because I don't want to import * from
        # units into this module.
        units = import_('simtk.unit')
        if not (node.id == 'units' or hasattr(units, node.id)):
            # also, let's take this opporunity to check that the node.id
            # (which supposed to be the name of the unit, like "nanometers")
            # is actually an attribute in simtk.unit
            raise ValueError('%s is not a valid unit' % node.id)

        return ast.Attribute(value=ast.Name(id='units', ctx=ast.Load()),
                             attr=node.id, ctx=ast.Load())

    def __init__(self, filename, mode='r', force_overwrite=False):
        self._closed = True
        self._mode = mode
        if StrictVersion(import_('scipy.version').short_version) < StrictVersion('0.12.0'):
            raise ImportError('MDTraj NetCDF support requires scipy>=0.12.0. '
                              'You have %s' % import_('scipy.version').short_version)
        netcdf = import_('scipy.io').netcdf_file

        if mode not in ('r', 'w'):
            raise ValueError("mode must be one of ['r', 'w']")

        if mode == 'w' and not force_overwrite and os.path.exists(filename):
            raise IOError('"%s" already exists' % filename)

        # AMBER uses the NetCDF3 format, with 64 bit encodings, which for
        # scipy.io.netcdf_file is "version=2"
        self._handle = netcdf(filename, mode=mode, version=2)
        self._closed = False
        if mode == 'w':
            self._needs_initialization = True
        elif mode == 'r':
            self._needs_initialization = False
        else:
            raise RuntimeError()

def mol2_to_dataframes(filename):
    """Convert a GAFF (or sybyl) mol2 file to a pair of pandas dataframes.


    Parameters
    ----------
    filename : str
        Name of mol2 filename

    Returns
    -------
    atoms_frame : pd.DataFrame
        DataFrame containing atom information
    bonds_frame : pd.DataFrame
        DataFrame containing bond information
    
    Notes
    -----
    These dataframes may contain force field information as well as the
    information necessary for constructing the coordinates and molecular
    topology.  This function has been tested for GAFF and sybyl-style 
    mol2 files but has been primarily tested on GAFF mol2 files. 
    This function does NOT accept multi-structure MOL2 files!!!    
    
    See Also
    --------
    If you just need the coordinates and bonds, use load_mol2(filename)
    to get a Trajectory object.
    """
    pd = import_("pandas")
    with open(filename) as f:
        data = dict((key, list(grp)) for key, grp in itertools.groupby(f, _parse_mol2_sections))

    csv = StringIO()
    csv.writelines(data["@<TRIPOS>BOND\n"][1:])
    csv.seek(0)
    bonds_frame = pd.read_table(csv, names=["bond_id", "id0", "id1", "bond_type"], index_col=0, header=None, sep="\s*")

    csv = StringIO()
    csv.writelines(data["@<TRIPOS>ATOM\n"][1:])
    csv.seek(0)
    atoms_frame = pd.read_csv(
        csv, sep="\s*", names=["serial", "name", "x", "y", "z", "atype", "code", "resName", "charge"], header=None
    )  # , usecols=range(1, 10))  # usecols not available in pandas 0.11
    return atoms_frame, bonds_frame

    def __init__(self, filename, mode="r", force_overwrite=True):
        self._closed = True  # is the file currently closed?
        self._mode = mode  # what mode were we opened in
        netcdf = import_("netCDF4")

        if mode not in ["r", "w", "a", "ws", "as"]:
            raise ValueError(
                (
                    "mode must be one of ['r', 'w', 'a', 'ws', 'as']"
                    " 'r' indicates read, 'w' indicates write, and 'a' indicates"
                    " append. 'a' and 'w' can be appended with 's', which turns "
                    " off buffering"
                )
            )

        if mode in ["w", "ws"] and not force_overwrite and os.path.exists(filename):
            raise IOError('"%s" already exists')

        # AMBER uses the NetCDF3 format, with 64 bit encodings
        self._handle = netcdf.Dataset(filename, mode=mode, format="NETCDF3_64BIT", clobber=force_overwrite)
        self._closed = False

        # self._frame_index is the current frame that we're at in the
        #     file
        # self._needs_initialization indicates whether we need to set the
        #     global properties of the file. This is required before the first
        #     write operation on a new file
        # self._n_atoms is the number of atoms in the file

        if mode in ["a", "as"]:
            self._frame_index = len(self._handle.dimensions["frame"])
            self._n_atoms = len(self._handle.dimensions["atom"])
            self._needs_initialization = False
        elif mode in ["w", "ws"]:
            self._frame_index = 0
            self._n_atoms = None
            # self._n_atoms will be set during _initialize_headers call
            self._needs_initialization = True
        elif mode == "r":
            self._frame_index = 0
            self._n_atoms = len(self._handle.dimensions["atom"])
            self._needs_initialization = False
        else:
            raise RuntimeError()

    def __init__(self, filename, mode="r", force_overwrite=True, compression="zlib"):
        self._open = False  # is the file handle currently open?
        self.mode = mode  # the mode in which the file was opened?

        if not mode in ["r", "w", "a"]:
            raise ValueError("mode must be one of ['r', 'w', 'a']")

        if mode == "w" and not force_overwrite and os.path.exists(filename):
            raise IOError('"%s" already exists' % filename)

        # import tables
        self.tables = import_("tables")

        if compression == "zlib":
            compression = self.tables.Filters(complib="zlib", shuffle=True, complevel=1)
        elif compression is None:
            compression = None
        else:
            raise ValueError('compression must be either "zlib" or None')

        self._handle = self._open_file(filename, mode=mode, filters=compression)
        self._open = True

        if mode == "w":
            # what frame are we currently reading or writing at?
            self._frame_index = 0
            # do we need to write the header information?
            self._needs_initialization = True
            if not filename.endswith(".h5"):
                warnings.warn("The .h5 extension is recommended.")

        elif mode == "a":
            try:
                self._frame_index = len(self._handle.root.coordinates)
                self._needs_initialization = False
            except self.tables.NoSuchNodeError:
                self._frame_index = 0
                self._needs_initialization = True
        elif mode == "r":
            self._frame_index = 0
            self._needs_initialization = False