# Transformation between two different sets of molecular vibrational normal coordinate systems

Lets assume we have $$N$$ Atoms and we treat them within the Born-Oppenheimer Approximation. We can calculate the adiabatic electronic groundstate potential. Lets assume we observe two local minima, denoted by $$\alpha$$ and $$\beta$$, on this potential energy surface. Lets denote these minima nuclear coordinates as $$R_\alpha$$ and $$R_{\beta}$$. At each minimum we can do an harmonic approximation of the potential, centered on the corresponding minimum, $$\begin{eqnarray} V_\alpha(R) = V(R_\alpha)+(R-R_\alpha)^T V'_{\alpha} + 1/2(R-R_\alpha)^TV''_{\alpha}(R-R_\alpha) \tag{1}\\ V_\beta(R) = V(R_\beta)+(R-R_\alpha)^T V'_{\beta} + 1/2(R-R_\beta)^TV''_{\beta}(R-R_\beta) \tag{2}\\ \end{eqnarray}$$

where the gradent is $$V'_\alpha=\nabla_RV(R)|_{R=R_\alpha}$$ and the Hessian $$V''_{\alpha,nm} = \frac{\partial^2 V}{\partial R_n \partial R_m}|_{R=R_\alpha}$$.

Since we are in a local minimum, we know that the gradients vanish,$$V'=0$$. For simplicities sake, we also assume that both potentials evaluated at the minimum strcucture are zero, $$V(R_{\alpha/\beta})=0$$.

The associated harmonic vibrational Hamiltonian in atomic units then takes on the form $$\tag{3} H_\alpha(R) = -\frac{1}{2}\nabla^T_RM^{-1}\nabla_R + 1/2(R-R_\alpha)^TV''_\alpha(R-R_\alpha)$$ where $$M$$ is a diagonal matrix with the nuclear masses.

Each harmonic potential allows the definition of vibrational normal modes and vibrational normal coordinates. The transformation to normal coordinates is defined as follows.

First we define the mass weighted displacement coordinates $$X_\alpha = M^{1/2}(R-R_\alpha)$$, in these coordinates we have $$\begin{eqnarray} H_\alpha(X_\alpha)= -\frac{1}{2}\nabla^T_{X_\alpha}\nabla_{X_\alpha} +1/2 X^T_\alpha M^{-1/2}V''_\alpha M^{-1/2} X_\alpha \tag{4}\\ H_\alpha(X_\alpha)= -\frac{1}{2}\nabla^T_{X_\alpha}\nabla_{X_\alpha} +1/2 X^T_\alpha W_\alpha X_\alpha \tag{5}\\ \end{eqnarray}$$ with the mass weighted Hesse matrix $$W=M^{-1/2}V''M^{-1/2}$$.

The last step to obtain normal coordinates is the transformation to coordinates that diagonalize the mass weighted Hesse matrix. Let $$WK = K\Lambda \tag{6}$$ where $$K$$ is the matrix with eigenvectors of $$W$$ as columns and $$\Lambda$$ a diagonal matrix with corresponding eigenvalues. The required transformation is then $$X= KQ\tag{7}$$ which leads to $$H_\alpha(Q_\alpha)= -\frac{1}{2}\nabla^T_{Q_\alpha}\nabla_{Q_\alpha} +1/2 Q^T_\alpha \Lambda_\alpha Q_\alpha \tag{8}$$

Summarized we have two different sets of transformations, one for each minimum, $$\begin{eqnarray} Q_\alpha = K^{-1}_\alpha M^{1/2}(R-R_\alpha)\tag{9}\\ Q_\beta = K^{-1}_\beta M^{1/2}(R-R_\beta) \tag{10} \end{eqnarray}$$

After this rather lengthy setup comes my question. How can I transform from one set of normal coordinates $$Q_\alpha$$ to the other $$Q_\beta$$ and how do I interpret this. Does that mean that I can express any molecules vibrational coordinates in the vibrational coordinates of another molecule(made from the same number and type of atoms)?

Intuitively I see problems when the geometries are translated and rotated with respect to each other assuming a global cartesian coordinate system for both molecules. There should be no simple internal vibrational coordinate that could connect two such geometries but on the other hand, internal vibrational coordinates should form a complete coordinate system and allow me to express any structure with the same amount of atoms in it. I think I am missing something very obvious to reconcile this "contradiction", which is why I need some input.

