Using fft to calculate vsdft
I believe the short answer is no.

The paper is behind a paywall so I didn't look at it.  In the abstract is the following blurb: " use interpolation techniques to resample the signal at constant angular increments, followed by a common discrete Fourier transform to shift from the angular domain to the order domain."  In which they are describing other techniques.

Looking at your equations, that is precisely what needs to be done to $\theta[n]$ in order to convert it into a regular DFT, suitable for a FFT.

The FFT exploits symmetries in the even spacing of the exponents in order to save computations.  With a varying $\theta[n]$, those symmetries aren't there so you aren't going to be able to implement a FFT algorithm.

You don't tell us what $x[n]$ is either.

How large is N?

Hope this helps at least a little.

Ced