Diffuse X-ray scattering from correlated motions in a protein crystal

Protein dynamics are integral to biological function, yet few techniques are sensitive to collective atomic motions. A long-standing goal of X-ray crystallography has been to combine structural information from Bragg diffraction with dynamic information contained in the diffuse scattering background. However, the origin of macromolecular diffuse scattering has been poorly understood, limiting its applicability. We present a finely sampled diffuse scattering map from triclinic lysozyme with unprecedented accuracy and detail, clearly resolving both the inter- and intramolecular correlations. These correlations are studied theoretically using both all-atom molecular dynamics and simple vibrational models. Although lattice dynamics reproduce most of the diffuse pattern, protein internal dynamics, which include hinge-bending motions, are needed to explain the short-ranged correlations revealed by Patterson analysis. These insights lay the groundwork for animating crystal structures with biochemically relevant motions.

To prevent dehydration at ambient temperature, each crystal was enclosed in a plastic capillary with crystallization solution in the tip. X-ray diffraction images were recorded by a photon-counting detector using shutterless data collection with fine phi-slicing (0.1 • per frame). (B) The total dose was distributed over 11 diffraction volumes from four large crystals (Supplementary Table 1). (C) At each spindle angle (phi), diffraction images were collected with the crystal in the beam (left) and with the crystal translated out of the beam (middle). As the background includes scattering from the instrument and narrow rings due to the plastic capillary, the pattern due to the sample only is estimated by subtraction (right). For clear illustration, the images shown are the sums of 10 sequential frames (1 • oscillation and 1 s exposure in total). Fig. 3: Rocking curve for the (2,-8,-1) reflection. In our diffraction data from triclinic lysozyme, the apparent crystal mosaicity of 0.02 to 0.03 • is smaller than the oscillation width of 0.1 • per frame, so most Bragg reflections appear only in one or two consecutive images and their rocking curves cannot be resolved. We thus identified an intense (2,-8,-1) Bragg reflection that is close to the spindle axis and passes through the Ewald sphere along c * at a rate approximately 3 times slower than the fastest route. The peak profile above the background was well-fit by a three-dimensional Gaussian model [2] with apparent mosaicity and beam divergence of 0.020 • and 0.024 • , respectively. (A) The total photon counts in a 25-pixel region around the Bragg peak is plotted vs. spindle angle from the diffraction maximum (black symbols). The sharp Bragg peak (Gaussian fit shown in blue) is easily distinguished from the more slowly-varying background. The vertical dashed lines show the voxel subdivisions along c * used in the construction of the fine diffuse map. The peak is confined to the central subdivision. (B) On the left is a zoom-in of the (2,-8,-1) Bragg reflection in the diffraction image (0.1 • oscillation) at its peak intensity. The x and y axes are distances from the beam center, and the spindle is parallel to the x-axis. On the right is the same oscillation simulated using the Gaussian model. In both images, the intersection of voxel edges in the a * -b * plane with the Ewald sphere are projected onto the detector plane, showing that the Bragg intensity is confined to the central voxel.  To improve the signal-to-noise ratio, the interpolation was performed by least-squares fitting a 2nd order polynomial to the 5x5x5 voxel region around each point in the target grid. The resulting "half-integer map" (gray points) was analyzed as a function of the resolution d. Within each resolution bin, the isotropic background (red curve) was defined as 1σ below the mean (dashed black curve). (B) The isotropic background was subtracted from the total intensity to aid in visualizing the cloudy variational pattern. Slices from left to right are perpendicular to a * , b * , and c * . (C-D) The 1-unit cell MD simulation was analyzed in the same way as the half-integer map. The isotropic scattering and variational scattering components are similar between the 1-unit cell MD simulation and the half-integer map.

