Is there a natural (suitable) definition for functional derivative in Curved space time If $$\delta S = \int \sqrt g F[\phi] \delta \phi\tag{1}$$
Then is it natural to define the functional derivative as follows,
$$\frac{\delta S}{\delta \phi} = F[\phi].\tag{2}$$
In particular does this definition satisfy the commutativity of the functional derivatives.
I understand that this is not the standard definition of a functional derivative, but if I define this way this makes a certain calculation I am doing much easier to control. So to Clarify what I want to know is that if I define the 'functional derivative' this way then if
$$S = \int \sqrt{g} L[\phi, g_{\mu \nu}]\tag{3}$$ then 
is the following true?
$$\frac{\delta}{\delta\phi} \frac{\delta}{\delta g_{\mu \nu}} S = \frac{\delta}{\delta g_{\mu \nu}} \frac{\delta}{\delta\phi} S \tag{4} $$
This is my procedure for computing the second functional derivative suppose
$$\frac{\delta}{\delta \phi}S = E[\phi, g_{\mu \nu}]\tag{5}$$ 
and 
$$\frac{\delta}{\delta g_{\mu \nu}}S = E^{\mu \nu}[\phi, g_{\mu \nu}]\tag{6}$$ 
(equations of motion) Then to compute second functional derivative we write, eg.
$$E[\phi, g_{\mu \nu}](x) = \int \sqrt{g} d^4 y E(y) \hat \delta(x-y)\tag{7}$$
where $$\hat \delta (x-y) = \frac{\delta(x-y)}{\sqrt{g}}\tag{8} $$ is the generalised delta function. After this we can use the same definition to compute
$$\frac{\delta} {\delta g_{\mu \nu}} \int \sqrt{g} d^4 y E(y) \hat \delta(x-y)\tag{9} .$$ Of course now the result would involve delta functions. 
 A: Comments to the question (v4):


*

*First a disclaimer. Note that even for a smooth local functional, the existence of a functional/variational derivatives is not guaranteed, but depends on appropriately chosen boundary conditions.

*Calculus of variations, functional/variational derivatives, Frechet derivatives,  Gateaux derivatives, etc, is a huge mathematical topic. If OP is looking for rigor, then we suggest to post the question on Math.SE or Mathoverflow.SE. In this answer we will just make some heuristic comments and ignore boundary terms. 

*Let there be given an $n$-dimensional pseudo-Riemannian manifold $(M,g)$.  

*We can consider a Dirac delta distribution
$$\tag{A} \delta_M(x,y)~=~\frac{\delta^n(x-y)}{\sqrt{|g(x)|}}, \qquad g(x)~:=~\det(g_{\mu\nu}(x)),$$
on the manifold, cf. OP's eq. (8). Here $x$ and $y$ denote two points in a coordinate neighborhood $U\subseteq M$, [and to simplify notation, they also denote the coordinates of the two points]. 

*The Dirac delta distribution (A) is independent of the choice of coordinates, i.e. transforms as a scalar.

*Next we would like to define the functional/variational derivative wrt. fields $\phi^{\alpha}(x)$ [which could include the metric field $g_{\mu\nu}(x)$].

*The functional derivative $\frac{\delta_M S}{\delta\phi^{\alpha}(x)}$ 
should satisfy
$$\tag{B}  \delta S~=~ \int \! d\mu(x)~\frac{\delta_M S}{\delta\phi^{\alpha}(x)}\delta\phi^{\alpha}(x), \qquad d\mu(x)~:=~\sqrt{|g(x)|}d^nx,$$ 
for appropriate infinitesimal variations $\delta\phi^{\alpha}(x)$, cf. OP's eq. (1).

*In particular, we have the rule
$$\tag{C} \frac{\delta_M \phi^{\beta}(y)}{\delta\phi^{\alpha}(x)}
~=~\delta^{\beta}_{\alpha}~\delta_M(x,y). $$ 
When $\phi^{\alpha}$ is the metric field $g_{\mu\nu}$, note this subtlety.

*When $\phi^{\alpha}$ and/or $\phi^{\beta}$ are the metric field $g_{\mu\nu}$, then commutativity of functional derivatives are generally not expected, since the manifold notion (B) of functional differentiation depends on the metric, cf OP's eq. (4). In fact, it is easy to find counterexamples.

*For the rest of this answer we will assume that $\phi^{\alpha}$ and $\phi^{\beta}$ are not the metric field $g_{\mu\nu}$. It is then expected that the functional derivatives commute
$$\tag{D} [\frac{\delta_M}{\delta\phi^{\alpha}}, \frac{\delta_M}{\delta\phi^{\beta}}]~=~0. $$ 
See also e.g. Ref. 1.

