I believe the probability is at least $\approx0.343$.

Let $\mu_n$ be a probability measure giving
$q_n = P_n[ X_1 + X_2 + X_3 < 2X_4]$.

Consider now $(Y_i)_{i\le 4}$ Bernoulli$(p)$. The $Y_i$'s produce the desired inequality __with an inclusive $\le$__ with probability
\begin{align}
I+II+III &= P[ \text{binomial}(3, p) \le 1] P[Y_4=1] 
\\&\quad+ P[ \text{binomial}(3, p) = 2] P[Y_4=1] 
\\&\quad+ (1-p)^4.
\end{align}
But the problem requires a strict inequality, for which the second and third terms are lost. We can boost this binomial scheme by using $X_1,X_2,X_3,X_4$ from $\mu_n$, independently $Y_i\sim$Bernoulli$(p)$ as above, and set
$$Z_i = t X_i + Y_i.$$
for some $t>0$ very small. Assuming that the support of $\mu_n$ is bounded, we can always choose $t>0$ small enough so that the first term is unchanged, while the second term and third terms become strict inequality with independent probability $q_n$:
$$
q_{n+1} = P[Z_1+Z_2+Z_3 < 2Z_4]
=
P[ \text{binomial}(3, p) \le 1] P[Y_4=1] 
+ q_n\Big(
P[ \text{binomial}(3, p) = 2] P[Y_4=1] 
+ (1-p)^4\Big).
$$
This defines two polynomial $f(p)$ and $g(p)$ of order $4$ in $p\in[0,1]$ such that
$$
q_{n+1} = f(p) + q_n g(p).
$$
Consider the fixed-point $q_\infty$ defined by
$$
q_{\infty} = 
f(p)
+ q_{\infty}
g(p).
$$
We can reach a large $q_\infty$ by maximizing
$$
M = \max_{p\in[0,1]} \frac{f(p)}{(1-g(p))_+}.
$$
If I am not mistaken, $f(p)= ((1-p)^3 + 3(1-p)^2p )p$
while $g(p) = 3p^2(1-p)p + (1-p)^4$.

The maximum is reached at $p^* \approx 0.404$ giving $M\approx 0.343$.

We now set $p=p^*$ for this maximizer, and perform this procedure, starting from some discrete distribution $q$ (e.g., Bernoulli($p^*$)). Since $g(p^*)\in(0,1)$, the sequence $q_n$ converges to $M\approx 0.343$.

Edit: using Will's explicit formula with $t=0.1$, with $n=30$ terms  confirm a probability of around $0.343$:
```
🖧[14]: import numpy as np
  ...: 
  ...: p = 0.404
  ...: t = 0.1
  ...: powers = np.arange(30)
  ...: ts = (t ** powers)[None, None, :]
  ...: 
  ...: bs = np.random.binomial(1, p, size=(100000, 4, 30))
  ...: ys = np.sum(ts * bs, axis=-1)
  ...: (np.sum(ys[:, :-1], axis=1) <  2 * ys[:, -1]).mean()
🖧[14]: 0.34353
```