Chapter 25.2. Programs and program systems in wide use

Heavy-atom-based phasing by the methods of isomorphous replacement (Green et al., 1954) and/or anomalous scattering (Pepinsky & Okaya, 1956) are initiated by reading one or more `scaled' files into the program PHASIT, along with estimates of the heavy-atom or anomalous-scatterer positional, occupancy and thermal parameters. Each input file can contain either isomorphous-replacement, derivative anomalous-scattering or native anomalous-scattering data. MAD (multiple-wavelength anomalous diffraction) data are treated as both isomorphous and anomalous-scattering data, in which case one simply inputs the scattering-factor differences (real part) appropriate for the wavelengths comprising the `isomorphous' data sets and the actual scattering factors (imaginary part) appropriate for the `anomalous' data sets. All possible combinations of isomorphous-replacement data, conventional anomalous-scattering data and MAD data are allowed and can be used simultaneously during phasing.

Solvent flattening and negative-density truncation are carried out following the strategy developed by Wang (1985); however, a reciprocal-space equivalent of the automated solvent-masking procedure is used (Furey & Swaminathan, 1997; Leslie, 1987). In addition, during solvent-mask construction all density near heavy-atom sites is automatically ignored, leading to more accurate masks. The complete process is fully automated and carries out three solvent-mask iterations with at least 16 solvent flattening and phase combination cycles. Optionally, an arbitrary number of additional cycles can be carried out for phase extension. A program is provided to interactively examine or edit the solvent mask or to create the mask by hand if desired.

Noncrystallographic (NC) symmetry averaging cases (Rossmann & Blow, 1963; Bricogne, 1974) are treated in direct space by operating on `submaps', i.e. arbitrary regions in an electron-density map encompassing all of the molecules to be averaged that are unique by true crystal symmetry. Many of the averaging programs were derived from routines originally written by W. Hendrickson & J. Smith and have been described earlier (Bolin et al., 1993), but they were substantially rewritten for incorporation into the PHASES package. Programs are supplied to: generate and examine the required submaps; refine the NC symmetry operators; interactively create averaging envelope masks; average density within the envelopes; convert the submaps to full-cell maps; invert the modified maps; and combine the phases with those from another source. An automated procedure is provided to carry out a specified number of averaging and phase combination cycles in addition to solvent flattening and negative-density truncation. This procedure allows for a gradual phase extension, if desired, extending by one reciprocal-lattice point in each direction for a given number of cycles.

Several programs are included to carry out partial structure phase combination with a variety of weighting options as an aid to structure completion. If density modification (solvent flattening or negative-density truncation and/or NC symmetry averaging) is performed, then phase (and amplitude) extension is also possible by manual or automated procedures.

The current package includes 44 Fortran programs and one C subroutine, with the C subroutine used only to provide an interface between the Fortran programs and standard X-Window graphics-library routines. All programs communicate only through files with a simple common format. For the major programs, memory is allocated from a single large one-dimensional array which gets partitioned as required for each problem at run time. This greatly simplifies redimensioning if needed for very large problems, since at most only two lines of code need to be changed. All source code is provided, along with compilation procedures or shell scripts appropriate for most workstations, including Silicon Graphics, Sun, IBM R6000, ESV and DEC Alpha AXP (both OSF and OpenVMS). A flow chart illustrating the major programs and data flow for common phasing procedures is given in Fig. 25.2.1.1.

Figure 25.2.1.1| top | pdf | Flow chart for the major phasing path encompassing native and derivative scaling, heavy-atom-based phasing, solvent flattening, negative-density truncation, and phase combination. Boxed entries represent programs while lines represent files. Optional paths for noncrystallographic symmetry averaging and phase extension are included by considering the additional programs offset from the main path by dashed lines.

Vital data common to nearly all calculations, such as the cell dimensions, lattice type and space-group symmetry, are entered only once in a single `parameter file'. All interactive programs prompt for the name of this file and for batch programs it is to be supplied on the first input line. The parameter file can also optionally contain the name of a `running' log file. If used, the running log file is opened in `append' mode by each program in the package, and a copy of all screen or printed output is added to the file along with a time and date stamp indicating what program was run and when. This allows the user to maintain a complete history of all calculations and results on a given problem in a single, chronologically accurate file.

With this method, the derivative scattering, on average, is made equal to the native scattering by plotting with the averages taken in corresponding resolution shells. A least-squares fit of a straight line to the plot yields a slope equal to (twice the difference between overall isotropic temperature parameters for derivative and native data sets) and an intercept of ln . From these values, the derivative data are put on the scale of the native by multiplying each derivative amplitude by

With this option, applied after relative Wilson scaling, the unique parameters of a symmetric 3 × 3 scaling tensor S are determined by two cycles of least-squares minimization of with respect to S, where is a weighting factor, and , , are direction cosines of the reciprocal-lattice vector expressed in an orthogonal system. The derivative data are then placed on the scale of the native by multiplying each derivative amplitude by the appropriate S.

With this option, again applied after relative Wilson scaling, a scale factor for each reflection is also determined by minimizing equation (25.2.1.3) with respect to S, but here S is a scalar and the summation is taken only over neighbouring reflections within a sphere centred on the reflection being scaled. The sphere radius is initially set to include roughly 125 neighbours, and the scale factor is accepted if at least 80 are actually present. If insufficient neighbours are available, then the sphere size is increased incrementally and the process repeated until a preset maximum radius is encountered. If the maximum is reached, the process terminates with the message that the data set is too sparse for local scaling. Scaling is achieved by multiplying each derivative amplitude by the appropriate S.

Rejection of outliers is often desirable, as erroneously large isomorphous or anomalous differences can lead to streaks in difference-Patterson maps and complicate identification of heavy-atom or anomalous-scatterer sites. The interactive program TOPDEL facilitates identification and rejection of such outliers while selecting reflections for use in difference-Patterson calculations. An input `scaled' file is read in, and user-supplied resolution and cutoffs are applied. The data are then sorted in descending order of magnitude of ΔF (either isomorphous or anomalous differences) and the largest differences are listed for examination. The user is then prompted to determine which, if any, of the large differences are to be rejected as outliers and to determine what percentage of the remaining largest differences are to be used in the Patterson-map synthesis. The appropriate Fourier-coefficient file is then created.

Selected regions of an electron-density map that are useful for NC symmetry applications can be extracted from the full-cell maps produced by FSFOUR with the programs EXTRMAP (batch) or MAPVIEW (interactive). The `submap' regions can cover any arbitrary volume and cross multiple cell edges if desired.

Programs MAPORTH and SKEW (both run in batch mode) are provided to modify submaps, as modification is sometimes useful or required with NC symmetry applications. MAPORTH simply converts the map to correspond to an orthogonal grid, which simplifies refinement of NC symmetry operators. SKEW also converts the map to an orthogonal grid, but changes the axis directions such that the new b axis can be arbitrarily oriented. This is useful in NC symmetry applications where one may want to examine maps looking directly down the NC symmetry rotation axis. Both programs compute density values at the new grid points by using a 64-point cubic spline interpolation and can also orthogonalize or skew masks to maintain correspondence with the modified submaps.

Program GMAP (interactive) is used to extract any region from a FSFOUR map, possibly crossing multiple cell edges, and convert it to a form directly readable by the external interactive graphics programs TOM [SGI version of FRODO (Jones, 1978)], O (Jones et al., 1991) or CHAIN (Sack, 1988). In addition to the output map file, one may also output a corresponding skeleton (Greer, 1974) file (for TOM) or skeleton data block (for O) to facilitate chain tracing.

Program PSRCH (batch) lists the largest peaks in a Fourier map and is useful in identifying additional heavy-atom or anomalous-scatterer sites from a map phased by a tentative model. Either positive or negative peaks can be listed, with the latter sometimes useful in MAD phasing applications, depending on the assignment of `native' and `derivative' data sets. Only unique peaks are listed, and the peak positions are interpolated from the map.

Phasing by heavy-atom-based methods (isomorphous replacement and/or anomalous scattering) begins when one or more `scaled' data sets are input to the program PHASIT (batch). User-specified rejection criteria are first applied to each data set, and structure factors corresponding to the heavy-atom or anomalous-scatterer substructure are computed from where is the occupancy, is the (possibly complex) scattering factor, is the isotropic temperature parameter and , , are the fractional coordinates of the jth atom. Sc is a scale factor relating the calculated structure factor (absolute scale) to the scale of the observed data. The summation is taken over all heavy atoms or anomalous scatterers in the unit cell. Alternatively, anisotropic temperature parameters can be used for each atom if desired. A subset of reflections comprising all centric data (plus the largest 25% of the isomorphous or anomalous differences if there are insufficient centric data) is selected and used to estimate Sc by a least-squares fit to the observed differences. Initial estimates of the `standard error' E (expected lack of closure) are determined from this subset as a function of F magnitude, treating centric and acentric data separately. SIR (single isomorphous replacement) or SAS (single-wavelength anomalous scattering) phase probability distributions are given by where the lack of closure is defined by for isomorphous-replacement data and for anomalous-scattering data, with the and − superscripts denoting members of a Bijvoet pair, and with φ denoting the protein phase, and and denoting the heavy-atom structure-factor amplitude and phase, respectively. The distributions, however, are cast in the A, B, C, D form (Hendrickson & Lattman, 1970). After all input data sets are processed in this manner, the individual phase probability distributions for common reflections are combined via with k as a normalization constant and the sums taken over all contributing data sets. The resulting combined distributions are then integrated to yield a centroid phase and figure of merit for each reflection. The standard error estimates, E, as a function of structure-factor magnitude are then updated for each data set, this time using all reflections and a probability-weighted average over all possible phase values for the contribution from each reflection (Terwilliger & Eisenberg, 1987a,b). With these updated standard error estimates, the individual SIR and/or SAS phase probability distributions are recomputed for all reflections and combined again to yield an improved centroid phase and figure of merit for each reflection. The resulting phases, figures of merit and probability distribution information are then available for use in map calculations or for further parameter or phase refinement. This method is used to produce MIR (multiple isomorphous replacement), SIRAS (single isomorphous replacement with anomalous scattering) MIRAS (multiple isomorphous replacement with anomalous scattering) and MAD phases as well as other possible phase combinations.

Structure-factor amplitudes and phases for a macromolecular structure can be computed directly from atomic coordinates corresponding to a tentative model with the programs PHASIT and GREF (both run as batch processes). This allows one to obtain structure-factor information from an input model typically derived from a partial chain trace or from a molecular-replacement solution. Equation (25.2.1.5) is used, but this time the sum is taken over all known atoms in the cell, and the scale factor is refined by least squares against the native amplitudes rather than against the magnitudes of isomorphous or anomalous differences. The computed structure factors may be used directly for map calculations, including `omit' maps, or for combination with other sources of phase information. One can output probability distribution information for the calculated phases, if desired, as well as coefficients for various Fourier syntheses, including those using sigma_A weighting (Read, 1986) for the generation of reduced-bias native or difference maps.

For the purpose of improving phases by density-modification methods, such as solvent flattening, negative-density truncation and/or NC symmetry averaging, one must compute structure factors by Fourier inversion of an electron-density map rather than from atomic coordinates. The program MAPINV (batch) is a companion program to FSFOUR and carries out this inverse Fourier transform. It accepts a full-cell map in FSFOUR format and inverts it to produce amplitudes and phases for a selected set of reflections when given the target range of Miller indices. A variable-radix 3D fast Fourier transform algorithm is used. Optionally, the program can modify the density prior to inversion by truncation below a cutoff and/or by squaring the density values. Other types of density modification are handled by different programs in the package and are carried out prior to running MAPINV. The indices, calculated amplitude and phase are written to a file for each target reflection.

The simplest procedure is to refine against the magnitudes of isomorphous or anomalous structure-factor amplitude differences, which can be carried out with the program GREF (batch mode). In this case, one minimizes for isomorphous-replacement data or for anomalous-scattering data with respect to the desired parameters contributing to , where is a weighting factor. For anomalous-scattering data, only the imaginary component of the scattering factors is used during the structure-factor calculation. For isomorphous-replacement data, the summation is taken only over centric reflections, plus the strongest 25% of differences for acentric reflections if insufficient centric data are present. For anomalous-scattering data, the summation is taken only over the strongest 25% of Bijvoet differences. An advantage of these methods is that only data from the derivative being refined are used (plus the native with isomorphous data), hence there is no possibility of feedback between other derivatives which may not be truly independent. A disadvantage is that, apart for the centric reflections, the target value in the minimization is only an approximation to the true . The accuracy of this approximation is improved by restricting the summations to the strongest differences.

An alternative procedure available in the program PHASIT (batch) is to refine against the observed derivative amplitudes. In this case, one minimizes the `lack of closure' (now based on amplitudes instead of intensities) with respect to the desired parameters contributing to , including the derivative-to-native scaling parameters. In all cases, the calculated derivative amplitudes are obtained from equation (25.2.1.9). To use this procedure, one must have an estimate of the protein phase φ. Several variations of this method, all available in PHASIT, are described below and are generally referred to as `phase refinement'.

