- 1. Introduction
- 2. Possible parallelity restraints and their gradient
- 3. Planarity targets calculation
- 4. Gradient calculation
- 5. Implementation of parallelity restraints in cctbx
- 6. Discussion
- A1. Planarity targets
- A2. Parallel-distance target
- D1. Parallelity targets
- D2. Gradient with respect to the normal vectors
- References

- 1. Introduction
- 2. Possible parallelity restraints and their gradient
- 3. Planarity targets calculation
- 4. Gradient calculation
- 5. Implementation of parallelity restraints in cctbx
- 6. Discussion
- A1. Planarity targets
- A2. Parallel-distance target
- D1. Parallelity targets
- D2. Gradient with respect to the normal vectors
- References

## research papers

## Programming new geometry restraints: parallelity of atomic groups

**Oleg V. Sobolev,**

^{a}^{*}Pavel V. Afonine,^{a}Paul D. Adams^{a,}^{b}and Alexandre Urzhumtsev^{c,}^{d}^{a}Lawrence Berkeley National Laboratory, One Cyclotron Road, MS64R0121, Berkeley, CA 94720, USA, ^{b}Department of Bioengineering, University of California Berkeley, Berkeley, CA 94720, USA, ^{c}Centre for Integrative Biology, IGBMC, CNRS-INSERM-UdS, 1 rue Laurent Fries, BP 10142, 67404 Illkirch, France, and ^{d}Département de Physique, Faculté des Sciences et des Technologies, Université de Lorraine, 54506 Vandoeuvre-lès-Nancy, France^{*}Correspondence e-mail: osobolev@lbl.gov

Improvements in structural biology methods, in particular crystallography and cryo-electron microscopy, have created an increased demand for the *a priori* information about model geometry that can be utilized in 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.

Keywords: restraints; atomic model refinement; parallel planes; *cctbx*; *PHENIX*; gradient calculation.

### 1. Introduction

The ) and Jack & Levitt (1978), this target is a weighted sum of several terms including *R*_{X} describing the model fit to the experimental data:

In the case of *R*_{n} 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; Headd *et al.*, 2012; Brown *et al.*, 2015) to complement the traditionally used restraints, such as bond lengths, bond angles, dihedral angles, non-bonded interactions, and planarity and restraints.

As an example, *cctbx* (Grosse-Kunstleve *et al.*, 2002), to illustrate the general approach of adding new restraints in *cctbx* and *PHENIX* (Adams *et al.*, 2010) and also to reinforce the general principles of crystallographic 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 will be presented elsewhere.

Each individual target in equation (1) 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 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, 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) 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 ):

programs (Lunin & Urzhumtsev, 1985(i)

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; Kim *et al.*, 1984). 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

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

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, **r**_{k,m} = (*x*_{k,m}, *y*_{k,m}, *z*_{k,m}), *k* = 1, …, *K _{m}*, in the same coordinate system. For each such group we need to define the corresponding `best plane'. A plane is defined by a vector

**N**

_{m}normal to it and by a point

**C**

_{m}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; Blanc & Paciorek, 2001; Brown

*et al.*, 2015), and these deviations may be weighted if desired. The `best plane' is the one that corresponds to the minimum of this value (

**N**

_{m}and

**C**

_{m}are variables),

