research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
APPLIED
CRYSTALLOGRAPHY
ISSN: 1600-5767

Programming new geometry restraints: parallelity of atomic groups

aLawrence Berkeley National Laboratory, One Cyclotron Road, MS64R0121, Berkeley, CA 94720, USA, bDepartment of Bioengineering, University of California Berkeley, Berkeley, CA 94720, USA, cCentre for Integrative Biology, IGBMC, CNRS-INSERM-UdS, 1 rue Laurent Fries, BP 10142, 67404 Illkirch, France, and dDépartement de Physique, Faculté des Sciences et des Technologies, Université de Lorraine, 54506 Vandoeuvre-lès-Nancy, France
*Correspondence e-mail: osobolev@lbl.gov

Edited by D. I. Svergun, European Molecular Biology Laboratory, Hamburg, Germany (Received 24 February 2015; accepted 31 May 2015; online 8 July 2015)

Improvements in structural biology methods, in particular crystallography and cryo-electron microscopy, have created an increased demand for the refinement of atomic models against low-resolution experimental data. One way to compensate for the lack of high-resolution experimental data is to use a priori information about model geometry that can be utilized in refinement in the form of stereochemical restraints or constraints. Here, the definition and calculation of the restraints that can be imposed on planar atomic groups, in particular the angle between such groups, are described. Detailed derivations of the restraint targets and their gradients are provided so that they can be readily implemented in other contexts. Practical implementations of the restraints, and of associated data structures, in the Computational Crystallography Toolbox (cctbx) are presented.

1. Introduction

The restrained refinement of atomic models is performed through minimization of a target function. First introduced in the macromolecular crystallography field by Konnert (1976[Konnert, J. H. (1976). Acta Cryst. A32, 614-617.]) and Jack & Levitt (1978[Jack, A. & Levitt, M. (1978). Acta Cryst. A34, 931-935.]), this target is a weighted sum of several terms including RX describing the model fit to the experimental data:

[R = {w_X}{R_X} + \textstyle\sum\limits_{n = 1}^N {{R_n}}. \eqno (1)]

In the case of refinement of atomic coordinates, the terms Rn describe the model geometry and are known as geometric or stereochemical restraints.

At low resolution there may be many models that may fit the experimental data equally well, but not all of these models may be chemically reasonable. This is why additional information needs to be used to restrict the models to those that conform to expectations about the correct chemical structure. These expectations may originate from various sources, such as small-molecule crystallographic studies, high-quality macromolecular structures or quantum chemistry calculations, forming the basis for libraries of a priori knowledge. At lower resolution, there are fewer experimental data available and therefore a larger number of restraints are required to ensure a high-quality refined model. Therefore, low-resolution structural studies, which have become more and more frequent, require additional restraints (see, for example, Smart et al., 2012[Smart, O. S., Womack, T. O., Flensburg, C., Keller, P., Paciorek, W., Sharff, A., Vonrhein, C. & Bricogne, G. (2012). Acta Cryst. D68, 368-380.]; Headd et al., 2012[Headd, J. J., Echols, N., Afonine, P. V., Grosse-Kunstleve, R. W., Chen, V. B., Moriarty, N. W., Richardson, D. C., Richardson, J. S. & Adams, P. D. (2012). Acta Cryst. D68, 381-390.]; Brown et al., 2015[Brown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). Acta Cryst. D71, 136-153.]) to complement the traditionally used restraints, such as bond lengths, bond angles, dihedral angles, non-bonded interactions, and planarity and chirality restraints.

As an example, nucleic acids have a number of base pairs, some of which are parallel or nearly parallel to each other. This knowledge may serve as a restraint that can be applied to the refinement of such models. Adding a new restraint into the refinement target means developing a corresponding algorithm to calculate its value, as well as an algorithm to calculate the gradient of this restraint function. We use the example of the parallelity restraint, implemented in cctbx (Grosse-Kunstleve et al., 2002[Grosse-Kunstleve, R. W., Sauter, N. K., Moriarty, N. W. & Adams, P. D. (2002). J. Appl. Cryst. 35, 126-136.]), to illustrate the general approach of adding new restraints in cctbx and PHENIX (Adams et al., 2010[Adams et al. (2010). Acta Cryst. D66, 213-221.]) and also to reinforce the general principles of crystallographic refinement programs. This article describes the algorithms in a ready-to-program manner, with sufficient mathematical and computational details to serve a didactic purpose as well. Practical examples of the application of this new restraint function to model refinement will be presented elsewhere.

Each individual target in equation (1)[link] is the sum of a large number of identical terms, each depending on a small number of parameters. For example, a diffraction target usually depends on a structure factor amplitude or the parameters of its phase distribution, a density target depends on the density values in a few grid points around a particular atom (see, for example, Afonine & Urzhumtsev, 2004[Afonine, P. V. & Urzhumtsev, A. (2004). Acta Cryst. A60, 19-32.], and references therein), a geometric restraint depends on the coordinates of a few atoms involved in particular restraint, etc. Such a composition of the target function (1)[link] is very important for several reasons. First, it allows the manipulation of each kind of restraint, e.g. selection of the atoms involved or an update of the function. Second, each restraint is relatively simple and easy to program. Third, each restraint can be weighted individually. Finally, the gradient of the crystallographic target used for the minimization can be calculated as a sum of gradients from the individual terms.

Several general points are important to note when designing refinement programs (Lunin & Urzhumtsev, 1985[Lunin, V. Yu. & Urzhumtsev, A. G. (1985). Acta Cryst. A41, 327-333.]):

(i) Refinement targets may be functions of the atomic parameters and values of electron density (or structure factors), in combination with the restraints and/or constraints imposed on the atomic model. Each component of the overall target is a relatively simple function when it is expressed through the parameters of the appropriate type.

(ii) These parameters are calculated from each other as a chain of consecutive transitions.

(iii) The gradient calculations are performed step by step using the chain rule, inverting the chain of steps used to calculate the target (Baur & Strassen, 1983[Baur, W. & Strassen, V. (1983). Theor. Comput. Sci. 22, 317-330.]; Kim et al., 1984[Kim, K. M., Nesterov, Yu. E. & Cherkassky, B. V. (1984). Dokl. Acad. Nauk. SSSR, 275, 1306-1309.]). This procedure results in exact gradient values that require approximately the same amount of time to calculate as a single value of the target, and this is independent of the number of parameters to be refined.

(iv) The gradient calculation does not require an analytical function; it may be any algorithm that provides a value for the corresponding target. At each refinement step, the target or gradient calculation is a local operation. It takes the numerical values obtained at the previous step, calculates the result of the step and passes the numerical values to the next step. Such a scheme allows one to avoid complicated mathematical equations connecting the final target or gradient value with the initial model parameters.

In order to introduce a parallelity restraint into refinement, the considerations stated above require that we propose an algorithm to calculate a numerical value describing by how much two atomic groups that form pseudo-planes (approximately belong to two planes) are not parallel to each other. This algorithm should be numerically robust enough to handle a broad range of model parameters, from good (towards the end of refinement) to distorted (at the beginning of refinement).

2. Possible parallelity restraints and their gradient

2.1. Definition of the `best plane'

Let us consider two pseudo-planar atomic groups (groups that should be planar but may not currently be so), m = 1, 2, containing atoms defined by their Cartesian coordinates, rk,m = (xk,m, yk,m, zk,m), k = 1, …, Km, in the same coordinate system. For each such group we need to define the corresponding `best plane'. A plane is defined by a vector Nm normal to it and by a point Cm that belongs to it. For a given atomic group, the plane can be characterized by the sum of the squared deviations of the atoms from it (Urzhumtsev, 1991[Urzhumtsev, A. G. (1991). Acta Cryst. A47, 723-727.]; Blanc & Paciorek, 2001[Blanc, E. & Paciorek, W. (2001). J. Appl. Cryst. 34, 480-483.]; Brown et al., 2015[Brown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). Acta Cryst. D71, 136-153.]), and these deviations may be weighted if desired. The `best plane' is the one that corresponds to the minimum of this value (Nm and Cm are variables),

[\eqalignno{ \sum\limits_{k = 1}^{K_m} w_{k,m} \left [ \left ( {\bf r}_{k,m} - {\bf C}_m \right ) \cdot {\bf n}_m \right ]^2 = & \, \sum\limits_{k = 1}^{K_m} w_{k,m} {{ \left [ \left ( {\bf r}_{k,m} - {\bf C}_{m} \right ) \cdot {\bf N}_{m} \right ]^2} \over {{\bf N}_m \cdot {\bf N}_m }} \cr & \, \to \min\limits_{{\bf N}_m , {\bf C}_m }, &(2)}]

