# Computing the bond spring constant from first principles

Given a set of atoms and their positions in space, is there a way to compute the bond spring constant for each pair of atoms that are considered to be bonded (I can already determine the bonds with OpenBabel)? I don't want to rely on empirical data, and it's okay if the spring constant is wrong by up to about 20%. I already know the bond length and bond energy, if that helps.

For background: my ultimate goal is to compute the normal modes of a molecule, by finding the generalized eigenvalues of an undamped mass-spring system. I already know the mass matrix, but I don't know the stiffness matrix.

• I guess one thing I could do is use OpenBabel to find the bond energy, then perturb the location of one atom in a bond, find the bond energy again, and divide the change in bond energy by the distance the atom was perturbed. – taktoa Feb 6 at 7:41

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.

• yes, I just put "first principals" in the title to distinguish from the myriad of questions about computing the bond spring constant from IR spectroscopy data and the like – taktoa Feb 7 at 3:01
• So if I understand correctly there is no transformation that will allow you to convert a (diagonal) matrix of force constants for a vector of normal modes into a (diagonal) matrix representing displacements along connected atoms. So the answer to the actual question should be "no" and then the second question (requesting accuracy bounds) would also be "no"? – Buck Thorn Feb 7 at 15:40
• @BuckThorn - it depends on the coordinate system for your Hessian. As I said above, if you evaluate the Hessian in internal coordinates, then you'll have force constants for the bonds in your z-matrix. – Geoff Hutchison Feb 7 at 19:11
• @GeoffHutchison Maybe I am mixing topics, but my understanding is that the weighted Hessian has cross-terms, it is not usually diagonal in any but normal coordinates (except maybe for a simple molecule where they coincide). Question becomes to what extent you can ignore those cross-terms if you transform into internal coordinates (a transformation which you can do starting from Cartesian if you assume small enough displacement ie linearize). I realize force fields often ignore cross terms, but the question is how big the error is. Perhaps this is totally unimportant. – Buck Thorn Feb 7 at 19:29
• The only way to know would be to do it for a molecule of interest. Indeed, for larger molecules, the normal modes typically involve multiple atoms and thus the internal coordinates won't help much - but my point is that the Cartesian coordinates aren't useful. – Geoff Hutchison Feb 7 at 20:46