• +1. Please do take a look at the edits I made, especially the change of 6 different instances of i to I. An academic journal wouldn't let you publish a paper with i not capitalized, and we hope to keep high standards here at MMSE. Commented May 22, 2021 at 7:39
• @NikeDattani Thank you for editing. I had the comfort of working with native English speaking coauthors in the past, who took care of the punctuation and grammar. But that is of course no reason to be sloppy. Commented May 22, 2021 at 8:42
• I think the answer to your question is yes. You just need to find the Jacobian that transforms between the two coordinates systems. And it is definitely possible to express the vibrations of one molecule in terms of the eigenvectors of another as long they span the same space. In all likelihood, the eigenvectors of one molecule are a bad basis for the eigenvectors of another molecule. I could imagine some interesting similarity analysis of configurational isomers being done in this way. Commented May 22, 2021 at 16:52
• @jheindel Equations 9 and 10 are linear with regard to the displacement vectors $\Delta R_\alpha=R-R_\alpha$ and $\Delta R_\beta = R-R_\beta$. But I think the problem is that it looks more like a affine transformation with respect to the molecular coordinates vector $R$, that connects both equations. I'm lacking experience with regard to affine transformation and the conversion between them. I think that complicates the problem but I'm not certain, hence the question. Commented May 22, 2021 at 17:35

I agree with Emil Zak's answer that normal mode coordinates from different minima should, in general, be put together with special care. However, relating the two different sets of normal mode origin from different minima (or transition states) is not of nonsense, see SVP model of Hua Guo and Bin Jiang and related theories.

To make things short, I will use the same symbols as in the question. The relation between the normal mode coordinates and Cartesian coordinates, are $$\mathbf{R}=\mathfrak{R}_{\alpha}\cdot\left(\mathbf{R}_{\alpha}+\sum_{i=7}^{3N}{{Q}_{\alpha}}_i\cdot {\mathbf{K}_{\alpha}}_i\right)+\mathbf{t},\tag{1}$$ and $$\tag{2} \mathbf{R}=\mathfrak{R}_{\beta}\cdot\left(\mathbf{R}_{\beta}+\sum_{i=7}^{3N}{{Q}_{\beta}}_i\cdot {\mathbf{K}_{\beta}}_i\right)+\mathbf{t}.$$ Here, we assume that the $$\mathbf{R}_{\alpha}$$ are originated from the COM of molecule. $$\mathfrak{R}_{\alpha}$$ mean a rotation matrix and $$\mathbf{t}$$ is the translation. And the mass-weighted part are omitted for simplicity.

Note that I separate rotation, translation, and vibration on purpose. Although in practice you can keep them together and compute the all $$3N$$ DOF based on normal modes fashion, like in Gaussian software you get 6 "low frequencies" for translations and rotations, which should be exact zeros in principle, when you use the normal modes vectors to represent the coordinate, it is not a good idea to use the rotation normal modes, more importantly, the molecular orientation must be the same as of the reference geometry, because the vibration normal modes are attached to the reference geometry, and that is why the rotation matrix is outside the big parentheses.

Using the formulae above, it is easy to convert the normal mode coordinate $$\mathbf{Q}_{\alpha}$$ to $$\mathbf{R}$$, and if we have an algorithm to get $$\mathbf{R}$$ to $$\mathbf{Q}_{\beta}$$, the problem will be solved.

The main idea is