where wk,m is a weight (e.g. a unit value or an atomic mass) for the contribution of an atom with coordinates rk,m to the target. Here and below, a centred dot `·' represents the scalar (dot) product of two vectors, and [{\bf n}_m = {\bf N}_m ({\bf N}_m \cdot {\bf N}_m)^{-1/2}] is the normalized vector Nm.

Two facts are useful (see, for example, Urzhumtsev, 1991[Urzhumtsev, A. G. (1991). Acta Cryst. A47, 723-727.]). First, the minimum of equation (2)[link] is reached when Cm is the weighted geometric center of the group

[{\bf C}_m = \left ( \textstyle\sum\limits_{k = 1}^{K_m} w_{k,m} \right )^{-1} \textstyle\sum\limits_{k = 1}^{K_m} w_{k,m} {\bf r}_{k,m} . \eqno(3)]

Second, if for a given atomic group we calculate the minimal eigenvalue of its inertia matrix (see §3.2[link]), the corresponding eigenvector is normal to the best plane defined by condition (2)[link] and therefore it can be taken as Nm.

The angle between two atomic groups can be described by the angle θ between two vectors normal to the best planes for the corresponding groups. This angle varies between the limits 0 ≤ θπ/2. This angle is expected to be zero if we require the planes to be parallel, which is the most frequent restraint of this kind. We therefore need to design a function which reaches its minimum when θ = 0. Note that there are multiple functions which could be constructed such that θ0 = 0 is a minimum. Also, the particular restraint function needs to describe the stereochemical knowledge near the point of interest (minimum) and have a desired behavior far from it. Additionally, the computation of the restraint and its derivatives must have no numerical irregularities. Below, we analyze some possible restraints that describe the angle between two planar groups, noting that different restraints can be used in different circumstances.

2.2. Analysis of various possible parallelity targets

The computationally simplest target to restrain the angle between normal vectors is

[\eqalignno{ R_{\rm parallelity} = & \, w \left [ 1 - \cos \left ( \theta - \theta _0 \right) \right ] \cr = & \, w \left [ 1 - \left ( \cos \theta \cos {\theta _0} + \sin \theta \sin {\theta _0} \right ) \right ] , &(4)}]

where θ = θ(n1, n2) is the angle between two unit vectors n1 and n2, θ0 is the target angle, and w is the individual weight factor for a particular pair of planes. By default, phenix.refine (Afonine et al., 2012[Afonine, P. V., Grosse-Kunstleve, R. W., Echols, N., Headd, J. J., Moriarty, N. W., Mustyakimov, M., Terwilliger, T. C., Urzhumtsev, A., Zwart, P. H. & Adams, P. D. (2012). Acta Cryst. D68, 352-367.]) uses the form w = 1/σ2, where σ is the root mean-square deviation from the ideal value obtained by some means. Introducing the angle θ0 in equation (4)[link] will force the planes to form the prescribed angle between them. For example, as suggested by Brown et al. (2015[Brown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). Acta Cryst. D71, 136-153.]), a restraint with θ0 = π/2 is important when working with aromatic groups. For simplicity, we continue to call this restraint `parallelity', even though this is not completely appropriate. We note that cosθ = n1·n2 and sinθ = ||n1 × n2||.

For small deviations from θ0, restraint (4)[link] behaves as the quadratic function [{1 \over 2}]w(θθ0)2. However, using the quadratic function of the angle itself [see in particular equation (14) of Brown et al. (2015[Brown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). Acta Cryst. D71, 136-153.])] requires its derivatives, i.e. the derivative of arccos(n1·n2), which is undefined for θ = θ0. Also, its calculation is computationally unstable for θθ0. For this reason, we have excluded this kind of target from our consideration.

Also, a seemingly plausible restraint function that one may wish to avoid is

[R_{\rm parallelity} = w \left ( \cos \theta - \cos {\theta _0} \right )^2 . \eqno(5)]

This function is easy to calculate but it has very different behavior across different ranges of the parameter θ0. For example, it is similar to a second-order polynomial near θ0 = π/2 and to a fourth-order polynomial near θ0 = 0. Such heterogeneity is generally disadvantageous for minimization.

Another known target for the parallelity restraint is

[R_{\rm parallelity} = \textstyle\sum\limits_{k = 1}^{K_1} w_{k,1} \left ( {\bf r}_{k,1} \cdot {\bf n}_1 - d_1 \right )^2 + \sum\limits_{k = 1}^{K_2} w_{k,2} \left ( {\bf r}_{k,2} \cdot {\bf n}_1 - d_2 \right )^2 , \eqno(6)]

[equation (13) of Brown et al. (2015[Brown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). Acta Cryst. D71, 136-153.]); the notation of the present article is used]. It imposes a rigid condition on the parallelity of two planes, n2 = n1, simultaneously with planarity of the corresponding groups. This constraint-type condition may be useful, for example, to simplify the parameterization of the model, but in practice nucleic acid bases may bend from one base pair to another, or they may be twisted inside a given base pair. As a consequence, it is more convenient to have planarity and parallelity restraints defined separately, allowing for more flexibility if desired, in particular allowing the interplanar angle to vary from one planar group to another.

When the function is far from the minimum (e.g. when the current structure is distorted with respect to the restraints), one may wish that the target be relatively insensitive to the variation in the angle θ. This can be achieved by using various kinds of `top-out' functions (Dennis & Welsch, 1987[Dennis, J. E. Jr & Welsch, R. E. (1987). Commun. Stat. Simul. Comput. 7, 345-359.]; Murshudov et al., 2011[Murshudov, G. N., Skubák, P., Lebedev, A. A., Pannu, N. S., Steiner, R. A., Nicholls, R. A., Winn, M. D., Long, F. & Vagin, A. A. (2011). Acta Cryst. D67, 355-367.]; Smart et al., 2012[Smart, O. S., Womack, T. O., Flensburg, C., Keller, P., Paciorek, W., Sharff, A., Vonrhein, C. & Bricogne, G. (2012). Acta Cryst. D68, 368-380.]; Headd et al., 2012[Headd, J. J., Echols, N., Afonine, P. V., Grosse-Kunstleve, R. W., Chen, V. B., Moriarty, N. W., Richardson, D. C., Richardson, J. S. & Adams, P. D. (2012). Acta Cryst. D68, 381-390.]). A direct application of this idea results in

[R_{\rm parallelity} = w {\Omega ^2} \left \{ 1 - \exp \left [ {{{\cos \left ( {\theta - {\theta _0}} \right ) - 1} \over {\Omega ^2}}} \right ] \right \} , \eqno(7)]

where Ω is a constant (Fig. 1[link]). Using (θθ0)2 or (cosθ − cosθ0)2 in equation (7)[link] instead of cos(θθ0) − 1 is not desirable for the reasons discussed above. In fact, we can use directly a sigmoid form of the cosine function and introduce the `top-out type' targets in a simpler form:

[R_{\rm parallelity} = w \left [ 1 - \cos 2 \left ( \theta - {\theta _0} \right ) \right ] = 2 w \left [ 1 - {\cos }^2 \left ( \theta - {\theta _0} \right ) \right ] \eqno(8)]

or

[R_{\rm parallelity} = \cases { w \left [ 1 - \cos n \left ( \theta - {\theta _0} \right) \right ] & if $ \left | \theta - {\theta _0} \right | \le \pi/n$, \cr 2w & if $ \left | \theta - {\theta _0} \right | \,\,\gt\,\, \pi/n $,} \eqno (9)]

where n > 2. These functions (Fig. 1[link]) have harmonic behavior when θθ0, while for large values of the argument they are constant [equation (9)[link]] or nearly constant [equation (8)[link]].

[Figure 1]
Figure 1
Examples of various parallelity targets, θ0 = 0. Heavy line with no markers: f(θ) = 1 − cosθ, equation (4)[link]. Line with filled triangles: f(θ) = 1 − cos2θ, equation (8)[link]. Line with open triangles: f(θ) = 1 − cos4θ for θπ/4, equation (9)[link]. Lines with black (Ω = 1) and open (Ω = 2) circles: f(θ) = Ω2{1 − exp[(cosθ − 1)/Ω2]}, equation (7)[link]. Line with filled squares: f(θ) = (1 − cosθ)4, equation (10)[link]. Line with asterisks: f(θ), slack function, equation (11)[link].

A further useful modification of target function (4)[link] can be achieved by introducing the `slack' parameter (square well-like potential), which scores small deviations from the target value equally. This can be achieved in a `soft' way (Fig. 1[link]) as

[R_{\rm parallelity} = w \left [ 1 - \cos \left ( \theta - {\theta _0} \right ) \right ]^n, \quad n \ge 2 , \eqno (10)]

or following explicitly the slack definition (Fig. 1[link]) as

