def tile_bonds(bonds, reps, num_solvent_atoms):
    """
    Like the Numpy tile function, tile the bonds array.
    The Numpy tile function cannot be used because the 
    values in the bonds array need to be incremented by
    the number of atoms in the solvent molecules.

    bonds : Psi-OMM Bonds array
        The bonds array for the solvent molecule.
    reps : int
        Number of repeats of the bonds array. For a simple
        array, arr=[1,2], tile(arr, 3) would result in
        [1,2,1,2,1,2].
    num_solvent_atoms : int
        Number of atoms in the solvent molecule. This value
        is used to increment every atom index in the bonds array
        (this is every value except bond orders).

    Returns
    -------
    List of lists of lists in the form of the bonds array.
    """
    ret_bonds0 = bonds[0].copy()
    ret_bonds1 = bonds[1].copy()

    for tile_ix in range(1, reps):
        for b0 in bonds[0]:
            working_bond = b0.copy()
            # Increment the indices at working_bond[0] and [1] by
            # the num_solvent_atoms. working_bond[2] is the bond
            # order and should be left alone.
            working_bond[0] += num_solvent_atoms * tile_ix
            working_bond[1] += num_solvent_atoms * tile_ix

            ret_bonds0.append(working_bond)

        for b1 in bonds[1]:
            working_bond = b1.copy()
            # Increment every value in the list by num_solvent_atoms.
            for ix in range(len(working_bond)):
                working_bond[ix] += num_solvent_atoms * tile_ix

            ret_bonds1.append(working_bond)

    # Return the new tiled bond array in the form of a typical bonds
    # array, (bonds0, bonds1)
    return (ret_bonds0, ret_bonds1)