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