I've recently written a simple code to numerically compute the Hessian of some function (at a point). Most electronic structure packages will compute the Hessian and then project out the translations and rotations. I have tried to look up how one actually does this, but I can't find any good information really.

I feel like this should be pretty simple, so I feel a bit silly for not knowing how to do it, but if someone can provide some mathematical detail on how to project out the translational and rotational modes, that would be greatly appreciated.

  • 1
    $\begingroup$ Related: mattermodeling.stackexchange.com/questions/441/…. The link to the Gaussian website talks about the projection $\endgroup$
    – Tyberius
    Feb 25 '21 at 20:13
  • 1
    $\begingroup$ @ShoubhikRMaiti I think this actually answers the question, so I might try to implement it from there. There's a lot of extra stuff being described in this paper, so maybe someone will summarize in a nice way. If I figure it out first, I'll give an answer. $\endgroup$
    – jheindel
    Feb 25 '21 at 22:08
  • 1
    $\begingroup$ While converting the Hessian to internal coordinates is an interesting one, it is not how it is done in practice. In practice the translation and rotations are removed by going to the Eckart frame, which separates vibrations from rotations and translations. Those modes are then easily identifiable because their frequency is zero. I can post some python code to see how this is done if you would like $\endgroup$
    – Cody Aldaz
    Feb 26 '21 at 7:51
  • 3
    $\begingroup$ I think this gaussian page answers the question very clearly: web.archive.org/web/20191229092611/https://gaussian.com/vib Basically, you do transform the hessian to internal coordinates. I will try this and answer the question here, though the gaussian link explains it as clearly as one could hope to really. $\endgroup$
    – jheindel
    Feb 26 '21 at 9:05
  • 1
    $\begingroup$ @CodyAldaz They are constructing the transformation matrix from cartesian coordinates to internal coordinates. This transformation essentially projects out the external degrees of freedom. They are not diagonalizing the Cartesian hessian, but the hessian which is transformed to internal coordinates which forces exactly the six zero eigenvalues. So, it is true they are doing the displacements in cartesian coordinates (makes finite differencing easier), but they are not diagonalizing the hessian made up those derivatives. $\endgroup$
    – jheindel
    Feb 26 '21 at 22:10

I've put some code online https://gist.github.com/craldaz/b38e1c951d515c807c67aac303406343

that demonstrates how one removes translations and rotations from the Hessian using the description you linked https://web.archive.org/web/20191229092611/https://gaussian.com/vib/

First you have to form the axes of inertia which is done in the function eckart_frame, then the basis of vibrations is formed in the function vibrational_basis and finally the hessian is projected into this vibrational basis in normal_modes to remove all the translations and rotations.

The axes of inertia are used in the vibrational_basis function to form the matrix TR (or D in the web archive) which contains all translations and rotation. The null space of TR is the vibrational basis, B, in the Eckart frame. The final projection to get the Hessian in internal coordinates is done using B

$$H_{int} = BHB$$

*Note that I am using the SVD and not the Gram-Schmidt algorithm the Gaussian archive uses. It's much easier to do.

  • 1
    $\begingroup$ Nice code. Does it handle linear molecules ? It seems that you always remove 6 external degrees of freedom, unless I am mistaken. $\endgroup$
    – Hans Wurst
    Sep 10 '21 at 10:12
  • 1
    $\begingroup$ Cool! This seems like it'd be a great addition to cclib (e.g., as one of the calculation methods) since they already process Hessians. $\endgroup$ Sep 16 '21 at 20:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.