You are asking the right question. And you can use kmeans!!! Despite what you may be told by some, you absolutely can cluster with kmeans. There is nothing about binary data that will cause kmeans to fail. However, you might want to consider the following:
1 - Zero-mean your matrix by column. This means that you compute the mean row vector, which now becomes a real valued vector, and then subtract that vector from each of the original binary vectors. Your 0/1 binary matrix of 650K row vectors now becomes a real valued matrix of 650K vectors. Note that this DOES NOT change the mutual distance (or similarity) between vectors. It is just a translation operation, applied identically to each vector.
2 - Apply the sign function to the matrix. The sign function forces each matrix element to -1 if it is negative, or to +1 otherwise. The result of this transformation, in steps 1 and 2, is that the new matrix is no longer sparse.
3 - Now apply kmeans. you can use the Euclidean metric, or experiment with other metrics that you kmeans implementation supports. No need to use a specific binary clustering algorithm. kmeans is simple and clustering 650K vectors should be easily feasible on a decent desktop.
4 - If you wish to have binary cluster vectors as the result, then apply the sign function to the final k clusters. You may also convert the final cluster vectors from +1/-1 representation to 0/1 representation (but only after applying the sign function).
Things to note:
Because you only have 62 dimensional vectors the range of 'similarity' values that are possible between vectors in the binary representation is 62 (corresponding to a Hamming distance between 0 and 62.) Since the range of distances between binary vectors is thus limited, any ranking by hamming distances will necessarily result in numerous ties. As you try to squeeze 650K vectors into only 62 possible distance buckets, the number of vectors per bucket will depend on the number of clusters, but will generally be large and you may need to resolve ties by going back to the original data from which you derived the initial binary matrix.