The files produced by PHASIT for isomorphous data sets contain the information needed to produce the Fourier coefficients where and are the calculated heavy-atom structure-factor amplitude and phase, respectively, and is computed from where and are the current derivative and native phases, respectively. These coefficients are more accurate than using simple isomorphous difference magnitudes to approximate and can be computed once phasing has begun, since estimates of the required phase differences are then available. Alternatively, the program MRGDF (interactive) can be used to produce Fourier coefficients of the form where m is the current figure of merit. This method suffers somewhat as phase differences are ignored, but it has the advantage that the amplitude difference does not necessarily involve any derivative previously used in the computation of . If amplitudes from a new derivative and from the native are used, then peaks in the resulting `cross-difference' Fourier synthesis for the new derivative will automatically correspond to the same origin and hand as prior sites used in the phasing process, although the hand may still be incorrect. Finally, GREF can be used to generate the Fourier coefficients with the second set producing a map similar to that obtained using equation (25.2.1.16). Both coefficient sets in equation (25.2.1.19) are lacking in that the phase difference is ignored, but the second set [and also those in equation (25.2.1.16)] has the advantage that heavy-atom sites already in the model are subtracted away, allowing any remaining minor sites to stand out in the resulting map.

The files produced by PHASIT for anomalous-scattering data sets contain the information needed to produce the Fourier coefficients or where is the protein phase used when computing . The coefficients in equation (25.2.1.20) correspond to a conventional Bijvoet difference Fourier map, which should show large positive peaks at the locations of anomalous-scattering sites when the hand is correct. The coefficients in equation (25.2.1.21) correspond to the case in which contributions from known anomalous scatterers are subtracted out. As in the isomorphous-replacement case, a program MRGBDF (interactive) is also provided to generate the Fourier coefficients where the Bijvoet difference does not necessarily have to come from a derivative used in the phasing. If the difference doesn't come from a derivative used in phasing, then a `cross-Bijvoet difference' Fourier map is obtained, which should produce large positive peaks at anomalous-scatterer locations in the new derivative when the original hand is correct. Additionally, GREF can be used to generate the Fourier coefficients where and are the heavy-atom structure-factor amplitude and phase, used when computing . These coefficients can also be used to identify additional anomalous-scatterer sites, but they are insensitive to the hand. As in equation (25.2.1.21), if the second set in equation (25.2.1.23) is used, then contributions from anomalous scatterers already included in the phasing will be subtracted out.

Finally, the program HNDCHK (interactive) is provided to determine the enantiomorph by examination of a Bijvoet difference Fourier map. One inputs the map along with the anomalous-scatterer positions used in the phasing. The program then uses a 64-point cubic spline interpolation algorithm to obtain the density precisely at the input coordinates and also at coordinates related to them by a centre of symmetry. If the input heavy-atom configuration had the correct hand, large positive peaks should occur exactly at the input locations. If the hand is incorrect, even larger negative peaks occur at the true positions, i.e. those related to the input positions by a centre of symmetry.

Solvent-mask construction follows the procedure suggested by Wang (1985), with the exception that electron density in the vicinity of heavy-atom sites is temporarily ignored during the mask-building process. This allows one to use a tight solvent mask, which maximizes the phasing power of the method while preventing artificial extension of the protein envelope into the solvent region in the vicinity of surface-bound heavy-atom sites. Failure to do this has occasionally been found to deplete the protein region elsewhere to compensate for the incorrectly extended region.

Once a solvent mask is constructed, solvent flattening and negative-density truncation is carried out using the program BNDRY. An electron-density map and corresponding mask are input along with an empirical constant S, which is used to estimate the value of on the scale of the input map. The estimation follows the procedure of Wang (1985) and is based on the assumption that for typical solvent conditions and proteins not containing heavy metals, the ratio of mean solvent electron density to maximum protein electron density is constant, although for phasing purposes the optimum values are resolution-dependent. Typical values of S are supplied in the package. One simply couples the value of S taken from known structures with density values obtained from the experimental maps to estimate on the appropriate (but unknown) map scale by solving the equation for . Once the estimate of is obtained, solvent flattening and negative-density truncation are carried out simultaneously by resetting all map values according to the relationships

The scaled data are sorted into bins according to d spacing, and a three-term polynomial is fitted to the mean values of as a function of resolution. For each reflection, a unimodal phase probability distribution is constructed using a modification (Bricogne, 1976) of the Sim (1959) weighting scheme via where the average in the appropriate resolution range is determined from the polynomial. This distribution is cast in the A, B, C, D form with Phase combination with the anchor set then proceeds according to equation (25.2.1.10), and the combined distributions are integrated to give a new phase and figure of merit for each reflection.

As an alternative to the procedure above, in the BNDRY program the weights, W, used when constructing the unimodal probability distributions in equations (25.2.1.30) can be computed according to where and are normalized structure-factor amplitudes for the observed and calculated structure factors, respectively, and is determined by the procedure described by Read (1986). For acentric reflections, equation (25.2.1.31) is used whereas for centric reflections, W is one half the value given by equation (25.2.1.31).

Normally, the distributions constructed for the calculated phases are combined with those for the anchor set with full weight in equation (25.2.1.10). However, in BNDRY, one can supply a damping factor in the range 0–1 to down-weight the contributions of the anchor set. The damping factor simply multiplies the distribution coefficients such that a factor of 1 (default) indicates no damping, and values less than one place more emphasis on the map-inverted or partial structure phases. If set to zero, the calculated phases are accepted as they are, since there is effectively no phase combination with the anchor set.

If phase extension is requested during the phase combination step, an additional file (prepared by the interactive program MISSNG) is also supplied to the BNDRY program. This file contains unique reflections absent from the anchor set but for which observed amplitudes (and possibly phase probability distribution coefficients) are available. Phase combination then proceeds exactly as above, except that for any extended reflections lacking phase probability information, the calculated phases are accepted as they are. Phase extension is required when phasing purely by SAS methods as it is the only way to phase centric reflections. As a final option, phase and amplitude extension is possible, in which case both the calculated amplitude and phase are accepted as they are for reflections having only indices provided on the extension file. This is sometimes desirable to include low-resolution reflections that may have been obscured by the beam stop.

NC symmetry operators are specified in terms of the parameters φ, ψ, χ, , , and t, which refer to a Cartesian coordinate system in Å, obtained by orthogonalization of the unit cell as in the Protein Data Bank (Bernstein et al., 1977). The angles φ and ψ determine the direction of the NC rotation axis, while χ determines the amount of rotation about it. , and are coordinates of a point through which the rotation axis passes, and t is a post-rotation translation parallel to the rotation axis. The relationships between the angles, orthogonal reference axes X, Y, Z and the unit cell are given in Fig. 25.2.1.2. Coordinates for a pair of points related by NC symmetry are then expressed in the orthogonal system by where and are three-element column vectors containing coordinates for the related points, is a 3 × 3 rotation matrix derived from the angles, O is a three-element column vector containing coordinates for a point through which the rotation axis passes, t is the post-rotation translation scalar in Å and is a three-element column vector containing direction cosines of the rotation axis. This type of parameterization simplifies transfer of information from self-rotation functions, which are usually calculated in spherical polar angles anyway, and also makes obvious pseudo-space-group symmetry type operations such as pseudo-screw axes. For convenience, a program O_TO_SP is provided to convert from a 3 × 3 rotation matrix and 1 × 3 column vector representation of the NC symmetry operation, as used in some programs, to the parameters described here.

Figure 25.2.1.2| top | pdf | Relationships between noncrystallographic symmetry rotation axis direction, orthogonal reference system axes X, Y, Z and crystallographic axes. The X axis is aligned with the crystal a. The Y axis is parallel to . The Z axis is parallel to , i.e. . ψ is the angle between the NC rotation and +Y axes. φ is the angle between the projection of the NC rotation axis in the XZ plane and the +X axis, with +φ counterclockwise when viewed from +Y toward the origin. χ is the amount of rotation about the directed axis, with +χ clockwise when viewed from the axis toward the origin.

Refinement of the NC symmetry operator parameters is achieved by least-squares minimization of the squared difference in electron density for all NC-symmetry-related points. Thus, one minimizes with respect to the operator parameters, where the sum is taken over all points within the appropriate averaging envelope(s). One starts refinement with low-resolution data (∼6 Å) on a coarse (∼2 Å) grid and monitors progress by following the correlation coefficient between the related electron-density values. Once convergence is obtained, the calculation is resumed with higher-resolution data on a finer grid. Typically, a correlation coefficient of around 0.4 or higher (for a 3 Å MIR map, 1 Å grid) indicates that the operator has been correctly located. The operator refinement is confined to submaps and is facilitated by use of an orthogonal grid. A submap containing the molecules to be averaged is obtained from the programs MAPVIEW or EXTRMAP and can be converted to an orthogonal grid, if needed, by the program MAPORTH, as described earlier.