1. Translate the molecule so that the COM located at $$(0,0,0)$$.
2. Rotate $$\mathbf{R}_{\beta}$$ to align with $$\mathbf{R}$$ with some algorithm (Wahba's problem). In my implement, Kabsch's algorithm is used. So the rotation matrix $$\mathfrak{R}_{\beta}$$ is get. (Edit: Here, one may minimize the mass-weighted RMSD $$\Delta=\sqrt{\frac{1}{N}\sum_{i=1}^{N}m_i|\mathfrak{R}_{\beta}\mathbf{R}_{\beta i}-\mathbf{R}_{\alpha i}|^2}, \tag{3}$$ where $$m_i$$ is the mass of the $$i$$-th atom. Kudin and Dymarsky has shown that minimization of mass-weighted RMSD is equivalent with the Eckart condition. I thank @Hans Wurst provided the information.)
3. Utilizing the orthogonality of mass-weighted normal mode vectors, get the coordinates $$\mathbf{Q}_{\beta}$$. Note that in this step you only need the vibration ones, and exclude the rotation and translations
4. Now rotate $$\mathbf{R}_{\beta}$$ to align with $$\mathbf{R}+{{Q}_{\beta}}_i\cdot {\mathbf{K}_{\beta}}_i$$, update $$\mathfrak{R}_{\beta}$$ and $$\mathbf{Q}_{\beta}$$.
5. Do it iteratively until the two molecules are perfectly aligned.

The algorithm above is implemented with Julia programming language by me as shown in gist.

Then you can do $$\tag{4} \mathbf{Q}_{\alpha} \leftrightarrow \mathbf{R} \leftrightarrow \mathbf{Q}_{\beta}$$

Problem solved.

I would suggest that the reader refer to some published work, thanks to the anonymous reviewers.

1. GÃ¡bor CzakÃ³, J. Phys. Chem. A 116, 116, 7467â€“7473 (2012)
2. Philip R. Bunker, Molecular Symmetry and Spectroscopy (1979), Academic Press, New York (see Ch. 7 therein)

There are some methods that use intramolecular curvilinear coordinates as the "bridge", although I think use Cartesian coordinates is easier to be understood.

• Why is the RMSD minimizing rotation the correct one ? I could also do a a rotation that aligns the tensor of inertia axes but these rotation are in general not identical and would then yield different values when doing the transformation. Is there a physical reason that the RMSD minimizing rotation is the only sensible rotation ? This is the only remaining point that is somewhat unclear to me, everything else is fine. My first (apparently wrong) intuition was that the transformation to normal coordinates should be independent of such a rotation but that clearly isn't the case. Commented Jul 1, 2021 at 12:25
• @Hans Wurst , Actually no special reason. I think you can used the principle axis of inertia tensor. What is most tricky, I think, is that I run the procedure of alignment-vibrate iteratively and finally the RMSD will be zero. Since the (3N-6) real vibrational NM are used to represent the intramolecular DOF and the rotation is represented by a rotational matrix, and by using the iterative procedure so the rotation-vibration DOF are self-consist. Commented Jul 1, 2021 at 13:24
• The reason why the RMSD approach works, is given here doi.org/10.1063/1.1929739 The RMSD minimizing rotation is the same as the one that transforms to the Eckart frame, defined by the Eckart conditions. These axes are the ones that should be used if one wants a local coordinate system in which rotation and vibration is (almost) decoupled. There should be no need to do it iteratively, you only have to make sure that you use the set of axes that is "continuous", as there can be mutliple solutions. Commented Aug 4, 2021 at 7:43
• Thank you, Hans. I did not know that work. What I missed in my answer is that I should use a mass-weighted RMSD approach, instead of the one I used (no-weighted) here. Commented Aug 4, 2021 at 9:00
• Since I really wrote a paper and submitted it. I am in a mess now. :-( Commented Aug 4, 2021 at 9:06

The paradox discovered by @HansWurst comes by the assumption that it is possible to treat both minima on the 1D PES separately and define a separate set of normal coordinates for them. 3N normal coordinates are functions of 3N Cartesian coordinates of the bodies in our system. It turns out that for a double-well potential in 1D, which you describe here, the normal-mode approximation is a very bad one. We only have one dynamical dimension, hence we can only define one normal coordinate for this system. Splitting into two independent coordinates is the 'mistake' here. To treat this problem we can define a variational basis composed of sums of 1D Harmonic Oscillator basis functions centred at the two local minima. Such a basis is likely to serve well.

In general, when more dimensions are involved it is rather customary to perform rotations on full sets of normal coordinates (so called Duschinsky rotation), such that, for instance, a set of normal coordinates for the electronic ground state can be expressed with normal coordinates of a displaced and deformed electronically excited state.

• Do you know a detailed reference to the Duschinsky rotation ? I find mention of it but haven't found a reference with a clear definition of it and all steps involved. For example when the ground state structure and the excited state structure are unaligned and quite different. Assume you have two optimized structures as xyz files with completely different orientations and origin for example and you want to transform from one set of normal coordinates to the other. Commented Jun 30, 2021 at 10:32
• @HansWurst, you may want to look into PhD thesis of Joonsuk Huh "Unified description of vibronic transitions with coherent states" (2010), sec. 2.2. Commented Jun 30, 2021 at 16:43
• Thank you, that is exactly what I was looking for. Commented Jun 30, 2021 at 21:01

Intuitively I would say that yes of course you can relate the two coordinate systems. But whereas in one of the two, the corresponding eigenvalue-matrix is diagonal {zero values at off-diagonal elements), if you express the other system in it, the eigenvalue-matrix will most likely not be diagonal.

Moreover, I think it might only work if you use the Cartesian equivalent of the normal mode vectors, as is often done for transformations.