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.

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 every detail; roughly speaking, the graph of the function $G(x)$ described below is an infinite spiral ramp around the unit disk.

[Edit: I re-wrote the following construction.] We suppose that $d = 2$. We use polar coordinates $(r, t)$ rather than Cartesian ones $(x, y)$. We also write $P_n$ (rather than $x_n$) for the sequence of points defined in the question: $P_{n+1} = P_n - \nabla G(P_n)$.

Step 1. We divide $\mathbb{R}^2$ into five regions:

• 'interior': the unit disk, $D_0 = \{r \le 1\}$;

• 'ramp': the area between two spirals, $D_1 = \{1 + e^{-t} \le r \le 1 + 2 e^{-t} , \, t \ge 0\}$;

• 'walls': everything between the consecutive turns of the ramp, $D_2 = \{1 + 2 e^{-t} < r < 1 + e^{-t + 2 \pi} , \, t \ge 2 \pi\}$;

• 'circus': $D_3 = \{r < 10\} \setminus (D_0 \cup D_1 \cup D_2 \cup D_3) = \{1 + 2 e^{-t} < r < 10 , \, 0 \le t < 2 \pi\}$;

• 'exterior': $D_4 = \{r \ge 10\}$.

Our goal is to define $G$ in such a way that: (A) the second order derivatives of $G$ are small; (B) a particle placed on the ramp and following the recurrence equation $P_{n+1} = P_n - \nabla G(P_n)$ stays on the ramp forever.

Step 2. First of all, we set $G = 0$ in the interior $D_0$ and the exterior $D_4$.

