Given an action
$$S[\psi]=\int d^4x \mathcal{L}(\psi,\nabla_\mu \psi),$$
how do we find the equations of motion? The principle of stationary action says that the equations of motion can be found by extremizing the action, i.e., by setting its first-order change with respect to a change in its argument to zero. That is, we demand that
$$\left.\frac{dS[\psi+\alpha \delta \psi]}{d\alpha}\right|_{\alpha=0}=0.$$
Note that $\delta \psi$ need not be small, it is an arbitrary variation. The only thing we demand from $\delta\psi$ is that it must vanish at infinity. Let us look at the action you provided:
$$
S\left[\varphi, g_{\mu\nu}\right]=\int d^4x \sqrt{-g} e^{-2 \varphi}\left[R+4 g_{\mu \nu} \nabla^\mu \varphi \nabla^\nu \varphi\right]
$$
Set $$\left.\frac{d S\left[\varphi+\alpha \delta \varphi, g_{\mu \nu}\right]}{d \alpha}\right|_{\alpha=0}$$ equal to zero to get
\begin{equation}
\begin{aligned}
& 0=\frac{d}{d\alpha}|_{\alpha=0}\int d^4x \sqrt{-g} e^{-2 (\varphi+\alpha \delta \varphi)}\left[R+4 g_{\mu \nu} \nabla^\mu (\varphi+\alpha \delta \varphi) \nabla^\nu (\varphi+\alpha \delta \varphi)\right]\\
&=\int d^4 x \sqrt{-g} e^{-2 \varphi}\left\{\left[R+4 \nabla_\mu \varphi \nabla^\mu \varphi\right](-2 \delta \varphi)+8 \nabla^\mu \delta \varphi \nabla_\mu \varphi\right\} \\
&=\int d^4 x \sqrt{-g} e^{-2 \varphi}\left\{-2 R+8 \nabla^\mu \varphi \nabla_\mu \varphi-8 \square \varphi\right\} \delta \varphi \\
& \Rightarrow R-4 \nabla^\mu \varphi \nabla_\mu \varphi+4 \square \varphi=0
\end{aligned}
\qquad(1)
\end{equation}
where $g_{\mu\nu}\nabla^{\mu}\nabla^{\nu}=\Box$. I have made use of
$$
\begin{aligned}
\nabla^\mu\left(\delta \varphi \nabla_\mu \varphi e^{-2 \varphi}\right)= & e^{-2 \varphi} \nabla^\mu \delta \varphi \nabla_\mu \varphi+\Box \varphi \delta \varphi e^{-2 \varphi} -2 \nabla^\mu \varphi \nabla_\mu \varphi \delta \varphi e^{-2 \varphi} \\
\Rightarrow e^{-2 \varphi} \nabla^\mu \delta \varphi \nabla_\mu \varphi= & \delta \varphi e^{-2 \varphi}\left(2 \nabla^\mu \varphi \nabla_\mu \varphi-\Box \varphi\right)+\underbrace{\nabla^\mu\left(\delta \varphi \nabla_\mu \varphi e^{-2\varphi}\right)}_{\text{Integrates to an inconsequential boundary term}}
\end{aligned}
$$
Next, let us vary the action with respect to the metric tensor. As you probably know, we usually vary with respect to the inverse metric rather than the metric itself because it's usually easier to do it this way. Without further ado, let us write our action as
$$
S\left[\varphi, g^{\mu\nu}\right]=\int d^4x \sqrt{-g} e^{-2 \varphi}\left[R+4 g^{\mu \nu} \nabla_\mu \varphi \nabla_\nu \varphi\right]
$$
and set
$$\left.\frac{d S\left[\varphi, g^{\mu \nu}+\alpha \delta g^{\mu\nu}\right]}{d \alpha}\right|_{\alpha=0}$$ equal to zero:
$$
\begin{aligned}
& 0=\left.\frac{d}{d \alpha}\right|_{\alpha=0} \int d^4x \sqrt{-\operatorname{det}\left(g(\alpha))\right.}\left[R\left(g^{\mu \nu}+\alpha \delta g^{\mu\nu}\right)+4\left(g^{\mu \nu}+\alpha \delta g^{\mu\nu}\right) \nabla_\mu \varphi \nabla_\nu \varphi\right] e^{-2 \varphi} \\
& =\int d^4x\left[-\frac{1}{2} \sqrt{-g} g_{\mu \nu} \delta g^{\mu \nu}\left[R+4 g^{\alpha \beta} \nabla_\alpha \varphi \nabla_\beta \varphi \right]e^{-2 \varphi}+\sqrt{-g}\left[g^{\mu\nu}\frac{d R_{\mu\nu}\left(g^{\mu \nu}+\alpha \delta g^{\mu \nu}\right)}{d \alpha}\right|_{\alpha=0}\right. \\
& \left.+\delta g^{\mu\nu}R_{\mu\nu}+4 \delta g^{\mu \nu} \nabla_\mu \varphi \nabla_\nu \varphi\right] e^{-2 \varphi} \\
& =\int d^4 x \sqrt{-g} e^{-2 \varphi}\left[R_{\mu \nu}-\frac{1}{2} R g_{\mu \nu}+2 g_{\mu \nu} \nabla_\alpha \varphi \nabla^\alpha \varphi-2 \square \varphi g_{\mu \nu}+2 \nabla_\mu \nabla_\nu \varphi \right]\delta g^{\mu\nu} \\
& \Rightarrow R_{\mu \nu}-\frac{1}{2} R g_{\mu \nu}=-2 g_{\mu \nu} \nabla_\alpha \varphi \nabla^\alpha \varphi+2 \square \varphi g_{\mu \nu}-2 \nabla_\mu \nabla_\nu \varphi. \\
&
\end{aligned}
$$
In going from the first equality to the second, I used
$$\frac{d}{d\alpha}\sqrt{-\operatorname{det}(g(\alpha))}|_{\alpha=0}=-\frac{1}{2}g_{\mu \nu}\delta g^{\mu\nu}\sqrt{-\operatorname{det}(g)}, $$
and
$$\begin{aligned}\left.\int d^4x \sqrt{-g}g^{\mu\nu}\frac{d R_{\mu\nu}\left(g^{\mu \nu}+\alpha \delta g^{\mu \nu}\right)}{d \alpha}\right|_{\alpha=0}e^{-2\varphi}=
& \int d^4 x \sqrt{-g} e^{-2 \varphi} \nabla_\sigma\left[g_{\mu \nu} \nabla^\sigma \delta g^{\mu \nu}-\nabla_\lambda \delta g^{\sigma \lambda}\right] \\
= & -\int d^4 x \sqrt{-g} \nabla_\sigma e^{-2 \varphi}\left[g_{\mu \nu} \nabla^\sigma \delta g^{\mu \nu}-\nabla_\lambda \delta g^{\sigma\lambda}\right] \\
= & \int d^4 x \sqrt{-g}\left[g_{\mu \nu} \nabla^\sigma \nabla_\sigma e^{-2 \varphi}-\nabla_\mu \nabla_{\nu} e^{-2 \varphi}\right] \delta g^{\mu \nu} \\
= & \int d^4 x \sqrt{-g}\delta g^{\mu\nu}\left[g_{\mu \nu}\left[4 \nabla_\alpha\varphi \nabla^\alpha \varphi-2 \Box \varphi\right]-\left[4 \nabla_\nu \varphi \nabla_\mu \varphi-2 \nabla_\mu \nabla_\nu \varphi\right]\right]e^{-2\varphi}
\end{aligned}
$$.
The identity $$\left.\int d^4x \sqrt{-g}g^{\mu\nu}\frac{d R_{\mu\nu}\left(g^{\mu \nu}+\alpha \delta g^{\mu \nu}\right)}{d \alpha}\right|_{\alpha=0}e^{-2\varphi}=
\int d^4 x \sqrt{-g} e^{-2 \varphi} \nabla_\sigma\left[g_{\mu \nu} \nabla^\sigma \delta g^{\mu \nu}-\nabla_\lambda \delta g^{\sigma \lambda}\right]$$
is proven in every general relativity textbook. For example, you can take a look at chapter 4 of Carroll's book. Now,
$$
\begin{aligned}
& R_{\mu \nu}-\frac{1}{2} R g_{\mu \nu}=-2 g_{\mu \nu} \nabla_\alpha \varphi \nabla^\alpha \varphi+2 \square \varphi g_{\mu \nu}-2 \nabla_\mu \nabla_\nu \varphi \\
&\Rightarrow R=8 \nabla_\alpha \varphi \nabla^\alpha \varphi-6 \square \varphi
\end{aligned}
$$
By the way, if we plug this back into the final equation in (1), we find that
\begin{equation}
\nabla^{\mu}\varphi\nabla_{\mu}\varphi=\frac{1}{2}\Box\varphi
\end{equation}
which means that the gravitational field equation can now be written as
$$R_{\mu\nu}=-2\nabla_\mu\nabla_\nu \varphi.$$
To find the energy-momentum tensor, we simply write the equation of motion resulting from varying the action with respect to the metric in the form
$$R_{\mu \nu}-\frac{1}{2}Rg_{\mu\nu}=T_{\mu\nu}.$$
In our case,
$$T_{\mu\nu}=\Box\varphi g_{\mu\nu}-2\nabla_\mu\nabla_\nu \varphi.$$
