Regarding the usage of 'classical potentials' in quantum mechanics I am familiar with basic quantum mechanics and I know that there is no concept of 'force' in quantum mechanics, unlike in classical mechanics. Problems in quantum mechanics are solved by writing down the Hamiltonian for a system, and trying to solve for the various eigenvalues. Some of the first problems that are taught to students learning Quantum Mechanics are the harmonic oscillator problem, and the Hydrogen atom problem, where the Hamiltonian takes the same form as a classical system. Since moving over to quantum mechanics requires one to lose several ideas that have been built up while learning classical mechanics, why is the potential found in quantum mechanics problem of the same form as classical mechanics?
The potential for the hydrogen atom, for example, is classical in origin and is derived from the Coulomb force. How is this direct usage of the potential, which is purely classical in origin, justified by the theory?  
 A: The Hamiltonian (and thus in particular the form of the potential) is not "justified by the theory". The form of the Hamiltonian is an input to the theory. However, there are various ways to see that using the classical Hamiltonian in its quantized version is really the right thing to do:


*

*The "classical" limit of large quantum numbers gives the classical behaviour, as demanded by the correspondence principle.

*Even without looking at the classical theory, the predicted spectral lines (for the hydrogen atom) with this Hamiltonian are empirically correct, as seen in experiments. This is the only "justification" the Hamiltonian/Lagrangian needs: The justification for the classical form is also, in the end, that it is the form that produces the correct predictions.

*The form of the potential itself can be derived from the underlying quantum field theory, quantum electrodynamics, yet there you give the Lagrangian of classical electrodynamics as an input, shifting the question to "Why does the Lagrangian have this form?". Answer, again: Because it predicts the correct results. (Although one might argue that there is a certain "naturality" to this minimally coupled Yang-Mills theory that is not as evident for the Coulomb potential.)
A: To justify this in theory, you must start from quantum electrodynamics which treats the electromagnetic field quantum mechanically. So, there is no classical potential at the fundamental level (although, in practice, one does make use the notion of classical external fields a lot in quantum field theory, but this notion is in principle flawed as a fundamental quantity).
In quantum electrodynamics there is only a coupling between the fields associated with charged particles and the field associated with photons. So, there is no Coulomb potential at this level, the classical Coulomb potential can be derived from first principles as an effective classical potential, and this comes with small corrections due to vacuum polarization which are important at small distances.
The use of classical potentials when justified dynamically may still lead to subtle artifacts, because quantum mechanics is fundamentally so different from classical mechanics. E.g. Lev Vaidman has recently suggested that the notion of non-locality in the Aharonov- Bohm effect may be such an artifact.
A: Well the usage of the potential is not purely classical: in quantum mechanics we promote $\hat x$ to an observable and then reinterpret $V$ as acting on square matrices via its Taylor expansion about some point. We can interpret these things by using Ehrenfest's theorem prescriptively, to obtain the rough equivalence. 
Some definitions will be helpful: an "observable" here is a Hermitian linear operator; the word Hermitian (equivalently self-adjoint) means the closest equivalent to "its matrix-transpose is equal to its term-by-term complex conjugate." The term anti-Hermitian is the same but inserts a minus sign; the transpose is equal to its negative term-by-term complex conjugate. We can represent this compactly by speaking of an operator's "dual" $\hat A^\dagger$ which is the conjugate-transpose of an operator; then Hermitian means $\hat A = \hat A^\dagger$ while anti-Hermitian means $\hat A = -\hat A^\dagger$. We care about Hermitian matrices because their eigenvalues are real numbers. So we draw the analogy between these eigenvalues and the quantized numbers that we actually observe, and to make these real our "observables" must be Hermitian operators.
Now suppose we start with the assumptions that we've got a vector space using the field $\mathbb C$ with observable quantities somehow represented by self-adjoint linear operators on that space, whose averages in a state $|\psi\rangle$ are equal to $\langle \psi| \hat A |\psi\rangle$ and whose evolution is determined by $|\partial_t \psi\rangle = \hat E |\psi\rangle$ for some $\psi$ and some special evolution operator $\hat E$. We want to incorporate Ehrenfest's theorem: this new mathematics must on average reproduce the results of classical mechanics. Let's step back a second and talk about how we can do that.
What are we emulating?
The target formulation is the Hamiltonian formulation, where there is a Hamiltonian function $H$ such that the evolution of any quantity over time is given by a Poisson bracket $$\frac{{\rm d}A}{{\rm d}t} = \frac{\partial A}{\partial t} - \{H, A\}.$$ The conditions which such a Poisson algebra must satisfy are a little complex; first you need to have an associative algebra on these quantities using the field of real numbers $\mathbb R$, and then you need to steal its product $\cdot$ and define your Poisson bracket $\{,\}$ as a Lie bracket (so $\{x,y\} = -\{y, x\};$ furthermore $x$, $y$, and $\{x, y\}$ have to all be the same kind of thing). This bracket must then obey two properties with respect to that product: 


