In general, this will not be possible. For instance, the first non-trivial eigenfunction will have two nodal domains, so $f_1$ cannot be arbitrary. Furthermore, you expect to have some sort of gradient-type estimate which prevent low-frequency eigenfunctions from oscillating too wildly. As such, you will need to restrict the functions $f_k$ quite a bit.
However, even if you don't specify the eigenfunctions, the eigenvalues $\lambda_k$ will need to satisfy a version of Weyl's law, and thus their asymptotic behavior cannot be arbitrary. For instance, there is a result of Hormander which essentially shows that if the potential is smooth
$$N_\mu(\lambda)=(2 \pi)^{-n} \omega_n \operatorname{Vol}_g(M) \lambda^n+O\left(\lambda^{n-1}\right) . $$
where $N_\mu(\lambda)$ is the number of eigenvalues (counted with multiplicity) less than $\lambda$.
In general, it is possible to take a list of finitely many (non-repeating) numbers $\{\lambda_k\}_{k=1}^N$ and find a Riemannian metric where the bottom of the spectrum is exactly the values desired. There are some bounds on the possible degeneracy of eigenvalues, but if the eigenvalues are simple there is no issue. Here is a recent paper of Xiang He which gives a good overview of the current state of this problem. https://arxiv.org/pdf/2308.04190.pdf