Masks encompassing the region(s) to be averaged are usually created by hand in the interactive program MAPVIEW. Here, one reads in a submap comprising the desired region of whatever type of map is available, usually an MIR map. An appropriate contour level and initial section are selected and the contoured electron density for that section appears on the screen. One then selects the `add next section' menu item two or three times to create a projection over several sections of the map, since in the projection the molecular boundary is usually more obvious. Selecting the `trace mask' menu item then allows the user to hand-contour the molecular envelope by directing the cursor tied to a mouse or other pointing device. One then moves to an adjacent section and repeats the process until the complete 3D mask is obtained. To simplify matters and speed up the process, there are `copy next mask' and `copy previous mask' menu items allowing one to take advantage of the fact that the mask is a slowly changing function, particularly when near the centre of the molecule. One can use this feature to copy a mask from the previous or following section and apply it to the current section. Up to twelve distinct 3D masks can be selected. Each mask is colour-coded and can be simultaneously displayed superimposed on the contoured electron-density section. Once the mask is completed, the `make asu' menu item is selected to apply crystallographic symmetry operations to all points within the generated envelope masks. If these operations generate a point also within the envelope masks, the point is flagged in red to indicate that it is redundant, indicating that when tracing the mask, one inadvertently strayed into a symmetry-related molecule. After this check for redundancy, all points within the submap distinct from but related to points within the molecular envelopes by crystal symmetry are flagged in green. This enables one to detect packing contacts and also to ensure that all significant electron density has been assigned to some envelope. Upon completion, the mask is written to a file suitable for use either in averaging, solvent flattening (after expansion by BLDCEL), or operator refinement. In cases of `proper' NC symmetry, it is often desirable to trace the averaging envelope mask in a `skewed' map, such that one is looking directly down the NC rotation axis. In this case, it is usually very obvious where the NC symmetry breaks down, simplifying identification of the averaging envelope. If the averaging mask is created in a skewed submap, then the batch program TRNMSK can be used to transform it so as to correspond to the original, unskewed submap for use in averaging calculations (which do not require skewing).

All averaging calculations are carried out by the program MAPAVG (batch), which requires the submap to be averaged along with the envelope masks and NC symmetry operators. A copy of the input submap is made and each grid point in the mask is examined in turn. If the grid point lies within any averaging envelope, then all points related to it by NC symmetry are generated from the operators and examined. If the generated points also lie within the appropriate envelope mask, the electron density there is interpolated, as described earlier, and the density values for all related points are summed. The average value of the electron density is then inserted at the original point in the submap copy. Upon completion, the averaged version of the submap is written to a file and correlation coefficients for regions related by the various NC symmetry operations are output. The averaged submap is then passed to the program BLDCEL along with the averaging mask and the original unaveraged FSFOUR map from which the submap was created. For all points within the averaging envelope(s), their electron-density values and those at points related by crystallographic symmetry are inserted into the full-cell map, and it is written to a file. This file then contains the NC symmetry averaged electron density expanded to a full-cell map that obeys space-group symmetry. As an option, the averaging mask can also be expanded in BLDCEL to a full-cell mask, which could then be used for solvent flattening.

During NC symmetry averaging, phase combination and extension is carried out precisely as described during solvent flattening and negative-density truncation. The only difference is that after generation of each electron-density map, the NC symmetry averaging is carried out on the appropriate submap region, which is then expanded back to a full-cell map prior to each solvent-flattening calculation.

A script to carry out a standard solvent-flattening run is provided along with a description of the expected input files, output files and examples. Not surprisingly, this DOALL procedure does it all. Execution of the script will create a map from an input `anchor' set of phases, typically obtained by MIR, SIR, or MAD methods, and will then create a solvent mask from the map after zeroing out density near heavy-atom sites. This solvent mask is used in four cycles of solvent flattening, combining the map-inverted phase information with the anchor phases. A new solvent mask is then generated, starting from a map produced with the phases after the first four cycles. Four cycles of solvent flattening using the second solvent mask are then carried out, restarting from the original map and combining with the anchor phases. These phases are then used to compute a new map from which a third solvent mask is built. The third mask is then used for eight cycles of solvent flattening, again restarting with the original map and combining with the anchor phases. Supplied in the script, but commented out, are instructions to carry out an arbitrary number of additional phase extension cycles, and then an arbitrary number of phase and amplitude extension cycles, all using the third solvent mask. The combined phases and distribution coefficients are written to a file after all cycles with a given mask are completed.

Additional scripts are provided to carry out phase extension and/or NC symmetry averaging iterations. These scripts are executed after completion of a normal solvent-flattening run with the DOALL procedure. With the EXTNDAVG script, an input number of additional solvent flattening and/or phase combination cycles are carried out, and phase (and possibly amplitude) extension may be requested. Initial and final d spacings are input to the program SLOEXT (batch) along with the number of map modification or phase combination iterations per step, where each step represents the extension by one reciprocal-lattice point in each direction if phase extension is to be carried out. The calculations proceed where the DOALL script leaves off, starting with a map made from the final phases and using the third solvent mask. If NC symmetry averaging is to be carried out, after each map calculation the appropriate submap is extracted from it and is passed to MAPAVG along with the averaging mask. The averaged submap is passed to BLDCEL, where it is expanded to a full-cell map, which then is passed to BNDRY for solvent flattening. Map inversion and phase combination then proceed normally (although possibly with phase extension). Note that, in general, separate masks are used for solvent flattening and averaging.

The EXTNDAVG_MC script carries out the same procedures and options as the EXTNDAVG script, except that it is used when carrying out NC symmetry averaging with multiple crystal forms. Starting phase files, anchor phases, solvent masks, averaging masks and control files are provided for each crystal form. For each form, the solvent flattening and phase combination steps are carried out independently with the appropriate data; however, during the averaging step, maps from all crystal forms are involved.

The interactive program PRECESS is provided to display diffraction data in the form of pseudo-precession photographs. One can display any zone or step through all zones, with the corresponding intensities mapped to a colour scheme. If a grey scale is selected, the image looks very much like a properly exposed precession photograph taken with Polaroid film. When the cursor is placed near a reciprocal-lattice point, the Miller indices, intensity, standard deviation and d spacing are displayed, allowing one to quickly confirm or identify space groups and Laue symmetry. If a scaled file is input containing isomorphous-replacement or anomalous-scattering data, one can display the corresponding intensity differences instead of the native intensities and quickly visualize the distribution of differences to help assess isomorphism.

The interactive program MAPVIEW can be used to examine contoured electron-density or Patterson maps, as well as to examine, create or edit solvent or averaging masks. Either full-cell FSFOUR maps or submaps (including skewed submaps) can be used, although only from the former can any arbitrary region be obtained and reordered interactively. The mask creation and editing functions have been described earlier. The program is very useful for Patterson analysis, evaluation of phasing results and to help decide which region is appropriate for isolating a molecule for use in model building. It is usually crucial for construction of averaging masks, but is also useful for examining or editing other masks.

While MAPVIEW is extremely useful, there are times when it is desirable to have individual plots available either for comparison, stereo viewing of electron density, or incorporation into documents. The program CTOUR (batch) handles these functions and accepts an input FSFOUR map or submap. The CTOUR program can create any number of plot files in a single run, with each consisting of either an individual section, a mono projection, or a stereo projection, with each projection over different multiple sections. If full-cell FSFOUR maps are input, any desired region may be selected, whereas if submaps (including skewed maps) are input, the accessible regions are limited by those present in the input map.

The plot files created by CTOUR are generic in nature and are not directly displayable. One needs a driver program to convert the generic files to the format appropriate for the desired display device. The appropriate drivers for several popular display devices are provided within the package and are described below.

Within the package, fractional atomic coordinates are used extensively, and the program is provided to convert from PDB (Protein Data Bank) to PHASES coordinate files and vice versa. The program prompts for input and output file names, the direction of the conversion, chain or residue ranges, and whether to reset occupancies and/or thermal factors to specified values. The coordinate ranges (both fractional and in PDB coordinates) spanned by the model are also listed.

The program O_TO_SP is provided to convert NC symmetry operators expressed in terms of a 3 × 3 rotation matrix and 1 × 3 translation vector to the PHASES-style spherical polar system described earlier. Although originally written to convert the transformation operator as defined in the O program (Jones et al., 1991), the procedure works for any rotation or translation operator expressed in this form, provided that the operator is applicable to Cartesian coordinates in Å orthogonalized as in the Protein Data Bank (Bernstein et al., 1977).

For efficiency, structure-factor files used within the package are binary; however, the program RD31 is provided to read these binary files and convert them to formatted files that can be examined and possibly edited by the user. The indices, amplitudes, phases, figures of merit, phase probability distribution coefficients, and markers indicating which reflections are centric along with the allowed phase values are thus made readily accessible. A corresponding program, MK31B, is also provided to reverse the process; it reads the formatted (and possibly edited) versions of the structure-factor files and generates the appropriate binary-file equivalents. Additionally, the program XPL_PHI is supplied to convert the binary structure-factor files to a form readable by the X-PLOR program (Brünger et al., 1987) in order to facilitate complete model refinement. Phase and figure-of-merit information are also passed to the output file, allowing refinement with phase restraints if desired.

The program IMPORT allows users to `import' phase information obtained from programs external to the package so it can be used for subsequent calculations within the package. For example, one can use phase and probability distribution information obtained elsewhere to initiate solvent flattening, negative-density truncation and/or NC symmetry averaging within PHASES, or simply to generate and display maps with MAPVIEW or the other graphics programs. Reflection indices, the observed structure-factor amplitude, figure of merit, phase and phase probability distribution coefficients must be supplied, although free format can be used.

The program PSTATS compares phases in two different structure-factor files. It lists mean phase differences as a function of d spacing for common reflections. The program is very useful for comparing results from different phasing strategies and for testing new procedures against error-free phases. It can also be used to check for convergence in iterative procedures or to assess the relative contributions of phase sets during phase combination.

The following density-modification modes (specified by the MODE keyword) are provided by DM:

The most commonly used modes are solvent flattening and histogram matching – these give a good first map in most cases. Recently, multi-resolution modification has been added to this list. Averaging is applied whenever possible. Skeletonization and Sayre's equation are generally only applied in special situations.

Density-modification calculations are somewhat prone to producing grossly overestimated figures of merit (Cowtan & Main, 1996). Users should be aware of this. In general the phases and figures of merit produced by density-modification calculations should only be used for the calculation of weighted maps. They should not be used for the calculation of difference maps or used in refinement or other calculations (the REFMAC program is an exception, containing a mechanism to deal with this form of bias). The use of -type maps should be avoided when the calculated phases are from density modification, since they are dependent on two assumptions, neither of which hold for density modification: that the current phases are very close to being correct and that the calculated amplitudes may only approach the observed values as the phase error approaches zero.

To limit the problems of overestimation, three phase-combination modes are provided (controlled by the COMBINE keyword):

In the case of a molecular-replacement calculation or high noncrystallographic symmetry, it may be desirable only to weight the modified phases and not to recombine them back with the initial phases so that any initial bias may be overcome. In the case of high noncrystallographic symmetry, it may also be possible to restore missing reflections in both amplitude and phase. Options are available for both these situations.

When performing phase extension, the order in which the structure factors are included will affect the final accuracy of the extended phases. The phases obtained from previous cycles of phase extension will be included in the calculation of new phases for the unphased structure factors in the next cycle. A reflection with more accurately determined phases might enhance the phase-extension power of the original set of reflections, whereas a reflection with less accurately determined phases might corrupt the phase-extension power of the original set of reflections and make the phase extension deteriorate quickly. The factors that might affect phase extension are the structure-factor amplitudes, the resolution shell and the figure of merit. Based on the above considerations, the following phase-extension schemes are provided in DM:

All forms of map modification are affected by the overall temperature factor of the data, and histogram matching in particular is critically dependent on the accurate determination of the scale factor. Wilson statistics have been found inadequate for scaling in this case, especially when the data resolution is worse than 3 Å, because of the dip in scattering below 5 Å.

More accurate estimates of the scale and temperature factors may be achieved by fitting the data to a semi-empirical scattering curve (Cowtan & Main, 1998). This curve is prepared using Parseval's theorem, which relates the sum of the intensities to the variance of the map: Thus, the sum of the intensities in a particular resolution shell is proportional to the difference in variance of maps calculated with and without that shell of data. The empirical curve is therefore calculated from the variance in the protein regions of a group of known structures, calculated as a function of resolution. The curve is scaled to the protein volume of the current structure, and a correction is made for the solvent, which is assumed to be flat.

The overall temperature factor is removed, and an absolute scale is imposed by fitting the data to this curve. The use of sharpened F's (with no overall temperature factor) is necessary for histogram matching and often increases the power of averaging for phase extension.

Since the solvent content is used in scaling the data, it is important that this value be entered correctly. However, the volume of the solvent mask may be varied independently of the true solvent content, as discussed in Section 25.2.2.3.

If the user does not supply a solvent mask, the solvent mask is calculated by Wang's (1985) method, using the reciprocal-space approach of Leslie (1987). A number of variants on this algorithm are implemented; however, the parameter that affects the quality of the solvent mask most dramatically is the radius of the smoothing function (Chapter 15.1 ). This parameter may be estimated empirically by where is the resolution limit of the observed amplitudes, and is the mean figure of merit over the same reflections (with w = 0 for unphased reflections).

Once the smoothed map has been determined, cutoff values are chosen to divide the map into protein and solvent regions. If the protein boundary is poorly defined, the user may specify protein, solvent and excluded volumes, in which case two cutoffs are specified and the intermediate region is marked as neither protein nor solvent.

If the user does not supply an averaging mask, it is determined by a local correlation method (Vellieux et al., 1995). A large region covering 27 unit cells is selected, and the local correlation between the maps before and after rotation by one of the noncrystallographic symmetry operators is calculated. The largest contiguous region that is in agreement among different NCS operators is isolated from the local correlation map, and a finer local correlation map is calculated over this volume. This process is iterated until a good mask with a detailed boundary is found.

This approach is fully automatic, except in the case where a noncrystallographic symmetry operator intersects a crystallographic symmetry operator, in which case the mask is not uniquely defined, and some user intervention may be required. The method is robust, and by increasing the radius of the sphere within which the local correlation is calculated, it may be used with very poor maps (Cowtan & Main, 1998). The method is easily extended to include information from multiple averaging operators.

For simplicity of coding, all Fourier transforms are performed in core using real-to-Hermitian and Hermitian-to-real fast Fourier transforms (FFTs). The data are expanded to space group P1 before calculating a map and averaged back to a reciprocal asymmetric unit after inverse transformation. Most of the map modifications preserve crystallographic symmetry, so restricted phases are not constrained except during phase combination.

The target histograms are calculated from the protein regions of several stationary-atom structures at resolutions from 6 to 1.5 Å, according to the method described by Zhang & Main (1990a). The histogram variances should be consistent with the map variances used in scaling the data. The resolution of the target histogram can be accurately matched to the data resolution by averaging the target histograms on either side of the current resolution.

Averaging is performed using a single-step approach (Rossmann et al., 1992), in which every copy of the molecule in a `virtual' asymmetric unit is averaged with every other copy. Density values are obtained at non-grid positions using a 27-point quadratic spectral spline interpolation. A sharpened map is first calculated by dividing by the Fourier transform of the quadratic spline function. The same spline function is then convoluted with the sharpened map to obtain the density value at an arbitrary coordinate (Cowtan & Main, 1998). This approach gives very accurate interpolation from a coarse grid map with relatively little computation and additionally provides gradient information for the refinement of averaging operators.