*

*$\{x, y\cdot z\} = \{x, y\}\cdot z + y\cdot\{x,z\},$ 

*$\{x, \{y, z\}\} + \{z, \{x, y\}\} + \{y, \{z, x\}\} = 0.$


That's a big mouthful; but if you can guarantee these things then there's going to be some canonical coordinates $\{q_i, q_j\} = \{p_i,p_j\} = 0$ with $\{q_i, p_j\} = \delta_{ij}$ within which that $H$ describes a total energy function for of a bunch of particles with positions $q_i$ and momenta $p_i,$ and the above evolution equation then describes how all of the physical quantities change according to Newton's laws acting within the system.
How do we do it?
Noticing the adjoint equation $\langle \partial_t \psi| = \langle \psi|\hat E^\dagger$ we can derive that the time-rate-of-change of any average quantity must be $$\frac{\rm d}{{\rm d}t} \big(\langle \psi |\hat A|\psi \rangle\big) = \langle \psi | \left(\frac{\partial \hat A}{\partial t} + \hat E^\dagger~\hat A + \hat A \hat E  \right)|\psi\rangle.$$ This looks suggestively like the above in how it starts with a partial time derivative, we see that there must be some connection then between $\hat E$ and $H$. As a nice freebie, the combination $\hat E^\dagger~\hat A + \hat A \hat E$ must also be Hermitian when $\hat A$ is, so this evolution is clearly preserving our real eigenvalues in a nice way. Clearly what this tells us, if we pursue this procedure, is that our implementation of $\{a, b\}$ is going to be "before you've taken the average values, you did something with the operators." So $\{a, b\} = \langle \psi | \{\hat a, \hat b\} |\psi\rangle.$
So our first requirement is that $\{a, b\}$ be a Lie bracket. That's actually harder than it sounds. The above equation naively suggests that our underlying operator-juggling should be, say, $-\hat a \hat b - \hat b^\dagger \hat a$ or so, but for Hermitian operators this is a symmetric bracket, not an antisymmetric one. In fact to be antisymmetric we need $\{a, b\} \propto \hat a \hat b - \hat b \hat a$ to provide that minus sign. Furthermore if we consider that bracket in isolation, called the "commutator," it actually maps two Hermitian operators to an anti-Hermitian operator, which means that our Lie algebra can't be restricted then to Hermitian operators! However, we can fix this: if you multiply the resulting anti-Hermitian operator by a pure-imaginary number, then you get back a Hermitian operator without destroying the antisymmetry. Let's call this pure-imaginary number $1/(i\lambda)$ for some real $\lambda$. We see that our underlying symbol-shuffling must be a bracket looking like $\{\hat a, \hat b\} = \frac{1}{i\lambda} (\hat a \hat b - \hat b \hat a).$
We can force the above to fit this mold by specifying that $\hat E = \frac{\hat H}{i \lambda},$ which remember because $\hat H$ is an observable hence Hermitian, causes $\hat E$ to be an anti-Hermitian operator, creating that pivotal $-$ sign when we write $\hat E^\dagger = -\hat E.$ The resulting equation then is simply:
$$\frac{\rm d}{{\rm d}t} \big(\langle \psi |\hat A|\psi \rangle\big) = \langle \psi | \left(\frac{\partial \hat A}{\partial t} - \{\hat H, \hat A\} \right)|\psi\rangle.$$
Follow-up maths.
You can also double-check the other two equations we have above and confirm that they hold if $\cdot$ is simply taken as the operator-composition of two Hermitian operators; however again this does not map (Hermitian, Hermitian) to Hermitian. So, the usual prescription is to use $a \cdot b \leftrightarrow (\hat a \hat b + \hat b \hat a)/2,$ to ensure this in an appropriate way. The Poisson algebra and the quantum algebra perfectly match on average if we use the modified commutator above for the Poisson bracket and the symmetric-product for the classical-variable-product.
This means that the $\lambda$ describes something which lives completely outside of these "on average" predictions and is particular to the quantum formalism; there are now a couple of ways to discover that in fact $\lambda = \hbar$. The first is simply to connect the Planck-Einsten equation for the energy levels of photons, $E = h f$, with the wavefunctions that we're obtaining: the equation $i\lambda |\partial_t \psi\rangle = \hat H|\psi\rangle$ is solved for an energy-eigenfunction $\hat H |\eta\rangle = \eta |\eta\rangle$ as $i \lambda |\partial_t \eta\rangle = \eta |\eta \rangle,$ which we can solve as $|\eta\rangle = e^{-it\eta/\lambda}|\eta_0\rangle,$ interpreting this as $e^{-i\omega t}$ gives $\eta = h f = h (\omega / 2 \pi) = h \eta / (2 \pi \lambda),$ ultimately showing $\lambda = h / (2 \pi).$ Or you can use the fact that $\{\hat x, \hat p\} = 1$ for the canonical coordinates implies $\hat x \hat p - \hat p \hat x = i \lambda$ which turns out to give a Heisenberg uncertainty principle of size $\lambda$, which you can scale precisely to $\hbar$. You can do a hybrid between these, embedding into the latter the De Broglie relation that $p = \hbar k,$ since the above commutation relation can be realized by pretending $\hat p = -i \lambda \partial_x.$ Applying this to some wave $e^{i(kx - \omega t)}$ we get $p = \lambda k$ and $\lambda = \hbar$ straight-away.
Perhaps the most honest way to phrase these is to say that $\lambda$ is the symbol $\hbar$ from the very start, and then look to experiments to pin down the values of $\lambda = \hbar$ in these particular cases which, on average, replicate classical mechanics. What I've just described is therefore ways to measure $\hbar,$ rather than ways to discover that $\lambda$ is $\hbar$.
A: Part of your trouble is the incorrect belief that QM doesn't contain force. Let's review how QM works. First, we take classical equations in however many variables we may have, including derivatives of variables with respect to others. We then replace every variable, except time and the variable in a derivative's "denominator", with an operator on the Hilbert space whose eigenvalues are the empirically possible values of the original classical variable. Then we replace each $f = 0$ equation with $\left\langle f\right\rangle = 0$, and we may be able to rearrange somewhat. For example, Newton's second law $\mathbf{F} = \frac{d}{dt}\mathbf{p}$ becomes $$\left\langle \mathbf{F}\right\rangle = \left\langle \frac{d}{dt}\mathbf{p}\right\rangle = \frac{d}{dt}\left\langle \mathbf{p}\right\rangle.$$ Or, if you prefer the potential-dependent form for conservative forces, we go from $-\frac{\partial U}{\partial\mathbf{x}} = \frac{d}{dt}\mathbf{p}$ to $$\left\langle-\frac{\partial U}{\partial\mathbf{x}}\right\rangle = \frac{d}{dt}\left\langle \mathbf{p}\right\rangle.$$ Suppose you knew classical mechanics inside out, but no QM. In that case, you'd know that Newtonian mechanics admits several equivalent formulations that rely less on force in their Hamilton-Jacobi equation. You'd know also that, although these approaches are equivalent, one may be more useful than another in some contexts. Each of them can be quantised; for example, the forces of Newtonian mechanics can be handled as described above. But QM typically used the Hamiltonian formalism because we're interested in energy eigenstates. The reason this is useful is because (i) when we write a state as a superposition of such eigenstates the coefficients have a complex exponential time-dependence and (ii) a quantity's mean is conserved iff its operator commutes with the quantised Hamiltonian (this is analogous to a Poisson brackets result). Similarly, in QFT it's worth starting from the field-theoretic Lagrangian formalism so we can identify the physical symmetries of a system.
