Score 100

{3, 7, 19, 23, 31, 37, 43, 47, 59, 67, 79, 101, 131, 137, 139, 191, 211, 239, 283, 431, 443, 463, 593, 599, 647, 751, 829, 859, 1039, 1103, 1129, 1237, 1597, 1871, 2029, 2251, 2293, 2311, 2341, 2347, 2377, 2521, 2531, 2591, 2731, 2927, 3061, 3121, 3169, 3221, 3361, 3457, 3571, 3911, 4003, 4523, 4831, 5279, 6091, 6271, 6553, 7393, 7481, 9001, 10501, 10847, 11551, 12421, 12973, 14281, 15541, 15661, 18061, 18481, 20749, 21121, 21313, 22541, 34217, 42331, 45289, 45361, 49201, 52541, 54721, 60601, 61813, 65701, 87721, 88609, 131581, 134677, 156157, 164341, 178537, 217561, 253381, 254041, 430921, 435481}

With a little more time I suspect I could get a much larger set.

I started by finding (using Mathematica) all 148 odd primes p such that p–1 is a product of distinct primes not exceeding 29. (There are 10 primes up to 29, and it's easy to check all 2^10 possibilities.) From this set I choose larger and larger subsets, eventually stopping at a subset of size 53, and fed them into an Egyptian fraction solver that a research student of mine, Yue Shi, coded up; the simple prime factorizations of the inputs made finding solutions much more likely. The program yielded a solution of size 44, which I used as a seed to the next step. (In hindsight, the next step worked so well that I suspect starting with a much simpler solution would have worked just fine.)

Next, I examined each prime \$p\$ in the current solution one by one, looking for ways to write\$\frac1 {p - 1} = \frac1 {q_ 1 - 1} + \frac1 {q_ 2 - 1}\$ with \$q_ 1\$ and \$q_ 2\$ both prime. It is known that the integer solutions to \$\frac1n = \frac1x + \ \frac1y\$ are in one-to-one correspondence with factorizations \$n^2 = f_ 1 f_ 2\$, the correspondence being given by \$x = n + f_ 1\$, \$y = n + f_ 2\$. I searched through the factorizations of \$ (p - 1)^2\$ to see if any of them yielded a solution where both new denominators were one less then primes; if so, I replaced \$p\$ by \$q_1, q_ 2\$ in the solution.

This process can be iterated, not just on each prime in the original solution but also on the primes used to replace old primes. I observed about 2 such splitting-iterations on average for every prime in the original solution, which leads me to suspect that the process would run quite a bit longer before terminating; I was doing the splitting by hand, and so I just stopped when I reached a solution of size 100.

Score 100 8605

I used an algorithm that starts with one solution and repeatedly tries to split a prime \$p\$ in the solution into two other primes \$q_ 1\$ and \$q_ 2\$ that satisfy \$\frac1{p-1} = \frac1{q_1-1}+\frac1{q_2-1}\$.

It is known (and can be quickly checked) that the positive integer solutions to \$\frac1n = \frac1x + \frac1y\$ are in one-to-one correspondence with factorizations \$n^2 = f_ 1 f_ 2\$, the correspondence being given by \$x = n + f_ 1\$, \$y = n + f_ 2\$. We can search through the factorizations of \$(p-1)^2\$ to see if any of them yielded a solution where both new denominators \$x,y\$ were one less then primes; if so, then \$p\$ can be replaced by \$q_1=x+1\$, \$q_2=y+1\$. (If there were multiple such factorizations, I used the one where the minimum of \$q_1\$ and \$q_2\$ was smallest.)

I started with this seed solution of length 44:

seed = {3, 7, 11, 23, 31, 43, 47, 67, 71, 79, 103, 131, 139, 191, 211, 239, 331, 419, 443, 463, 547, 571, 599, 647, 691, 859, 911, 967, 1103, 1327, 1483, 1871, 2003, 2311, 2347, 2731, 3191, 3307, 3911, 4003, 4931, 6007, 6091, 8779}

This seed was found using an Egyptian fraction solver that a former research student of mine, Yue Shi, coded in ChezScheme. (The primes \$p\$ involved all have the property that \$p-1\$ is the product of distinct primes less than 30, which increased the likelihood of Shi's program finding a solution.)

The following Mathematica code continually updates a current solution by looking at its primes one by one, trying to split them in the manner described above. (The number 1000 in the third line is an arbitrary stopping point; in principle one could let the algorithm run forever.)

solution = seed;
j = 1; (* j is the index of the element of the solution that we'll try to split *)

While[j <= 1000 && j <= Length[solution], 
  currentP = solution[[j]];
  allDivisors = Divisors[(currentP - 1)^2];
  allFactorizations = {#, (currentP - 1)^2/#} & /@
    Take[allDivisors, Floor[Length[allDivisors]/2]];
  allSplits = currentP + allFactorizations;
  goodFactorizations = Select[allSplits, 
    And @@ PrimeQ[#] && Intersection[#, solution] == {} &];
  If[goodFactorizations == {}, 
    j++, 
    solution = Union[Complement[solution, {currentP}], First@goodFactorizations]
  ]
]

The code above yields a solution of length 4126, whose largest element is about \$8.7\times10^{20}\$; by the end, it was factoring integers \$(p-1)^2\$ of size about \$8.8\times10^{21}\$.

In practice, I ran the code several times, using the previous output as the next seed in each case and increasing the cutoff for j each time; this allowed for the recovery of some small prime splits that had become non-redundant thanks to previous splitting, which somewhat mitigated the size of the integers the algorithm factored.

The final solution, which took about an hour to obtain, is too long to fit in this answer but has been posted online. It has length 8605 and largest element about \$4.62\times10^{19}\$.

Various runs of this code consistently found that the length of the solution was about 3–4 times as long as the set of primes that had been examined for splitting. In other words, the solution was growing much faster than the code scanned through the initial elements. It seems likely that this behavior would continue for a long time, yielding some gargantuan solutions.