The multi-crystal averaging calculation in DMMULTI is equivalent to several single-crystal averaging calculations running simultaneously, with the exception that during the averaging step, the molecule density is averaged across every copy in every crystal form. This average is weighted by the mean figure of merit of each crystal form; this allows the inclusion of unphased crystal forms, since in the first cycle they will have zero weight and therefore not disrupt the phasing that is already present. In subsequent cycles, the unphased form contains phase information from the back-transformed density.

This technique can be extremely useful, since adding a new crystal form usually provides considerably more phase information than adding a new derivative if the cross-rotation and translation functions can be solved.

In the multi-crystal case, averaging is performed using a two-step approach, first building an averaged molecule from all the copies in all crystal forms, then replacing the density in each crystal form with the averaged values. This approach is computationally more efficient when there are many copies of the molecule.

The conceptual flow chart of simultaneous density-modification calculations across multiple crystal forms is shown in Fig. 25.2.2.3(a); in practice, this scheme is implemented using a single process and looping over every crystal form at each stage (Fig. 25.2.2.3b). Maps are reconstructed from a large data object containing all the reflection data in every crystal form. Averaging is performed using a second data object containing maps of each averaging domain. By this means, an arbitrary number of domains may be averaged across an arbitrary number of crystal forms.

Multi-crystal averaging has been particularly successful in solving structures from very weak initial phasing, since the data redundancy is usually higher than for single-crystal problems.

The user should not be burdened with the choice of input parameters that they may not be qualified to choose. They also should not be forced to construct an input file that is obscure and difficult to understand. It is hard now to remember what most computer programs were like in the 1970s. Usually, the input to the program was a block of numbers and flags where the meaning of each item was defined by its line and column numbers. This block not only contained information the programmer could never anticipate, like the cell constants, but defined how the computer's memory should be allocated and obscure parameters that could only be estimated after careful reading of research papers.

TNT was one of the first programs in crystallography to have its input introduced with keywords and to allow input statements to come in any order. As an example of the difference, consider the resolution limits. Usually, a crystallographic program would have a line in its input similar to One had to recognize this line amongst many as the line containing the resolution limits. (In many programs, a value of 99 was used to indicate that no lower-resolution limit was to be applied.) In TNT the same data would be entered as The keyword identifies the data as the resolution limit(s). If the statement contained two numbers, they were considered the upper and lower limits of the diffraction data.

The preceding example also shows how default values can be implemented by a program much more safely with keyword-based input. In the previous scheme, if a value was ever to be changed by the user, its place had to be allocated in the input block. This often left numbers floating in the block which were almost never changed, and because they were so infrequently referred to, they were usually unrecognized by the user. It was quite possible for one of these numbers to be accidentally changed and the error unnoticed for quite some time. When the data are introduced with keywords, a data item is not mentioned if the default value is suitable.

The most time-consuming calculations in refinement are the calculation of structure factors from atomic coordinates and the calculation of derivatives of the part of the residual dependent upon the diffraction data with respect to the atomic parameters. The quickest means of performing these calculations requires the use of space-group-optimized fast Fourier transforms (FFTs). The initial implementation of TNT used FFTs to calculate structure factors, but the much slower direct summation method to calculate the derivatives. Within a few years, Agarwal's method (Agarwal, 1978; Agarwal et al., 1981) was incorporated into TNT and from then on all crystallographic calculations were performed with FFTs.

