We intend to find values of the $n$-the point $t$ such that:
$\Re\left(\zeta \left(i t+\frac{1}{2}\right)\right) \neq 0$
and
$\Im\left(\zeta \left(i t+\frac{1}{2}\right)\right) = 0$
which should be what are known as the Gram points.
We start with $(1)$ the Franca-LeClair asymptotic (A slightly changed version of formula (163) page 47 in "A theory for the zeros of Riemann ζ and other L-functions" by Guilherme Franca and Andre LeClair, July 2014) - entered by Peter Bala in the OEIS, to the Riemann zeta zeros:
$(1)$ $$x=2 \pi \exp (1) \exp \left(W\left(\frac{n-\frac{11}{8}}{\exp (1)}\right)\right)$$
where $W(z)$ above is the Lambert W function.
We then feed the right hand side (RHS) into:
$(2)$ $$\frac{\vartheta \left(x\right)}{\pi}$$
where $\vartheta \left(x\right)$ is the Riemann–Siegel theta function and see what happens, and notice that it appears that approximately:
$(3)$ $$k-\frac{\vartheta \left(x\right)}{\pi }-\frac{3}{2} \approx 0$$
(4):
$$k-\frac{\vartheta \left(2 \pi \exp (1) \exp \left(W\left(\frac{n-\frac{11}{8}}{\exp (1)}\right)\right)\right)}{\pi }-\frac{3}{2} \approx 0$$
with the left hand side $\text{(LHS)}$ being numerically:
$\text{LHS} =-0.000456796..., -0.000321091..., -0.000260155..., -0.000223000..., -0.000197231...,...$
for $k=1,2,3,4,5,...$ a natural number.
to remedy the situation and make it approximately equal to $n$ again, we now take formula (226) at page 65 in "A theory for the zeros of Riemann ζ and other L-functions" by Guilherme Franca and Andre LeClair, July 2014, which is:
$(5)$
$$\frac{\widetilde{y}_{n}}{2\pi}\log\left(\frac{\widetilde{y}_{n}}{2\pi e}\right)=n-\frac{11}{8}$$
(This of course resembles the Riemann–von Mangoldt formula.) We then change sides of some of the terms and get:
$(6)$ $$n=\frac{\widetilde{y}_{n}}{2\pi}\log\left(\frac{\widetilde{y}_{n}}{2\pi e}\right)+\frac{11}{8}$$
The searched for remedy is now that we add the LHS of $(4)$ to the RHS of $(6)$ to make the error smaller:
$(7)$ $$n \approx \frac{x}{2 \pi }\log \left(\frac{x}{2 \pi e}\right)+ \frac{11}{8}+k-\frac{\vartheta (x)}{\pi }-\frac{3}{2}$$
where $\widetilde{y}_{n}$ is substituted with $x$. And $\frac{11}{8}-\frac{3}{2} = \frac{1}{8}$ so $(7)$ becomes:
$(8)$ $$n \approx \frac{x}{2 \pi }\log \left(\frac{x}{2 \pi e}\right)+ \frac{1}{8}+k-\frac{\vartheta (x)}{\pi }$$
For easier mouse scrolling we repeat $(1)$:
$(1)$ $$x=2 \pi \exp (1) \exp \left(W\left(\frac{n-\frac{11}{8}}{\exp (1)}\right)\right)$$
and guess that the RHS of $(8)$, could fit into the place of $n$ in $(1)$ so that we get $(9)$:
$(9)$ $$x=2 \pi \exp (1) \exp \left(W\left(\frac{\frac{x}{2 \pi }\log \left(\frac{x}{2 \pi e}\right)+ \frac{1}{8}+k-\frac{\vartheta (x)}{\pi }-\frac{11}{8}}{\exp (1)}\right)\right)$$
Now: $\frac{1}{8}-\frac{11}{8}=-\frac{5}{4}$
$(10)$ $$x=2 \pi \exp (1) \exp \left(W\left(\frac{\frac{x}{2 \pi }\log \left(\frac{x}{2 \pi e}\right) -\frac{5}{4}+k-\frac{\vartheta (x)}{\pi }}{\exp (1)}\right)\right)$$
But: This should, but does not converge to a Franca-LeClair point. As in the comment below, the definition of Franca-LeClair point is:
$$\Re\left(\zeta \left(i t+\frac{1}{2}\right)\right) = 0$$ $$\Im\left(\zeta \left(i t+\frac{1}{2}\right)\right) \neq 0$$
and the iterative formula $(11)$ converges to a Franca-LeClair point (my own wording), for $c=\frac{1}{2}$.
However, and to the main question, we do find Gram point experimentally by setting the constant $-\frac{5}{4}$ to zero $c=0$. This then gives us the iterative formula:
$(11)$ $$x=2 \pi \exp (1) \exp \left(W\left(\frac{\frac{x}{2 \pi }\log \left(\frac{x}{2 \pi e}\right) -c+k-\frac{\vartheta (x)}{\pi }}{\exp (1)}\right)\right)$$
where as said $c=0$, and where $x$ converges to a Gram point.
A demonstration with Mathematica shows that formula $(11)$ indeed converges to a Gram point:
Clear[x, k, c, numberOfIterations, i, t]
numberOfIterations = 20;
(* Interesting values of c are: c = 0, c = 1/2, c = 1/4, c = 3/4 *)
(* c = 0 gives Gram points *)
(* c = 1/2 gives Franca-LeClair points *)
(* c = 1/4 gives non-zero self intersections: Re[Zeta[1/2+I*t]] = Im[Zeta[1/2+I*t]] *)
(* c = 3/4 gives: Re[Zeta[1/2+I*t]] = -Im[Zeta[1/2+I*t]] *)
c = 0;
t = Table[x = 1;
Do[x = N[Round[
2*Pi*E^1*
E^LambertW[((x/(2*Pi))*Log[x/(2*Pi*E)] - c + k -
RiemannSiegelTheta[x]/Pi)/E^1], 10^-20], 20];, {i, 1,
numberOfIterations}]; x, {k, 0, 12}]
Zeta[1/2 + I*t]
So my question is, what is the correct way to handle the constant $c$ so that one arrives at $c=0$ and thereby the Gram points in the iterative formula $(11)$?
$$x=2 \pi \exp (1) \exp \left(W\left(\frac{\frac{x}{2 \pi}\log \left(\frac{x}{2 \pi e}\right) -c+k-\frac{\vartheta (x)}{\pi}}{\exp (1)}\right)\right)$$
This is not a complete answer either, but it might be worth noticing that the equation in the iterative formula can be put in a different form as follows:
When $x=$ the $k$-th Gram point: $$x=2 \pi \exp (1) \exp \left(W\left(\frac{\frac{x}{2 \pi}\log \left(\frac{x}{2 \pi e}\right) -c+k-\frac{\vartheta (x)}{\pi}}{\exp (1)}\right)\right) $$ the equation can be written: $$\frac{x}{2 \pi \exp (1)}= \exp \left(W\left(\frac{\frac{x}{2 \pi}\log \left(\frac{x}{2 \pi e}\right) -c+k-\frac{\vartheta (x)}{\pi}}{\exp (1)}\right)\right) $$
$$\log\left(\frac{x}{2 \pi \exp (1)}\right)= W\left(\frac{\frac{x}{2 \pi}\log \left(\frac{x}{2 \pi e}\right) -c+k-\frac{\vartheta (x)}{\pi}}{\exp (1)}\right) $$
$$\frac{x}{2 \pi e}\log\left(\frac{x}{2 \pi \exp (1)}\right)= \frac{\frac{x}{2 \pi}\log \left(\frac{x}{2 \pi e}\right) -c+k-\frac{\vartheta (x)}{\pi}}{\exp (1)}$$
$$\frac{x}{2 \pi e}\log\left(\frac{x}{2 \pi e}\right) = \frac{x}{2 \pi e}\log \left(\frac{x}{2 \pi e}\right) + \frac{-c+k-\frac{\vartheta (x)}{\pi}}{e}$$
$$\frac{x}{2 \pi e}\log\left(\frac{x}{2 \pi e}\right) = \frac{x}{2 \pi e}\log \left(\frac{x}{2 \pi e}\right) + \frac{-c+k}{e}-\frac{\vartheta (x)}{\pi e}$$
Series expansion in Mathematica of RiemannSiegelTheta[x] at infinity;
gives:
$$\frac{\vartheta (x)}{\pi e} \approx -\frac{x \left(\log \left(\frac{2 \pi }{x}\right)+1\right)}{2 (e \pi )}-\frac{1}{8 e}+\frac{1}{48 e \pi x}+\frac{7}{5760 e \pi x^3}+\frac{31}{80640 e \pi x^5}+O\left(\left(\frac{1}{x}\right)^7\right)$$
which is:
$$\frac{\vartheta (x)}{\pi e} \approx \frac{x}{2 \pi e}\log\left(\frac{x}{2 \pi e}\right)-\frac{1}{8 e}+\frac{1}{48 e \pi x}+\frac{7}{5760 e \pi x^3}+\frac{31}{80640 e \pi x^5}+O\left(\left(\frac{1}{x}\right)^7\right)$$
at $x= \infty$ this simplifies to:
$$\frac{\vartheta (x)}{\pi e} \approx \frac{x}{2 \pi e}\log\left(\frac{x}{2 \pi e}\right)-\frac{1}{8 e}$$
But that is as far as I can get.