def __init__(self, topology):
        self._ref_topology = topology.copy()

        self._N_coeff = np.array([0.48318, 0.70328, -0.18643])  # CA_i-1, CA_i, O_i-1
        self._C_coeff = np.array([0.44365, 0.23520, 0.32115])   # CA_i, CA_i+1, O_i
        self._H_coeff = np.array([0.84100, 0.89296, -0.73389])  # CA_i-1, CA_i, O_i-1

        newTopology = Topology()
        CACBO_idxs = []; N_idxs = []; C_idxs = []; H_idxs = []
        res_idx = 1
        atm_idx = 0
        prev_ca = None
        prev_o = None
        for chain in topology._chains:
            newChain = newTopology.add_chain()
            for residue in chain._residues:
                newTopology, atm_idx, res_idx = self._add_residue(newTopology, 
                        newChain, residue, chain, res_idx, atm_idx, 
                        N_idxs, C_idxs, H_idxs, CACBO_idxs, prev_ca, prev_o)

        self._CACBO_idxs = np.array(CACBO_idxs)
        self._N_idxs = np.array(N_idxs)
        self._C_idxs = np.array(C_idxs)
        self._H_idxs = np.array(H_idxs)
        self.topology = newTopology

    def _initialize(self, simulation):
        """Deferred initialization of the reporter, which happens before
        processing the first report.

        At the time that the first report is processed, we now have access
        to the simulation object, which we don't have at the point when the
        reporter is instantiated

        Parameters
        ----------
        simulation : simtk.openmm.app.Simulation
            The Simulation to generate a report for
        """
        if self._atomSubset is not None:
            if not min(self._atomSubset) >= 0:
                raise ValueError('atomSubset must be zero indexed. '
                                 'the smallest allowable value is zero')
            if not max(self._atomSubset) < simulation.system.getNumParticles():
                raise ValueError('atomSubset must be zero indexed. '
                                 'the largest value must be less than the number '
                                 'of particles in the system')
            if not all(a==int(a) for a in self._atomSubset):
                raise ValueError('all of the indices in atomSubset must be integers')

            self._atomSlice = self._atomSubset
            if hasattr(self._traj_file, 'topology'):
                self._traj_file.topology = _topology_from_subset(
                    Topology.from_openmm(simulation.topology), self._atomSubset)
        else:
            self._atomSlice = slice(None)
            if hasattr(self._traj_file, 'topology'):
                self._traj_file.topology = Topology.from_openmm(simulation.topology)


        system = simulation.system
        if self._temperature:
            # Compute the number of degrees of freedom.
            dof = 0
            for i in range(system.getNumParticles()):
                if system.getParticleMass(i) > 0*units.dalton:
                    dof += 3
            dof -= system.getNumConstraints()
            if any(type(system.getForce(i)) == mm.CMMotionRemover for i in range(system.getNumForces())):
                dof -= 3
            self._dof = dof

        if simulation.topology.getUnitCellDimensions() is None:
            self._cell = False