The FFT programs of Ten Eyck (1973, 1977) made very efficient use of computer memory. Another means of saving memory was to recognize that the code for calculating stereochemical restraints did not need to be in the memory when the crystallographic calculations were being performed and vice versa. There were two ways to save memory using this information. One could create a series of `overlays' or one could break the calculation into a series of separate programs. The means for defining an overlay structure were never standardized and could not be ported from one type of computer to another and were, therefore, never attempted in TNT. For this reason, and a number of others mentioned here, TNT is not a single program but a collection of programs, each with a well defined and specialized purpose.

The need to state this goal seems remarkable in these modern times, but the truth is that most computer programs in the 1970s required specific customizations before they could be used. The simplest modifications were the definitions of the maximum number of atoms, residues, atom types etc. accepted by the program. These modifications are still required in Fortran77 programs because that language does not allow the dynamic allocation of memory. However, in most programs today the limits are set high enough that the standard configuration does not present a problem.

The most difficult modification required for programs like PROLSQ was to adapt the calculations to the space group in hand. Their authors usually included code for the space groups they were particularly interested in, leaving all others to be implemented by the user. Writing code for a new space group was often a daunting task for someone who was not an expert programmer and had no tools for testing the modifications.

It is too burdensome to require the user to understand sufficiently the internal workings of a complex calculation that they can code and debug central subroutines of a refinement program. In its initial implementation, TNT avoided this problem, to an extent, by performing the space-group-specific calculations in separate programs. At least the user did not need to modify an existing program. All that was required was the construction of a program that read the proper format file, performed the calculation and wrote its answer in the proper format. The user was required to supply both a program that could calculate structure factors from the model and another program that could calculate the derivative of the diffraction component of the residual function with respect to the atomic parameters of the model.

While a structure-factor program could usually be located, either by finding an existing program or by expanding the model to a lower-symmetry space group for which a program did exist, the requirement of creating a derivative program proved too great a burden. The derivation of the space-group-specific calculation, its implementation and debugging proved too difficult for almost everyone, and this design was quickly abandoned. Instead, an implementation of Agarwal's (1978) algorithm was created. In this method, the derivatives are calculated with a series of convolutions with an map. The calculation of the map is the only space-group-specific part of the calculation, and this was done with a separate program for calculating Fourier syntheses. Such programs were as easy to come by as structure-factor calculation programs and could be replaced by a lower-symmetry program if required.

While it is easier to find or write a program that only calculates a Fourier transform and much easier to debug one than to debug a modification to a larger and more complex program, it is still difficult. The lack of availability of programs for the space group of a crystal often prevented the use of TNT. Over time, programs for more space groups were written and distributed with TNT. Eventually, a method was developed by one of TNT's authors in which FFTs could be calculated using a single program as efficiently as the original space-group-specific programs. Once this program existed, there was no longer the need for isolated structure-factor and Fourier synthesis programs. These calculations have disappeared into the heart of TNT, and TNT consists of many fewer programs today than in the past.

Without a search for symmetry-related bad contacts, it was quite common to build atoms into the same density from two different sides of the molecule. A number of models in the PDB contain these types of errors because neither the refinement nor the graphics programs available at that time would indicate this type of error.

Prior to TNT, one often started a refinement with rigid-body refinement using CORELS and then switched to another program. TNT was the first refinement package to allow both styles of refinement. One was not required to learn about two different packages when running TNT.

This innovation allowed TNT to run efficiently in all space groups available to macromolecular crystals.

With a simple and quick model of the scattering of the bulk solvent in the crystal (Tronrud, 1997), the low-resolution data could be used in refinement for the first time. The inclusion of these data in the calculation of maps greatly improved their appearance.

This method of minimization (Axelsson & Barker, 1984; Tronrud, 1992) allows the direct inclusion of the diagonal elements of the second-derivative matrix and the indirect inclusion of its off-diagonal elements. An additional benefit is that it allows both positional parameters and B factors to be optimized in each cycle. Previously, one was required to hold one class of parameter fixed while the other was optimized. It is much more efficient and simpler for the user to optimize all parameters at once. This method, because it incorporates the diagonal elements directly, produces sets of B factors that agree with the diffraction data better than those from the simple conjugate-gradient method.

It is not uncommon for crystallization enhancers to be found on a special position in the crystal. In addition, cross-linking the molecules in a crystal is often done for various reasons. In both cases, the model contains chemical bonds to a molecule or atoms in another asymmetric unit of the crystal. In order for the stereochemistry of these links to be properly restrained, it must be possible to describe such a link to the refinement program.

When the resolution of the diffraction data set is less than about 2 Å, the individual B factors of a refined model are observed to vary wildly from atom to atom, even when the atoms are bonded to one another. This pattern is not reasonable if one interprets the B factor as a measure of the vibrational motion of the atom. Traditionally, one applies an additional restraint on the B factors of the model, where the ideal value for the difference in B factor for two bonded atoms is zero.

Since it is clear from examinations of higher-resolution models that the B factors generally increase from one side of a bond to the other (e.g. moving from the main chain to the end of a side chain), the traditional restraint is flawed. A restraint library was generated (Tronrud, 1996) where each bond in a residue is assigned a preferred increment in B factor and a confidence (standard deviation) in that increment.

With this enhancement, TNT's minimizer treats the second-derivative matrix as a collection of 5 × 5 element blocks along its diagonal, one block for each atom. While this method improves the rate of convergence for noncrystallographic symmetry restraints, its most significant feature is that it allows the refinement of atoms located on special positions without special handling by the user.

It is rather common in crystals containing multiple copies of a molecule in the asymmetric unit for one or more molecules to have a higher B factor than the others. If the transformation that generates each copy of the molecule consists only of a rotation and translation of the positions of the atoms, the difference in B factors cannot be modelled. The transformations used in TNT now consist of a rotation, translation, a B-factor shift and an occupancy shift.

At Florida State University, Chapman has implemented a real-space refinement package, principally intended for the refinement of virus models, using TNT. He was able to use TNT's minimizer and stereochemical restraints unchanged along with programs he developed to implement his method. More information about this package can be found at http://www.sb.fsu.edu/∼rsref .

Bricogne & Irwin (1996) have developed a maximum-likelihood refinement package using TNT. Not only are TNT's minimizer and stereochemical restraints used, but many of the calculations of the maximum-likelihood residual's derivatives are performed using TNT programs. While Bricogne and co-workers have not needed to modify TNT programs to implement their ideas, there is ongoing collaboration between them and TNT's authors on the development of commands that allow access to some previously internal calculations. More information about Buster can be found at http://babinet.globalphasing.com/buster/ .

When Navraj Pannu wanted to implement Read's maximum-likelihood refinement functions (Pannu & Read, 1996b) in TNT, he choose not to implement it as a separate program, but modified TNT's source code to create a new version of the program Rfactor, named Maxfactor.

Abrahams (1996) conceived the idea that because some amino-acid side chains can be expected to violate the noncrystallographic symmetry (NCS) of the crystal more than others, one could develop a library of the relative strength with which each atom of each residue type would be held by the NCS restraint. He chose to determine these strengths from the average of the current agreement to the NCS of all residues of the same type. For example, if the lysine side chains do not agree well with their NCS mates, the NCS will be loosely enforced for those side chains. On the other hand, if almost all the valine side chains agree well with their mates, then the NCS will be strongly enforced for the few that do not agree well.

He chose to implement this idea by modifying the source code for the TNT program NCS. Since the calculations involved in implementing this idea are simple, the extent of the modifications were not large.

The quickest way to change the position of an atom substantially is not to move it, but rather involves a two-step procedure to remove it from its current (probably wrong) site and to add a new atom at a new (hopefully right) position. Such updating of the model does not imply that all rejected atoms are immediately repositioned in a new site, so the number of atoms to be added does not have to be equal to the number rejected.

Atom rejection in ARP is primarily based on the interpolated or electron density at its atomic centre and the agreement of the atomic density distribution with a target shape. Applied together, these criteria offer powerful means of identifying incorrectly placed atoms, but can suggest false positives. However, a correctly located atom that happens to be rejected should be selected again and put back in the model. Developments of further, perhaps more elegant, criteria may be expected in the future development of the technique.

Atom addition uses the difference or Fourier synthesis. The selection is based on grid points rather than peaks, as the latter are often poorly defined and may overlap with neighbouring peaks or existing atoms, especially if the resolution and phases are poor. The map grid point with the highest electron density satisfying the defined distance constraints is selected as a new atom, grid points within a defined radius around this atom are rejected and the next highest grid point is selected. This is iterated until the desired number of new atoms is found and reciprocal-space minimization is used to optimize the new atomic parameters.

Real-space refinement based on density shape analysis around an atom can be used for the definition of the optimum atomic position. Atoms are moved to the centre of the peak using a target function that differs from that employed in reciprocal-space minimization. The function used is the sphericity of the site, which keeps an atom in the centre of the density cloud but has little influence on the R factor and phase quality. It is only applicable for well separated atoms and is mainly used for solvent atoms at high resolution.

Geometrical constraints are based on a priori chemical knowledge of the distances between covalently linked carbon, nitrogen and oxygen atoms (1.2 to 1.6 Å) and hydrogen-bonded atoms (2.2 to 3.3 Å). Such constraints are applied in rejection and addition of atoms.

The main problem in automatically reconstructing a protein model from electron-density maps is in achieving an initial tracing of the polypeptide chain, even if the result is only partially complete. Subsequent building of side chains and filling of possible gaps is a relatively straightforward task. The complexity of the autotracing can be nicely illustrated as the well known travelling-salesman problem. Suppose one is faced with 100 trial peptide units possessing two incoming and two outgoing connections on average, which is close to what happens in a typical ARP refinement of a 10 kDa protein. Assuming that one of the chain ends is known and that it is possible to connect all the points regardless of the chosen route, then one is faced with the problem of choosing the best chain out of 2⁹⁸. In practice, the situation is even more complex, as not all trial peptides are necessarily correctly identified in the first iteration and some may be missing – analogous to the correctness or incorrectness of the atomic positions described above.

If the connections can be assigned a probability of the peptide being correct, then only the path that visits each node exactly once and maximizes the total probability remains to be identified. Automatic density-map interpretation is based on the location of the atoms in the current model and consists of several steps. Firstly, each atom of the free-atom model is assigned a probability of being correct. Secondly, these weighted atoms are used for identification of patterns typical for a protein. The method utilizes the fact that all residues that comprise a protein, with the exception of cis peptides, have chemically identical main-chain fragments which are close to planar: the structurally identical Cα—C—O—N—Cα trans peptide units.

The problem of searching for possible peptide units and their connections thus becomes straightforward. The most crucial factor is that proteins are composed of linear non-branching polypeptide chains, allowing sets of connected peptides to be obtained from an initial list of all possible tracings. Choosing the direction of a chain path is carried out on the basis of the electron density and observed backbone conformations. The set of peptide units and the list of how they are interconnected do not, however, allow unambiguous tracing of a full-length chain in most cases.

Taken together, the probabilistic identification of the peptide units, the naturally high conformational flexibility of the connections of the peptide units and the limited quality of the X-ray data and/or phases introduce large enough errors to cause density breaks in the middle of the chains or result in density overlaps. Thus, the result of such a tracing is usually a set of several main-chain fragments. The less accurate the starting maps (i.e. initial phases) and the lower the resolution and quality of the X-ray data, the more breaks there will be in the tracing and the greater the number of peptide units which will be difficult to identify.

Residues are differentiated only as glycine, alanine, serine and valine, and complete side chains are not built at this stage. For every polypeptide fragment, a side-chain type can be assigned with a defined probability, using connectivity criteria from the free-atom models and the α-carbon positions of the main-chain fragments. Given these guesses for the side chains and provided the sequence is known, the next step employs docking of the polypeptide fragments into the sequence. Each possible docking position is assigned a score, which allows automated inspection of the side-chain densities, search for expected patterns and building of the most probable side-chain conformations.

An electron-density map can be used to create a free-atom atomic model, with equal atoms placed in regions of high density (Perrakis et al., 1997). To build this model, only the molecular weight of the protein is required, without any sequence information. In brief, a map covering a crystallographic asymmetric unit on a fine grid of about 0.25 Å is constructed. The model is slowly expanded from a random seed by the stepwise addition of atoms in significant electron density and at bonding distances from existing atoms. All atoms in this model and in all subsequent steps are considered to be of the same type. As ARP proceeds, the geometrical criteria remain the same, but the density threshold is gradually reduced, allowing positioning of atoms in lower-density areas of the map. The procedure continues until the number of atoms is about three times that expected. This number is then reduced to about n + 20% atoms by removing atoms in weak density. This method of map parameterization has the advantage that it puts atoms at protein-like distances while covering the whole volume of the protein.

A free-atom model can describe almost every feature of an electron-density map, but this interpretation rarely resembles a conventional conception of a protein. Nevertheless, information from parts of the improved map and the free-atom model can be automatically recognized as containing elements of protein structure by applying the algorithms briefly described for model reconstruction, and at least a partial atomic protein model can be built. Combination of this partial protein model with a free-atom set (a hybrid model) allows a considerably better description of the current map. The protein model provides additional information (in the form of stereochemical restraints), while prominent features in the electron density (unaccounted for by the current model) are described by free atoms.

The procedure of real-space manipulation is coupled to least-squares or maximum-likelihood optimization of the model's parameters against the X-ray data. This is the scheme that we generally refer to as ARP refinement, though there are two distinct modes of ARP: In the unrestrained mode, all atoms in reciprocal-space refinement are treated as free atoms with unknown connectivity and are refined against the experimental data alone. This mode has a higher radius of convergence but needs high-resolution diffraction data to perform effectively. In the restrained mode, a model or a hybrid model is required, i.e. the atoms must belong to groups of known stereochemistry. This stereochemical information, in the form of restraints, can then be utilized during the reciprocal-space minimization, allowing it to proceed with less data, presuming that the connectivity of the input atoms is basically correct.

The hybrid models described above are used as the main tool for obtaining as full a protein model as possible from the map calculated with the initial phases.

Given the information contained in the hybrid model in the traditional form of stereochemical restraints, reciprocal-space refinement can work more efficiently, new improved phases can be obtained and a more accurate and complete protein model can be constructed. The new hybrid model can be re-input to refinement and these steps can be iterated so that improved phases result in construction of ever larger parts of the protein. An almost complete protein model can be obtained in a fully automated way.

Starting from a molecular-replacement solution implies that a search model positioned in the new lattice is already available. The model can be directly incorporated in restrained ARP refinement. If the starting model is very incomplete or different, its atoms can be regarded as free atoms and the solution can be treated as starting from just initial phases. This increases the radius of convergence and minimizes the bias introduced by the search model.

Slightly varying the protocol described for generating models from maps results in a set of slightly different free-atom models. Each model is then submitted to ARP. In protein crystallography, there are generally insufficient data for convergence of free-atom refinement to a global minimum and different starting models result in final models with small differences, i.e. containing different errors. Averaging of these models can be utilized to minimize the overall error. The procedure in effect imposes a random noise, small enough to be eliminated during the subsequent averaging, but large enough to overcome at least some of the systematic errors.

Structure factors are calculated for all the refined models and a vector average of the calculated structure factors is derived. The phase of the vector average is more accurate than that from any of the individual models. A weight, W_wARP, is assigned to each structure factor on the basis of the variance of the two-dimensional distribution of the individual structure factors around the mean. The mean value of W_wARP over all reflections and the R factor after averaging can be used to judge the progress of the averaging procedure.

If the coordinates of one or a few heavy atoms are known, initial phases can be calculated. The problem of solving the structure of such a metalloprotein from the sites of the metal alone can be considered in the same framework as for heavy-atom-replacement solutions. Maps calculated from the phases of heavy atoms alone often have the best defined features within a defined radius of the heavy atom(s). Thus protocols that do not place all atoms at the start but instead perform a slow building while extending the model in a growing sphere around the heavy atom are preferred. When such a model is essentially complete, it can be used for automated tracing and completion of the model.

In this application, the protein (or nucleic acid) model is not rebuilt during refinement, and only the solvent structure is continuously updated, allowing the construction of a solvent model without iterative manual map inspection.

The viewpoint in the MolScript right-handed coordinate system is always located on the positive z axis, looking towards the origin, with the positive x axis to the right. The user obtains the desired view of the structure by specifying rotations and translations of the atomic coordinates; it is not possible to change the location or the direction of the viewpoint.

There are two main benefits with this scheme: The first is that it is similar to the way we handle objects in everyday life: we do not normally fly around the object, but rather move it about with our hands. The other benefit is that together with the coordinate copy feature it can be used to compose an image containing several geometrically related subunits.

The disadvantage is that the atomic coordinates must be transformed before the creation of the graphics objects. This may complicate the composition of an image where another structure or geometric object is to be included. For example, if two separate structures have been aligned structurally by some external procedure, then the user must take care not to destroy the alignment in the process of setting the viewing transformation.

The graphics state consists of the parameters that determine the exact visual appearance of the graphics objects. The default values of the graphics state parameters are reasonable, so that an image of acceptable quality can be produced quickly. However, to obtain high-quality images which emphasize the relevant structural features, the user must usually fine-tune the rendering by modifying the graphics state parameters appropriately for the various graphics objects.

The graphics state may be modified by using the `set' command at any point in the script file. The change can have an effect only on graphics commands below that point in the script. When a graphics command is processed, the object is created according to the current values of the graphics state parameters at that exact point in the script file. It is this property of the graphics state that gives the user a very high degree of control over the composition and appearance of the graphics objects.