Step 3. We now define $G$ on the ramp $D_1$ in such a way that condition (B) is satisfied: for a fixed $t \ge 0$, if $1 + e^{-t} \le r \le 1 + 2 e^{-t}$, we let $G$ to be a linear function of $r$: $$ G = c (e^{-2 t} + a(t) (r - 1 - e^{-t})) , $$ where $a(t) > 0$ is extremely small. More precisely, we require that if the point $P$ lies on the inner edge of the ramp: $r = 1 + e^{-t}$, then the point $$ P' = P - \nabla G(P) $$ also lies on the same curve, a little bit down the ramp. Denote the polar coordinates of $P$ by $t$ and $r = 1 + e^{-t}$, and the polar coordinates of $P'$ by $t'$ and $r' = 1 + e^{-t'}$. At $r = 1 + e^{-t}$ we have $\partial_r G = c a(t)$ and $\partial_t G = -2 c e^{-2 t} + c a(t) e^{-t}$. Our condition reads $$ r' \cos(t' - t) = r - \partial_r G(P) , \quad r' \sin(t' - t) = -r \partial_t G(P) , $$ that is, $$ (1 + e^{-t'}) \cos(t' - t) = (1 + e^{-t}) - c a(t) , \quad (1 + e^{-t'}) \sin(t' - t) = c (1 + e^{-t}) (2 e^{-2 t} - a(t) e^{-t}) . $$ Given a small $c > 0$, by the inverse mapping theorem, it looks like we can choose $t'$ and $a(t)$ in such a way that the above condition holds, and $$ \begin{gathered} t' - t \sim 2 c e^{-2 t} , \\ a(t) \sim 2 e^{-3t} , \\ a'(t) \sim -6 e^{-3 t} , \\ a''(t) \sim 18 e^{-3 t} . \end{gathered} $$ (I did not work out the details here). The above properties of $a(t)$ imply that the second order derivatives of $G$ are bounded by a constant times $c$ on $D_1$.

Step 4. We now build the walls for our ramp: we define $G$ in $D_2$. For $t \ge 2 \pi$ and $1 + 2 e^{-t} < r < 1 + e^{-t + 2 \pi}$ we let $G$ to be the cubic polynomial in $r$ which matches the values and the derivatives of $G$ with respect to $r$ at the boundary values $r_0 = 1 + 2 e^{-t}$ and $r_1 = 1 + e^{-t + 2 \pi}$, defined already in Step 3. The formula for $G$ can be given explicitly in terms of $t$, $a(t)$ and $a(t - 2 \pi)$. The important thing is that the second derivative of $G$ with respect to $r$ is bounded by the ratio of the 'height of the wall' (which is at most $c e^{-2 t} (e^{2 \pi} - 1)$) to the square of the 'width of the wall' (which is $e^{-t} (e^{2 \pi} - 2)$) plus the ratio of $a(t)$ to the 'width of the wall'. Thus, it is bounded by a constant times $c$. Similarly, one can bound second order derivatives of $G$ in any direction at all points of $D_2$ (and again I did not work out the details).

Step 5. We need to extend $G$ to the circus $D_3$. Since this is a nice domain, and we have already defined $G$ elsewhere so that it has second order derivtives bounded by a constant times $c$, we can define $G$ in $D_3$ in such a way that the bound on the second order derivatives remains the same (possibly with a larger constant).

Step 6. 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^2/2 - G$.

Now $F$ is convex (because the second derivatives of $F$ 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 of course does not affect convexity of $F$.

End of the construction. By property $(B)$, if $P_0$ corresponds to $t = 0$ and $r = 1 + e$ (which lies on the inner edge of the ramp $r = 1 + e^{-t}$), then $P_n$ corresponds to $t_n$ and $r_n = 1 + e^{-t_n}$ for some $t_n > 0$. Additionally, $t_{n+1} - t_n = t_n' - t_n \sim 2 c e^{-2 t_n}$, and hence $t_{n+1} - t_n \to 0$ and $t_n \to \infty$. Therefore, $(P_n)$ is not convergent.

I am still not satisfied with the above description (and the quality of the image). If you are able to simplify it, feel free to edit this answer.

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:

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 = 1}^\infty |x_{n+1} - x_n|^2 < \infty$. In particular, $$\lim_{n \to \infty} |\nabla G(x_n)| = \lim_{n \to \infty} |x_{n+1} - x_n| = 0.$$ This implies that the set of accumulation points of $(x_n)$ is connected (see below). Additionally, assuming continuity of $\nabla F$, for every accumulation point $x^* = \lim_{n \to \infty} x_{k_n}$ we have $\nabla G(x^*) = \lim_{n \to \infty} \nabla G(x_{k_n}) = 0$. To summarise:

Let me remark why the set of accumulation points of a bounded sequence $(x_n)$, satisfying the condition $\lim_{n \to \infty} |x_{n+1} - x_n| = 0$, is connected. Suppose that it is not, that is, it can be divided into two closed parts $F_1$ and $F_2$, whose distance is positive. Let $3 \delta$ denote the distance between $F_1$ and $F_2$. Choose and arbitrary $N > 0$ large enough, so that $|x_{n + 1} - x_n| < \delta$ when $n > N$. There exist elements $x_i, x_k$ such that $N < i < k$, $\operatorname{dist}(x_i, F_1) < \delta$ and $\operatorname{dist}(x_k, F_2) < \delta$. For $j = i, i+1, \ldots, k$ we have $\operatorname{dist}(x_j, F_1) + \operatorname{dist}(x_j, F_2) \ge 3 \delta$, each summand changes by at most $\delta$ when $k$ is replaced by $k + 1$. It is therefore easy to see that for some $k \in \{i + 1, i + 2, \ldots, k - 1\}$ we must have $\operatorname{dist}(x_k, F_1) \ge \delta$ and $\operatorname{dist}(x_k, F_2) \ge \delta$. In other words, for any $N$ there is $k > N$ such that $\operatorname{dist}(x_k, F_1) \ge \delta$ and $\operatorname{dist}(x_k, F_2) \ge \delta$. By choosing a convergent subsequence of these elements $x_k$, we find an accumulation point of $(x_n)$ away from $F_1 \cup F_2$, a contradiction.