[R_{\rm parallelity} = \left\{\matrix { 0 \hfill & {\rm if} \,\left | \theta - {\theta _0} \right | \le {\rm slack}, \hfill \cr w [ 1 - \cos ( \theta - {\theta _0} - {\rm slack} ) ] & {\rm if} \, {\rm slack} \, \lt \, \theta - {\theta _0}, \hfill \cr w [ 1 - \cos ( \theta - {\theta _0} + {\rm slack} ) ] & {\rm if }\,\theta - {\theta _0} \,\lt\, - {\rm slack}. }\right. \eqno (11)]

2.3. Practical calculation of the parallelity targets and their gradient

All the restraint functions mentioned above are similar in the sense that they require calculation of the normal vectors to the optimal planes of the atomic groups; then their dot product and finally a vector product of these vectors are calculated to obtain the target value.

A number of other restraints, such as group planarity or the distance between parallel groups, are naturally expressed through the same types of parameters used for the parallelity restraint (Appendix A[link]). These parameters are the center of the atomic group, and the eigenvalues and eigenvectors of its inertia matrix (generally speaking, all three eigenvalues and all three eigenvectors may be required); a particular target may use only some of the parameters mentioned. Thus, the calculation of these parameters from the initial set of atomic coordinates is common to all such restraints on planar groups (Fig. 2[link]). For this reason, in §3[link] we provide details of these common steps.

[Figure 2]
Figure 2
The overall calculation scheme. (a) The steps used to calculate parameters common to all targets. The letters af from top to bottom refer to the steps for calculating the target described in §3[link]. The same steps in the direction from bottom to top refer to the steps for calculating the derivatives described in §4[link]. The steps are applied independently to each atomic group required to be planar. (b) The particular parallelity targets that can be constructed on the basis of the calculated parameters of an atomic group.

Assuming we have a function that we can calculate numerically, the principle of efficient gradient calculation (Baur & Strassen, 1983[Baur, W. & Strassen, V. (1983). Theor. Comput. Sci. 22, 317-330.]; Kim et al., 1984[Kim, K. M., Nesterov, Yu. E. & Cherkassky, B. V. (1984). Dokl. Acad. Nauk. SSSR, 275, 1306-1309.]) suggests inverting the steps of the function's calculation and obtaining the gradient step by step using the chain rule. This would result in the exact value of the gradient with respect to the original parameters, and this calculation will take no more than four times the CPU time required to calculate a single value of the function. This factor of four is invariant with respect to the number of parameters and the function type. In crystallographic practice, this factor is typically closer to one than to four (Lunin & Urzhumtsev, 1985[Lunin, V. Yu. & Urzhumtsev, A. G. (1985). Acta Cryst. A41, 327-333.]). §4[link] illustrates this principle by calculating the gradient of the parallelity restraints. It is easy to see that the calculation of the gradient of other restraints for pseudo-plane groups differs only in the first step.

3. Planarity targets calculation

This section describes the calculation of the target value for the suggested restraint. Fig. 2[link] shows the overall scheme, where the letters a to f refer to the corresponding steps described below.

3.1. Centered coordinate system (Step a)

We start by calculating the atomic positions in the coordinate system with the origin shifted into the geometric center Cm of the group, equation (3)[link]. This operation is similar for all pseudo-plane groups, and in what follows we omit the group index m to simplify the formulae. The new parameters are the Cartesian coordinates (Xk, Yk, Zk) of the atomic positions:

[{\bf q}_k = {\bf r}_k - {\left ( \textstyle\sum\limits_{k = 1}^K w_k \right )^{-1} \textstyle\sum\limits_{k = 1}^K w_k {\bf r}_k}, \quad k = 1, \ldots, K , \eqno (12)]

and the Cartesian coordinates (Cx, Cy, Cz) of the weighted geometric center of the group.

3.2. Inertia matrix (Step b)

For a given atomic group, its inertia matrix is

[\eqalignno{ {\bf S} = & \, \left ( { \matrix { {\sum\limits_{k = 1}^K {w_k X_k X_k} } & {\sum\limits_{k = 1}^K {w_k X_k Y_k} } & {\sum\limits_{k = 1}^K {w_k X_k Z_k} } \cr {\sum\limits_{k = 1}^K {w_k Y_k X_k} } & {\sum\limits_{k = 1}^K {w_k Y_k Y_k} } & {\sum\limits_{k = 1}^K {w_k Y_k Z_k} } \cr {\sum\limits_{k = 1}^K {w_k Z_k X_k} } & {\sum\limits_{k = 1}^K {w_k Z_k Y_k} } & {\sum\limits_{k = 1}^K {w_k Z_k Z_k} } \cr } } \right) \cr = & \, \left ( { \matrix { S_{XX} & S_{XY} & S_{XZ} \cr S_{YX} & S_{YY} & S_{YZ} \cr S_{ZX} & S_{ZY} & S_{ZZ} \cr } } \right ) . & (13)}]

This matrix is symmetric, SXY = SYX, SYZ = SZY, SZX = SXZ, but for convenience in subsequent calculations we refer to all nine of its elements. Matrix (13)[link] is positive semidefinite, having always three non-negative eigenvalues, 0 ≤ λ+λλ0, with the corresponding eigenvectors orthogonal to each other. The result of this step contains nine elements of matrix (13)[link] and three coordinates of the center [equation (12)[link]] taken from the previous step.

It is important to remember that at each step we calculate the values of the new parameters and do not keep their analytical expression.

3.3. Characteristic equation (Step c)

The eigenvalues 0 ≤ λ+λλ0 of the inertia matrix (13)[link] are the roots of the cubic characteristic equation

[\lambda ^3 + a_S \lambda ^2 + b_S \lambda + c_S = 0 , \eqno (14)]

with the coefficients

[a_S = - S_{XX} - S_{YY} - S_{ZZ} , \eqno (15a)]

[\eqalignno{ b_S = & \, \det \left ( {\matrix{ S_{XX} & S_{XY} \cr S_{YX} & S_{YY} \cr } } \right ) + \det \left ( {\matrix{ S_{YY} & S_{YZ} \cr S_{ZY} & S_{ZZ} \cr } } \right ) \cr & \, + \det \left ( {\matrix{ S_{ZZ} & S_{ZX} \cr S_{XZ} & S_{XX} \cr } } \right ) , & (15b)}]

[c_S = - \det {\bf S} . \eqno (15c)]

Further steps of the procedure use both the coefficients of the cubic equation (14)[link] and the elements of the inertia matrix (13)[link]. To avoid confusion when calculating the gradients [see, for example, equation (18)[link] below], it is more convenient to introduce new variables formally:

[\eqalign{ \hat S_{XX} = & \, S_{XX}, \quad \hat S_{XY} = S_{XY}, \quad \hat S_{XZ} = S_{XZ}, \cr \hat S_{YX} = & \, S_{YX}, \quad \hat S_{YY} = S_{YY}, \quad \hat S_{YZ} = S_{YZ}, \cr \hat S_{ZX} = & \, S_{ZX}, \quad \hat S_{ZY} = S_{ZY}, \quad \hat S_{ZZ} = S_{ZZ}. \cr} \eqno (16)]

Thus, the new set of parameters contains 15 values: three coordinates (Cx, Cy, Cz) of the center of the atomic group that are simply taken from the previous step, three equation coefficients aS, bS, cS, and nine elements of the inertia matrix [\hat S_{XX},\hat S_{XY},...,\hat S_{ZZ}].

3.4. Eigenvalues of the inertia matrix (Step d)

Different algorithms may be applied to calculate the eigenvalues and choose the minimum of them; we use the one described by Urzhumtsev (1991[Urzhumtsev, A. G. (1991). Acta Cryst. A47, 723-727.]), as it gives analytically the minimal real eigenvalue that is required for the planarity and parallelity restraints. Appendix B[link] gives an extended description in the notation of the present paper. The coordinates of the center of the atomic group [equation (12)[link]] and the matrix elements [equation (16)[link]] are not used at this step and pass to the next steps without change.

3.5. Eigenvectors of the inertia matrix (Step e)

Different algorithms may be applied to calculate the eigenvectors N of the inertia matrix that are solutions of the matrix equation

[\left ( {\hat {\bf S}} - \lambda {\bf I} \right ) {\bf N} = 0 , \eqno(17)]

with the eigenvalues λ obtained in the previous step. Appendix B[link] gives technical details of the method used in our program. For an atomic group that is approximately planar, λ+ < λ, the solution of equation (17)[link] with λ = λ+ defines unambiguously a single direction N+ of the normal to the best plane.

Note that other restraints on the plane groups may require knowledge of other eigenvalues and/or eigenvectors (Appendix A[link]). Therefore, in the general case, the result of this step is a set of parameters fully describing a pseudo-planar atomic group, namely the three coordinates (Cx, Cy, Cz) of its center, all three eigenvalues 0 ≤ λ+λλ0 and all three eigenvectors N+, N, N0 of the inertia matrix.

3.6. Parallelity targets (Step f)

To calculate a planarity restraint for two pseudo-planar atomic groups, each described by the Cartesian coordinates of the vector N1 = N1+ and N2 = N2+ normal to their respective best planes, we may need to follow several intermediate steps (Appendix D[link]). We start (Step f1) from normalization of the normal vectors, {N1, N2} [\to] {n1, n2}. This normalization includes choosing the correct direction of the normalized vectors so that the value of the angle θ between the normal vectors is not larger than π/2. Once the normalized vectors are known, we calculate (Step f2) the cosine and sine of the angle θ between them, which we denote Cθ and Sθ, respectively. Finally (Step f3), we calculate the value of the chosen target.

4. Gradient calculation

4.1. Gradient with respect to the plane parameters

We start by inverting the formulae to calculate the chosen target (Step f3 in Fig. 2[link]) and obtain [{{\partial {R_{\rm parallelity}}} / {\partial {C_\theta }}}] and [{{\partial {R_{\rm parallelity}}} / {\partial {S_\theta }}}] (Appendix D[link]).

Inverting the previous step (Step f2), we take as input the values of [{{\partial {R_{\rm parallelity}}} / {\partial {C_\theta }}}] and [{{\partial {R_{\rm parallelity}}} / {\partial {S_\theta }}}] (and not their analytical expressions), whatever the target Rparallelity is. As a result, we obtain ([\partial R_{\rm parallelity} / \partial n_{1X}], [\partial R_{\rm parallelity} / \partial n_{1Y}], [\partial R_{\rm parallelity} / \partial n_{1Z}]) and ([\partial R_{\rm parallelity} / \partial n_{2X}], [{{\partial {R_{\rm parallelity}}} / {\partial {n_{2Y}}}}], [\partial R_{\rm parallelity} / \partial n_{2Z}]). Finally, inverting the normalization (Step f1) we obtain ([\partial R_{\rm parallelity} / \partial N_{1X}], [\partial R_{\rm parallelity} / \partial N_{1Y}], [\partial R_{\rm parallelity} / \partial N_{1Z}]) and ([\partial R_{\rm parallelity}/ \partial N_{2X}], [\partial R_{\rm parallelity} / \partial N_{2Y}], [\partial R_{\rm parallelity} / \partial N_{2Z}]).

Working with other restraints on plane groups, we can apply the same kind of procedure as Steps f3 to f1 above. For all these targets, we arrive at a set of partial derivatives of a target with respect to the following:

(i) The three eigenvalues for each pseudo-planar group.

(ii) The coordinates of the three eigenvectors for each pseudo-planar group.

(iii) The coordinates of the center of each pseudo-planar group.

For the parallelity restraints discussed above, all these derivatives are equal to zero except those for the coordinates of N1+ and N2+. Further steps of the gradient calculation use these partial derivatives as input data whatever the target is. These steps are common to all restraints and all atomic groups. Therefore, starting with the next section we omit the index used to indicate an atomic group and the index of a particular plane-group restraint for which all the partial derivatives above have been calculated. Partial derivatives with respect to the coordinates (Cx, Cy, Cz) will not change until the last step.

4.2. Gradient with respect to the coefficients of the characteristic equation (Steps e and d)

As mentioned in §§3.4[link] and 3.5[link], different algorithms exist for calculating the eigenvalues and eigenvectors of the inertia matrix from the coefficients of the characteristic equation (14)[link]. According to the main scheme, we invert these algorithms to calculate the corresponding derivatives. Appendix C[link] provides these for the particular method of calculation chosen in our program. The output of this step contains the partial derivatives of the target R with respect to the three coefficients aS, bS, cS of equation (14)[link], to the nine elements [\hat S_{XX}, \hat S_{XY}], …, [\hat S_{ZZ}] of equation (16)[link], and to the three coordinates (Cx, Cy, Cz) of the center of the atomic group that are taken, with no change, from the previous step.

4.3. Gradient with respect to the elements of the inertia matrix (Step c)

Using the partial derivatives of aS, bS, cS [equation (50)[link] in Appendix C[link]] and of the [{\hat {\bf S}}] elements [equation (43)[link] in Appendix C[link]] with respect to the elements of the matrix S [equation (13)[link]], we obtain

[\eqalign{{{\partial R} \over {\partial {S_{XX}}}} & = - {{\partial R} \over {\partial {a_S}}} + \left({{S_{YY}} + {S_{ZZ}}} \right){{\partial R} \over {\partial {b_S}}} + \left({{S_{YZ}}{S_{ZY}} - {S_{YY}}{S_{ZZ}}} \right){{\partial R} \over {\partial {c_S}}} \cr & \quad + {{\partial R} \over {\partial {{\hat S}_{XX}}}}, \cr {{\partial R} \over {\partial {S_{XY}}}} & = - {S_{YX}}{{\partial R} \over {\partial {b_S}}} + \left({{S_{YX}}{S_{ZZ}} - {S_{ZX}}{S_{YZ}}} \right){{\partial R} \over {\partial {c_S}}} + {{\partial R} \over {\partial {{\hat S}_{XY}}}}, \cr {{\partial R} \over {\partial {S_{XZ}}}} & = - {S_{ZX}}{{\partial R} \over {\partial {b_S}}} + \left({{S_{ZX}}{S_{YY}} - {S_{YX}}{S_{ZY}}} \right){{\partial R} \over {\partial {c_S}}} + {{\partial R} \over {\partial {{\hat S}_{XZ}}}}, } \eqno(18)]

and similar expressions for the other six derivatives easily obtained from equation (18)[link] by the corresponding cyclic substitution [X \to Y \to Z \to X] of the indices, e.g. to obtain [{{\partial R} / {\partial {S_{YZ}}}}] we make the following substitutions in the second equation: [XY \to YZ], [YX \to ZY], [ZZ \to XX], [ZX \to XY], [YZ \to ZX].

4.4. Gradient with respect to the coordinates in the centered system (Step b)

We obtain the derivatives with respect to the atomic coordinates, k = 1, …, K,

[\eqalign {{{\partial R} \over {\partial {X_K}}} = & \, 2{X_K}{{\partial R} \over {\partial {S_{XX}}}} + {Y_K}\left({{{\partial R} \over {\partial {S_{XY}}}} + {{\partial R} \over {\partial {S_{YX}}}}} \right) + {Z_K}\left({{{\partial R} \over {\partial {S_{XZ}}}} + {{\partial R} \over {\partial {S_{ZX}}}}} \right), \cr {{\partial R} \over {\partial {Y_K}}} = & \, 2{Y_K}{{\partial R} \over {\partial {S_{YY}}}} + {Z_K}\left({{{\partial R} \over {\partial {S_{YZ}}}} + {{\partial R} \over {\partial {S_{ZY}}}}} \right) + {X_K}\left({{{\partial R} \over {\partial {S_{YX}}}} + {{\partial R} \over {\partial {S_{XY}}}}} \right), \cr {{\partial R} \over {\partial {Z_K}}} = & \, 2{Z_K}{{\partial R} \over {\partial {S_{ZZ}}}} + {X_K}\left({{{\partial R} \over {\partial {S_{ZX}}}} + {{\partial R} \over {\partial {S_{XZ}}}}} \right) + {Y_K}\left({{{\partial R} \over {\partial {S_{ZY}}}} + {{\partial R} \over {\partial {S_{YZ}}}}} \right). } \eqno(19)]

4.5. Gradient of restraint with respect to the original coordinates (Step a)

Using the partial derivatives of equations (3)[link] and (12)[link] and the notation Wj = [w_j {({\sum _{k = 1}^K {w_k} })^{-1}}],

[\eqalign{ {{\partial {C_x}} \over {\partial {x_j}}} = & \, {{\partial {C_y}} \over {\partial {y_j}}} = {{\partial {C_z}} \over {\partial {z_j}}} = {W_j}, \cr {{\partial {X_j}} \over {\partial {x_j}}} = & \, {{\partial {Y_j}} \over {\partial {y_j}}} = {{\partial {Z_j}} \over {\partial {z_j}}} = 1 - {W_j}, \cr {{\partial {X_k}} \over {\partial {x_j}}} = & \, {{\partial {Y_k}} \over {\partial {y_j}}} = {{\partial {Z_k}} \over {\partial {z_j}}} = - {W_j} \quad {\rm if } \quad j \ne k, } \eqno (20)]

and knowing that all cross-term derivatives are equal to 0, we obtain the required gradient with respect to the original atomic coordinates

[\eqalign{ {{\partial R} \over {\partial x_j}} = & \, {{\partial C_x} \over {\partial x_j}}{{\partial R} \over {\partial C_x}} + \sum\limits_{k = 1}^K {{{\partial X_k} \over {\partial x_j}}{{\partial R} \over {\partial X_k}}} = {{\partial R} \over {\partial X_j}} \cr & \, + {W_j} \left ({{{\partial R} \over {\partial C_x}} - \sum\limits_{k = 1}^K {{{\partial R} \over {\partial X_k}}} } \right ), \cr {{\partial R} \over {\partial y_j}} = & \, {{\partial R} \over {\partial Y_j}} + {W_j} \left ({{{\partial R} \over {\partial C_y}} - \sum\limits_{k = 1}^K {{{\partial R} \over {\partial Y_k}}} } \right), \cr {{\partial R} \over {\partial z_j}} = & \, {{\partial R} \over {\partial Z_j}} + {W_j} \left ({{{\partial R} \over {\partial C_z}} - \sum\limits_{k = 1}^K {{{\partial R} \over {\partial Z_k}}} } \right). } \eqno(21)]

5. Implementation of parallelity restraints in cctbx

5.1. Computational Crystallography Toolbox

The Computational Crystallography Toolbox (cctbx) is an open-source library with numerous functionalities necessary to construct crystallography-oriented programs. In fact, it finds applications beyond crystallography, such as small-angle scattering (sastbx; Liu et al., 2012[Liu, H., Hexemer, A. & Zwart, P. H. (2012). J. Appl. Cryst. 45, 587-593.]) or cryo-EM (Afonine et al., 2013[Afonine, P. V., Headd, J. J., Terwilliger, T. C. & Adams, P. D. (2013). Comput. Crystallogr. Newsl. 4, 43-44.]). cctbx contains tools to handle various objects such as crystal symmetry, diffraction data, maps (crystallographic or cryo-EM) and atomic models, and provides an array of building blocks necessary to construct complex procedures such as model refinement (Afonine et al., 2012[Afonine, P. V., Grosse-Kunstleve, R. W., Echols, N., Headd, J. J., Moriarty, N. W., Mustyakimov, M., Terwilliger, T. C., Urzhumtsev, A., Zwart, P. H. & Adams, P. D. (2012). Acta Cryst. D68, 352-367.]) or diffraction data analysis (Sauter et al., 2013[Sauter, N. K., Hattne, J., Grosse-Kunstleve, R. W. & Echols, N. (2013). Acta Cryst. D69, 1274-1282.]; Parkhurst et al., 2014[Parkhurst, J. M., Brewster, A. S., Fuentes-Montero, L., Waterman, D. G., Hattne, J., Ashton, A. W., Echols, N., Evans, G., Sauter, N. K. & Winter, G. (2014). J. Appl. Cryst. 47, 1459-1465.]). In particular, to support macromolecular refinement algorithms cctbx includes tools to handle various geometry restraints (Grosse-Kunstleve & Adams, 2004[Grosse-Kunstleve, R. W. & Adams, P. D. (2004). IUCr Comput. Commission. Newsl. 4, 19-36.]). The main set of geometry restraints includes those on covalent bonds and angles, chirality, dihedral (torsion) angles, non-bonded interactions, planar groups, and now parallelity. The overall restraint target, the sum in equation (1)[link], is the sum of contributions from all restraint types

[\eqalignno{ \sum\limits_{n = 1}^N w_n R_n = & \, \textstyle\sum\limits_{i = 1}^{N _{\rm bonds}} w_i R_{{\rm bonds},i} \, + \sum\limits_{i = 1}^{N _{\rm angles}} w_i R_{{\rm angles},i} \cr & \, +\cdots + \textstyle\sum\limits_{i = 1}^{N _{\rm parallelities}} w_i R_{{\rm parallelities},i} . & (22)}]

If additional information about the molecule or its environment is available and deemed necessary to use, additional terms can be added to equation (22)[link] (Echols et al., 2010[Echols, N., Headd, J. J., Afonine, P. & Adams, P. D. (2010). Comput. Crystallogr. Newsl. 1, 12-17.]; Headd et al., 2012[Headd, J. J., Echols, N., Afonine, P. V., Grosse-Kunstleve, R. W., Chen, V. B., Moriarty, N. W., Richardson, D. C., Richardson, J. S. & Adams, P. D. (2012). Acta Cryst. D68, 381-390.], 2014[Headd, J. J., Echols, N., Afonine, P. V., Moriarty, N. W., Gildea. R. J. & Adams, P. D. (2014). Acta Cryst. D70, 1346-1356.]). Such additional information may be the regular conformations of a protein's secondary structure, non-crystallographic symmetry or similarity to a previously solved structure. Technically, adding additional restraints means adding additional terms to equation (22)[link], which in turn means that an algorithm for calculating the additional restraint target and its gradient with respect to atomic coordinates, and the data structures for holding information about the atoms involved, are required.

The implementation of the parallelity restraint in cctbx is very similar to other restraints (angle, dihedral angle, planarity etc.) and therefore may serve as an example of the common framework for restraints.

5.2. General approach to defining stereochemical restraints in cctbx

The source code for commonly used restraints is located at http://sourceforge.net/p/cctbx/code/HEAD/tree/trunk/cctbx/geometry_restraints in the files bond.h, bond_bpl.cpp, angle.h, angle_bpl.cpp, dihedral.h, dihedral_bpl.cpp, chirality.h, chirality_bpl.cpp, nonbonded.h, nonbonded_bpl.cpp, planarity.h, planarity_bpl.cpp, parallelity.h and parallelity_bpl.cpp.

Files with the name <type of restraint>_bpl.cpp contain boost.python (Abrahams & Grosse-Kunstleve, 2003[Abrahams, D. & Grosse-Kunstleve, R. W. (2003). C/C++ Users J. 21, 29-36.]) descriptions of classes and functions for the restraint, and the bpl in the file names stands for `boost.python library'. Generally, these files contain boost.python wrappers for classes describing the restraint and the associated data structure (termed proxy; see below) holding the information necessary to calculate the restraint. Files with the name <type of restraint>.h contain the actual implementation of classes and functions for the particular type of restraint.

5.3. Proxy class

The calculation of a particular restraint requires information about the atoms involved, such as the atomic coordinates, and the restraint target and weight. Since the restraints are built once at the beginning of refinement, while the atomic coordinates change during refinement, information about the atoms that participate in a restraint and the actual atomic coordinates is decoupled. The data structure that holds the restraint information is called a proxy, and it contains integer indices that can be used to access a particular atom in an array of atoms and therefore individual coordinates, along with some other restraint-specific information.

For example, the bond proxy contains a pair of integer indices (atom numbers in an atom list), the restraint target value and weight, the slack and limit values, a flag that specifies which function to use (least-squares or top-out) and a rotation matrix operator in case bonded atoms belong to different crystallographic symmetry copies, and an indicator to distinguish the origin of a particular restraint (e.g. covalent or hydrogen bond). Since each proxy object contains all the parameters necessary to define a particular restraint, each restraint can be parameterized individually.

The parallelity proxy is organized in a very similar manner to the planarity proxy, the difference being that the atom numbers are split and contained in two groups because they define two separate planes with an arbitrary number of atoms in each one.

Usually, one constructs multiple restraints of one type for a structure, so it is convenient to have a structure to store an array of proxies of one type and perform various operations on them. Similarly to other restraint types, for the parallelity restraints there is an array that can contain multiple proxy objects, called the shared_parallelity_proxy. Operations on proxy arrays include selection (this allows selection of a subset of restraints corresponding to a part of the structure), deletion (restraints can be deleted for a part of structure), and calulation of root mean-square deviations of the model from target values, absolute differences between the restraint target and actual model values, the restraint target value, gradients and so on.

5.4. Class to handle parallelity restraint

The parallelity restraint itself [equation (4)[link]] and its modifications [equations (7)[link] and (11)[link]] are implemented as a C++ class parallelity available in Python via boost.python from cctbx.geometry_restraints.parallelity. It has two constructors. One requires providing all parameters to define a parallelity proxy (such as target value, type of function etc.) and two sets of Cartesian atomic coordinates corresponding to the two planes. The other constructor takes an array with Cartesian coordinates of atoms and an instance of the parallelity proxy class that contains the rest of the required information.

All constructors contain a call to the init_deltas() function which calculates all necessary values, such as the overall restraint value and the gradients. This means that instantiation of a parallelity object leads to the immediate calculation of the target value and gradients for the restraint.

To simplify the handling of the slack functionality and to make it universal for harmonic and top-out potentials, the conditions of equation (11)[link] are evaluated first and θ0 adjusted accordingly before any other calculations. This makes it possible to avoid the gradient calculation if it is not necessary, for example when |θθ0| ≤ slack, and also to use all of the derivative equations for harmonic and top-out potentials without change.

In addition to the two constructors, the parallelity class has member functions residual() and gradients() that return the value of the restraint function and gradients with respect to the atomic coordinates. Since all necessary calculations are already performed in the constructor, these functions return precalculated values.

6. Discussion

Applying the general principles of refinement program construction, it is straightforward to build various stereochemical restraints. However, a careful analysis is required to choose an appropriate function for the restraint target. This choice is based on a function's ability to provide the desired behavior, both near the minimum and far from it. Possible singularities when calculating the target and its derivatives must also be taken into account. With a view to improving the refinement of atomic models, especially at low resolution, we have developed new restraints to control the angle between two approximately planar atomic groups. As a particular and more frequent case, this restraint imposes the groups to be parallel or near parallel. We have proposed a practical algorithm to calculate this target and its gradient. This algorithm has been implemented in the cctbx open-source library, making it available in PHENIX and to the broader crystallographic and macromolecular structure community. All calculation details are provided, making it easy to add other types of restraint in the future.

APPENDIX A

Other restraints on plane groups and their gradients

A1. Planarity targets

This restraint requires an atomic group to be as planar as possible. Calculation of the target uses exactly the same Steps ae described in §3[link] and Fig. 2[link] and differs only in Step f. As for the other steps, in order to avoid confusion below when calculating the gradients we introduce new parameters [{\hat \lambda _+} = {\lambda _+ }], [{\hat \lambda _- } = {\lambda _- }], [{\hat \lambda _0} = {\lambda _0}] instead of the eigenvalues of the inertia matrix (Appendix C[link]).

(Step f) The simplest target can be expressed as (Urzhumtsev, 1991[Urzhumtsev, A. G. (1991). Acta Cryst. A47, 723-727.])

[R_{\rm planarity} \left( {{\hat \lambda }_+ } \right ) = w {\hat \lambda _+ } , \eqno(23)]

which takes the value of the minimal eigenvalue already calculated at Step d. Another approach used currently (Blanc & Paciorek, 2001[Blanc, E. & Paciorek, W. (2001). J. Appl. Cryst. 34, 480-483.]) uses a different method of calculation but is essentially the same.

The restraint of equation (23)[link] does not take into account the number of atoms in the restrained group. This leads to the effect that the more populated group (i.e. the group with more atoms) will be forced to be more planar than the less populated group given the same weight. This may be corrected by normalizing by the number of atoms in the group:

[R_{\rm planarity} \left ( {{\hat \lambda }_ + } \right ) = {w \over {K_{\rm atoms}}} {\hat \lambda _+ } . \eqno(24)]

While equation (23)[link] expresses a type of absolute value of the `flatness' of the ellipsoid of inertia, other targets may estimate a relative value of such a flatness using several eigenvalues, for example

[R_{\rm planarity} \left ( {\hat \lambda _+ }, {\hat \lambda _0} \right ) = {{w {\hat \lambda _+ }} \over {\hat \lambda _0}} . \eqno(25)]

For this target the normalization with respect to the number of atoms is performed automatically.

A2. Parallel-distance target

Another geometric restraint that can be imposed on atomic models, especially those of nucleic acids, is requiring the distance between two parallel planes to be equal to a prescribed value ltarget. Its calculation follows the same algorithm presented above. We first pass through the same Steps ae, finishing as follows.

(Step f1) We repeat Step f1 of the parallelity restraint (§3.6[link]) and calculate normalized eigenvectors.

(Step f2) We then calculate the median (subscript med) direction between two normal vectors,

[{\bf N}_{\rm med} = {{{\bf n}_1 + {\bf n}_2} \over 2} , \eqno(26)]

(Step f3) and normalize it,

[{\bf n}_{\rm med} = {{{\bf N}_{\rm med}} \over {\left ( {\bf N}_{\rm med} \cdot {\bf N}_{\rm med} \right)^{1/2}}} . \eqno(27)]

(Step f4) Finally, we calculate the `oriented' distance between the best planes that pass through the two atomic groups as

[l = \left ( {\bf C}_2 - {\bf C}_1 \right) \cdot {\bf n}_{\rm med} , \eqno(28)]

(Step f5) and the target itself,

[R_{{\rm par} - {\rm dist}} \left ( {\bf N}_1, {\bf N}_2, {\bf C}_1, {\bf C}_2, \right) = w \left ( l^2 - l_{\rm target}^2 \right)^2 . \eqno(29)]

Target (29)[link] uses the square of the parameter from equation (28)[link] to remove the influence of the sign of the projection of (C2C1) on the unit vector nmed.

APPENDIX B

Eigenvalues and eigenvectors of the inertia matrix

Below, when describing the intermediate transformations of parameters at Steps d and e of the principal algorithm used in our program, the output parameters of each transformation are the input parameters for the next one. We start from the following input parameters:

[\eqalignno{ & C_x, C_y, C_z, \hat S_{XX}, \hat S_{XY}, \hat S_{XZ}, \hat S_{YX}, \hat S_{YY}, \hat S_{YZ}, \hat S_{ZX}, \hat S_{ZY}, \hat S_{ZZ}, a_S, \cr &b_S, c_S . &(30)}]

At Step d, all these parameters except aS, bS and cS remain unchanged and we omit them in the explicit list of parameters in this appendix.

(Step d1) We replace the initial equation (14)[link] by

[\left ( \lambda + {\textstyle{1 \over 3}} a_S \right )^3 + g_S \left ( \lambda + {\textstyle{1 \over 3}} a_S \right ) + h_S = 0 , \eqno(31)]

with the new parameters replacing aS, bS and cS

[\eqalign {g_S = & \, - \left (- {\textstyle{1 \over 3}} a_S^2 + b_S \right) = \textstyle {1 \over 3} a_S^2 - b_S, \cr {h_S} = & \, {\textstyle{2 \over 27}} a_S^3 - {\textstyle{1 \over 3}} a_S b_S + c_S, \cr {\hat a_S} = & \, {\textstyle{1 \over 3}} a_S , } \eqno(32)]

If the value of gS in equation (32)[link] is close to zero, or even slightly negative because of rounding errors, we assign gS = 0; when an atomic group is already more or less planar this will not be the case.

(Step d2) We calculate the cosine ξS of the characteristic angle 3γS as

[\eqalign{ \xi _S = & \, - {\textstyle{1 \over 2}} h_S \left ( {\textstyle{1 \over 3}} g_S \right )^{-3/2}, \cr {\hat g_S} = & \, 2 \left ( {\textstyle{1 \over 3}} g_S \right )^{1/2}, }\eqno(33)]

resulting in the output parameters [{\hat a_S}], [{\hat g_S}] and [{\xi _S}].

(Step d3) First, we calculate

[{\gamma _S} = {\textstyle{1 \over 3}} \arccos \,{\xi _S}, \eqno(34)]

and then

[{\tau _{S+}} = \cos \left ( {\gamma _S} + {\textstyle{2 \over 3}} \pi \right ), \,{\tau _{S-}} = \cos \left ( {\gamma _S} - {\textstyle{2 \over 3}} \pi \right ),\, {\tau _{S0}} = \cos \left ({\gamma _S} \right ). \eqno(35)]

The output parameters are [{\hat a_S}], [{\hat g_S}], [{\tau _{S+}}], [{\tau _{S-}}] and [{\tau _{S0}}].

(Step d4) The eigenvalues are

[{\lambda _+} = \hat g_S {\tau _{S+}} - \hat a_S, \,{\lambda _-} = \hat g_S{\tau _{S-}} - \hat a_S, \,{\lambda _0} = \hat g_S {\tau _{S0}} - \hat a_S, \eqno(36)]

where λ+λλ0. We then calculate the corresponding eigenvectors starting from the set of parameters

[\eqalignno{ & {C_x}, {C_y}, {C_z}, \hat S_{XX}, \hat S_{XY}, \hat S_{XZ}, \hat S_{YX}, \hat S_{YY}, \hat S_{YZ}, \hat S_{ZX}, \hat S_{ZY}, \hat S_{ZZ}, {\lambda _+}, \cr & {\lambda _-}, {\lambda _0} . &(37)}]

(Step e1) The eigenvector N+ corresponding to the minimal eigenvalue λ+ satisfies the condition

[\left ( {\hat {\bf S}} - {\lambda _+} {\bf I} \right) {\bf N}_+ = 0 . \eqno(38)]

This means that it is orthogonal to the three vectors tX = [( {\hat S_{XX}} - {\lambda _+}, {\hat S}_{XY}, {\hat S_{XZ}} )], tY = [( {\hat S_{YX}}, \hat S_{YY} - {\lambda _+}, \hat S_{YZ} )], tZ = [( {\hat S_{ZX}}, {\hat S}_{ZY} ,\hat S_{ZZ} - {\lambda _+} )]. In other words, N+ can be taken as the vector product of any two of them. To minimize rounding errors and avoid singular situations, for example when two chosen vectors are collinear or practically collinear to each other, we check all permutations

[\left ( {\bf t}_1 = {\bf t}_X, {\bf t}_2 = {\bf t}_Y \right ), \left ( {\bf t}_1 = {\bf t}_Y, {\bf t}_2 = {\bf t}_Z \right ), \left ( {\bf t}_1 = {\bf t}_Z, {\bf t}_2 = {\bf t}_X \right ) . \eqno(39)]

For each of these permutations we calculate the vector

[\eqalignno{ {\bf N}_+ = & \, {\bf t}_1 \times {\bf t}_2 \cr = & \, \left ( {t_{1Y}} {t_{2Z}} - {t_{1Z}} {t_{2Y}}, \,{t_{1Z}} {t_{2X}} - {t_{1X}} {t_{2Z}},\, {t_{1X}} {t_{2Y}} - {t_{1Y}} {t_{2X}} \right ) , &(40)}]

and accept the permutation for which the norm of equation (40)[link] is maximal.

(Step e2) After selecting vectors t1 and t2 we take the corresponding product of equation (40)[link] as the eigenvector to define the `best plane'.

The calculation of other eigenvectors (if necessary) follows the same procedure. Currently known geometric restraints, including parallelity ones, do not require more than one eigenvector that corresponds to the minimal eigenvalue λ+. However, considering other possible new restraints, the complete transition from atomic coordinates to plane parameters in the general case should include as output all three eigenvalues and all three corresponding eigenvectors

[{\lambda _+}, {\lambda _-}, {\lambda _0}, {\bf N}_+, {\bf N}_-, {\bf N}_0, {C_x}, {C_y}, {C_z} . \eqno(41)]

APPENDIX C

Gradient with respect to the coefficients of the characteristic equation

Recalculation of the gradient is applied to each pseudo-planar group, one by one, and does not depend on the particular plane-group restraint R. The input to this computational module is a set of partial derivatives of R with respect to the center of the group, to the three eigenvalues and to the coordinates of the three eigenvectors, as mentioned above. For a particular restraint, some of these derivatives, if not most of them, may be equal to zero.

(Step e2) Derivation of the eigenvector with respect to the t vectors.

For the eigenvector N = N+ calculated by equation (40)[link] we have

[\eqalign{ {{\partial R} \over {\partial {t_{1X}}}} = & \, {{\partial {N_X}} \over {\partial {t_{1X}}}} {{\partial R} \over {\partial {N_X}}} + {{\partial {N_Y}} \over {\partial {t_{1X}}}} {{\partial R} \over {\partial {N_Y}}} + {{\partial {N_Z}} \over {\partial {t_{1X}}}} {{\partial R} \over {\partial {N_Z}}} \cr = & \,0 - {t_{2Z}} {{\partial R} \over {\partial {N_Y}}} + {t_{2Y}} {{\partial R} \over {\partial {N_Z}}} , \cr {{\partial R} \over {\partial {t_{1Y}}}} = & \, - {t_{2X}} {{\partial R} \over {\partial {N_Z}}} + {t_{2Z}} {{\partial R} \over {\partial {N_X}}} , \cr \quad {{\partial R} \over {\partial {t_{1Z}}}} = & \, - {t_{2Y}} {{\partial R} \over {\partial {N_X}}} + {t_{2X}} {{\partial R} \over {\partial {N_Y}}} . } \eqno(42)]

The derivatives [{{\partial R} / {\partial {t_{2X}}}}], [{{\partial R} / {\partial {t_{2Y}}}}] and [{{\partial R} / {\partial {t_{2Z}}}}] differ from equation (42)[link] by swapping the indices 1 and 2 and inverting the sign.

(Step e1) Derivation of the eigenvector with respect to the eigenvalues and [{\hat {\bf S}}] elements.

The partial derivatives with respect to the [{\hat {\bf S}}] elements depend on the choice of vectors for t1 and t2 when calculating N, equation (40)[link]. For example, when (t1 = tX, t2 = tY)

[\eqalign{ {{\partial R} \over {\partial {{\hat S}_{XX}}}} = & \, {{\partial R} \over {\partial {t_{1X}}}}, \cr {{\partial R} \over {\partial {{\hat S}_{XY}}}} = & \, {{\partial R} \over {\partial {t_{1Y}}}}, \cr {{\partial R} \over {\partial {{\hat S}_{XZ}}}} = & \, {{\partial R} \over {\partial {t_{1Z}}}}, \cr {{\partial R} \over {\partial {{\hat S}_{YX}}}} = & \, {{\partial R} \over {\partial {t_{2X}}}}, \cr {{\partial R} \over {\partial {{\hat S}_{YY}}}} = & \, {{\partial R} \over {\partial {t_{2Y}}}}, \cr {{\partial R} \over {\partial {{\hat S}_{YZ}}}} = & \, {{\partial R} \over {\partial {t_{2Z}}}}, \cr {{\partial R} \over {\partial {{\hat S}_{ZX}}}} = & \, {{\partial R} \over {\partial {{\hat S}_{ZY}}}} = {{\partial R} \over {\partial {{\hat S}_{ZZ}}}} = 0, \cr {{\partial R} \over {\partial {\lambda _ + }}} = & \, - {{\partial R} \over {\partial {t_{1X}}}} - {{\partial R} \over {\partial {t_{2Y}}}} + {{\partial R} \over {\partial {{\hat \lambda }_ + }}} . } \eqno(43)]

The equations are similar for other choices of t1 and t2. It is important not to forget the direct partial derivative [{{\partial R} / {\partial {{\hat \lambda }_+}}}] if the target depends directly on the eigenvalue. The partial derivatives of (43)[link] with respect to other eigenvalues, for example when using targets such as equation (24)[link], are calculated similarly. In this work we do not consider the targets depending on eigenvectors other than N+; equations (42)[link] and (43)[link] may easily be generalized for more complicated targets in the future, but here we avoid giving these simple but lengthy expressions.

(Step d4) Derivation of eigenvalues with respect to [{\hat a_S}], [{\hat g_S}], [{\tau _{S+}}], [{\tau _{S-}}], [{\tau _{S0}}]. For each eigenvalue we obtain

[\eqalign{ {{\partial R} \over {\partial {{\hat a}_S}}} = & \, {{\partial {\lambda _ + }} \over {\partial {{\hat a}_S}}}{{\partial R} \over {\partial {\lambda _ + }}} + {{\partial {\lambda _ - }} \over {\partial {{\hat a}_S}}}{{\partial R} \over {\partial {\lambda _ - }}} + {{\partial {\lambda _0}} \over {\partial {{\hat a}_S}}}{{\partial R} \over {\partial {\lambda _0}}} \cr = & \, - {{\partial R} \over {\partial {\lambda _ + }}} - {{\partial R} \over {\partial {\lambda _ - }}} - {{\partial R} \over {\partial {\lambda _0}}}, \cr {{\partial R} \over {\partial {{\hat g}_S}}} = & \, {\tau _{S + }}{{\partial R} \over {\partial {\lambda _ + }}} + {\tau _{S - }}{{\partial R} \over {\partial {\lambda _ - }}} + {\tau _{S0}}{{\partial R} \over {\partial {\lambda _0}}}, \cr {{\partial R} \over {\partial {\tau _{S + }}}} = & \, {\hat g_S}{{\partial R} \over {\partial {\lambda _ + }}}, \cr {{\partial R} \over {\partial {\tau _{S - }}}} = & \, {\hat g_S}{{\partial R} \over {\partial \lambda _ - ^{}}}, \cr {{\partial R} \over {\partial {\tau _{S0}}}} = & \, {\hat g_S}{{\partial R} \over {\partial {\lambda _0}}}. } \eqno(44)]

Partial derivatives with respect to the matrix [{\hat {\bf S}}] elements were calculated above [equation (43)[link]].

(Step d3) Derivation of eigenvalues with respect to intermediate variables. Using the property of a derivative of the inverse function, we have

[\eqalignno { {{d \tau _{S + }} \over {d \xi _S}} = & \, {\left ({{d \xi _S} \over {d \tau _{S+}}} \right)^{-1} } = {\left [{{d \cos (3 \gamma)} \over {d \cos (\gamma + {{2\pi /3}})}} \right] ^{-1}} \cr = & \, {\left [{{d (4 \tau _{S+}^3 - 3 \tau _{S+})} \over {d \tau _{S+}}} \right] ^{-1}} = {1 \over {12 \tau _{S+}^2 - 3}}, &(45)}]

and similar equations for τS and τS0. Therefore

[{{\partial R} \over {\partial {\xi _S}}} = {1 \over {\left({12\tau _{S + }^2 - 3} \right)}}{{\partial R} \over {\partial {\tau _{S + }}}} + {1 \over {\left({12\tau _{S - }^2 - 3} \right)}}{{\partial R} \over {\partial {\tau _{S - }}}} + {1 \over {\left({12\tau _{S0}^2 - 3} \right)}}{{\partial R} \over {\partial {\tau _{S0}}}}. \eqno(46)]

(Step d2) To obtain the derivatives with respect to gS and hS, we first calculate

[\eqalign{{{{\rm d}{{\hat g}_S}} \over {{\rm d}{g_S}}} & = {\left ({ {1 \over 3}} g_S \right)^{ - 1/2}}, \quad{{{\rm d}{\xi _S}} \over {{\rm d}{g_S}}} = { {1 \over 4}} h_S {\left ({ {1 \over 3}} g_S \right)^{ - 5/2}}, \cr {{{\rm d}{\xi _S}} \over {{\rm d}{h_S}}} &= \left ({ - { {1 \over 2}}} \right) {\left ({ {1 \over 3}} g_S \right)^{ - 3/2}},} \eqno(47)]

with which

[{{\partial R} \over {d{g_S}}} = {{\partial {\xi _S}} \over {\partial {g_S}}}{{\partial R} \over {\partial {\xi _S}}} + {{\partial {{\hat g}_S}} \over {\partial {g_S}}}{{\partial R} \over {\partial {{\hat g}_S}}} , \quad {{\partial R} \over {\partial {h_S}}} = {{\partial {\xi _S}} \over {\partial {h_S}}}{{\partial R} \over {\partial {\xi _S}}}. \eqno (48)]

(Step d1) Derivation of eigenvalues with respect to the coefficients of the original characteristic equation. We can use

[\eqalign{ {{\partial {h_S}} \over {\partial {a_S}}} = & \, { {2 \over 9}} a_S^2 - { {1 \over 3}} b_S, \quad {{\partial {h_S}} \over {\partial {b_S}}} = - { {1 \over 3}} a_S, \quad {{\partial {h_S}} \over {\partial {c_S}}} = 1 \cr {{\partial {g_S}} \over {\partial {a_S}}} = & \, { {2 \over 3}} a_S, \quad {{\partial {g_S}} \over {\partial {b_S}}} = - 1, \quad {{\partial {{\hat a}_S}} \over {\partial {a_S}}} = { {1 \over 3}}, } \eqno(49)]

This results in

[\eqalign{ {{\partial R} \over {\partial {a_S}}} = & \, {{\partial {h_S}} \over {\partial {a_S}}}{{\partial R} \over {\partial {h_S}}} + {{\partial {g_S}} \over {\partial {a_S}}}{{\partial R} \over {\partial {g_S}}} + {{\partial {{\hat a}_S}} \over {\partial {a_S}}}{{\partial R} \over {\partial {{\hat a}_S}}} \cr = & \, \left({{{2a_S^2 - 3b_S^{}} \over 9}} \right){{\partial R} \over {\partial {h_S}}} + {{2a_S^{}} \over 3}{{\partial R} \over {\partial {g_S}}} + {1 \over 3} {{\partial R} \over {\partial {{\hat a}_S}}}, \cr {{\partial R} \over {\partial {b_S}}} = & \, {{\partial {h_S}} \over {\partial {b_S}}}{{\partial R} \over {\partial {h_S}}} + {{\partial {g_S}} \over {\partial {b_S}}}{{\partial R} \over {\partial {g_S}}} = - { {1 \over 3}} a_S {{\partial R} \over {\partial {h_S}}} - {{\partial R} \over {\partial {g_S}}}, \cr {{\partial R} \over {\partial {c_S}}} = & \, {{\partial {h_S}} \over {\partial {c_S}}}{{\partial R} \over {\partial {h_S}}} = {{\partial R} \over {\partial {h_S}}}, } \eqno(50)]

APPENDIX D

Calculation of the parallelity targets and their gradients with respect to the normal vectors

D1. Parallelity targets

Let each of the two pseudo-planar atomic groups be described by the Cartesian coordinates of the vector normal to its respective best plane (N1 and N2), not necessarily of unit length.

(Step f1) First we calculate the unit normal vectors as

[{\bf n}_1 = {{{\bf N}_1} \over {\left ( {\bf N}_1 \cdot {\bf N}_1 \right )^{1/2}}}, \quad {\bf n}_2 = {{\sigma {\bf N}_2} \over {\left ( {\bf N}_2 \cdot {\bf N}_2 \right)^{1/2}}}, \quad \sigma = \pm 1. \eqno(51)]

The correct sign for σ must be chosen so that the absolute value of the angle θ between the normal vectors is not larger than π/2.

(Step f2) Once the normalized vectors are known, we calculate the cosine and sine of the angle θ between them, considering that sinθ ≥ 0 because 0 ≤ θπ/2:

[C_\theta = \cos \left ( {\bf n}_1, {\bf n}_2 \right ) = {\bf n}_1 \cdot {\bf n}_2 = n_{1X} n_{2X} + n_{1Y} n_{2Y} + n_{1Z} n_{2Z}, \eqno(52)]

[\eqalignno{ {\bf n} = & \, \left ( n_X, n_Y, n_Z \right ) = {\bf n}_1 \times {\bf n}_2 \cr = & \, \left ( n_{1Y} n_{2Z} - n_{1Z} n_{2Y}, \,n_{1Z} n_{2X} - n_{1X} n_{2Z},\, n_{1X} n_{2Y} - n_{1Y} n_{2X} \right ), \cr &&(53)}]

[S_\theta = \left\| {\bf n} \right\| = \left ( n_X^2 + n_Y^2 + n_Z^2 \right )^{1/2} . \eqno(54)]

(Step f3) Finally, we calculate the value of the chosen target. Generally, one needs first to calculate the intermediate parameters

[\eqalign{ C_{D\theta} = & \, C_\theta \cos{\theta _0} + S_\theta \sin {\theta _0}, \cr S_{D\theta} = & \, S_\theta \cos {\theta _0} - C_\theta \sin {\theta _0} , } \eqno(55)]

which lead to the targets of equations (4)[link], (7)[link] and (8)[link] as

[R_{\rm parallelity} = w (1 - C_{D\theta}) , \eqno(56)]

[R_{\rm parallelity} = w {\Omega ^2} \left [ 1 - \exp {\left ( {{C_{D\theta} - 1} \over {\Omega ^2}} \right )} \right ] , \eqno(57)]

[R_{\rm parallelity} = 2w \left ( 1 - C_{D\theta}^2 \right ). \eqno(58)]

Equation (10)[link] for n > 2 may require, for example, using the De Moivre formula (Lial et al., 2008[Lial, M. L., Hornsby, J., Schneider, D. I., Callie, J. & Daniels (2008). College Algebra and Trigonometry, 4th ed., p.792. Boston: Pearson/Addison Wesley.]) CnDθ = Re(CDθ + iSDθ)n; we do not provide these expressions here.

Obviously, the particular and most frequent case θ0 = 0 can be considered as an independent target

[R_{\rm parallelity} = w \left ( 1 - {C_\theta } \right ), \eqno(59)]

which does not require equations (53)[link]–(55)[link].

D2. Gradient with respect to the normal vectors

(Step f3) Inversion of the last step of the target calculation provides the values of the gradients of the targets of equations (4)[link], (7)[link] and (8)[link] with respect to the intermediate parameters as

[{{\partial {R_{\rm parallelity}}} \over {\partial {C_{D\theta }}}} = - w, \quad {{\partial {R_{\rm parallelity}}} \over {\partial {S_{D\theta }}}} = 0, \eqno(60)]

[{{\partial {R_{\rm parallelity}}} \over {\partial {C_{D\theta }}}} = - w \exp { \left ( {{C_{D\theta} - 1} \over {\Omega ^2}} \right )} , \quad {{\partial {R_{\rm parallelity}}} \over {\partial {S_{D\theta }}}} = 0, \eqno(61)]

[{{\partial {R_{\rm parallelity}}} \over {\partial C_{D\theta }^{}}} = - 4wC_{D\theta } , \quad {{\partial {R_{\rm parallelity}}} \over {\partial S_{D\theta}}} = 0. \eqno(62)]

Then, using the chain rule,

[\eqalign{ {{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} = & \, {{\partial {R_{\rm parallelity}}} \over {\partial {C_{D\theta }}}}{{\partial {C_{D\theta }}} \over {\partial {C_\theta }}} + {{\partial {R_{\rm parallelity}}} \over {\partial {S_{D\theta }}}}{{\partial {S_{D\theta }}} \over {\partial {C_\theta }}} \cr = & \, {{\partial {R_{\rm parallelity}}} \over {\partial {C_{D\theta }}}} \cos{\theta _0} - {{\partial {R_{\rm parallelity}}} \over {\partial {S_{D\theta }}}} \sin{\theta _0}, \cr {{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}} = & \, {{\partial {R_{\rm parallelity}}} \over {\partial {C_{D\theta }}}}{{\partial {C_{D\theta }}} \over {\partial {S_\theta }}} + {{\partial {R_{\rm parallelity}}} \over {\partial {S_{D\theta }}}}{{\partial {S_{D\theta }}} \over {\partial {S_\theta }}} \cr = & \, {{\partial {R_{\rm parallelity}}} \over {\partial {C_{D\theta }}}} \sin{\theta _0} + {{\partial {R_{\rm parallelity}}} \over {\partial {S_{D\theta }}}} \cos{\theta _0}. } \eqno(63)]

(Step f2) At the next step, applying the chain rule, we obtain

[\eqalignno{ {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1X}}}}& = {{\partial {C_\theta }} \over {\partial {n_{1X}}}}{{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} + {{\partial {S_\theta }} \over {\partial {n_{1X}}}}{{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}} \cr = & \, {n_{2X}}{{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} + \left({{{\partial {n_X}} \over {\partial {n_{1X}}}}{{\partial {S_\theta }} \over {\partial {n_X}}} + {{\partial {n_Y}} \over {\partial {n_{1X}}}}{{\partial {S_\theta }} \over {\partial {n_Y}}} + {{\partial {n_Z}} \over {\partial {n_{1X}}}}{{\partial {S_\theta }} \over {\partial {n_Z}}}} \right) \cr & \, \times {{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}} \cr = & \, {n_{2X}}{{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} + \left({{{\partial {n_X}} \over {\partial {n_{1X}}}}{{{n_X}} \over {{S_\theta }}} + {{\partial {n_Y}} \over {\partial {n_{1X}}}}{{{n_Y}} \over {{S_\theta }}} + {{\partial {n_Z}} \over {\partial {n_{1X}}}}{{{n_Z}} \over {{S_\theta }}}} \right) \cr & \, \times {{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}} \cr = & \, {n_{2X}}{{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} + \left({{n_{2Y}}{{{n_Z}} \over {{S_\theta }}} - {n_{2Z}}{{{n_Y}} \over {{S_\theta }}}} \right){{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}}, & (64)}]

and similarly

[\eqalign{ {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Y}}}} = & \, {n_{2Y}}{{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} + \left({{n_{2Z}}{{{n_X}} \over {{S_\theta }}} - {n_{2X}}{{{n_Z}} \over {{S_\theta }}}} \right){{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}}, \cr {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Z}}}} = & \, {n_{2Z}}{{\partial {R_{\rm parallelity}}} \over {\partial {C_\theta }}} + \left({{n_{2X}}{{{n_Y}} \over {{S_\theta }}} - {n_{2Y}}{{{n_X}} \over {{S_\theta }}}} \right) {{\partial {R_{\rm parallelity}}} \over {\partial {S_\theta }}}. } \eqno(65)]