A new feature in version 2 of MolScript is the ability to set the colour of residues on a residue-by-residue basis in schematic secondary-structure representations. It is also possible to set the colour of atoms and residues according to a linear function of the B factors.

The graphics commands create the graphics objects to be rendered in the final image. The commands need an atom or residue selection (see below) as argument. The visual attributes, and in some cases the dimensions, of the objects are determined by the graphics state parameters.

The graphics objects include the most common ways to represent atoms, such as simple line drawings, ball-and-stick models and CPK (Corey–Pauling–Koltun) spheres approximating the van der Waals radii of the atoms. The graphics objects representing high-level structures are mainly designed for protein structures, and comprise arrows for β-sheet strands, cylinders or helical ribbons for α-helices, and coils for non-repetitive peptide chain structures. The coil object can also be used to represent oligonucleotide backbone structures.

A set of basic atom and residue selections are provided in MolScript for use as arguments to the graphics commands. Arbitrary subsets of atoms or residues can be specified by joining together the basic selections using a form of Boolean operators. Unfortunately, the Boolean expression feature may sometimes be difficult to understand for the non-expert user. One should consider the entire expression as a test to be applied to every atom or residue. Any atom or residue for which the Boolean expression evaluates to `true' will be selected as argument for the command.

Externally defined objects described by points, lines or triangular surfaces may be included in the image. The objects may optionally be transformed by the most recent transformation applied to the coordinate data. This feature allows import of arbitrary geometry created by some external software, e.g. molecular surfaces, electron-density representations or electrostatic field lines. The graphics state parameters apply to the rendering of the external objects in the image.

PostScript (Adobe Systems Inc., 1985) is a page description language for controlling high-quality printers. More information can be found at the Adobe Inc. web site, http://www.adobe.com/print/ .

The PostScript output mode relies on the painter's algorithm for hidden-surface removal. The most distant graphics segments are output first, continuing with the segments closer to the viewpoint, which may obliterate previously rendered segments. The implementation of this procedure is straightforward, and gives good results provided that the graphics objects are subdivided into sufficiently small segments. The PostScript mode allows more than one plot (image) to be rendered on a single page.

The Raster3D suite of programs (Merritt & Bacon, 1997) produces high-quality images using a ray-tracing algorithm. MolScript can produce the input file required for the `render' program, which is the core program of the Raster3D suite. The web site for Raster3D is http://www.bmsc.washington.edu/raster3d/ .

The Raster3D mode features highlighting, transparency and shadows to produce MolScript images of very high visual quality.

The VRML97 standard (Virtual Reality Modeling Language, formerly VRML 2.0) allows storage and transmission of 3D scenes in a system-independent manner over the web. Software to view VRML97 files is typically included in any modern web browser. A web site containing more information on VRML97 is http://www.web3d.org/x3d/specifications/vrml/ .

