Using fft to calculate vsdft We are currently implementing the velocity synchronous discrete Fourier transform (vsdft) for order
tracking rotating machinery in Python. 

The formula is similar to dft. Simplified, it is given by

$$ \text{VSDFT}(\Omega) = C\sum_{n=0}^{N-1} x[n] \omega[n] e^{j\Omega \theta[n]} $$

where $C$ is a constant and 

$$ \omega = \frac{d\theta}{dt} $$

Since dft is given by 

$$ F_k = \sum_{n=0}^{N-1} f_n\cdot e^{-i 2 \pi k n / N} $$

and we wanted to use the optimized fft implementing of numpy, and was thinking if there is a way to write the input f so that we could
use fft directly. However, my math skills are failing me.

Can any of you see a way or any pointer to implement the above VSDFT without having to reimplement the Cooley-Tukey FFT algorithm? 

Refrence to the original paper can be found in here.
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