Reid Barton's approximation
$$
f(n,2) \simeq R(n) := \frac12 \sum_{k\in{\bf Z}} \phantom. 2^k n \exp(-2^k n) \phantom{\infty\infty}(n\rightarrow\infty)
$$
is corroborated by further numerical computation (see plot below for $2^6 \leq n \leq 2^{13}$), and accounts for both the averaged limit of $.72134752\ldots$, which arises as $1/(2 \log 2)$, and the residual oscillation, whose limiting form as a function of $t := \log_2 n$ is a nearly perfect sine wave of period $1$ and amplitude
$$
\frac1{\log 2}\left|\phantom.\Gamma\bigl(1 + \frac{2\pi i}{\log 2}\bigr)\right| \phantom.=\phantom.
\frac1{\log 2}\left[\frac{(2\pi^2/ \log2)}{\sinh(2\pi^2/ \log2)}\right]^{1/2} = \phantom. 7.130117\ldots \cdot 10^{-6}.
$$
This is obtained by applying Poisson summation to
$$
\varphi(x) := 2^x \exp(-2^x),
$$
which as Reid observed decays rapidly both as $x \rightarrow \infty$ and as $x \rightarrow -\infty$.  We have
$$
R(n) = \frac12 \sum_{k \in {\bf Z}} \phantom. \varphi(t+k),
$$
a smooth function of period $1$ in $t = \log_2 n$, which thus has a Fourier expansion
$$
R(n) = \sum_{m\in{\bf Z}} \phantom. a_m e^{2\pi i m t}
$$
where
$$
a_m = \frac12 \int_0^1 R(2^t) \phantom. e^{-2\pi i m t} dt
= \frac12 \int_{-\infty}^\infty \varphi(t) \phantom. e^{-2\pi i m t} dt
= \frac12\hat\varphi(-m).
$$
Changing the variable of integration from $t$ to $2^t$ lets us recognize the Fourier transform $\varphi$ as $1/\log(2)$ times a Gamma integral:
$$
\varphi(y) = \frac1{\log 2} \Gamma\Bigl(1 + \frac{2 \pi i y} {\log 2}\Bigr).
$$ 
This gives us the coefficients $a_m$ in closed form.  The constant coefficient $a_0 = 1 / (2\log 2)$ can be interpreted as the approximation of the Riemann sum $R(n)$ by an integral; the oscillating terms $a_m e^{2\pi i m t}$ for $m \neq 0$ are the corrections to this approximation, and are small due to the exponential decay of the Gamma function on vertical strips; indeed we can compute the magnitude $|a_m|$ in elementary closed form using the formula
$|\Gamma(1+i\tau)| = (\pi\tau / \sinh(\pi\tau))^{1/2}$.  Taking $m = \pm 1$ recovers the amplitude $2|a_1|$ exhibited above; the $m = \pm 2$ terms yield an oscillation twice as fast but with magnitude $6.6033857\ldots \cdot 10^{-12}$, and further terms are smaller yet.

The plot below shows $R(n)$ for $2^6 \leq n \leq 2^{13}$ and compares it with actual values $f(n,2)$.  (See http://math.harvard.edu/~elkies/mo11255.pdf for the original PDF plot, which can be "zoomed in" to view details.)  The horizontal coordinate is $\log_2 n$; the vertical coordinate is centered at $1/(2 \log 2)$, showing also the lines $1/(2 \log 2) \pm 2|a_1|$; the red contour shows $R(n)$; and the black cross-hairs, eventually merging visually into a continuous curve, show the numerical values of $f(n,2)$.

<img src="http://math.harvard.edu/~elkies/mo11255.png">

Numerical computation of the difference $f(n,2) - R(n)$ suggests that it oscillates with the same period but with decreasing amplitude asymptotic to $C/n$ for some constant $C$; the blue curves show this with the numerical estimate $2.947 \cdot 10^{-4}$ for $C$.

Much the same behavior should hold for other $k$, with even tighter oscillations; e.g. for $k=6$ the Gamma factor in $a_1$ would be $\Gamma(1 + \frac{2\pi i}{\log 1.2})$ which has magnitude about $4.55 \cdot 10^{-23}$.