# Diffraction Measurements and Equilibrium Parameters

Authors: Victor A. Sipachev
Journal: Advances in Physical Chemistry
 (2011)

---

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


---

## Body

## 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 CO2 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 CO2 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 systemR(t) can in the first approximation be represented as(1)R(t)=Re+∑i=1nXicos⁡(ωit+θi),
where Re is the vector of equilibrium geometric parameters, ωi and θi are the frequencies and phases of vibrational motions, Xi is the vector of geometric parameter changes at ωit+θi=2kπ and ωj≠it+θ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 ωit+θi values simultaneously satisfy the condition ωit+θ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-to-one correspondence between the configuration determined by (1) and parameters measured experimentally, and the O⋯O distance in CO2 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 (C3O2, 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 C3O2 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(2)ddt(∂L∂q̇)-(∂L∂q)=0
(ℒ is the Lagrange function) written for the nuclear subsystem takes the form(3)ddt(∂T(q,q̇)∂q̇)-∂∂q(T(q,q̇)-U(q))=0
(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(4)T=(12)∑i=1Nmi(ẋi2+ẏi2+żi2),U=(12)∑i=13N-6(5)(∂2U∂qi∂qj)qiqj+⋯,
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)Mẋiẋj and U=(1/2)Fqiqj. 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(5)T=(12)(B-1q̇j)†MB-1q̇i=G-1q̇iq̇j
(† 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(6)qi=∑j3N(∂qi∂xj)xj+(12)∑j,k3N(∂2qi∂xj∂xk)xjxk+⋯,
that is, B=B(0)+(1/2)(∂B/∂x), and(7)U=(12)∑i,j3N-6(5)(∂2U∂qi∂qj)qiqj+(16)∑i,j,k3N-6(5)(∂3U∂qi∂qj∂qk)qiqjqk+⋯,
that is, F=F(0)+(1/3)(∂F/∂q), allow (3) to be rewritten as(8)G-1(0)q̈j+(12)(∂G-1∂q)(qiq̈j+q̇iq̇j)-(12)(∂G-1∂q)q̇iq̇j+F(0)qi+(12)(∂F∂q)qiqj=0,
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 l-dimensional vector b is understood to be a set of m×n×l(∂aij/∂bk) values. If the Cartesian system of coordinates is used, G-1(q)=G-1(0)=M.In the approximation of infinitesimal amplitudes ((∂G-1/∂q)q=0 and (∂F/∂q)q=0), (8) is rewritten as(9)G-1q̈j+Fqj=0.
Solutions to this system of equations are sought in the form q=lcos(ωt+θ) (q̈=-lω2cos(ωt+θ)), that is,(10)G-1lαωα2-Flα=0,GFlα=lαωα2,
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 theG-1 matrix (and, simultaneously, to the F matrix: V1†G-1V1=D1 and V1†FV1=F1, where D1 is some diagonal matrix and F1 remains symmetrical (Hermitian) by construction. Since G-1 is positive definite, all the elements of the D1 matrix are positive numbers. We can therefore construct the D1-1/2 diagonal matrix such that D1-1/2D1D1-1/2=I (I is the identity matrix). As a result, the kinetic energy matrix transforms into the identity matrix whereas the potential energy matrix D1-1/2F1D1-1/2=F2 remains symmetrical. It can therefore be diagonalized with the use of a unitary (orthogonal) transformation, V2†D1-1/2V1†FV1D1-1/2V2=D2. As to the identity matrix, its unitary transformation leaves it unchanged, V2†IV2=I.Let theL matrix be constructed as described above (L=V1D1-1/2V2). We have(11)L-1G(L-1)†L†FL=L-1GFL=Ω,GFL=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, theL 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 (12)B=L-1B,G=BM-1B†,F=L†FL.
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π2cνα)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=s0+s1 (s0,α=eαcos(ωαt+θα)) and ω=ω0+ω1. Let us rewrite (9) and (8) using the notation introduced above and subtract the former from the latter,
(13)Is̈1j+Ωs1j=-(12)(∂G-1∂s)(s0is̈0j+ṡ0iṡ0j)+(12)(∂G-1∂s)ṡ0iṡ0j-(12)(∂F∂s)s0is0j
(no corrections should be introduced into sisj-type products, because this would cause the appearance of terms smaller in magnitude than s1). Here, s0=∑α=13N-6(5)eαcos(ωαt+θα) is a linear combination of “unit’’ vectors eα (the superposition principle), which are also basis vectors for the construction of s1.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-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:(14)Is̈1,kj+Ωs1,kj=-(12)(∂G-1∂s)(s0is̈0j+ṡ0iṡ0j),(15)Is̈1,cj+Ωs1,cj=(12)(∂G-1∂s)ṡ0iṡ0j,(16)Is̈1,aj+Ωs1,aj=-(12)(∂F∂s)s0is0j.
Clearly, s1=s1,k+s1,c+s1,a 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 G elements) on the instantaneous configuration of the system. We have(17)Is̈1,k,α+Ωs1,k,α=12ωα2cos⁡(ωαt+θα)∑β=13N-6(5)cos⁡(ωβt+θβ)ΔGβ-1eα-12ωαsin⁡(ωαt+θα)∑β=13N-6(5)ωβsin⁡(ωβ+θβ)ΔGβ-1eα+2Iω0,αω1,αcos⁡(ωαt+θα)eα.
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), Δ𝔊β-1 is L†[(∂G-1/∂q)qβ]L, that is, the increment of the 𝔊-1 matrix caused by vibration at an ωβ frequency. The eα 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(18)Is̈1,k,α+Ωs1,k,α=14ωα∑β=13N-6(5)(ωα+ωβ)cos⁡[(ωα+ωβ)t+(θα+θβ)]ĝβ(α)+14ωα∑β=13N-6(5)(ωα-ωβ)cos⁡[(ωα-ωβ)t+(θα-θβ)]ĝβ(α).
The solution is sought in the form of a linear combination of the periodic functions present on the right-hand side of (18),(19)s1,k,α=∑β=13N-6(5)aβcos⁡[(ωα+ωβ)t+(θα+θβ)]+bβcos⁡[(ωα-ωβ)t+(θα-θβ)],s̈1,k,α=-∑β=13N-6(5)aβ(ωα+ωβ)2cos⁡[(ωα+ωβ)t+(θα+θβ)]+bβ(ωα-ωβ)2cos⁡[(ωα-ωβ)t+(θα-θβ)],
where a and b are the sought vector coefficients. Substituting (19) into (18) and equating the corresponding harmonic terms, we obtain(20)aβ=14ωα(ωα+ωβ)P+(β)ĝβ(α),bβ=14ωα(ωα-ωβ)P-(β)ĝβ(α),
where 𝔓±(β) is the diagonal matrix with the elements 𝔭γγ(β)=1/[ωγ2-(ωα±ωβ)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:(21)s1,k,α=14∑β=13N-6(5)ωα(ωα+ωβ)cos⁡[(ωα+ωβ)t+(θα+θβ)]P+(β)ĝβ(α)+ωα(ωα-ωβ)cos⁡[(ωα-ωβ)t+(θα-θβ)]P-(β)ĝβ(α).
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 cos2α and sin2α are equal in magnitude and opposite in sign. The s1,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.


## 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 Figure2. Let the AB bond rotate with respect to the A point at an angular velocity q̇φ. The linear velocity of the B point is then vB=rq̇φ, where r is the distance between points A and B, and the kinetic energy of the material point with mass mB is TB=(1/2)mBr2q̇φ2. Derivative ∂TB/∂r calculations yield mBrq̇φ2, which exactly equals the fc=mBrq̇φ2 centrifugal force that acts on the material point with mass mB which moves at an angular velocity q̇φ with respect to point A.Figure 2
Bond rotation about a center.Similar considerations apply to wagging and torsional vibrations, when the distance from the axis of rotation changes (Figure3). We then have fc=mB(rsin(∠CAB)q̇φ2) and TB=(1/2)mB(rsin(∠CAB)2q̇φ2). The derivative of TB with respect to r, ∂TB/∂r=fcsin(∠CAB) (Figure 3, fcbond), has the meaning of the projection of the centrifugal force onto the AB bond, which rotates about the CA axis at an angular velocity q̇φ (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.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 Figure3 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, ∂TB/∂∠CAB=fcrcos(∠CAB) (Figure 3, fcangle), 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 fc 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(22)Is̈1,c,α+Ωs1,c,α=12ωαsin⁡(ωαt+θα)×∑β=13N-6(5)ωβsin⁡(ωβ+θβ)ΔGβ-1eα
(the resonance term is excluded). Transforming the products of sine functions as in (18), we obtain(23)Is̈1,c,α+Ωs1,c,α=-14ωα∑β=13N-6(5)ωβcos⁡[(ωα+ωβ)t+(θα+θβ)]ĝβ(α)+14ωα∑β=13N-6(5)ωβcos⁡[(ωα-ωβ)t+(θα-θβ)]ĝβ(α).
Substituting (19) into the left-hand side of (23) yields(24)aβ=-(14)ωαωβP+(β)ĝβ(α),bβ=-(14)ωαωβP-(β)ĝβ(α)
(the notation is the same as before). It follows that(25)s1,c,α=14∑β=13N-6(5)-ωαωβcos⁡[(ωα+ωβ)t+(θα+θβ)]P+(β)ĝβ(α)+ωαωβcos⁡[(ωα-ωβ)t+(θα-θβ)]P-(β)ĝβ(α).
Averaging over time and ensemble gives the contribution of the “centrifugal’’ term to the shrinkage effect,(26)〈s1,c〉=14∑α=13N-6(5)ωα2P-(α)ĝα(α)=14∑α=13N-6(5)ωα2Ω-1ĝα(α)
(clearly, 𝔓-(α)=Ω-1). 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’’∂𝔉/∂s=ℌ. ℌ 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):(27)Is̈1,a,α+Ωs1,a,α=-12ωαcos⁡(ωαt+θα)∑β=13N-6(5)ωβcos⁡(ωbt+θβ)Heαeβ=-14∑β=13N-6(5){cos⁡[(ωα+ωβ)+(θα+θβ)]Heαeβ+cos⁡[(ωα-ωβ)+(θα-θβ)]×Heαeβ}
(the resonance term is excluded). Substituting (19) into the left-hand side of this equation yields(28)aβ=-(14)P+(β)Heαeβ,bβ=-(14)P-(β)Heαeβ.
This gives(29)s1,a,α=-14∑β=13N-6(5){cos⁡[(ωα+ωβ)+(θα+θβ)]P+(β)+cos⁡[(ωα-ωβ)+(θα-θβ)]×P-(β)}Heαeβ.
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(30)〈s1,a,α〉=-14Ω-1Heαeβ.
The eαeβ 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 e (eα=[(h/4π2cνα)coth(hcνα/2kT)]1/2). These vectors become longer than one starting from ω≈84.5cm-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:(31)sα=s0,α(t)+14∑β=13N-6(5){cos⁡[(ωα-ωβ)+(θα-θβ)]P-(β)ĥβ(α)[ωα(ωα+ωβ)-ωαωβ]×cos⁡[(ωα+ωβ)t+(θα+θβ)]×P+(β)ĝβ(α)-cos⁡[(ωα+ωβ)+(θα+θβ)]×P+(β)ĥβ(α)+[ωα(ωα-ωβ)+ωαωβ]×cos⁡[(ωα-ωβ)t+(θα-θβ)]P-(β)ĝβ(α)-cos⁡[(ωα-ωβ)+(θα-θβ)]P-(β)ĥβ(α)}
for the vector of displacements with respect to the equilibrium configuration and(32)〈sα〉=14Ω-1{ωα2ĝα(α)-ĥα(α)}
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 coordinatesq is trivial,(33)lα=l0,α+Ls1,α,
where s1,α is the sum of small terms in (31) (the sum of corrections to eigenvectors obtained by solving the kinematic, centrifugal, and anharmonic problems), and(34)〈lα〉=L〈s1,α〉.
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 𝔓± 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) coordinatesq. 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(35)xα=B-1lαcos⁡(ωαt+θα),
or, according to (33),(36)xα=B-1l0,αcos⁡(ωαt+θα)+B-1Ls1,α.
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 B-1=B(0)-1+(1/2)∑β=13N-6(5)(∂B-1/∂x)xβcos(ωβt+θβ), we have(37)xαtr⁡=12∑β=13N-6(5)ΔBβ-1l0,αcos⁡(ωαt+θα)cos⁡(ωβt+θβ),
where xtr 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 cos2(ωαt+θα) does not vanish in averaging. For this reason,(38)〈xαtr⁡〉=〈12∂B-1∂xxαl0,αcos⁡2(ωαt+θα)〉=14σαΔBα-1l0,α,
where σα=(h/4π2cνα)coth(hcνα/2kT) [6]. The 〈xαtr〉 value can be obtained without comparatively laborious calculations of the ΔBα-1 matrices. For this purpose, let us rewrite (38) as follows. Since B-1=(B(0)+(1/2)ΔB)-1=(I+(1/2)B(0)-1ΔB)-1B(0)-1≈(I-(1/2)B(0)-1ΔB)B(0)-1, we have ΔB-1≈B(0)-1ΔBB(0)-1. On the other hand, B(0)-1l0,α=x0,α. Therefore,(39)〈xαtr⁡〉≈(14)σαB(0)-1ΔBαB(0)-1l0,α=(14)σαB(0)-1ΔBαx0,α.
Here, x0,α is the eigenvector in Cartesian coordinates obtained by solving the vibrational problem in the first approximation. Let us calculate the l0,α+ and l0,α- vectors for the R0+x0,α and R0-x0,α configurations, respectively (R0 is the vector of the Cartesian coordinates of atoms in the equilibrium configuration). Clearly,(40)l0,α+=(B(0)+(12)ΔBα)x0,α,l0,α-=(B(0)-(12)ΔBα)(-x0,α).
Summing these equations yields(41)l0,α++l0,α-=ΔBαx0,α.
Substituting the result into (39) yields(42)〈xαtr⁡〉=(14)σαB(0)-1(lα++lα-).
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 changelα 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(43)〈xα〉=14B(0)-1(σαΔBα-1l0,α+L[Ωωα2ĝα(α)-ĥα(α)])
(frequency factors σ are already contained in the 𝔤̂ and 𝔥̂ vectors by virtue of the definition of eigenvectors in Section 2) and(44)〈x〉=∑α=13N-6(5)〈xα〉
(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 Xi, Yi, and Zi, and the corresponding components of the 〈x〉 and 〈xα〉 vectors, by ξi, υi, ζi and ξα,i, υα,i, ζα,i. The Rav vector with the components(45)Xi+∑α=13N-6(5)ξα,i,Yi+∑α=13N-6(5)υα,i,Zi+∑α=13N-6(5)ζα,i
then determines the mean configuration of the molecule and, therefore, the mean values of its rotational constants.On the other hand, for theRij distance between atoms i and j, we have(46)ΔRij=∑α=13N-6(5){([(Xi-Xj)+(ξα,i-ξα,j)]2+[(Yi-Yj)+(υα,i-υα,j)]2+[(Zi-Zj)+(ζα,i-ζα,j)]2)1/2-R0,ij[(Xi-Xj)+(ξα,i-ξα,j)]2},
where R0,ij is the distance between atoms i and j in the equilibrium configuration. The ΔRij values determine the shrinkage effect for internuclear distances; in diffraction experiments, we measure the R0,ij+ΔRij 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(47)xαtr⁡=14∑β=13N-6(5)ΔBβ-1l0,α(cos⁡[(ωα+ωβ)t+(θα+θβ)]+cos⁡[(ωα-ωβ)t+(θα-θβ)]).
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)-1L matrix (the s0,α(t) value then transforms into x0,α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 cos2α-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,(48)〈xα2〉=σαSq(x0,α)+12∑β=13N-6(5)(Sq(aβ+)+Sq(aβ-))
(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(49)aβ+=14{(12)[ωα(ωα+ωβ)-ωαωβ]B(0)-1LP+(β)ĝβ(α)+B(0)-1LP+(β)ĥβ(α)+(12)σασβΔBβ-1l0,α},aβ-=14{(12)[ωα(ωα-ωβ)+ωαωβ]B(0)-1LP-(β)ĝβ(α)+B(0)-1LP-(β)ĥβ(α)+(12)σασβΔBβ-1l0,α}.
The mean-square amplitudes of changes in internuclear distances 〈xij2〉 are given by(50)〈xij2〉=12∑α=13N-6(5){(x¯α,i,x-x¯α,j,x)2+(x¯α,i,y-x¯α,j,y)2+(x¯α,i,z-x¯α,j,z)2},
where x¯α,i,x, x¯α,i,y, and x¯α,i,z are the components of the x0,α+∑β=13N-6(5)(aβ++aβ-) vector corresponding to the displacements of the ith atom along the x, y, and z axes; again, 〈cos2(at+b)〉 gives the coefficient 1/2. To pass to the central moment from the moment about zero, we must subtract 〈Rα,ij〉2 from 〈xα,ij2〉. This is equivalent to excluding terms with cos(α-α) from amplitude calculations.If theΔRα,ij and 〈xα,ij2〉 values are known, asymmetry parameter 〈xij3〉/〈xij2〉3/2 calculations present no difficulties. Indeed,(51)〈xij3〉=∑α=13N-6(5)ΔRα,ij3.


## 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 a3N×3N×3N cubic table (N is the number of atoms in the molecule). Let us consider how the kth “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 theHC (Cartesian) symbol. The cut with the elements hk,ijC (the kth cut) is formed as follows. First, the kth Cartesian coordinate of the molecule xk changes by a small value δ, and the Hessian H of the molecule is calculated in this intentionally distorted nonequilibrium configuration,(52)Hk+=H(x1,…,xk+δ,…,x3N).
The same is made after the subtraction of δ from xk,(53)Hk-=H(x1,…,xk-δ,…,x3N).
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,(54)HkC=Hk+-Hk-2δ.
It follows that the hk,ijC element of the table of cubic force constants is calculated as(55)hk,ijC=hij+-hij-2δ.
The equation for calculating the error in the hk,ijC element obtained this way is well known [11], (56)Δhk,ijC=16δ2|∂3hij∂xk3|μ+ɛ(ij)δ,
where μ∈[xk-δ,xk+δ]. Here, the first term is the so-called truncation error. Its origin becomes clear if the difference hij+(xk+δ)-hij-(xk-δ) is expanded in powers of δ. The second term is determined by the error in the hij values themselves.It follows from this equation that, first, the error inhk,ijC is not a monotonic function of δ: the smaller δ, the smaller the truncation error, but the larger the contribution of the error in hij. The δ values can, of course, be optimized, and the corresponding algorithms are well known [10], but the use of δ of its own for each hk,ijC 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 inhk,ijC elements that differ only by the order of indices are different. Indeed,(57)Δhk,ijC=16δ2|∂3hij∂xk3|μ+ɛ(ij)δ,Δhj,kiC=16δ2|∂3hki∂xk3|μ+ɛ(ki)δ.
It follows that we cannot expect to obtain the theoretical equality hijkC=hikjC=hjikC=hjkiC=hkijC=hkjiC 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 i,j,k indices. The mean hαC¯ value is then(58)hαC¯=16∑αhαC.
Let us put(59)Δα=(16∑α(hαC-hαC¯)2)1/2,Δα′=max⁡hαC-min⁡hαC2.
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(60)Δ=max⁡Δα,Δ′=max⁡Δα′,
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 HC matrix elements, and Δ′, as a (conventional) estimate of errors in hiiiC. 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.


### 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)(∂F∂q)ijk=HijkC(BjiBkjBjk)-1,
where (∂F/∂q)ijk is the three-dimensional matrix of the third derivatives of potential energy with respect to internal coordinates, Hijk is, as before, the three-dimensional matrix of the third derivatives of potential energy with respect to Cartesian coordinates, and Bij is the transformation matrix between internal and Cartesian coordinates, qj=Bijxi (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{x1,…,xk+δ,…,x3N} configuration (see the preceding subsection) will be labeled by the superscript “+,’’ and the values related to the xk-δ configuration, by the superscript “−.’’ Let us introduce the denotations(62)B(k)+=B0+Δk,B(k)-=B0-Δk,
where B0 is the B matrix calculated for the equilibrium configuration and Δk=(∂B/∂xk)δ.Calculations of thekth cut of the HC matrix begin with calculations of the gradient vector for the given configuration and then its derivatives with respect to the atomic coordinates,(63)∂U∂x=(B(k)+)†∂U∂q,∂2U∂x2=(B(k)+)†∂2U∂q2B(k)++∂∂x(B(k)+)†∂U∂q
(because ∂/∂x=(∂/∂q)(∂q/∂x)). Substituting this result and a similar equation for the configuration with xk-δ into (54) yields(64)Hk+=B0†F+B0+B0†F+Δk+Δk†F+B0+(∂Δk†∂x)+(∂U∂q)+,Hk-=B0†F-B0-B0†F-Δk-Δk†F+B0-(∂Δk†∂x)-(∂U∂q)-,
where F is the matrix of second-order force constants in internal coordinates. It follows that the numerator in (54) takes the form(65)Hk+-Hk-=B0†(F+-F-)B0+2B0†FΔk+2Δk†FB0+(∂Δk†∂x)(∂U∂q+-∂U∂q-)
(indeed, F++F-=2F and (∂Δk/∂x)-=-(∂Δk/∂x)+). We therefore have(66)B0†(F+-F-)B02δ=Hk+-Hk-2δ-1δ[B0†F(∂B∂xk)+(∂B∂xk)†FB0+(∂Δk∂x)(∂2U∂q∂xk)].
Here, ∂2U/∂q∂xk=B(0)-1(∂2U/∂x∂xk); 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 HI (internal) symbol,(67)HijkI=Hijk′C(BjiBkjBjk)-1.
Above, the algorithm of calculations is described. It corresponds to the following analytic equations:(68)∂U∂xk=(Bkl)(∂U∂ql),∂2U(∂xk∂xl)=[∂2U(∂qk∂ql)]BlkBkl+(∂Bkl∂xl)(∂U∂ql),∂3U(∂xk∂xl∂xm)=[∂3U(∂qk∂ql∂qm)]BklBmkBlm+[∂2U(∂qk∂ql)]×Blk(∂Bkl∂xm)+[∂2U(∂qk∂ql)](∂Blk∂xm)Bkl+[∂2U(∂ql∂qm)]Bmk(∂Bkl∂xl)
(the last term differs from the preceding ones in the order of indices),(69)∂3U(∂qk∂ql∂qm)={[∂3U(∂xk∂xl∂xm)]-[∂2U(∂qk∂ql)]×Blk(∂Bkl∂xm)-[∂2U(∂qk∂ql)](∂Blk∂xm)Bkl-[∂2U(∂ql∂qm)]Bmk(∂Bkl∂xl)}×(BklBmkBlm)-1.


## 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 a3N×3N×3N cubic table (N is the number of atoms in the molecule). Let us consider how the kth “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 theHC (Cartesian) symbol. The cut with the elements hk,ijC (the kth cut) is formed as follows. First, the kth Cartesian coordinate of the molecule xk changes by a small value δ, and the Hessian H of the molecule is calculated in this intentionally distorted nonequilibrium configuration,(52)Hk+=H(x1,…,xk+δ,…,x3N).
The same is made after the subtraction of δ from xk,(53)Hk-=H(x1,…,xk-δ,…,x3N).
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,(54)HkC=Hk+-Hk-2δ.
It follows that the hk,ijC element of the table of cubic force constants is calculated as(55)hk,ijC=hij+-hij-2δ.
The equation for calculating the error in the hk,ijC element obtained this way is well known [11], (56)Δhk,ijC=16δ2|∂3hij∂xk3|μ+ɛ(ij)δ,
where μ∈[xk-δ,xk+δ]. Here, the first term is the so-called truncation error. Its origin becomes clear if the difference hij+(xk+δ)-hij-(xk-δ) is expanded in powers of δ. The second term is determined by the error in the hij values themselves.It follows from this equation that, first, the error inhk,ijC is not a monotonic function of δ: the smaller δ, the smaller the truncation error, but the larger the contribution of the error in hij. The δ values can, of course, be optimized, and the corresponding algorithms are well known [10], but the use of δ of its own for each hk,ijC 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 inhk,ijC elements that differ only by the order of indices are different. Indeed,(57)Δhk,ijC=16δ2|∂3hij∂xk3|μ+ɛ(ij)δ,Δhj,kiC=16δ2|∂3hki∂xk3|μ+ɛ(ki)δ.
It follows that we cannot expect to obtain the theoretical equality hijkC=hikjC=hjikC=hjkiC=hkijC=hkjiC 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 i,j,k indices. The mean hαC¯ value is then(58)hαC¯=16∑αhαC.
Let us put(59)Δα=(16∑α(hαC-hαC¯)2)1/2,Δα′=max⁡hαC-min⁡hαC2.
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(60)Δ=max⁡Δα,Δ′=max⁡Δα′,
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 HC matrix elements, and Δ′, as a (conventional) estimate of errors in hiiiC. 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.


## 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)(∂F∂q)ijk=HijkC(BjiBkjBjk)-1,
where (∂F/∂q)ijk is the three-dimensional matrix of the third derivatives of potential energy with respect to internal coordinates, Hijk is, as before, the three-dimensional matrix of the third derivatives of potential energy with respect to Cartesian coordinates, and Bij is the transformation matrix between internal and Cartesian coordinates, qj=Bijxi (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{x1,…,xk+δ,…,x3N} configuration (see the preceding subsection) will be labeled by the superscript “+,’’ and the values related to the xk-δ configuration, by the superscript “−.’’ Let us introduce the denotations(62)B(k)+=B0+Δk,B(k)-=B0-Δk,
where B0 is the B matrix calculated for the equilibrium configuration and Δk=(∂B/∂xk)δ.Calculations of thekth cut of the HC matrix begin with calculations of the gradient vector for the given configuration and then its derivatives with respect to the atomic coordinates,(63)∂U∂x=(B(k)+)†∂U∂q,∂2U∂x2=(B(k)+)†∂2U∂q2B(k)++∂∂x(B(k)+)†∂U∂q
(because ∂/∂x=(∂/∂q)(∂q/∂x)). Substituting this result and a similar equation for the configuration with xk-δ into (54) yields(64)Hk+=B0†F+B0+B0†F+Δk+Δk†F+B0+(∂Δk†∂x)+(∂U∂q)+,Hk-=B0†F-B0-B0†F-Δk-Δk†F+B0-(∂Δk†∂x)-(∂U∂q)-,
where F is the matrix of second-order force constants in internal coordinates. It follows that the numerator in (54) takes the form(65)Hk+-Hk-=B0†(F+-F-)B0+2B0†FΔk+2Δk†FB0+(∂Δk†∂x)(∂U∂q+-∂U∂q-)
(indeed, F++F-=2F and (∂Δk/∂x)-=-(∂Δk/∂x)+). We therefore have(66)B0†(F+-F-)B02δ=Hk+-Hk-2δ-1δ[B0†F(∂B∂xk)+(∂B∂xk)†FB0+(∂Δk∂x)(∂2U∂q∂xk)].
Here, ∂2U/∂q∂xk=B(0)-1(∂2U/∂x∂xk); 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 HI (internal) symbol,(67)HijkI=Hijk′C(BjiBkjBjk)-1.
Above, the algorithm of calculations is described. It corresponds to the following analytic equations:(68)∂U∂xk=(Bkl)(∂U∂ql),∂2U(∂xk∂xl)=[∂2U(∂qk∂ql)]BlkBkl+(∂Bkl∂xl)(∂U∂ql),∂3U(∂xk∂xl∂xm)=[∂3U(∂qk∂ql∂qm)]BklBmkBlm+[∂2U(∂qk∂ql)]×Blk(∂Bkl∂xm)+[∂2U(∂qk∂ql)](∂Blk∂xm)Bkl+[∂2U(∂ql∂qm)]Bmk(∂Bkl∂xl)
(the last term differs from the preceding ones in the order of indices),(69)∂3U(∂qk∂ql∂qm)={[∂3U(∂xk∂xl∂xm)]-[∂2U(∂qk∂ql)]×Blk(∂Bkl∂xm)-[∂2U(∂qk∂ql)](∂Blk∂xm)Bkl-[∂2U(∂ql∂qm)]Bmk(∂Bkl∂xl)}×(BklBmkBlm)-1.


## 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 ofq 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 Section7, we can write(70)HkI=Fk+-Fk-2δ,hijkI=fij+-fij-2δ,
where δ is the shift along the kth internal coordinate, kth coordinate according to Pulay, or kth eigenvector. If the transformation(71)Fscaled=D1/2FtheorD1/2
(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,(72)(fij+-fij-)scaled=(didj)1/2(fij+-fij-)theor.
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(73)B†{D1/2FtheorD1/2}B={B†D1/2}Ftheor{D1/2B},
that is, qitheor=qscaled/d1/2. This means that the 2δ denominator (shift along the kth coordinate) should be divided by dk1/2,(74)hijkI,scaled=(didj)1/2(fij+-fij-)theor2δ/dk1/2=(didjdk)1/2hijkI,theor.
The equation hijkscaled=(didjdk)1/3hijktheor 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 forexperimentally 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 (Cs symmetry) was considered. Torsional vibration of the NO2 group in this molecule transforms under the A′′ representation of the Cs 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 anAB2 triatomic molecule. Clearly, the ∠BAB angle should change along the corresponding minimum energy path. But the angular coordinate transforms under the A1 representation, whereas the antisymmetric stretching coordinate, under the B1 (or B2) representation. In addition, it is impossible to make the angle vibrate at the antisymmetric stretching vibration frequency.One more example is a linear CO2-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 Figure5.Figure 5
Potential energy curve along a torsional coordinate.The horizontal lines in Figure5 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 NO2 group rotations about the N–C bond, 26.4 cm−1 (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 NO2 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,(75)V(ϕ)=(12)∑m=1∞Vms[1-cos⁡(msϕ)],
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 theB(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. Eigenvectorlα 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](76)Z=12((q̈α-q̈β)A(q̈α-q̈β))|t0,
where A=∂2ℒ/(∂q̇∂q̇), qα is the free path, and qβ is an admissible path. At time t0, 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 remain zero in the kinematic approximation. The application of (76) then gives(77)Z=ωα4(Δlα)G-1(Δlα)=ωα2(Δlα)F(Δlα),
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 matrixFdiag with the |fii|1/2 elements, where fii is the ith diagonal element of the matrix of force constants in internal coordinates,(78)Δxα=Fdiag⁡(Fdiag⁡B(0))-1Δqα.
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.


## 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 Figure6, all the values below except skewness are in angstrom units; skewness is, of course, dimensionless):(79)Distancel0l1Skewnessre-raC1–C20.05060.05240.0063-0.0196C1–H40.07680.07760.0025-0.0145C1⋯N30.06940.0771-0.03540.0041C1⋯O90.10170.1274-0.3512-0.1605,
and so forth (here, re and ra are the equilibrium and experimentally observed internuclear distances; column ℓ0 contains amplitudes calculated in the approximation of infinitesimal amplitudes, in which skewness and re-ra 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(80)Distancel0l1Skewnessre-raC1–C20.05060.0527-29.5846-0.01967C1–H40.08980.0819-0.1031-0.0532C1⋯N30.06940.0754-0.4085-0.0242C1⋯O90.10170.11170.5720-0.2607,
and so forth. The exaggerated re-ra (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.



---

*Source: 864714-2011-12-28.xml*