def bakan_bahar_ensemble_align(coords, tolerance = 0.001, verbose = False ):
    '''
    input: a list of coordinates in the format:
    [
      [ (x, y, z), (x, y, z), (x, y, z) ], # atoms in model 1
      [ (x, y, z), (x, y, z), (x, y, z) ], # atoms in model 2
      # etc.
    ]
    '''
    rmsd_tolerance = float("inf")
    average_struct_coords = np.array( random.choice(coords) )
    cycle_count = 0
    while( rmsd_tolerance > tolerance ):
        if verbose:
            print 'Cycle %d alignment, current tolerance: %.4f (threshold %.4f)' % (cycle_count, rmsd_tolerance, tolerance)
        old_average_struct_coords = average_struct_coords
        coords = np.array([ kabsch( average_struct_coords, x ) for x in coords ])
        average_struct_coords = np.mean( coords, axis = 0 )
        rmsd_tolerance = rmsd( old_average_struct_coords, average_struct_coords )
        cycle_count += 1
    return coords