where *w*_{k,m} is a weight (*e.g.* a unit value or an atomic mass) for the contribution of an atom with coordinates **r**_{k,m} to the target. Here and below, a centred dot `·' represents the scalar (dot) product of two vectors, and is the normalized vector **N**_{m}.

Two facts are useful (see, for example, Urzhumtsev, 1991). First, the minimum of equation (2) is reached when **C**_{m} is the weighted geometric center of the group

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

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

where θ = θ(**n**_{1}, **n**_{2}) is the angle between two unit vectors **n**_{1} and **n**_{2}, θ_{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) 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) will force the planes to form the prescribed angle between them. For example, as suggested by Brown *et al.* (2015), 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θ = **n**_{1}·**n**_{2} and sinθ = ||**n**_{1} × **n**_{2}||.

For small deviations from θ_{0}, restraint (4) behaves as the quadratic function *w*(θ − θ_{0})^{2}. However, using the quadratic function of the angle itself [see in particular equation (14) of Brown *et al.* (2015)] requires its derivatives, *i.e.* the derivative of arccos(**n**_{1}·**n**_{2}), 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

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

[equation (13) of Brown *et al.* (2015); the notation of the present article is used]. It imposes a rigid condition on the parallelity of two planes, **n**_{2} = **n**_{1}, 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; Murshudov *et al.*, 2011; Smart *et al.*, 2012; Headd *et al.*, 2012). A direct application of this idea results in

where Ω is a constant (Fig. 1). Using (θ − θ_{0})^{2} or (cosθ − cosθ_{0})^{2} in equation (7) 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:

or

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

A further useful modification of target function (4) 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) as

#### 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*). 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). For this reason, in §3 we provide details of these common steps.

Assuming we have a function that we can calculate numerically, the principle of efficient gradient calculation (Baur & Strassen, 1983; Kim *et al.*, 1984) 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). §4 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 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 **C**_{m} of the group, equation (3). 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 (*X*_{k}, *Y _{k}*,

*Z*) of the atomic positions:

_{k}and the Cartesian coordinates (*C _{x}*,

*C*,

_{y}*C*) of the weighted geometric center of the group.

_{z}#### 3.2. Inertia matrix (Step *b*)

For a given atomic group, its inertia matrix is

This matrix is symmetric, *S _{XY}* =

*S*,

_{YX}*S*=

_{YZ}*S*,

_{ZY}*S*=

_{ZX}*S*, but for convenience in subsequent calculations we refer to all nine of its elements. Matrix (13) is positive semidefinite, having always three non-negative eigenvalues, 0 ≤ λ

_{XZ}_{+}≤ λ

_{−}≤ λ

_{0}, with the corresponding eigenvectors orthogonal to each other. The result of this step contains nine elements of matrix (13) and three coordinates of the center [equation (12)] 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) are the roots of the cubic characteristic equation

with the coefficients

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

Thus, the new set of parameters contains 15 values: three coordinates (*C _{x}*,

*C*,

_{y}*C*) of the center of the atomic group that are simply taken from the previous step, three equation coefficients

_{z}*a*,

_{S}*b*,

_{S}*c*, and nine elements of the inertia matrix .

_{S}#### 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), as it gives analytically the minimal real eigenvalue that is required for the planarity and parallelity restraints. Appendix *B* gives an extended description in the notation of the present paper. The coordinates of the center of the atomic group [equation (12)] and the matrix elements [equation (16)] 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

with the eigenvalues λ obtained in the previous step. Appendix *B* gives technical details of the method used in our program. For an atomic group that is approximately planar, λ_{+} < λ_{−}, the solution of equation (17) 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*). 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 (*C _{x}*,

*C*,

_{y}*C*) of its center, all three eigenvalues 0 ≤ λ

_{z}_{+}≤ λ

_{−}≤ λ

_{0}and all three eigenvectors

**N**

_{+},

**N**

_{−},

**N**

_{0}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 **N**_{1} = **N**_{1+} and **N**_{2} = **N**_{2+} normal to their respective best planes, we may need to follow several intermediate steps (Appendix *D*). We start (Step *f*1) from normalization of the normal vectors, {**N**_{1}, **N**_{2}} {**n**_{1}, **n**_{2}}. 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 *f*2) the cosine and sine of the angle θ between them, which we denote *C*_{θ} and *S*_{θ}, respectively. Finally (Step *f*3), 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 *f*3 in Fig. 2) and obtain and (Appendix *D*).

Inverting the previous step (Step *f*2), we take as input the values of and (and not their analytical expressions), whatever the target *R*_{parallelity} is. As a result, we obtain (, , ) and (, , ). Finally, inverting the normalization (Step *f*1) we obtain (, , ) and (, , ).

Working with other restraints on plane groups, we can apply the same kind of procedure as Steps *f*3 to *f*1 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 **N**_{1+} and **N**_{2+}. 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 (*C _{x}*,

*C*,

_{y}*C*) will not change until the last step.

_{z}#### 4.2. Gradient with respect to the coefficients of the characteristic equation (Steps *e* and *d*)

As mentioned in §§3.4 and 3.5, different algorithms exist for calculating the eigenvalues and eigenvectors of the inertia matrix from the coefficients of the characteristic equation (14). According to the main scheme, we invert these algorithms to calculate the corresponding derivatives. Appendix *C* 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 *a _{S}*,

*b*,

_{S}*c*of equation (14), to the nine elements , …, of equation (16), and to the three coordinates (

_{S}*C*,

_{x}*C*,

_{y}*C*) of the center of the atomic group that are taken, with no change, from the previous step.

_{z}#### 4.3. Gradient with respect to the elements of the inertia matrix (Step *c*)

Using the partial derivatives of *a _{S}*,

*b*,

_{S}*c*[equation (50) in Appendix

_{S}*C*] and of the elements [equation (43) in Appendix

*C*] with respect to the elements of the matrix

**S**[equation (13)], we obtain

and similar expressions for the other six derivatives easily obtained from equation (18) by the corresponding cyclic substitution of the indices, *e.g.* to obtain we make the following substitutions in the second equation: , , , , .

### 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) or cryo-EM (Afonine *et al.*, 2013). *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 (Afonine *et al.*, 2012) or diffraction data analysis (Sauter *et al.*, 2013; Parkhurst *et al.*, 2014). In particular, to support macromolecular algorithms *cctbx* includes tools to handle various geometry restraints (Grosse-Kunstleve & Adams, 2004). The main set of geometry restraints includes those on covalent bonds and angles, dihedral (torsion) angles, non-bonded interactions, planar groups, and now parallelity. The overall restraint target, the sum in equation (1), is the sum of contributions from all restraint types

If additional information about the molecule or its environment is available and deemed necessary to use, additional terms can be added to equation (22) (Echols *et al.*, 2010; Headd *et al.*, 2012, 2014). Such additional information may be the regular conformations of a protein's secondary structure, or similarity to a previously solved structure. Technically, adding additional restraints means adding additional terms to equation (22), 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) 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

while the atomic coordinates change during 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 *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)] and its modifications [equations (7) and (11)] 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) 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 *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 *a*–*e* described in §3 and Fig. 2 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 , , instead of the eigenvalues of the inertia matrix (Appendix *C*).

(Step *f*) The simplest target can be expressed as (Urzhumtsev, 1991)

which takes the value of the minimal eigenvalue already calculated at Step *d*. Another approach used currently (Blanc & Paciorek, 2001) uses a different method of calculation but is essentially the same.

The restraint of equation (23) 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:

While equation (23) 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

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 *l*_{target}. Its calculation follows the same algorithm presented above. We first pass through the same Steps *a*–*e*, finishing as follows.

(Step *f*1) We repeat Step *f*1 of the parallelity restraint (§3.6) and calculate normalized eigenvectors.

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

(Step *f*3) and normalize it,

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

(Step *f*5) and the target itself,

Target (29) uses the square of the parameter from equation (28) to remove the influence of the sign of the projection of (**C**_{2} − **C**_{1}) on the unit vector **n**_{med}.

### 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:

At Step *d*, all these parameters except *a _{S}*,

*b*and

_{S}*c*remain unchanged and we omit them in the explicit list of parameters in this appendix.

_{S}(Step *d*1) We replace the initial equation (14) by

with the new parameters replacing *a _{S}*,

*b*and

_{S}*c*

_{S}If the value of *g _{S}* in equation (32) is close to zero, or even slightly negative because of rounding errors, we assign

*g*

_{S}= 0; when an atomic group is already more or less planar this will not be the case.

(Step *d*2) We calculate the cosine ξ_{S} of the characteristic angle 3γ_{S} as

resulting in the output parameters , and .

(Step *d*3) First, we calculate

and then

The output parameters are , , , and .

(Step *d*4) The eigenvalues are

where λ_{+} ≤ λ_{−} ≤ λ_{0}. We then calculate the corresponding eigenvectors starting from the set of parameters

(Step *e*1) The eigenvector **N**_{+} corresponding to the minimal eigenvalue λ_{+} satisfies the condition

This means that it is orthogonal to the three vectors **t**_{X} = , **t**_{Y} = , **t**_{Z} = . 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

For each of these permutations we calculate the vector

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

(Step *e*2) After selecting vectors **t**_{1} and **t**_{2} we take the corresponding product of equation (40) 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

### 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 *e*2) Derivation of the eigenvector with respect to the **t** vectors.

For the eigenvector **N** = **N**_{+} calculated by equation (40) we have

The derivatives , and differ from equation (42) by swapping the indices 1 and 2 and inverting the sign.

(Step *e*1) Derivation of the eigenvector with respect to the eigenvalues and elements.

The partial derivatives with respect to the elements depend on the choice of vectors for **t**_{1} and **t**_{2} when calculating **N**, equation (40). For example, when (**t**_{1} = **t**_{X}, **t**_{2} = **t**_{Y})

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

(Step *d*4) Derivation of eigenvalues with respect to , , , , . For each eigenvalue we obtain

Partial derivatives with respect to the matrix elements were calculated above [equation (43)].

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

and similar equations for τ_{S−} and τ_{S0}. Therefore

(Step *d*2) To obtain the derivatives with respect to *g _{S}* and

*h*, we first calculate

_{S}with which

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

This results in

### 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 (**N**_{1} and **N**_{2}), not necessarily of unit length.

(Step *f*1) First we calculate the unit normal vectors as

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 *f*2) Once the normalized vectors are known, we calculate the cosine and sine of the angle θ between them, considering that sinθ ≥ 0 because 0 ≤ θ ≤ π/2:

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

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

Equation (10) for *n* > 2 may require, for example, using the De Moivre formula (Lial *et al.*, 2008) *C*_{nDθ} = Re(*C*_{Dθ} + *iS*_{Dθ})^{n}; we do not provide these expressions here.

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

#### D2. Gradient with respect to the normal vectors

(Step *f*3) Inversion of the last step of the target calculation provides the values of the gradients of the targets of equations (4), (7) and (8) with respect to the intermediate parameters as

Then, using the chain rule,

(Step *f*2) At the next step, applying the chain rule, we obtain

and similarly

The expressions for , and are similar to equations (64) and (65) with *n*_{1X}, *n*_{1Y}, *n*_{1Z} instead of *n*_{2X}, *n*_{2Y}, *n*_{2Z} in the right-hand parts and with a − sign instead of a + before the parentheses.

(Step *f*1) Taking the values obtained in equations (64) and (65), we calculate

where ||**N**_{1}|| = (**N**_{1}·**N**_{1})^{1/2}. In the same way, we obtain

Expressions for , and are obtained from equations (66)–(68) by substituting σ||**N**_{2}||^{−3} for ||**N**_{1}||^{−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

Abrahams, D. & Grosse-Kunstleve, R. W. (2003). *C/C++ Users J.* **21**, 29–36. Google Scholar

Adams *et al.* (2010). *Acta Cryst.* D**66**, 213–221. Google Scholar

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.* D**68**, 352–367. Web of Science CrossRef CAS IUCr Journals Google Scholar

Afonine, P. V., Headd, J. J., Terwilliger, T. C. & Adams, P. D. (2013). *Comput. Crystallogr. Newsl.* **4**, 43–44. Google Scholar

Afonine, P. V. & Urzhumtsev, A. (2004). *Acta Cryst.* A**60**, 19–32. Web of Science CrossRef CAS IUCr Journals Google Scholar

Baur, W. & Strassen, V. (1983). *Theor. Comput. Sci.* **22**, 317–330. CrossRef Web of Science Google Scholar

Blanc, E. & Paciorek, W. (2001). *J. Appl. Cryst.* **34**, 480–483. Web of Science CrossRef CAS IUCr Journals Google Scholar

Brown, A., Long, F., Nicholls, R. A., Toots, J., Emsley, P. & Murshudov, G. (2015). *Acta Cryst.* D**71**, 136–153. Web of Science CrossRef IUCr Journals Google Scholar

Dennis, J. E. Jr & Welsch, R. E. (1987). *Commun. Stat. Simul. Comput.* **7**, 345–359. CrossRef Web of Science Google Scholar

Echols, N., Headd, J. J., Afonine, P. & Adams, P. D. (2010). *Comput. Crystallogr. Newsl.* **1**, 12–17. Google Scholar

Grosse-Kunstleve, R. W. & Adams, P. D. (2004). *IUCr Comput. Commission. Newsl.* **4**, 19–36. Google Scholar

Grosse-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

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.* D**68**, 381–390. Web of Science CrossRef CAS IUCr Journals Google Scholar

Headd, J. J., Echols, N., Afonine, P. V., Moriarty, N. W., Gildea. R. J. & Adams, P. D. (2014). *Acta Cryst.* D**70**, 1346–1356. Web of Science CrossRef IUCr Journals Google Scholar

Jack, A. & Levitt, M. (1978). *Acta Cryst.* A**34**, 931–935. CrossRef CAS IUCr Journals Web of Science Google Scholar

Kim, K. M., Nesterov, Yu. E. & Cherkassky, B. V. (1984). *Dokl. Acad. Nauk. SSSR*, **275**, 1306–1309. Google Scholar

Konnert, J. H. (1976). *Acta Cryst.* A**32**, 614–617. CrossRef CAS IUCr Journals Web of Science Google Scholar

Lial, 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

Liu, H., Hexemer, A. & Zwart, P. H. (2012). *J. Appl. Cryst.* **45**, 587–593. Web of Science CrossRef CAS IUCr Journals Google Scholar

Lunin, V. Yu. & Urzhumtsev, A. G. (1985). *Acta Cryst.* A**41**, 327–333. CrossRef CAS Web of Science IUCr Journals Google Scholar

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.* D**67**, 355–367. Web of Science CrossRef CAS IUCr Journals Google Scholar

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. Web of Science CrossRef CAS IUCr Journals Google Scholar

Sauter, N. K., Hattne, J., Grosse-Kunstleve, R. W. & Echols, N. (2013). *Acta Cryst.* D**69**, 1274–1282. Web of Science CrossRef CAS IUCr Journals Google Scholar

Smart, O. S., Womack, T. O., Flensburg, C., Keller, P., Paciorek, W., Sharff, A., Vonrhein, C. & Bricogne, G. (2012). *Acta Cryst.* D**68**, 368–380. Web of Science CrossRef CAS IUCr Journals Google Scholar

Urzhumtsev, A. G. (1991). *Acta Cryst.* A**47**, 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.