$M$ is a geodesically convex Riemannian manifold, that is, for any two points $p,q$ on $M$, there is a unique minimizing geodesic connecting them. Can we conclude that for any $p \in M$, the function $f(x)=(d(x,p))^2$ is smooth where $d$ is the distance function ?
Yes. As noted on Wikipedia
"The distance function from $p$ is a smooth function except at the point $p$ itself and the cut locus."
The cut locus is the set of all points $q \in M$ such that there is more than one distinct minimizing geodesic between $q$ and $p$. Therefore, $d$ will be smooth at all points other than $p$, and hence $d^2$ will be smooth on all of $M$.