Determination of the Eckart molecule

Determination of the Eckart molecule-fixed frame by use of the apparatus of quaternion
algebra
Sergey V. Krasnoshchekov, Elena V. Isayeva, and Nikolay F. Stepanov
Citation: The Journal of Chemical Physics 140, 154104 (2014); doi: 10.1063/1.4870936
View online: http://dx.doi.org/10.1063/1.4870936
View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/15?ver=pdfcov
Published by the AIP Publishing
Articles you may be interested in
Two dimensional symmetric correlation functions of the S operator and two dimensional Fourier transforms:
Considering the line coupling for P and R lines of linear molecules
J. Chem. Phys. 140, 104304 (2014); 10.1063/1.4867417
Mixed quantum/classical theory of rotationally and vibrationally inelastic scattering in space-fixed and body-fixed
reference frames
J. Chem. Phys. 139, 174108 (2013); 10.1063/1.4827256
Molecular structure calculations: A unified quantum mechanical description of electrons and nuclei using
explicitly correlated Gaussian functions and the global vector representation
J. Chem. Phys. 137, 024104 (2012); 10.1063/1.4731696
Derivation of the electronic nonadiabatic coupling field in molecular systems: An algebraic-vectorial approach
J. Chem. Phys. 120, 8420 (2004); 10.1063/1.1691394
Non-IPR fullerenes: C 36 and C 72
AIP Conf. Proc. 486, 179 (1999); 10.1063/1.59851
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
THE JOURNAL OF CHEMICAL PHYSICS 140, 154104 (2014)
Determination of the Eckart molecule-fixed frame by use of the apparatus
of quaternion algebra
Sergey V. Krasnoshchekov, Elena V. Isayeva, and Nikolay F. Stepanov
Chemistry Department, Lomonosov Moscow State University, Moscow 119991, Russian Federation
(Received 10 January 2014; accepted 31 March 2014; published online 16 April 2014)
The problem of determining the Eckart molecule-fixed frame for an arbitrary molecule with nuclei
displaced from the equilibrium positions is considered. The solution of the problem is formulated by
minimizing the sum of mass-weighted squared deviations (MWSD) of the nuclei of a displaced configuration from the nuclei of the equilibrium configuration. A mathematical proof of the equivalence
of Eckart conditions and the minimum of MWSD is given. It is shown that the extrema of the sum
of MWSD coincide with eigenvalues of a special 4 × 4 symmetric matrix. Its minimal eigenvalue
corresponds to the desired solution, and the respective eigenvector can be treated as the quaternion
containing the necessary information for rotating the original coordinate system and aligning its axes
with the molecule-fixed coordinate system. A detailed scheme for an efficient numerical implementation of the method is provided, and a numerical example is given. © 2014 AIP Publishing LLC.
[http://dx.doi.org/10.1063/1.4870936]
I. INTRODUCTION
The construction of the vibration-rotation Hamiltonian
operator plays a central role both in molecular spectroscopy
and in a substantial number of other problems of quantum and
physical chemistry. Introduction of an appropriate coordinate
system for a molecule under study is one of the key elements
in the theoretical treatment.1 It is a non-trivial problem to find
the proper placement of coordinate axes for a molecule where
one or more nuclei are shifted from the equilibrium positions
by finite distances. Since in the general case the vibrational
and rotational movements of a molecule as a whole cannot
be fully separated, a customary compromise is to choose a
molecule-fixed coordinate system (MFCS) where vibration
and rotation can be separated to the maximum extent. It is
well known that in the framework of the usual approximation of small displacements of nuclei an efficient separation
of motions can be achieved when the Eckart conditions are
satisfied.2, 3 The Eckart frame chosen this way is the one most
universally used, and its efficient determination is the subject of this article. For large amplitude motion, however, the
Eckart embedding may not be the best option. In this case a
flexible Eckart-Sayvetz embedding can be introduced.4 A recent good review on the body-fixed reference frames can be
found in Ref. 5.
Determination of the orientation of the Eckart frame is
also required in many applied problems: for instance, numerical differentiation of vector and tensor properties (such as
dipole moment and polarizability) in Cartesian coordinates.
As was shown in Ref. 6, sometimes the choice of the coordinate system for the representation of surfaces of the dipole
moment components substantially affects calculated intensities of vibrational transitions. In a simple case of a triatomic molecule, an exact (within the Born-Oppenheimer approximation), body-fixed Hamiltonian for the nuclear motions
0021-9606/2014/140(15)/154104/7/$30.00
can be constructed.7 A variational solution of the vibrationrotation problem of a more complex molecule using the curvilinear coordinate Meyer-Günthard-Pickett8, 9 kinetic energy
operator with the Eckart frame requires evaluation of the rovibrational G matrix for configurations displaced along curvilinear coordinates. For this purpose it is necessary to know the
Eckart frame Cartesian coordinates for displaced molecular
shapes.10 There are many other chemical applications of the
Eckart frame, for example, in the molecular dynamics simulations for large biomolecules.11
A voluminous literature is devoted to methods for determining the Eckart frame in special cases, for example, for
triatomic,12–15 planar,16–18 and highly symmetric molecules.19
Such approaches are convenient in particular situations, but
they do not solve the problem in general. Although the problem of finding the Euler angles to satisfy Eckart conditions
was posed a long time ago and its importance for molecular applications was recognized, the solution is still actively
discussed in the literature.20–33 In this paper we, on the basis
of existing work and new developments, provide all the necessary proofs and details of the most general and convenient
method for finding the Eckart frame. This method is based on
the apparatus of quaternions, which is a particular case of the
geometric algebra.
A convenient and standard method of specifying the orientation of the axes of an Eckart frame is a set of Euler
angles: ζ¯ = (θ, ϕ, χ ), 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π, 0 ≤ χ < 2π .
These angles unambiguously define the form of the orthogonal matrix U = U(θ , ϕ, χ ), also called “the direction
cosine matrix,” for rotation of the coordinate system of equilibrium configuration of an N-atomic molecule.1 After a rotation specified by matrix U, the Cartesian coordinates of
atoms of the distorted molecular configuration r¯ = (x, y, z)
are transformed into a new representation r¯ ∗ = (x ∗ , y ∗ , z∗ )
so that the new coordinate system coincides with the Eckart
140, 154104-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
154104-2
Krasnoshchekov, Isayeva, and Stepanov
J. Chem. Phys. 140, 154104 (2014)
frame:
⎛
xa∗
⎞
⎛
cos θ · cos ϕ · cos χ − sin ϕ · sin χ
⎜ ∗⎟ ⎜
⎜ y ⎟ = ⎜ −cos θ · cos ϕ · sin χ − sin ϕ · cos χ
⎝ a⎠ ⎝
sin θ · cos ϕ
za∗
cos θ · sin ϕ · cos χ + cos ϕ · sin χ
−cos θ · sin ϕ · sin χ + cos ϕ · cos χ
sin θ · sin ϕ
The Euler-angle representation of the rotation matrix is
not the only one possible. The U matrix can alternatively be
represented34 in a form of an exponential series, U = exp (K),
where the argument of the exponent is a 3 × 3 antisymmetric
matrix containing three rotational parameters k¯ = (k1 , k2 , k3 );
specifically
⎞
⎛
0
−k3 k2
⎟
⎜
0
−k1 ⎟
(2)
K=⎜
⎠.
⎝ k3
−k2 k1
0
The Eckart conditions (we use this term as an abbreviation
for “the rotational Eckart conditions”) can be expressed in
the form of three equations (for each axis of the coordinate
system):2–4
N
ma (R¯ a × r¯a∗ ) = 0,
(3)
a=1
where ma are the masses of atoms, and R¯ a , r¯a∗ are the radiusvectors of atoms in the equilibrium configuration and of distorted atoms after rotation to the desired Eckart frame.
Since the elements of the rotational matrix U, the direction cosine matrix, are expressed through Euler angles in the
form of complex trigonometric relationships in Eq. (1), the
Eckart conditions in Eq. (3) gain the form of transcendental
equations, which are not directly solvable in analytic form in
a general case. One possible method of solving this system of
equations is direct minimization of the sum of discrepancies
using one of the existing numerical methods for finding an extremum of a function of several variables. However, in such
a case, ambiguity in the solution can be encountered. Indeed,
it is well known that when the Euler angle θ is equal to zero,
the two remaining angles can vary in a coordinated way, ϕ
= −χ , so that the direction cosine matrix will remain intact.
Besides, numerical minimization of a function requires rather
complicated computer code.
In general, two approaches have been developed for finding the Eckart frame. The first one, which we shall call the
matrix method, was proposed in 1970 by Pickett and Strauss
in Ref. 20 and was further developed in Refs. 21, 22, 24, 25,
27, 28, and 33. It is based on the construction of the following
matrix:
Aij =
N
ma Rai raj ,
(4)
a=1
where Rai and raj are the components of the vectors R¯ a , r¯a . It
can be shown that the Eckart conditions are satisfied if ma-
−sin θ · cos χ
⎞⎛
xa
⎞
⎟⎜ ⎟
⎜ ⎟
sin θ · sin χ ⎟
⎠ ⎝ ya ⎠ .
cos θ
za
(1)
trix A is symmetric.35 If not symmetric, this matrix can be
symmetrized by means of an orthogonal matrix that coincides
with the direction cosine matrix U,
AU = S ≡ S† .
(5)
Since the matrix A in Eq. (4) can become singular or have
its eigenvalues close to each other, unwanted numerical instability and ambiguity can arise. These effects can complicate the method and make it subject to potential errors.
As stated in Refs. 22, 25, and 33, there are eight possible
raw solutions produced by the matrix method, of which four
belong to SO(3)—the group of proper rotations. Dymarsky
and Kudin25 proposed a convenient way to choose the required transformation matrix U that is the closest to I ∈ O(3)
(the identity matrix). This method was successfully employed recently in Ref. 33 for the numerical construction of
the curvilinear internal coordinate Hamiltonian for the NH3
molecule.
The second method of finding the Eckart frame, what we
shall call the method of minimal deviation,36,26, 29 is less actively discussed in the literature of rotation-vibration spectroscopy in comparison to the matrix method. The main
point of this method consists in the minimization of the
sum of mass-weighted squared deviations (MWSD) between
pairs of equivalent atoms for shifted nuclear and equilibrium configurations by means of an appropriate rotation of
the coordinate system through the choice of the orthogonal
matrix U,
min
U ∈SO(3)
N
ma |R¯ a − (U¯ra )|2 .
(6)
a=1
As was shown in Refs. 23, 26, and 37, imposing of
the condition of Eq. (6) ensures satisfying the Eckart
conditions.
Kudin and Dymarsky26 pointed out that the method can
be conveniently implemented by using quaternion algebra, referring to the early (1989) work of Kearsley36 and to the software code of Heisterberg, which was available on-line but not
described in a publication. The quaternion implementation of
the minimal deviation method was referred to as a universal
one due to the absence of special cases.26 Our analysis shows
that Heisterberg’s computer code differs from the Kearsley’s
algorithm in spite of similar features.
Generally speaking, the concept of the root-mean-square
deviation (RMSD) between molecular configurations of similar molecules has a much wider field of application than
theoretical studies of vibration-rotation spectra. This concept is used when superimposing structures for biochemistry,
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
154104-3
Krasnoshchekov, Isayeva, and Stepanov
J. Chem. Phys. 140, 154104 (2014)
crystallography, and macromolecular dynamics studies.38–40
The minimization of the RMSD can be performed using
quaternion-based algorithms. Several approaches have been
proposed and discussed, including the method of Lagrange
multipliers or finding roots of the characteristic quartic
polynomials.29, 36, 39–46 Good reviews of rotational superposition and of using quaternions in molecular modeling can be
found in Refs. 47 and 48, respectively. Our analysis shows
that none of the previously proposed methods is the optimum
one. In addition, we could not find any publication where the
quaternion-based approach for finding Eckart frame is systematically applied in the context of theoretical studies of
vibration-rotation spectra.
Kearsley’s approach36 is probably one of the best existing
methods for minimizing RMSD. It is based on quaternion algebra and is targeted at the solution of a crystallographic problem, where it is unnecessary to take into account the atomic
masses. Besides, there are a few misprints in Ref. 36. For
determining the Eckart frame, the theoretical background of
Kearsley’s approach must be clarified, and its optimum algorithmic implementation must be developed.
There is a general need for a universal method for finding
the MFCS that would effectively work with molecules of
arbitrary size and symmetry properties. Such a method must
be convenient for software implementation, be compact and
free of numerical instability, ambiguity and dependence on
special cases, and be able to provide a solution with full
numerical accuracy. In the present work, we consider the
necessary theoretical foundations of the problem of finding
the Eckart frame using the method of minimal deviation
and give proofs and intermediate formula derivations. The
algorithmic scheme presented is targeted for the most general
case, is specified in detail, and is illustrated by a numerical
example.
⎛
x∗
⎞
⎛
q02 + q12 − q22 − q32
⎜ ∗⎟
⎜
⎝ y ⎠ = U¯r ≡ ⎝ 2(q1 q2 − q0 q3 )
z∗
2(q1 q3 + q0 q2 )
a=1
ma |R¯ a − r¯a∗ |2
N
i
[0, r¯ ∗ ] = Q−1 [0, r¯ ]Q
¯ q¯ + 2q0 (¯r × q)
¯ − q¯ × r¯ × q¯ ,
= 0, q02 r¯ + (¯r · q)
(8)
where the vectors are treated as quaternions with a zero scalar
component. Using this definition, it is easy to show that the
rotation of vector r¯ in the matrix form is given by the components of the quaternion in the following way:
⎞⎛ ⎞
x
⎟⎜ ⎟
2(q2 q3 + q0 q1 ) ⎠ ⎝ y ⎠ .
z
q02 − q12 − q22 + q32
2(q1 q3 − q0 q2 )
(9)
where the summation over the index i is conducted along three
coordinate axes. In this expression only the last term in square
brackets depends on the matrix of the coordinate transformation U. Let us write down this term in explicit form (accurate to a constant multiplier) taking into account the equality
r¯ ∗ = U r¯ :
N
i
∗
ma (Rai rai
)=
a=1
N
a=1
=
a=1
=
¯
¯ 0 , q]
¯ = [p0 q0 − p¯ · q,
¯ p0 q¯ + q0 p¯ + p¯ × q].
P Q = [p0 , p][q
(7)
The squared norm of the quaternion is equal to |Q|2 = (q02
¯ while the norm of the product of quaternions is equal
+ q¯ · q),
to the product of their norms. The inverse quaternion is de2
¯
.
fined as follows: Q−1 = [q0 , −q]/|Q|
¯ |Q|2 = 1,
The quaternion with unit norm, Q = [q0 , q],
can be used for rotation of the vector r¯ = (x, y, z) into a new
position r¯ ∗ = (x ∗ , y ∗ , z∗ ) expressed in the original axes system in the following manner:
2(q2 q3 − q0 q1 )
With the designation r¯ ∗ = U r¯ the sum of mass-weighted
squared deviations in Eq. (6) can be transformed using the
equality |¯r ∗ |2 = |¯r |2 ,
N
A quaternion is a generalization of an imaginary number
for the four-dimensional case. A quaternion can be considered as a row-vector with four components, the first of which
is a scalar value, while the other three are components of
¯ A
a three-dimensional vector: Q = (q0 , q1 , q2 , q3 ) = [q0 , q].
quaternion can also be considered as a rotation of the threedimensional space about the axis, specified by the vector
part of the quaternion, by the angle, specified by its scalar
part. The essential properties of quaternions, necessary for the
problem under consideration, are given below.
The product of two quaternions is a quaternion, and it is
given by the following expression:
q02 − q12 + q22 − q32
B. The equivalence of rotational Eckart conditions
and the minimum of the MWSD
ma |R¯ a − (U¯ra )|2 =
A. Essential properties of quaternions
2(q1 q2 + q0 q3 )
The determinant of the matrix of proper rotation U is positive,
thereby avoiding inversion of the coordinate system.
N
II. THEORETICAL FOUNDATIONS
N
ma
Rai
i
ma R¯ a† U r¯a .
Uij raj
j
(11)
a=1
∗
ma [|Rai |2 + |rai |2 − 2Rai rai
],
a=1
(10)
The extremum of this expression with respect to variations
of three variables ζ¯ (see Eqs. (1) and (2)) defining the form
of a matrix U(ζ¯ ) is specified by equating its derivatives to
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
154104-4
Krasnoshchekov, Isayeva, and Stepanov
J. Chem. Phys. 140, 154104 (2014)
zero,
where indices i, j, and an additional index k are obtained from
the cyclic permutation of Cartesian axes x, y, z. Three conditions from Eq. (17) together with the definition r¯ ∗ = U r¯ can
finally be written in the following vector form:
N
∂ ma R¯ a† U(ζ¯ )¯ra
∂ ζ¯ a=1
=
N
∂
ma R¯ a† U ζ¯ r¯a = 0. (12)
(U(ζ¯ )) r¯a ≡
∂ ζ¯
N
ma R¯ a†
a=1
a=1
This expression is further transformed into the following form
after inserting the identity matrix U† U = I and using the abbreviation Uζ¯ U† = X,
N
ma R¯ a† Uζ¯ U† U¯ra
a=1
=
N
ma R¯ a† (Uζ¯ U† )¯ra∗ =
a=1
N
ma R¯ a† X¯ra∗ = 0. (13)
∂
(UU† ) = Uζ¯ U† + U(U† )ζ¯ = X + X† = 0.
∂ ζ¯
(14)
Hence, the matrix X is skew-symmetric, X = −X† , and xij
= −xji .
Now, once this important property of the matrix X is
established, we prove that the matrix U that minimizes the
functional in Eq. (10) coincides with the direction cosine
matrix of the MFCS so that the Eckart conditions are also
satisfied.
Since the result of summation does not depend on the
designation of indices of summation, the following equality
will be true:
∗
ma Rai Xij raj
=
N 3 3
a=1 i=1 j =1
∗
ma Raj Xj i rai
. (15)
a=1 i=1 j =1
Both sides of this equality are equal to the left side of the conditions in Eq. (13) with independent derivatives with respect
to variables ζ¯ . Therefore, putting both versions together and
also taking into account the skew-symmetric property of the
matrix X, one obtains the following condition:
N 3 3
N 3 3
ma R¯ a × (U¯ra ) = 0,
(18)
a=1
C. Equivalence of the condition of zeroing the
derivative of the residual MWSD functional and the
eigenvalue problem formulation
As was shown in Kearsley’s paper,36 the search for the
proper rotation matrix ensuring the minimum of the sum of
squared deviations between corresponding atoms can be reduced to the eigenvalue problem using the quaternion algebra. In this method, quaternions are used for representing the
vector positions of atoms and for specifying the orthogonal
matrices for coordinate system rotations.
Let us represent the radius-vectors of nuclei in the equilibrium and shifted configurations as purely vector quaternions, i.e., quaternions with a zero scalar component: [0, R¯ a ]
and [0, r¯a ], respectively. Then, let Q = [q1 , q2 , q3 , q4 ] be the
quaternion, corresponding to the required rotation of the coordinate system. Since the equilibrium and shifted configurations of a molecule cannot be aligned in a general case, after the rotation of the coordinate system, each pair of equivalent atoms will correspond to the “residual” purely vector
quaternion,
[0, δ¯a ] = [0, R¯ a ] − Q−1 [0, r¯a ]Q.
∗
∗
ma (Rai raj
− Raj rai
)Xij = 0.
∗
∗
ma (Rai raj
− Raj rai
)=
N
a=1
ma (R¯ a × r¯a∗ )k = 0,
(19)
In accordance with the principle of minimum total MWSD,
the quaternion Q is defined as the rotation that minimizes the
mass-weighted sum of squared “residual” vectors,
N
a=1
ma |[0, δ¯a ]|2 =
N
ma |δ¯a |2 = min .
(20)
a=1
(16)
Since Xij = −Xji , Xii = 0, and the elements of matrix X are
independent for j ≥ i + 1, this condition can be written in the
following form using the vector product notation:
a=1
N
which can be immediately recognized as the rotational Eckart
conditions.
We have just proved that the matrix U minimizing the
sum in Eq. (10) also ensures fulfillment of the Eckart conditions. In Subsection II C we will show that the problem of
minimization of the sum in Eq. (10) can be formulated in
closed form in terms of the eigenvalues and eigenvectors of
a special 4 × 4 matrix.
ε(Q) =
a=1 i=1 j =1
N
a=1
∗
∗
ma (Rai Xij raj
+ Raj Xj i rai
)
a=1 i=1 j =1
=
ma (R¯ a × r¯a∗ ) =
a=1
For any orthogonal matrix U, (U† U = UU† = I) the following
equalities are obeyed:
N 3 3
N
(17)
The norm of the quaternion Q must be equal to unity because
of orthogonality of the transformation that rotates the coordinate system. Since the norm of a product of quaternions is
equal to the product of their norms, the condition of the minimum in Eq. (20) will not be violated if elementary quaternions
are multiplied from the left side by Q:
Q[0, δ¯a ] = Q[0, R¯ a ] − [0, r¯a ]Q.
(21)
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
154104-5
Krasnoshchekov, Isayeva, and Stepanov
J. Chem. Phys. 140, 154104 (2014)
In this case the sum of mass-weighted squared residual deviations ε(Q) will take the form,
ε(Q) =
N
ma |Q[0, δ¯a ]|2
a=1
=
N
It can be shown in a straightforward manner that the functional ε(Q) can be represented by the following quadratic
form:
⎛ ⎞† ⎛
⎞⎛ ⎞
q1
q1
C11 C12 C13 C14
⎜
⎜q ⎟ ⎜C
⎟
⎟
⎜ 2 ⎟ ⎜ 21 C22 C23 C24 ⎟ ⎜ q2 ⎟
ε(Q) = ⎜ ⎟ ⎜
⎟ ⎜ ⎟ , (23)
⎝ q3 ⎠ ⎝ C31 C32 C33 C34 ⎠ ⎝ q3 ⎠
ma |[−q¯ · (R¯ a − r¯a ), q1 (R¯ a − r¯a )
q4
a=1
+q¯ × (R¯ a + r¯a )]|2 .
C11 =
C13 =
C22 =
C24 =
N
a=1
N
a=1
N
a=1
N
(22)
2
2
2
,
ma x−a
+ y−a
+ z−a
ma (x−a z+a − x+a z−a ),
2
2
2
ma x−a
+ y+a
+ z+a
,
ma (x−a z−a − x+a z+a ),
C34 =
C42
C43
C44
q4
where the matrix elements of the matrix C are given by the
following formulas:
C12 =
C14 =
C23 =
C33 =
a=1
N
C41
N
a=1
N
a=1
N
a=1
N
ma (y+a z−a − y−a z+a ),
ma (x+a y−a − x−a y+a ),
ma (x−a y−a − x+a y+a ),
(24)
2
2
2
,
ma x+a
+ y−a
+ z+a
a=1
ma (y−a z−a − y+a z+a ),
C44 =
a=1
2
2
2
ma x+a
+ y+a
+ z−a
.
a=1
The following notation is introduced here: ξ +a = a
+ ξ a , ξ −a = a − ξ a , = X, Y, Z, where ξ = x, y, z. Here capital Greek letters signify the atomic coordinates of the original
equilibrium system and lowercase letters signify the displaced
atomic coordinates in the new center-of-mass system.
If an eigenvector V of the matrix C is chosen as a Qvector in Eq. (23), the value of this quadratic form is reduced
to the corresponding eigenvalue λ multiplied by the squared
norm of the eigenvector, i.e.,
ε(V ) = V † CV = V + λV = λ|V |2 .
N
(25)
Since eigenvectors are usually normalized to unity, the value
of the functional ε(V ) is merely the corresponding eigenvalue. The eigenvalues of the matrix C are all non-negative,
which can be easily proven. It is evident that the minimum
value of the sum of the mass-weighted squared residual deviations is equal to the minimum eigenvalue of the matrix C.
The corresponding eigenvector can be treated as the quaternion, specifying the orthogonal transformation to the Eckart
frame.
The solution of the problem of minimum mass-weighted
deviation between two molecular configurations reduces to
the solution of a numerically stable and trivial eigenvalue
problem. The eigenvalue problem for the 4 × 4 symmetric
matrix can be solved through finding the roots of the characteristic quartic polynomial.39 These roots can be found analytically, so that the whole problem of finding the Eckart frame
can, in principle, be solved in a purely analytic closed form.
III. NUMERICAL EXAMPLE
As an illustration of the feasibility of the proposed technique and in the interest of simplifying debugging of its future software implementations, we will present the numerical
example for the water molecule, H2 O. All calculations were
performed using Fortran code with double precision for floating point numbers (64 bits). The geometric parameters of the
molecule were calculated quantum-mechanically49 using the
MP2 method with the aug-cc-pVTZ basis set: the bond length
re (OH) = 0.96137027 Å, and the valence angle (HOH)
= 127.94547◦ . Atomic coordinates in the principal axes system (the Z-axis is aligned with the C2 symmetry axis) are presented in Table I. In the search for the Eckart frame, the atom
TABLE I. Cartesian coordinates (x, y, z, in Å units) of water molecule in the
equilibrium, the shifted, the center-of-mass, and the molecule-fixed (Eckart
frame) coordinate systems. (The mass of oxygen atom is 15.99491502 amu,
and the mass of hydrogen atoms is 1.00782522 amu.)
Atom
Equilibrium
O (1)
0.0
0.0
0.06615931
H (2)
0.0
0.75813308
− 0.52499806
H (3)
0.0
− 0.75813308
− 0.52499806
Shift
0.05000000
0.05000000
0.0
Center of mass
Eckart frame
0.00559574
0.00559574
0.06615931
0.0
0.00362754
0.06653210
0.0
0.0
0.0
− 0.04440426
0.71372882
− 0.52499806
0.0
0.72901492
− 0.50551041
0.0
0.0
0.0
− 0.04440426
− 0.80253733
− 0.52499806
0.0
− 0.78658654
− 0.55040211
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
154104-6
Krasnoshchekov, Isayeva, and Stepanov
J. Chem. Phys. 140, 154104 (2014)
of oxygen was arbitrarily displaced by the vector (0.05, 0.05,
0.0). Eight significant figures are used for comparing results
of different implementations of the method.
The optimization of the target function as the sum of
the discrepancies in three Eckart conditions by the BroydenFletcher-Goldfarb-Shanno numerical method50 gave the following values of Euler angles (in degrees): θ = 5.1229451, ϕ
= 19.363965, χ = −19.292339. These values correspond to
the residual root mean square discrepancy in Eq. (3) of the order of ≈10−13 . The direction cosine matrix, obtained from the
⎛
optimized Euler angles using the definition of Eq. (1), takes
the following form:
⎞
⎛
0.99644220
0.0
−0.08427893
⎟
⎜
U = ⎝ −0.00249522 0.99956162 −0.02950141 ⎠ .
0.08424198
0.02960675
(26)
The direct method of calculation of the matrix of direction
cosines begins from the calculation of the symmetric matrix
of Eq. (24) and leads to the following numerical result:
⎞
0.00895035
⎜
⎜ 0.10582126
⎜
C=⎜
⎜ −0.10582126
⎝
0.0
0.99600539
⎟
⎟
⎟
⎟.
⎟
⎠
7.14533839
0.0
2.51123114
−0.10582126
−0.10582126
The diagonalization of this matrix yields the following
eigenvalues in descending order: λ1 = 7.15137244, λ2
= 4.64385230, λ3 = 2.51043643, and λ4 = 0.00291629.
The minimum eigenvalue (λ4 ) corresponds to the eigenvector
V4 = (0.99900065, −0.01479182, 0.04217237, 0.00062443).
The application of Eq. (9) for construction of the direction
cosine matrix using the corresponding quaternion produces
the result, which is numerically identical to the one obtained
above from the optimized Euler angles. The residual discrepancy in the Eckart conditions with these values is a quantity
of the order of ≈10−17 .
(27)
4.64305759
plies to larger molecules, are seemingly underappreciated in
vibration-rotation spectroscopy studies. This method is not
only a universal one, but it is also fast and well suited for numerical computations. The scope of its applications is wider
than molecular spectroscopy; it can be equally used in biochemistry and crystallography.
In the current paper, all necessary theoretical considerations for the implementation of the quaternion method are
assembled, a number of theoretical relationships are proven
and all steps necessary for practical solution of the problem
are specified. The method has been coded in Fortran and illustrated by a numerical example.
IV. CONCLUSIONS
The choice of the Eckart frame attached to a molecule
with atoms shifted from the equilibrium positions is one of the
important problems of the quantum mechanics of molecules
and molecular spectroscopy. The proper choice of the Eckart
frame ensures the optimum separation of vibration and rotation and can be fulfilled in accordance with the Eckart conditions. The orthogonal transformation, responsible for the rotation of the equilibrium coordinate system and aligning its
axes to the ones of the Eckart frame, can be described by a set
of three Euler angles or simply by the elements of the 3 × 3
orthogonal direction cosine matrix. In a general case, specifying the Eckart frame is a non-trivial problem.
The Eckart frame can be found analytically in a number of simple cases (such as triatomic, planar, and selected
highly-symmetric molecules). In many studies of vibrational
and vibration-rotation spectra of small molecules there is no
need for developing the most efficient and general computational method for finding the direction cosine matrix specifying the Eckart frame. In this situation a matrix method is
commonly used that is based on symmetrization of a special 3 × 3 matrix.10, 33 Certain advantages of the alternative minimal-deviation, quaternion-based method, which ap-
ACKNOWLEDGMENTS
The authors are deeply indebted to Professor Norman C.
Craig for his general invaluable help in reading the manuscript
and in making useful suggestions.
1 P.
R. Bunker and P. Jensen, Molecular Symmetry and Spectroscopy (NRC
Research Press, Ottawa, 1998).
2 C. Eckart, Phys. Rev. 47, 552 (1935).
3 J. D. Louck and H. W. Galbraith, Rev. Mod. Phys. 48, 69 (1976).
4 A. Sayvetz, J. Chem. Phys. 7, 383 (1939).
5 A. V. Meremianin, J. Chem. Phys. 120, 7861 (2004).
6 C. R. Le Sueur, S. Miller, J. Tennyson, and B. T. Sutcliffe, Mol. Phys. 76,
1147 (1992).
7 B. T. Sutcliffe and J. Tennyson, Int. J. Quantum Chem. 39, 183 (1991).
8 R. Meyer and H. H. Günthard, J. Chem. Phys. 49, 1510 (1968).
9 H. M. Pickett, J. Chem. Phys. 56, 1715 (1972).
10 X.-G. Wang and T. Carrington, Jr., J. Chem. Phys. 138, 104106
(2013).
11 G. Chevrot, P. Calligari, K. Hinsen, and G. R. Kneller, J. Chem. Phys. 135,
084110 (2011).
12 G. A. Natanson, Mol. Phys. 66, 129 (1989).
13 H. Wei and T. Carrington, J. Chem. Phys. 107, 2813 (1997).
14 H. Wei and T. Carrington, J. Chem. Phys. 107, 9493 (1997).
15 H. Wei and T. Carrington, Chem. Phys. Lett. 287, 289 (1998).
16 G. A. Natanson, Chem. Phys. Lett. 121, 343 (1985).
17 H. Wei, J. Chem. Phys. 118, 7202 (2003).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46
154104-7
Krasnoshchekov, Isayeva, and Stepanov
J. Chem. Phys. 140, 154104 (2014)
18 H.
36 S.
19 F.
37 N.
Wei, J. Chem. Phys. 118, 7208 (2003).
O. Meyer III and R. W. Redding, J. Mol. Spectrosc. 70, 410 (1978).
20 H. M. Pickett and H. L. Strauss, J. Am. Chem. Soc. 92, 7281 (1970).
21 R. W. Redding and J. T. Hougen, J. Mol. Spectrosc. 37, 366 (1971).
22 N. J. D. Lucas, J. Phys. B: At. Mol. Phys. 6, 155 (1973).
23 F. Jorgensen, Int. J. Quantum Chem. 14, 55 (1978).
24 R. W. Redding and F. O. Meyer III, J. Mol. Spectrosc. 74, 486 (1979).
25 A. Y. Dymarsky and K. N. Kudin, J. Chem. Phys. 122, 124103 (2005).
26 K. N. Kudin and A. Y. Dymarsky, J. Chem. Phys. 122, 224105 (2005).
27 M. Dierksen, J. Chem. Phys. 122, 227101 (2005).
28 A. Y. Dymarsky and K. N. Kudin, J. Chem. Phys. 122, 227102 (2005).
29 G. R. Kneller, J. Chem. Phys. 128, 194101 (2008).
30 M. E. Castro, A. Niño, and C. Muñoz-Caro, Comput. Phys. Commun. 181,
967 (2010).
31 T. Szidarovszky, C. Fábri, and A. G. Császár, J. Chem. Phys. 136, 174112
(2012).
32 A. V. Meremianin, J. Math. Chem. 51, 1376 (2013).
33 C. Fábri, E. Mátyus, and A. G. Császár, Spectrochim. Acta, Part A 119, 84
(2014).
34 W. D. Allen and A. G. Császár, J. Chem. Phys. 98, 2983 (1993).
35 See supplementary material at http://dx.doi.org/10.1063/1.4870936 for the
proof.
K. Kearsley, Acta Crystallogr., Sect. A 45, 208 (1989).
F. Stepanov and V. I. Pupyshev, Quantum Mechanics of Molecules
and Quantum Chemistry (Moscow State University Publishing House,
Moscow, 1991).
38 M. Praprotnik and D. Janežiˇ
c, J. Chem. Inf. Model. 45, 1571 (2005).
39 D. L. Theobald, Acta Crystallogr., Sect. A 61, 478 (2005).
40 P. Liu, D. K. Agrafiotis, and D. L. Theobald, J. Comput. Chem. 31, 1561
(2010).
41 T. F. Havel and I. Najfeld, J. Mol. Struct. THEOCHEM 308, 241 (1994).
42 E. A. Coutsias, C. Seok, and K. A. Dill, J. Comput. Chem. 25, 1849 (2004).
43 G. R. Kneller, J. Comput. Chem. 26, 1660 (2005).
44 E. A. Coutsias, C. Seok, and K. A. Dill, J. Comput. Chem. 26, 1663 (2005).
45 G. R. Kneller, J. Comput. Chem. 32, 183 (2011).
46 P. Liu, D. K. Agrafiotis, and D. L. Theobald, J. Comput. Chem. 32, 185
(2011).
47 D. R. Flower, J. Mol. Graphics Modell. 17, 238 (1999).
48 C. F. F. Karney, J. Mol. Graphics Modell. 25, 595 (2007).
49 M. J. Frisch, G. W. Trucks, H. B. Schlegel et al., Gaussian 03, Revision
C.02, Gaussian, Inc., Wallingford, CT, 2004.
50 W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes. The Art of Scientific Computing (Cambridge University Press,
Cambridge, 1986).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP:
85.140.189.227 On: Thu, 17 Apr 2014 05:10:46

Download Report