def generate_bond_indices(natoms):
    """
    natoms: int
        The number of atoms
    Finds the array of bond indices of the upper triangle of an interatomic distance matrix, in column wise order
    ( or equivalently, lower triangle of interatomic distance matrix in row wise order):
    [[0,1], [0,2], [1,2], [0,3], [1,3], [2,3], ...,[0, natom], ...,[natom-1, natom]]
    """ 
    # 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