Entering a cutoff radius and max_nieghbors does not produce the expected number of edges.

With a simple cubic structure (lattice parameter = 1 Angstrom) and setting cutoff=1.001 and max_nieghbors=12, nearest_neighbor_edges() produces 18 edges - which is both more edges than are within the specified cutoff and more than specified max_nieghbors.

Under these conditions the expected result should be 6 edges.

Cause: the if min_nbrs < max_neighbors: branch is activated since there are less than 12 neighbors within the radius.
This results in the parameter max_nieghbors actually becoming a minimum on the number of neighbors.

Potential solution: include an additional extend_radius flag as a parameter of nearest_neighbor_edges() that when True would activate this branch and when False would constrain the graph to the specified cutoff radius.

Example

In [0]: atoms
Out [0]: Cu
1.0
1 0 0
0 1 0
0 0 1
Cu
1
Cartesian
0.0 0.0 0.0

In [1]:  atoms.get_all_neighbors(r = 1.001)
Out [1]:  array([[[0, 0, 1.0, (-1.0, 0.0, 0.0)],
        [0, 0, 1.0, (0.0, -1.0, 0.0)],
        [0, 0, 1.0, (0.0, 0.0, -1.0)],
        [0, 0, 1.0, (0.0, 0.0, 1.0)],
        [0, 0, 1.0, (0.0, 1.0, 0.0)],
        [0, 0, 1.0, (1.0, 0.0, 0.0)]]], dtype=object)

In [2]: edges = nearest_neighbor_edges(atoms, cutoff=1.001, use_canonize=True)
    len(edges[(0,0)])
Out [2]: 18 

Should be 6 edges.