Assume, that FFT tells you that signal is

SIGNAL = a1 x sin (wt) + b1 x cos (wt) + a2 x sin (2wt) + b2 x cos (2wt) + a3 x sin (3wt) + b3 x cos (3wt) + ...

then

SIGNALRMS = SQRT(1/2 x (a1^2 + b1^2 + a2^2 + b2^2 + a3^2 + b3^2 + ...))

and

FUNDRMS = SQRT(1/2 x (a1^2 + b1^2))

So, as you intend to calculate SIGNALRMS numerically, means not by using the FFT, then you should calculate FUNDRMS as above and

THD = SQRT(SIGNALRMS^2-FUNDRMS^2)*100/SIGNALRMS

will finally do the trick.

By the way, if you synchronize the taking of samples with the zero crossings, "b1" will vanish eventually. Check this with a simulation.

Kai