(Skip the initial two – or even three – paragraphs for the actual solution).
The equation $$x_{n+1} - x_n = \nabla F(x_n) - x_n$$ seems to be a discretised version of $$x'(t) = \nabla F(x(t)) - x(t).$$ The vector field $\nabla F(x) - x$ is conservative, with scalar potential $|x|^2/2 - F(x)$ which is confining (that is, it goes to $\infty$ at infinity). Thus, the gradient flow will converge to the set of critical points of the potential. In a typical situation, this set is discrete, in which case the answer for the continuous problem is affirmative (although the limit need not be unique: it may depend on the initial value). If, however, you choose a sufficiently nasty potential, the gradient flow may fail to converge.
Now let us try a similar approach for the discrete equation: write $G(x) = |x|^2/2 - \nabla F(x)$, so that $x_{n+1} - x_n = -\nabla G(x_n)$. This is the gradient descent sequence for the function $G$, with a constant parameter. It is well-known that gradient decent may fail to converge if the parameter is too large, compared to the second derivatives of $G$. However, in our case the second derivative of $G$ in any direction is bounded by $1$, which is precisely what we need.
Here the actual solution begins: Write $G(x) = |x|^2/2 - F(x)$, so that $x_{n+1} = x_n - \nabla G(x_n)$. By convexity of $F$, $$ G(x) \le G(x_n) + \nabla G(x_n) \cdot (x - x_n) + |x - x_n|^2/2 . $$ Setting $x = x_{n+1}$ leads to $$ \begin{aligned} G(x_{n+1}) & \le G(x_n) + \nabla G(x_n) \cdot (-\nabla G(x_n)) + |\nabla G(x_n)|^2/2 \\ & = G(x_n) - |\nabla G(x_n)|^2/2 . \end{aligned} $$ It follows that $(G(x_n))$ is a (stricly) decreasing sequence, and therefore it converges to a certain limit. In other words, $\sum_n |\nabla G(x_n)|^2 < \infty$, that is, $\sum_n |x_{n+1} - x_n|^2 < \infty$. In particular, the set of accumulation points of $(x_n)$ is connected. Additionally, assuming continuity of $\nabla F$, for every accumulation point $x^*$ we have $\nabla G(x^*) = 0$. To summarise:
If $\nabla F$ is continuous and the set of critical points of $|x|^2/2 - F(x)$ is discrete, then the answer is affirmative: $(x_n)$ converges to a critical point of $F$.
OK, I cheated: the actual solution starts here. What happens if the set of critical points of $G(x) = |x|^2/2 - F(x)$ is not discrete? It turns out that $(x_n)$ may fail to converge! It is difficult to describe a counterexample, and I did not check the details; roughly speaking, the graph $G(x)$ described below is an infinite spiral ramp around the unit disk.
We suppose that $d = 2$, we set $c > 0$ small enough (to be determined later), and we identify $(x, y)$ with $x + i y$. The curves $(1 + e^{-t}) e^{i t}$ and $(1 + 2 e^{-t}) e^{i t}$ are two spirals that will describe the edges of the ramp. The height of the ramp will be roughly $e^{-2 t}$, so the slope will decrease rapidly as the ramp approaches the unit circle.
We define $G(r e^{i t})$ in the following way. First of all, we set $G(r e^{i t}) = 0$ when $r \le 1$. For a fixed $t > 0$, if $1 + e^{-t} \le r < 1 + 2 e^{-t}$, we let $G$ to be a linear function of $r$: $$ G(r e^{i t}) = c (e^{-2 t} + a(t) (r - 1 - e^{-t})) , $$ where $a(t)$ (extremely small) is chosen in such a way that for every $t > 0$, the point $$ G((1 + e^{-t}) e^{i t}) - \nabla G((1 + e^{-t}) e^{i t}) $$ also lies on the same spiral $(1 + e^{-s}) e^{i s}$. A quick-and-dirty calculation suggests that there is a solution $a(t)$ such that the corresponding $s$ satisfies $s \approx t + c e^t a(t)$, and we have $a(t) \sim 2 e^{-3t}$ (and similarly for $a'(t)$ and $a''(t)$); but here I must admit that I did not really check the details.
For $t > 2 \pi$ we define $G(r e^{i t})$ when $1 + 2 e^{-t} < r < 1 + e^{-t + 2 \pi}$ as a cubic polynomial in $r$ that matches the values and the derivatives of $G$ with respect to $r$ at the boundary values $r = 1 + 2 e^{-t}$ and $r = 1 + e^{-t + 2 \pi}$. Another quick-and-dirty calculation shows that wherever we have defined $G$ so far, the second derivatives of $G$ in any direction are bounded by a constant multiple of $c$. We finally extend $G$ to all of $\mathbb{R}^2$ so that $G(r e^{i t}) = 0$ for $r \ge 10$, and the second derivatives of $G$ are everywhere bounded by a constant multiple of $c$.
Now we are ready to eventually choose $c$: we make it small enough, so that the second derivatives of $G$ are in fact bounded by $1$. And we eventually define $F(r e^{i t}) = r^2/2 - G(r e^{i t})$ and $F(r e^{i t}) = F_0(r e^{i t})$ when $r \le 10$.
Now $F$ is convex (because the second derivatives in any direction belong to the interval $(0, 2)$) and $F(r e^{i t}) = r^2 / 2$ for $r \ge 10$. There is one more step needed to assert that $\nabla F$ is bounded: we modify $F(r e^{i t})$ for $r \ge 10$ and we set $F(r e^{i t}) = 10 (r - 10)$ for $r \ge 10$ (this does not affect convexity of $F$).
By the above construction, if $x_0 = 1 + e$, then $x_n = (1 + e^{-t_n}) e^{i t_n}$ for sequence $(t_n)$ such that $t_{n+1} - t_n \to 0$. On the other hand, we must have $t_n \to \infty$, because $\nabla F((1 + e^{-t}) e^{i t}) \ne 0$ for all $t > 0$. Therefore, $(x_n)$ is not convergent.
The above construction is not clear, I know. If you are able to simplify it, feel free to edit this answer.