POLARRFN (CCP4: Supported Program)


polarrfn - fast rotation function that works in polar angles, updated to work with up to the 100th order of the spherical harmonics


polarrfn hklin foo_in_Fobs.mtz hklin2 foo_in_Fmodel.mtz mapin bar.map mapout foo.map PLOT fred.plt
[Keyworded input]


This is a fast rotation function that works in polar angles, written by W. Kabsch. The program produces sections of constant rotation angle kappa for different axis directions defined by omega (angle from pole) and phi (angle around equator). The direction cosines of the axis is given by:
        (  sin omega  cos phi  )
        (  sin omega  sin phi  )
        (     cos omega        )

These sections may be plotted as stereograms (PLOT option). For example, the kappa = 180 section can be used to identify 2-fold axes in a self-rotation function.

The program searches for peaks in the rotation function (FIND option), and these are listed in the log output in both polar and Eulerian angles. For each peak, peaks related by crystallographic symmetry are also listed. If the peak, or one of its symmetry relations, corresponds to the identity operation (i.e. zero rotation), then this is indicated by "Origin peak" in the output. For a self-rotation function, peaks corresponding to genuine non-crystallographic symmetry will appear after the origin peaks.

Options are provided in the program for writing the calculated rotation function to a map file, that can be read by the program on another occasion for plotting, peak searching etc (options MAP, READ). The results of the first stage of the calculation (ALMN) may also be saved for subsequent rotation function summations on different grids (SAVE, SUM). A list of peaks may be read in, and all symmetry related peaks generated (option PEAK).

Compared to the Crowther program (ALMN), this program is a little slower but a self rotation function in polar angles is easier to understand. This program can also be used to calculate a small part of the rotation function on a fine grid (down to 1 degree steps). However, the crystal symmetry is expressed less fully in polar angle space than in Eulerian angle space, so for cross-rotation functions the Eulerian program is normally more appropriate.

This version of POLARRFN was modified for CCP4 5.0 to use up to the 100th order of spherical harmonics (previously 30). This allows a larger radius of integration and/or a higher resolution cutoff to be used.


Only the first 4 letters of each keyword are necessary. All input is free format. Most data cards are optional. The various data control lines are identified by keywords, those available being:
The following keywords are compulsory: SELF/CROSS, CRYSTAL, LABIN.

TITLE <title>

The rest of the line is taken as a title.

SELF <arad> <eps>

Set flag for self rotation (cf CROSS)
integration radius in Patterson space (default 20Å)
Radius (Angstrom) of sphere within which the Patterson is removed. Default Eps set = 0.7*<resmax>.

CROSS <arad> <eps>

Set flag for cross rotation (cf SELF)
integration radius in Patterson space (default 20Å)
Radius (Angstrom) of sphere within which the Patterson is removed. Default Eps set = 0.7*<resmax>.

RESOLUTION <resmin> <resmax>

Read resolution limits in Å, in either order. If only one treat as high resolution limit. Note that because of the limit of 100 in the order of the spherical harmonics in this version, the high resolution limit cannot be numerically less than the integration radius (arad) / 17.4: if it is, the program resets <resmax> to <arad> / 17.4 .

CRYSTAL FILE <nfile> ORTH <ncode> BFAC <badd> FLIMITS <fmin> <fmax> [CELL SYMMETRY]

This card is COMPULSORY for each hklin set (one for self-rotation, two for cross-rotation). The first subkeyword FILE (or NUMBER) indicates whether this CRYSTAL card refers to crystal 1 or to crystal 2.

followed by crystal number <nfile>. This sets the crystal number, 1 or 2 (default 1). The syntax is any one of FILE 1 or FILE 2.
followed by orthogonalisation code <ncode> (default 1). If the PEAK option is used ORTH should be given. Ncode defaults to 1 if unspecified.
            ncode      orthogonalisation code for this crystal
                  = 1, orthogonal x y z  along  a,c*xa,c*  (Brookhaven, default)
                  = 2                           b,a*xb,a*
                  = 3                           c,b*xc,b*
                  = 4                           a+b,c*x(a+b),c*
                  = 5                           a*,cxa*,c   (Rollett)
followed by temperature factor <badd> (default = 0). This can be used to sharpen the data (e.g. = -20). Fused = Finput * exp(-badd*ssq/lsq). A better alternative is to use normalised amplitudes (from the program ECALC), in which case the BFAC option must NOT be used.
followed by <fmin> and <fmax> (defaults to no cutoffs). These tests are done after application of the temperature factor, but before squaring to intensities.
followed by symmetry operation or crystal space group name or number. Default: take from MTZ file.
followed by cell dimensions. Default: take from MTZ file.

