Interchanging integration and summation vs. convergence

I'm trying to calculate the series (all numerical data is calculated in Mathematica)

$\displaystyle S=\sum_{n=0}^{\infty}{\frac{1}{(n+1)(2n+1)(3n+1)}} \approx 1.05964$

I proceeded as follows. First the partial fraction decomposition of the summand,

$\displaystyle {\frac{1}{(n+1)(2n+1)(3n+1)}=\frac{1}{2}\frac{1}{n +1}-4\frac{1}{2n+1}+\frac{9}{2}\frac{1}{3n+1}$,

then introducing the variable x:

$\displaystyle \frac{1}{(n+1)(2n+1)(3n+1)}=\frac{1}{2}\frac{x^{n+ 1}}{n+1}-4\frac{x^{2n+1}}{2n+1}+\frac{9}{2}\frac{x^{3n+1}}{ 3n+1}|_0^1$.

This is then recognized as the integral

$\displaystyle S=\sum_{n=0}^{\infty}{\int_0^1{(\frac{1}{2}x^n-4 x^{2n}+\frac{9}{2}x^{3n})dx}}$

This is still numerically equal to the first expression. Now, I want to interchange the summation and integration and sum the geometric series,

$\displaystyle S=\int_0^1{(\frac{1}{2}\frac{1}{1-x}-4 \frac{1}{1-x^2}+\frac{9}{2}\frac{1}{1-x^3})dx}$,

but this integral is numerically equal to $\displaystyle 0.7980$.

Why aren't I allowed to sum the geometric series?