This can be evaluated numerically using the Euler-Maclaurin formula. Let $t(n) = (x\log x)^{-1}$ be the $n$-th term, pick a threshold $a$ and write the sum as $ \sum_{2\leq k < a}t(n) - \log\log a+\frac12t(a) + \sum_{k\geq 1}R_k, $ $ R_k = \frac{B_{2k}}{(2k)!}\partial_n^{2k-1}\big|_{n=a}t(n). $ The sum of $R_k$ is asymptotic in nature and diverges as $k\to\infty$, so it is necessary to sum $R_k$ only until the point $k$ when $|R_k|$ starts increasing.
Here is the value to what should be about $200$ digits: $\begin{array}{rrrrr} 0.7946786454 & 5289940220 & 3897962065 & 1495140649 & 9959088280 \\ 4968901512 & 0950148178 & 5896068756 & 6696614777 & 6273344714 \\ 3933647273 & 2836173295 & 5585193860 & 6587583982 & 7169312934 \\ 2858106252 & 2997064914 & 8599822651 & 4777778609 & 7884785692 \end{array}$ This was produced with the following code:
from mpmath import * @memoize @monitor def fd(n, x, a=1, b=1): # n-th derivative of 1/(x**a*log(x)**b) if n < 0: raise ValueError('fd must have n >= 0') dp = [[1/(x**(i+a)*log(x)**(j+b)) for j in range(n+1)] for i in range(n+1)] for k in range(1,n+1): for i in range(n+1-k): for j in range(n+1-k): dp[i][j] = -(b+j)*dp[i+1][j+1]-(a+i)*dp[i+1][j] return dp[0][0] @monitor def qq(a = None): if a is None: a = int(2 + ceil(mp.dps / 2.7)) ax = mpf(a) def term(x): x = mpf(x) return 1/(x*log(x)) s = sum(term(x) for x in range(2, a)) - log(log(ax)) + term(a)/2 b, bPrev, bCount = None, None, 0 for k in range(1, 1000): b = bernoulli(2*k)/gamma(2*k+1) * fd(2*k-1, ax) if (bPrev is not None and abs(b) > abs(bPrev)) or abs(b) < mpf(10)**(-mp.dps): break bPrev = b s -= b bCount += 1 print('Precision %.3f (%.3f) from %d terms at a=%d' % (log10(abs(bPrev)), log10(abs(b)), bCount, a)) return s