There is no such operation, because (for any domain $D \subset \Bbb R^n$) the linear operator $$\nabla : C^{\infty}(D) \to \{\textrm{conservative vector fields on $D$}\}$$
is not injective: $$\ker \nabla = \{\textrm{constant functions $D \to \Bbb R$}\} \cong \Bbb R.$$ (For simplicity we take all objects to be smooth, i.e., infinitely differentiable, but this can be weakened considerably.)
On the other hand, recognizing $\nabla$ as linear tells us how to find the next-best thing to an inverse. Since $\nabla$ is surjective onto the set of conservative vector fields (this is just the definition of conservative), taking the quotient by the kernel induces an isomorphism
$$C^{\infty}(D) / \{\textrm{constant functions $D \to \Bbb R$}\} \stackrel{\cong}{\to} \{\textrm{conservative vector fields on $D$}\} .$$ By construction, we can interpret its inverse, $$\tilde\beta : \{\textrm{conservative vector fields}\} \stackrel{\cong}{\to} C^{\infty}(D) / \{\textrm{constant functions $D \to \Bbb R$}\} ,$$ as the map that assigns to any conservative vector field its $1$-parameter family of potentials.
Unwinding definitions shows that computing $\tilde\beta({\bf X})$ for a conservative vector field $${\bf X} = \pmatrix{P\\Q\\R}$$ (here taking $n = 3$ for notational simplicity) just amounts to integrating: A function $f$ is a representative of $\tilde\beta({\bf X})$ iff it satisfies $f_x = P$, $f_y = Q$, $f_z = R$, so
$$f(x, y, z) = \int_{\gamma} {\bf X} \cdot d{\bf s} = \int_{\gamma} P \,dx + Q \,dy + R \,dz,$$ where $\gamma$ is a path from some fixed reference point $(x_0, y_0, z_0) \in D$ to $(x, y, z)$. The conservativeness of $\bf X$ guarantees precisely that the above integrals are independent of the choice of path $\gamma$.
Let $\partial_l=\frac{\partial}{\partial x_l}$, $a=(a_1,a_2,a_3);r=(r_1,r_2,r_3)$
$$\nabla(\Phi)=\begin{pmatrix}\partial_1\cos(a_1r_1+a_2r_2+a_3r_3)\\
\partial_2\cos(a_1r_1+a_2r_2+a_3r_3)\\
\partial_3\cos(a_1r_1+a_2r_2+a_3r_3)\end{pmatrix}=\begin{pmatrix}
-a_1\sin(a\cdot r)\\
-a_2\sin(a\cdot r)\\
-a_3\sin(a\cdot r)\end{pmatrix}$$
Best Answer
Yes, it is possible. I will give an example to demonstrate the general procedure. Consider $f(x, y) = x^3 - 3xy^2 + x^2 + y^2 + \log x$. Then $\nabla f = \dfrac{\partial f}{\partial x}\hat{i} + \dfrac{\partial f}{\partial y}\hat{j} = \left(3x^2 - 3y^2 + 2x + \dfrac{1}{x}\right)\hat{i} + (2y -6xy)\hat{j}$.
To reverse this, we look at each component individually. We know that
$\dfrac{\partial f}{\partial x} = 3x^2 - 3y^2 + 2x + \dfrac{1}{x}\\ \dfrac{\partial f}{\partial y} = 2y - 6xy $
Therefore:
$\begin{align} \displaystyle f(x,y) & = \int \dfrac{\partial f}{\partial x}\, dx\\ & = \int 3x^2 - 3y^2 + 2x + \dfrac{1}{x}\, dx\\ & = x^3 - 3xy^2 + x^2 + \log x + u(y) \end{align}$
What is that $u(y)$? It's the "constant" of integration, of course.When we differentiate $f$ with respect to $x$ partially, any term of $f$ not containing $x$ is a constant - this includes terms containing only $y$.
Now to determine $u(y)$, we look at $\dfrac{\partial f}{\partial y}$. We could integrate this with respect to $y$, similar to what we did with $\dfrac{\partial f}{\partial x}$. Then some of the terms after the integration will be common. The terms that are not common are those that constitute $u(y)$. So we need to look for terms of $\dfrac{\partial f}{\partial y}$ that do not contain $x$, and integrate them. Here, $\dfrac{\partial f}{\partial y} = 2y - 6xy$, and the only term not containing $x$ is $2y$. Therefore:
$\displaystyle u(y) = \int 2y\, dy = y^2 + C$.
Thus, $\boxed{f(x, y) = x^3 - 3xy^2 + x^2 + \log x + y^2 + C}$.
In general: $$f(x, y) = \int \dfrac{\partial f}{\partial x} dx + \int \left[\text{terms of $\dfrac{\partial f}{\partial y}$ that do not contain $x$}\right]\, dy$$