def load_mol2(filename):
    """Load a TRIPOS mol2 file from disk.

    Parameters
    ----------
    filename : str
        Path to the prmtop file on disk.

    Returns
    -------
    traj : md.Trajectory
        The resulting topology, as an md.Topology object.

    Notes
    -----
    This function should work on GAFF and sybyl style MOL2 files, but has
    been primarily tested on GAFF mol2 files.
    This function does NOT accept multi-structure MOL2 files!!!
    The elements are guessed using GAFF atom types or via the atype string.

    Examples
    --------
    >>> traj = md.load_mol2('mysystem.mol2')
    """
    from mdtraj.core.trajectory import Trajectory
    from mdtraj.core.topology import Topology

    atoms, bonds = mol2_to_dataframes(filename)

    atoms_mdtraj = atoms[["name", "resName"]].copy()
    atoms_mdtraj["serial"] = atoms.index

    #Figure out 1 letter element names

    # IF this is a GAFF mol2, this line should work without issues
    atoms_mdtraj["element"] = atoms.atype.map(gaff_elements)
    # If this is a sybyl mol2, there should be NAN (null) values
    if atoms_mdtraj.element.isnull().any():
        # If this is a sybyl mol2, I think this works generally.
        atoms_mdtraj["element"] = atoms.atype.apply(lambda x: x.strip(".")[0])

    atoms_mdtraj["resSeq"] = np.ones(len(atoms), 'int')
    atoms_mdtraj["chainID"] = np.ones(len(atoms), 'int')

    if bonds is not None:
        bonds_mdtraj = bonds[["id0", "id1"]].values
        offset = bonds_mdtraj.min()  # Should this just be 1???
        bonds_mdtraj -= offset
    else:
        bonds_mdtraj = None

    top = Topology.from_dataframe(atoms_mdtraj, bonds_mdtraj)

    xyzlist = np.array([atoms[["x", "y", "z"]].values])
    xyzlist /= 10.0  # Convert from angstrom to nanometer

    traj = Trajectory(xyzlist, top)

    return traj