Supplementary
Supplementary Fig. 10: Diffuse streaks attributed to intense features in the molecular transform. In the experimental maps of variational scattering (left panel), certain pairs of halos are connected by streaks of diffuse intensity, such as (0,-5,12) and (0,-6,12), while others are not, such as (0,-4,12) and (0,-5,12). Similar features are seen in the one-phonon scattering simulation of the lattice dynamics model (center panel). These anisotropic features occur because the halo intensities are modulated by the molecular transform, the continuous diffraction intensity of an isolated unit cell (right); the molecular transform is intense between (0,-5,12) and (0,-6,12) but is close to zero between (0,-4,12) and (0,-5,12), explaining the presence of a diffuse streak in the former case and its absence in the latter.
Supplementary Fig. 11: Contribution of halos to the total variational scattering. The mean variational intensity in each resolution shell was calculated for the full map, which includes the halo features, and for the half-integer interpolated map (described in Supplementary Fig. 9A), which does not. The halo contribution was estimated as the difference between the two (shaded region). In most resolution shells, the halos account for about half of the total variational intensity (solid curve).
Supplementary Fig. 12: Correlations between Bragg and diffuse intensities. The total variational intensity, which includes halos and cloudy scattering, was averaged in neighborhoods (h,k,l) around each Bragg peak . The CC is strong at low resolution, but decays significantly at high resolution. In contrast, the CC between the experiment and the refined crystal structure model (black curve) remains excellent at high resolution.
Supplementary Fig. 13: Refinement of lattice dynamics model against a small subset of diffuse halos. (A) The spring constants in the lattice dynamics model were refined iteratively to minimize the least-squares difference (reduced χ 2 , blue symbols) between the calculated one-phonon intensity and the measured variational intensity around 400 intense halos. To avoid over-fitting, the refinement was carried out in four stages (i-iv) as restraints were removed (see Methods in the Main Text). The first stage (i) is a Gaussian model with equal springs, while in the final stage (iv) all springs are refined individually. (B) The ability of the model to explain halo anisotropy (see Fig. 3C in the Main Text) was assessed using the Pearson correlation coefficient for the 400 halos (CC aniso ). Initially, the anisotropy parameters for model and experiment show little correlation (CC aniso = 0.45, scatter plot in lower left). After refinement, the correlation improves (CC aniso = 0.92, scatter plot in lower right).
Supplementary Fig. 14: Correlation of experiment and simulated one-phonon scattering from the refined lattice dynamics model calculated with intensities from the full 13x11x11 grid. The data-model correlation (solid blue curve) is excellent (CC ∼ 0.9) between 2 and 5Å resolution (solid), but decays rapidly beyond ∼ 2Å resolution. At high resolution, the correlation is limited by noise in the experimental map, as estimated using the CC* statistic (dashed curve). Interpolating the maps on a 7x7x7 grid improves the signal-to-noise (see Fig. 2D in the Main Text).      Two elastic network models for internal protein motion were fit to the residual ADPs. The number of normal modes for the elastic network is 129 × 6 − 6 = 771, since there are 129 rigid groups (one per residue) with 6 degrees of freedom each, and the restraint that there is no rigid-body motion of the entire protein reduces the number of modes by 6. In the model with domain motion suppressed, similar restraints applied to the 3 rigid groups (α, β , and hinge) result in 129 × 6 − 3 × 6 = 756 normal modes. Thus, by restraining 3 groups instead of 1, the number of normal modes is very moderately reduced, which in principle might make it more difficult to fit the ADPs. However, we did not adjust the normal mode amplitudes to fit the ADPs as was done in prior studies [7,8]. Instead, we refined the physical model (spring constants) through a conservative parameterization with one coupling constant per residue (or 129 free parameters, see Equation 12 in the Main Text). Thus, the ability of each model to fit the ADPs is not limited by the number of normal modes, but is instead limited primarily by the parameterization, which is the same in both cases.

Supplementary Methods 1 General Theory
The kinematical theory of X-ray diffraction applicable to diffuse scattering from proteins has been described previously [9,10]. Here we provide an overview of the main results as relevant to this study.

