By a straight application of Pick's lemma, $|f'(0)|\le1-e^{-2}$. However, the fact that $f(z)\ne0$ throws a wrench in the search for a best upper bound. A lower bound on the supremum can be observed in the solution $\varphi(z)=\exp\big(\!\frac{z-1}{z+1}\!\big)$, which satisfies $\varphi(0)=e^{-1}$, $\varphi(1)=1$, and has an essential singularity at $z=-1$, and $\varphi'(0)=2/e$, so this proves that $\sup_f|f'(0)|\ge2/e$. There is still a gap between $2/e$ and $1-e^{-2}$, though, so there is still work to be done.
However, we can use the same trick used by Twiceler below, and consider the function $F(z)$ defined by $f=\varphi\circ F$. Since $f(z)\ne0$, it is the exponential of an analytic function, and the Mobius transform is invertible, so such an $F$ must exist. (A little extra detail in this step: By the Weierstrass factor theorem, if $f$ is nonzero, then $f(z)=e^{g(z)}$ for some $g$. Thus, if we define $F(z):=\frac{1+g(z)}{1-g(z)}$, then $f(z)=\exp\big(\!\frac{F(z)-1}{F(z)+1}\!\big)=\varphi(F(z))$.) Then
$$\begin{align}
f(0)=\exp\bigg(\!\frac{F(0)-1}{F(0)+1}\!\bigg)=e^{-1}&\Rightarrow F(0):=c_n=\frac{i\pi n}{1-i\pi n}\Rightarrow |c_n|^2=1-\big(1+(n\pi)^2\big)^{-1} \\
&\Rightarrow \varphi'(c_n)=\frac2e(1-i\pi n)^2\Rightarrow |\varphi'(c_n)|=\frac2e\big(1+(n\pi)^2\big)
\end{align}$$
and so by Pick's lemma, $|f'(0)|=|\varphi'(c_n)F'(0)|\le|\varphi'(c_n)|\cdot(1-|c_n|^2)=2/e$. (This proves that $\sup_f|f'(0)|=2/e$, and $\varphi$ saturates the bound.) Moreover, if the bound is saturated, then $F$ is a Mobius transformation that preserves the unit circle, i.e. $F(z)=\theta\dfrac{c_n+\theta z}{\theta-\overline{c_n}z}$ for unimodular $\theta$, so
$$f(z)=\exp\!\bigg[\!\!\frac{c_n+\theta z-(1-\overline{c_n}\theta z)}{c_n+\theta z+(1-\overline{c_n}\theta z)}\!\!\bigg]
=\exp\!\bigg[\!\!\frac{(1-in\pi)\theta z-(2n^2\pi^2-in\pi+1)}{(1+in\pi)+(2n^2\pi^2+in\pi+1)\theta z}\!\!\bigg],$$
where $|\theta|=1$ and $n\in\mathbb Z$.
Some explanatory remarks: The moral of the story is that if we have to avoid zero, we can transition from the problem of finding a nonzero function $f:\mathbb D\to\mathbb D-\{0\}$ to a function $g:\mathbb D\to\{x\in\mathbb C:\operatorname{Re}x<0\}$ whose range is the left half-plane (this is so that $|e^g|<1$), and via a Mobius transform similar to the Cayley transform, we can map that range back to $\mathbb D$ to get a function $F:\mathbb D\to\mathbb D$ (note that there are no longer any holes in the range), to which we can apply Twiceler's analysis using Pick's lemma.
Although I have defined $F$ via $f=\varphi\circ F$, I have been careful not to say $F:=\varphi^{-1}\circ f$, because $\varphi$ is not injective, and indeed takes on every point in $\mathbb D-\{0\}$ infinitely many times (as you approach $z=-1$ from different directions). In particular, the value $\varphi(z)=e^{-1}$ happens at an infinite number of points, indexed as $\{c_n\}_{-\infty}^\infty$, so $f(0)=e^{-1}\Rightarrow g(0)=-1+2i\pi n\Rightarrow F(0)=\frac{i\pi n}{1-i\pi n}:=c_n$. But after this, Pick's lemma applies as normal on $F$, so $|F'(0)|\le1-|c_n|^2=1/(1+n^2\pi^2)\Rightarrow |f'(0)|\le2/e$.
The remaining question is how to characterize the solutions that saturate the inequality. The algebra is not fun, so I have just shown the solution here, but the point is that all solutions are in the form $\varphi\circ F$ where $F$ is a Mobius transformation which preserves the unit circle and satisfies $F(0)=c_n$, so that $f(0)=e^{-1}$ is satisfied.
