Diffraction Measurements and Equilibrium Parameters

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).


Introduction
Diffraction measurements yield some time-and ensembleaveraged 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 2 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 2 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 R(t) can in the first approximation be represented as where R e is the vector of equilibrium geometric parameters, ω i and θ i are the frequencies and phases of vibrational motions, X i is the vector of geometric parameter changes at ω i t + θ i = 2kπ and ω j / = i t + θ j = (2k + 1)π/2 (k = 0, . . .), and n 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 ω i t + θ i values simultaneously satisfy the condition ω i t + θ i = (2k + 1)π/2 (k = 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-toone correspondence between the configuration determined by (1) and parameters measured experimentally, and the O· · · O distance in CO 2 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 R 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 3 O 2 , O=C=C=C=O), the observed 2 Advances in Physical Chemistry 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 3 O 2 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.
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.

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 (L is the Lagrange function) written for the nuclear subsystem takes the form (U is assumed to be independent of the velocities of the particles constituting the system). Here, T(q,q) is the additive part of the Lagrange function (the kinetic energy), U(q) is the potential energy (the nonadditive part), and q stands for some "generalized" coordinates describing the configuration of the system. The kinetic and potential energies are written as where N is the number of vibrational system particles and 3N − 6(5) is the number of vibrational degrees of freedom (all the q and x coordinates are taken to be zero in the equilibrium configuration). In the matrix form, T = (1/2) Mẋ iẋ j and U = (1/2) Fq i q j . Here, M is the diagonal matrix of atomic masses containing three masses of each atom corresponding to motions along the x, y, and z axes and F is the potential energy operator written in generalized coordinates q. The F operator is symmetrical by construction. In the q = Bx generalized coordinates (B is the transformation matrix between the generalized and Cartesian coordinates), the kinetic energy takes the form ( † is the symbol of transposition). The G −1 matrix is symmetrical by construction. This matrix is positive definite, because any motion with a nonzero velocity increases the kinetic energy. The trivial equalities that is, B = B(0) + (1/2)(∂B/∂x), and that is, F = F(0) + (1/3)(∂F/∂q), allow (3) to be rewritten as where B(0), G −1 (0), and F(0) are the matrices corresponding to the approximation of infinitesimal vibrational amplitudes. The derivative of an m × n matrix A with respect to an ldimensional vector b is understood to be a set of m × n × l(∂a i j /∂b k ) values. If the Cartesian system of coordinates is used, In the approximation of infinitesimal amplitudes ((∂G −1 /∂q)q = 0 and (∂F/∂q)q = 0), (8) is rewritten as Advances in Physical Chemistry 3 Solutions to this system of equations are sought in the form q = l cos(ωt + θ) (q = −lω 2 cos(ωt + θ)), that is, where l α 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 G −1 and F matrices are symmetrical and the G −1 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 G −1 matrix (and, simultaneously, to the F matrix: where D 1 is some diagonal matrix and F 1 remains symmetrical (Hermitian) by construction. Since G −1 is positive definite, all the elements of the D 1 matrix are positive numbers. We can therefore construct the D −1/2 . As a result, the kinetic energy matrix transforms into the identity matrix whereas the potential energy matrix As to the identity matrix, its unitary transformation leaves it unchanged, V † 2 IV 2 = I. Let the L matrix be constructed as described above (L = where Ω is the diagonal matrix of the squares of vibrational frequencies and L is the matrix whose columns are the eigenvectors of the vibrational problem, l α . 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 L matrix is independent of the instantaneous configuration of the system. The columns of this matrix are mere linear combinations of generalized coordinates q selected to describe the geometry of the molecule. In what follows, it will be more convenient to use the matrices In the corresponding system of coordinates, G(0) transforms into the identity matrix (I); F(0), into the diagonal matrix of eigenvalues Ω; eigenvectors, into unit vectors l α → e α (e α 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 [(h/4π 2 cν α ) coth(hcν α /2kT)] 1/2 rather than one [6]; this is implied in what follows. Here, h is the Planck constant, c is the velocity of light, k is the Boltzmann constant, and ν α is the αth frequency).
Solutions to (8) will be sought in the form s = s 0 + s 1 (s 0,α = e α cos(ω α t + θ α )) and ω = ω 0 + ω 1 . Let us rewrite (9) and (8) using the notation introduced above and subtract the former from the latter, (no corrections should be introduced into s i s j -type products, because this would cause the appearance of terms smaller in magnitude than s 1 ). Here, s 0 = 3N−6(5) α=1 e α cos(ω α t + θ α ) is a linear combination of "unit" vectors e α (the superposition principle), which are also basis vectors for the construction of s 1 .
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 G matrix. If the model of molecular motions with the F 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 R α structure [6] calculated in a "one-anda-half " approximation (taking into account terms secondorder 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: Is j 1,c + Ωs 4

Kinematic Problem: Equation
There should be no resonance terms containing cos(ω α t +θ α ) in a conservative system. For this reason, ω 1 = 0 for all α (the known result [7]). In (17), L, that is, the increment of the G −1 matrix caused by vibration at an ω β frequency. The e α vector "cuts" the αth column, g (α) β , from the ΔG −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 The solution is sought in the form of a linear combination of the periodic functions present on the right-hand side of (18), where a and b are the sought vector coefficients. Substituting (19) into (18) and equating the corresponding harmonic terms, we obtain where P (β) ± is the diagonal matrix with the elements p Of course, there are no problems with P (β) matrices. As to ΔG −1 β matrices, they are easy to obtain by numerical differentiation.
To summarize, in the "kinematic" approximation, we have the following: 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 2 α and sin 2 α are equal in magnitude and opposite in sign. The s 1,k 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 q 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. (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.

"Centrifugal" Problem: Equation
Let us consider the triatomic fragment shown in Figure 2. Let the AB bond rotate with respect to the A point at an angular velocityq ϕ . The linear velocity of the B point is then v B = rq ϕ , where r is the distance between points A and B, and the kinetic energy of the material point with mass m B is ϕ . Derivative ∂T B /∂r calculations yield m B rq 2 ϕ , which exactly equals the f c = m B rq 2 ϕ centrifugal force that acts on the material point with mass m B which moves at an angular velocityq ϕ with respect to point A.
Similar considerations apply to wagging and torsional vibrations, when the distance from the axis of rotation changes ( Figure 3). We then have f c = m B (r sin(∠CAB)q 2 ϕ ) and , has the meaning of the projection of the centrifugal force onto the AB bond, which rotates about the CA axis at an angular velocityq ϕ (the sine of the ∠CAB angle equals the cosine of the angle between the AB bond and the normal to the axis of rotation). Clearly, an increase in the AB distance always increases the kinetic energy, and the corresponding derivatives are always positive.
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 ∠CAB angle decreases the kinetic energy of torsional motion if this angle is larger than π/2 (the B point then approaches the axis of rotation) and increases it if ∠CAB < π/2 (the B point moves from the axis). The derivative of the kinetic energy of torsional motion with respect to the angular coordinate, ∂T B /∂∠CAB = f c r cos(∠CAB) (Figure 3, f angle c ), can be treated as the centrifugal force acting on the angle. This force equals the product of the length r of the "lever" by the projection of the f c force onto the direction normal to the AB 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 (the resonance term is excluded). Transforming the products of sine functions as in (18), we obtain Substituting (19) into the left-hand side of (23) yields (the notation is the same as before). It follows that Averaging over time and ensemble gives the contribution of the "centrifugal" term to the shrinkage effect, . This, for instance, corresponds to elongation of bonds caused by bending vibrations (bond rotations about a center, see Figure 2). (16) This problem is quite similar to the preceding one. Let us introduce the denotation for the table of "cubic force constants" ∂F/∂s = H. H is a table (three-dimensional matrix) of the third derivatives of the potential energy with 6 Advances in Physical Chemistry 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):

Anharmonic Problem: Equation
(the resonance term is excluded). Substituting (19) into the left-hand side of this equation yields This gives 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 The e α e β vectors "cut" a column with the h αβ... elements from the H 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 g and h vectors, which are, in turn, determined by the lengths of "unit" vectors e (e α = [(h/4π 2 cν α ) coth(hcν α /2kT)] 1/2 ). These vectors become longer than one starting from ω ≈ 84.5 cm −1 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: for the vector of displacements with respect to the equilibrium configuration and for the contribution of the vibration with the ω α frequency to the shrinkage effect. Note that g (α) α (these vectors only differ in the order of differentiation), and, naturally, h (α) α . For this reason, the summation in β can be performed from β = 1 to β = α.
The transition to vectors in internal coordinates q is trivial, where s 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 Here, L is the matrix defined by (11). The solution becomes meaningless if (ω α ±ω β ) 2 −ω 2 γ = 0, for instance, if ω β = 2ω α (the matrix with the [(ω α ± ω β ) 2 − ω 2 γ ] elements then has no inverse; it follows, that the P ± 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.

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 q. For the transition to Advances in Physical Chemistry 7 internuclear distances and the spatial configuration of the system, the results should be transformed into Cartesian coordinates. This transformation is performed as or, according to (33), Because of the smallness of the second term, we can assume B −1 = B(0) −1 in it and exclude it from consideration in this section. Since where x tr is the displacement that appears because of nonlinearity of the transition between internal and Cartesian coordinates and ΔB −1 β is the increment of the B −1 matrix when the configuration of the system changes by the x β vector. In this sum, only the term with cos 2 (ω α t + θ α ) does not vanish in averaging. For this reason, where σ α = (h/4π 2 cν α ) coth(hcν α /2kT) [6]. The x tr α value can be obtained without comparatively laborious calculations of the ΔB Here, x 0,α is the eigenvector in Cartesian coordinates obtained by solving the vibrational problem in the first approximation. Let us calculate the l + 0,α and l − 0,α vectors for the R 0 + x 0,α and R 0 −x 0,α configurations, respectively (R 0 is the vector of the Cartesian coordinates of atoms in the equilibrium configuration). Clearly, Summing these equations yields Substituting the result into (39) yields On the other hand, the x + α and x − α vectors satisfying the conditions l + α = Bx + α and l − α = Bx − α can be calculated by the method of successive approximations by fitting the B(x) matrix and x vector components to obtain q α = B(x)x α . Clearly, x α = (x + α + x − α )/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 l α eigenvector components in generalized coordinates and, therefore, mean bond lengths, valence angles, and so forth. The refinement of the x α Cartesian displacements of atoms and their mean values x , however, leads to substantial shrinkage effect corrections for distances between nonbonded atoms measured by diffraction methods.
Combining (32) and (38), we obtain (frequency factors σ are already contained in the g and h vectors by virtue of the definition of eigenvectors in Section 2) and (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 ith atom by X i , Y i , and Z i , and the corresponding components of the x and x α vectors, by ξ i , υ i , ζ i and ξ α,i , υ α,i , ζ α,i . The R av vector with the components then determines the mean configuration of the molecule and, therefore, the mean values of its rotational constants.

Advances in Physical Chemistry
On the other hand, for the R i j distance between atoms i and j, we have where R 0,i j is the distance between atoms i and j in the equilibrium configuration. The ΔR i j values determine the shrinkage effect for internuclear distances; in diffraction experiments, we measure the R 0,i j + ΔR i j 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 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 B(0) −1 L matrix (the s 0,α (t) value then transforms into x 0,α cos(ω α t + θ α )), (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 2 α-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 x 2 α vector, (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 a + β and a − β vectors are The mean-square amplitudes of changes in internuclear distances x 2 i j are given by where x α,i,x , x α,i,y , and x α,i,z are the components of the x 0,α + 3N−6(5) β=1 (a + β + a − β ) vector corresponding to the displacements of the ith atom along the x, y, and z axes; again, cos 2 (at + b) gives the coefficient 1/2. To pass to the central moment from the moment about zero, we must subtract R α,i j 2 from x 2 α,i j . This is equivalent to excluding terms with cos(α − α) from amplitude calculations.
If the ΔR α,i j and x 2 α,i j values are known, asymmetry parameter x 3 i j / x 2 i j 3/2 calculations present no difficulties. Indeed,

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 quantummechanical programs, first, by a numerical method (the method of finite differences) and, secondly, in Cartesian coordinates. The first circumstance requires estimating possible calculation errors.  (Figure 4). This table can be denoted by the H C (Cartesian) symbol. The cut with the elements h C k,i j (the kth cut) is formed as follows. First, the kth Cartesian coordinate of the molecule x k changes by a small value δ, and the Hessian H of the molecule is calculated in this intentionally distorted nonequilibrium configuration, H(x 1 , . . . , x k + δ, . . . , x 3N ).

Errors in Cubic Force
The same is made after the subtraction of δ from x k , H(x 1 , . . . , x k − δ, . . . , x 3N ). Next, we subtract the second result from the first one and divide the difference by 2δ. As a result, we obtain the kth cut of the table of cubic constants, It follows that the h C k,i j element of the table of cubic force constants is calculated as The equation for calculating the error in the h C k,i j element obtained this way is well known [11], where μ ∈ [x k −δ, x k +δ]. 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 h i j values themselves.
It follows from this equation that, first, the error in h C k,i j is not a monotonic function of δ: the smaller δ, the smaller the truncation error, but the larger the contribution of the error in h i j . The δ values can, of course, be optimized, and the corresponding algorithms are well known [10], but the use of δ of its own for each h C k,i j 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 h C
k,i j elements that differ only by the order of indices are different. Indeed, (57) It follows that we cannot expect to obtain the theoretical equality The following procedure for estimating calculation errors can be suggested. Let index α run over all the permutations of the i, j, k indices. The mean h C α value is then Let us put In other words, Δ is the root-mean-square deviation for the given α (for the given set of i, j, and k indices), and Δ is the half-spread of the h C α values. We put where maxima are sought over all α (over all different sets of i, j, and k indices). Δ can be suggested as an estimate of errors in off-diagonal H C matrix elements, and Δ , as a (conventional) estimate of errors in h C iii . 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 −4 au and Δ = 0.984625 × 10 −4 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.

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 where (∂F/∂q) i jk is the three-dimensional matrix of the third derivatives of potential energy with respect to internal coordinates, H i jk is, as before, the three-dimensional matrix of the third derivatives of potential energy with respect to Cartesian coordinates, and B j i is the transformation matrix between internal and Cartesian coordinates, q j = B j i x i (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 {x 1 , . . . , x k + δ, . . . , x 3N } configuration (see the preceding subsection) will be labeled by the superscript "+," and the values related to the x k − δ configuration, by the superscript "−." Let us introduce the denotations where B 0 is the B matrix calculated for the equilibrium configuration and Δ k = (∂B/∂x k )δ.
Calculations of the kth cut of the H C matrix begin with calculations of the gradient vector for the given configuration and then its derivatives with respect to the atomic coordinates, (because ∂/∂x = (∂/∂q)(∂q/∂x)). Substituting this result and a similar equation for the configuration with x k − δ into (54) yields where F is the matrix of second-order force constants in internal coordinates. It follows that the numerator in (54) takes the form (indeed, F + + F − = 2F and (∂Δ k /∂x) − = −(∂Δ k /∂x) + ). We therefore have Here, ∂ 2 U/∂q∂x k = B(0) −1 (∂ 2 U/∂x∂x k ); that is, this vector is obtained from the kth column of the H matrix, and the other matrices are calculated by the finite-difference method. The right-hand side of this equation gives the kth cut of some H C 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 H I (internal) symbol, Above, the algorithm of calculations is described. It corresponds to the following analytic equations: (the last term differs from the preceding ones in the order of indices),

Scaling of Cubic Force Constants
A popular trend of recent years is to scale quantummechanical 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 q 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 l α 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 where δ is the shift along the kth internal coordinate, kth coordinate according to Pulay, or kth eigenvector. If the transformation (D is a diagonal matrix) scales the force field in the equilibrium configuration, it is easy to accept that the same D matrix should scale the force fields of configurations slightly changed with respect to equilibrium, that is, 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 suggested in [15] without any justification is, of course, absolutely incorrect. Conversely, in [16][17][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]).

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.

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 (C s symmetry) was considered. Torsional vibration of the NO 2 group in this molecule transforms under the A representation of the C s group. However, the minimum energy path calculated by the authors included changes in the C-C and C-N distances and other coordinates of A 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 AB 2 triatomic molecule. Clearly, the ∠BAB angle should change along the corresponding minimum energy path. But the angular coordinate transforms under the A 1 representation, whereas the antisymmetric stretching coordinate, under the B 1 (or B 2 ) representation. In addition, it is impossible to make the angle vibrate at the antisymmetric stretching vibration frequency.
One more example is a linear CO 2 -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.
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 2 group rotations about the N-C bond, 26.4 cm −1 (the B3LYP/aug-cc-pVTZ data).
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 2 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, 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.

The Problem of Redundant Coordinates
Attention should be given to calculations with the use of the B(0) −1 matrix, which is one way or another necessary, in particular, for the determination of the x vector. Calculations of vibrational spectra are usually performed with redundant coordinates. Of course, we can always construct a Moore-Penrose pseudoinverse of the B 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 l α components compatible at x = 0 then become incompatible at x / = 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 (Z) is the quadratic form [21] where A = ∂ 2 L/(∂q∂q), q α is the free path, and q β is an admissible path. At time t 0 , the states of the system (q,q) on both paths are identical. The admissible path becomes real if Z is minimum. For our problem, the free path is the path along which corrections to bond lengths, valence angles, and so forth Advances in Physical Chemistry 13 remain zero in the kinematic approximation. The application of (76) then gives that is, the problem reduces to the minimization of the (Δl α )F(Δl α ) 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 F diag with the | f ii | 1/2 elements, where f ii is the ith diagonal element of the matrix of force constants in internal coordinates, It can be assumed that the Δl α = B(0)Δx matrix should minimize constraint. In any event, the difference between the potential energy along the l α + Δl α 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.

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][23][24][25][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): and so forth (here, r e and r a are the equilibrium and experimentally observed internuclear distances; column 0 contains amplitudes calculated in the approximation of infinitesimal amplitudes, in which skewness and r e − r a are, naturally, zero; column 1 , 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, obtain and so forth. The exaggerated r e − r a (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.