Total scattering cross-section
The macroscopic differential scattering cross-section, dΣ/dΩ, describes the probability of observing a scattering event in the solid angle dΩ when the sample is illuminated by X-rays. The cross-section depends on the respective wavevectors of the incident and scattered X-rays, s 0 and s . The wavevectors have magnitudes of inverse wavelength: s 0 = |s 0 | = 1/λ and s = |s | = 1/λ , where λ and λ are the wavelengths of the incident and scattered X-ray photons. Given an incident flux, J 0 , the total scattered flux into a solid angle ∆Ω in the directionŝ is as follows: In the X-ray experiments described here, both coherent and incoherent processes contribute to X-ray scattering: The coherent X-ray scattering has a cross-section where s = R T (s − s 0 ) is the scattering vector in the sample frame which is rotated with respect to the lab frame by a rotation matrix R, I sample is the coherent intensity of the sample in per-electron units, and I e is the Thomson scattering cross-section of a free electron, I e = r 2 e P(ŝ 0 ,ŝ ), where r e is the classical electron radius and P(ŝ 0 ,ŝ ) is the polarization factor. Incoherent, or Compton, X-ray scattering has a cross-section [11]: where S sample is the sum of atomic incoherent scattering functions over all illuminated atoms [12]. The Compton scattering process is inelastic: the scattered photon has less energy than the incident photon by an amount that increases with scattering angle. However, the wavelength shift, which is at most twice the Compton wavelength of the electron λ C = hm −1 e c −1 ≈ 0.02426Å, is small compared with the λ ∼ 1Å wavelengths generally employed in X-ray crystallography. Therefore, in the following, we assume λ = λ .
The sample in this case is a crystal composed of unit cells. The crystal is assumed to be homogeneous, meaning that the unit cells are equivalent to each other on average, and therefore one can write: and where V 0 is the volume of the sample illuminated by the X-ray beam, v c is the unit cell volume (V 0 /v c is the number of unit cells illuminated), and I(s) and S(s) are the coherent and incoherent intensities per unit cell. Finally, we can write the scattered flux in terms of the intensities per unit cell by combining Equations 1-7, as follows:

Coherent scattering from crystals
Coherent scattering of X-rays depends on the spatial distribution of electrons in the sample. If ρ(r) is the instantaneous number density of electrons at a position r, the scattering is described by the Fourier transform of the electron density, a quantity known as the structure factor, F(s): The intensity per unit cell is proportional to the mean square structure factor, as follows: where N is the number of unit cells under consideration and brackets denote the ensemble (or time) average. In general, the electron density is time-dependent, but the system is in equilibrium so that we can define an average electron density ρ(r) with corresponding structure factor F(s) . If the structure factor is separated into the average and fluctuating terms, F = F +F, the total intensity (Equation 10) becomes We can further define an average unit-cell electron density, in which the densities of all unit cells are superimposed, ρ cell (r) , and corresponding structure factor F cell (s) . If the unit cells under consideration are equivalent, ensemble averaging is equivalent to spatial averaging, and the two average structure factors are related as follows: where the sum runs over all unit cells, N, and r n is the origin of each unit cell. If s happens to coincide with a node of the reciprocal lattice g h (h is the vector of Miller indices, h, k, and l), the phase factor in the sum is equal to 1 for all n, and the sum evaluates to N. When N is large, the amplitude of the lattice term falls off very quickly with |s − g h | , and as N → ∞, it can be written as a sum of delta functions as follows: where v * c is the volume of a unit cell of the reciprocal lattice. Thus, for the large crystal, we can combine Equations 11-13 to write the intensity per unit cell as follows: The second term of Equation 14 is the diffuse scattering per unit cell, The first term of Equation 14 contains the Bragg peaks. We can define the Bragg intensity as follows: which is obtained for s = g h upon integration of the Bragg peaks in reciprocal space. With these definitions (Equations 15-16), the coherent intensity from a macroscopic crystal (Equation 14) is as follows:

Diffuse Patterson function
The diffuse Patterson function, or three-dimensional difference pair distance distribution function (3D-∆PDF), is defined as the mean autocorrelation of the electron density fluctuations per unit cell, as follows: where ∆ρ = ρ − ρ is the difference between the instantaneous electron density and the ensemble average. The diffuse intensity per unit cell (Equation 15) can be rearranged as: Plugging in the definition of the structure factor (Equation 9), we see that I D is the Fourier transform of the ∆PDF, as follows: Thus, I D and ∆PDF are a Fourier transform pair, and the inverse transform is:

Scattering simulation
In order to simulate total scattering without having to model the entire crystal, it is useful to simulate a small number of unit cells and to impose periodic boundary conditions by what is known as the supercell method [13][14][15]. Given a supercell consisting of a three-dimensional array of N scell = N 1 × N 2 × N 3 unit cells, a "supercell reciprocal lattice" can be defined, which subdivides the regular reciprocal lattice by N 1 along a * , N 2 along b * , and N 3 along c * . The total intensity per unit cell (Equation 10) is calculated in the same manner as usual (Equation 11), except that the structure factor F scell is now a periodic function whose Fourier transform is defined only at nodes of the reciprocal lattice of the supercell. The mean of F scell , evaluated at the Bragg peak locations (the nodes of the unit cell reciprocal lattice, g h ) is proportional to the regular unit cell structure factor F cell : Thus, we can calculate the Bragg intensity per unit cell (Equation 16) from a supercell simulation as follows: where F scell is evaluated only at nodes of the regular, unit-cell reciprocal lattice. The diffuse intensity per unit cell is defined in the same way as the macroscopic crystal case (Equation 15), except that it is no longer a continuous function of s and is evaluated only at the nodes of the supercell reciprocal lattice.