The expressions for [{{\partial {R_{\rm parallelity}}} / {\partial {n_{2X}}}}], [{{\partial {R_{\rm parallelity}}} / {\partial {n_{2Y}}}}] and [{{\partial {R_{\rm parallelity}}} / {\partial {n_{2Z}}}}] are similar to equations (64)[link] and (65)[link] with n1X, n1Y, n1Z instead of n2X, n2Y, n2Z in the right-hand parts and with a − sign instead of a + before the parentheses.

(Step f1) Taking the values obtained in equations (64)[link] and (65)[link], we calculate

[\eqalign{ {{\partial {R_{\rm parallelity}}} \over {\partial {N_{1X}}}} = & \, {{\partial {n_{1X}}} \over {\partial {N_{1X}}}} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1X}}}} + {{\partial {n_{1Y}}} \over {\partial {N_{1X}}}} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Y}}}} \cr & \, + {{\partial {n_{1Z}}} \over {\partial {N_{1X}}}} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Z}}}} \cr = & \, {{N_{1Y}^2 + N_{1Z}^2} \over {\left ( {\bf N}_1 \cdot {\bf N}_1 \right )^{3/2} }} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1X}}}} - {{N_{1X} N_{1Y}} \over {\left ( {\bf N}_1 \cdot {\bf N}_1 \right )^{3/2} }} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Y}}}} \cr & \, - {{N_{1X} N_{1Z}} \over {\left ( {\bf N}_1 \cdot {\bf N}_1 \right )^{3/2} }} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Z}}}} \cr = & \, \left \| {\bf N}_1 \right \|^{-3} \biggl [ \left ( N_{1Y}^2 + N_{1Z}^2 \right ) {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1X}}}} \cr & \, - N_{1X} N_{1Y} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Y}}}} - N_{1X} N_{1Z} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Z}}}} \biggr ] , } \eqno(66)]