LABIN <data assignments>

LABIN FILE <nfile> F=<label> SIGF=<label2>
File column data assignments for F, SIGF . If you use the output from ECALC (recommended), set F=E and SIGF=SIGE .

FILE sets the crystal number, 1 or 2 (default 1). The syntax is any one of FILE 1 or FILE 2.

LIMITS <iphi1> <iphi2> <nphi> <iomega1> <iomega2> <nomega> <ikappa1> <ikappa2> <nkappa>

i.e. limits and steps on phi, omega, kappa
      iphi1 iomega1 ikappa1   start points in degrees
      iphi2 iomega2 ikappa2   stop  points in degrees
      nphi  nomega  nkappa    intervals in degrees
The default start & stop points are 0 & 180 which covers the basic asymmetric unit for the self-rotation function; the default step size is set numerically equal to the high resolution cutoff of the data, or that specified by the RESOLUTION keyword, whichever is numerically the greater. It is recommended to omit specification of the angle limits, at least for self-rotation functions.

[NO]PRINT <lprint>

Switch on or off printing of map as a figure-field in the logfile, and set print threshold LPRINT. The rotation function is scaled to a maximum of 100. Values less than 0 are printed as ' -' and values between 0 and LPRINT are printed as ' .'. NOPRINT (= PRINT 0) suppresses printing of the figure-field. The default is to print the map (PRINT +1), except in the case of the READ option, for which the default is NOPRINT.


Write rotation function to a map file. This may be read back later for plotting, peak searching, etc. using the READ option.

MAXORD <maxord>

This is primarily for test purposes, and should not normally be set. If you set MAXORD 30, you should get the same result as the original 30-order version, all else being equal. This will have the effect of placing a numerically greater limit on the high resolution cutoff. Note however that the old version does not set the default angle step equal to the high resolution cutoff of the data, and it does not take the latter into account when determining the maximum order to be used. What this means in practice is that to get the same result you should set the resolution cutoffs and angular steps explicitly, as the default values used by the two versions of the program are likely to be different.

PLOT <cstart> <cint> [ NOTITLE ]

Plot rotation function as a stereographic projection of each kappa section. These sections have omega = 0 or 180 at the centre, omega = 90 around the edge, and phi as marked around the periphery.
contour start level (default = 10)
contour interval (default = 10)
It is recommended to set the interval equal to the RMS deviation of the map, and the start level to twice this (this means you have to run the program twice, the first time to get the value of the RMS deviation printed near the end of the log file and then again to get the correctly contoured plot).

If the keyword NOTITLE is present, the maps are plotted with no labels or titles: this may be useful for publication purposes. The plots are not normally useful if a small part of the function is being calculated.

JOIN <kappa> <omega1> <phi1> <omega2> <phi2> [DASHED [<repeat> <dash>]]

Draw arc connecting two peaks at (omega1,phi1) and (omega2,phi2) on section kappa.

If the keyword DASHED is present, the arc will be dashed with repeat and dash length as specified (in multiples of the radius of the stereogram). Default values for <repeat> and <dash> are 0.05, 0.025. Up to 30 JOIN commands may be given.

Warning: this option is primarily intended for self-rotation functions, and probably will not work sensibly if the kappa sections are split into two sections or half sections.

FIND <peak> <maxpeak> [RMS] [OUTPUT <filename>]

Read peak threshold and maximum number of peaks. Peak searching may be suppressed by PEAK 0 0.

Note that although the program attempts to interpolate the positions of peaks, this may not be very accurate in regions where the peaks are very spread out (kappa near 0) or on the edges of the map.

threshold for peaks (default = 10).
maximum number of peaks to find (default = 10)
Up to <maxpeak> peaks above <peak> will be found, and all symmetry related peaks generated

If the keyword RMS is present, then the peak threshold is <peak> * RMS density, otherwise <peak> is the absolute threshold in the scaled map.

If the keyword OUTPUT is present, the unique peaks will be output to a file whose name may also be given (default filename or logical name PEAKS). The keyword RMS must precede OUTPUT if both are present.

READ <kappa1> <kappa2>

Instead of calculating a rotation function, read a previously calculated one (written out using the MAP flag). Sections <kappa1> to <kappa2> will be processed (if no values are given, the whole map will be used), printed (if PRINT is on, NOPRINT is default), plotted (if a PLOT card is present) and searched for peaks (see FIND).

SAVE [ <fname> ]

save output from first stage of rotation (calculation of Almn coefficients for each crystal) in files <fname>1.dat and <fname>2.dat. <fname> defaults to 'COEFFS'. These files can be used to save time on a later run of the program (see SUM card) to calculate a rotation function on a different grid (e.g. on a fine grid around a peak).

