7
$\begingroup$

I have been working to implement optimisation in internal coordinates (for a transition state search software I am working on). However, the implementation details are almost always glossed over in the paper.

I have been trying to optimise with the Delocalised Internal Coordinates (Baker). In this case, "primitive" internal coordinates which include bonds, angles, dihedrals, out-of-plane bending etc. are built up. Then the B matrix (partial derivatives of primitives vs Cartesian components of the atoms) is formed, and then the G matrix ($B B^T$) is diagonalised and only the non-zero components are taken to obtain the non-redundant set of coordinates. For a non-linear molecule the number of non-redundant delocalised should be 3N-6. If the number is less than this, then it means not enought primitives were considered (e.g. there were two fragments not connected with any bond, or angles were not added).

This part is covered in papers. However, what is not covered in detail is how to deal with linear molecule. For vibrational analysis linear molecules have 3N-5 normal modes. However, I am not sure how many internal coordinates are required for optimising linear molecules. Because with a linear molecule like $\ce{CO2}$ for example, there are two degenerate linear bends. However, in terms of optimisation, it is obvious to me that only one bend mode is needed for optimisation, since as soon as the structure deviates from linearity, there is only one bending mode, the other one vanishes. I also know that for linear molecules, two orthogonal primitive linear bends are used by Gaussian optimisers for example.

So how many non-redundant coordinates should be generated for linear molecules? How to handle nearly linear molecules? Is there a threshold in terms of numbers that can be used for removing the redundant eigenvectors from G matrix for these cases?

$\endgroup$
4
  • 3
    $\begingroup$ +1. You might find this paper interesting: sciencedirect.com/science/article/abs/pii/S0022286020304129. Also papers about these internal coordinates often say that the algorithm is able to take advantage of point group symmetry, so CO2 should not require optimization of as many coordinates as COS for example. $\endgroup$ Oct 14, 2023 at 19:53
  • $\begingroup$ @NikeDattani Thanks! I will have a look. Although it would be possible to use symmetry, I don't want to try to code symmetry right now, and want to set up a general code that works reasonably for everything. $\endgroup$
    – S R Maiti
    Oct 14, 2023 at 20:00
  • 1
    $\begingroup$ This page has lots of useful information. I think you would indeed just have to set up a dummy atom for near-linear optimisation to break the degeneracies, switching from your normal optimisation to that mode when you detect the angle is close enough to 180 degrees. $\endgroup$ Oct 14, 2023 at 21:12
  • $\begingroup$ @ShernRenTee Yes, most of the papers on linear angles (Hoy et al) mention using two orthogonal bends with dummy atom, or a dummy axis. I wonder how many nonredundant degrees of freedom you would get with that method though if you have multiple linear bonds (carbon suboxide), or if the molecule is nearly linear. $\endgroup$
    – S R Maiti
    Oct 15, 2023 at 9:45

1 Answer 1

7
$\begingroup$

Indeed, Baker's delocalized coordinates routines generate 3N-6 (or 3N-5 for linear molecules) coordinates automagically. Note our small contribution where we distinguish between strong and weak coordinates.

For linear molecules ABC, you can indeed use symmetry (if the program has C(lin) symmetry), and if there are no other molecules present. Therefore, within QUILD we use dummy atoms (D). These are placed 1.0 Bohrs from the central atom B, with an angle of 90 degrees to one of the end point atoms (A). The angle DBC is then close to 90 degrees as well. All of the derivatives for angle ABC are then replaced by the dihedral ABDC, which is well defined.

Three small details: the choice of where to place the D atom is chosen within QUILD to point as close as possible to the center of mass if there are more molecules in the system. I.e. for D we have a bond of 1.0 to B, an angle of 90 degrees to A, but we need a dihedral angle to some other atom. And if only three atoms are present, we make an arbitrary choice for this reference atom.

Second detail: up to an deviation of (I think) 5 degrees from linearity in the starting structure, a dummy atom is added. During the optimization this may go beyond this value, and no dummy atom is any longer needed. Then the coordinates setup is reset, without dummy atom for this angle, and new delocalized coordinates are generated. Normally, the delocalized coordinates are only set once, and updated after e.g. 50 or 100 geometry steps.

Third detail: a dummy atom is only added if the middle atom has only two connections. If there are more atoms connected, the linear angle is either included with low weight (depending on deviation from linearity), such that the other "normal" angles are used, or the almost-linear angles are ignored.

$\endgroup$
8
  • $\begingroup$ Thank you very much for the detailed answer! Is the dummy atom moved during optimisation or is it stationary? Also, I understand that in linear cases there are two orthogonal bend directions, in the dummy atom case is there only one direction? Also, what is the threshold for the DLC generation (i.e. minimum eigenvalue for that is kept in U matrix) in almost linear cases where I imagine there would be some almost zero but not completely eigenvalues of the G matrix? $\endgroup$
    – S R Maiti
    Oct 16, 2023 at 8:46
  • $\begingroup$ At every geometry step it is repositioned, always in the same manner. In truly linear cases, yes, there are two bend directions, but by positioning the dummy atom, one specific one is chosen (and kept fixed). The threshold is dynamic. It is adapted automatically, such that 3N-6 coordinates are obtained (or 3N-5 for linear molecules). This is in particular needed for weakly bound systems. $\endgroup$
    – MSwart
    Oct 19, 2023 at 10:53
  • $\begingroup$ Also, I was wondering if there is detorioration of the hessian update if the dummy atom is repositioned? Because the coordinate system of that angle would effectively change every step. $\endgroup$
    – S R Maiti
    Oct 26, 2023 at 12:57
  • 2
    $\begingroup$ No, I don't see that happening. In particular for linear angles, where the same position is chosen every time. $\endgroup$
    – MSwart
    Oct 29, 2023 at 15:16
  • 1
    $\begingroup$ Probably I didn't explain myself well enough. The cos(theta) will suffer the same fate as a regular angle representation. For non-linear angles the displacement vectors are perpendicular to the bonds. With linear angles they are parallel to it (along the z-axis if all three atoms lie at the z-axis). So, you don't have a representation of the angle, only a "bond" between the end points. Adding a dummy atom corrects this. $\endgroup$
    – MSwart
    Dec 20, 2023 at 12:54

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .