About this Journal Submit a Manuscript Table of Contents
Advances in Physical Chemistry
Volume 2011 (2011), Article ID 864714, 14 pages
Research Article

Diffraction Measurements and Equilibrium Parameters


Received 6 June 2011; Accepted 20 September 2011

Academic Editor: Jan Skov Pedersen




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 OO 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 system 𝑅(𝑡) can in the first approximation be represented as𝑅(𝑡)=𝑅𝑒+𝑛𝑖=1𝑋𝑖𝜔cos𝑖𝑡+𝜃𝑖,(1) where 𝑅𝑒 is the vector of equilibrium geometric parameters, 𝜔𝑖 and 𝜃𝑖 are the frequencies and phases of vibrational motions, 𝑋𝑖 is the vector of geometric parameter changes at 𝜔𝑖𝑡+𝜃𝑖=2𝑘𝜋 and 𝜔𝑗𝑖𝑡+𝜃𝑗=(2𝑘+1)𝜋/2(𝑘=0,), and 𝑛 is the number of vibrational degrees of freedom. Since the ratios between the frequencies and phases of various vibrations must not necessarily be rational numbers, it is clear that a vibrational system has a chance to assume its equilibrium configuration only once during the whole time of its existence. This requires that all the 𝜔𝑖𝑡+𝜃𝑖 values simultaneously satisfy the condition 𝜔𝑖𝑡+𝜃𝑖=(2𝑘+1)𝜋/2(𝑘=0,). For this reason, “measurements’’ of equilibrium parameters are out of the question. It seems to follow from (1) that the mean (measured) geometric parameters of a vibrational system coincide with equilibrium. However, this is not the case even in the first approximation. Indeed, there is no one-to-one correspondence between the configuration determined by (1) and parameters measured experimentally, and the OO 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 𝑅 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𝑑𝑑𝑡𝜕𝜕̇𝑞𝜕𝜕𝑞=0(2) ( is the Lagrange function) written for the nuclear subsystem takes the form𝑑𝑑𝑡𝜕𝑇(𝑞,̇𝑞)𝜕𝜕̇𝑞𝜕𝑞(𝑇(𝑞,̇𝑞)𝑈(𝑞))=0(3) (𝑈 is assumed to be independent of the velocities of the particles constituting the system). Here, 𝑇(𝑞,̇𝑞) is the additive part of the Lagrange function (the kinetic energy), 𝑈(𝑞) is the potential energy (the nonadditive part), and 𝑞 stands for some “generalized’’ coordinates describing the configuration of the system.

The kinetic and potential energies are written as1𝑇=2𝑁𝑖=1𝑚𝑖̇𝑥2𝑖+̇𝑦2𝑖+̇𝑧2𝑖,1𝑈=23𝑁6(5)𝑖=1𝜕2𝑈𝜕𝑞𝑖𝜕𝑞𝑗𝑞𝑖𝑞𝑗+,(4) where 𝑁 is the number of vibrational system particles and 3𝑁6(5) is the number of vibrational degrees of freedom (all the 𝑞 and 𝑥 coordinates are taken to be zero in the equilibrium configuration). In the matrix form, 𝑇=(1/2)𝑀̇𝑥𝑖̇𝑥𝑗 and 𝑈=(1/2)𝐹𝑞𝑖𝑞𝑗. Here, 𝑀 is the diagonal matrix of atomic masses containing three masses of each atom corresponding to motions along the 𝑥, 𝑦, and 𝑧 axes and 𝐹 is the potential energy operator written in generalized coordinates 𝑞. The 𝐹 operator is symmetrical by construction. In the 𝑞=𝐵𝑥 generalized coordinates (𝐵 is the transformation matrix between the generalized and Cartesian coordinates), the kinetic energy takes the form1𝑇=2𝐵1̇𝑞𝑗𝑀𝐵1̇𝑞𝑖=𝐺1̇𝑞𝑖̇𝑞𝑗(5) ( is the symbol of transposition). The 𝐺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𝑞𝑖=3𝑁𝑗𝜕𝑞𝑖𝜕𝑥𝑗𝑥𝑗+123𝑁𝑗,𝑘𝜕2𝑞𝑖𝜕𝑥𝑗𝜕𝑥𝑘𝑥𝑗𝑥𝑘+,(6) that is, 𝐵=𝐵(0)+(1/2)(𝜕𝐵/𝜕𝑥), and1𝑈=23𝑁6(5)𝑖,𝑗𝜕2𝑈𝜕𝑞𝑖𝜕𝑞𝑗𝑞𝑖𝑞𝑗+163𝑁6(5)𝑖,𝑗,𝑘𝜕3𝑈𝜕𝑞𝑖𝜕𝑞𝑗𝜕𝑞𝑘𝑞𝑖𝑞𝑗𝑞𝑘+,(7) that is, 𝐹=𝐹(0)+(1/3)(𝜕𝐹/𝜕𝑞), allow (3) to be rewritten as𝐺1(0)̈𝑞𝑗+12𝜕𝐺1𝑞𝜕𝑞𝑖̈𝑞𝑗+̇𝑞𝑖̇𝑞𝑗12𝜕𝐺1𝜕𝑞̇𝑞𝑖̇𝑞𝑗+𝐹(0)𝑞𝑖+12𝜕𝐹𝑞𝜕𝑞𝑖𝑞𝑗=0,(8) where 𝐵(0), 𝐺1(0), and 𝐹(0) are the matrices corresponding to the approximation of infinitesimal vibrational amplitudes. The derivative of an 𝑚×𝑛 matrix 𝐴 with respect to an 𝑙-dimensional vector 𝑏 is understood to be a set of 𝑚×𝑛×𝑙(𝜕𝑎𝑖𝑗/𝜕𝑏𝑘) values. If the Cartesian system of coordinates is used, 𝐺1(𝑞)=𝐺1(0)=𝑀.

