The posted code doesn't run. The frequency vector f is defined in terms of numel(X), but X isn't defined yet.
I thing that introducing the zero padding with nextpow2 on only one side of the signal introduces some phase distortion. The code below doesn't use zero padding.
The values of X won't be exactly real because of computational innacuruacy. However, you can check that the imaginary part is very, very small relative to the real part in frequencies of interest, so the phase will be either 0 or +-180. But then unwrapping that causes problems. Don't unwrap.
x = exp(-(1+1i*0)/2*(t/10).^6);
f = linspace(-fs/2*1e3,fs/2*1e3,nfft);
X = fftshift(fft(x,nfft));
If you're confident that X should be real, then maybe consider