While smoothness is not yet clear to me, the local existence and uniqueness of a $C^2$ solution near $p$ satisfying $\phi(p)=0$ and $\phi\ge0$ near $p$ is not hard. Here is the argument (and an indication of what will need to be done to prove smoothness, if it does actually follow).
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 that $X_H$ has a $C^1$, $n$-dimensional unstable submanifold $N_+$ given by $p_i = f_i(x) = L_{ij}x^j + O(|x|^2)$ and a $C^1$, $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)$ 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$ evidently does give a $C^2$ solution to the original equation.
Whether this $\phi$ has higher regularity is completely determined by whether the unstable manifold of $X_H$ is more regular at $x=p=0$ than $C^1$. I haven't had time to think about this, so I don't have an opinion at the moment.