Diffraction Measurements and Equilibrium Parameters

Abstract

Structural studies are largely performed without taking into account vibrational effects or with incorrectly taking them into account. The paper presents a first-order perturbation theory analysis of the problem. It is shown that vibrational effects introduce errors on the order of 0.02 Å or larger (sometimes, up to 0.1-0.2 Å) into the results of diffraction measurements. Methods for calculating the mean rotational constants, mean-square vibrational amplitudes, vibrational corrections to internuclear distances, and asymmetry parameters are described. Problems related to low-frequency motions, including torsional motions that transform into free rotation at low excitation levels, are discussed. The algorithms described are implemented in the program available from the author (free).

1. Introduction

Diffraction measurements yield some time- and ensemble-averaged parameter values which do not necessarily coincide (generally, almost never coincide) with equilibrium geometric characteristics. This is a consequence of vibrational motions, which remain unfrozen even at absolute zero. For instance, the CO₂ molecule is bent on average because of bending vibrations, and the mean O⋯O distance in this linear molecule is therefore smaller than the sum of C–O bond lengths (“shrinkage” effect [1, 2]). The properties of the CO₂ molecule (first of all, optical, for instance, laser properties) are, however, determined by its equilibrium linear rather than average bent configuration. For this reason, the problem of the reconstruction of equilibrium geometry from diffraction measurement data is very important.

The instantaneous configuration of a vibrational system 𝑅(𝑡) can in the first approximation be represented as𝑅(𝑡)=𝑅_𝑒+^𝑛_𝑖=1𝑋_𝑖cos𝜔_𝑖𝑡+𝜃_𝑖,(1) where 𝑅_𝑒 is the vector of equilibrium geometric parameters, 𝜔_𝑖 and 𝜃_𝑖 are the frequencies and phases of vibrational motions, 𝑋_𝑖 is the vector of geometric parameter changes at 𝜔_𝑖𝑡+𝜃_𝑖=2𝑘𝜋 and 𝜔_𝑗≠𝑖𝑡+𝜃_𝑗=(2𝑘+1)𝜋/2(𝑘=0,…), and 𝑛 is the number of vibrational degrees of freedom. Since the ratios between the frequencies and phases of various vibrations must not necessarily be rational numbers, it is clear that a vibrational system has a chance to assume its equilibrium configuration only once during the whole time of its existence. This requires that all the 𝜔_𝑖𝑡+𝜃_𝑖 values simultaneously satisfy the condition 𝜔_𝑖𝑡+𝜃_𝑖=(2𝑘+1)𝜋/2(𝑘=0,…). For this reason, “measurements’’ of equilibrium parameters are out of the question. It seems to follow from (1) that the mean (measured) geometric parameters of a vibrational system coincide with equilibrium. However, this is not the case even in the first approximation. Indeed, there is no one-to-one correspondence between the configuration determined by (1) and parameters measured experimentally, and the O⋯O distance in CO₂ decreases irrespective of the sign of changes in the ∠OCO angle. The set of mean (measured by diffraction methods) internuclear distances can be inconsistent with any vibrational system configuration. For instance, the set of mean internuclear distances for a tetrahedron corresponds to angles smaller than tetrahedral, which is, of course, geometrically impossible. This, naturally, increases the 𝑅 factor, because refinements are performed for certain structural models.

The shrinkage effect for the carbon dioxide molecule is very small (~0.005 Å). However, this is not the case with more complex systems. For instance, in the carbon suboxide molecule (C₃O₂, O=C=C=C=O), the observed distance between the terminal oxygen atoms is shortened with respect to its equilibrium value by 0.198 Å at 508 K [3], 0.150 Å at 293 K [4], 0.140 Å at 237 K [3] (all these values are experimental), and 0.0325 Å at 0 K (calculated value) [5]. (The experimental radial distribution curve for C₃O₂ from [4] is reproduced in Figure 1 to show that measurement errors are minimum for this molecule.) Naturally, such corrections to diffractionally measured parameters cannot be ignored.

Figure 1 

Experimental radial distribution curve for C3O2 from [4].

This makes it necessary to analyze the influence of vibrational motions on the ensemble-averaged geometric parameters of vibrational systems, primarily internuclear distances, determined by experimental diffraction methods. The spectra of condensed phases are exceedingly complex, and we shall not touch them but shall concentrate on molecular vibrational systems. Nevertheless, the algorithms described below are likely necessary to use in structural analyses of molecular crystals. We shall proceed using the classical formalism and the harmonic approximation as a zero-order step.

2. Equation of Motion

In the adiabatic approximation, the nuclear subsystem of a vibrational system is considered separately. After minor simplifications, the Lagrange equation of motion𝑑𝑑𝑡𝜕𝜕ℒ𝑞̇−𝜕𝜕ℒ𝑞=0(2) (ℒ is the Lagrange function) written for the nuclear subsystem takes the form𝑑𝑑𝑡𝜕𝑇𝜕(𝑞𝑞,̇𝑞̇)−𝜕𝜕𝑞(𝑇(𝑞,𝑞̇)−𝑈(𝑞))=0(3) (𝑈 is assumed to be independent of the velocities of the particles constituting the system). Here, 𝑇(𝑞,𝑞̇) is the additive part of the Lagrange function (the kinetic energy), 𝑈(𝑞) is the potential energy (the nonadditive part), and 𝑞 stands for some “generalized’’ coordinates describing the configuration of the system.

The kinetic and potential energies are written as𝑇=12^𝑁_𝑖=1𝑈=1𝑚_𝑖𝑥̇^2𝑖+𝑦̇^2𝑖+𝑧̇^2𝑖,2^3𝑁−6(5)_𝑖=1𝜕²𝑈𝜕𝑞^𝑖𝜕𝑞^𝑗𝑞^𝑖𝑞^𝑗+⋯,(4) where 𝑁 is the number of vibrational system particles and 3𝑁−6(5) is the number of vibrational degrees of freedom (all the 𝑞 and 𝑥 coordinates are taken to be zero in the equilibrium configuration). In the matrix form, 𝑇=(1/2)𝑀𝑥̇^𝑖𝑥̇^𝑗 and 𝑈=(1/2)𝐹𝑞^𝑖𝑞^𝑗. Here, 𝑀 is the diagonal matrix of atomic masses containing three masses of each atom corresponding to motions along the 𝑥, 𝑦, and 𝑧 axes and 𝐹 is the potential energy operator written in generalized coordinates 𝑞. The 𝐹 operator is symmetrical by construction. In the 𝑞=𝐵𝑥 generalized coordinates (𝐵 is the transformation matrix between the generalized and Cartesian coordinates), the kinetic energy takes the form𝑇=12𝐵⁻¹𝑞̇^𝑗^†𝑀𝐵⁻¹𝑞̇^𝑖=𝐺⁻¹𝑞̇^𝑖𝑞̇^𝑗(5) († is the symbol of transposition). The 𝐺⁻¹ matrix is symmetrical by construction. This matrix is positive definite, because any motion with a nonzero velocity increases the kinetic energy.

The trivial equalities𝑞_𝑖=^3𝑁_𝑗𝜕𝜕𝑥𝑞^𝑖_𝑗𝑥^𝑗+12^3𝑁_𝑗,𝑘𝜕²𝑞_𝑖𝜕𝑥_𝑗𝜕𝑥_𝑘𝑥^𝑗𝑥_𝑘+⋯,(6) that is, 𝐵=𝐵(0)+(1/2)(𝜕𝐵/𝜕𝑥), and𝑈=12^3𝑁−6(5)_𝑖,𝑗𝜕²𝑈𝜕𝑞_𝑖𝜕𝑞_𝑗𝑞^𝑖𝑞_𝑗+16^3𝑁−6(5)_{𝑖,𝑗,𝑘}𝜕³𝑈𝜕𝑞_𝑖𝜕𝑞_𝑗𝜕𝑞_𝑘𝑞^𝑖𝑞^𝑗𝑞^𝑘+⋯,(7) that is, 𝐹=𝐹(0)+(1/3)(𝜕𝐹/𝜕𝑞), allow (3) to be rewritten as𝐺⁻¹(0)𝑞̈^𝑗+12𝜕𝐺𝜕𝑞⁻¹𝑞^𝑖𝑞̈^𝑗+𝑞̇^𝑖𝑞̇^𝑗−12𝜕𝐺𝜕𝑞⁻¹𝑞̇^𝑖𝑞̇^𝑗+𝐹(0)𝑞^𝑖+12𝜕𝜕𝐹𝑞𝑞^𝑖𝑞^𝑗=0,(8) where 𝐵(0), 𝐺⁻¹(0), and 𝐹(0) are the matrices corresponding to the approximation of infinitesimal vibrational amplitudes. The derivative of an 𝑚×𝑛 matrix 𝐴 with respect to an 𝑙-dimensional vector 𝑏 is understood to be a set of 𝑚×𝑛×𝑙(𝜕𝑎_𝑖𝑗/𝜕𝑏_𝑘) values. If the Cartesian system of coordinates is used, 𝐺⁻¹(𝑞)=𝐺⁻¹(0)=𝑀.

