International
Tables for
Crystallography
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 3.2, pp. 416-417   | 1 | 2 |
doi: 10.1107/97809553602060000769

Appendix A3.2.1. 

R. E. Marsha* and V. Schomakerb

aThe Beckman Institute–139–74, California Institute of Technology, 1201 East California Blvd, Pasadena, California 91125, USA, and  bDepartment of Chemistry, University of Washington, Seattle, Washington 98195, USA
Correspondence e-mail:  rem@xray.caltech.edu

Consider n atoms at observed vector positions [\sf r] (expressed in Cartesians), n constraints (each adjusted position [{\sf r}_{\rm a}] – `a' for `adjusted' – must be on the plane) and [3n + 3] adjustable parameters (3n [{\sf r}_{\rm a}] components and the 3 components of the vector [\sf a] of reciprocal intercepts of the plane), so that the problem has [n - 3] degrees of freedom. The weight matrices [\sspecialfonts{\bsf P}] may be differently anisotropic for each atom, but there are no interatomic correlations. As before, square brackets, `[[\ldots]]', represent the Gaussian sum over all atoms, usually suppressing the atom indices. We also write [\lambda], not the [\lambda] of Section 3.2.2[link], for the Lagrange multipliers (one for each atom); [\sf m] for the direction cosines of the plane normal; and d for the perpendicular origin-to-plane distance.

As before, [\sspecialfonts{\bsf P}_{k}] is the reciprocal of the atomic error matrix: [\sspecialfonts{\bsf P}_{k} = {\bsf M}_{k}^{-1}] (correspondingly, [{\bi P} \equiv {\bi M}^{-1})] but `[\sspecialfonts{\bsf M}]' is no longer the `[\sspecialfonts{\bsf M}]' of Section 3.2.2.[link] The appropriate least-squares sum is[\sspecialfonts S = [({\sf r} - \hbox{\sf r}_{\rm a})^{T} {\bsf P}({\sf r} - {\sf r}_{\rm a})]]and the augmented sum for applying the method of Lagrange multipliers is[\Phi = S/2 - [\lambda {\sf a}^{T} {\sf r}_{\rm a}].][\Phi] is to be minimized with respect to variations of the adjusted atom positions [{\sf r}_{ka}] and plane reciprocal intercepts [b_{i}], leading to the equations[\sspecialfonts\eqalign{ 0 &= {-\partial \phi\over \partial {\sf r}_{\rm a}^{T}} = {\bsf P}({\sf r} - {\sf r}_{\rm a}) + \lambda {\sf b} \,\,{\rm and}\,\, \cr 0 &= {-\partial \Phi\over \partial {\sf b}^{T}} = [\lambda {\sf r}_{\rm a}],}]subject to the plane conditions [{\sf b}^{T}{\sf r}_{\rm a} = 1], each atom, with [d^{2} = 1/({\sf b}^{T}{\sf b})], [{\sf m} = {\sf b}/\sqrt{{\sf b}^{T}{\sf b}}]. These equations are nonlinear.

A convenient solution runs as follows: first multiply the first equation by [\sspecialfonts{\bsf M}] and solve for the adjusted atom positions in terms of the Lagrange multipliers [\lambda] and the reciprocal intercepts [\sf b] of the plane; then multiply that result by [{\sf b}^{T}] applying the plane conditions, and solve for the [\lambda]'s[\sspecialfonts\eqalign{ {\sf r}_{\rm a} &= {\sf r} + \lambda {\bsf M}{\sf b},\quad {\bsf M} \equiv {\bsf P}^{-1}\cr 1 &= {\sf b}^{T} {\sf r} + \lambda {\sf b}^{T} {\bsf M}{\sf b},\cr \lambda &= (1 - {\sf b}^{T}{\sf r})/({\sf b}^{T} {\bsf M}{\sf b}).}]Next insert these values for the [\lambda]'s and [{\sf r}_{\rm a}]'s into the second equation:[\sspecialfonts 0 = {\partial \Phi\over \partial {\sf b}^{T}} = [\lambda {\sf r}_{\rm a}] = \left[{1 - {\sf b}^{T}{\sf r}\over {\sf b}^{T} {\bsf M}{\sf b}} \left({\sf r} + {1 - {\sf b}^{T}{\sf r}\over {\sf b}^{T}{\bsf M}{\sf b}} {\bsf M}{\sf b}\right)\right].]

This last equation, [F({\sf b}) = 0], is to be solved for [{\sf b}]. It is highly nonlinear: [F({\sf b}) = O({\sf b}^{3}) / O({\sf b}^{4})]. One can proceed to a first approximation by writing [0 = [(1 - {\sf b}^{T}{\sf r}) {\sf r}]]; i.e., [\sspecialfonts[{\bsf r}{\bsf r}]\cdot {\sf b} = [{\bsf r}]], in dyadic notation. [\sspecialfonts[{\bsf M} = {\bsf I}], all atoms; [1 - {\sf b}^{T} {\sf r} = 0] in the multiplier of [\sspecialfonts{\bsf M}{\sf b}/({\sf b}^{T}{\bsf M}{\sf b}).]] A linear equation in [\sf b], this approximation usually works well.6 We have also used the iterative Frazer, Duncan & Collar eigenvalue solution as described by SWMB (1959[link]), which works even when the plane passes exactly through the origin. To continue the solution of the nonlinear equations, by linearizing and iterating, write [F({\sf b}) = 0] in the form[\sspecialfonts\eqalign{ 0 &= [\lambda {\sf r} + \lambda^{2} {\bsf M}{\sf b}] \cr &\approx \left[{\sf r}{\partial \lambda\over \partial {\sf b}} + 2{\bsf M}{\sf b} \lambda {\partial \lambda\over \partial {\sf b}} + \lambda^{2}{\bsf M}\right] \delta {\sf b} + [\lambda ({\sf r} + \lambda {\bsf M}{\sf b})]_{0},}]solve for [\delta{\sf b}], reset [\sf b] to [{\sf b} + \delta {\sf b}], etc., until the desired degree of convergence of [|\delta {\sf b}|/|{\sf b}|] toward zero has been attained. By [\partial \lambda / \partial {\sf b}] is meant the partial derivative of the above expression for [\lambda] with respect to [\sf b], as detailed in the next paragraph.

In the Fortran program DDLELSP (double precision Deming Lagrange, with error estimates, least-squares plane, written to explore this solution) the preceding equation is recast as[\sspecialfonts\eqalign{ {\bsf B} \delta {\sf b} &\equiv \left[\lambda^{2} {\bsf M} + ({\sf r} + 2\lambda {\bsf M}{\sf b}) {\partial \lambda\over \partial {\sf b}}\right] \delta {\sf b}\cr &\equiv {\bsf Y} = - [\lambda ({\sf r} + \lambda {\bsf M}{\sf b})]_{0},}]with[\sspecialfonts\eqalign{ {\partial \lambda\over \partial {\sf b}} &= {\sf r}^{T} / {\sf b}^{T} {\bsf M}{\sf b} - [2(1 - {\sf b}^{T} {\sf r}) / ({\sf b}^{T} {\bsf M}{\sf b})^{2}]\cr &= -({\sf r}^{T} + 2\lambda {\sf b}^{T} {\bsf M}) / {\sf b}^{T} {\bsf M}{\sf b}.}]

The usual goodness of fit, GOF2 in DDLELSP, evaluates to[\sspecialfonts\eqalign{ \hbox{GOF2} &= \left({S_{\min}\over n - 3}\right)^{1/2} = \left({1\over n - 3} [\lambda^{2}{\sf b}^{T} {\bsf M}{\bsf P}{\bsf M}{\sf b}]\right)^{1/2}\cr &= \left({1\over n - 3} [\lambda^{2}{\sf b}^{T}{\bsf M}{\sf b}]\right)^{1/2}\cr &= \left({1\over n - 3} \left[{(1 - {\sf b}^{T} {\sf r})^{2}\over {\sf b}^{T}{\bsf M}{\sf b}}\right]\right)^{1/2}.}]This is only an approximation, because the residuals [1-{\sf b}^{T}{\sf r}] are not the differences between the observations and appropriate linear functions of the parameters, nor are their variances (the [\sspecialfonts{\sf b}^{T} {\bsf M}{\sf b}]'s) independent of the parameters (or, in turn, the errors in the observations).

We ask also about the perpendicular distances, e, of atoms to plane and the mean-square deviation [\overline{(\delta e)^{2}}] to be expected in e.[\eqalign{ e &= (1 - {\sf b}^{T} {\sf r})/\sqrt{{\sf b}^{T} {\sf b}} = d(1 - {\sf b}^{T} {\sf r})\cr \delta e &= -d ({\sf b}^{T} \eta + {\sf r}^{T} \varepsilon) + O (\eta^{2}).}]Here [\eta] and [\varepsilon] are the errors in [\sf r] and [\sf b]. Neglecting `[O(\eta^{2})]' then leads to[\overline{(\delta e)^{2}} = d^{2} ({\sf b}^{T} \overline{\eta \eta^{T}} {\sf b} + 2{\sf b}^{T} \overline{\eta \varepsilon^{T}} {\sf r} + {\sf r}^{T} \overline{\varepsilon \varepsilon^{T}} {\sf r}).]We have [\sspecialfonts\overline{\eta \eta^{T}} = {\bsf M} = {\bsf P}^{-1}], but [\overline{\eta \varepsilon^{T}}] and [\overline{\varepsilon \varepsilon^{T}}] perhaps still need to be evaluated.

References

Schomaker, V., Waser, J., Marsh, R. E. & Bergman, G. (1959). To fit a plane or a line to a set of points by least squares. Acta Cryst. 12, 600–604.








































to end of page
to top of page