The VRML97 mode allows hyperlinking of objects. The MolScript implementation is optimized to produce output files that are as small as possible, but the file size is strongly dependent on the value of the `segments' parameter in the graphics state.

OpenGL is a standard API (Applications Programming Interface) for interactive 3D graphics. It is available on most current computer systems. For more information, see the web site http://www.opengl.org/ .

The OpenGL output mode allows a certain degree of interactivity, in contrast to the other output modes. It is possible to initiate execution of the MolScript program in OpenGL mode in one window on the screen, while keeping the script file in a separate text-editor window. The image is rotatable in 3D in the OpenGL window. The script can be edited in its window, and the modified script can be re-read and displayed directly by the MolScript program in its OpenGL window. This simplifies to some extent the iterative fine-tuning of the script.

Raster image files in several different formats can be created by MolScript. Currently these include SGI RGB, encapsulated PostScript (EPS), JPEG, PNG and GIF image formats. The JPEG, PNG and GIF formats require that external software libraries are available during the compilation and linking of the MolScript program. Software libraries for several of these image formats are available on the web; links are given at the official MolScript web site http://www.avatar.se/molscript/ .

The image file formats essentially capture the raster image created by the OpenGL implementation. The EPS format was a variant of the PostScript output mode in version 1 of MolScript, but for various reasons this has changed in version 2 to an encoding of the OpenGL raster image.

XDS is organized into eight steps (major subroutines) which are called in succession by the main program. Information is exchanged between the steps by files (see Table 25.2.9.1), which allows repetition of selected steps with a different set of input parameters without rerunning the whole program. ASCII files can be inspected and modified using a text editor, whereas types DIR and BIN indicate binary random access and unformatted sequential access files, respectively. All files have a fixed name defined by XDS, which makes it mandatory to process each data set in a newly created directory. Clearly, one should not run more than one XDS job at a time in any given directory. Output files affected by rerunning selected steps (see Table 25.2.9.1) should also first be given another name if their original contents are meant to be saved.

Data processing begins by copying an appropriate input file into the new directory. Input-file templates are provided with the XDS package for a number of frequently used data-collection facilities. The copied input file must be renamed XDS.INP and edited to provide the correct parameter values for the actual data-collection experiment. All parameters in XDS.INP are named by keywords containing an equal sign as the last character, and many of them will be mentioned here in context to clarify their meaning. Execution of XDS (JOB= XDS) invokes each of the eight program steps as described below. Results and diagnostics from each step are saved in files with the extension LP attached to the program step name. These files should always be studied carefully to see whether processing was satisfactory or – in case of failure – to find out what could have gone wrong.

XYCORR calculates a lookup table of additive spatial corrections at each detector pixel and stores it in the file XYCORR.TBL. The data images are often already corrected for geometrical distortions, in which case XYCORR produces a table of zeros or – as for spiral read-out imaging plate detectors – computes the small corrections resulting from radial (ROFF=) and tangential (TOFF=) offset errors of the scanner. For some multiwire and CCD detectors that deliver geometrically distorted images, corrections are derived from a calibration image (BRASS_PLATE_IMAGE= file name). This image displays the response to a brass plate containing a regular grid of holes which is mounted in front of the detector and illuminated by an X-ray point source, e.g. ⁵⁵Fe. Clearly, the source must be placed exactly at the location to be occupied by the crystal during the actual data collection, as photons emanating from the calibration source are meant to simulate all possible diffracted beam directions. For visual control using the VIEW program, spots that have been located and accepted from the brass-plate image by XYCORR are marked in the file FRAME.pck.

Problems : (a) A misplaced calibration source leads to an incorrect lookup table, impairing the correct prediction of the observed diffraction pattern in subsequent program steps. (b) Underexposure of the calibration image results in an incomplete and unreliable list of calibration spots.

INIT estimates the initial background at each pixel and determines the trusted region of the detector surface. The total background at each pixel is the sum of the detector noise and the X-ray background. The detector noise, saved in the lookup table BLANK.TBL, is determined from a specific image recorded in the absence of X-rays (DARK_CURRENT_IMAGE=) or is assumed to be a constant derived from the mean recorded value in each corner of the data images. A lookup table of the X-ray background, saved on the file BKGPIX.TBL, is obtained from the first few data images by the following two-pass procedure. To exclude diffraction spots in the data image, the minimum of the five values at (x, y), is used as a lower background estimate at pixel (x, y) in the first pass. In the second pass, the background is taken as the maximum of the lower estimates at these five locations. Ideally, the parameters SPOT_WIDTH_ALONG_X= 2*dx + 1, SPOT_WIDTH_ALONG_Y= 2*dy + 1 are chosen to match the extent of a spot. The lookup table is obtained by adding the X-ray background from each image. Shaded regions on the detector (i.e. from the beam stop), pixels outside a user-defined circular region (RMAX=) or pixels with an undefined spatial correction value are classified as untrustworthy and marked by −3. The table should be inspected using the VIEW program.

Problems : (a) The addition of background from too many data images may exceed 262 144 at some pixels, which are removed from the trusted detector region due to internal number overflow. (b) Some detectors with insufficient protection from electromagnetic pulses may generate badly spoiled images whose inclusion leads to a completely wrong X-ray background table. These images can be identified in INIT.LP by their unexpected high mean pixel contents, and this step should be repeated with a different set of images.

COLSPOT locates, at most, 500 000 strong diffraction spots occurring in a subset of the data images and saves their centroids on the file SPOT.XDS. Up to ten ranges of contiguous images (SPOT_RANGE=) may be specified explicitly; otherwise, spots are taken from the first few data images, covering a total rotation range of 5°. Spots are located automatically by comparing each pixel value with the mean value and standard deviation of surrounding pixels, as described in Chapter 11.3 . A lower threshold for accepting pixels and a minimum required number of such pixels within a spot can be defined in XDS.INP by the parameters MINIMUM_SIGNAL_TO_NOISE_FOR_LOCATING_SPOTS= and MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT=, respectively.

Problem : Sharp edges like ice rings in the images can lead to an excessive number of pixels erroneously classified as contributing to a diffraction spot which extends over many adjacent images, thereby causing a hash-table overflow. The problem can be avoided by specifying non-adjacent images for spot search.

IDXREF uses the initial parameters describing the diffraction experiment as provided by XDS.INP and the observed centroids of the spots occurring in the file SPOT.XDS to find the orientation, metric and symmetry of the crystal lattice, and it refines all or a specified subset of these parameters. On return, the complete set of parameters are saved in the file XPARM.XDS, and the original file SPOT.XDS is replaced by a file of identical name – now with indices attached to each observed spot. Spots not belonging to the crystal lattice are given indices 0, 0, 0. XDS considers the run successful if at least 70% of the given spots can be explained with reasonable accuracy; otherwise, XDS will stop with an error message. Alien spots often arise because of the presence of ice or small satellite crystals, and continuation of data processing may still be meaningful. In this case, XDS is called again with an explicit list of the subsequent steps specified in XDS.INP.

Using and understanding the results reported in IDXREF.LP requires a knowledge of the concepts employed by this step, as described in Chapter 11.3 . First, a reciprocal-lattice vector, referring to the unrotated crystal, is computed from each observed spot centroid. Differences between any two reciprocal-lattice vectors that are above a specified minimal length (SEPMIN=) are accumulated in a three-dimensional histogram. These difference vectors will form clusters in the histogram, since there are many different pairs of reciprocal-lattice vectors of nearly identical vector difference. The clusters are found as maxima in the smoothed histogram (CLUSTER_RADIUS=), and a basis of three linearly independent cluster vectors is selected that allows all other cluster vectors to be expressed as nearly integral multiples of small magnitude with respect to this basis. The basis vectors and the 60 most populated clusters with attached indices are listed in IDXREF.LP. If many of the indices deviate significantly from integral values, the program is unable to find a reasonable lattice basis and all further processing will be meaningless.

If the space group and cell constants are specified, a reduced cell is derived, and the reciprocal-basis vectors found above are reinterpreted accordingly; otherwise, a reduced cell is determined directly from the reciprocal basis. Parameters of the reduced cell, coordinates of the reciprocal-basis vectors and their indices with respect to the reduced cell are reported.

Based on the orientation and metric of the reduced cell now available, IDXREF indexes up to 3000 of the strongest spots by the local-indexing method. This method considers each spot as a node of a tree and identifies the largest subtree of nodes which can be assigned reliable indices. The number of reflections in the ten largest subtrees is reported and usually shows a dominant first tree corresponding to a single lattice, whereas alien spots are found in small subtrees. Reflections in the largest subtree are used for initial refinement of the basis vectors of the reduced cell, the incident beam wave vector and the origin of the detector, which is the point in the detector plane nearest to the crystal. Experience has shown that the detector origin and the direction of the incident beam are often specified with insufficient accuracy, which could easily lead to a misindexing of the reflections by a constant offset. For this reason, IDXREF considers alternative choices for the index origin and reports their likelihood for being correct. Parameters controlling the local indexing are INDEX_ERROR=, INDEX_MAGNITUDE=, INDEX_QUALITY= (corresponding to ɛ, δ, in Chapter 11.3 ) and INDEX_ORIGIN= h₀, k₀, l₀, which is added to the indices of all reflections in the tree. After initial refinement based on the reflections in the largest subtree, all spots that can now be indexed are included. Usually, the detector distance and the direction of the rotation axis are not refined, but if the spots were extracted from images covering a large range of total crystal rotation, better results are obtained by including these parameters in the refinement (REFINE=).

The refined metric parameters of the reduced cell are used for testing each of the 44 possible lattice types, as described in Chapter 11.3 . For each lattice type, IDXREF reports the likelihood of being correct, the conventional cell parameters and the linear transformation relating original indices to the new indices with respect to the conventional cell. However, no automatic decisions for space-group assignment are made by XDS. If the space group and cell constants are provided by the user, the reduced-cell vectors are reinterpreted accordingly; otherwise, data processing continues with the crystal being described by its reduced-cell basis vectors and triclinic symmetry. On completion, when integrated intensities are available, the user chooses any plausible space group according to the rated list of the 44 possible lattice types and repeats only the CORRECT and GLOREF steps with the appropriate conventional cell parameters and reindexing transformation (see below).

Problems : (a) Indices of many difference-vector clusters deviate significantly from integral values. This can be caused by incorrect input parameters, such as rotation axis, oscillation angle or detector position, by a large fraction of alien spots in SPOT.XDS, by placing the detector too close to the crystal, or by inappropriate choice of parameters SEPMIN= and CLUSTER_RADIUS= in densely populated images. (b) Indexing and refinement is unsatisfactory despite well indexed difference-vector clusters. This is probably caused by selection of an incorrect index origin, and IDXREF should be rerun with plausible alternatives for INDEX_ORIGIN= after a visual check of a data image with the VIEW program.

COLPROF extracts the three-dimensional profile of each reflection predicted to occur in the rotation images within the trusted region of the detector surface and saves all profiles on the file XREC.XDS. A scaling factor is determined for each image, derived by comparing its background region (after subtraction of the detector noise) with the current X-ray background table. This table, initially obtained from the file BKGPIX.TBL, is updated by the background from each data image at a rate defined by the input parameter BFRAC=. For visual control, the contents of the updated X-ray background table are saved on the file BKGPIX.IMG at the end of this program step. Information for predicting reflection positions is initially provided by the file XPARM.XDS. These parameters are either kept constant or refined periodically using centroids of the most recently found strong diffraction spots as data reduction proceeds (REFINE=, NUMBER_OF_FRAMES_BETWEEN_REFINEMENT_IN_COLPROF=, NUMBER_OF_REFLECTIONS_USED_FOR_REFINEMENT_IN_COLPROF=, WEAK=).

In order to include all pixels contributing to the intensity of a spot, approximate values describing their extension and form must be specified, as defined in Chapter 11.3 by the parameters δ_M, σ_M, δ_D, σ_D. The value for BEAM_DIVERGENCE= δ_D = arctan(spot diameter/detector distance) is found by measuring the diameter of a strong spot in a data image displayed by the VIEW program and should include a few adjacent background pixels. The form of a spot is roughly described as a Gaussian and its standard deviation is specified by the parameter BEAM_DIVERGENCE_E.S.D.= σ_D, which is usually about one-sixth to a tenth of δ_D. Similarly, REFLECTING_RANGE= δ_M is the approximate rotation angle required for a strong spot recorded perpendicular to the rotation axis to pass completely through the Ewald sphere. The standard deviation of the intensity distribution is given by the mosaicity REFLECTING_RANGE_E.S.D.= σ_M. Thus, a three-dimensional domain of pixels belonging to each reflection is defined by the above parameters, and the program automatically removes pixels contaminated by neighbouring reflections. It determines and subtracts the background, corrects for spatial distortions, and maps each pixel content into a reflection-specific coordinate system centred on the Ewald sphere (see Chapter 11.3 ). The form of these profiles is then similar for all reflections, and their mean obtained from superimposition of strong reflections is reported at regular intervals. On return from this step, the data image last processed with all expected spots encircled is saved in the file FRAME.pck for inspection using the VIEW program.

Problems : (a) Off-centred profiles indicate incorrectly predicted reflection positions by using the parameters provided by the file XPARM.XDS (i.e. misindexing by using a wrong origin of the indices), crystal slippage, or change in the incident beam direction. (b) Profiles extending to the borders of the box indicate too-small values for BEAM_DIVERGENCE= or REFLECTING_RANGE=. This leads to incorrect integrated intensities because of truncated reflection profiles and unreliable background determination. (c) Display of the file FRAME.pck showing spots which are not encircled. If these unexpected reflections are not close to the spindle and are not ice reflections, it is likely that the parameters provided by the file XPARM.XDS are wrong.

PROFIT estimates intensities from the three-dimensional profiles of the reflections stored in the input file XREC.XDS and saves the results in the file PROFIT.HKL. In the first pass, templates are generated by superimposing profiles of fully recorded strong reflections, and all grid points with a value above a minimum percentage of the maximum in the template (CUT=) are defined as elements of the integration domain. To allow for variations of their shape, profile templates are generated from reflections located at nine regions of equal size covering the detector surface and additional sets of nine to cover equally sized batches of images. Standard deviations, REFLECTING_RANGE_E.S.D.= and BEAM_DIVERGENCE_E.S.D.=, observed for each profile template are reported and – in the case of large discrepancies – could be used for rerunning COLPROF with better values for these parameters. In the second pass, intensities and their standard deviations are estimated by fitting the reflection profile to its template, as described in Chapter 11.3 . Overloaded (OVERLOAD=) or incomplete reflections covering less than a minimum percentage of the template volume (MINPK=) or reflections with unreliable background are excluded from further processing.

Problem : The program stops because there are no strong spots for learning profile templates. It is likely that parameters REFLECTING_RANGE=, BEAM_DIVERGENCE= etc., which define the box dimensions, have been incorrectly chosen. After correction, both the COLPROF and PROFIT step should be repeated.

CORRECT applies Lorentz and polarization correction factors as well as factors that partially compensate for radiation damage and absorption effects to intensities and standard deviations of all reflections found in the file PROFIT.HKL, and saves the results on the files XDS.HKL and either NORMAL.HKL or ANOMAL.HKL (if Friedel's law is broken, as specified by a positive value for the input parameter DELFRM=). These factors are determined from many symmetry-equivalent reflections usually found in the data images such that their integrated intensities become as similar as possible. The residual scatter of these intensities is a more realistic measure of their errors and is used to determine a correction factor for the standard deviations previously estimated from profile fitting. An initial guess for this factor (WFAC1=) is provided in XDS.INP and is used to identify outliers, which are collected in the file MISFITS for separate analysis.

Data quality as a function of resolution is described by the agreement of the intensities of symmetry-related reflections and quantified by the R factors R_sym and the more robust indicator R_meas (Diederichs & Karplus, 1997). These R factors as well as the intensities of all reflections with indices of type h00, 0k0 and 00l and those expected to be systematically absent are important indicators for identification of the correct space group. Clearly, large R factors or many rejected reflections (MISFITS) or large observed intensities for systematically absent reflections suggest that the assumed space group or the indexing is incorrect. It is easy to test other possible space groups (SPACE_GROUP_NUMBER=) by simply repeating the CORRECT and GLOREF steps after copying the appropriate reindexing transformation (REIDX=) and conventional cell constants (UNIT_CELL_CONSTANTS=) found in the rated table of the 44 possible lattice types in IDXREF.LP to XDS.INP. One should remember, however, that the final choice to be kept should be run last, as XDS overwrites earlier versions of the output files.

Another useful feature is the possibility of comparing the new data with those from a previously measured crystal (REFERENCE_DATA_SET= file name). For some space groups, like P4₂, possessing an ambiguity in the choice of axes, comparison with the reference data set allows one to identify the consistent solution from the complete set of alternatives already listed in IDXREF.LP together with their required index transformation. Reference data are also quite useful for recognizing misindexing or for testing potential heavy-atom derivatives.

Problems : (a) Incomplete data sets may lead to wrong conclusions about the space group, as some of its symmetry operators might not be involved in the R-factor calculations. (b) Conventional cell parameters, as listed in IDXREF.LP, often violate constraints imposed by the space group and must be edited accordingly after copying to XDS.INP.

GLOREF refines the diffraction parameters by using the observed positions of all strong spots contained in the file PROFIT.HKL. It reports the root-mean-square error between calculated and observed positions along with the refined unit-cell constants. Again, for testing possible space groups, the crystallographer consults the table printed by the IDXREF step and selects the appropriate reindexing transformation and starting values for the conventional cell constants. The refined diffraction parameters (after possible reindexing) are saved on the file GXPARM.XDS, which is identical in format to XPARM.XDS. Replacing XPARM.XDS with the new file offers a convenient way for repeating COLPROF, now with a better set of parameters.

Problem : GLOREF will fail if the crystal slips during data collection.

XPLAN supports the planning of data collection. It is based upon information provided by XPLAN.INP and the input files XPARM.XDS and BKGPIX.TBL, both of which become available on processing a few test images with XDS. XPLAN estimates the completeness of new reflection data, expected to be collected for each given starting angle and total crystal rotation, and reports the results for a number of selected resolution shells in the file XPLAN.LP. To minimize recollection of data, the name of a file containing already measured reflections can be specified in XPLAN.INP.

Problems : (a) Incorrect results may occur for some space groups, e.g. P4₂, if the unit cell determined by XDS from processing a few test images implicates reflection indices inconsistent with those from the already-measured data. The correct cell choice can be found, however, by using the old data as a reference and repeating CORRECT and GLOREF with the appropriate reindexing transformation, followed by copying GXPARM.XDS to XPARM.XDS. The same applies if IDXREF was run for an unknown space group and then reindexed in CORRECT and GLOREF. (b) XPLAN ignores potential reflection overlap due to the finite oscillation range covered by each image.

XSCALE accepts one or more files of type XDS.HKL, NORMAL.HKL or ANOMAL.HKL, as obtained from data processing with XDS or merged output files from DENZO (see Chapter 11.4 ), determines scaling factors for the reflection intensities, merges symmetry-equivalent observations into a unique set, and reports data completeness and quality in the file XSCALE.LP. The desired program action is specified in the file XSCALE.INP. It consists of a definition of shells used for analysing the resolution dependency of data quality and completeness, space-group number and cell constants, and one line for each set of input reflections. Each set has a file name, type identifier, a resolution window for accepting data, a weighting factor for the standard deviation of the reflection intensity, a decision constant for accepting Bijvoet pairs, a number controlling the degree of smoothness of the scaling function and an optional file name for specifying the output file.

Problem : All reflections beyond the highest shell specified for analysing the resolution dependency of data quality and completeness will be ignored regardless of the resolution window given for each data set.

VIEW is used for visualizing data images as well as control images produced by XDS. It responds to navigation commands entered by movements of a mouse, and reports the corresponding image coordinates and their pixel contents upon activation of the mouse buttons. VIEW also allows magnification of selected image portions and changes in colour.

Problem : For many detector image formats as well as XDS produced images, the true pixel value is stored in a coded form which is interpreted by VIEW as a signed integer. Numbers less than −4095 displayed by VIEW correspond to large positive pixel values.

XDSCONV accepts reflection-intensity data files as produced by XSCALE or CORRECT and converts them into a format required by software packages for structure determination. XDSCONV estimates structure-factor moduli based on the assumption that the intensity data set obeys Wilson's distribution and uses a Bayesian approach to statistical inference as described by French & Wilson (1978). For anomalous intensity data, both structure-factor amplitudes and are simultaneously estimated from the Bijvoet intensity pair by a method similar to that described by Lewis & Rees (1983) – which accounts for the correlation between and . The output file generated may inherit test reflections previously used for calculating a free R factor (Brünger, 1992b) or may contain new test reflections selected by XDSCONV.

Space-group-general automatic Patterson map interpretation was introduced in the program SHELXS86 (Sheldrick, 1985); completely different algorithms are employed in the current version of SHELXS, based on the Patterson superposition minimum function (Buerger, 1959, 1964; Richardson & Jacobson, 1987; Sheldrick, 1991, 1998a; Sheldrick et al., 1993). The algorithm used in SHELXS is as follows:

The program SHELXD (Sheldrick & Gould, 1995; Sheldrick, 1997, 1998b) is now part of the SHELX system. It is designed both for the ab initio solution of macromolecular structures from atomic resolution native data alone and for the location of heavy-atom sites from ΔF or F_A values at much lower resolution, in particular for the location of larger numbers of anomalous scatterers from MAD data. The dual-space approach of SHELXD was inspired by the Shake and Bake philosophy of Miller et al. (1993, 1994) but differs in many details, in particular in the extensive use it makes of the Patterson function that proves very effective in applications involving ΔF or F_A data. The ab initio applications of SHELXD have been described in Chapter 16.1 , so only the location of heavy atoms will be described here. An advantage of the Patterson function is that it provides a good noise filter for the ΔF or F_A data: negative regions of the Patterson function can simply be ignored. On the other hand, the direct-methods approach is efficient at handling a large number of sites, whereas the number of Patterson peaks to analyse increases with the square of the number of atoms. Thus, for reasons of efficiency, the Patterson function is employed at two stages in SHELXD: at the beginning to obtain starting atom positions (otherwise random starting atoms would be employed) and at the end, in the form of the triangular table described above, to recognize which atoms are correct. In between, several cycles of real/reciprocal space alternation are employed as in the ab initio structure solution, alternating between tangent refinement, E-map calculation and peak search, and possibly random omit maps, in which a specified fraction of the potential atoms are left out at random.

Since the input files for the direct and Patterson methods in SHELXS and the integrated method in SHELXD are almost identical (usually only one instruction needs to be changed), it is easy to try all three methods for difficult problems. The Patterson map interpretation in SHELXS is a good choice if the heavy atoms have variable occupancies and it is not known how many heavy-atom sites need to be found; the direct-methods approaches work best with equal atoms. In general, the conventional direct methods in SHELXS will tend to perform best in a non-polar space group that does not possess special positions; however, for more than about a dozen sites, only the integrated approach in SHELXD is likely to prove effective; the SHELXD algorithm works best when the number of sites is known. Especially for the MAD method, the quality of the data is decisive; it is essential to collect data with a high redundancy to optimize the signal-to-noise ratio and eliminate outliers. In general, a resolution of 3.5 Å is adequate for the location of heavy-atom sites. At the time of writing, SHELXD does not include facilities for the further calculations necessary to obtain maps. Experience indicates that it is only necessary to refine the B values of the heavy atoms using other programs; their coordinates are already rather precise.

Excellent accounts of the theory of direct and Patterson methods with extensive literature references have been presented in IT B Chapter 2.2 by Giacovazzo (2001) and Chapter 2.3 by Rossmann & Arnold (2001).

In refining macromolecular structures, it is almost always necessary to supplement the diffraction data with chemical information in the form of restraints. A typical restraint is the condition that a bond length should approximate to a target value with a given estimated standard deviation; restraints are treated as extra experimental data items. Even if the crystal diffracts to 1.0 Å, there may well be poorly defined disordered regions for which restraints are essential to obtain a chemically sensible model (the same can be true of small molecules too!). SHELXL is generally not suitable for refinements at resolutions lower than about 2.5 Å because it cannot handle general potential-energy functions, e.g. for torsion angles or hydrogen bonds; if noncrystallographic symmetry restraints can be employed, this limit can be relaxed a little.

For some purposes (e.g. riding hydrogen atoms, rigid-group refinement, or occupancies of atoms in disordered side chains), constraints, exact conditions that lead to a reduction in the number of variable parameters, may be more appropriate than restraints; SHELXL allows such constraints and restraints to be mixed freely. Riding hydrogen atoms are defined such that the C—H vector remains constant in magnitude and direction, but the carbon atom is free to move; the same shifts are applied to both atoms, and both atoms contribute to the least-squares derivative sums. This model may be combined with anti-bumping restraints that involve hydrogen atoms, which helps to avoid unfavourable side-chain conformations. SHELXL also provides, e.g., methyl groups that can rotate about their local threefold axes; the initial torsion angle may be found using a difference-electron-density synthesis calculated around the circle of possible hydrogen-atom positions.

The original SHELX refinement algorithms were modelled closely on those described by Cruickshank (1970). For macromolecular refinement, an alternative to (blocked) full-matrix refinement is provided by the conjugate-gradient solution of the least-squares normal equations as described by Hendrickson & Konnert (1980), including preconditioning of the normal matrix that enables positional and displacement parameters to be refined in the same cycle. The structure-factor derivatives contribute only to the diagonal elements of the normal matrix, but all restraints contribute fully to both the diagonal and non-diagonal elements, although neither the Jacobian nor the normal matrix itself are ever generated by SHELXL. The parameter shifts are modified by comparison with those in the previous cycle to accelerate convergence whilst reducing oscillations. Thus, a larger shift is applied to a parameter when the current shift is similar to the previous shift, and a smaller shift is applied when the current and previous shifts have opposite signs.

SHELXL refines against F² rather than F, which enables all data to be used in the refinement with weights that include contributions from the experimental uncertainties, rather than having to reject F values below a preset threshold; there is a choice of appropriate weighting schemes. Provided that reasonable estimates of σ(F²) are available, this enables more experimental information to be employed in the refinement; it also allows refinement against data from twinned crystals.

Inversion of the full normal matrix (or of large matrix blocks, e.g. for all positional parameters) enables the precision of individual parameters to be estimated (Rollett, 1970), either with or without the inclusion of the restraints in the matrix. The standard uncertainties in dependent quantities (e.g. torsion angles or distances from mean planes) are calculated in SHELXL using the full least-squares correlation matrix. These standard uncertainties reflect the data-to-parameter ratio, i.e. the resolution and completeness of the data and the percentage of solvent, and the quality of the agreement between the observed and calculated F² values (and the agreement of restrained quantities with their target values when restraints are included).

Full-matrix refinement is also useful when domains are refined as rigid groups in the early stages of refinement (e.g. after structure solution by molecular replacement), since the total number of parameters is small and the correlation between parameters may be large.

The motion of macromolecules is clearly anisotropic, but the data-to-parameter ratio rarely permits the refinement of the six independent anisotropic displacement parameters (ADPs) per atom; even for small molecules and data to atomic resolution, the anisotropic refinement of disordered regions requires the use of restraints. SHELXL employs three types of ADP restraint (Sheldrick 1993; Sheldrick & Schneider, 1997). The rigid bond restraint, first suggested by Rollett (1970), assumes that the components of the ADPs of two atoms connected via one (or two) chemical bonds are equal within a specified standard deviation. This has been shown to hold accurately (Hirshfeld, 1976; Trueblood & Dunitz, 1983) for precise structures of small molecules, so it can be applied as a `hard' restraint with small estimated standard deviation. The similar ADP restraint assumes that atoms that are spatially close (but not necessarily bonded, because they may be different components of a disordered group) have similar U^ij components. An approximately isotropic restraint is useful for isolated solvent molecules. These two restraints are only approximate and so should be applied with low weights, i.e. high estimated standard deviations.