In the approximation of infinitesimal amplitudes ((𝜕𝐺⁻¹/𝜕𝑞)𝑞=0 and (𝜕𝐹/𝜕𝑞)𝑞=0), (8) is rewritten as𝐺⁻¹𝑞̈^𝑗+𝐹𝑞^𝑗=0.(9) Solutions to this system of equations are sought in the form 𝑞=𝑙cos(𝜔𝑡+𝜃) (𝑞̈=−𝑙𝜔²cos(𝜔𝑡+𝜃)), that is,𝐺⁻¹𝑙_𝛼𝜔^2𝛼−𝐹𝑙_𝛼=0,𝐺𝐹𝑙_𝛼=𝑙_𝛼𝜔^2𝛼,(10) where 𝑙_𝛼 is the vector determining the change in the configuration of the vibrational system vibrating at the 𝜔_𝛼 frequency (eigenvectors and frequencies are indexed by Greek letters). Since the 𝐺⁻¹ and 𝐹 matrices are symmetrical and the 𝐺⁻¹ matrix is positive definite, they can simultaneously be reduced to the diagonal form by one similarity transformation using the procedure that transforms the positive definite into identity matrix. This is the simplest matrix analysis theorem.

Indeed, it is known that Hermitian matrices can be diagonalized by a similarity transformation with unitary matrices. Let us apply such a transformation to the 𝐺⁻¹ matrix (and, simultaneously, to the 𝐹 matrix: 𝑉^†1𝐺⁻¹𝑉₁=𝐷₁ and 𝑉^†1𝐹𝑉₁=𝐹₁, where 𝐷₁ is some diagonal matrix and 𝐹₁ remains symmetrical (Hermitian) by construction. Since 𝐺⁻¹ is positive definite, all the elements of the 𝐷₁ matrix are positive numbers. We can therefore construct the 𝐷^1−1/2 diagonal matrix such that 𝐷^1−1/2𝐷₁𝐷^1−1/2=𝐼 (𝐼 is the identity matrix). As a result, the kinetic energy matrix transforms into the identity matrix whereas the potential energy matrix 𝐷^1−1/2𝐹₁𝐷^1−1/2=𝐹₂ remains symmetrical. It can therefore be diagonalized with the use of a unitary (orthogonal) transformation, 𝑉^†2𝐷^1−1/2𝑉^†1𝐹𝑉₁𝐷^1−1/2𝑉₂=𝐷₂. As to the identity matrix, its unitary transformation leaves it unchanged, 𝑉^†2𝐼𝑉₂=𝐼.

Let the 𝐿 matrix be constructed as described above (𝐿=𝑉₁𝐷^1−1/2𝑉₂). We have𝐿⁻¹𝐺𝐿⁻¹^†𝐺𝐿𝐹^†𝐿𝐹𝐿==𝐿Ω𝐿,⁻¹𝐺𝐹𝐿=Ω,(11) where Ω is the diagonal matrix of the squares of vibrational frequencies and 𝐿 is the matrix whose columns are the eigenvectors of the vibrational problem, 𝑙_𝛼. We obtained a not very convenient (non-Hermitian) Hamilton function, which is more simple to analyze by classical methods, bearing in mind that classical and quantum results coincide at the level of first-order perturbation theory.

As distinct from all the matrices introduced above, the 𝐿 matrix is independent of the instantaneous configuration of the system. The columns of this matrix are mere linear combinations of generalized coordinates 𝑞 selected to describe the geometry of the molecule. In what follows, it will be more convenient to use the matrices 𝔊𝔅==𝔅𝑀𝐿⁻¹₋₁𝐵,𝔉=𝐿_†𝔅^†,𝐹𝐿.(12) In the corresponding system of coordinates, 𝐺(0) transforms into the identity matrix (𝐼); 𝐹(0), into the diagonal matrix of eigenvalues Ω; eigenvectors, into unit vectors 𝑙_𝛼→𝑒_𝛼 (𝑒_𝛼 is the vector whose all components except the 𝛼th component are zero, and the 𝛼th component is one; for the transition to the quantum normalization in the transformations described below, it is more convenient to set the 𝛼th element equal to [(ℎ/4𝜋²𝑐𝜈_𝛼)coth(ℎ𝑐𝜈_𝛼/2𝑘𝑇)]^1/2 rather than one [6]; this is implied in what follows. Here, ℎ is the Planck constant, 𝑐 is the velocity of light, 𝑘 is the Boltzmann constant, and 𝜈_𝛼 is the 𝛼th frequency).

Solutions to (8) will be sought in the form 𝑠=𝑠₀+𝑠₁ (𝑠_0,𝛼=𝑒_𝛼cos(𝜔_𝛼𝑡+𝜃_𝛼)) and 𝜔=𝜔₀+𝜔₁. Let us rewrite (9) and (8) using the notation introduced above and subtract the former from the latter, 𝐼𝑠^𝑗1̈+Ω𝑠^𝑗1=−12𝜕𝔊𝜕𝑠⁻¹𝑠^𝑖0+1𝑠^𝑗0̈+𝑠^𝑖0̇𝑠^𝑗0̇2𝜕𝔊𝜕𝑠⁻¹𝑠^𝑖0̇𝑠^𝑗0̇−12𝜕𝜕𝔉𝑠𝑠^𝑖0𝑠^𝑗0(13) (no corrections should be introduced into 𝑠^𝑖𝑠^𝑗-type products, because this would cause the appearance of terms smaller in magnitude than 𝑠₁). Here, 𝑠₀=∑^{3𝑁−6(5)𝛼=1}𝑒_𝛼cos(𝜔_𝛼𝑡+𝜃_𝛼) is a linear combination of “unit’’ vectors 𝑒_𝛼 (the superposition principle), which are also basis vectors for the construction of 𝑠₁.

A very important point should be mentioned in relation to (13). The first two terms on its right-hand side explicitly depend on atomic masses; this dependence is determined by the 𝐺 matrix. If the model of molecular motions with the 𝐹 matrix written in Cartesian coordinates is used, both these terms vanish. Of course, this cannot be balanced by any terms that appear in the expansion of the potential function into a series. It follows that the solution to the problem depends on the selected model of molecular motions.

In the model based on the use of Cartesian coordinates, atoms are “tied’’ to their equilibrium positions and move rectilinearly with respect to them. Naturally, no centrifugal effects (the second term on the right-hand side of (13), see below) can appear in such systems, because there are no bonds between atoms. This model is used as a basis of the so-called 𝑅_𝛼 structure [6] calculated in a “one-and-a-half’’ approximation (taking into account terms second-order in atomic displacements but ignoring terms of the same order of smallness that appear in the expansion of the potential energy function into a series). The model based on the use of Cartesian coordinates cannot be considered satisfactory, because the potential energy is determined by the mutual arrangement of vibrational system particles (valence interactions, angles between bonds, etc.) rather than their positions in space.

Problem (13) can be divided into three independent problems:𝐼𝑠̈^𝑗1,𝑘+Ω𝑠^𝑗1,𝑘=−12𝜕𝔊𝜕𝑠⁻¹𝑠^𝑖0𝐼𝑠̈_𝑗1,𝑐𝑠^𝑗0̈+𝑠^𝑖0̇𝑠^𝑗0̇,(14)+Ω𝑠^𝑗1,𝑐=12𝜕𝔊𝜕𝑠⁻¹𝑠^𝑖0̇𝐼𝑠̈_𝑗1,𝑎𝑠^𝑗0̇,(15)+Ω𝑠^𝑗1,𝑎=−12𝜕𝜕𝔉𝑠𝑠^𝑖0𝑠^𝑗0.(16) Clearly, 𝑠₁=𝑠_1,𝑘+𝑠_1,𝑐+𝑠_1,𝑎 is the solution to (13). Let us consider problems (14)–(16) sequentially.

3. Kinematic Problem: Equation (14)

Problems (14)–(16) are solved following similar schemes [7]. Problem (14) arises because of the dependence of kinematic coefficients (matrix 𝐺 elements) on the instantaneous configuration of the system. We have𝐼𝑠̈_1,𝑘,𝛼+Ω𝑠_1,𝑘,𝛼=12𝜔^2𝛼cos𝜔_𝛼𝑡+𝜃_𝛼^3𝑁−6(5)_𝛽=1cos𝜔_𝛽𝑡+𝜃_𝛽Δ𝔊^𝛽−1𝑒_𝛼−12𝜔^𝛼sin𝜔_𝛼𝑡+𝜃_𝛼^3𝑁−6(5)_𝛽=1𝜔_𝛽sin𝜔_𝛽+𝜃_𝛽Δ𝔊^𝛽−1𝑒_𝛼+2𝐼𝜔_0,𝛼𝜔_1,𝛼cos𝜔_𝛼𝑡+𝜃_𝛼𝑒_𝛼.(17) There should be no resonance terms containing cos(𝜔_𝛼𝑡+𝜃_𝛼) in a conservative system. For this reason, 𝜔₁=0 for all 𝛼 (the known result [7]). In (17), Δ𝔊^𝛽−1 is 𝐿^†[(𝜕𝐺⁻¹/𝜕𝑞)𝑞_𝛽]𝐿, that is, the increment of the 𝔊⁻¹ matrix caused by vibration at an 𝜔_𝛽 frequency. The 𝑒_𝛼 vector “cuts’’ the 𝛼th column, 𝔤̂^𝛽(𝛼), from the Δ𝔊^𝛽−1 matrix. Transforming the products of cosine and sine functions as cos𝛼cos𝛽=(1/2)[cos(𝛼+𝛽)+cos(𝛼−𝛽)] and sin𝛼sin𝛽=(1/2)[cos(𝛼−𝛽)−cos(𝛼+𝛽)], we obtain𝐼𝑠̈_1,𝑘,𝛼+Ω𝑠_1,𝑘,𝛼=14𝜔^{𝛼3𝑁−6(5)}_𝛽=1𝜔_𝛼+𝜔_𝛽cos𝜔_𝛼+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽𝔤̂^𝛽(𝛼)+14𝜔^{𝛼3𝑁−6(5)}_𝛽=1𝜔_𝛼−𝜔_𝛽cos𝜔_𝛼−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽𝔤̂^𝛽(𝛼).(18) The solution is sought in the form of a linear combination of the periodic functions present on the right-hand side of (18),𝑠_1,𝑘,𝛼=^3𝑁−6(5)_𝛽=1+𝑏_𝛽𝑎_𝛽cos𝜔_𝛼+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽𝑠̈_1,𝑘,𝛼cos𝜔_𝛼−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽,=−^3𝑁−6(5)_𝛽=1+𝑏_𝛽𝑎_𝛽𝜔_𝛼+𝜔_𝛽²cos𝜔_𝛼+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽𝜔_𝛼−𝜔_𝛽²cos𝜔_𝛼−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽,(19) where 𝑎 and 𝑏 are the sought vector coefficients. Substituting (19) into (18) and equating the corresponding harmonic terms, we obtain𝑎_𝛽=14𝜔^𝛼𝜔_𝛼𝑏_𝛽+𝜔_𝛽𝔓^+(𝛽)𝔤̂^𝛽(𝛼),=14𝜔^𝛼𝜔_𝛼−𝜔_𝛽𝔓^−(𝛽)𝔤̂^𝛽(𝛼),(20) where 𝔓^±(𝛽) is the diagonal matrix with the elements 𝔭^{(𝛽)𝛾𝛾}=1/[𝜔^2𝛾−(𝜔_𝛼±𝜔_𝛽)²]. Of course, there are no problems with 𝔓^(𝛽) matrices. As to Δ𝔊^𝛽−1 matrices, they are easy to obtain by numerical differentiation.

To summarize, in the “kinematic’’ approximation, we have the following:𝑠_1,𝑘,𝛼=14^3𝑁−6(5)_𝛽=1𝜔_𝛼𝜔_𝛼+𝜔_𝛽cos𝜔_𝛼+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽𝔓^+(𝛽)𝔤̂^𝛽(𝛼)+𝜔_𝛼𝜔_𝛼−𝜔_𝛽cos𝜔_𝛼−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽𝔓^−(𝛽)𝔤̂^𝛽(𝛼).(21) Of course, negative frequencies are meaningless, but cos(−𝛼)=cos𝛼. According to (21), the coefficient of the term with zero frequency cos(𝜔_𝛼−𝜔_𝛼) is zero; that is, there is no constant displacement related to the shrinkage effect. This is, however, already clear from (14), in which the coefficients of cos²𝛼 and sin²𝛼 are equal in magnitude and opposite in sign. The 𝑠_1,𝑘 values, however, make a contribution (not very significant) to vibrational amplitudes, which is important for the interpretation of electron diffraction experiments.

Although in problem (14), the 𝑞 coordinates retain their equilibrium values after averaging, because all the cosine-containing terms then vanish, the shrinkage effect for nonbonded distances can be quite substantial. Its origin is explained in Section 6.

4. “Centrifugal’’ Problem: Equation (15)

This problem (“bond-on-a-block problem’’) was casually discussed by Bartell [8]. It involves the derivatives of the kinetic energy with respect to coordinates. I call this problem centrifugal from the following considerations.

Let us consider the triatomic fragment shown in Figure 2. Let the 𝐴𝐵 bond rotate with respect to the 𝐴 point at an angular velocity 𝑞_𝜑̇. The linear velocity of the 𝐵 point is then 𝑣_𝐵=𝑟𝑞_𝜑̇, where 𝑟 is the distance between points 𝐴 and 𝐵, and the kinetic energy of the material point with mass 𝑚_𝐵 is 𝑇_𝐵=(1/2)𝑚_𝐵𝑟²𝑞^2𝜑̇. Derivative 𝜕𝑇_𝐵/𝜕𝑟 calculations yield 𝑚_𝐵𝑟𝑞^2𝜑̇, which exactly equals the 𝑓_𝑐=𝑚_𝐵𝑟𝑞^2𝜑̇ centrifugal force that acts on the material point with mass 𝑚_𝐵 which moves at an angular velocity 𝑞_𝜑̇ with respect to point 𝐴.

Figure 2 

Bond rotation about a center.

Similar considerations apply to wagging and torsional vibrations, when the distance from the axis of rotation changes (Figure 3). We then have 𝑓_𝑐=𝑚_𝐵(𝑟sin(∠𝐶𝐴𝐵)𝑞^2𝜑̇) and 𝑇_𝐵=(1/2)𝑚_𝐵(𝑟sin(∠𝐶𝐴𝐵)²𝑞^2𝜑̇). The derivative of 𝑇_𝐵 with respect to 𝑟, 𝜕𝑇_𝐵/𝜕𝑟=𝑓_𝑐sin(∠𝐶𝐴𝐵) (Figure 3, 𝑓^bond_𝑐), has the meaning of the projection of the centrifugal force onto the 𝐴𝐵 bond, which rotates about the 𝐶𝐴 axis at an angular velocity 𝑞_𝜑̇ (the sine of the ∠𝐶𝐴𝐵 angle equals the cosine of the angle between the 𝐴𝐵 bond and the normal to the axis of rotation). Clearly, an increase in the 𝐴𝐵 distance always increases the kinetic energy, and the corresponding derivatives are always positive.

Figure 3 

Bond rotation about an axis.

Valence angle changes do not influence the kinetic energy of stretching and bending vibrations, because they only rotate the corresponding velocity vectors but do not change their lengths. On the other hand, it is clear from Figure 3 that an increase in the ∠𝐶𝐴𝐵 angle decreases the kinetic energy of torsional motion if this angle is larger than 𝜋/2 (the 𝐵 point then approaches the axis of rotation) and increases it if ∠𝐶𝐴𝐵<𝜋/2 (the 𝐵 point moves from the axis). The derivative of the kinetic energy of torsional motion with respect to the angular coordinate, 𝜕𝑇_𝐵/𝜕∠𝐶𝐴𝐵=𝑓_𝑐𝑟cos(∠𝐶𝐴𝐵) (Figure 3, 𝑓^angle_𝑐), can be treated as the centrifugal force acting on the angle. This force equals the product of the length 𝑟 of the “lever” by the projection of the 𝑓_𝑐 force onto the direction normal to the 𝐴𝐵 bond. It is easy to imagine what happens to a system comprising three flexibly connected rods when the central rod is rotated (torsional motion): two end rods tend to assume the orientation normal with respect to the axis of rotation.

Similarly, the derivative of the contribution of a wagging coordinate to the kinetic energy with respect to the angle between the bonds lying on one side of the axis of rotation should be negative, and the derivative with respect to two other angles, positive. This is easy to understand considering a simple mechanical model and centrifugal forces that arise as a result of wagging vibrations: bonds that lie on the one side of the axis of rotation tend to approach each other.

For problem (15), the equation similar to (17) has the form𝐼𝑠̈_1,𝑐,𝛼+Ω𝑠_1,𝑐,𝛼=12𝜔^𝛼×sin_3𝑁−6(5)𝜔_𝛼𝑡+𝜃_𝛼_𝛽=1𝜔_𝛽sin𝜔_𝛽+𝜃_𝛽Δ𝔊^𝛽−1𝑒_𝛼(22) (the resonance term is excluded). Transforming the products of sine functions as in (18), we obtain𝐼𝑠̈_1,𝑐,𝛼+Ω𝑠_1,𝑐,𝛼=−14𝜔^{𝛼3𝑁−6(5)}_𝛽=1𝜔_𝛽cos𝜔_𝛼+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽𝔤̂^𝛽(𝛼)+14𝜔^{𝛼3𝑁−6(5)}_𝛽=1𝜔_𝛽cos𝜔_𝛼−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽𝔤̂^𝛽(𝛼).(23) Substituting (19) into the left-hand side of (23) yields𝑎_𝛽=−14𝜔^𝛼𝑏_𝛽𝜔_𝛽𝔓^+(𝛽)𝔤̂^𝛽(𝛼),=−14𝜔^𝛼𝜔_𝛽𝔓^−(𝛽)𝔤̂^𝛽(𝛼)(24) (the notation is the same as before). It follows that𝑠_1,𝑐,𝛼=14^3𝑁−6(5)_𝛽=1−𝜔_𝛼𝜔_𝛽cos𝜔_𝛼+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽𝔓^+(𝛽)𝔤̂^𝛽(𝛼)+𝜔_𝛼𝜔_𝛽cos𝜔_𝛼−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽𝔓^−(𝛽)𝔤̂^𝛽(𝛼).(25) Averaging over time and ensemble gives the contribution of the “centrifugal’’ term to the shrinkage effect,𝑠_1,𝑐=14^3𝑁−6(5)_𝛼=1𝜔^2𝛼𝔓^−(𝛼)𝔤̂^𝛼(𝛼)=14^3𝑁−6(5)_𝛼=1𝜔^2𝛼Ω⁻¹𝔤̂^𝛼(𝛼)(26) (clearly, 𝔓^−(𝛼)=Ω⁻¹). This, for instance, corresponds to elongation of bonds caused by bending vibrations (bond rotations about a center, see Figure 2).

5. Anharmonic Problem: Equation (16)

This problem is quite similar to the preceding one. Let us introduce the denotation for the table of “cubic force constants’’ 𝜕𝔉/𝜕𝑠=ℌ. ℌ is a table (three-dimensional matrix) of the third derivatives of the potential energy with respect to the coordinates. In this system, the eigenvectors of the vibrational problem of the first approximation are “unit" vectors. Let us stage-by-stage reproduce solution to problem (15):𝐼𝑠̈_1,𝑎,𝛼+Ω𝑠_1,𝑎,𝛼=−12𝜔^𝛼cos𝜔_𝛼𝑡+𝜃_𝛼^3𝑁−6(5)_𝛽=1𝜔_𝛽cos𝜔_𝑏𝑡+𝜃_𝛽ℌ𝑒_𝛼𝑒_𝛽=−14^3𝑁−6(5)_𝛽=1cos𝜔_𝛼+𝜔_𝛽+𝜃_𝛼+𝜃_𝛽ℌ𝑒_𝛼𝑒_𝛽+×ℌco𝑒s_𝛼𝜔_𝛼−𝜔_𝛽+𝜃_𝛼−𝜃_𝛽𝑒_𝛽(27) (the resonance term is excluded). Substituting (19) into the left-hand side of this equation yields𝑎_𝛽=−14𝔓^+(𝛽)𝑏_𝛽ℌ𝑒_𝛼𝑒_𝛽,=−14𝔓^−(𝛽)ℌ𝑒_𝛼𝑒_𝛽.(28) This gives𝑠_1,𝑎,𝛼=−14^3𝑁−6(5)_𝛽=1cos𝜔^𝛼+𝜔_𝛽+𝜃_𝛼+𝜃_𝛽𝔓^+(𝛽)+×𝔓cos_−(𝛽)𝜔_𝛼−𝜔_𝛽+𝜃_𝛼−𝜃_𝛽ℌ𝑒_𝛼𝑒_𝛽.(29) We find that the contribution of the anharmonic component to the shrinkage effect, which is determined by the term containing cos[(𝜔_𝛼−𝜔_𝛼)+(𝜃_𝛼−𝜃_𝛼)](𝛽=𝛼), can be written as𝑠_1,𝑎,𝛼=−41Ω⁻¹ℌ𝑒_𝛼𝑒_𝛽.(30) The 𝑒_𝛼𝑒_𝛽 vectors “cut’’ a column with the 𝔥̂_𝛼𝛽… elements from the ℌ matrix, and calculations present no difficulties.

Note that (21), (25), and (27) present the first members of the series describing atomic displacements from equilibrium positions. The convergence of these series depends, in particular, on the lengths of the 𝔤̂ and 𝔥̂ vectors, which are, in turn, determined by the lengths of “unit’’ vectors 𝑒 (𝑒_𝛼=[(ℎ/4𝜋²𝑐𝜈_𝛼)coth(ℎ𝑐𝜈_𝛼/2𝑘𝑇)]^1/2). These vectors become longer than one starting from 𝜔≈84.5cm⁻¹ at 300 K, and their length increases as the frequency decreases. The convergence of the series can then be fairly slow, and, in the presence of low frequencies, atomic displacements and, therefore, amplitudes are calculated very approximately. This does not relate to shrinkage corrections, because only odd series terms contribute to them, and we have every reason to hope that fifth-order contributions will nevertheless be much smaller than third-order ones.

To summarize, we obtained the following result:𝑠_𝛼=+𝑠1_0,𝛼(𝑡)4^3𝑁−6(5)_𝛽=1𝜔^𝛼𝜔_𝛼×cos𝜔+_𝛼𝜔_𝛽−𝜔_𝛼𝜔_𝛽×𝔓_+(𝛽)+𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽×𝔓_+(𝛽)𝔤̂^𝛽(𝛼)−cos𝜔_𝛼+𝜔_𝛽+𝜃_𝛼+𝜃_𝛽×cos𝔥̂𝜔^𝛽(𝛼)_𝛼+𝜔_𝛼𝜔_𝛼−𝜔_𝛽+𝜔_𝛼𝜔_𝛽−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽𝔓^−(𝛽)𝔤̂^𝛽(𝛼)−cos𝜔_𝛼−𝜔_𝛽+𝜃_𝛼−𝜃_𝛽𝔓^−(𝛽)𝔥̂^𝛽(𝛼)(31) for the vector of displacements with respect to the equilibrium configuration and⟨𝑠_𝛼⟩=14Ω⁻¹𝜔^2𝛼𝔤̂^𝛼(𝛼)−𝔥̂^𝛼(𝛼)(32) for the contribution of the vibration with the 𝜔_𝛼 frequency to the shrinkage effect. Note that 𝔤̂^𝛽(𝛼)=𝔤̂^𝛼(𝛽) (these vectors only differ in the order of differentiation), and, naturally, 𝔥̂^𝛽(𝛼)=𝔥̂^𝛼(𝛽). For this reason, the summation in 𝛽 can be performed from 𝛽=1 to 𝛽=𝛼.

The transition to vectors in internal coordinates 𝑞 is trivial,𝑙_𝛼=𝑙_0,𝛼+𝐿𝑠_1,𝛼,(33) where 𝑠_1,𝛼 is the sum of small terms in (31) (the sum of corrections to eigenvectors obtained by solving the kinematic, centrifugal, and anharmonic problems), and⟨𝑙_𝛼⟩=𝐿𝑠_1,𝛼.(34) Here, 𝐿 is the matrix defined by (11).

The solution becomes meaningless if (𝜔_𝛼±𝜔_𝛽)²−𝜔^2𝛾=0, for instance, if 𝜔_𝛽=2𝜔_𝛼 (the matrix with the [(𝜔_𝛼±𝜔_𝛽)²−𝜔^2𝛾] elements then has no inverse; it follows, that the 𝔓_± matrix cannot be constructed). This is a situation of the type of parametric resonance known in classical mechanics [7], when, for instance, the mass of a particle changes at twice the vibrational frequency. An analysis of such situations is outside the scope of first-order perturbation theory used in this work.

6. Calculations of Corrections for the Shrinkage Effect and Amplitudes of Changes in Internuclear Distances

Above, the algorithm was described for calculations of shrinkage effect corrections and vibrational amplitudes for internal (generalized) coordinates 𝑞. For the transition to internuclear distances and the spatial configuration of the system, the results should be transformed into Cartesian coordinates. This transformation is performed as𝑥_𝛼=𝐵⁻¹𝑙_𝛼cos𝜔_𝛼𝑡+𝜃_𝛼,(35) or, according to (33),𝑥_𝛼=𝐵⁻¹𝑙_0,𝛼cos𝜔_𝛼𝑡+𝜃_𝛼+𝐵⁻¹𝐿𝑠_1,𝛼.(36) Because of the smallness of the second term, we can assume 𝐵⁻¹=𝐵(0)⁻¹ in it and exclude it from consideration in this section. Since 𝐵⁻¹=𝐵(0)⁻¹+(1/2)∑^{3𝑁−6(5)𝛽=1}(𝜕𝐵⁻¹/𝜕𝑥)𝑥_𝛽cos(𝜔_𝛽𝑡+𝜃_𝛽), we have𝑥^𝛼tr=12^3𝑁−6(5)_𝛽=1Δ𝐵^𝛽−1𝑙_0,𝛼cos𝜔_𝛼𝑡+𝜃_𝛼cos𝜔_𝛽𝑡+𝜃_𝛽,(37) where 𝑥^tr is the displacement that appears because of nonlinearity of the transition between internal and Cartesian coordinates and Δ𝐵^𝛽−1 is the increment of the 𝐵⁻¹ matrix when the configuration of the system changes by the 𝑥_𝛽 vector. In this sum, only the term with cos²(𝜔_𝛼𝑡+𝜃_𝛼) does not vanish in averaging. For this reason,𝑥^𝛼tr=12𝜕𝐵𝜕𝑥⁻¹𝑥^𝛼𝑙_0,𝛼cos²𝜔_𝛼𝑡+𝜃_𝛼=41𝜎^𝛼Δ𝐵^𝛼−1𝑙_0,𝛼,(38) where 𝜎_𝛼=(ℎ/4𝜋²𝑐𝜈_𝛼)coth(ℎ𝑐𝜈_𝛼/2𝑘𝑇) [6]. The ⟨𝑥^𝛼tr⟩ value can be obtained without comparatively laborious calculations of the Δ𝐵^𝛼−1 matrices. For this purpose, let us rewrite (38) as follows. Since 𝐵⁻¹=(𝐵(0)+(1/2)Δ𝐵)⁻¹=(𝐼+(1/2)𝐵(0)⁻¹Δ𝐵)⁻¹𝐵(0)⁻¹≈(𝐼−(1/2)𝐵(0)⁻¹Δ𝐵)𝐵(0)⁻¹, we have Δ𝐵⁻¹≈𝐵(0)⁻¹Δ𝐵𝐵(0)⁻¹. On the other hand, 𝐵(0)⁻¹𝑙_0,𝛼=𝑥_0,𝛼. Therefore,𝑥^𝛼tr≈41𝜎^𝛼𝐵(0)⁻¹Δ𝐵_𝛼𝐵(0)⁻¹𝑙_0,𝛼=14𝜎^𝛼𝐵(0)⁻¹Δ𝐵_𝛼𝑥_0,𝛼.(39) Here, 𝑥_0,𝛼 is the eigenvector in Cartesian coordinates obtained by solving the vibrational problem in the first approximation. Let us calculate the 𝑙^+0,𝛼 and 𝑙^−0,𝛼 vectors for the 𝑅₀+𝑥_0,𝛼 and 𝑅₀−𝑥_0,𝛼 configurations, respectively (𝑅₀ is the vector of the Cartesian coordinates of atoms in the equilibrium configuration). Clearly,𝑙^+0,𝛼=𝐵(0)+12Δ𝐵^𝛼𝑥^0,𝛼𝑙_−0,𝛼,=𝐵(0)−12Δ𝐵^𝛼−𝑥^0,𝛼.(40) Summing these equations yields𝑙^+0,𝛼+𝑙^−0,𝛼=Δ𝐵_𝛼𝑥_0,𝛼.(41) Substituting the result into (39) yields𝑥^𝛼tr=41𝜎^𝛼𝐵(0)⁻¹𝑙^+𝛼+𝑙^−𝛼.(42) On the other hand, the 𝑥^+𝛼 and 𝑥^−𝛼 vectors satisfying the conditions 𝑙^+𝛼=𝐵𝑥^+𝛼 and 𝑙^−𝛼=𝐵𝑥^−𝛼 can be calculated by the method of successive approximations by fitting the 𝐵(𝑥) matrix and 𝑥 vector components to obtain 𝑞_𝛼=𝐵(𝑥)𝑥_𝛼. Clearly, ⟨𝑥_𝛼⟩=(𝑥^+𝛼+𝑥^−𝛼)/2. This method gives somewhat more accurate results, especially when vibrational amplitudes are indeed large. It is implemented in the Shrink09 program [9].

Clearly, the procedure described in this section does not change 𝑙_𝛼 eigenvector components in generalized coordinates and, therefore, mean bond lengths, valence angles, and so forth. The refinement of the 𝑥_𝛼 Cartesian displacements of atoms and their mean values ⟨𝑥⟩, however, leads to substantial shrinkage effect corrections for distances between nonbonded atoms measured by diffraction methods.

Combining (32) and (38), we obtain⟨𝑥_𝛼⟩=14𝐵(0)⁻¹𝜎^𝛼Δ𝐵^𝛼−1𝑙_0,𝛼+𝐿Ω𝜔^2𝛼𝔤̂^𝛼(𝛼)−𝔥̂^𝛼(𝛼)(43) (frequency factors 𝜎 are already contained in the 𝔤̂ and 𝔥̂ vectors by virtue of the definition of eigenvectors in Section 2) and⟨𝑥⟩=^3𝑁−6(5)_𝛼=1⟨𝑥_𝛼⟩(44) (clearly, we must add Hedberg's corrections for centrifugal distortions caused by rotations of a molecule as a whole to this result [10]). This result allows us to determine four parameters important for structural studies: the mean configuration of the molecule necessary for the introduction of corrections into experimentally observed rotational constants (microwave experiments), shrinkage effect values for the experimental set of internuclear distances (diffraction experiments), mean-square amplitudes of changes in internuclear distances (gas-phase diffraction), and asymmetry parameters (skewness; gas-phase diffraction). Let us denote the Cartesian coordinates of the 𝑖th atom by 𝑋_𝑖, 𝑌_𝑖, and 𝑍_𝑖, and the corresponding components of the ⟨𝑥⟩ and ⟨𝑥_𝛼⟩ vectors, by 𝜉_𝑖, 𝜐_𝑖, 𝜁_𝑖 and 𝜉_𝛼,𝑖, 𝜐_𝛼,𝑖, 𝜁_𝛼,𝑖. The 𝑅_av vector with the components𝑋_𝑖+^3𝑁−6(5)_𝛼=1𝑌_𝑖𝜉_𝛼,𝑖,+^3𝑁−6(5)_𝛼=1𝑍_𝑖𝜐_𝛼,𝑖,+^3𝑁−6(5)_𝛼=1𝜁_𝛼,𝑖(45) then determines the mean configuration of the molecule and, therefore, the mean values of its rotational constants.

On the other hand, for the 𝑅_𝑖𝑗 distance between atoms 𝑖 and 𝑗, we haveΔ𝑅_𝑖𝑗=∑^3𝑁−6(5)_𝛼=1𝑋^𝑖−𝑋_𝑗+𝜉_𝛼,𝑖−𝜉_𝛼,𝑗²+𝑌_𝑖−𝑌_𝑗+𝜐_𝛼,𝑖−𝜐_𝛼,𝑗²+𝑍_𝑖−𝑍_𝑗+𝜁_𝛼,𝑖−𝜁_𝛼,𝑗²^1/2−𝑅_0,𝑖𝑗,(46) where 𝑅_0,𝑖𝑗 is the distance between atoms 𝑖 and 𝑗 in the equilibrium configuration. The Δ𝑅_𝑖𝑗 values determine the shrinkage effect for internuclear distances; in diffraction experiments, we measure the 𝑅_0,𝑖𝑗+Δ𝑅_𝑖𝑗 distances, and, for the transition to the equilibrium configuration, the corresponding corrections should be introduced into these values.

Two other important parameters are mean square amplitudes of changes in internuclear distances and skewness. Let us rewrite (37) in the form𝑥^𝛼tr=14^3𝑁−6(5)_𝛽=1Δ𝐵^𝛽−1𝑙_0,𝛼+cocsos𝜔_𝛼𝜔+_𝛼𝜔_𝛽𝑡+𝜃_𝛼+𝜃_𝛽−𝜔_𝛽𝑡+𝜃_𝛼−𝜃_𝛽.(47) To calculate the mean square amplitudes of the displacement of atoms from equilibrium positions, it is necessary: (i) to perform the transition to the Cartesian coordinates in (31) through multiplying by the 𝐵(0)⁻¹𝐿 matrix (the 𝑠_0,𝛼(𝑡) value then transforms into 𝑥_0,𝛼cos(𝜔_𝛼𝑡+𝜃_𝛼)), (ii) sum (31) and (47), (iii) raise the result to a power of two, and (iv) perform averaging over time and phases (the latter for degenerate vibrations). In averaging, only cos²𝛼-type terms do not vanish (give 1/2) whereas all the cross terms (of the cos𝛼cos𝛽 type) reduce to zero. For this reason, it is sufficient to calculate the coefficients of the harmonic terms in (31) (written in Cartesian coordinates) and (47) to obtain the ⟨𝑥^2𝛼⟩ vector,𝑥^2𝛼=𝜎_𝛼Sq𝑥_0,𝛼+12^3𝑁−6(5)_𝛽=1Sq𝑎^+𝛽+Sq𝑎^−𝛽(48) (the 1/2 coefficient appears in the averaging of squared cosine functions over time and phases). Here the Sq symbol denotes the multiplication by the diagonal matrix whose components are the components of the vector following this symbol (squaring of vector components), and the 𝑎^+𝛽 and 𝑎^−𝛽 vectors are𝑎^+𝛽=14𝜔^𝛼𝜔_𝛼+𝜔_𝛽−𝜔_𝛼𝜔_𝛽𝐵(0)⁻¹𝐿𝔓^+(𝛽)𝔤̂^𝛽(𝛼)+𝐵(0)⁻¹𝐿𝔓^+(𝛽)𝔥̂^𝛽(𝛼)+12√𝑎_−𝛽𝜎_𝛼𝜎_𝛽Δ𝐵^𝛽−1𝑙_0,𝛼,=14𝜔^𝛼𝜔_𝛼−𝜔_𝛽+𝜔_𝛼𝜔_𝛽𝐵(0)⁻¹𝐿𝔓^−(𝛽)𝔤̂^𝛽(𝛼)+𝐵(0)⁻¹𝐿𝔓^−(𝛽)𝔥̂^𝛽(𝛼)+12√𝜎_𝛼𝜎_𝛽Δ𝐵^𝛽−1𝑙_0,𝛼.(49) The mean-square amplitudes of changes in internuclear distances ⟨𝑥^2𝑖𝑗⟩ are given by𝑥^2𝑖𝑗=12^3𝑁−6(5)_𝛼=1𝑥_{𝛼,𝑖,𝑥}−𝑥_{𝛼,𝑗,𝑥}²+𝑥_{𝛼,𝑖,𝑦}−𝑥_{𝛼,𝑗,𝑦}²+𝑥_{𝛼,𝑖,𝑧}−𝑥_{𝛼,𝑗,𝑧}²,(50) where 𝑥_{𝛼,𝑖,𝑥}, 𝑥_{𝛼,𝑖,𝑦}, and 𝑥_{𝛼,𝑖,𝑧} are the components of the 𝑥_0,𝛼+∑^{3𝑁−6(5)𝛽=1}(𝑎^+𝛽+𝑎^−𝛽) vector corresponding to the displacements of the 𝑖th atom along the 𝑥, 𝑦, and 𝑧 axes; again, ⟨cos²(𝑎𝑡+𝑏)⟩ gives the coefficient 1/2. To pass to the central moment from the moment about zero, we must subtract ⟨𝑅_{𝛼,𝑖𝑗}⟩² from ⟨𝑥^{2𝛼,𝑖𝑗}⟩. This is equivalent to excluding terms with cos(𝛼−𝛼) from amplitude calculations.

If the Δ𝑅_{𝛼,𝑖𝑗} and ⟨𝑥^{2𝛼,𝑖𝑗}⟩ values are known, asymmetry parameter ⟨𝑥^3𝑖𝑗⟩/⟨𝑥^2𝑖𝑗⟩^3/2 calculations present no difficulties. Indeed,𝑥^3𝑖𝑗=^3𝑁−6(5)_𝛼=1Δ𝑅^{3𝛼,𝑖𝑗}.(51)

7. Calculations of the Third Derivatives of Potential Energy with respect to Internal Coordinates

The above equations are fairly cumbersome, but the corresponding algorithm is easy to write. Certain difficulties arise with cubic force constants, which are formed by quantum-mechanical programs, first, by a numerical method (the method of finite differences) and, secondly, in Cartesian coordinates. The first circumstance requires estimating possible calculation errors.

7.1. Errors in Cubic Force Constant Calculations

The set of cubic force constants found in the results of quantum-mechanical calculations can conveniently be represented in the form of a 3𝑁×3𝑁×3𝑁 cubic table (𝑁 is the number of atoms in the molecule). Let us consider how the 𝑘th “cut’’ of this table is obtained (Figure 4).

Figure 4 

Cubic table of third-order potential energy derivatives with respect to Cartesian coordinates.

This table can be denoted by the 𝐻^𝐶 (Cartesian) symbol. The cut with the elements ℎ^{𝐶𝑘,𝑖𝑗} (the 𝑘th cut) is formed as follows. First, the 𝑘th Cartesian coordinate of the molecule 𝑥_𝑘 changes by a small value 𝛿, and the Hessian 𝐻 of the molecule is calculated in this intentionally distorted nonequilibrium configuration,𝐻^+𝑘=𝐻𝑥₁,…,𝑥_𝑘+𝛿,…,𝑥_3𝑁.(52) The same is made after the subtraction of 𝛿 from 𝑥_𝑘,𝐻^−𝑘=𝐻𝑥₁,…,𝑥_𝑘−𝛿,…,𝑥_3𝑁.(53) Next, we subtract the second result from the first one and divide the difference by 2𝛿. As a result, we obtain the 𝑘th cut of the table of cubic constants,𝐻^𝐶𝑘=𝐻^+𝑘−𝐻^−𝑘2𝛿.(54) It follows that the ℎ^{𝐶𝑘,𝑖𝑗} element of the table of cubic force constants is calculated asℎ^{𝐶𝑘,𝑖𝑗}=ℎ^+𝑖𝑗−ℎ^−𝑖𝑗2𝛿.(55) The equation for calculating the error in the ℎ^{𝐶𝑘,𝑖𝑗} element obtained this way is well known [11], Δℎ^{𝐶𝑘,𝑖𝑗}=16𝛿²||||𝜕³ℎ_𝑖𝑗𝜕𝑥^3𝑘||||_𝜇+𝜀^(𝑖𝑗)𝛿,(56) where 𝜇∈[𝑥_𝑘−𝛿,𝑥_𝑘+𝛿]. Here, the first term is the so-called truncation error. Its origin becomes clear if the difference ℎ^+𝑖𝑗(𝑥_𝑘+𝛿)−ℎ^−𝑖𝑗(𝑥_𝑘−𝛿) is expanded in powers of 𝛿. The second term is determined by the error in the ℎ_𝑖𝑗 values themselves.

It follows from this equation that, first, the error in ℎ^{𝐶𝑘,𝑖𝑗} is not a monotonic function of 𝛿: the smaller 𝛿, the smaller the truncation error, but the larger the contribution of the error in ℎ_𝑖𝑗. The 𝛿 values can, of course, be optimized, and the corresponding algorithms are well known [10], but the use of 𝛿 of its own for each ℎ^{𝐶𝑘,𝑖𝑗} element would make the problem of calculations of cubic force constants quite unrealistic. This means that calculations are performed with 𝛿 values that are far from optimum.

Secondly, errors in ℎ^{𝐶𝑘,𝑖𝑗} elements that differ only by the order of indices are different. Indeed,Δℎ^{𝐶𝑘,𝑖𝑗}=16𝛿²||||𝜕³ℎ_𝑖𝑗𝜕𝑥^3𝑘||||_𝜇+𝜀^(𝑖𝑗)Δℎ_{𝐶𝑗,𝑘𝑖}𝛿,=16𝛿²||||𝜕³ℎ_𝑘𝑖𝜕𝑥^3𝑘||||_𝜇+𝜀^(𝑘𝑖)𝛿.(57) It follows that we cannot expect to obtain the theoretical equality ℎ^{𝐶𝑖𝑗𝑘}=ℎ^{𝐶𝑖𝑘𝑗}=ℎ^{𝐶𝑗𝑖𝑘}=ℎ^{𝐶𝑗𝑘𝑖}=ℎ^{𝐶𝑘𝑖𝑗}=ℎ^{𝐶𝑘𝑗𝑖} with the use of finite-difference schemes.

The following procedure for estimating calculation errors can be suggested. Let index 𝛼 run over all the permutations of the 𝑖,𝑗,𝑘 indices. The mean ℎ^𝐶𝛼 value is thenℎ^𝐶𝛼=16_𝛼ℎ^𝐶𝛼.(58) Let us putΔ_𝛼=16_𝛼ℎ^𝐶𝛼−ℎ^𝐶𝛼²^1/2,Δ^𝛼=maxℎ^𝐶𝛼−minℎ^𝐶𝛼2.(59) In other words, Δ is the root-mean-square deviation for the given 𝛼 (for the given set of 𝑖, 𝑗, and 𝑘 indices), and Δ^ is the half-spread of the ℎ^𝐶𝛼 values. We putΔ=maxΔ_𝛼,Δ^=maxΔ^𝛼,(60) where maxima are sought over all 𝛼 (over all different sets of 𝑖, 𝑗, and 𝑘 indices). Δ can be suggested as an estimate of errors in off-diagonal 𝐻^𝐶 matrix elements, and Δ′, as a (conventional) estimate of errors in ℎ^{𝐶𝑖𝑖𝑖}. The Gaussian program cannot be used to obtain the required information. The corresponding calculations were therefore performed using the Gamess [12] package and a special program external with respect to Gamess (the Gamess package is not intended for cubic constant calculations). Calculations for three molecules gave Δ=0.825737×10⁻⁴ au and Δ′=0.984625×10⁻⁴ au (maximum cubic constant values for these molecules were on the order of 1.5–2 au).

These results were obtained using options that imposed severe requirements on the accuracy of calculations. The use of default options increased errors by an order of magnitude. Nevertheless, they remained within the limits acceptable for practical purposes.

However, note that, in high-accuracy calculations, errors increased as the number of atoms in molecules grew, which seems to be a natural result. In calculations with options “by default,” the largest errors were obtained for the smallest molecule containing chlorine (the other molecules contained C, H, N, O, and Si). It may well be that, for molecules with atoms of different chemical natures or a large number of atoms, errors can be outside admissible limits.

Note that cubic constant tables in outputs of the Gaussian program often contain inexplicable misprints. It is much safer to extract the required information from the last (inconveniently formatted) table printed out.

7.2. The Transition from Cartesian to Internal Coordinates

The three-dimensional matrix of cubic constants in Cartesian coordinates is not a tensor [13], and cannot be transformed linearly, as𝜕𝜕𝐹𝑞_𝑖𝑗𝑘=𝐻^{𝐶𝑖𝑗𝑘}𝐵^𝑖𝑗𝐵^𝑗𝑘𝐵^𝑘𝑗⁻¹,(61) where (𝜕𝐹/𝜕𝑞)_𝑖𝑗𝑘 is the three-dimensional matrix of the third derivatives of potential energy with respect to internal coordinates, 𝐻_𝑖𝑗𝑘 is, as before, the three-dimensional matrix of the third derivatives of potential energy with respect to Cartesian coordinates, and 𝐵^𝑗𝑖 is the transformation matrix between internal and Cartesian coordinates, 𝑞^𝑗=𝐵^𝑗𝑖𝑥^𝑖 (the matrix of second derivatives can be transformed this way only because the corresponding gradient vector is zero at equilibrium).

Let us consider the problem in more detail. As before, the values related to the {𝑥₁,…,𝑥_𝑘+𝛿,…,𝑥_3𝑁} configuration (see the preceding subsection) will be labeled by the superscript “+,’’ and the values related to the 𝑥_𝑘−𝛿 configuration, by the superscript “−.’’ Let us introduce the denotations𝐵^(𝑘)+=𝐵₀+Δ_𝑘,𝐵^(𝑘)−=𝐵₀−Δ_𝑘,(62) where 𝐵₀ is the 𝐵 matrix calculated for the equilibrium configuration and Δ_𝑘=(𝜕𝐵/𝜕𝑥_𝑘)𝛿.

Calculations of the 𝑘th cut of the 𝐻^𝐶 matrix begin with calculations of the gradient vector for the given configuration and then its derivatives with respect to the atomic coordinates,𝜕𝑈𝜕𝑥=𝐵^(𝑘)+^†𝜕𝑈𝜕𝜕𝑞₂,𝑈𝜕𝑥²=𝐵^(𝑘)+^†𝜕²𝑈𝜕𝑞²𝐵^(𝑘)++𝜕𝜕𝑥𝐵^(𝑘)+^†𝜕𝑈𝜕𝑞(63) (because 𝜕/𝜕𝑥=(𝜕/𝜕𝑞)(𝜕𝑞/𝜕𝑥)). Substituting this result and a similar equation for the configuration with 𝑥_𝑘−𝛿 into (54) yields𝐻^+𝑘=𝐵^†0𝐹⁺𝐵₀+𝐵^†0𝐹⁺Δ_𝑘+Δ^†𝑘𝐹⁺𝐵₀+𝜕Δ^†𝑘𝜕𝑥⁺𝜕𝜕𝑈𝑞⁺𝐻_−𝑘,=𝐵^†0𝐹⁻𝐵₀−𝐵^†0𝐹⁻Δ_𝑘−Δ^†𝑘𝐹⁺𝐵₀−𝜕Δ^†𝑘𝜕𝑥⁻𝜕𝜕𝑈𝑞⁻,(64) where 𝐹 is the matrix of second-order force constants in internal coordinates. It follows that the numerator in (54) takes the form𝐻^+𝑘−𝐻^−𝑘=𝐵^†0𝐹⁺−𝐹⁻𝐵₀+2𝐵^†0𝐹Δ_𝑘+2Δ^†𝑘𝐹𝐵₀+𝜕Δ^†𝑘𝜕𝑥𝜕𝜕𝑞𝑈⁺−𝜕𝑈𝜕𝑞⁻(65) (indeed, 𝐹⁺+𝐹⁻=2𝐹 and (𝜕Δ_𝑘/𝜕𝑥)⁻=−(𝜕Δ_𝑘/𝜕𝑥)⁺). We therefore have𝐵^†0𝐹⁺−𝐹⁻𝐵₀2𝛿=𝐻^+𝑘−𝐻^−𝑘2𝛿−𝛿1𝐵^†0𝐹𝜕𝐵𝜕𝑥_𝑘+𝜕𝜕𝑥𝐵_𝑘^†𝐹𝐵₀+𝜕Δ^𝑘𝜕𝑥𝜕²𝑈𝜕𝑞𝜕𝑥_𝑘.(66) Here, 𝜕²𝑈/𝜕𝑞𝜕𝑥_𝑘=𝐵(0)⁻¹(𝜕²𝑈/𝜕𝑥𝜕𝑥_𝑘); that is, this vector is obtained from the 𝑘th column of the 𝐻 matrix, and the other matrices are calculated by the finite-difference method. The right-hand side of this equation gives the 𝑘th cut of some 𝐻^𝐶 matrix, which can now safely be transformed into the matrix of third derivatives of potential energy with respect to internal coordinates. Let us denote this matrix by the 𝐻^𝐼 (internal) symbol,𝐻^{𝐼𝑖𝑗𝑘}=𝐻^{𝐶𝑖𝑗𝑘}𝐵^𝑖𝑗𝐵^𝑗𝑘𝐵^𝑘𝑗⁻¹.(67) Above, the algorithm of calculations is described. It corresponds to the following analytic equations:𝜕𝑈𝜕𝑥^𝑘=𝐵^𝑙𝑘𝜕𝑈𝜕₂𝜕𝑞^𝑙,𝑈𝜕𝑥^𝑘𝜕𝑥^𝑙=𝜕²𝑈𝜕𝑞^𝑘𝜕𝑞^𝑙𝐵^𝑘𝑙𝐵^𝑙𝑘+𝜕𝐵^𝑙𝑘𝜕₃𝜕𝑥^𝑙𝜕𝜕𝑈𝑞^𝑙,𝑈𝜕𝑥^𝑘𝜕𝑥^𝑙𝜕𝑥^𝑚=𝜕³𝑈𝜕𝑞^𝑘𝜕𝑞^𝑙𝜕𝑞^𝑚𝐵^𝑙𝑘𝐵^𝑘𝑚𝐵^𝑚𝑙+𝜕²𝑈×𝐵_𝑘𝑙𝜕𝑞^𝑘𝜕𝑞^𝑙𝜕𝜕𝐵𝑥^𝑙𝑘𝑚+𝜕²𝑈𝜕𝑞^𝑘𝜕𝑞^𝑙𝜕𝜕𝐵𝑥^𝑘𝑙𝑚𝐵^𝑙𝑘+𝜕²𝑈𝜕𝑞^𝑙𝜕𝑞^𝑚𝐵^𝑘𝑚𝜕𝜕𝐵𝑥^𝑙𝑘𝑙(68) (the last term differs from the preceding ones in the order of indices),𝜕³𝑈𝜕𝑞^𝑘𝜕𝑞^𝑙𝜕𝑞^𝑚=𝜕³𝑈𝜕𝑥^𝑘𝜕𝑥^𝑙𝜕𝑥^𝑚−𝜕²𝑈×𝐵_𝑘𝑙𝜕𝑞^𝑘𝜕𝑞^𝑙𝜕𝜕𝐵𝑥^𝑙𝑘𝑚−𝜕²𝑈𝜕𝑞^𝑘𝜕𝑞^𝑙𝜕𝜕𝐵𝑥^𝑘𝑙𝑚𝐵^𝑙𝑘−𝜕²𝑈×𝐵𝜕𝑞^𝑙𝜕𝑞^𝑚𝐵^𝑘𝑚𝜕𝜕𝐵𝑥^𝑙𝑘𝑙_𝑙𝑘𝐵^𝑘𝑚𝐵^𝑚𝑙⁻¹.(69)

8. Scaling of Cubic Force Constants

A popular trend of recent years is to scale quantum-mechanical force fields to approximate the calculated normal vibration frequencies to the experimental values. Scaling is usually performed for so-called pseudosymmetry coordinates (for certain linear combinations of 𝑞 coordinates) suggested by Pulay et al. [14]. I prefer to perform scaling which leaves quantum-mechanical eigenvectors unchanged (i.e., constants for the linear combinations of internal coordinates corresponding to 𝑙_𝛼 eigenvectors are scaled).

In any event, the question of how cubic force constants should be scaled arises. Using the notation introduced in Section 7, we can write𝐻^𝐼𝑘=𝐹^+𝑘−𝐹^−𝑘2𝛿,ℎ^{𝐼𝑖𝑗𝑘}=𝑓^+𝑖𝑗−𝑓^−𝑖𝑗2𝛿,(70) where 𝛿 is the shift along the 𝑘th internal coordinate, 𝑘th coordinate according to Pulay, or 𝑘th eigenvector. If the transformation𝐹^scaled=𝐷^1/2𝐹^theor𝐷^1/2(71) (𝐷 is a diagonal matrix) scales the force field in the equilibrium configuration, it is easy to accept that the same 𝐷 matrix should scale the force fields of configurations slightly changed with respect to equilibrium, that is,𝑓^+𝑖𝑗−𝑓^−𝑖𝑗^scaled=𝑑_𝑖𝑑_𝑗^1/2𝑓^+𝑖𝑗−𝑓^−𝑖𝑗^theor.(72) On the other hand, force field scaling is equivalent to the scaling of coordinates. Indeed, the transition from internal to Cartesian coordinates is performed for a scaled force field as𝐵^†𝐷^1/2𝐹^theor𝐷^1/2𝐵=𝐵^†𝐷^1/2𝐹^theor𝐷^1/2𝐵,(73) that is, 𝑞^theor_𝑖=𝑞^scaled/𝑑^1/2. This means that the 2𝛿 denominator (shift along the 𝑘th coordinate) should be divided by 𝑑^𝑘1/2,ℎ^𝐼,scaled_𝑖𝑗𝑘=𝑑_𝑖𝑑_𝑗^1/2𝑓^+𝑖𝑗−𝑓^−𝑖𝑗^theor2𝛿/𝑑^𝑘1/2=𝑑_𝑖𝑑_𝑗𝑑_𝑘^1/2ℎ^𝐼,theor_𝑖𝑗𝑘.(74) The equation ℎ^scaled_𝑖𝑗𝑘=(𝑑_𝑖𝑑_𝑗𝑑_𝑘)^1/3ℎ^theor_𝑖𝑗𝑘 suggested in [15] without any justification is, of course, absolutely incorrect. Conversely, in [16–18], the equation given above was used.

The scaling of third potential energy derivatives is an important problem. The point is that their high-level calculations take too much time. On the other hand, scaled cubic constants determined in low-level calculations give results almost indistinguishable from those obtained using high-level calculations (see Table 4 in [19]).

9. The Problem of Low Frequencies

The experimental values of low frequencies are often inaccessible, which causes difficulties in the scaling of theoretical force fields. The simplest way out is to scan the low frequency changing it with a certain step, but leaving the corresponding eigenvector unchanged. These calculations do not take much time. We can then select the value that best describes the diffraction experiment. It may well be that this is the only method for experimentally determining the position of low frequencies if they cannot be extracted from vibrational of vibronic spectra.

10. Internal Rotations and Similar Problems

The situation is possible when internal rotation or inversion become free at a fairly low excitation level (when torsional or inversion vibrational frequencies are very low). In such a situation, parameters for treatment of diffraction data are sometimes determined by calculating the “minimum energy path,’’ that is, quantum-mechanical calculations are performed with the optimization of the geometries that arise in scanning the system along, for instance, the torsional coordinate [20]. It may well be (almost inevitable) that the “torsional vibration’’ along the minimum energy path will then include coordinates inconsistent with it by symmetry. For instance, in the cited work, the nitroethane molecule (𝐶_𝑠 symmetry) was considered. Torsional vibration of the NO₂ group in this molecule transforms under the 𝐴^ representation of the 𝐶_𝑠 group. However, the minimum energy path calculated by the authors included changes in the C–C and C–N distances and other coordinates of 𝐴^ symmetry. This is hardly possible. No matter what order of perturbation theory is used, because potential function symmetry always coincides with geometric configuration symmetry of the vibrational system, matrix elements between coordinates with different symmetries are always zero.

The approach under consideration is unsatisfactory for the following reasons.(1)The resulting system of vibrational motions violates selection rules.(2)The requirement that one of the coordinates (for instance, torsional) change in the direction of increasing energy, and the other coordinates (including the components of the same eigenvector), in the direction of decreasing it is physically ungrounded. Generally, the minimum energy path for any coordinate is the absence of any vibrations.(3)The approach based on the search for a minimum energy path ignores the kinetic component, in particular, contributions to the kinetic energy from changes in the coordinates that minimize potential energy. It should be borne in mind that we deal with molecular vibrations rather than equilibrium configurational transformations, that is, with a system for which dynamic effects cannot be ignored.

Let us turn to some obvious examples. Let us, for instance, consider antisymmetric stretching vibration of an 𝐴𝐵₂ triatomic molecule. Clearly, the ∠𝐵𝐴𝐵 angle should change along the corresponding minimum energy path. But the angular coordinate transforms under the 𝐴₁ representation, whereas the antisymmetric stretching coordinate, under the 𝐵₁ (or 𝐵₂) representation. In addition, it is impossible to make the angle vibrate at the antisymmetric stretching vibration frequency.

One more example is a linear CO₂-type molecule. The minimum energy path along the bending coordinate will necessarily include changes in bond lengths. We again have inconsistency by symmetry. In addition, bonds characterized by force constants that are much larger than bending vibration constants cannot vibrate at the corresponding frequency.

In my view, it is more reasonable to approach the problem as follows. Let the potential energy curve along the eigenvector corresponding to the vibration under consideration have the form shown in Figure 5.

Figure 5 

Potential energy curve along a torsional coordinate.

The horizontal lines in Figure 5 are vibrational levels, and the potential well contains five of them. It is always possible to solve the vibrational problem for five levels. We can determine their populations and select a coefficient for recalculating the corresponding eigenvector components into the Cartesian displacements of atoms. Such calculations actually divide the system into two subsystems with known populations. For instance, let us consider the nitroethane molecule (Figure 6), which has a very low torsional frequency corresponding to NO₂ group rotations about the N–C bond, 26.4 cm⁻¹ (the B3LYP/aug-cc-pVTZ data).

Figure 6 

Nitroethane molecule.

The potential well along the torsional coordinate contains ten vibrational levels. At 325 K (the temperature of electron diffraction measurements), calculations taking into account level populations give a frequency factor 𝜎 of 4.5088 (rather than 10.9868 as calculated following the usual scheme [6]) for 68.85% of molecules with “in well” vibrations. As concerns 31.15% of molecules with freely rotating NO₂ groups, the vibrational amplitudes increase enormously for them, to 1.25–1.29 Å for the O10–C1, O9–H5, and O10–H5 distances. In any event, the corresponding parameters should be used as free variables in the refinement of the electron diffraction structure for 31.15% of molecules with free rotations. In all probability, the corresponding coordinates will not give an appreciable contribution to the diffraction picture.

Also note that the amplitudes related to free internal rotations do not change as the vibrational quantum number changes (except small centrifugal effects), only the frequency of rotations increases. For this reason, the results of potential energy surface scan calculations can be used as staring vibrational amplitude values for molecules with free internal rotations.

True, the potential for hindered rotation is strongly anharmonic,𝑉(𝜙)=12^∞_𝑚=1𝑉_𝑚𝑠[1−cos(𝑚𝑠𝜙)],(75) but this is an even function of the torsional angle 𝜙, which cannot contribute to the shrinkage effect. As to the amplitude of this vibration, it can be measured directly as the distance between potential barriers obtained in quantum-mechanical calculations.

11. The Problem of Redundant Coordinates

Attention should be given to calculations with the use of the 𝐵(0)⁻¹ matrix, which is one way or another necessary, in particular, for the determination of the ⟨𝑥⟩ vector. Calculations of vibrational spectra are usually performed with redundant coordinates. Of course, we can always construct a Moore-Penrose pseudoinverse of the 𝐵 matrix after introducing the Eckart coordinates into it. No problems then arise with open systems containing dependences of the type of angles at nodal atoms.

The case is, however, somewhat different with cyclic systems containing nonlinear dependences between bond lengths and ring angles. Eigenvector 𝑙_𝛼 components compatible at 𝑥=0 then become incompatible at 𝑥≠0, and the use of the procedure described above gives corrections (though very insignificant) for bond lengths and valence angles even in the kinematic approximation, in which corrections for these parameters should be zero (Section 3). The following simplified approach to the problem seems to be quite reasonable. According to the Gauss least constraint principle, the difference between the actual and free motion (motion without constraints) should be minimum. Constraint (𝑍) is the quadratic form [21]𝑍=12𝑞^𝛼̈−𝑞_𝛽̈𝐴𝑞_𝛼̈−𝑞_𝛽̈||_𝑡0,(76) where 𝐴=𝜕²ℒ/(𝜕𝑞̇𝜕𝑞̇), 𝑞_𝛼 is the free path, and 𝑞_𝛽 is an admissible path. At time 𝑡₀, the states of the system (𝑞,𝑞̇) on both paths are identical. The admissible path becomes real if 𝑍 is minimum.

For our problem, the free path is the path along which corrections to bond lengths, valence angles, and so forth remain zero in the kinematic approximation. The application of (76) then gives𝑍=𝜔^4𝛼Δ𝑙_𝛼𝐺⁻¹Δ𝑙_𝛼=𝜔^2𝛼Δ𝑙_𝛼𝐹Δ𝑙_𝛼,(77) that is, the problem reduces to the minimization of the (Δ𝑙_𝛼)𝐹(Δ𝑙_𝛼) quadratic form.

In the Shrink09 program, the real path is sought by the introduction of weight factors, the role of which is played by the diagonal matrix 𝐹^diag with the |𝑓_𝑖𝑖|^1/2 elements, where 𝑓_𝑖𝑖 is the 𝑖th diagonal element of the matrix of force constants in internal coordinates,Δ𝑥_𝛼=𝐹^diag𝐹^diag𝐵(0)⁻¹Δ𝑞_𝛼.(78) It can be assumed that the Δ𝑙_𝛼=𝐵(0)Δ𝑥 matrix should minimize constraint. In any event, the difference between the potential energy along the 𝑙_𝛼+Δ𝑙_𝛼 and free paths should be minimum. The use of this procedure considerably decreases corrections to bond lengths obtained in the kinematic approximation for cyclic structures and slightly increases corrections for valence angles, which seems reasonable.

Note that, for open systems, such scaling does not introduce any changes into calculation results.

12. Conclusions

The calculations described above solve the vibrational problem and the problem of the search for parameters necessary for the interpretation of diffraction data at the level of first-order perturbation theory. The solution was obtained using the classical formalism, the harmonic approximation serving as a zero order step. The paper generalizes the results obtained in the preceding works [19, 22–26]. The derivations given in those works, however, contained several inaccuracies (my fault). They had therefore to be refined and repeated here. The initial calculation scheme suggested, however, remained unchanged on the whole. The same scheme can be used in higher perturbation theory orders, but the corresponding quantum-mechanical data are unavailable. All that concerns vibrational amplitudes is entirely new.

In conclusion and by way of illustration, let us consider a fragment of the results obtained for quite a trivial molecule, nitroethane (see Figure 6, all the values below except skewness are in angstrom units; skewness is, of course, dimensionless):Distanceℓ₀ℓ₁Skewness𝑟_𝑒−𝑟_𝑎CCCC1111⋯⋯−−HCNO24390000....00015760069168470000....0100725777271446−−0000....0003003562513542−−−0000....1000601104945165,(79) and so forth (here, 𝑟_𝑒 and 𝑟_𝑎 are the equilibrium and experimentally observed internuclear distances; column ℓ₀ contains amplitudes calculated in the approximation of infinitesimal amplitudes, in which skewness and 𝑟_𝑒−𝑟_𝑎 are, naturally, zero; column ℓ₁, amplitudes calculated as described above).

If a different model of molecular motions is used (calculations are performed in Cartesian coordinates and the quantum-mechanical three-dimensional matrix of cubic constants is used as is), we, for instance, obtainDistanceℓ₀ℓ₁Skewness𝑟_𝑒−𝑟_𝑎CCCC1111⋯⋯−−HCNO24390000....00108605991084760000....0010871515129477−−−0200.95..14.7500283804156−−−−0000...02.001652903467227,(80) and so forth. The exaggerated 𝑟_𝑒−𝑟_𝑎 (and, therefore, skewness) values can hardly be considered realistic.

Attempts at the introduction of corrections for vibrational motions of atoms into X-ray data on molecular crystals were made long ago [27], but were abandoned since for reasons unknown.