*Example. If the functional is of the form
$$\tag{E}S~=~\int \! d\mu(x)~L(\phi(x),\partial\phi(x),x),$$ 
and we define 
$$\tag{F}\nabla^{(x)}_{\mu}
~:=~ \frac{1}{\sqrt{|g(x)|}}d^{(x)}_{\mu}\sqrt{|g(x)|} ,\qquad 
d^{(x)}_{\mu}~:=~\frac{d}{dx^{\mu}},$$
then the 1st functional derivative reads
$$\frac{\delta_M S}{\delta\phi^{\alpha}(x)}
~=~\frac{\partial L(x)}{\partial \phi^{\alpha}(x)}
-\nabla^{(x)}_{\mu}\frac{\partial L(x)}{\partial \phi^{\alpha}_{,\mu}(x)}$$
$$\tag{G}~=~\int \! d\mu(y)\left[\frac{\partial L(y)}{\partial \phi^{\alpha}(y)}\delta_M(x,y)
+ \frac{\partial L(y)}{\partial \phi^{\alpha}_{,\mu}(y)}d^{(y)}_{\mu}\delta_M(x,y)\right],$$
so that the 2nd functional derivative becomes symmetric
$$\frac{ \delta^2_M S}{\delta\phi^{\beta}(y)\delta\phi^{\alpha}(x)}
~=~\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y) 
-\nabla^{(y)}_{\nu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y)\right]$$
$$+\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}_{,\mu}(y)}d^{(y)}_{\mu}\delta_M(x,y)
-\nabla^{(y)}_{\nu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}_{,\mu}(y)}d^{(y)}_{\mu}\delta_M(x,y)\right]$$
$$~=~\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y) 
-\nabla^{(y)}_{\nu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}(y)}\delta_M(x,y)\right]$$
$$-\nabla^{(x)}_{\mu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}(y)\partial \phi^{\alpha}_{,\mu}(y)}\delta_M(x,y)\right]
-\nabla^{(y)}_{\nu}\nabla^{(x)}_{\mu}\left[\frac{\partial^2 L(y)}{\partial\phi^{\beta}_{,\nu}(y)\partial \phi^{\alpha}_{,\mu}(y)}\delta_M(x,y)\right]$$
$$\tag{H}~=~\left[(x,\alpha) \leftrightarrow (y,\beta) \right].$$
Here we have used that
$$\tag{I} [f(x)-f(y)]\delta_M(x,y)~=~0, \qquad [\nabla^{(x)}_{\mu}+d^{(y)}_{\mu}]\delta_M(x,y)~=~0.$$
References:


*

*Bryce DeWitt, Supermanifolds, Cambridge Univ. Press, 1992; eq. (5.1.3).

A: When I'm confused about all the functional derivatives, I'm doing that:
$$ S(\phi_1,\phi_2,\phi_3,g_1,g_2,g_3) =\sqrt{g_1}\,L(\phi_1,g_1) + \sqrt{g_2}\,L(\phi_2,g_2) + \sqrt{g_3}\,L(\phi_3,g_3)$$
I hope you'll recognize your $\int\sqrt{g}\, L[\phi,g]$ in the expression above.
The functional derivatives then become more familiar gradients:
$$\frac{\partial S}{\partial \phi_i} = \pmatrix{
\sqrt{g_1}\,\frac{dL}{d\phi}\left(\phi_1,g_1\right)\\
\sqrt{g_2}\,\frac{dL}{d\phi}\left(\phi_2,g_2\right)\\
\sqrt{g_2}\,\frac{dL}{d\phi}\left(\phi_3,g_3\right)\\
}\qquad
\frac{\partial S}{\partial g_i} = \pmatrix{
\sqrt{g_1}\,\frac{dL}{dg}\left(\phi_1,g_1\right)\\
\sqrt{g_2}\,\frac{dL}{dg}\left(\phi_2,g_2\right)\\
\sqrt{g_2}\,\frac{dL}{dg}\left(\phi_3,g_3\right)\\
}+\pmatrix{
\frac{d\sqrt{\phantom{l}g_1}}{dg_1}\,L\left(\phi_1,g_1\right)\\
\frac{d\sqrt{\phantom{l}g_2}}{dg_2}\,L\left(\phi_2,g_2\right)\\
\frac{d\sqrt{\phantom{l}g_3}}{dg_3}\,L\left(\phi_3,g_3\right)\\
}$$
The mixed derivative will become a diagonal matrix:
$$\frac{\partial^2 S}{\partial \phi_i\partial g_j} = \pmatrix{
    \sqrt{g_1}\,\left.\frac{d^2L}{d\phi dg}\right|_{\phi_1g_1} + \frac{d\sqrt{\phantom{l}g_1}}{dg_1}\,\left.\frac{dL}{d\phi}\right|_{\phi_1g_1}&0&0\\
0&  \sqrt{g_2}\,\left.\frac{d^2L}{d\phi dg}\right|_{\phi_2g_2} + \frac{d\sqrt{\phantom{l}g_2}}{dg_2}\,\left.\frac{dL}{d\phi}\right|_{\phi_2g_2}&0\\
0&0&...\\
}$$
From that point of view your "convenience definition" is (a) incomplete, (b) doesnt really look sensible and (c) doesn't look too convenient after all. 
I mean -- are you going to divide the gradients componentwise (b) by the $\sqrt{g_i}$?  Okay, but what about (a) the extra terms in the derivatives over $g$?  Do you want to subtract the extra terms, or leave them? In any case that "new definitions" will definitely not commute (c).
I suggest that you make up your mind in those simple terms, and then generalize from 3D up to the functional spaces...
A: Assuming that g is constant, the only way your equations "hold," is if $F[\phi] = e^\phi$. 
I hope this is helpful.
