def generate_bond_indices(natoms):
    """
    Finds the array of bond indices of an interatomic distance matrix, in row wise order:
    [[0,1], [0,2], [1,2], [0,3], [1,3], [2,3], ..., [0, natoms], [1, natoms], ...,[natoms-1, natoms]]
    
    Parameters
    ----------
    natoms: int
        The number of atoms

    Returns
    ----------
    bond_indices : list
        A list of lists, where each sublist is the subscripts of an interatomic distance
        from an interatomic distance matrix representation of a molecular system.
        e.g. r_12, r_01, r_05 
    """
    # initialize j as the number of atoms
    j = natoms - 1
    # now loop backward until you generate all bond indices 
    bond_indices = []
    while j > 0:
        i = j - 1
        while i >= 0:
            new = [i, j]
            bond_indices.insert(0, new)
            i -= 1
        j -= 1
    return bond_indices