This code computes just the numerical data for A321005. 

CountPairs := proc (p::prime) local b, a, n; option remember; a := 2; n := 0; while a < p do b := `mod`(1/a, p); if a <= b and isprime(b) then n := n+1 end if; a := nextprime(a) end do; n end proc; seq(CountPairs(ithprime(k)), k = 1 .. 1000)


This is the code for generating a graph (It takes a little while to run).

CountPairs:= proc(p::prime)
option remember;
local b, a:= 2, n:= 0;
   while a < p do
      b:= 1/a mod p; 
      if a <= b and isprime(b) then n:= n+1 fi;
      a:= nextprime(a)
   od;
   n
end proc:

#Compute as many terms as we can within a certain time limit:
try
   timelimit(
      999, #seconds
      #This infinite loop is intentional:
      proc() local p:= 1; do p:= nextprime(p); CountPairs(p) od end proc() 
   )
catch:
end try:

#Extract remember table and convert to a Matrix:
PP:= Matrix(
   [lhs,rhs]~([entries(op(4, eval(CountPairs)), pairs, indexorder)]),
    datatype= float[8]
):
N:= max(PP):
Asy:= p-> Li(p)^2/2/p: #1st asymptotic term
Res:= p-> 2*sqrt(p)/Pi^2: #estimated residual bound
plot(
   [PP, Asy(n), Asy(n)+Res(n), Asy(n)-Res(n)], n= 2..N, 
   style= [point,line$3], color= [grey,red,blue$2], thickness= [0,1,0,0], 
   symbolsize= 1, symbol= point
);