def create_water_topology_on_disc(n):
    topfile = tempfile.mktemp('.pdb')
    top = Topology()
    chain = top.add_chain()

    for i in xrange(n):
        res = top.add_residue('r%i' % i, chain)
        h1 = top.add_atom('H', hydrogen, res)
        o = top.add_atom('O', oxygen, res)
        h2 = top.add_atom('H', hydrogen, res)
        top.add_bond(h1, o)
        top.add_bond(h2, o)

    xyz = np.zeros((n * 3, 3))
    Trajectory(xyz, top).save_pdb(topfile)
    return topfile

    def mutate(self, mut_res_idx, mut_new_resname):
        """Mutate residue

        Parameters 
        ----------
        mut_res_idx : int
            Index of residue to mutate.
        mut_new_resname : str
            Three-letter code of residue to mutate to.
        """

        assert (self.topology.residue(mut_res_idx).name != mut_new_resname), "mutating the residue to itself!"

        # Build new topology
        newTopology = Topology()
        for chain in self.topology.chains:
            newChain = newTopology.add_chain()
            for residue in chain._residues:
                res_idx = residue.index
                if res_idx == mut_res_idx:
                    # create mutated residue
                    self._add_mutated_residue(mut_new_resname, newTopology, newChain, res_idx, residue)
                else:
                    # copy old residue atoms directly
                    newResidue = newTopology.add_residue(residue.name, newChain, res_idx)
                    for atom in residue.atoms:
                        newTopology.add_atom(atom.name, 
                                    md.core.element.get_by_symbol(atom.element.symbol), 
                                    newResidue, serial=atom.index)

        # The bond connectivity should stay identical
        for atm1, atm2 in self.topology._bonds:
            new_atm1 = newTopology.atom(atm1.index)
            new_atm2 = newTopology.atom(atm2.index)
            newTopology.add_bond(new_atm1, new_atm2)

        self._prev_topology = self.topology.copy()
        self.topology = newTopology

    def topology(self):
        """Get the topology out from the file

        Returns
        -------
        topology : mdtraj.Topology
            A topology object
        """
        try:
            raw = self._get_node('/', name='topology')[0]
            if not isinstance(raw, string_types):
                raw = raw.decode()
            topology_dict = json.loads(raw)
        except self.tables.NoSuchNodeError:
            return None

        topology = Topology()

        for chain_dict in sorted(topology_dict['chains'], key=operator.itemgetter('index')):
            chain = topology.add_chain()
            for residue_dict in sorted(chain_dict['residues'], key=operator.itemgetter('index')):
                try:
                    resSeq = residue_dict["resSeq"]
                except KeyError:
                    resSeq = None
                    warnings.warn('No resSeq information found in HDF file, defaulting to zero-based indices')
                try:
                    segment_id = residue_dict["segmentID"]
                except KeyError:
                    segment_id = ""
                residue = topology.add_residue(residue_dict['name'], chain, resSeq=resSeq, segment_id=segment_id)
                for atom_dict in sorted(residue_dict['atoms'], key=operator.itemgetter('index')):
                    try:
                        element = elem.get_by_symbol(atom_dict['element'])
                    except KeyError:
                        element = elem.virtual
                    topology.add_atom(atom_dict['name'], element, residue)

        atoms = list(topology.atoms)
        for index1, index2 in topology_dict['bonds']:
            topology.add_bond(atoms[index1], atoms[index2])

        return topology

    def topology(self):
        """Get the topology out from the file

        Returns
        -------
        topology : mdtraj.Topology
            A topology object
        """
        try:
            raw = self._get_node("/", name="topology")[0]
            if not isinstance(raw, string_types):
                raw = raw.decode()
            topology_dict = json.loads(raw)
        except self.tables.NoSuchNodeError:
            return None

        topology = Topology()

        for chain_dict in sorted(topology_dict["chains"], key=operator.itemgetter("index")):
            chain = topology.add_chain()
            for residue_dict in sorted(chain_dict["residues"], key=operator.itemgetter("index")):
                try:
                    resSeq = residue_dict["resSeq"]
                except KeyError:
                    resSeq = None
                    warnings.warn("No resSeq information found in HDF file, defaulting to zero-based indices")
                residue = topology.add_residue(residue_dict["name"], chain, resSeq=resSeq)
                for atom_dict in sorted(residue_dict["atoms"], key=operator.itemgetter("index")):
                    try:
                        element = elem.get_by_symbol(atom_dict["element"])
                    except KeyError:
                        element = None
                    topology.add_atom(atom_dict["name"], element, residue)

        atoms = list(topology.atoms)
        for index1, index2 in topology_dict["bonds"]:
            topology.add_bond(atoms[index1], atoms[index2])

        return topology

    def _read_models(self):
        if not self._mode == 'r':
            raise ValueError('file not opened for reading')

        self._topology = Topology()

        pdb = PdbStructure(self._file, load_all_models=True)

        atomByNumber = {}
        for chain in pdb.iter_chains():
            c = self._topology.add_chain()
            for residue in chain.iter_residues():
                resName = residue.get_name()
                if resName in PDBTrajectoryFile._residueNameReplacements:
                    resName = PDBTrajectoryFile._residueNameReplacements[resName]
                r = self._topology.add_residue(resName, c, residue.number)
                if resName in PDBTrajectoryFile._atomNameReplacements:
                    atomReplacements = PDBTrajectoryFile._atomNameReplacements[resName]
                else:
                    atomReplacements = {}
                for atom in residue.atoms:
                    atomName = atom.get_name()
                    if atomName in atomReplacements:
                        atomName = atomReplacements[atomName]
                    atomName = atomName.strip()
                    element = atom.element
                    if element is None:
                        element = self._guess_element(atomName, residue)

                    newAtom = self._topology.add_atom(atomName, element, r, serial=atom.serial_number)
                    atomByNumber[atom.serial_number] = newAtom

        # load all of the positions (from every model)
        _positions = []
        for model in pdb.iter_models(use_all_models=True):
            coords = []
            for chain in model.iter_chains():
                for residue in chain.iter_residues():
                    for atom in residue.atoms:
                        coords.append(atom.get_position())
            _positions.append(coords)

        if not all(len(f) == len(_positions[0]) for f in _positions):
            raise ValueError('PDB Error: All MODELs must contain the same number of ATOMs')

        self._positions = np.array(_positions)

        ## The atom positions read from the PDB file
        self._unitcell_lengths = pdb.get_unit_cell_lengths()
        self._unitcell_angles = pdb.get_unit_cell_angles()
        self._topology.create_standard_bonds()
        self._topology.create_disulfide_bonds(self.positions[0])

        # Add bonds based on CONECT records.
        connectBonds = []
        for connect in pdb.models[0].connects:
            i = connect[0]
            for j in connect[1:]:
                if i in atomByNumber and j in atomByNumber:
                    connectBonds.append((atomByNumber[i], atomByNumber[j]))
        if len(connectBonds) > 0:
            # Only add bonds that don't already exist.
            existingBonds = set(self._topology.bonds)
            for bond in connectBonds:
                if bond not in existingBonds and (bond[1], bond[0]) not in existingBonds:
                    self._topology.add_bond(bond[0], bond[1])
                    existingBonds.add(bond)