Experimental Corrections
Here we consider geometric and experimental corrections that relate the photon counts observed by a detector to the coherent intensity per unit cell, I(s) (see Section 1 and Equation 10). The first section describes the factors depending on the scattering geometry. The second section describes corrections to the Bragg intensity for a crystal that rotates continuously during the exposure. The final section describes a method for placing the total intensity on an absolute scale and subtracting the incoherent contribution.

Geometric factors
We consider corrections applicable to the geometry normally used in macromolecular crystallography, which consists of a planar X-ray detector with square pixels and a sample that is rotated about a fixed spindle axis. Corrections are derived assuming that detection occurs in the far field from a small sample, so that the illuminated portion of the sample acts as a point-like source. In Section 1, Equation 8 relates the scattered flux J(ŝ ) in the directionŝ to the coherent and incoherent intensities. Two prefactors in that equation, the polarization P and solid angle ∆Ω, depend on the scattering geometry. The solid angle of the pixel can be calculated as follows: where d is the distance between the sample and detector pixel, ∆x and ∆y are the pixel width and height, and ω is the angle between the detector surface normal and the path of the incident photon. The polarization factor depends on the scattering geometry as well as X-ray source properties: the fraction of polarization, p, and the vector normal to the plane of polarization,n. In terms of these parameters, the polarization factor can be written as follows [16]: In addition to solid angle and polarization, we must also correct for the efficiency of detection. In an ideal situation, the expected number of photons recorded by a detector pixel, n, would be equal to the scattered flux (Equation 8) times the exposure time: n = J∆t. However, two non-ideal effects may be important. First is the possibility that the detector pixel does not absorb all of the photons it intercepts. For a pixel of thickness δ and absorption coefficient of the sensor κ, the fraction of photons absorbed, E, is as follows: Second is the possibility that the material between the sample and detector absorbs or scatters some of the X-rays before they can be detected. For example, if a uniform material (such as air) with attenuation coefficient µ occupies the space between the sample and the detector (the path length is d), the probability that a photon is transmitted, A, is as follows: The sample rotates during the exposure, so that each pixel samples a range of scattering vectors in the reference frame of the sample. Letm be the spindle axis, φ the spindle angle, ∆t the exposure time, and ∆φ /∆t the rate of rotation. Then, for small ∆φ and ∆Ω, the volume element swept by the pixel in reciprocal space is as follows [10]: The measured intensity thus depends on the average flux (Equation 8) for this volume, where the brackets signify averaging over the volume element ∆v * , as follows: Finally, there may be other photons detected during the exposure that come from background sources, such as air scatter. The background scattering rate for each pixel, r b , will be sample and instrument-dependent, but it can be estimated experimentally, for example by performing a measurement with the sample removed.
With all corrections put together, the expected number of photons detected is related to the intensity per unit cell as follows:

