tldr, you want the force constant from the Hessian matrix .. ideally in internal coordinates
First off, you're saying 'first principles,' which most computational chemists would take to be some level of quantum chemical method, whether density functional or wavefunction-based. None of that is available in Open Babel - it only implements a very small number of classical force fields.
Even if you have optimal bond lengths and energies, that doesn't directly give you the force constants. For that, you want the second derivative.
In principle, you could start perturbing atoms and evaluating the energy. In that way you could build up a stochastic approximation of the PES and evaluate numerical second derivatives. I don't recommend that.
A more accurate way would be at least to get the forces / gradients at each of those points, which you can use to build up the Hessian matrix. There are some notes about that as part of the psi4numpy notebooks which can be used to drive the Psi4 quantum program.
The problem is that for any molecule over a few atoms, the Cartesian Hessian matrix won't directly give you a bond force constant. Most vibrations include multiple atoms, and are rarely along any particular Cartesian axis.
The psi4numpy notebooks give some indication of how you can evaluate the Hessian in internal coordinates, which should be closer to what you want.