By the [Euler–Maclaurin formula][1] (with $p=4$, $m=0$, and $g(x):=\frac1n\,f^2(\frac xn)$ in place of $f(x)$ there in the formula),
$$d_n(f):=\int_0^1 f^2(y) \, dy - \frac{1}{n}\sum_{i=1}^n f^2(i/n) \\
=-\frac1{2n}\,f^2(1)-\frac1{6n^3}\,f(1)f'(1)+O(n^{-4})\le O(n^{-4})$$
for each $f\in\mathcal F$. 

---

However, one cannot get a constant factor $O(n^{-4})$ good for all $f$ at once. Specifically,  
$$\sup_{f\in\mathcal F}d_n(f)\ge d_n(f_n)
=\frac{1}{\pi ^4 (2 n+1)^2}\sim\frac1{4\pi^4 n^2}$$
as $n\to\infty$, where 
$$f_n(x):=\frac1{\pi^2}\,
\Big(\sin (\pi  x)-\frac{\sin (\pi  (2 n+1) x)}{(2 n+1)^2}\Big).$$

(Note that $f_n\notin\mathcal F$, since $f_n'(0)\ne0$ and $f_n'(1)\ne0$. However, $f_n$ can be approximated however closely by functions in $\mathcal F$ with respect to the norm given by the formula 
$$\|f\|^2=\max_{[0,1]}(f^2)+\int_0^1(f'')^2.$$ 
More generally, this approximation shows that $\sup_{f\in\mathcal F}d_n(f)$ will not change if the conditions $f'(0)=0=f'(1)$ are removed from the definition of $\mathcal F$.) 

  [1]: https://en.wikipedia.org/wiki/Euler%E2%80%93Maclaurin_formula#The_formula