SUM [ <fname> ]

read Almn coefficients from files generated in a previous run of the program with the SAVE option. The files are <fname>1.dat and <fname>2.dat <fname> defaults to 'COEFFS'


read a list of peak positions (omega, phi, kappa in degrees) from subsequent lines, terminated by blank line, 'END' or end of file. The program will not calculate a rotation function, but will generate all symmetry related peaks from the peaks given, and print out the corresponding matrices. This option requires cards SELF or CROSS, CRYSTAL or CELL, and SYMMETRY (ie the number of crystals, their symmetry and orthogonalisation codes must be defined). If the keyword MATRICES is present, the rotation matrices corresponding to each peak will be printed.


Terminates input. If present must be the last card.


The polar angles omega, phi, kappa have an intrinsic symmetry operation 180-omega, 180+phi, -kappa by definition. This gives an asymmetric unit in the absence of other symmetry of omega 0 to 180, phi 0 to 360, kappa 0 to 180. Additional symmetry relations between omega and phi for a given kappa are created in self rotation functions, and by parallel dyad axes in the two crystals (ie any dyad axis in a self rotation function). These are as follows
  1. Self-rotation function, and also kappa=180 in all cases:
    kappa = -kappa, hence 180-omega, 180+phi
  2. dyad axes in both crystals parallel to orthogonal x axes (after any permutation produced by the orthogonalisation code NCODE):
    180 - omega, -phi
  3. dyad axes parallel to orthogonal y axes:
    180 - omega, 180 - phi
  4. dyad axes parallel to orthogonal z axes:
    omega, 180 + phi

If more than one of these symmetries is present, additional symmetry is generate by their combination. The asymmetric units required are given in the following table:

Symmetry                            Asymmetric unit
--------                           omega   phi  kappa

      no symmetry                   180    360   180

1.   180 - omega, 180 + phi          90    360   180
                                 or 180    180   180

2.   180 - omega, -phi               90    360   180
                                 or 180    180   180

3.   180 - omega, 180 - phi          90    360   180
                                 or 180    -90)  180
                                        to +90)

4.   omega, 180 + phi               180    180   180

In the program, the PLOT option will always try to plot a complete stereogram of all kappa sections. In the absence of symmetry, each kappa section is split into two hemispheres, omega 0 to 90 and omega 90 to 180. For symmetries 1, 2 and 3, only the hemisphere omega 0 to 90 is plotted. For symmetry 4, the two quarter spheres phi 0 to 180 for omega 0 to 90 and omega 90 to 180 are plotted together in the same circle. The PLOT option is not normally useful if only a small part of the asymmetric unit is calculated.


Input amplitudes for crystal 1, or for sole crystal in the case of self-rotation function.
Input amplitudes for crystal 2, for cross-rotation function.
Input map for READ option.
Output map containing rotation function, for MAP option.
Plotter output containing stereographic projection of each kappa section, for PLOT option.
Almn coefficients for crystal 1. Default created as scratch file.
Almn coefficients for crystal 2. Default created as scratch file.



polarrfn hklin tpfktrunc.mtz 
             mapout tself.map 
             plot self.plt 
             << eof
TITLE T-PFK self-rotation  R=29A, 5 - 7 A
PLOT 10 10
FIND 2 20 RMS OUTPUT peaks.dat

#!/bin/csh -f
#  Self-rotation function
set resolution = 15 3
set radius = 20
set u1af = {$scr0}/u1afobs_fc117

#  Calculate whole map, & search for peaks
polarrfn hklin {$u1af} mapout u1aself << eof-1
title U1A Self-Rotation function, resolution ${resolution} radius ${radius}
self  ${radius}
resolution  ${resolution}
crystal  file 1 bfac -20 
labin  file 1  F=FP SIGF=SIGFP
find  5 100

# Now plot just the kappa = 180 deg section
polarrfn hklin {$u1af} mapin u1aself plot u1aself << eof-2
title U1A Self-Rotation function, resolution ${resolution} radius ${radius}
self  ${radius}
resolution  ${resolution}
crystal  file 1 bfac -20 
labin  file 1  F=FP SIGF=SIGFP
read 180 180
plot  10 10

Cross rotation

polarrfn hklin  tpfktrunc 
             hklin2 uniqsfpfk.mtz
             mapout tcross.map 
             << eof
TITLE T-PFK cross-rotation  R=29A, 5 - 7 A
LIMITS 0 180 5 0 360 5 0 180 5
FIND 10 20


Author:    W. Kabsch

Contacts:  Phil Evans   (pre@mrc-lmb.cam.ac.uk)
           Ian Tickle   (i.tickle@astex-technology.com)