class PDBTrajectoryFile(object):
    """Interface for reading and writing Protein Data Bank (PDB) files

    Parameters
    ----------
    filename : str
        The filename to open. A path to a file on disk.
    mode : {'r', 'w'}
        The mode in which to open the file, either 'r' for read or 'w' for write.
    force_overwrite : bool
        If opened in write mode, and a file by the name of `filename` already
        exists on disk, should we overwrite it?

    Attributes
    ----------
    positions : np.ndarray, shape=(n_frames, n_atoms, 3)
    topology : mdtraj.Topology
    closed : bool

    Notes
    -----
    When writing pdb files, mdtraj follows the PDB3.0 standard as closely as
    possible. During *reading* however, we try to be more lenient. For instance,
    we will parse common nonstandard atom names during reading, and convert them
    into the standard names. The replacement table used by mdtraj is at
    {mdtraj_source}/formats/pdb/data/pdbNames.xml.

    See Also
    --------
    mdtraj.load_pdb : High-level wrapper that returns a ``md.Trajectory``
    """
    distance_unit = 'angstroms'
    _residueNameReplacements = {}
    _atomNameReplacements = {}
    _chain_names = [chr(ord('A') + i) for i in range(26)]

    def __init__(self, filename, mode='r', force_overwrite=True):
        self._open = False
        self._file = None
        self._topology = None
        self._positions = None
        self._mode = mode
        self._last_topology = None

        if mode == 'r':
            PDBTrajectoryFile._loadNameReplacementTables()

            if _is_url(filename):
                self._file = urlopen(filename)
                if filename.lower().endswith('.gz'):
                    if six.PY3:
                        self._file = gzip.GzipFile(fileobj=self._file)
                    else:
                        self._file = gzip.GzipFile(fileobj=six.StringIO(
                            self._file.read()))
                if six.PY3:
                    self._file = six.StringIO(self._file.read().decode('utf-8'))
            else:
                if filename.lower().endswith('.gz'):
                    self._file = gzip.open(filename, 'r')
                    self._file = six.StringIO(self._file.read().decode('utf-8'))                    
                else:
                    self._file = open(filename, 'r')

            self._read_models()
        elif mode == 'w':
            self._header_written = False
            self._footer_written = False
            if os.path.exists(filename) and not force_overwrite:
                raise IOError('"%s" already exists' % filename)
            self._file = open(filename, 'w')
        else:
            raise ValueError("invalid mode: %s" % mode)

        self._open = True

    def write(self, positions, topology, modelIndex=None, unitcell_lengths=None, 
              unitcell_angles=None, bfactors=None):
        """Write a PDB file to disk

        Parameters
        ----------
        positions : array_like
            The list of atomic positions to write.
        topology : mdtraj.Topology
            The Topology defining the model to write.
        modelIndex : {int, None}
            If not None, the model will be surrounded by MODEL/ENDMDL records
            with this index
        unitcell_lengths : {tuple, None}
            Lengths of the three unit cell vectors, or None for a non-periodic system
        unitcell_angles : {tuple, None}
            Angles between the three unit cell vectors, or None for a non-periodic system
        bfactors : array_like, default=None, shape=(n_atoms,)
            Save bfactors with pdb file. Should contain a single number for
            each atom in the topology
        """
        if not self._mode == 'w':
            raise ValueError('file not opened for writing')
        if not self._header_written:
