Programming new geometry restraints: parallelity of atomic groups

Details are described of the calculation of new parallelity restraints recently introduced in cctbx and PHENIX.


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) 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 refinement of atomic coordinates, the terms 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 ISSN 1600-5767 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), to illustrate the general approach of adding new restraints in cctbx and PHENIX  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) 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, 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 refinement programs (Lunin & Urzhumtsev, 1985): (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;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 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, 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 n m ¼ N m ðN m Á N m Þ À1=2 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 x3.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.

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  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 1 2 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)].
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 x3 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). x4 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 pseudoplane groups differs only in the first step.

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.

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 k ) of the atomic positions: The overall calculation scheme. (a) The steps used to calculate parameters common to all targets. The letters a-f from top to bottom refer to the steps for calculating the target described in x3. The same steps in the direction from bottom to top refer to the steps for calculating the derivatives described in x4. 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.

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 XZ , 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 + À 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.

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 z ) of the center of the atomic group that are simply taken from the previous step, three equation coefficients a S , b S , c S , and nine elements of the inertia matrix S S XX ;Ŝ S XY ; :::;Ŝ S ZZ .

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.

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 z ) of its center, all three eigenvalues 0 + À 0 and all three eigenvectors N + , N À , N 0 of the inertia matrix.

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 f1) 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.

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) and obtain @R parallelity =@C and @R parallelity =@S (Appendix D).
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 z ) will not change until the last step.

Gradient with respect to the coefficients of the characteristic equation (Steps e and d)
As mentioned in xx3.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 S of equation (14), to the nine elementsŜ S XX ;Ŝ S XY , . . . ,Ŝ S ZZ of equation (16), and to the three coordinates (C x , C y , C z ) of the center of the atomic group that are taken, with no change, from the previous step.

Gradient with respect to the elements of the inertia matrix (Step c)
Using the partial derivatives of a S , b S , c S [equation (50) and similar expressions for the other six derivatives easily obtained from equation (18) by the corresponding cyclic substitution X ! Y ! Z ! X of the indices, e.g. to obtain @R=@S YZ we make the following substitutions in the second equation: XY ! YZ, YX ! ZY, ZZ ! XX, ZX ! XY, YZ ! ZX.

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,

Gradient of restraint with respect to the original coordinates (Step a)
Using the partial derivatives of equations (3) and (12) and the notation W j = w j ð and knowing that all cross-term derivatives are equal to 0, we obtain the required gradient with respect to the original atomic coordinates 5. Implementation of parallelity restraints in cctbx

Computational Crystallography Toolbox
The Computational Crystallography Toolbox (cctbx) is an open-source library with numerous functionalities necessary research papers 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 refinement  or diffraction data analysis (Sauter et al., 2013;Parkhurst et al., 2014). In particular, to support macromolecular refinement 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, chirality, 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., 2012Headd et al., , 2014. 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), 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.

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.

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.

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 paral-research papers lelity 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.

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 a-e described in x3 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 parameterŝ þ ¼ þ , À ¼ À , 0 ¼ 0 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 nucleic acids, is requiring the distance between two parallel planes to be equal to a prescribed value 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 (x3.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 f4) 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 .

research papers 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 S and c S remain unchanged and we omit them in the explicit list of parameters in this appendix. (Step d1) We replace the initial equation (14) by with the new parameters replacing a S , b S and 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 d2) We calculate the cosine S of the characteristic angle 3 S as resulting in the output parametersâ a S ,ĝ g S and S . ( Step d3) First, we calculate and then The output parameters areâ a S ,ĝ g S , Sþ , SÀ and S0 . (Step d4) The eigenvalues are þ ¼ĝ g S Sþ Àâ a S ; À ¼ĝ g S SÀ Àâ a S ; 0 ¼ĝ g S S0 Àâ a S ; ð36Þ where + À 0 . We then calculate the corresponding eigenvectors starting from the set of parameters C x ; C y ; C z ;Ŝ S XX ;Ŝ S XY ;Ŝ S XZ ;Ŝ S YX ;Ŝ S YY ;Ŝ S YZ ;Ŝ S ZX ;Ŝ S ZY ;Ŝ S ZZ ; þ ; Step e1) The eigenvector N + corresponding to the minimal eigenvalue + satisfies the condition This means that it is orthogonal to the three vectors t X = ðŜ S XX À þ ;Ŝ S XY ;Ŝ S XZ Þ, t Y = ðŜ S YX ;Ŝ S YY À þ ;Ŝ S YZ Þ, t Z = ðŜ S ZX ;Ŝ S ZY ;Ŝ S ZZ À þ Þ. 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 e2) 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 e2) Derivation of the eigenvector with respect to the t vectors.
For the eigenvector N = N + calculated by equation (40) we have The derivatives @R=@t 2X , @R=@t 2Y and @R=@t 2Z differ from equation (42) by swapping the indices 1 and 2 and inverting the sign. (Step e1) Derivation of the eigenvector with respect to the eigenvalues andŜ S elements.

research papers
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 f3) 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 R parallelity ¼ w 2 1 À exp C D À 1 2 ! ; ð57Þ 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 which does not require equations (53)-(55).