Integration of a rapidly oscillating function

A bit of experimentation indicates that "LevinRule" does not work well, and that "MaxErrorIncreases" and MaxRecursion must be increased. Most importantly, WorkingPrecision must be increased substantially. To do so, first increase the precision of the only finite precision constant.

alpha = SetPrecision[7.2973525664*^-3, 65];

Then, set the options of NIntegrate as follows.

NIntegrate[rho[r], {r, 0, 55000}, Method -> {"GlobalAdaptive", 
    "SymbolicProcessing" -> 0, "MaxErrorIncreases" -> 2000}, MaxRecursion -> 20, 
    WorkingPrecision -> 60, PrecisionGoal -> 8, AccuracyGoal -> Infinity]
(* 0.999999999862974654189196482270003235458935716234694791163667 *)

as desired. AbsoluteTiming for the computation is 1.4 seconds.

The source of the problem

The problem is the expansion of LaguerreL[], with degree nr == 118. If it is evaluated with a symbolic argument r, it will expand to power series form (the usual linear combination of powers of r); in this form, it is impossible to evaluate accurately at machine precision. A very high precision is needed, and in addition, the machine precision coefficients, which derive from the definition of alpha, need to have their precision increased. I set alpha = 7.2973525664`100*^-3 for the following. For example, compare the plots:

Plot[rho[r], {r, 0, 55000}, PlotRange -> All, 
 PlotLabel -> "Evaluated on a numeric argument"]

Plot[Evaluate@rho[r], {r, 0, 55000}, PlotRange -> All, 
 PlotLabel -> "Evaluated on a symbolic argument, with low precision"]

Plot[Evaluate@Rationalize[rho[r], 0], {r, 0, 55000}, PlotRange -> All,
  WorkingPrecision -> 50, 
 PlotLabel -> "Evaluated on a symbolic argument, with insufficient precision"]

Plot[Evaluate@Rationalize[rho[r], 0], {r, 0, 55000}, PlotRange -> All, 
 WorkingPrecision -> 60, 
 PlotLabel -> "Evaluated on a symbolic argument, with high precision"]

Here we can see why WorkingPrecision -> 60 was successful in @bbgodfrey's answer.

An easy fix to the problem

The trick is to evaluate LaguerreL[] only on numeric arguments using ?NumericQ to control the evaluation. The built-in routine will evaluate LaguerreL[] accurately on machine precision inputs. This may be used closest to the source of the problem in up[] and um[] or simply on the top-level function rho[]:

(* with alpha = 7.2973525664`*^-3, machine precision as in the OP *)
ClearAll[rho];
rho[r_?NumericQ] := Abs[g[r]]^2 + Abs[f[r]]^2; 
NIntegrate[rho[r], {r, 0, Infinity}, MinRecursion -> 3, MaxRecursion -> 20] // AbsoluteTiming
(*  {0.727109, 1.0000000000001033`}   *)

I wouldn't consider this an oscillatory integral, since it eventually stops oscillating. (The Levin rule expects an oscillatory function like Sin[] or BesselJ[] etc.) However it does have many oscillation, so a sufficiently fine sampling is necessary, which is the purpose of MinRecursion -> 3 above. This could also be achieved with a high number of sample points, such as Method -> {"GaussKronrodRule", "Points" -> 61}.

Splitting the interval up can speed things up in this case, since the function effectively has support over a finite interval and is analytic. Something like this takes only 0.32 sec.:

NIntegrate[rho[r], {r, 0, 50000, Infinity}, Method -> {"GaussKronrodRule", "Points" -> 21}]

One final remark: The integral can be evaluated symbolically; but it takes a really long time, and the numerics issue of the expanded polynomial persists. Below is a relatively fast way, since the integrand, after removing the unnecessary Abs, expands into a linear combination of terms of the form Exp[-cnt r] * r^s for various powers s.

alpha = 7.2973525664`200*^-3;
(* ... re-execute OP's defs ...*)

i[s_] = Integrate[E^(-2 cnk r) r^s, {r, 0, Infinity}, Assumptions -> s > 0];

List @@ (rho[r] /. Abs -> Identity // Expand) /. E^(-2 cnk r)*r^s_ :> i[s] // Total
(*1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000*)