#.........这里部分代码省略.........

    def _to_topology(self, atom_list, chain_types=None, residue_types=None):
        """Create a mdtraj.Topology from a Compound.

        Parameters
        ----------
        atom_list :
        chain_types :
        residue_types :

        Returns
        -------
        top : mtraj.Topology

        """

        if isinstance(chain_types, Compound):
            chain_types = [Compound]
        if isinstance(chain_types, (list, set)):
            chain_types = tuple(chain_types)

        if isinstance(residue_types, Compound):
            residue_types = [Compound]
        if isinstance(residue_types, (list, set)):
            residue_types = tuple(residue_types)
        top = Topology()
        atom_mapping = {}

        default_chain = top.add_chain()
        default_residue = top.add_residue('RES', default_chain)

        last_residue_compound = None
        last_chain_compound = None
        last_residue = None
        last_chain = None

        for atom in atom_list:
            # Chains
            for parent in atom.ancestors():
                if chain_types and isinstance(parent, chain_types):
                    if parent != last_chain_compound:
                        last_chain_compound = parent
                        last_chain = top.add_chain()
                        last_chain_default_residue = top.add_residue('RES', last_chain)
                        last_chain.compound = last_chain_compound
                    break
            else:
                last_chain = default_chain
                last_chain.compound = last_chain_compound

            # Residues
            for parent in atom.ancestors():
                if residue_types and isinstance(parent, residue_types):
                    if parent != last_residue_compound:
                        last_residue_compound = parent
                        last_residue = top.add_residue(parent.__class__.__name__, last_chain)
                        last_residue.compound = last_residue_compound
                    break
            else:
                if last_chain != default_chain:
                    last_residue = last_chain_default_residue
                else:
                    last_residue = default_residue
                last_residue.compound = last_residue_compound

            # Add the actual atoms
            try:
                elem = get_by_symbol(atom.name)
            except KeyError:
                elem = get_by_symbol("VS")
            at = top.add_atom(atom.name, elem, last_residue)
            at.charge = atom.charge
            atom_mapping[atom] = at

        # Remove empty default residues.
        chains_to_remove = [chain for chain in top.chains if chain.n_atoms == 0]
        residues_to_remove = [res for res in top.residues if res.n_atoms == 0]
        for chain in chains_to_remove:
            top._chains.remove(chain)
        for res in residues_to_remove:
            for chain in top.chains:
                try:
                    chain._residues.remove(res)
                except ValueError:  # Already gone.
                    pass

        for atom1, atom2 in self.bonds():
            # Ensure that both atoms are part of the compound. This becomes an
            # issue if you try to convert a sub-compound to a topology which is
            # bonded to a different subcompound.
            if all(a in atom_mapping.keys() for a in [atom1, atom2]):
                top.add_bond(atom_mapping[atom1], atom_mapping[atom2])
        return top

