2.24 **********************
* LSQROTGEN WRITE-UP *
**********************
This program refines the orientation and position of a general
noncrystallographic symmetry operator by least squares. It is
applicable to any noncrystallographic symmetry transformation,
including arbitrary rotation angles and post rotation translations.
The operator being refined may relate regions of density within the
same crystal or in different crystals, so that cross-crystal averaging
is possible. The input map(s) (which MUST be orthogonal) are read in
along with control information specifying initial values for the
operator, the number of crystals, and what areas in the map(s) to
consider. The map areas considered can be all points within a given
distance from arbitrary input points, all points within input masks,
or all points simultaneously satisfying both conditions.
The input map(s) usually encompasses a dimer, trimer etc and are
extracted from FSFOUR maps via programs MAPVIEW or EXTRMAP. The input
masks, if used, are prepared by MAPVIEW or MDLMSK and must correspond
precisely to the input maps. In the single crystal case if masks are
used, different mask values corresponding to the different molecules
present allow a single mask file to be used (see MAPVIEW). One merely
specifies which molecules (mask numbers) are to be used in refinement.
If the input map (and its corresponding mask) do not correspond to an
orthogonal system, program MAPORTH should be used to convert them
before they can be used for operator refinement here. After each cycle
of refinement the correlation coefficient is printed along with the
new parameters. In the two crystal case a scale factor relating the
density within the appropriate envelopes in each crystal is
automatically refined. In most cases one should start with a low
resolution map (roughly 6 angstrom data, on a 2 angstrom grid) in case
the initial operator parameters are inaccurate. As the refinement
converges higher resolution data and finer map grids should be used.
It is often sufficient to refine considering only density within
spheres of radius 15-25 angstroms, centered on either the aggregate
centroid (if pure rotational symmetry is present), or centered on
each molecules center of gravity (the general case). This enables one
to refine the operator without the need for an input mask (although a
mask will be needed later for averaging). Usually correlation
coefficients of about 0.4 or higher (in a 4 angstrom map) indicate the
noncrystallographic symmetry axis is well positioned, and that
averaging will be useful.
INPUT DATA (UNIT 5)
CARD I PAMFIL (free format)
PAMFIL = Input file specifying cell and symmetry
parameters, used only to get "running log" file
CARD II NCYCLE, NCRYST (free format)
NCYCLE = No. of cycles of least squares refinement.
NCRYST = No. of crystals (1 or 2). Normally 1 but for
cross-crystal averaging it should be 2.
CARD III MAPFILE1 (free format)
MAPFILE1 = Name of file containing map for crystal 1.
***** Include card IIIA ONLY if NCRYST=2 *****
CARD IIIA MAPFILE2 (free format)
MAPFILE2 = Name of file containing map for crystal 2.
CARD IV PHI, PSI, CHI, OX, OY, OZ, T (free format)
Spherical polar angles defining direction and rotational
PHI = order of noncrystallographic symmetry axis, oriented with
respect to orthogonal frame with X along a, Y along c*
cross a, and Z along x cross y (i.e. c*). PSI = angle
PSI = between NC symmetry axis and +Y axis. PHI = angle between
projection of NC symmetry axis on XZ plane and
CHI = +X axis. +PHI = CCW rotation about +Y axis as measured
from +X axis. +CHI = CW rotation about the directed axis,
when viewed from the +axis toward the origin
OX =
Origin of NC symmetry rotation axis, in angstroms with
OY = respect to the orthogonal axes. The axis passes through
this point.
OZ =
T = Post rotation translational shift (in angstroms) parallel
to the rotation axis.
Note that the transformation operator refined is defined as that which
moves molecule "A" to molecule "B" via
Xb = (Rm) (Xa - Xo) + Xo + T*Rx
where Rm is a 3x3 rotation matrix expressed in terms of the spherical
polar angles, Xb, Xa are 3 element column vectors containing new and
old coordinates, respectively, Xo is a 3 element column vector
containing coordinates of the origin point for the rotation axis, T is
a post rotation translation shift scalar (screw like, in angstroms)
and Rx is a 3 element column vector containing direction cosines of
the rotation axis. All components of the operator are given in terms
of orthogonal coordinates in Angstroms, and the operator is applied
to (and yields) orthogonal coordinates.
For cross-crystal averaging applications molecule "A" is assumed to
reside in crystal 1 and molecule "B" in crystal 2. The ORTHOGONAL
coordinate systems in both crystals are then simply superimposed.
Note that the operator is defined relative to only the (common)
ORTHOGONAL axes, so one need not be concerned about its orientation
realatve to each set of CRYSTAL axes.
CARD V ISPHERE_A, MSK_A (free format)
ISPHERE_A = 0 consider all points in map (subject only to
MSK_A criteria below) for molecule A.
= 1 consider only grid points within specified
sphere for molecule A. (also subject to
MSK_A criteria below).
MSK_A = 0 no mask input for molecule A.
= 1 consider only grid points within envelopes
specified by input mask for molecule A. (also
subject to ISPHERE_A criteria)
Note that ISPHERE_A and MSK_A should not BOTH
be 0, although both can be 1 in which case
both criteria are applied for grid point
selection.
***** Include cards VA and VB ONLY if MSK_A = 1 *****
CARD VA MASK_A (free format)
MASK_A = Mask number (from 1-12) identifying points within
molecular envelope for "A". The value should
correspond to that used during mask creation.
CARD VB MASKFILE1 (free format)
MASKFILE1 = Name of file containing mask for crystal 1.
***** Include card VC ONLY if ISPHERE_A = 1 *****
CARD VC XCENA,YCENA,ZCENA,RADA (free format)
XCENA =
Sphere center, in Angstroms, with respect to orthogonal
YCENA =
coordinate system, situated in molecule "A".
ZCENA =
RADA = Sphere radius, in Angstroms, for molecule "A"
CARD VI ISPHERE_B, MSK_B (free format)
ISPHERE_B = 0 consider all points in map (subject only to
MSK_B criteria below) for molecule B.
= 1 consider only grid points within specified
sphere for molecule B. (also subject to MSK_B
criteria below).
MSK_B = 0 no mask input for molecule B.
= 1 consider only grid points within envelopes
specified by input mask for molecule B. (also
subject to ISPHERE_B criteria)
Note that ISPHERE_B and MSK_B should not BOTH
be 0, although both can be 1 in which case
both criteria are applied for grid point
selection.
***** Include cards VIA and VIB ONLY if MSK_B = 1 *****
CARD VIA MASK_B (free format)
MASK_B = Mask number (from 1-12) identifying points within
molecular envelope for "B". The value should
correspond to that used during mask creation.
CARD VIB MASKFILE2 (free format)
MASKFILE2 = Name of file containing mask for crystal 2.
(may be same as MASKFILE1, if NCRYST=1).
***** Include card VIC ONLY if ISPHERE_B = 1 *****
CARD VIC XCENB,YCENB,ZCENB,RADB (free format)
XCENB =
Sphere center, in Angstroms, with respect to orthogonal
YCENB =
coordinate system, situated in molecule "B".
ZCENB =
RADB = Sphere radius, in Angstroms, for molecule "B"
CARD VII ( IVAR(I), I=1,7 ) (free format)
Variable selection information
IVAR(1) = 1 to refine PHI, 0 to hold fixed
IVAR(2) = 1 to refine PSI, 0 to hold fixed
IVAR(3) = 1 to refine CHI, 0 to hold fixed
IVAR(4) = 1 to refine OX, 0 to hold fixed
IVAR(5) = 1 to refine OY, 0 to hold fixed
IVAR(6) = 1 to refine OZ, 0 to hold fixed
IVAR(7) = 1 to refine T, 0 to hold fixed
NOTES: Input maps (and masks, if used) must be orthogonal. If the
crystal systems do not have orthogonal axes, program MAPORTH must be
run to orthogonalize the maps (and masks, if used).
If masks are input, they must coincide exactly with their
corresponding maprs.
Normally all parameters are refined, but occasionally one must use
the IVAR selection flags to hold one or more parameters fixed.
An example would be the case where PSI is close to 0, in which case
PHI is then indeterminate (and irrelevant!). One could then hold PHI
fixed to avoid matrix singularities. Also, in cases where pure
translations are involved, one could hold all of the angles fixed.
The translation shift T is for a translation parallel to the
rotation axis (screw like) as translations in any other direction
can be achieved simply by changing the rotation axis origin. An
initial estimate of T can be obtained from two points P1, P2
related by the NC symmetry from
T = DX cos(PHI)sin(PSI) + DY cos(PSI) -DZ sin(PHI)sin(PSI)
where DX = P2x-P1x, DY = P2y-P1y, DZ = P2z-P1z
and the P's are expressed in the orthogonal axial system.
Note the directionality of the transformation (P1 going to P2 as
opposed to P2 going to P1) affects the sign of T (and CHI). If
the transformation operator is available as a simple 3x3 matrix
and 1x3 vector WHICH OPERATES ON ORTHOGONAL COORDINATES AS IN THE
PROTEIN DATA BANK FRAMEWORK, then the program O_to_sp can be used
to convert that representation to the spherical polar angles, axis
offset and post rotation translation needed here. An example of
this usage would be to get the matrix and vector from one of the
"lsq" options in the graphics program "O", and use o_to_sp to
convert the information to PHASES style.
******** FILES ********
INPUT MAP FILES (BINARY)
record 1) A,B,C,AL,BE,GA,NX,NY,NZ,IXMN,IYMN,IZMN,IXMX,IYMX,IZMX
with first 6 values REAL*4, next 9 INTEGER*4, lengths in Angstroms,
angles in degrees.
NX =
Number of grid points defining one "cell length" along
NY = respective axis. Implicitly defines grid spacing as
del x = A/NX, del y = B/NY and del z = C/NZ
NZ =
IXMN, IXMX =
Minimum, maximum grid index defining map region such
IYMN, IYMX = that x (fractional) = IX * (del x) / A etc.
There are no restrictions on magnitudes or signs.
IZMN, IZMX =
The map follows as (IYMX-IYMN+1)*(IZMX-IZMN+1) records, with
each containing one row (IXMX-IXMN+1 REAL*4 values) along X,
starting at IXMN. Y is slowest varying, i.e. the file could have
been created with the following FORTRAN code:
DO 30 IY=IYMN,IYMX
DO 20 IZ=IZMN,IZMX
20 WRITE(LU)(RHO(IX,IY,IZ),IX=IXMN,IXMX)
30 CONTINUE
INPUT MASK FILES (BINARY), if needed
Header record identical to map file.
Mask records similar to normal map records except that the mask values
are written as FORTRAN type "BYTE" (INTEGER*1). Only grid points with
mask values of 0, 10, 20, 30, 40 etc will be used (i.e. inside
envelope masks 1,2,3,4,5 etc, respectively).