Note: I have realized that, using the Stable Manifold Theorem, one can prove the smoothness of the solution $\phi$ that I describe below. Thus, I am modifying my answer to incorporate that.
Local existence and uniqueness of a smooth solution near $p$ satisfying $\phi(p)=0$ and $\phi\ge0$ near $p$ follows from an application of the Stable Manifold Theorem. Here is the argument.
This is a local question, so we might as well assume that $M=\mathbb{R}^n$, that $p=0$, that $g = g_{ij}dx^idx^j$ satisfies $g_{ij}(0)=\delta_{ij}$, and that the function $V$ has a Taylor expansion $V = h_{ij}x^ix^j + O(|x|^3)$, where $(h_{ij})$ is a symmetric, positive definite matrix. We are looking for a closed $1$-form $d\phi = f_i\ dx^i$ near $x=0$ so that $\phi$ satisfies the equation $g^{ij}f_if_j = V$ and, at the same time, satisfies $\phi(0)=0$ and $\phi\ge0$ near $x=0$.
Let $p_i$ be the coordinates on $T^*\mathbb{R}^n$ such that the canonical $1$-form has the expression $p_i\ dx^i$. Then the graph of $d\phi$, described by equations $p_i = f_i(x)$, will be a Lagrangian submanifold for the $2$-form $dp_i\wedge dx^i$ and will lie in the zero locus of the Hamiltonian $H(x,p) = g^{ij}(x)p_ip_j - V(x)$. Therefore, it will be a union of integral curves of the Hamiltonian vector field
$$
X_H = 2g^{ij}p_i\frac{\partial\ \ }{\partial x^j}
+ \left(\frac{\partial V}{\partial x^k}
- \frac{\partial g^{ij}}{\partial x^k}p_ip_j\right)\frac{\partial\ \ }{\partial p_k}.
$$
This graph will have to pass through the unique singular point of $X_H$, i.e., $x = p = 0$ (since $\phi$ clearly must have a critical point at $x=0$ because it vanishes there and is nonnegative nearby), and the linear part of $X_H$ at $x=p=0$ is
$$
Y_H = 2\delta^{ij}p_i\frac{\partial\ \ }{\partial x^j}
+ 2h_{ij}x^i\frac{\partial\ \ }{\partial p_j}.
$$
The unstable manifold of $Y_H$ is the $n$-dimensional submanifold defined by $p_i = L_{ij} x^j$, where $L$ is the (unique) symmetric positive definite square root of $(h_{ij})$. The stable manifold of $Y_H$ is defined by $p_i = -L_{ij}x^j$. It follows from the Stable Manifold Theorem that $X_H$ has a smooth $n$-dimensional unstable submanifold $N_+$ given by $p_i = f_i(x) = L_{ij}x^j + O(|x|^2)$ and a smooth $n$-dimensional stable submanifold $N_-$ given by $p_i = -L_{ij}x^j + O(|x|^2)$.
From the dynamics of $X_H$, it is clear that the only possibility for $d\phi$, when $\phi$ satisfies the above conditions and is at least $C^2$, is to have its graph be $N_+$. Conversely, taking $d\phi = f_i(x)\ dx^i$ where $p_i = f_i(x)$ is the (necessarily Lagrangian) unstable manifold of $X_H$ and fixing the additive constant by requiring that $\phi(0)=0$ does give a smooth solution to the original equation.