where ||N1|| = (N1·N1)1/2. In the same way, we obtain

[\eqalignno{ {{\partial {R_{\rm parallelity}}} \over {\partial {N_{1Y}}}} = & \, {\left\| {\bf N}_1 \right\|^{-3}} \biggl [ \left (N_{1Z}^2 + N_{1X}^2 \right ) {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Y}}}} \cr & \, - N_{1Y} N_{1Z} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Z}}}} - N_{1Y} N_{1X} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1X}}}} \biggr ] , &(67)}]

[\eqalignno{ {{\partial {R_{\rm parallelity}}} \over {\partial {N_{1Z}}}} = & \, {\left\| {\bf N}_1 \right\|^{-3}} \biggl [ \left ( N_{1X}^2 + N_{1Y}^2 \right ) {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Z}}}} \cr & \, - N_{1Z} N_{1X} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1X}}}} - N_{1Z} N_{1Y} {{\partial {R_{\rm parallelity}}} \over {\partial {n_{1Y}}}} \biggr ] . &(68)}]

Expressions for [{{\partial {R_{\rm parallelity}}} /{\partial {N_{2X}}}}], [{{\partial {R_{\rm parallelity}}} / {\partial {N_{2Y}}}}] and [{{\partial {R_{\rm parallelity}}} / {\partial {N_{2Z}}}}] are obtained from equations (66)[link]–(68)[link] by substituting σ||N2||−3 for ||N1||−3 and replacing the index 1 by 2 for the rest.

