Hi,
some ideas.
is.dag() (in R) tells you whether there are loops at all. But probably you know that there are some.
Start graph.dfs() from a vertex with in-degree 0, and then look at the DFS order of vertices, and look for edges that are pointing backwards according to the DFS order. I believe that these must signal loops.
Alternatively, calculate the strongly connected components of the graphs. All vertices of a loop belong to the same strongly connected component. Also, each non-trivial strongly connected component must contain at least one directed loop. This is also a way to parallelize the problem, because you can just work with the individual strongly connected components.
I cannot say for sure that these methods work for your graph, so try on a smaller (snowball) sample first, and see how your solution scales with the number of edges and vertices.
Best,
Gabor