The "simplest" classical explanation I know is the van der Waals interaction described by Keesom between two permanent dipoles. Let us consider two permanent dipoles $\vec{p}_1$ (located at $O_1$) and $\vec{p}_2$ located at $O_2$. Their potential energy of interaction is:
\begin{equation}
U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2}) = -\vec{p}_1\cdot \vec{E}_2 = -\vec{p}_2 \cdot \vec{E}_1
= -\frac{[3(\vec{p}_1\cdot\vec{u})(\vec{p}_2\cdot\vec{u})-\vec{p}_1\cdot \vec{p}_2]}{4\pi \epsilon_0 ||\vec{O_1 O_2}||^3}
\end{equation}
where $\vec{u}= \vec{O_1 O_2}/||\vec{O_1O_2}||$ and $\vec{O_1 O_2}$ is a vector going from $O_1$ to $O_2$.
Now, if the system is subject to thermal fluctuations then the dipoles can change their orientations according to a Boltzmann weight. If the time scale associated to the motion of the particles is much longer than the one associated to the dipole orientations then for each distance $r = ||\vec{O_1O_2}||$ on can trace over all the possible orientations.
The effective interaction between the dipoles is then characterized by the free energy of the system i.e. by:
\begin{equation}
\mathcal{F}(r|p_1,p_2) \equiv -k_B T \ln \int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2} \:e^{-\beta U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2})}
\end{equation}
where $d\Omega_i$ is the solid angle element associated with the direction of dipole $\vec{p_i}$ and $\beta=1/k_B T$.
This quantity is super difficult to compute. However, if the dipoles are sufficiently far apart, their interaction energy is small and one may expand the exponential in powers of $\beta U$:
\begin{eqnarray}
\mathcal{F}(r|p_1,p_2) \approx -k_B T \ln \int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2} && \times\\
&& \:\left(1 \right.\\
&& -\beta U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2}) \\
&& \left.+\frac{1}{2}\beta^2 U(\vec{p}_1,\vec{p}_2,\vec{O_1 O_2})^2\right)
\end{eqnarray}
The first integral $\int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2} 1 $ is trivial and equal to $1$. The second integral of the potential is zero by symmetry (each dipole is as likely to be at either a positive or negative angle from $\vec{u}$ and the integral of a scalar product of two independent vectors is zero). The only remaining term is of order 2 in $\beta U$ which decays with $1/r^6$. At the end of the day we end up with something like:
\begin{equation}
\mathcal{F}(r|p_1,p_2) \approx -k_B T \ln \left(1+\beta^2\frac{C}{r^6} \right) \approx -\frac{C}{k_B T}\frac{1}{r^6}
\end{equation}
where
\begin{equation}
C = \int \frac{d\Omega_1 d\Omega_2}{(4\pi)^2}\left( \frac{[3(\vec{p}_1\cdot\vec{u})(\vec{p}_2\cdot\vec{u})-\vec{p}_1\cdot \vec{p}_2]}{4\pi \epsilon_0}\right)^2
\end{equation}
is a number not important in terms of physical insight.
What matters here is that the average value of the dipole-dipole interaction is zero while the first non zero term not to be zero is owing to fluctuations. This is roughly the rule with van der Waals like interactions which is that they are fluctuation driven (the same is true in the quantum case).
Remark: the previous derivation assumes that the interactions travel at infinite speed which is not true. When the limitation of the speed of light is introduced then at large distances, the van der Waals interaction decays as $1/r^7$.
There is really just one underlying principle here: Any equation we write down should not depend on any arbitrary choices we made in order to define the quantities. All the examples you can discuss can be understood in this principle.
Can't add a vector and a scalar. Well, of course a vector is three numbers and a scalar is one number, so, for example, $\mathbf{v} + v$, , where $v$ is a speed, i.e. a scalar with units of velocity, doesn't even make mathematical sense. But we could add imagine adding one component of a vector to a scalar, i.e. $v_z + v$. But, this quantity shouldn't appear in a fundamental law of physics because our choice of what axis to call the $z$ axis is completely arbitrary, and if we made a different choice our equations would look different. But this is situation-dependent. For example, if we were discussing physics in a background uniform gravitational field, then we can use a convention where $z$ points along the direction of gravitational field. This is not arbitrary because the gravitational field sets a preferred direction. By declaring that we are going to call that particular direction the "z-direction'', it makes sense that any equations we then write down will only hold for that particular choice of $z$-axis. That is why the equation for the gravitational potential energy in a gravitational field, $U = mgz$, is valid even though $z$ is a component of the displacement vector $\textbf{r}$. However, you can still translate this equation into one that is valid for any choice of axes, namely $U = m \textbf{g} \cdot \textbf{r}$, where $\textbf{g}$ is the gravitational field vector.
Can't add numbers with different units. The point is that we normally work in units of physics which are chosen completely arbitrarily. If time $t$ is measured in seconds and position $x$ is measured in meters, then it
makes no sense to write down an equation involving $x+t$ because this equation would depend on our definition of "second" and "meter", and there is no reason why the laws of physics should depend on the second being defined to be 9,192,631,770 times the period of some radiation mode of a cesium atom. But, if we choose to work in natural units, then this is not an arbitrary choice because, as the name suggests, natural units are uniquely determined given fundamental constants of physics. In natural units, there is nothing wrong with writing an equation involving $x+t$, because we remember that we have made a special choice of units, and the equation will hold only in those units.
Of course, any equation that you can write in natural units can still be translated into arbitrary units. Take Einstein's famous mass-energy equivalence. In natural units ($c=1$) it states that $E = m$. Obviously, in arbitrary units, this is a bad equation because if $E$ is measured in Joules, and $m$ is measured in kg, then it would depend on the definitions of Joules and kg. But that's fine, because this equation only holds in natural units. Its translation into arbitrary units is $E = mc^2$, and the units now match up.
Can't add covariant 4-vectors and contravariant 4-vectors. Again, this is because in special relativity, in order to write down components of vectors we have to make an arbitrary choice of coordinate directions in space-time. Equations we write down in special relativity shouldn't depend on this choice, and this prevents us from adding covariant 4-vectors and contravariant 4-vectors because they transform differently when you change coordinate directions in space-time.
Can't add things in different vector spaces. This is just because, if $\mathbf{v}$ is in one vector space and $\mathbf{w}$ is in a totally different vector space, then you wrote an equation involving $\mathbf{v} + \mathbf{w}$ then it would depend on how you relate the bases between the two vector spaces, which -- given that they are totally different spaces -- there is no way to do non-arbitrarily.
Best Answer
I guess to obtain a potential "per one molecule" you should divide your potential by the Avogadro number 6x$10^{23}$/mol, you obtain the potential in (kilo)Joules, and after you differentiate it with respect to the coordinate you'll obtain the force in Newtons.