The traditional way of getting the right geometry for your system is by taking a fixed atomic structure, perform an electronic structure minimization, calculate the forces acting on the atoms from that, use them to alter the atomic structure according to the forces and then use the altered structure to start over again. This is repeated until you converge to the equilibrium geometry.
However there is an alternative (quite elegant) way for DFT calculations using pseudopotentials or the PAW method with plane-wave basis sets: the ab-initio molecular dynamics method introduced by Roberto Car and Michele Parrinello. The "killer feature" of the Car-Parrinello method is that it allows the simultaneous solution of the electronic structure problem and the equations of motion for the nuclei which allows for the relaxation of the nuclei to find stable structures, i.e. the equilibrium geometry, as well as for thermal simulations of solids and liquids - so this method does offer much more than a simple geometry optimization. How does it achieve that? In the Car-Parrinello approach, the total Kohn-Sham (KS) energy is the potential energy as a function of the positions of the nuclei. The forces from this energy then determine the molecular dynamics (MD) for the nuclei. The special feature of the Car-Parrinello algorithm is that it also solves the quantum electronic problem using MD. This is accomplished by deriving the equations of motion from a fictitious Lagrangian which contains a fictitious kinetic energy of the electronic (KS) wave functions and some constaints to ensure that the wave functions remain orthonormal:
\begin{align}
\mathcal{L} &= \sum_{n} f_{n} m_{\psi} \langle \dot{\psi}_{n} | \dot{\psi}_{n} \rangle + \frac{1}{2} \sum_{i} M_{i} \dot{\vec{R}}_{i} \! {}^{2} + E_{\mathrm{DFT}}[\psi_{n}, \vec{R}_{i}] - \sum_{m, n} \Lambda_{n,m} \left( \langle \psi_{m} | \psi_{n} \rangle - \delta_{n,m} \right)
\end{align}
where $f_{n}$ is the occupation of the $n^{\mathrm{th}}$ Kohn-Sham eigenstate $\psi_{n}$, $M_{i}$ and $\vec{R}_{i}$ are the mass and the position of the $i^{\mathrm{th}}$ nucleus, $\delta_{n,m} = \begin{cases}1 & n=m \\ 0 & n \neq m \end{cases}$ is the Kronecker delta, $\dot{\psi}_{n} = \frac{\mathrm{d} \psi_{n}}{\mathrm{d}t}$ and $\dot{\vec{R}}_{i} = \frac{\mathrm{d} \vec{R}_{i}}{\mathrm{d}t}$ are the first derivatives with respect to time and $\Lambda_{n,m}$ is a Lagrange multiplier ensuring the constraint of wave function orthonormality.
The first term in this Lagrangian represents the fictitious kinetic energy of the wave functions. There is no physical correspondence for this term. Ideally one would choose the fictitious mass $m_{\psi}$ of the wave functions equal to zero. This introduction of this unphysical quantity also lead to naming the Lagrangian fictitious.
The second term describes the classical kinetic energy of the nuclei.
The third term $E_{\mathrm{DFT}}[\psi_{n}, \vec{R}_{i}]$ is the density functional total energy, which is a functional of the electronic wave functions and the atomic positions.
The last term is the constraint of orthonormal wave functions $| \psi_{n} \rangle$.
Applying the principle of least action to this Lagrangian leads to the following Euler-Lagrange equations (equations of motion for the ab initio MD simulation):
\begin{align}
M_{i} \ddot{\vec{R}}_{i} &= \underbrace{- \nabla_{\vec{R}_{i}} E_{\mathrm{DFT}}}_{= \, \vec{F}_{i}} \\
m_{\psi} | \ddot{\psi}_{n} \rangle &= - \underbrace{ \frac{1}{f_{n}} \frac{\delta E_{\mathrm{DFT}}}{\delta \langle \psi_{n} |} }_{= \, \hat{H} | \psi_{n} \rangle } + \underbrace{ \sum_{m} | \psi_{m} \rangle \frac{\Lambda_{n,m}}{f_{n}} }_{\tilde \, | \psi_{n} \rangle \epsilon_{n}}
\end{align}
where $\ddot{\psi}_{n} = \frac{\mathrm{d}^{2} \psi_{n}}{\mathrm{d}t^{2}}$ and $\ddot{\vec{R}}_{i} = \frac{\mathrm{d}^{2} \vec{R}_{i}}{\mathrm{d}t^{2}}$ are the second derivatives with respect to time, $\hat{H}$ is the Kohn-Sham hamiltonian, $\epsilon_{n}$ is $n^{\mathrm{th}}$ Kohn-Sham eigenvalue and $\vec{F}_{i}$ is the force acting on the $i^{\mathrm{th}}$ nucleus.
The first set of equations is just Newton's equations of motion for the nuclei moving under the forces derived from $E_{\mathrm{DFT}}$.
The stationary solution of the second set of equations is equivalent to the Kohn-Sham equations, since for a steady state all time derivatives vanish so that:
\begin{align}
\hat{H} | \psi_{n} \rangle &= \sum_{m} | \psi_{m} \rangle \frac{\Lambda_{n,m}}{f_{n}}
\end{align}
Constructing a matrix $\mathbf{\Lambda}$ from the Lagrangian multipliers $\Lambda_{n,m}$ shows that $\mathbf{\Lambda}$ is proportional to the transpose of the matrix representation of $\hat{H}$, i.e. $\Lambda_{n,m} \propto H_{m,n}$. Diagonalizing $\mathbf{\Lambda}$ leads to the eigenvalues $\epsilon_{n}$ of the Kohn-Sham equations.
So, this Lagrangian does not lead to the time-dependent Schroedinger equation for the electronic wave functions. Rather it creates a dynamic that (in the ideal case) keeps the electrons always in their ground state.