The transition from isotropic to anisotropic roughly doubles the number of parameters and almost always results in an appreciable reduction in the R factor. However, this represents an improvement in the model only when it is accompanied by a significant reduction in the free R factor (Brünger, 1992b). Since the free R factor is itself subject to uncertainty because of the small sample used, a drop of at least 1% is needed to justify anisotropic refinement. There should also be a reduction in the goodness of fit, and the resulting displacement ellipsoids should make chemical sense and not be `non-positive-definite'!

When there are several identical chemical moieties in the asymmetric unit, a very effective restraint is to assume that the chemically equivalent 1,2 and 1,3 distances are the same, but unknown. This technique is easy to apply using SHELXL and is often employed for small-molecule structures and, in particular, for oligosaccharides. Similarly, the terminal P—O bond lengths in DNA structures can be assumed to be the same (but without a target value), i.e. it is assumed that the whole crystal is at the same pH. For proteins, the method is less suitable because of the different abundance of the different amino acids, and, in any case, good target distances are available (Engh & Huber, 1991).

Local noncrystallographic symmetry (NCS) restraints (Usón et al., 1999) may be applied to restrain corresponding 1,4 distances and isotropic displacement parameters to be the same when there are several identical macromolecular domains in the asymmetric unit; usually, the 1,2 and 1,3 distances are restrained to standard values in such cases and so do not require NCS restraints. Such local NCS restraints are more flexible than global NCS constraints and – unlike the latter – do not require the specification of a transformation matrix and mask.

There are many ways of modelling disorder using SHELXL, but for macromolecules the most convenient is to retain the same atom and residue names for the two or more components and assign a different `part number' (analogous to the PDB alternative site flag) to each component. With this technique, no change is required to the input restraints etc. Atoms in the same component will normally have a common occupancy that is assigned to a `free variable'. If there are only two components, the sum of their occupancies can be constrained to be unity; if there are more than two components, the sum of their free variables may be restrained to be unity. Since any linear restraint may be applied to the free variables, they are very flexible, e.g. for modelling complicated disorder. By restraining distances to be equal to a free variable, a standard deviation of the mean distance may be calculated rigorously using full-matrix least-squares algebra.

Babinet's principle is used to define a bulk solvent model with two refinable parameters (Moews & Kretsinger, 1975), and global anisotropic scaling (Usón et al., 1999) may be applied using a parameterization proposed by Parkin et al. (1995). An auxiliary program, SHELXWAT, allows automatic water divining by iterative least-squares refinement, rejection of waters with high displacement parameters, difference-electron-density calculation, and a peak search for potential water molecules that make at least one good hydrogen bond and no bad contacts; this is a simplified version of the ARP procedure of Lamzin & Wilson (1993).

SHELXL provides facilities for refining against data from merohedral, pseudo-merohedral and non-merohedral twins (Herbst-Irmer & Sheldrick, 1998). Refinement against data from merohedrally twinned crystals is particularly straightforward, requiring only the twin law (a 3 × 3 matrix) and starting values for the volume fractions of the twin components. Failure to recognize such twinning not only results in high R factors and poor quality maps, it can also lead to incorrect biochemical conclusions (Luecke et al., 1998). Twinning can often be detected by statistical tests (Yeates & Fam, 1999), and it is probably much more widespread in macromolecular crystals than is generally appreciated!

Least-squares refinement as implemented in SHELXL and other programs is appropriate for structural models that are relatively complete, but when an appreciable fraction of the structure is still to be located, maximum-likelihood refinement (Bricogne, 1991; Pannu & Read, 1996a; Murshudov et al., 1997) is likely to be more effective, especially when experimental phase information can be incorporated (Pannu et al., 1998). Within the least-squares framework, there are still several possible ways of improving the radius of convergence. SHELXL provides the option of gradually extending the resolution of the data during the refinement; a similar effect may be achieved by a resolution-dependent weighting scheme (Terwilliger & Berendzen, 1996). Unimodal restraints, such as target distances, are less likely to result in local minima than are multimodal restraints, such as torsion angles; multimodal functions are better used as validation criteria. It is fortunate that validation programs, such as PROCHECK (Laskowski et al., 1993), make good use of multimodal functions such as torsion angles and hydrogen-bonding patterns that are not employed as restraints in SHELXL refinements.