Corrections to experimental data
Equation 31 relates the measured photon counts to the signal of interest: the coherent intensity per unit cell, I, which consists of the Bragg and diffuse components. Here we describe the data processing steps that are necessary to convert photon counts into estimates of I B and I D on an absolute scale. In a standard macromolecular crystallography experiment, intensity is determined on an arbitrary scale (this is because the incident flux, the illuminated volume, or both, are typically unknown), and the background scattering is not measured. This is adequate for Bragg intensity determination because the intensity is integrated relative to a local background, and the overall scale factor is fit during model refinement. For diffuse scattering, it is necessary to measure and subtract the background. As a first step in data analysis, the background and geometric corrections are applied to the measured photon counts, so that one obtains an intensity on an arbitrary scale, which we call I meas , as follows: where the factors in the denominator are defined in Equations 24-27. This intensity is proportional to true intensity per unit cell, with the expected theoretical relationship (Equations 31and 32): For a given exposure, the prefactor J 0 r 2 e V 0 /v c is a constant that is the same for all pixels. The next step is to separate I meas into Bragg and continuous scattering components. This is simplest in the case where the coherent intensity I has a diffuse component that varies gradually on the scale of ∆v * and Bragg peaks that are delta function-like on this scale (Equation 17). Then, I meas can be treated differently depending on whether the integration volume ∆v * contains a Bragg peak. If it does not, the intensity can be factored out of the integral (Equation 30), and we have I + S ∆v * ≈ I D + S. If the volume contains a Bragg peak (at s = g h ), we have I + S ∆v * ≈ I D + S + (v * c /∆v * )I B . Thus, the Bragg intensity can be estimated from I meas by first subtracting the continuous scattering beneath the Bragg peak and then scaling by ∆v * /v * c (the Lorentz correction), as follows: where I bkg is an estimate for the continuous scattering contribution at h.

Incoherent scattering and absolute intensity
A convenient method for determining the absolute scale factor based on the integral of the total intensity was first described by Krogh-Moe [17] and further elucidated by Norman [18]. The total intensity includes the coherent scattering (Bragg and diffuse) as well as incoherent Compton scattering. The total integrated intensity (per unit cell) is defined as follows: where I and S are the coherent and incoherent scattering per unit cell. Compton scattering is insensitive to molecular structure and can be calculated given the atomic inventory of the unit cell as follows: where the sum runs over all atoms in the unit cell and S n (s) are the atomic incoherent scattering functions [11,12]. Application of Parseval's theorem from harmonic analysis allows us to write the integral of the elastic scattering in terms of independent contributions from each atom, as follows: where f n are the coherent atomic scattering factors [19]. Thus, both the coherent and incoherent contributions to the total intensity can be calculated from the atomic scattering functions, as follows: Alternatively, the total intensity can be found by directly integrating the measured three-dimensional reciprocal space map. Care must be taken however to place Bragg and diffuse intensities on the same (arbitrary) scale. If Bragg peaks are modeled as delta functions (Equation 17), the intensity integral is as follows: As discussed in the previous section, the Bragg intensity and the continuous scattering (proportional to I D + S) are measured on an arbitrary scale. To determine the unknown scale factor, the total measured intensity (Equation 39) can be compared with the calculated theoretical intensity (Equation 38). There are two complications, however. The first is that the Bragg peak at h=0 is unmeasurable. The second is that the limits of integration in Equations 38 and 39 go to infinity, while the measurement ends at some finite maximum scattering vector, s max . To overcome these complications in the theoretical scattering, the h=0 Bragg intensity is subtracted and the integral truncated at s max , as follows: where Z = ∑ n f n (0) is the number of electrons in the unit cell. This is compared with the integral of the measured data, over the same region of reciprocal space: where I B and I C are the measured Bragg and continuous intensities on an arbitrary scale. The scale factor is calculated from the ratio of predicted and measured total scattering (Equations 40 and 41).

Lattice Dynamics Simulation
Lattice dynamics techniques are used to calculate the vibrational frequencies and normal modes (dispersion relations) of atomic or molecular crystals [20,21], often in the context of diffuse scattering [22,23]. Here we consider a class of lattice models where the forces between atoms are modeled as a network of springs, and the atoms are assigned to groups that move collectively as rigid bodies.

Equations of motion in a rigid-body coordinate system
The equations of motion for the network can be expressed using a generalized coordinate system for rigid-body displacement. The instantaneous displacement of a rigid group from its resting position is described by a sixcoordinate vector, w = (t 1 ,t 2 ,t 3 , λ 1 , λ 2 , λ 3 ) , where t and λ are the vectors of translation and libration [6]. For small rotations, the displacement u of an atom in the group is related to w by a linear operator: u = A(r)w, where A depends on the coordinates of the atom at rest, r = (r 1 , r 2 , r 3 ), as follows: In the derivation below, we introduce a compact notation where the cartesian displacements of all n atoms in a unit cell l are represented by a 3n-dimensional column vector u (l) , and similarly the generalized coordinates of all m groups in the unit cell are represented by a 6m-dimensional column vector w (l) . A 3n-by-6m matrix A is constructed so that u (l) = A w (l) . (43) Using this notation, the harmonic part of the total potential energy is where V (l,l ) is the portion of the Hessian matrix (or the force constants matrix) pertaining to the interaction between a unit cell l and l , defined as follows: where V is the total potential energy. The total kinetic energy is where M is the matrix of generalized masses (which includes the moments of inertia), defined as follows: where diag(m) is a square matrix with the atomic masses along the diagonal, I 3 is the 3 × 3 identity matrix, and ⊗ is the Kronecker product. The equations of motion can be derived from the kinetic and potential energies (Equations 44 and 46) using Lagrangian mechanics, with the following result: The generalized mass matrix is positive definite and can be decomposed as M = LL T , where L is a lower triangular matrix (found by Cholesky decomposition).

