Perhaps there's some smart way to count the number of permutations of the 52 cards such that each player gets at least one card of each suite, which is the complement of what is asked for. But sometimes we do not need the exact solution and a simulation can give a pretty good approximation. In this case, with $10^6$ samples, we can estimate the probability being asked for is about $0.184171$. I bet this is pretty close the exact result.
Here's the code in Mathematica
In[2]:= n = 52
Out[2]= 52
In[13]:= n0 = n/4
Out[13]= 13
In[3]:= l = Range[n]
Out[3]= {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, \
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, \
35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52}
In[14]:= ll = Partition[l, n0]
Out[14]= {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13}, {14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25, 26}, {27, 28, 29, 30, 31, 32,
33, 34, 35, 36, 37, 38, 39}, {40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51, 52}}
In[39]:= containAll[sample_] := Module[{i, sl, crossl},
sl = Partition[sample, n0];
crossl = Tuples[{sl, ll}];
If[And @@ Map[ContainsAny[#[[1]], #[[2]]] &, crossl], 1, 0]
]
In[53]:= s0 =
Total[Parallelize[
Map[(Permute[l, #] // containAll) &,
RandomPermutation[n, 1000000]]]] // AbsoluteTiming
Out[53]= {191.22, 815829}
In[59]:= 1 - 815829/1000000 // N
Out[59]= 0.184171