If we're lucky, the $f_j$, $1\le j\le n$ (I use index notation for partial derivatives) are already linearly independent. If not, then $v\cdot\nabla f=0$ for some vector $v\not= 0$, but then also $Pv\cdot\nabla f=0$ for all permutations of the components of $v$. I believe there are three possibilities, though I haven't proved this carefully: the $Pv$ will typically span $\mathbb R^n$ (then $f=0$), or they can span an $(n-1)$ dimensional subspace, so that $f$ can then only depend on $w\cdot x$; this corresponds to the exceptional example you mentioned. The only way to avoid such a large span is to have specifically $v=(1,1,\ldots , 1)$.
In this case, we can write
$$
f(x_1,x_2,\ldots ,x_n)=f(0,x_2-x_1, \ldots, x_n-x_1) = g(x_2-x_1,\ldots, x_n-x_1) ,
$$
and $g(t_2,\ldots , t_n)=f(0,t_2,\ldots ,t_n)$ is still a symmetric homogeneous polynomial of its arguments.
Now suppose we have an additional linear relation $\sum a_jf_j=\sum x_j f_{jj}$. We can rewrite this as
$$
\sum_{j\ge 2} b_j g_j(t_2,\ldots ,t_n) - \sum_{j\ge 2} t_j g_{jj}(t_2,\ldots ,t_n) = x_1\Delta f(0,t_2,\ldots ,t_n) ,
$$
with $t_j=x_j-x_1$ and $b_j=a_j-a_1$. Since $x_1$ and the $t_j$'s can be varied independently, both sides must be (identically) equal to zero. Now $g$ is symmetric, so we can as above manipulate for example the first term as follows:
$$
\left( \sum b_j D_j\right) g(t_2,\ldots, t_n) = \left( \sum b_j D_j\right) g(t_{\pi 2},\ldots, t_{\pi n}) = \sum b_{\pi k} g_k (t_{\pi 2},\ldots ,t_{\pi n})
$$
The same procedure applied to $\sum t_j g_{jj}$ produces
$$
\sum t_{\pi k} g_{kk}(t_{\pi 2}, \ldots , t_{\pi n}) ,
$$
so by comparing these we obtain that in fact $\sum b_{\pi k} g_k(t) -\sum t_k g_{kk}(t)=0$ for all permutations $\pi$ (this worked because the $b_j$'s, unlike the $t_j$'s, do not get reordered when we relabel the variables at the end). Thus we also have that
$$
\sum (b_{\pi k}-b_k) g_k = 0 . \quad\quad\quad\quad (1)
$$
If (1) actually had non-zero coefficients, then, since the permutation $\pi$ is arbitrary, I again obtain $v\cdot\nabla g$ on an at least $(n-2)$ dimensional subspace of vectors $v$, so $g$ could only depend on $\sum c_j t_j$, and, being a symmetric function, it really would have to depend on $\sum t_j$ only, but this isn't working.
So the conclusion is that we must have $b_2=\ldots =b_n$, and thus also $a_2=\ldots =a_n$. I can repeat the whole argument with another variable taking the role of $x_1$ to see that in fact (if $n\ge 3$, but $n=2$ is easy to do directly) $a_1=\ldots =a_n$. So, to sum this up, the only possibility left open at this point is $E(f)=\Delta f=\sum f_j=0$. I don't think there can be such an $f$, but it's apparently not quite as easy as I originally thought and I'm running out of steam now, so I'll leave it at that for now. (Perhaps the way to go is to exploit the extra symmetry of $g$ that comes from swapping $x_1$ with another variable in $f$; this produces lots of identities.)
Addendum: I have an argument now, but it's very messy and there must be a better, more insightful way of doing this. Here's a very quick outline: as suggested above, swap $x_1$ and $x_2$, say. This shows that
$$
g(t_2,\ldots, t_n) = g(-t_2,t_3-t_2,\ldots , t_n-t_2) .
$$
Use this on both sides of $\sum t_j g_{jj} = c\sum g_j$ to produce identities for $g_k$. After some fooling around with these and with the help of Euler's identity, I finally arrived at
$$
\left( |x|^2 - n\overline{x} \right)^2 f_j = d\left( x_j-\overline{x}\right) f ,
$$
with $\overline{x} = (x_1+\ldots + x_n)/n$, and $d$ is the degree of $f$ and $g$. Take the derivative wrt $x_j$ on both sides and then sum over $j$. Use that $\Delta f=0$ and $\sum x_j f_j = df$. After a calculation, this shows that $f=0$.
