Short answer: no, of course not. Perpetual motion machines do not exist. There are deep theoretical reasons for that, as well as less deep but nevertheless convincing experimental data. That's why we care so much ../about renewable energy and stuff: we know that energy supply is limited in any real physical system.
From the perspective of mechanics (classical, relativistic, or quantum) when you take some energy from a system, it jumps to a state with lower energy. If there were a system with no lower bound on energy levels, we could withdraw energy forever, forcing the system to fall lower and lower on the energy ladder (this is roughly what forced Dirac to introduce the concept of the Dirac sea). Thus, a sufficiently good mechanical system should have a lower bound on possible energy states.
We now focus on atoms and molecules. What would bound their energy levels from below? To begin, we need a good framework to describe these energy levels. Unfortunately, classical mechanics fails here: the virial theorem implies that for a potential of the form $V(r) \sim \frac{1}{r}$ the time average of kinetic and potential energies, denoted $\langle T\rangle$ and $\langle V\rangle$, respectively, are related by $$2\langle T \rangle = -\langle V \rangle,$$ and the total energy $E$ (which is conserved) is given by $$E = \langle E \rangle = \langle T \rangle + \langle V \rangle = \frac{1}{2} \langle V \rangle.$$
For a single electron (which has charge $-e$) moving in a Coulomb field of a nucleus with charge $Ze$, the potential energy is $V(r) = -\frac{kZe^2}{r}$, where $k$ is a Couloumb constant. Thus, if the distance $r$ to the nucleus tends to zero, the potential energy tends to $-\infty$, and so does the total energy $E$. This, together with the fact that an electron orbiting a nucleus would lose energy in the form of electromagnetic waves, was one of the major problems of classical mechanics that needed a "quantum treatment".
Let's dive into the (non-relativistic) quantum world. Unfortunately, we are immediately forced to accept an approximation; as of this writing, I'm unaware of a solution that doesn't use one. The Born-Oppenheimer approximation suggests that since the protons and neutrons are ../about 2000 times heavier than electrons, they move much slower and can be considered static. Thus, we must now prove that a system of $N$ electrons moving in a field of $K$ fixed nuclei has a lower bound on possible values of energy.
Under this framework, a molecule is a carefully chosen separable Hilbert space of possible states together with a "nice" essentially self-adjoint operator acting on it. The Hilbert space we are will use is denoted $(L^2(\mathbb{R}^3))^{\otimes N} \cong L^2(\mathbb{R}^{3N})$, where $\cong$ denotes an isomorphism. The notation indicates that we can represent an element of our Hilbert space in two equivalent ways; we'll use the second, which is a square-integrable function of $3N$ real variables which are coordinates of the electrons - 3 coordinates for each of the $N$ electrons. The inner product of two functions $\psi$ and $\phi$ is defined as an integral over the whole $3N$-dimensional space:
$$\langle\psi,\phi\rangle = \int\limits_{\mathbb{R}^{3N}} \overline\psi(r_1\dots r_N) \phi(r_1 \dots r_N) dr_1 \dots dr_N$$
Remark: The actual space should be tensored by $\mathbb C^{2^N}$ to account for spin, but the effect of spin on energy is too small for us to care in this case. Nevertheless, it exists. Furthermore, the elements of the space are actually not functions;they are equivalence classes of functions modulo function having zero norm. Following the usual colloquial terminology, we'll call them functions nevertheless; we'll refer to them as wavefunctions.
A physical state is not the same as a wavefunction: you can multiply the wavefunction by an arbitrary non-zero complex constant, and the result represents the same state. States form some peculiar kind of space; we shall stick to wavefunctions, which are elements of a vector space so we can use the full power of linearity. Thus, it is common to assume wavefunctions to be normalized such that[note]Here we use <,> to denote inner product.[/note] $\langle \psi,\psi\rangle=1$. We shall assume this too.
Next, we need to define the Hamiltonian. Here, the Hamiltonian is similar to that in non-quantum mechanics, with the exception that we must quantize it: turn numbers into operators. The general solution to this problem is unknown. (See Groenewold's theorem on the inconsistency of canonical quantization and John Baez's article series on some modern approaches). In our case, however, the recipe is simple: turn the kinetic energy $\frac{p_i^2}{2m}$ into $-\frac{\hbar^2}{2m} \Delta_i$ and turn the potential energy $V(r_1,r_2,\dots,r_N)$ into an operator that multiplies a wavefunction by $V$.
The energy consists of
The usual quantum-mechanical way to write this is: $$H = -\sum\limits_i \frac{\hbar^2}{2m} \Delta_i - \sum\limits_{i,A} \frac{kZ_Ae^2}{|r_i - R_A|} + \sum\limits_{i,j} \frac{ke^2}{|r_i - r_j|}$$ We now detail each item in the above equation
Nobody likes constants, right? Let's stick to atomic units instead, where $m=k=e=\hbar=1$, and the speed of light is $137$ (this is exactly $\frac{1}{\alpha}$, where $\alpha$ is the fine-structure constant). The Hamiltonian becomes a bit cleaner: $$H = -\sum\limits_i \frac{1}{2} \Delta_i-\sum\limits_{i,A} \frac{Z_A}{|r_i - R_A|} + \sum\limits_{i,j} \frac{1}{|r_i - r_j|}$$
What does it mean? An operator should act on wavefunctions. For a function $\psi$, the operator takes some second derivatives of $\psi$, multiply $\psi$ by other functions, and adds all this up. The result is another function - the result of operator acting on $\psi$, called $H\psi$. One important thing is that this operator is essentially self-adjoint, though we shall treat it as if it's adjoint: this follows from the Laplacian $\Delta_i$ being self-adjoint, and multiplication byreal-valued function being a self-adjoint operator as well.
The energy (well, the average or expectation value energy), of a state $\psi$ is $\langle \psi, H\psi\rangle$ (the inner product of $\psi$ and $H\psi$), which happens to be equal to $\langle H\psi , \psi\rangle$, thanks to self-adjointness. Physicists love to use bra-ket notation for expressions like this, and write it as $\langle \psi | H | \psi\rangle$.
We have a space of functions, a differential operator acting on these functions, and an integral that we wish to bound from below whenever the functions are normalized. I'll leave a huge number of unpleasant technical issues under the hood - we are going to presume that all equations make sense, that no domain issues occur, etc. Look into "Hall, Quantum Mechanics for Mathematicians" for all the details.
Our way to prove that the expression $\langle \psi, H\psi\rangle$ is indeed bounded below will be broken up into three parts.
A good thing ../about the average energy is that it is linear in the energy operator. This means that, for a sum of two operators, we have $$\langle \psi, (H_1+H_2)\psi \rangle = \langle \psi, H_1\psi \rangle + \langle \psi, H_2\psi \rangle,$ which means that we can analyze the terms in the Hamiltonian separately.
Now, the last term is electron-electron repulsion. Let's write it down again: $$\sum\limits_{i,j} \frac{1}{|r_i - r_j|}$$ The only thing that matters for our purposes is that this operator is multiplication by a nonnegative function, which we'll call $f(r_1,\dots,r_N)=\sum\limits_{i,j} \frac{1}{|r_i - r_j|}$ for now. Non-negativeness implies that the expectation value of this term is never below zero: $$\langle \psi, f \cdot \psi \rangle = \int\limits_{\mathbb R^{3N}} |\psi|^2 f dr_1 \dots dr_N \geq 0$$
Thus, since the expectation value is linear, we can simply drop the electron-electron repulsion term: this action can only lower the energy; since we want to find a lower bound on energy levels, this is okay.
Now the simplified Hamiltonian looks like this: $$H' = -\sum\limits_i \frac{1}{2} \Delta_i - \sum\limits_{i,A} \frac{Z_A}{|r_i - R_A|}$$ The electrons do not interact anymore, which makes the total Hamiltonian a sum of single-particle Hamiltonians of individual electrons: $$H' = \sum\limits_i \left[ -\frac{1}{2} \Delta_i - \sum\limits_{A} \frac{Z_A}{|r_i - R_A|} \right] = \sum\limits_i H_i,$$
where each of $H_i$ depends only on the $i$-th electron coordinates. In the tensor-product-version of the description of our Hilbert space of states, this means that $H_i$ acts on the $i$-th component of the tensor product $(L^2(\mathbb R^3))^{\otimes N}$. Whatever Hilbert space description you use, this implies that the whole system Hamiltonian $H'$ is separable, so analyzing the expectation values of $H_i$ is, in essence, a single-particle problem: the wavefunction $\psi(r_1, \dots, r_N)$ does still depend on all electrons' coordinates, but they are ignored by $H_i$. So, if we find a lower bound for each of $H_i$, we are done.
Everything before this section was actually pretty straightforward; it is now that we'll do a simple trick. Look at $H_i$ again: $$H_i = -\frac{1}{2} \Delta_i - \sum\limits_{A} \frac{Z_A}{|r_i - R_A|}$$ It would be awesome if we could put the Laplacian under the sum, for the resulting summands would be just the Hamiltonians of a hydrogen-like atom, --- a quantum system with any charged particle orbiting around another charged particle. This is one of the most important quantum systems and one of the few that has an exact and known analytical solution.
So, let's actually put the Laplacian under the sum!
$$H_i = \sum\limits_{A} \left[ -\frac{1}{2M} \Delta_i - \frac{Z_A}{|r_i - R_A|} \right] = \sum\limits_A H_{i,A},$$where $M$ is the number of nuclei. Now, each term of the sum is the Hamiltonian of an atom-like system, with a nucleus of charge $Z_A$ and an electron orbiting it, except that the electron has mass $M$ instead of $1$. Thankfully, this is not a problem: plug $m_e=M$ into the exact solutions of the hydrogen-like atom, and we get solutions for each $H_{i,A}$ term of the sum above. We don't need the actual solutions, though; all we need is the fact that expectation value of energy is bounded below in this solution.
So, what we've done is shown $$H \geq H' = \sum\limits_i H_i = \sum\limits_{i,A} H_{i,A},$$ and the expectation value of each $H_{i,A}$ is bounded below, therefore the expectation value of $H$ itself is bounded below.
This whole idea was inspired by the paper of Tosio Kato "Fundamental properties of Hamiltonian operators of Shrödinger type". You'll find many more details on the topic there.
All this is just a tiny bit of the whole story of why the world around us doesn't explode/freeze/break apart. For a full exposition, check out "The stability of matter: from atoms to stars" by Elliott H. Lieb. I should thank Valentin Fadeev for introducing me to this awesome work.