In the approximation of infinitesimal amplitudes ((𝜕𝐺1/𝜕𝑞)𝑞=0 and (𝜕𝐹/𝜕𝑞)𝑞=0), (8) is rewritten as𝐺1̈𝑞𝑗+𝐹𝑞𝑗=0.(9) Solutions to this system of equations are sought in the form 𝑞=𝑙cos(𝜔𝑡+𝜃) (̈𝑞=𝑙𝜔2cos(𝜔𝑡+𝜃)), that is,𝐺1𝑙𝛼𝜔2𝛼𝐹𝑙𝛼=0,𝐺𝐹𝑙𝛼=𝑙𝛼𝜔2𝛼,(10) where 𝑙𝛼 is the vector determining the change in the configuration of the vibrational system vibrating at the 𝜔𝛼 frequency (eigenvectors and frequencies are indexed by Greek letters). Since the 𝐺1 and 𝐹 matrices are symmetrical and the 𝐺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 𝐺1 matrix (and, simultaneously, to the 𝐹 matrix: 𝑉1𝐺1𝑉1=𝐷1 and 𝑉1𝐹𝑉1=𝐹1, where 𝐷1 is some diagonal matrix and 𝐹1 remains symmetrical (Hermitian) by construction. Since 𝐺1 is positive definite, all the elements of the 𝐷1 matrix are positive numbers. We can therefore construct the 𝐷11/2 diagonal matrix such that 𝐷11/2𝐷1𝐷11/2=𝐼 (𝐼 is the identity matrix). As a result, the kinetic energy matrix transforms into the identity matrix whereas the potential energy matrix 𝐷11/2𝐹1𝐷11/2=𝐹2 remains symmetrical. It can therefore be diagonalized with the use of a unitary (orthogonal) transformation, 𝑉2𝐷11/2𝑉1𝐹𝑉1𝐷11/2𝑉2=𝐷2. As to the identity matrix, its unitary transformation leaves it unchanged, 𝑉2𝐼𝑉2=𝐼.

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

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

Solutions to (8) will be sought in the form 𝑠=𝑠0+𝑠1 (𝑠0,𝛼=𝑒𝛼cos(𝜔𝛼𝑡+𝜃𝛼)) and 𝜔=𝜔0+𝜔1. Let us rewrite (9) and (8) using the notation introduced above and subtract the former from the latter, 𝐼̈𝑠𝑗1+Ω𝑠𝑗11=2𝜕𝔊1𝑠𝜕𝑠𝑖0̈𝑠𝑗0+̇𝑠𝑖0̇𝑠𝑗0+12𝜕𝔊1𝜕𝑠̇𝑠𝑖0̇𝑠𝑗012𝜕𝔉𝑠𝜕𝑠𝑖0𝑠𝑗0(13) (no corrections should be introduced into 𝑠𝑖𝑠𝑗-type products, because this would cause the appearance of terms smaller in magnitude than 𝑠1). Here, 𝑠0=3𝑁6(5)𝛼=1𝑒𝛼cos(𝜔𝛼𝑡+𝜃𝛼) is a linear combination of “unit’’ vectors 𝑒𝛼 (the superposition principle), which are also basis vectors for the construction of 𝑠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 𝐺 matrix. If the model of molecular motions with the 𝐹 matrix written in Cartesian coordinates is used, both these terms vanish. Of course, this cannot be balanced by any terms that appear in the expansion of the potential function into a series. It follows that the solution to the problem depends on the selected model of molecular motions.

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

Problem (13) can be divided into three independent problems:𝐼̈𝑠𝑗1,𝑘+Ω𝑠𝑗1,𝑘1=2𝜕𝔊1𝑠𝜕𝑠𝑖0̈𝑠𝑗0+̇𝑠𝑖0̇𝑠𝑗0,(14)𝐼̈𝑠𝑗1,𝑐+Ω𝑠𝑗1,𝑐=12𝜕𝔊1𝜕𝑠̇𝑠𝑖0̇𝑠𝑗0,(15)𝐼̈𝑠𝑗1,𝑎+Ω𝑠𝑗1,𝑎1=2𝜕𝔉𝑠𝜕𝑠𝑖0𝑠𝑗0.(16) Clearly, 𝑠1=𝑠1,𝑘+𝑠1,𝑐+𝑠1,𝑎 is the solution to (13). Let us consider problems (14)–(16) sequentially.

3. Kinematic Problem: Equation (14)

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

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

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

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

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

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

Figure 2: Bond rotation about a center.

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

Figure 3: Bond rotation about an axis.

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

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

For problem (15), the equation similar to (17) has the form𝐼̈𝑠1,𝑐,𝛼+Ω𝑠1,𝑐,𝛼=12𝜔𝛼𝜔sin𝛼𝑡+𝜃𝛼×3𝑁6(5)𝛽=1𝜔𝛽𝜔sin𝛽+𝜃𝛽Δ𝔊𝛽1𝑒𝛼(22) (the resonance term is excluded). Transforming the products of sine functions as in (18), we obtain𝐼̈𝑠1,𝑐,𝛼+Ω𝑠1,𝑐,𝛼1=4𝜔𝛼3𝑁6(5)𝛽=1𝜔𝛽𝜔cos𝛼+𝜔𝛽𝜃𝑡+𝛼+𝜃𝛽̂𝔤𝛽(𝛼)+14𝜔𝛼3𝑁6(5)𝛽=1𝜔𝛽𝜔cos𝛼𝜔𝛽𝜃𝑡+𝛼𝜃𝛽̂𝔤𝛽(𝛼).(23) Substituting (19) into the left-hand side of (23) yields𝑎𝛽1=4𝜔𝛼𝜔𝛽𝔓+(𝛽)̂𝔤𝛽(𝛼),𝑏𝛽1=4𝜔𝛼𝜔𝛽𝔓(𝛽)̂𝔤𝛽(𝛼)(24) (the notation is the same as before). It follows that𝑠1,𝑐,𝛼=143𝑁6(5)𝛽=1𝜔𝛼𝜔𝛽𝜔cos𝛼+𝜔𝛽𝜃𝑡+𝛼+𝜃𝛽𝔓+(𝛽)̂𝔤𝛽(𝛼)+𝜔𝛼𝜔𝛽𝜔cos𝛼𝜔𝛽𝜃𝑡+𝛼𝜃𝛽𝔓(𝛽)̂𝔤𝛽(𝛼).(25) Averaging over time and ensemble gives the contribution of the “centrifugal’’ term to the shrinkage effect,𝑠1,𝑐=143𝑁6(5)𝛼=1𝜔2𝛼𝔓(𝛼)̂𝔤𝛼(𝛼)=143𝑁6(5)𝛼=1𝜔2𝛼Ω1̂𝔤𝛼(𝛼)(26) (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’’ 𝜕𝔉/𝜕𝑠=. is a table (three-dimensional matrix) of the third derivatives of the potential energy with respect to the coordinates. In this system, the eigenvectors of the vibrational problem of the first approximation are “unit" vectors. Let us stage-by-stage reproduce solution to problem (15):𝐼̈𝑠1,𝑎,𝛼+Ω𝑠1,𝑎,𝛼1=2𝜔𝛼𝜔cos𝛼𝑡+𝜃𝛼3𝑁6(5)𝛽=1𝜔𝛽𝜔cos𝑏𝑡+𝜃𝛽𝑒𝛼𝑒𝛽1=43𝑁6(5)𝛽=1𝜔cos𝛼+𝜔𝛽+𝜃𝛼+𝜃𝛽𝑒𝛼𝑒𝛽𝜔+cos𝛼𝜔𝛽+𝜃𝛼𝜃𝛽×𝑒𝛼𝑒𝛽(27) (the resonance term is excluded). Substituting (19) into the left-hand side of this equation yields𝑎𝛽1=4𝔓+(𝛽)𝑒𝛼𝑒𝛽,𝑏𝛽1=4𝔓(𝛽)𝑒𝛼𝑒𝛽.(28) This gives𝑠1,𝑎,𝛼1=43𝑁6(5)𝛽=1𝜔cos𝛼+𝜔𝛽+𝜃𝛼+𝜃𝛽𝔓+(𝛽)𝜔+cos𝛼𝜔𝛽+𝜃𝛼𝜃𝛽×𝔓(𝛽)𝑒𝛼𝑒𝛽.(29) We find that the contribution of the anharmonic component to the shrinkage effect, which is determined by the term containing cos[(𝜔𝛼𝜔𝛼)+(𝜃𝛼𝜃𝛼)](𝛽=𝛼), can be written as𝑠1,𝑎,𝛼1=4Ω1𝑒𝛼𝑒𝛽.(30) The 𝑒𝛼𝑒𝛽 vectors “cut’’ a column with the ̂𝔥𝛼𝛽 elements from the matrix, and calculations present no difficulties.

Note that (21), (25), and (27) present the first members of the series describing atomic displacements from equilibrium positions. The convergence of these series depends, in particular, on the lengths of the ̂𝔤 and ̂𝔥 vectors, which are, in turn, determined by the lengths of “unit’’ vectors 𝑒 (𝑒𝛼=[(/4𝜋2𝑐𝜈𝛼)coth(𝑐𝜈𝛼/2𝑘𝑇)]1/2). These vectors become longer than one starting from 𝜔84.5cm1 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:𝑠𝛼=𝑠0,𝛼+1(𝑡)43𝑁6(5)𝛽=1𝜔𝛼𝜔𝛼+𝜔𝛽𝜔𝛼𝜔𝛽𝜔×cos𝛼+𝜔𝛽𝜃𝑡+𝛼+𝜃𝛽×𝔓+(𝛽)̂𝔤𝛽(𝛼)𝜔cos𝛼+𝜔𝛽+𝜃𝛼+𝜃𝛽×𝔓+(𝛽)̂𝔥𝛽(𝛼)+𝜔𝛼𝜔𝛼𝜔𝛽+𝜔𝛼𝜔𝛽𝜔×cos𝛼𝜔𝛽𝜃𝑡+𝛼𝜃𝛽𝔓(𝛽)̂𝔤𝛽(𝛼)𝜔cos𝛼𝜔𝛽+𝜃𝛼𝜃𝛽𝔓(𝛽)̂𝔥𝛽(𝛼)(31) for the vector of displacements with respect to the equilibrium configuration and𝑠𝛼1=4Ω1𝜔2𝛼̂𝔤𝛼(𝛼)̂𝔥𝛼(𝛼)(32) for the contribution of the vibration with the 𝜔𝛼 frequency to the shrinkage effect. Note that ̂𝔤𝛽(𝛼)=̂𝔤𝛼(𝛽) (these vectors only differ in the order of differentiation), and, naturally, ̂𝔥𝛽(𝛼)=̂𝔥𝛼(𝛽). For this reason, the summation in 𝛽 can be performed from 𝛽=1 to 𝛽=𝛼.

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

The solution becomes meaningless if (𝜔𝛼±𝜔𝛽)2𝜔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) coordinates 𝑞. For the transition to internuclear distances and the spatial configuration of the system, the results should be transformed into Cartesian coordinates. This transformation is performed as𝑥𝛼=𝐵1𝑙𝛼𝜔cos𝛼𝑡+𝜃𝛼,(35) or, according to (33),𝑥𝛼=𝐵1𝑙0,𝛼𝜔cos𝛼𝑡+𝜃𝛼+𝐵1𝐿𝑠1,𝛼.(36) Because of the smallness of the second term, we can assume 𝐵1=𝐵(0)1 in it and exclude it from consideration in this section. Since 𝐵1=𝐵(0)1+(1/2)3𝑁6(5)𝛽=1(𝜕𝐵1/𝜕𝑥)𝑥𝛽cos(𝜔𝛽𝑡+𝜃𝛽), we have𝑥𝛼tr=123𝑁6(5)𝛽=1Δ𝐵𝛽1𝑙0,𝛼𝜔cos𝛼𝑡+𝜃𝛼𝜔cos𝛽𝑡+𝜃𝛽,(37) where 𝑥tr is the displacement that appears because of nonlinearity of the transition between internal and Cartesian coordinates and Δ𝐵𝛽1 is the increment of the 𝐵1 matrix when the configuration of the system changes by the 𝑥𝛽 vector. In this sum, only the term with cos2(𝜔𝛼𝑡+𝜃𝛼) does not vanish in averaging. For this reason,𝑥𝛼tr=12𝜕𝐵1𝑥𝜕𝑥𝛼𝑙0,𝛼cos2𝜔𝛼𝑡+𝜃𝛼=14𝜎𝛼Δ𝐵𝛼1𝑙0,𝛼,(38) where 𝜎𝛼=(/4𝜋2𝑐𝜈𝛼)coth(𝑐𝜈𝛼/2𝑘𝑇) [6]. The 𝑥𝛼tr value can be obtained without comparatively laborious calculations of the Δ𝐵𝛼1 matrices. For this purpose, let us rewrite (38) as follows. Since 𝐵1=(𝐵(0)+(1/2)Δ𝐵)1=(𝐼+(1/2)𝐵(0)1Δ𝐵)1𝐵(0)1(𝐼(1/2)𝐵(0)1Δ𝐵)𝐵(0)1, we have Δ𝐵1𝐵(0)1Δ𝐵𝐵(0)1. On the other hand, 𝐵(0)1𝑙0,𝛼=𝑥0,𝛼. Therefore,𝑥𝛼tr14𝜎𝛼𝐵(0)1Δ𝐵𝛼𝐵(0)1𝑙0,𝛼=14𝜎𝛼𝐵(0)1Δ𝐵𝛼𝑥0,𝛼.(39) Here, 𝑥0,𝛼 is the eigenvector in Cartesian coordinates obtained by solving the vibrational problem in the first approximation. Let us calculate the 𝑙+0,𝛼 and 𝑙0,𝛼 vectors for the 𝑅0+𝑥0,𝛼 and 𝑅0𝑥0,𝛼 configurations, respectively (𝑅0 is the vector of the Cartesian coordinates of atoms in the equilibrium configuration). Clearly,𝑙+0,𝛼=1𝐵(0)+2Δ𝐵𝛼𝑥0,𝛼,𝑙0,𝛼=1𝐵(0)2Δ𝐵𝛼𝑥0,𝛼.(40) Summing these equations yields𝑙+0,𝛼+𝑙0,𝛼=Δ𝐵𝛼𝑥0,𝛼.(41) Substituting the result into (39) yields𝑥𝛼tr=14𝜎𝛼𝐵(0)1𝑙+𝛼+𝑙𝛼.(42) On the other hand, the 𝑥+𝛼 and 𝑥𝛼 vectors satisfying the conditions 𝑙+𝛼=𝐵𝑥+𝛼 and 𝑙𝛼=𝐵𝑥𝛼 can be calculated by the method of successive approximations by fitting the 𝐵(𝑥) matrix and 𝑥 vector components to obtain 𝑞𝛼=𝐵(𝑥)𝑥𝛼. Clearly, 𝑥𝛼=(𝑥+𝛼+𝑥𝛼)/2. This method gives somewhat more accurate results, especially when vibrational amplitudes are indeed large. It is implemented in the Shrink09 program [9].

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

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

On the other hand, for the 𝑅𝑖𝑗 distance between atoms 𝑖 and 𝑗, we haveΔ𝑅𝑖𝑗=3𝑁6(5)𝛼=1𝑋𝑖𝑋𝑗+𝜉𝛼,𝑖𝜉𝛼,𝑗2+𝑌𝑖𝑌𝑗+𝜐𝛼,𝑖𝜐𝛼,𝑗2+𝑍𝑖𝑍𝑗+𝜁𝛼,𝑖𝜁𝛼,𝑗21/2𝑅0,𝑖𝑗,(46) where 𝑅0,𝑖𝑗 is the distance between atoms 𝑖 and 𝑗 in the equilibrium configuration. The Δ𝑅𝑖𝑗 values determine the shrinkage effect for internuclear distances; in diffraction experiments, we measure the 𝑅0,𝑖𝑗+Δ𝑅𝑖𝑗 distances, and, for the transition to the equilibrium configuration, the corresponding corrections should be introduced into these values.

Two other important parameters are mean square amplitudes of changes in internuclear distances and skewness. Let us rewrite (37) in the form𝑥𝛼tr=143𝑁6(5)𝛽=1Δ𝐵𝛽1𝑙0,𝛼𝜔cos𝛼+𝜔𝛽𝜃𝑡+𝛼+𝜃𝛽𝜔+cos𝛼𝜔𝛽𝜃𝑡+𝛼𝜃𝛽.(47) To calculate the mean square amplitudes of the displacement of atoms from equilibrium positions, it is necessary: (i) to perform the transition to the Cartesian coordinates in (31) through multiplying by the 𝐵(0)1𝐿 matrix (the 𝑠0,𝛼(𝑡) value then transforms into 𝑥0,𝛼cos(𝜔𝛼𝑡+𝜃𝛼)), (ii) sum (31) and (47), (iii) raise the result to a power of two, and (iv) perform averaging over time and phases (the latter for degenerate vibrations). In averaging, only 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 𝑥2𝛼 vector,𝑥2𝛼=𝜎𝛼𝑥Sq0,𝛼+123𝑁6(5)𝛽=1𝑎Sq+𝛽𝑎+Sq𝛽(48) (the 1/2 coefficient appears in the averaging of squared cosine functions over time and phases). Here the Sq symbol denotes the multiplication by the diagonal matrix whose components are the components of the vector following this symbol (squaring of vector components), and the 𝑎+𝛽 and 𝑎𝛽 vectors are𝑎+𝛽=14𝜔𝛼𝜔𝛼+𝜔𝛽𝜔𝛼𝜔𝛽𝐵(0)1𝐿𝔓+(𝛽)̂𝔤𝛽(𝛼)+𝐵(0)1𝐿𝔓+(𝛽)̂𝔥𝛽(𝛼)+12𝜎𝛼𝜎𝛽Δ𝐵𝛽1𝑙0,𝛼,𝑎𝛽=14𝜔𝛼𝜔𝛼𝜔𝛽+𝜔𝛼𝜔𝛽𝐵(0)1𝐿𝔓(𝛽)̂𝔤𝛽(𝛼)+𝐵(0)1𝐿𝔓(𝛽)̂𝔥𝛽(𝛼)+12𝜎𝛼𝜎𝛽Δ𝐵𝛽1𝑙0,𝛼.(49) The mean-square amplitudes of changes in internuclear distances 𝑥2𝑖𝑗 are given by𝑥2𝑖𝑗=123𝑁6(5)𝛼=1𝑥𝛼,𝑖,𝑥𝑥𝛼,𝑗,𝑥2+𝑥𝛼,𝑖,𝑦𝑥𝛼,𝑗,𝑦2+𝑥𝛼,𝑖,𝑧𝑥𝛼,𝑗,𝑧2,(50) where 𝑥𝛼,𝑖,𝑥, 𝑥𝛼,𝑖,𝑦, and 𝑥𝛼,𝑖,𝑧 are the components of the 𝑥0,𝛼+3𝑁6(5)𝛽=1(𝑎+𝛽+𝑎𝛽) vector corresponding to the displacements of the 𝑖th atom along the 𝑥, 𝑦, and 𝑧 axes; again, cos2(𝑎𝑡+𝑏) gives the coefficient 1/2. To pass to the central moment from the moment about zero, we must subtract 𝑅𝛼,𝑖𝑗2 from 𝑥2𝛼,𝑖𝑗. This is equivalent to excluding terms with cos(𝛼𝛼) from amplitude calculations.

If the Δ𝑅𝛼,𝑖𝑗 and 𝑥2𝛼,𝑖𝑗 values are known, asymmetry parameter 𝑥3𝑖𝑗/𝑥2𝑖𝑗3/2 calculations present no difficulties. Indeed,𝑥3𝑖𝑗=3𝑁6(5)𝛼=1Δ𝑅3𝛼,𝑖𝑗.(51)

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

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

7.1. Errors in Cubic Force Constant Calculations

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

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

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

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

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

The following procedure for estimating calculation errors can be suggested. Let index 𝛼 run over all the permutations of the 𝑖,𝑗,𝑘 indices. The mean 𝐶𝛼 value is then𝐶𝛼=16𝛼𝐶𝛼.(58) Let us putΔ𝛼=16𝛼𝐶𝛼𝐶𝛼21/2,Δ𝛼=max𝐶𝛼min𝐶𝛼2.(59) In other words, Δ is the root-mean-square deviation for the given 𝛼 (for the given set of 𝑖, 𝑗, and 𝑘 indices), and Δ is the half-spread of the 𝐶𝛼 values. We putΔ=maxΔ𝛼,Δ=maxΔ𝛼,(60) where maxima are sought over all 𝛼 (over all different sets of 𝑖, 𝑗, and 𝑘 indices). Δ can be suggested as an estimate of errors in off-diagonal 𝐻𝐶 matrix elements, and Δ, as a (conventional) estimate of errors in 𝐶𝑖𝑖𝑖. The Gaussian program cannot be used to obtain the required information. The corresponding calculations were therefore performed using the Gamess [12] package and a special program external with respect to Gamess (the Gamess package is not intended for cubic constant calculations). Calculations for three molecules gave Δ=0.825737×104 au and Δ=0.984625×104 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𝜕𝐹𝜕𝑞𝑖𝑗𝑘=𝐻𝐶𝑖𝑗𝑘𝐵𝑖𝑗𝐵𝑗𝑘𝐵𝑘𝑗1,(61) where (𝜕𝐹/𝜕𝑞)𝑖𝑗𝑘 is the three-dimensional matrix of the third derivatives of potential energy with respect to internal coordinates, 𝐻𝑖𝑗𝑘 is, as before, the three-dimensional matrix of the third derivatives of potential energy with respect to Cartesian coordinates, and 𝐵𝑗𝑖 is the transformation matrix between internal and Cartesian coordinates, 𝑞𝑗=𝐵𝑗𝑖𝑥𝑖 (the matrix of second derivatives can be transformed this way only because the corresponding gradient vector is zero at equilibrium).

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

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

8. Scaling of Cubic Force Constants

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

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

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

9. The Problem of Low Frequencies

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

10. Internal Rotations and Similar Problems

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

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

Let us turn to some obvious examples. Let us, for instance, consider antisymmetric stretching vibration of an 𝐴𝐵2 triatomic molecule. Clearly, the 𝐵𝐴𝐵 angle should change along the corresponding minimum energy path. But the angular coordinate transforms under the 𝐴1 representation, whereas the antisymmetric stretching coordinate, under the 𝐵1 (or 𝐵2) 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 Figure 5.

Figure 5: Potential energy curve along a torsional coordinate.

The horizontal lines in Figure 5 are vibrational levels, and the potential well contains five of them. It is always possible to solve the vibrational problem for five levels. We can determine their populations and select a coefficient for recalculating the corresponding eigenvector components into the Cartesian displacements of atoms. Such calculations actually divide the system into two subsystems with known populations. For instance, let us consider the nitroethane molecule (Figure 6), which has a very low torsional frequency corresponding to 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,1𝑉(𝜙)=2𝑚=1𝑉𝑚𝑠[],1cos(𝑚𝑠𝜙)(75) but this is an even function of the torsional angle 𝜙, which cannot contribute to the shrinkage effect. As to the amplitude of this vibration, it can be measured directly as the distance between potential barriers obtained in quantum-mechanical calculations.

11. The Problem of Redundant Coordinates

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

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

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

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

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

12. Conclusions

The calculations described above solve the vibrational problem and the problem of the search for parameters necessary for the interpretation of diffraction data at the level of first-order perturbation theory. The solution was obtained using the classical formalism, the harmonic approximation serving as a zero order step. The paper generalizes the results obtained in the preceding works [19, 2226]. 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):Distance01Skewness𝑟𝑒𝑟𝑎C1C20.05060.05240.00630.0196C1H40.07680.07760.00250.0145C1N30.06940.07710.03540.0041C1O90.10170.12740.35120.1605,(79) and so forth (here, 𝑟𝑒 and 𝑟𝑎 are the equilibrium and experimentally observed internuclear distances; column 0 contains amplitudes calculated in the approximation of infinitesimal amplitudes, in which skewness and 𝑟𝑒𝑟𝑎 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, obtainDistance01Skewness𝑟𝑒𝑟𝑎C1C20.05060.052729.58460.01967C1H40.08980.08190.10310.0532C1N30.06940.07540.40850.0242C1O90.10170.11170.57200.2607,(80) and so forth. The exaggerated 𝑟𝑒𝑟𝑎 (and, therefore, skewness) values can hardly be considered realistic.

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


  1. O. Bastiansen and M. Traettebrg, “The influence of thermal motion on structure determination of linear molecules using the electron diffraction method,” Acta Crystallographica, vol. 13, p. 1108, 1960. View at Google Scholar
  2. Y. Morino, “The effect of thermal vibrations of gaseous molecules in electron diffraction studies,” Acta Crystallographica, vol. 13, p. 1107, 1960. View at Google Scholar
  3. A. Almenningen, S. P. Arnesen, O. Bastiansen, H. M. Seip, and R. Seip, “The effect of temperature variatin on the amplitudes of vibration and shrinkage effects in carbon suboxide studied by gas electron diffraction,” Chemical Physics Letters, vol. 1, no. 12, pp. 597–598, 1968. View at Publisher · View at Google Scholar · View at Scopus
  4. M. Tanimoto, K. Kuchitsu, and Y. Morino, “Molecular structure and the effect of large-amplitude vibration of carbon suboxide as studied by gas electron diffraction,” Bulletin of the Chemical Society of Japan, vol. 43, no. 9, pp. 2776–2785, 1970. View at Publisher · View at Google Scholar
  5. A. Clark and H. M. Seip, “The potential function for the CCC bending in carbon suboxide,” Chemical Physics Letters, vol. 6, no. 5, pp. 452–456, 1970. View at Publisher · View at Google Scholar
  6. S. J. Cyvin, Molecular Vibrations and Mean Square Amplitudes, Elsevier, Amsterdam, The Netherlands, 1968.
  7. L. D. Landau and E. M. Lifshits, Mekhanika (Mechanics), Fizmatgiz, Moscow, Russia, 1958.
  8. L. S. Bartell, “Calculation of mean atomic positions in vibrating polyatomic molecules,” Journal of Chemical Physics, vol. 38, p. 1827, 1963. View at Google Scholar
  9. The program is available from the author.
  10. M. Iwasaki and K. Hedberg, “Centrifugal distortion of bond distances and bond angles,” The Journal of Chemical Physics, vol. 36, no. 11, pp. 2961–2963, 1962. View at Publisher · View at Google Scholar · View at Scopus
  11. Ph. E. Gill, W. Murray, and M. H. Wright, Practical Optimization, Academic, New York, NY, USA, 1981.
  12. A. A. Granovskii, http://classic.chem.msu.su/gran/gamess/index.html.
  13. I. S. Sokolnikoff, Tensor Analysis, John Wiley & Sons, New York, NY, USA, 1951.
  14. P. Pulay, G. Fogarasi, F. Pang, and J. E. Boggs, “Systematic ab initio gradient calculation of molecular geometries, force constants, and dipole moment derivatives,” Journal of the American Chemical Society, vol. 101, no. 10, pp. 2550–2560, 1979. View at Publisher · View at Google Scholar · View at Scopus
  15. I. V. Kochikov, Y. I. Tarasov, V. P. Spiridonov, et al., “The equilibrium structure of thiophene by the combined use of electron diffraction, vibrational spectroscopy, and microwave spectroscopy guided by theoretical calculations,” Journal of Molecular Structure, vol. 567-568, p. 29, 2001. View at Google Scholar · View at Scopus
  16. S. Kondo, Y. Koga, and T. Nakanaga, “Ab initio MO calculation of the anharmonic force field of methyl fluoride and methyl chloride,” The Journal of Chemical Physics, vol. 81, no. 4, pp. 1951–1959, 1984. View at Publisher · View at Google Scholar · View at Scopus
  17. S. Kondo, “Empirical improvement of the anharmonic ab initio force field of methyl fluoride,” The Journal of Chemical Physics, vol. 81, no. 12, pp. 5945–5951, 1984. View at Publisher · View at Google Scholar · View at Scopus
  18. G. Klatt, A. Willetts, N. C. Handy, R. Tarroni, and P. Palmieri, “An improved anharmonic potential for CHF3,” Journal of Molecular Spectroscopy, vol. 176, no. 1, pp. 64–74, 1996. View at Publisher · View at Google Scholar · View at Scopus
  19. V. A. Sipachev, “The use of quantum-mechanical third-order force constants in structural studies,” Journal of Molecular Structure, vol. 693, no. 1–3, pp. 235–240, 2004. View at Publisher · View at Google Scholar · View at Scopus
  20. Y. I. Tarasov, I. V. Kochikov, N. Vogt et al., “Electron diffraction and quantum chemical study of the structure and internal rotation in nitroethane,” Journal of Molecular Structure, vol. 872, no. 2-3, pp. 150–165, 2008. View at Publisher · View at Google Scholar · View at Scopus
  21. V. I. Arnol’d, V. V. Kozlov, and A. I. Neishtadt, “Mathematical aspects of classical and celestial mechanics,” in Advances in Science and Technology. Modern Problems of Mathematics. Basic Research, V. I. Arnol’d, Ed., Dynamic Systems-3, VINITI, Moscow, Russia, 1985. View at Google Scholar
  22. V. A. Sipachev, “Calculation of shrinkage corrections in harmonic approximation,” Journal of Molecular Structure, vol. 22, no. 1–5, pp. 143–151, 1985. View at Publisher · View at Google Scholar · View at Scopus
  23. V. A. Sipachev, “Vibrational effects in diffraction and microwave experiments: a start on the problem,” in Advances in Molecular Structure Research, I. Hargittai and M. Hargittai, Eds., vol. 5, JAI, Greenwich, UK, 1999. View at Google Scholar
  24. V. A. Sipachev, “Anharmonic corrections to structural experiment data,” Structural Chemistry, vol. 11, no. 2-3, pp. 167–172, 2000. View at Publisher · View at Google Scholar · View at Scopus
  25. V. A. Sipachev, “Local centrifugal distortions caused by internal motions of molecules,” Journal of Molecular Structure, vol. 567-568, pp. 67–72, 2001. View at Publisher · View at Google Scholar · View at Scopus
  26. V. A. Sipachev, Strength from Weakness: Structural Consequences of Weak Interactions in Molecules, Supermolecules, and Crystals, Kluwer, Dordrecht, The Netherlands, 2002, Edited by A. Domenicano and I. Hargittai.
  27. A. Tulinsky and C. R. Worthington, “Basic beryllium acetate. I. The collection of intensity data,” Acta Crystallographica, vol. 12, pp. 623–626, 1959. View at Google Scholar