def load_hoomdxml(filename, top=None):
    """Load a single conformation from an HOOMD-Blue XML file.

    For more information on this file format, see:
    http://codeblue.umich.edu/hoomd-blue/doc/page_xml_file_format.html
    Notably, all node names and attributes are in all lower case.
    HOOMD-Blue does not contain residue and chain information explicitly. 
    For this reason, chains will be found by looping over all the bonds and 
    finding what is bonded to what. 
    Each chain consisists of exactly one residue. 

    Parameters
    ----------
    filename : string
        The path on disk to the XML file

    Returns
    -------
    trajectory : md.Trajectory
        The resulting trajectory, as an md.Trajectory object, with corresponding 
        Topology.

    Notes
    -----
    This function requires the NetworkX python package.
    """
    from mdtraj.core.trajectory import Trajectory
    from mdtraj.core.topology import Topology
    topology = Topology()
    tree = cElementTree.parse(filename)
    config = tree.getroot().find('configuration')
    position = config.find('position')
    bond = config.find('bond')
    atom_type = config.find('type')  # MDTraj calls this "name"

    box = config.find('box')
    box.attrib = dict((key.lower(), val) for key, val in box.attrib.items())
    # be generous for case of box attributes
    lx = float(box.attrib['lx'])
    ly = float(box.attrib['ly'])
    lz = float(box.attrib['lz'])
    try:
        xy = float(box.attrib['xy'])
        xz = float(box.attrib['xz'])
        yz = float(box.attrib['yz'])
    except:
        xy = 0.0
        xz = 0.0
        yz = 0.0
    unitcell_vectors = np.array([[[lx,  xy*ly, xz*lz],
                                  [0.0, ly,    yz*lz],
                                  [0.0, 0.0,   lz   ]]])

    positions, types = [], {}
    for pos in position.text.splitlines()[1:]:
        positions.append((float(pos.split()[0]),
                          float(pos.split()[1]),
                          float(pos.split()[2])))

    for idx, atom_name in enumerate(atom_type.text.splitlines()[1:]):
        types[idx] = str(atom_name.split()[0])
    if len(types) != len(positions):
        raise ValueError('Different number of types and positions in xml file')

    # ignore the bond type
    bonds = [(int(b.split()[1]), int(b.split()[2])) for b in bond.text.splitlines()[1:]]
    chains = _find_chains(bonds)
    ions = [i for i in range(len(types)) if not _in_chain(chains, i)]

    # add chains, bonds and ions (each chain = 1 residue)
    for chain in chains:
        t_chain = topology.add_chain()
        t_residue = topology.add_residue('A', t_chain)
        for atom in chain:
            topology.add_atom(types[atom], 'U', t_residue)
    for ion in ions:
        t_chain = topology.add_chain()
        t_residue = topology.add_residue('A', t_chain)
        topology.add_atom(types[atom], 'U', t_residue)
    for bond in bonds:
        atom1, atom2 = bond[0], bond[1]
        topology.add_bond(topology.atom(atom1), topology.atom(atom2))

    traj = Trajectory(xyz=np.array(positions), topology=topology)
    traj.unitcell_vectors = unitcell_vectors

    return traj