Acknowledgements

This work was supported by the NIH (project 1P01 GM063210), by the PHENIX Industrial Consortium, and in part by the US Department of Energy under contract No. DE-AC02-05CH11231. AU thanks the French Infrastructure for Integrated Structural Biology (FRISBI) (grant No. ANR-10-INSB-05-01) and Instruct, part of the European Strategy Forum on Research Infrastructures (ESFRI).

References

First citationAbrahams, D. & Grosse-Kunstleve, R. W. (2003). C/C++ Users J. 21, 29–36.  Google Scholar
First citationAdams et al. (2010). Acta Cryst. D66, 213–221.  Google Scholar
First citationAfonine, P. V., Grosse-Kunstleve, R. W., Echols, N., Headd, J. J., Moriarty, N. W., Mustyakimov, M., Terwilliger, T. C., Urzhumtsev, A., Zwart, P. H. & Adams, P. D. (2012). Acta Cryst. D68, 352–367.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationAfonine, P. V., Headd, J. J., Terwilliger, T. C. & Adams, P. D. (2013). Comput. Crystallogr. Newsl. 4, 43–44.  Google Scholar
First citationAfonine, P. V. & Urzhumtsev, A. (2004). Acta Cryst. A60, 19–32.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationBaur, W. & Strassen, V. (1983). Theor. Comput. Sci. 22, 317–330.  CrossRef Web of Science Google Scholar
First citationBlanc, E. & Paciorek, W. (2001). J. Appl. Cryst. 34, 480–483.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationBrown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). Acta Cryst. D71, 136–153.  Web of Science CrossRef IUCr Journals Google Scholar
First citationDennis, J. E. Jr & Welsch, R. E. (1987). Commun. Stat. Simul. Comput. 7, 345–359.  CrossRef Web of Science Google Scholar
First citationEchols, N., Headd, J. J., Afonine, P. & Adams, P. D. (2010). Comput. Crystallogr. Newsl. 1, 12–17.  Google Scholar
First citationGrosse-Kunstleve, R. W. & Adams, P. D. (2004). IUCr Comput. Commission. Newsl. 4, 19–36.  Google Scholar
First citationGrosse-Kunstleve, R. W., Sauter, N. K., Moriarty, N. W. & Adams, P. D. (2002). J. Appl. Cryst. 35, 126–136.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationHeadd, J. J., Echols, N., Afonine, P. V., Grosse-Kunstleve, R. W., Chen, V. B., Moriarty, N. W., Richardson, D. C., Richardson, J. S. & Adams, P. D. (2012). Acta Cryst. D68, 381–390.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationHeadd, J. J., Echols, N., Afonine, P. V., Moriarty, N. W., Gildea. R. J. & Adams, P. D. (2014). Acta Cryst. D70, 1346–1356.  Web of Science CrossRef IUCr Journals Google Scholar
First citationJack, A. & Levitt, M. (1978). Acta Cryst. A34, 931–935.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationKim, K. M., Nesterov, Yu. E. & Cherkassky, B. V. (1984). Dokl. Acad. Nauk. SSSR, 275, 1306–1309.  Google Scholar
First citationKonnert, J. H. (1976). Acta Cryst. A32, 614–617.  CrossRef CAS IUCr Journals Web of Science Google Scholar
First citationLial, M. L., Hornsby, J., Schneider, D. I., Callie, J. & Daniels (2008). College Algebra and Trigonometry, 4th ed., p.792. Boston: Pearson/Addison Wesley.  Google Scholar
First citationLiu, H., Hexemer, A. & Zwart, P. H. (2012). J. Appl. Cryst. 45, 587–593.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationLunin, V. Yu. & Urzhumtsev, A. G. (1985). Acta Cryst. A41, 327–333.  CrossRef CAS Web of Science IUCr Journals Google Scholar
First citationMurshudov, G. N., Skubák, P., Lebedev, A. A., Pannu, N. S., Steiner, R. A., Nicholls, R. A., Winn, M. D., Long, F. & Vagin, A. A. (2011). Acta Cryst. D67, 355–367.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationParkhurst, J. M., Brewster, A. S., Fuentes-Montero, L., Waterman, D. G., Hattne, J., Ashton, A. W., Echols, N., Evans, G., Sauter, N. K. & Winter, G. (2014). J. Appl. Cryst. 47, 1459–1465.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationSauter, N. K., Hattne, J., Grosse-Kunstleve, R. W. & Echols, N. (2013). Acta Cryst. D69, 1274–1282.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationSmart, O. S., Womack, T. O., Flensburg, C., Keller, P., Paciorek, W., Sharff, A., Vonrhein, C. & Bricogne, G. (2012). Acta Cryst. D68, 368–380.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationUrzhumtsev, A. G. (1991). Acta Cryst. A47, 723–727.  CrossRef CAS Web of Science IUCr Journals Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoJOURNAL OF
APPLIED
CRYSTALLOGRAPHY
ISSN: 1600-5767
Follow J. Appl. Cryst.
Sign up for e-alerts
Follow J. Appl. Cryst. on Twitter
Follow us on facebook
Sign up for RSS feeds