Normal modes
The Born / von-Karman approach [20][21][22] can be used to solve the equations of motion (Equation 48) by transforming from real space to reciprocal space (k-space, where k is the wavevector). First, displacement-wave solutions are proposed with the following form: where e is a (complex) polarization vector, r l is the origin of unit cell l, ω is the angular frequency, and t is time. Substituting these solutions into the equations of motion, This can be cast in the form of an eigenvalue equation, ω 2 e = D (k) e, where D (k) is known as the dynamical matrix.
Let ω 2 (k, j) and e (k, j) be the j th eigenvalue and eigenvector of D (k) . These satisfy Here we consider finite supercells consisting of N primitive unit cells and periodic boundary conditions so that k is restricted to a discrete set. It can be shown [21] that the general solution is a sum of all allowed waves, or normal modes, as follows: where σ (k, j) are complex, time-dependent amplitudes of the normal modes.

Displacement covariances for thermally-excited vibrations
In thermodynamic equilibrium, the mean-squared amplitudes of the normal coordinates can be found by applying the equipartition theorem (high-temperature limit): Using this result, the covariance matrix can be derived for the generalized coordinates. It can be expressed in the following compact form: where K (k) is the Fourier transform of the force constants matrix (K (k) = LD (k) L T ) and K + (k) is its generalized inverse (in the generalized inverse, normal modes corresponding to rotation and translations of the entire crystal have ω = 0 and are excluded). The covariance of individual atomic displacements can be calculated by projection (Equation 42), as follows: where the index κ refers to the rigid group to which atom j belongs (w (κl) is a 6-element vector). The self-terms of the covariance matrix are the familiar atomic displacement parameters (ADPs) in protein crystallography: The equivalent isotropic B-factor is defined as Because the generalized coordinates of the lattice model w are the same as those used in TLS refinement [6], the T, L, and S matrices are 3 × 3 blocks within the larger covariance matrix for w. For a particular rigid group κ, the matrices are related as follows:

Diffuse scattering and the one-phonon approximation
It is customary in phonon scattering theory to represent the reciprocal space coordinate by the momentum transfer vector q rather than the scattering vector s, with the relationship q = 2πs. We follow this convention below. Since the lattice vibrations are harmonic, the diffuse scattering intensity per unit cell can be calculated in the harmonic approximation: l,l e iq·(r l −r l ) ∑ j, j f j f j e iq·(r j −r j ) T j T j T l j,l j − 1 .
The first summation is over all pairs of unit cells indexed by l and l , and r l is the origin of unit cell l. The second summation is over all atom pairs in unit cell l (indexed by j) and l (indexed by j ). Here, f j is the atomic scattering factor and r j is the average position of the atom relative to the unit cell origin. The factors T j (the Debye-Waller factor) and T l j,l j depend on the self-and cross-terms of the displacement covariance matrix, as follows: and T l j,l j = exp 1 2 q T u j l u T jl + u jl u T j l q , where u jl is the instantaneous displacement of atom j in unit cell l and U j is defined in Equation 56. When the cross-terms of the covariance matrix are small, Equation 61 can be expanded as a Taylor series, and the intensity separated into terms of increasing order (I = I (1) + I (2) + . . . .) where I (1) is the "one-phonon" diffuse scattering, I (2) is the "two-phonon" diffuse scattering, and so on. In the "one-phonon" term, a phonon with wavevector k contributes only to points in reciprocal space that are displaced by k from the Bragg peak locations g h , as follows: This can be written in a more compact form that is computationally convenient by defining the vector-valued "onephonon structure factor", Then, the one-phonon intensity can be written as follows: where K + (k) is the generalized inverse of the force constants matrix (see Equation 54).