Algorithms for computing FFT in parallel
I think your main problem isn't how to parallel the algorithm (which actually can be done) but it's numeric precision. 
FFTs of a size that large are numerically quite tricky. The FFT coefficients are of the form
$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$ and if N is very large the coefficient calculation gets noisy. Lets say you have $N = 2^{40}$ and you use 64 bit double precision arithmetic. The first 1000 coefficients have a real part that's exactly unity (although it shouldn't be that way) so you will need higher precision math, which is very inefficient and cumbersome to use. 

You will also amass lots of rounding and truncation errors as the sheer number of operations that go into a single output number is very large too. Due to the "every output depends on every input" nature of the FFT, error propagation is rampant.

I'm not aware of an easy way to work around that. Your request is an unusual one. Most applications that do spectral analysis of large data sets do a running analysis where you don't have that problem. Perhaps if you can describe your application and it's constraints a but more, we can point you to a more suitable solution.