def load_mol2(filename):
    """Load a TRIPOS mol2 file from disk.

    Parameters
    ----------
    filename : str
        Path to the prmtop file on disk.

    Returns
    -------
    traj : md.Trajectory
        The resulting topology, as an md.Topology object.

    Notes
    -----
    This function should work on GAFF and sybyl style MOL2 files, but has
    been primarily tested on GAFF mol2 files.
    This function does NOT accept multi-structure MOL2 files!!!
    The elements are guessed using GAFF atom types or via the atype string.

    Examples
    --------
    >>> traj = md.load_mol2('mysystem.mol2')
    """
    from mdtraj.core.trajectory import Trajectory
    from mdtraj.core.topology import Topology, Single, Double, Triple, Aromatic, Amide

    atoms, bonds = mol2_to_dataframes(filename)

    atoms_mdtraj = atoms[["name", "resName"]].copy()
    atoms_mdtraj["serial"] = atoms.index

    #Figure out 1 letter element names

    # IF this is a GAFF mol2, this line should work without issues
    atoms_mdtraj["element"] = atoms.atype.map(gaff_elements)
    # If this is a sybyl mol2, there should be NAN (null) values
    if atoms_mdtraj.element.isnull().any():
        # If this is a sybyl mol2, I think this works generally.
        # Argument x is being passed as a list with only one element.
        def to_element(x):
            if isinstance(x, (list, tuple)):
                assert len(x) == 1
                x = x[0]

            if '.' in x:  # orbital-hybridizations in SYBL
                return x.split('.')[0]
            try:
                # check if we can convert the whole str to an Element,
                # if not, we only pass the first letter.
                from mdtraj.core.element import Element
                Element.getBySymbol(x)
            except KeyError:
                return to_element(x[0])
            return x
        atoms_mdtraj["element"] = atoms.atype.apply(to_element)

    atoms_mdtraj["resSeq"] = np.ones(len(atoms), 'int')
    atoms_mdtraj["chainID"] = np.ones(len(atoms), 'int')

    bond_type_map = {
        '1': Single,
        '2': Double,
        '3': Triple,
        'am': Amide,
        'ar': Aromatic
    }
    if bonds is not None:
        bonds_mdtraj = bonds[["id0", "id1"]].values
        offset = bonds_mdtraj.min()  # Should this just be 1???
        bonds_mdtraj -= offset
        # Create the bond augment information
        n_bonds = bonds_mdtraj.shape[0]
        bond_augment = np.zeros([n_bonds, 2], dtype=float)
        # Add bond type information
        bond_augment[:, 0] = [float(bond_type_map[str(bond_value)]) for bond_value in bonds["bond_type"].values]
        # Add Bond "order" information, this is not known from Mol2 files
        bond_augment[:, 1] = [0.0 for _ in range(n_bonds)]
        # Augment array, dtype is cast to minimal representation of float
        bonds_mdtraj = np.append(bonds_mdtraj, bond_augment, axis=-1)
    else:
        bonds_mdtraj = None

    top = Topology.from_dataframe(atoms_mdtraj, bonds_mdtraj)

    xyzlist = np.array([atoms[["x", "y", "z"]].values])
    xyzlist /= 10.0  # Convert from angstrom to nanometer

    traj = Trajectory(xyzlist, top)

    return traj

    def __init__(self, topology, use_chains=None):
        if use_chains is None:
            use_chains = range(len(topology._chains))

        self._ref_topology = topology.copy()

        # Build new topology
        newTopology = Topology()
        new_atm_idx = 0
        res_idx = 1
        prev_ca = None
        ca_idxs = []
        self._sidechain_idxs = []
        self._sidechain_mass = []
        self._chain_indices = []
        for chain_count, chain in enumerate(topology._chains):
            if chain_count in use_chains:
                newChain = newTopology.add_chain()
                for residue in chain._residues:
                    #resSeq = getattr(residue, 'resSeq', None) or residue.index
                    newResidue = newTopology.add_residue(residue.name, newChain, res_idx)
                    # map CA
                    new_ca = newTopology.add_atom('CA', md.core.element.get_by_symbol('C'),
                                        newResidue, serial=new_atm_idx)
                    self._chain_indices.append(chain_count)
                    if prev_ca is None:
                        prev_ca = new_ca
                    else:
                        # only bond atoms in the same chain.
                        if new_ca.residue.chain.index == prev_ca.residue.chain.index:
                            newTopology.add_bond(prev_ca, new_ca)
                        prev_ca = new_ca
                    try:
                        ca_idxs.append([[ atm.index for atm in residue.atoms if \
                                (atm.name == "CA") ][0], new_atm_idx ])
                    except:
                        print(residue)
                        print(chain)
                        for atm in residue.atoms:
                            atm.name
                        raise
                    new_atm_idx += 1

                    if residue.name == 'GLY':
                        self._sidechain_idxs.append([])
                        self._sidechain_mass.append([])
                    else:
                        # map CB
                        cb_name = "CB%s" % atom_types.residue_code[residue.name]
                        new_cb = newTopology.add_atom(cb_name, md.core.element.get_by_symbol('C'),
                                            newResidue, serial=new_atm_idx)
                        self._chain_indices.append(chain_count)

                        newTopology.add_bond(new_cb, new_ca)

                        self._sidechain_idxs.append([[ atm.index for atm in residue.atoms if \
                                    (atm.is_sidechain) and (atm.element.symbol != "H") ], new_atm_idx ])
                        self._sidechain_mass.append(np.array([ atm.element.mass for atm in residue.atoms if \
                                    (atm.is_sidechain) and (atm.element.symbol != "H") ]))
                        new_atm_idx += 1
                    res_idx += 1

        self._ca_idxs = np.array(ca_idxs)
        self.topology = newTopology
        assert self.topology.n_atoms == len(self._chain_indices)