Distortion is something very, very difficult to predict and compute analytically because it involves non-linear equations and models. What you are asking is not trivial; it's something that has been keeping many great minds occupied for decades. Some very important mathematical schemes have been developed, for example, by using Volterra series it is possible to model the non-linearities of amplifiers including their memory effects (the inclusion of the effect that capacitive and/or inductive elements have in the signal). They are quite cumbersome to do by hand though. From Volterra Series an interesting tool has been developed, which has applications in high-frequency amplifiers, they are called the "X Parameters". RF/Microwave power amplifiers are purposely heavily driven into their non-linear range (aka make them distort), so analyzing their distortion behavior is really important.
Audio doesn't work this way since you don't want to purposely drive the amp into distortion as in the RF world (unless you are building a Fuzz pedal or something), so, perhaps it could be more interesting and relevant to analyze audio amplifiers in what some call "weak non-linearities", that is, if you think of 'non-linear' as a clipped sine wave almost turned into a square wave, and 'linear' as a pristine sine wave with no distortion, 'weakly non-linear' is somewhere in between, but much more tilted towards the pristine sine wave side response than to the square wave. Those sort of models make use of polynomials which have low-valued coefficients for the quadratic and higher order terms (hence its name weakly non-linear, since the linear terms dominate).
In my opinion, as Abbey suggests, all calculations must be made in the computer world in tandem with real measurements. You can propose a weak non-linear model, measure a particular op-amp vs loading and try to fit your model into the measurements. It is more of an academic exercise to study how a particular op-amp distorts rather than something practical for everyday DIYing. In the case of loading, the non-linear components would appear primarily due to the output stage distortion. Again, it is very interesting if you are someone like me who is into academia or research, don't know how useful this would be to a DIYer though.
That being said, some op-amp datasheets do provide measurements of THD vs loading, they are not very comprehensive, usually 600 ohms and 2K will be the reported results. As a general rule, you don't want to go below 600 ohms, and 600 ohms is where an op-amp will distort more when considering full voltage swing. Of course, if you are not outputting full voltage swing, you can go much lower than 600 ohms.
I guess that, for some raw calculations, what you can do is to check the THD measurements and look at the swing values they used, for example, if they used 600 ohms and +/- 12V swing, you can divide 12/600 = 20 mA, and extrapolate to different values of loading/voltage. For example, if you use 1 V into 50 ohms, 1/50 = 20 mA, you are using the same current, and the distortion values MIGHT be similar to the values reported for 12V and 600 ohms. I emphasize MIGHT since I already mentioned that there are other factors being involved and the linear relationship I just described might not be valid. Measurements must be carried out to check if your predictions were right or completely off-target.