PROGRAM REFIX C ============= C C C REFIX C C C C This program is for academic distribution only C ============================================== C C C Determines unit cell orientation and refines diffraction C parameters using a list of observed reflection positions C C C C WOLFGANG KABSCH C MAX-PLANCK-INSTITUTE FOR MEDICAL RESEARCH C DEPARTMENT OF BIOPHYSICS C JAHNSTRASSE 29 C 6900 HEIDELBERG FRG C C To whom acknowledgement should be made C C C C REFERENCES C C "AUTOMATIC INDEXING OF ROTATION DIFFRACTION PATTERNS" C Kabsch,W. (1988),J.Appl.Cryst.21,67-71. C C C C MODIFIED BY HOWARD TERRY & ZBYSZEK DAUTER (EMBL HAMBURG) SEP.88 C TO READ OUTPUT FROM MOSCO STILLS PROGRAM (SPOTS.DAT) C C C C C 8-11-88 Modified for keyworded input. C Allow threshold intensity cutoff for spot selection. C Correct for film tilt on scanner. C A.G.W. Leslie, LMB, Cambridge C C 17-3-89 Bug in cassette tilt correction (notified by Duilio Cascio) C corrected. C C 8-1-90 Modified to read in a defining U matrix. C C In high symmetry space-groups there are a number C of possible auto-indexing solutions which are just C permutations of axes. C C This has the annoying consequence that misetting C angles close to zero may come out (for instance) C as close to 90 -180 -270 . We can obviate this by C permuting the REFIX solution under the Laue C Group operators and choose the one which gives C a rotation matrix closest to the given U matrix. C C This is useful in three respects-(1) a suitable C U Matrix can be chosen to make the missetting C angles small (necessary for POSTCHK at the moment) C (2) if we have two REFIX solutions from different C images we can impose the U Matrix of the first C on the second and so see how much the missetting C angles have changed (ie rather than C being c. 90 ,180, -270 degrees etc. C away from each other, the missetting angles will C be small)(3)If a crystal has slipped, when we C impose a defining U matrix we can keep the indexing C consistent throughout. C C The defining U matrix is read in with the C UMAT card (q.v.) or C it can be calculated from a given A matrix. C The A matrix is read C with the MATRIX card (q.v) C either as nine values or read from C a MATRIX.DAT file produced by a different run of REFIX C Paul McLaughlin, LMB, Cambridge C C 13-Mar-1990 Change default pixel size ("Q") back to 0.05 mm C A.G.W.L. C C 14-Mar-1990 Merge the film and image plate versions, C can now choose mode of operation by keyword FILM or IMAGE. C Add option of supplying camera constants. A.G.W. Leslie C C 21-Mar-1990 Very minor mods to write more output to SYS$OUTPUT, C since as it was originally the final refined cell (based C on using all reflections) was not given. A.G.W.L. C C 23-Jun-1990 Renormalise the ROT matrix before working out missetting C angles using UNORM routine.This stops INVALID ARGUMENT TO C MATH LIBRARY errors. PMCL C C 25-Jun-1990 Correct error in how work out the rotation axis and angle C of matrices in RTUMAT. Use DCOSFD. PMCL C C 26-Jul-1990 For image plates, the detector origin was not corrected C for ccomega (when given on CCOM card). This led to failure C when ccomega was given as -0.3. A.G.W.L. C C -Nov-1990 Toolpacked by K. Henrick C C 28-Nov-1990 Correct infinite format loop if detector origin not C specified for image plate. A.G.W.L. C Add keyword "LMB" (equivalent to SCRn keywords) to define scanner C origin for LMB scanner. A.G.W.L. C C 8-Jan-1991 Change the format of the spots file produced by IMSTILLS C so that it has "pseudo-fiducials" present. This eliminates the necessity C of differentiating between film and image plates. The keyword FILM C or IMAGE now only affects the default pixel size. Write the UMAT, C cell parameters and missetting angles to the output file MATRIX C in addition to AMAT and missets. Then the user has the option of C either using the AMAT or the UMAT when this file is read by IDXREF C or by OSCGEN. Reading UMAT has some advantages, as it is then easier C to update cell parameters explicitly. Also it allows the UMAT option C described above to be carried over into IDXREF and OSCGEN without C having to edit the matrices into the command files. Note that C this version of the program MUST be used in conjunction with the new C version of IMSTILLS, which adds the "pseudo fiducial" coordinates. C C It also means that REFIX no longer needs to know about the direct C beam position, as this is worked out from the "pseudo fiducial" C coordinates, so the keyword CENTRE is now redundant and has been C removed. Corrections to the image centre should be made using C the CAMCON keyword. A.G.W.L. C C 14-Nov-1991 Change so that output is written to unit 6 when running C in batch (previously only went to .LP file) A.G.W.L. C C FILES USED BY THIS ROUTINE C C FILE NAME DESCRIPTION STATUS C ------------------------------------------------------------------- C C REFIX.DAT - data cards for the program (given) C SPOT.DAT - a list of spots. on return indices (modified) C are attached to each spot. C REFIX.PARM - refined diffraction parameters (result) C REFIX.LP - printed messages and results (if runnign online) (result) C MATRIX.DAT - orientation matrix for oscgen C C Note SPOT and MATRIX are logical filenames and may be assigned to C other files. C C C DESCRIPTION OF INPUT PARAMETERS C C The input parameters necessary to run this program are assumed to C be present in the sequential formatted file "REFIX.DATA" in the C current directory. The laboratory frame can be any (!) orthonormal C right-handed coordinate system chosen by the user. C C **** INPUT IS NOW KEYWORDED, see below *** C C line # contents C 1 FLAMDA- x-ray wavelength (Angstroem) C S0 - x,y,z components of the direct beam wavevector with C respect to the laboratory coordinate system. C Direction of "S0" is from X-ray source towards the C crystal. C ACHSE - x,y,z components of rotation axis with respect to C laboratory system. C 2 IGROUP- Space-group number of crystal. The numbers C corresponding to each possible space-group C are defined in the "INTERNATIONAL TABLES I". C CELL - Unit cell parameters a,b,c (Angstroem) and C alpha,beta,gamma (degrees). C 3 Q - Length of a detector pixel (mm). C F - Signed distance between crystal and detector (mm). C Note that "F" may be negative depending upon "ED". C ORGX - Detector X-coordinate (pixels) of origin, C ORGY - Detector Y-coordinate (pixels) of origin. C ED - Orientation matrix of the detector with respect to C the laboratory system. The lab coordinates of C a point IX,IY (pixels) on the detector are : C x=Q*(IX-ORGX)*ED(1,1)+Q*(IY-ORGY)*ED(1,2)+F*ED(1,3) C y=Q*(IX-ORGX)*ED(2,1)+Q*(IY-ORGY)*ED(2,2)+F*ED(2,3) C z=Q*(IX-ORGX)*ED(3,1)+Q*(IY-ORGY)*ED(3,2)+F*ED(3,3) C Only the following values must be given: C ED(1,1),ED(2,1),ED(3,1),ED(1,2),ED(2,2),ED(3,2) C The third unit vector is then calculated by the C program as ED(.,3)=ED(.,1) X ED(.,2) C 4 FIXF - Refinement control variable C 0: refine detector to crystal distance, direct C beam direction and unit cell orientation. C The cell parameters are fixed. C 1: fix detector to crystal distance and refine C direct beam direction, unit cell orientation C and independent cell parameters. C DPHI - Expected error (degrees) of spindle position of C the spots given on file "SPOT.LIST" C C KEYWORDED INPUT C================ C C Keyworded input has been added using the PARSER from the MOSFLM C package. Refer to above for meaning of variables. Defaults for C the direct beam vector, rotation axis orientation and detector C orientation matrix are those appropriate to the MOSFLM package, C namely X along X-ray beam, Z along rotation axis. C C obligatory keywords C =================== C C SPACEGROUP Spacegroup number C CELL Unit cell parameters (real) C DISTANCE Crystal to detector distance (mm) C RUN autoindex. C C Optional keywords C ================= C C WAVE Wavelength (Default 1.5418) C PIXEL Pixel size (mm) (Default 0.05 for film) C Pixel size (mm) (Default 0.150mm if C processing image plates, C THRESHOLD Intensity threshold for spot selection C (Default 0) C DELPHI Expected error in PHI (Default 0.5 degrees) C XRAY Direct beam vector S0 C ROTAXIS Rotaion axis vector C DETMAT Detector orientation matrix ED, 6 components C as specified above C FIXC If present, cell parameters are fixed and C detector distance is refined C UMAT Defining U matrix (9 values).The misetting C angles are referred to this matrix. C MATRIX (9 values or file-name).This allows input C of an A matrix from which a defining U C Matrix can be calculated.It can be typed C explicitly or read from a file (such as C MATRIX.DAT which was produced by a different C run of REFIX.) C C The defining U matrix is read in with the UMAT card (q.v.) or C it can be calculated from a given A matrix. The A matrix is read C with the MATRIX card (q.v) either as nine values or read from C a MATRIX.DAT file produced by a different run of REFIX C C FILM Sets default pixel size to 0.05mm C C IMAGE Sets default pixel size to 0.15mm C C CAMCON Allows camera constants to be supplied. C CCX x These values will override C CCY x those present in the SPOTS file. Any combination C CCOM x of CCX,CCY and CCOMEGA may be given. C C C C DESCRIPTION OF RESULT FILES C C FILE NAME DESCRIPTION C C SPOT.LIST - Sequential formatted file. Each record consists of C 3 real and 3 integer numbers as follows: C X - X-coordinate (pixels) of observed spot C Y - Y-coordinate (pixels) of observed spot C PHI - Spindle position (degrees) of observed spot C On return indices are attached to each spot C IH - h-index of spot C IK - k-index of spot C IL - l-index of spot C A spot which could not be indexed has IH=IK=IL=0. C On input only the first 3 numbers need to be present. C A minimum of 25 neighbouring spots is required but C more spots will lead to better parameter estimates. C The program accepts up to 3000 spots and ignores the C rest. The spots should all be selected from a narrow C oscillation range of about 5 degrees. In case the C spots are from a single oscillation film set all C values "PHI" equal to the spindle angle at the center C of the oscillation range. C REFIX.PARM - Refined diffraction parameters written as a single C unformatted sequential record. C ACHSE - real(3) lab coordinates of rotation axis C ED - real(3,3) lab coordinates of detector C orientation matrix. C ORGX - real X-origin on film C ORGY - real Y-origin on film C F - real crystal to film distance C S0 - real(3) direct beam wavevector C RCELL - real(6) reciprocal cell parameters C U - real(3,3) orientation matrix of crystal C SDU - real standard deviation of U C SDXY - real standard deviation of spot C locations (pixels) C SDPHI - real standard deviation of spindle C positions (degrees) C SDCELL - real(6) standard deviations of reciprocal C unit cell parameters C REFIX.LP - Printed messages and results (if running online) C C MATRIX.DAT - Contains orientation matrix AMAT and missetting C angles, followed by UMAT, cell parameters and C missetting angles in format suitable for reading C into IDXREF or OSCGEN (keywords MATRIX or UMAT) C C C C SPECIAL CASE OF A SUPPER-ROTATION CAMERA C C The lab coordinate system used is described in the reference. C C The rotation axis has coordinates: 0.0 , 1.0 , 0.0 C C Using data frames as produced by the HARVARD software C we have for a detector swing angle CHI C ED(1,1)=-cos(CHI) ED(2,1)= 0 ED(3,1)= sin(CHI) C ED(1,2)= 0 ED(2,2)=-1 ED(3,2)= 0 C C Since the vector ED(.,3) is pointing from the crystal towards the C detector we usually get a positive sign for "F". C C C C SAMPLE INPUT DATA (old version) C C 1.5418 0.0 0.0 1.0 0.0 1.0 0.0 C 19 133.0 56.4 110.0 90.0 90.0 90.0 C 0.2 120.0 256.5 237.3 -1.0 0.0 0.0 0.0 -1.0 0.0 C 1 0.45 C C KEYWORDED VERSION C CELL 107.6 107.6 123.4 90. 90. 120. C SPACEGROUP 155 C FILM C DISTANCE 80.0 C DELPHI 0.4 C THRESH 500 C RUN C C C Subroutines required: ABSHKL,AUTOIX,CRYSYS,INDIX,INVCEL, C INVERS,MATMUL,REFINE,SETMAT,UNORM C C C C .. Parameters .. INTEGER MXRFL,MXY,MXSPOT,MXBEST,NCYCLE,SPOT PARAMETER (MXRFL=5000,MXY=101,MXSPOT=200,MXBEST=4,NCYCLE=6, + SPOT=13) C .. C .. Scalars in Common .. REAL CCOM,CCX,CCY,DPHI,F,FLAMDA,ORGX,ORGY,Q,THRESH INTEGER FIXF,IGROUP LOGICAL DCOMP,FILM,LCAMC,TARGET INTEGER LP,LINOUT LOGICAL ONLINE C .. C .. Arrays in Common .. REAL ACHSE,CELL,ED,S0,TARMAT C .. C .. Local Scalars .. REAL ALPHA,BETA,CBETA,COSOM,D,DE,DTOR,DTR,FINT,GAMMA,OMEGA,ORGXN, + ORGYN,PHI,PHIXYZ,PI,SDPHI,SDU,SDXY,SINOM,SXOFF,SYOFF,VOLUME, + X,XF,XFN,XN,XOFF,XOFFN,Y,YF,YFN,YN,YOFF,YOFFN,Z INTEGER ABSTYP,I,IBEST,ICS,ICS0,IER,IFAIL,ISPOT,J,J1,JKIM, + MAT,NACC,NBEST,NFAR,NOVLAP,NRFL,NSPOT INTEGER NRFL0,NVEC C .. C .. Local Arrays .. REAL A(3,3),ASTV(3),B(3,3),BSTV(3),C(6),CD(3,3),CR(3,3),CSTV(3), + DUM(3,3),EDB(3,3,MXBEST),ESD(MXBEST),ESDCEL(6,MXBEST), + ESDPHI(MXBEST),ESDU(MXBEST),ESDXY(MXBEST),FB(MXBEST), + FHKL(4,MXSPOT,MXBEST),ORGXB(MXBEST),ORGYB(MXBEST),RCELL(6), + RCELLB(6,MXBEST),RMS(MXBEST),ROT(3,3),S0B(3,MXBEST), + SDCELL(6),TARINV(3,3),U(3,3),UB(3,3,MXBEST),UBM(3,3),UU(3,3), + V(3,MXRFL),VC(3,MXSPOT),XFID(3),YFID(3),UDEF(3,3) INTEGER NACCB(MXBEST) INTEGER*2 IH(MXRFL),IK(MXRFL),IL(MXRFL),INTY(MXSPOT),IPHI(MXRFL), + IX(MXRFL),IY(MXRFL),TABLE(MXY,MXY) CHARACTER MSG(11)*62 C .. C .. External Functions .. LOGICAL CCPONL EXTERNAL CCPONL C .. C .. External Subroutines .. EXTERNAL ABSHKL,AUTOIX,CCPOPN,CROSS,CRYSYS,INDIX,INVCEL,INVERS, + MATCOP,MATMUL,READDAT,REFINE,RTUMAT,SETMAT,UNORM C .. C .. External Block Data .. EXTERNAL REFBLK C .. C .. Intrinsic Functions .. INTRINSIC ASIN,ATAN,ATAN2,COS,NINT,SIN,SQRT C .. C .. Common blocks .. COMMON /REFCOM/IGROUP,FLAMDA,CELL(6),F,THRESH,Q,DPHI,ED(3,3), + ACHSE(3),S0(3),ORGX,ORGY,FIXF,TARMAT(3,3),TARGET,DCOMP, + LCAMC,CCX,CCY,CCOM,FILM COMMON /RINOUT/ ONLINE,LP,LINOUT C .. SAVE C .. Data statements .. DATA UDEF/1.0,0.0,0.0, 0.0,1.0,0.0, 0.0,0.0,1.0/ DATA (MSG(JKIM),JKIM=1,6)/ + ' No error. solution is unique. ', + ' !!! WARNING !!! Solution may not be unique.', + ' !!! ERROR !!! Best solution cannot explain data.', + ' !!! ERROR !!! Value of NBEST is ZERO.', + ' !!! ERROR !!! Insufficient number of accepted spots.', + ' !!! ERROR !!! Dimension of difference vector set less than 2.'/ DATA (MSG(JKIM),JKIM=7,11)/ + ' !!! ERROR !!! Illegal space group number.', + ' !!! ERROR !!! Illegal cell parameters.', + ' !!! ERROR !!! Illegal incident beam vector or rotation axis.', + ' !!! ERROR !!! Crystal to detector distance is zero.', + ' !!! ERROR !!! Undefined detector orientation or pixel size.'/ C .. C C C---- values for COMMON /RINOUT/ C LINOUT = 6 LP = 6 ONLINE = (CCPONL()) C C---- Open .lp file if running online C IF (ONLINE) THEN C IFAIL = 1 LP = 4 C C ********************************** CALL CCPOPN(LP,'refix.lp',4,1,80,IFAIL) C ********************************** C IF (IFAIL.EQ.-1) THEN WRITE (LINOUT,FMT=*) + ' !!! ERROR !!! Cannot open PRINTER' STOP END IF END IF C REWIND LP C PI = ATAN(1.0)*4.0 DTOR = PI/180.0 C C---- Open matrix output C IFAIL = 1 C C ********************************** CALL CCPOPN(MAT,'MATRIX',4,1,80,IFAIL) C ********************************** C IF (IFAIL.NE.-1) THEN C C---- read data cards C C ******* CALL READDAT C ******* C C---- Get lattice type of crystal C C ************************** CALL ABSHKL(IGROUP,CELL,ABSTYP) C ************************** C C---- get crystal system C C *********************** CALL CRYSYS(IGROUP,CELL,ICS) C *********************** C C---- Get reciprocal cell parameters C C ************************* CALL INVCEL(CELL,RCELL,VOLUME) C ************************* C WRITE (LP,FMT=6030) IGROUP,CELL,RCELL IF (ONLINE) WRITE (LINOUT,FMT=6030) IGROUP,CELL,RCELL C C---- Calculate detector orientation matrix C ED(1,3) = ED(2,1)*ED(3,2) - ED(2,2)*ED(3,1) ED(2,3) = ED(3,1)*ED(1,2) - ED(3,2)*ED(1,1) ED(3,3) = ED(1,1)*ED(2,2) - ED(1,2)*ED(2,1) C C *********** CALL UNORM(ED,J) C *********** C C---- Check parameters C IF ((Q.LE.0.0) .OR. (J.NE.0)) THEN WRITE (LP,FMT=6032) MSG(11) IF (ONLINE) WRITE (LINOUT,FMT=6032) MSG(11) STOP C C ELSE IF (F.NE.0.0) THEN C C X = S0(1)**2 + S0(2)**2 + S0(3)**2 Y = ACHSE(1)**2 + ACHSE(2)**2 + ACHSE(3)**2 IF ((FLAMDA.GT.0.0) .AND. (X.GT.0.0) .AND. (Y.GT.0.0)) THEN C C X = SQRT(X)*FLAMDA Y = SQRT(Y) C C DO 10 J = 1,3 ACHSE(J) = ACHSE(J)/Y S0(J) = S0(J)/X 10 CONTINUE C C IF (VOLUME.LE.0.0) THEN WRITE (LP,FMT=6032) MSG(8) IF (ONLINE) WRITE (LINOUT,FMT=6032) MSG(8) STOP C C ELSE IF (ICS.GT.0) THEN C C---- Read spot list (rather spots.dat) C IFAIL = 1 C C ********************************** CALL CCPOPN(SPOT,'SPOTS',3,1,80,IFAIL) C ********************************** C IF (IFAIL.NE.-1) THEN C C NRFL = 0 REWIND SPOT C C---- Read in camera constants ccx,ccy C and the fiducial coordinates C READ (SPOT,FMT=*,END=60,ERR=150) XOFF,YOFF XF = 0.0 YF = 0.0 C C DO 20 I = 1,3 READ (SPOT,FMT=*,END=60,ERR=150) XFID(I), + YFID(I) IF (I.NE.2) THEN XF = XFID(I) + XF YF = YFID(I) + YF END IF 20 CONTINUE C C---- Calculate tilt of cassette C OMEGA = ATAN2((YFID(3)-YFID(2)), (XFID(3)-XFID(2))) WRITE (LP,FMT=6000) OMEGA/DTOR IF (ONLINE) WRITE (LINOUT,FMT=6000) OMEGA/DTOR C C---- Allow for camera constants supplied by user C IF (LCAMC) THEN OMEGA = CCOM*DTOR + OMEGA SXOFF = XOFF SYOFF = YOFF C C---- Reset CCX, CCY if they have been supplied C IF (CCX.GT.-998.0) XOFF = CCX IF (CCY.GT.-998.0) YOFF = CCY IF (ONLINE) WRITE (LINOUT,FMT=6002) WRITE (LP,FMT=6002) WRITE (LP,FMT=6004) XOFF,YOFF,OMEGA/DTOR WRITE (LP,FMT=6006) SXOFF,SYOFF IF (ONLINE) WRITE (LINOUT,FMT=6006) SXOFF,SYOFF END IF C C COSOM = COS(OMEGA) SINOM = SIN(OMEGA) C C---- Correct film centre and camera constants for cassette tilt. C XOFFN = XOFF*COSOM + YOFF*SINOM YOFFN = -XOFF*SINOM + YOFF*COSOM XFN = XF*COSOM + YF*SINOM YFN = -XF*SINOM + YF*COSOM C C ORGXN = ORGX*COSOM + ORGY*SINOM ORGYN = -ORGX*SINOM + ORGY*COSOM C C---- Define film centre C ORGX = 0.5*XFN/Q ORGY = 0.5*YFN/Q 30 CONTINUE C C READ (SPOT,FMT=*,END=60,ERR=150) X,Y,Z,PHI,FINT IF (X.GE.-100.0) THEN IF (X.LT.-98.0) THEN GO TO 30 ELSE IF (FINT.LT.THRESH) THEN GO TO 30 ELSE C C---- Correct for CCOMEGA and tilt C XN = X*COSOM + Y*SINOM YN = -X*SINOM + Y*COSOM C C---- Correct for camera constants CCX, CCY C X = (XN-XOFFN)/Q Y = (YN-YOFFN)/Q IF (NRFL.LT.MXRFL) THEN NRFL = NRFL + 1 IX(NRFL) = NINT(X*10.0) IY(NRFL) = NINT(Y*10.0) 40 CONTINUE IF (PHI.GT.180.0) THEN PHI = PHI - 360.0 GO TO 40 END IF 50 CONTINUE IF (PHI.LT.-180.0) THEN PHI = PHI + 360.0 GO TO 50 END IF IPHI(NRFL) = NINT(PHI*100.0) IH(NRFL) = 0 IK(NRFL) = 0 IL(NRFL) = 0 GO TO 30 END IF END IF END IF 60 WRITE (LP,FMT=6042) NRFL,THRESH IF (ONLINE) WRITE (LINOUT,FMT=6042) NRFL,THRESH C C---- Get orientation of crystal and refined parameters C SDXY = 2.0 SDPHI = DPHI C C *********************************************** CALL AUTOIX(IGROUP,CELL,FIXF,Q,ORGX,ORGY,F,ED,ACHSE, + S0,RCELL,U,SDU,SDXY,SDPHI,SDCELL,MXBEST, + NBEST,ORGXB,ORGYB,FB,EDB,S0B,RCELLB,UB, + ESDU,ESDXY,ESDPHI,ESDCEL,ESD,RMS,NACCB, + MXSPOT,INTY,VC,FHKL,NSPOT,MXRFL,NRFL,IX, + IY,IPHI,IH,IK,IL,V,MXY,TABLE,NVEC,IER) C *********************************************** C IF ((IER.LT.1) .OR. (IER.GT.11)) THEN WRITE (LP,FMT=*) ' !!! FATAL PROGRAMMING ERROR !!!' IF (ONLINE) WRITE (LINOUT, + FMT=*) ' !!! FATAL PROGRAMMING ERROR !!!' STOP ELSE C C WRITE (LP,FMT=6032) MSG(IER) IF (ONLINE) WRITE (LINOUT,FMT=6032) MSG(IER) C C ********************** CALL INVCEL(RCELL,C,VOLUME) CALL SETMAT(RCELL,A) CALL MATMUL(U,A,B) CALL INVERS(B,A,VOLUME) C ********************** C WRITE (LP,FMT=6034) NACCB(1),SDXY,Q*SDXY,SDPHI, + ORGX,ORGY,F WRITE (LP,FMT=6036) (ED(J1,1),J1=1,3), + (ED(J1,2),J1=1,3),S0, ((A(J1,J),J=1,3),J1=1,3),C, + RCELL,SDCELL IF (ONLINE) THEN WRITE (LINOUT,FMT=6034) NACCB(1),SDXY,Q*SDXY, + SDPHI,ORGX,ORGY,F WRITE (LINOUT,FMT=6036) (ED(J1,1),J1=1,3), + (ED(J1,2),J1=1,3),S0, ((A(J1,J),J=1,3),J1=1,3), + C,RCELL,SDCELL END IF C C IF ((NSPOT.LT.1) .OR. (NBEST.LT.1)) THEN STOP ELSE C C WRITE (LP,FMT=6008) IF (ONLINE) WRITE (LINOUT,FMT=6008) C C DO 110 IBEST = 1,NBEST WRITE (LP,FMT=6038) IBEST,ESD(IBEST) IF (ONLINE) WRITE (LINOUT,FMT=6038) IBEST, + ESD(IBEST) C C DO 70 ISPOT = 1,NSPOT IF (ONLINE) WRITE (LINOUT,FMT=6040) ISPOT, + (VC(J,ISPOT),J=1,3),INTY(ISPOT), + (FHKL(J,ISPOT,IBEST),J=1,4) WRITE (LP,FMT=6040) ISPOT, + (VC(J,ISPOT),J=1,3),INTY(ISPOT), + (FHKL(J,ISPOT,IBEST),J=1,4) 70 CONTINUE C C DO 80 J = 1,6 C(J) = RCELLB(J,IBEST) 80 CONTINUE C C ********************* CALL INVCEL(C,CELL,VOLUME) C ********************* C DO 100 J = 1,3 DO 90 J1 = 1,3 UU(J,J1) = UB(J,J1,IBEST) 90 CONTINUE 100 CONTINUE C C ****************** CALL SETMAT(C,A) CALL MATMUL(UU,A,B) CALL INVERS(B,A,VOLUME) C ****************** C WRITE (LP,FMT=6034) NACCB(IBEST),ESDXY(IBEST), + Q*ESDXY(IBEST),ESDPHI(IBEST),ORGXB(IBEST), + ORGYB(IBEST),FB(IBEST) WRITE (LP,FMT=6036) + (EDB(J1,1,IBEST),J1=1,3), + (EDB(J1,2,IBEST),J1=1,3), + (S0B(J,IBEST),J=1,3), ((A(J1,J),J=1,3),J1=1, + 3),CELL, (RCELLB(J,IBEST),J=1,6), + (ESDCEL(J,IBEST),J=1,6) C IF (ONLINE) THEN WRITE (LINOUT,FMT=6034) NACCB(IBEST), + ESDXY(IBEST),Q*ESDXY(IBEST),ESDPHI(IBEST), + ORGXB(IBEST),ORGYB(IBEST),FB(IBEST) WRITE (LINOUT,FMT=6036) + (EDB(J1,1,IBEST),J1=1,3), + (EDB(J1,2,IBEST),J1=1,3), + (S0B(J,IBEST),J=1,3), ((A(J1,J),J=1,3),J1=1, + 3),CELL, (RCELLB(J,IBEST),J=1,6), + (ESDCEL(J,IBEST),J=1,6) END IF 110 CONTINUE C C---- Index all spots C IF (NRFL.GT.0) THEN C C IF (FIXF.GT.0) THEN ICS0 = ICS ELSE ICS0 = FIXF END IF C C NRFL0 = 0 C C ******************************************* CALL INDIX(ABSTYP,Q,ORGX,ORGY,F,ED,ACHSE,S0, + RCELL,U,MXRFL,NRFL0,NRFL,IX,IY,IPHI, + IH,IK,IL,NACC,NOVLAP,NFAR,SDPHI,SDXY) CALL REFINE(NCYCLE,NRFL,IH,IK,IL,IX,IY,IPHI, + ACHSE,ICS0,Q,ORGX,ORGY,F,S0,ED,U, + RCELL,SDU,SDCELL,SDPHI,SDXY,IER) C ******************************************** C IF (IER.LE.0) THEN STOP ' Error in REFINE' ELSE C C ************************************* CALL INVCEL(RCELL,CELL,VOLUME) CALL SETMAT(RCELL,A) CALL MATMUL(U,A,B) CALL INVERS(B,A,VOLUME) CALL INVCEL(RCELL,CELL,VOLUME) CALL CROSS(A,ASTV,BSTV,CSTV,VOLUME,FLAMDA) C ************************************* C C---- If a target matrix was given, find permutation of the UMAT C which is closest to it. C C *************************** IF (TARGET) CALL RTUMAT(ASTV,BSTV,CSTV) C *************************** C C WRITE (LP,FMT=6010) IF (ONLINE) WRITE (LINOUT,FMT=6010) WRITE (LP,FMT=6034) NACC,SDXY,Q*SDXY,SDPHI, + ORGX,ORGY,F WRITE (LP,FMT=6036) (ED(J1,1),J1=1,3), + (ED(J1,2),J1=1,3),S0, + ((A(J1,J),J=1,3),J1=1,3),CELL,RCELL,SDCELL IF (ONLINE) THEN WRITE (LINOUT,FMT=6034) NACC, + SDXY,Q*SDXY,SDPHI,ORGX,ORGY,F WRITE (LINOUT,FMT=6036) + (ED(J1,1),J1=1,3), (ED(J1,2),J1=1,3),S0, + ((A(J1,J),J=1,3),J1=1,3),CELL,RCELL,SDCELL END IF C WRITE (LP,FMT=*) + '"A" Matrix (equivalent to UB) for OSCGEN' IF (ONLINE) WRITE (LINOUT,FMT=*) + '"A" Matrix (equivalent to UB) for OSCGEN' WRITE (LP,FMT=*) + ' (This matrix is also written to logical file MATRIX)' IF (ONLINE) WRITE (LINOUT,FMT=*) + ' (This matrix is also written to logical file MATRIX)' WRITE (LP,FMT=6014) (ASTV(I),BSTV(I),CSTV(I), + I=1,3) IF (ONLINE) WRITE (LINOUT,FMT=6014) (ASTV(I), + BSTV(I),CSTV(I),I=1,3) C C PHIXYZ = 0.0 WRITE (MAT,FMT=6012) (ASTV(I),BSTV(I), + CSTV(I),I=1,3),PHIXYZ,PHIXYZ,PHIXYZ C DO 120 I = 1,3 UBM(I,1) = ASTV(I) UBM(I,2) = BSTV(I) UBM(I,3) = CSTV(I) 120 CONTINUE C C DTR = 0.0174533 CD(1,1) = RCELL(1) CD(1,2) = COS(RCELL(6)*DTR)*RCELL(2) CD(1,3) = COS(RCELL(5)*DTR)*RCELL(3) CD(2,1) = 0.0 CD(2,2) = SIN(RCELL(6)*DTR)*RCELL(2) CD(2,3) = -SIN(RCELL(5)*DTR)*RCELL(3)* + COS(CELL(4)*DTR) CD(3,1) = 0.0 CD(3,2) = 0.0 CD(3,3) = 1.0/CELL(3) C C **************** CALL INVERS(CD,CR,DE) C **************** C C---- If target U MAT was given, take this out of ROT C IF (TARGET) THEN C C *********************** CALL INVERS(TARMAT,TARINV,D) CALL MATMUL(CR,TARINV,DUM) CALL MATCOP(DUM,CR) C ************************ C END IF C C ****************** CALL MATMUL(UBM,CR,ROT) C ****************** C DO 140 I = 1,3 DO 130 J = 1,3 ROT(I,J) = ROT(I,J)/FLAMDA 130 CONTINUE 140 CONTINUE IF (ONLINE) WRITE (LINOUT,FMT=6018) WRITE (LP,FMT=6018) IF (ONLINE) WRITE (LINOUT, + FMT=6016) ((ROT(I,J),J=1,3),I=1,3) WRITE (LP,FMT=6016) ((ROT(I,J),J=1,3),I=1,3) C C---- Renormalise ROT C C ************ CALL UNORM(ROT,J) C ************ C BETA = ASIN(-ROT(3,1)) CBETA = COS(BETA) ALPHA = ASIN(ROT(3,2)/COS(BETA)) IF (ROT(3,3).LT.0.0) ALPHA = 3.14159 - ALPHA GAMMA = ASIN(ROT(2,1)/COS(BETA)) IF (ROT(1,1).LT.0.0) GAMMA = 3.14159 - GAMMA ALPHA = 57.29578*ALPHA BETA = 57.29578*BETA GAMMA = 57.29578*GAMMA IF (ONLINE) WRITE (LINOUT,FMT=6020) ALPHA, + BETA,GAMMA IF (.NOT.TARGET) WRITE (MAT,FMT=6023) + ((UDEF(I,J),J=1,3),I=1,3) IF (TARGET) WRITE (MAT,FMT=6023) + ((TARMAT(I,J),J=1,3),I=1,3) WRITE (MAT,FMT=6024) CELL WRITE (MAT,FMT=6024) ALPHA,BETA,GAMMA ccc endfile mat CLOSE (MAT) C WRITE (LP,FMT=6034) NACC,SDXY,Q*SDXY,SDPHI, + ORGX,ORGY,F WRITE (LP,FMT=6036) (ED(J1,1),J1=1,3), + (ED(J1,2),J1=1,3),S0, + ((A(J1,J),J=1,3),J1=1,3),CELL,RCELL,SDCELL IF (ONLINE) THEN WRITE (LINOUT,FMT=6034) NACC, + SDXY,Q*SDXY,SDPHI,ORGX,ORGY,F WRITE (LINOUT,FMT=6036) + (ED(J1,1),J1=1,3), (ED(J1,2),J1=1,3),S0, + ((A(J1,J),J=1,3),J1=1,3),CELL,RCELL,SDCELL END IF END IF END IF C C C C---- Save refined parameters C C---- Check if solution is acceptable C IF ((SDXY.GT.1.5) .OR. (SDPHI.GT.DPHI) .OR. + (NACC.LT. (3*NRFL/4))) THEN WRITE (LP,FMT=6028) SDXY,SDPHI,DPHI,NACC, + NRFL IF (ONLINE) WRITE (LINOUT,FMT=6028) + SDXY,SDPHI,DPHI,NACC,NRFL END IF C C CLOSE (LP) STOP END IF END IF GO TO 160 END IF C C 150 WRITE (LP,FMT=*) ' !!! ERROR !!! Cannot read SPOT.LIST' IF (ONLINE) WRITE (LINOUT,FMT=* + ) ' !!! ERROR !!! Cannot read SPOT.LIST' STOP ELSE WRITE (LP,FMT=6032) MSG(7) IF (ONLINE) WRITE (LINOUT,FMT=6032) MSG(7) STOP END IF ELSE WRITE (LP,FMT=6032) MSG(9) IF (ONLINE) WRITE (LINOUT,FMT=6032) MSG(9) STOP END IF ELSE WRITE (LP,FMT=6032) MSG(10) IF (ONLINE) WRITE (LINOUT,FMT=6032) MSG(10) STOP END IF END IF C C C C---- Format statements C 6000 FORMAT (/1X,'FILM TILT',F6.1,' DEGREES') 6002 FORMAT (/1X,'The camera constants supplied will override those p', + 'resent in the SPOTS file') 6004 FORMAT (/1X,'Camera constants used:',/1X,'CCX ',F6.3,'mm CCY', + F6.3,'mm combined tilt and ccomega',F7.3,' degrees') 6006 FORMAT (1X,'(Camera constants in SPOTS file are: CCX ',F6.3,'mm ', + ' CCY',F6.3,'mm)') 6008 FORMAT (//1X,'*** List of possible solutions ***',//) 6010 FORMAT (//1X,'** Final refinement using ALL reflections **',/) 6012 FORMAT (3F12.7) 6014 FORMAT (3 (3 (5X,F9.7),/)) 6016 FORMAT (5X,3F10.5) 6018 FORMAT (' U- Matrix') 6020 FORMAT (' Misset angles:',/' Phi(X) = ',F10.5,' Phi(Y) = ',F10.5, + ' Phi(Z) = ',F10.5) 6022 FORMAT ('UMAT 1 0 0 -',/' 0 1 0 -',/' 0 0 1') 6023 FORMAT (3F12.7) 6024 FORMAT (6F12.4) 6028 FORMAT (//,1X,'******** WARNING ********',/,1X, + ' One or more of the following tests failed:',/,1X, + ' Positional residual in pixels SDXY .gt. 1.5 ', + ' ; with SDXY = ',F10.5,/' or ',/,1X,'Standard deviation', + ' in PHI, SDPHI .gt. DPHI ; with SDPHI = ',F10.5,/,1X, + ' and DPHI = ', + F10.5,/' or ',/,1X,'Number of accepted reflections NACC ', + '.lt. (0.75*NRFL); with NACC = ', + I8,/' and NRFL = ',I8, + /,1X,'The determined orientation MAY be incorrect') 6030 FORMAT (' ***** REFIX ***** ',//' Space-group number ',I5, + /' Initial unit cell parameters ',3F8.2,3F8.3,/' Reciproc', + 'al cell parameters ',3F8.6,3F8.3,/) 6032 FORMAT (/1X,A62,/) 6034 FORMAT (/' Refined values of diffraction parameters derived from', + I5,' indexed spots',/' Crystal orientation is given at sp', + 'indle=0',/' Standard deviation of spot position (pixe', + 'ls) ',F8.2,/' or ',F8.2,' mm',/' Standard deviation of s', + 'pindle position (degrees)',F8.2,/' Detector origin (pixe', + 'ls) at ',2F10.2,/' Crystal to detector distance (mm) ', + F10.2,/) 6036 FORMAT (' Lab coordinates of detector X-axis ',3F10.6,/' Lab coo', + 'rdinates of detector Y-axis ',3F10.6,/' Direct beam coor', + 'dinates (rec. Angstroem) ',3F10.6,/' Coordinates of unit', + ' cell A-axis ',3F10.3,/' Coordinates of unit cell B-axis ' + ,3F10.3,/' Coordinates of unit cell C-axis ',3F10.3,/' Un', + 'it cell parameters ',3F8.2,3F8.3,/' Reciprocal cel', + 'l parameters ',3F8.6,3F8.3,/ + ' Standard deviations ',6F8.6,/) 6038 FORMAT (//' ***** Solution #',I3,' Quality ',F8.5,//' Interpreta', + 'tion of observed difference vector clusters',/' Coordina', + 'tes(rec.Angstroem) Frequency Logical Indices ', + ' rms',/' X Y Z ', + 'H K L (rec.Ang.)',/) 6040 FORMAT (I5,3F9.6,I7,1X,3F9.3,F10.6) 6042 FORMAT (/' Autoindexing is based on ',I6,' Observed spot positio', + 'ns above threshold intensity of ',F7.0) C C 160 CONTINUE END C C C BLOCK DATA REFBLK C ================= C C C C C .. Scalars in Common .. REAL CCOM,CCX,CCY,DPHI,F,FLAMDA,ORGX,ORGY,Q,THRESH INTEGER FIXF,IGROUP LOGICAL DCOMP,FILM,LCAMC,TARGET C .. C .. Arrays in Common .. REAL ACHSE,CELL,ED,S0,TARMAT C .. C .. Common blocks .. COMMON /REFCOM/IGROUP,FLAMDA,CELL(6),F,THRESH,Q,DPHI,ED(3,3), + ACHSE(3),S0(3),ORGX,ORGY,FIXF,TARMAT(3,3),TARGET,DCOMP, + LCAMC,CCX,CCY,CCOM,FILM C .. SAVE C .. Data statements .. C C---- Data values are appropriate for MOSFLM axis convention, C with X along X-ray beam, and Z along rotation axis. C C C---- For hamburg image plates on synchrotron C RASTER (Q) = 0.187, beam centres are orgx,orgy values C depending on scanner program version C C C ccc Q/0.050/ DATA IGROUP/0/,FLAMDA/1.5418/,CELL/3*0.0,3*90.0/,F/0.0/, + THRESH/0.0/,S0/1.0,0.0,0.0/,ACHSE/0.0,0.0,1.0/,ED/0.0,1.0, + 0.0,0.0,0.0,1.0,3*0.0/,FIXF/1/,DPHI/1.0/,Q/0.150/, + ORGX/-999.0/,ORGY/-999.0/ DATA TARGET/.FALSE./,DCOMP/.FALSE./,FILM/.FALSE./,LCAMC/.FALSE./ C .. C C END C C C SUBROUTINE READDAT C ================== C C C .. Parameters .. INTEGER NPARS PARAMETER (NPARS=200) C .. C .. Scalars in Common .. REAL CCOM,CCX,CCY,DPHI,F,FLAMDA,ORGX,ORGY,Q,THRESH INTEGER FIXF,IGROUP LOGICAL DCOMP,FILM,LCAMC,TARGET INTEGER LP,LINOUT LOGICAL ONLINE C .. C .. Arrays in Common .. REAL ACHSE,CELL,ED,S0,TARMAT C .. C .. Local Scalars .. REAL ORGX1,ORGX2,ORGX3,ORGY1,ORGY2,ORGY3 INTEGER I,ICOMM,IFAIL,IOUT,ITIN,ITINS,J,J1,K,NTOK, + NWRN LOGICAL COMREAD,LEND,RCELL,RDIST,RDPHI,RGROUP,RRAST,RWAVE CHARACTER K1*4,KEY*4,MATFIL*255,LINE*400 C .. C .. Local Arrays .. REAL VALUE(NPARS) INTEGER IBEG(NPARS),IDEC(NPARS),IEND(NPARS),ITYP(NPARS) CHARACTER CVALUE(NPARS)*4 C .. C .. External Subroutines .. EXTERNAL CCPOPN,CCPUPC,KEYERR,KEYNUM,PARSER C .. C .. Intrinsic Functions .. INTRINSIC NINT C .. C .. Common blocks .. COMMON /REFCOM/IGROUP,FLAMDA,CELL(6),F,THRESH,Q,DPHI,ED(3,3), + ACHSE(3),S0(3),ORGX,ORGY,FIXF,TARMAT(3,3),TARGET,DCOMP, + LCAMC,CCX,CCY,CCOM,FILM COMMON /RINOUT/ ONLINE,LP,LINOUT C .. SAVE C .. Data statements .. C C---- data statements moved to BLOCK DATA C C DATA RCELL/.FALSE./,RGROUP/.FALSE./,RDPHI/.FALSE./,RRAST/.FALSE./, + RDIST/.FALSE./ C .. C C C---- Logical FILM simply determines default pixel size C ITIN = 5 IOUT = 6 ICOMM = 7 COMREAD = .FALSE. 10 CONTINUE C IF (.NOT.COMREAD) WRITE (LP,FMT=6000) IF (.NOT.COMREAD .AND. ONLINE) WRITE (LINOUT,FMT=6000) 20 CONTINUE C C---- Read next keyword C C LINE = ' ' KEY = ' ' NTOK = NPARS C C ******************************************************* CALL PARSER(KEY,LINE,IBEG,IEND,ITYP,VALUE,CVALUE,IDEC,NTOK,LEND, + .TRUE.) C ******************************************************** C C---- End of file? C IF (LEND) THEN GO TO 340 ELSE C C---- First 4 chars C KEY = LINE(IBEG(1) :IEND(1)) C C---- Convert to upper case C C *********** CALL CCPUPC(KEY) C *********** C C---- SPAC*EGROUP C =========== C IF (KEY.EQ.'SPAC') THEN GO TO 320 C C---- WAVE*LENGTH C =========== C ELSE IF (KEY.EQ.'WAVE') THEN GO TO 310 C C---- CELL C ==== C ELSE IF (KEY.EQ.'CELL') THEN GO TO 290 C C---- DIST*ANCE C ========= C ELSE IF (KEY.EQ.'DIST') THEN GO TO 280 C C---- THRE*SHOLD C ========== C ELSE IF (KEY.EQ.'THRE') THEN GO TO 270 C C---- PIXE*L C ====== C ELSE IF (KEY.EQ.'PIXE') THEN GO TO 260 C C---- DELP*HI C ======= C ELSE IF (KEY.EQ.'DELP') THEN GO TO 250 C C---- FIXC*ELL C ======== C ELSE IF (KEY.EQ.'FIXC') THEN GO TO 200 C C---- XRAY C ==== C ELSE IF (KEY.EQ.'XRAY') THEN GO TO 180 C C---- ROTA*XIS C ======== C ELSE IF (KEY.EQ.'ROTA') THEN GO TO 160 C C---- DETM*AT C ======= C ELSE IF (KEY.EQ.'DETM') THEN GO TO 130 C C---- UMAT C ==== C ELSE IF (KEY.EQ.'UMAT') THEN GO TO 100 C C---- MATRix x x x x x x x x x C ========================= C ELSE IF (KEY.EQ.'MATR') THEN GO TO 70 C C---- CAMCON C ====== C ELSE IF (KEY.EQ.'CAMC') THEN GO TO 50 C C---- FILM or IMAGE plate C ==================== C ELSE IF (KEY.EQ.'FILM') THEN GO TO 40 ELSE IF (KEY.EQ.'IMAG') THEN GO TO 30 C C---- RUN C === C ELSE IF (KEY.EQ.'RUN ') THEN C IF (COMREAD) ITIN = ITINS C C---- check that spacegoup,cell and distance have been supplied C IF (.NOT.RGROUP) THEN IF (ONLINE) WRITE (LINOUT,FMT=6038) WRITE (LP,FMT=6038) IF (ONLINE) THEN GO TO 20 ELSE GO TO 380 END IF C C ELSE IF (.NOT.RCELL) THEN IF (ONLINE) WRITE (LINOUT,FMT=6040) WRITE (LP,FMT=6040) IF (ONLINE) THEN GO TO 20 ELSE GO TO 370 END IF C C ELSE IF (RDIST) THEN GO TO 350 ELSE IF (ONLINE) WRITE (LINOUT,FMT=6042) WRITE (LP,FMT=6042) IF (ONLINE) THEN GO TO 20 ELSE GO TO 360 END IF END IF END IF END IF C C---- not recognised C NWRN = NWRN + 1 IF (ONLINE) WRITE (LINOUT,FMT=6058) WRITE (LP,FMT=6058) C C---- end of control cards C GO TO 330 30 FILM = .FALSE. GO TO 330 40 FILM = .TRUE. GO TO 330 50 LCAMC = .TRUE. I = 1 60 CONTINUE I = I + 1 C C---- Skip if no more tokens on line C IF (I.GT.NTOK) THEN GO TO 330 ELSE K1 = LINE(IBEG(I) :IEND(I)) C C ********** CALL CCPUPC(K1) C ********** C C---- CAMCon CCX C =========== C IF (K1.EQ.'CCX') THEN I = I + 1 C C ************************************ CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C CCX = VALUE(I) C C---- CAMCon CCY C ========== C ELSE IF (K1.EQ.'CCY') THEN I = I + 1 C C ************************************ CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C CCY = VALUE(I) C C---- CAMCon CCOM C ============ C ELSE IF (K1.EQ.'CCOM') THEN I = I + 1 C C ************************************ CALL KEYNUM(1,I,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C CCOM = VALUE(I) ELSE GO TO 390 END IF C C GO TO 60 END IF 70 DCOMP = .TRUE. C C---- Check if second token is alphanumeric, if so it should C be the file containing AMAT. C IF (ITYP(2).EQ.1) THEN MATFIL = LINE(IBEG(2) :IEND(NTOK)) IF (ONLINE) WRITE (LINOUT,FMT=*) ' MATFIL',MATFIL WRITE (LP,FMT=*) ' MATFIL',MATFIL C C IFAIL = 1 C C ***************************** CALL CCPOPN(3,MATFIL,3,1,80,IFAIL) C ***************************** C READ (3,FMT=6034,END=400) ((TARMAT(I,J),J=1,3),I=1,3) CLOSE (UNIT=3) ELSE C C ************************************ CALL KEYNUM(6,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C K = 1 C C DO 90 J = 1,3 DO 80 I = 1,3 K = K + 1 TARMAT(I,J) = VALUE(K) 80 CONTINUE 90 CONTINUE C C END IF C C IF (ONLINE) WRITE (LINOUT,FMT=6036) + ((TARMAT(J1,J),J1=1,3),J=1,3) WRITE (LP,FMT=6036) + ((TARMAT(J1,J),J1=1,3),J=1,3) TARGET = .TRUE. GO TO 330 C C ************************************ 100 CALL KEYNUM(6,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C K = 1 C C DO 120 J = 1,3 DO 110 I = 1,3 K = K + 1 TARMAT(J,I) = VALUE(K) 110 CONTINUE 120 CONTINUE C C IF (ONLINE) WRITE (LINOUT,FMT=6032) + ((TARMAT(J1,J),J1=1,3),J=1,3) WRITE (LP,FMT=6032) ((TARMAT(J1,J),J1=1,3),J=1,3) TARGET = .TRUE. GO TO 330 C C ************************************ 130 CALL KEYNUM(6,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C K = 1 C C DO 150 J = 1,2 DO 140 I = 1,3 K = K + 1 ED(I,J) = VALUE(K) 140 CONTINUE 150 CONTINUE C C IF (ONLINE) WRITE (LINOUT,FMT=6030) + ((ED(J1,J),J1=1,3),J=1,2) WRITE (LP,FMT=6030) ((ED(J1,J),J1=1,3),J=1,2) GO TO 330 C C ************************************ 160 CALL KEYNUM(3,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C DO 170 I = 1,3 ACHSE(I) = VALUE(I+1) 170 CONTINUE C C IF (ONLINE) WRITE (LINOUT,FMT=6028) ACHSE WRITE (LP,FMT=6028) ACHSE GO TO 330 C C ************************************ 180 CALL KEYNUM(3,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C DO 190 I = 1,3 S0(I) = VALUE(I+1) 190 CONTINUE C C IF (ONLINE) WRITE (LINOUT,FMT=6026) S0 WRITE (LP,FMT=6026) S0 GO TO 330 200 FIXF = 0 IF (ONLINE) WRITE (LINOUT,FMT=6024) WRITE (LP,FMT=6024) GO TO 330 C C ************************************ 250 CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C DPHI = VALUE(2) IF (ONLINE) WRITE (LINOUT,FMT=6014) DPHI WRITE (LP,FMT=6014) DPHI RDPHI = .TRUE. GO TO 330 C C ************************************ 260 CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C Q = VALUE(2) IF (ONLINE) WRITE (LINOUT,FMT=6012) Q WRITE (LP,FMT=6012) Q RRAST = .TRUE. GO TO 330 C C ************************************ 270 CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C THRESH = VALUE(2) IF (ONLINE) WRITE (LINOUT,FMT=6010) THRESH WRITE (LP,FMT=6010) THRESH GO TO 330 C C ************************************ 280 CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C F = VALUE(2) IF (ONLINE) WRITE (LINOUT,FMT=6008) F WRITE (LP,FMT=6008) F RDIST = .TRUE. GO TO 330 C C ************************************ 290 CALL KEYNUM(3,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C CELL(1) = VALUE(2) CELL(2) = VALUE(3) CELL(3) = VALUE(4) C C IF (NTOK.GT.4) THEN C C ****************************************** CALL KEYNUM((NTOK-4),5,LINE,IBEG,IEND,ITYP,NTOK) C ****************************************** C DO 300 I = 4,NTOK - 1 CELL(I) = VALUE(I+1) 300 CONTINUE END IF C C IF (ONLINE) WRITE (LINOUT,FMT=6006) CELL WRITE (LP,FMT=6006) CELL RCELL = .TRUE. GO TO 330 C C ************************************ 310 CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C FLAMDA = VALUE(2) IF (ONLINE) WRITE (LINOUT,FMT=6004) FLAMDA WRITE (LP,FMT=6004) FLAMDA RWAVE = .TRUE. GO TO 330 C C ************************************ 320 CALL KEYNUM(1,2,LINE,IBEG,IEND,ITYP,NTOK) C ************************************ C IGROUP = NINT(VALUE(2)) IF (ONLINE) WRITE (LINOUT,FMT=6002) IGROUP WRITE (LP,FMT=6002) IGROUP RGROUP = .TRUE. C C---- Read next control card C 330 CONTINUE GO TO 10 340 IF (COMREAD) THEN COMREAD = .FALSE. ITIN = ITINS CLOSE (UNIT=ICOMM) GO TO 10 ELSE GO TO 410 END IF C C 350 IF (.NOT.RWAVE .AND. ONLINE) WRITE (LINOUT,FMT=6044) FLAMDA IF (.NOT.RDPHI .AND. ONLINE) WRITE (LINOUT,FMT=6046) DPHI IF (.NOT.RWAVE ) WRITE (LP,FMT=6044) FLAMDA IF (.NOT.RDPHI ) WRITE (LP,FMT=6046) DPHI C C---- Set default pixel size for film (default raster is 0.150) C IF (.NOT.RRAST .AND. FILM) Q = 0.05 C C---- Check pixel size C IF (FILM .AND. (Q.NE.0.05)) THEN IF (ONLINE) WRITE (LINOUT,FMT=6048) Q WRITE (LP,FMT=6048) Q ELSE IF ((.NOT.FILM) .AND. (Q.NE.0.150)) THEN IF (ONLINE) WRITE (LINOUT,FMT=6050) Q WRITE (LP,FMT=6050) Q END IF RETURN 360 STOP 370 STOP 380 STOP C C---- CAMCon ?????? C ============= C 390 CONTINUE C C ******************************* CALL KEYERR(I,0,LINE,IBEG,IEND,ITYP) C ******************************* C STOP C C---- eof on input 'matrix' file C 400 IF (ONLINE) WRITE (LINOUT,FMT=6060) MATFIL WRITE (LP,FMT=6060) MATFIL C C---- eof on input C 410 IF (ONLINE) WRITE (LINOUT,FMT=6062) WRITE (LP,FMT=6062) STOP C C---- Format statements C 6000 FORMAT (1X,'REFIX=> ',$) 6002 FORMAT (/1X,'Spacegroup number',I5) 6004 FORMAT (/1X,'WAVELENGTH',F7.4) 6006 FORMAT (/1X,'Unit cell parameters ',6F8.2) 6008 FORMAT (/1X,'Crystal to detector distance ',F7.2,'mm') 6010 FORMAT (/1X,'Intensity threshold for spot selection ',F8.1) 6012 FORMAT (/1X,'Pixel size ',F6.3,'mm') 6014 FORMAT (/1X,'Expected phi error ',F7.2,' degrees') 6024 FORMAT (/1X,'Cell parameters will be fixed, crystal to detector ', + 'distance will be refined') 6026 FORMAT (/1X,'Direct beam vector ',3F7.1) 6028 FORMAT (/1X,'X,Y,Z Components of rotation axis ',3F7.4) 6030 FORMAT (/1X,'Detector orientation matrix',/1X,'values of ED(1,1)', + ',ED(2,1),ED(3,1),ED(1,2),ED(2,2),ED(3,2)',/1X,6F8.5) 6032 FORMAT (/1X,'TARGET U MATRIX',3 (/1X,3F8.5)) 6034 FORMAT (3F12.6) 6036 FORMAT (/1X,'A matrix read (from which target u matrix will be c', + 'alculated)',/3 (1X,3F8.5,/)) 6038 FORMAT (//1X,'*** Spacegroup not given ***') 6040 FORMAT (//1X,'*** Unit cell not given ***') 6042 FORMAT (//1X,'*** Detector distance not given ***') 6044 FORMAT (/1X,'Wavelength (Default) is ',F7.4) 6046 FORMAT (/1X,'Phi error (Default) is ',F4.2) 6048 FORMAT (//1X,'****** WARNING *******',/2X,'Program expects film ', + 'data but pixel size is',F5.3,' is this correct ?') 6050 FORMAT (//1X,'****** WARNING *******',/2X,'Program expects image', + ' plate data but pixel size is ',F5.3,'mm') 6058 FORMAT (2X,'****** Keyword not recognised ***** ') 6060 FORMAT (//1X,'******* End of file when reading matrixfrom ',A, + '*******') 6062 FORMAT (//1X,'*** Warning ** NO "RUN" card given') C C END C C C SUBROUTINE ABSHKL(IGROUP,CELL,ABSTYP) C ===================================== C C C C C---- Obtain constraint number for general systematic absences C from space-group number and unit cell parameters C C VERSION 10-1986 C C---- Description of calling parameters C C IGROUP - space-group number from international tables (given) C CELL - unit cell parameter (given) C ABSTYP - constraint number for general conditions limiting (result) C possible reflections C 0: no conditions on h,k,l C 1: h+k must be even C 2: h+l must be even C 3: k+l must be even C 4: h+k+l must be even C 5: -h+k+l must be a multiple of three C 6: h+k and h+l must be even C C .. Scalar Arguments .. INTEGER ABSTYP,IGROUP C .. C .. Array Arguments .. REAL CELL(6) C .. C .. Local Scalars .. INTEGER G,J C .. C .. Local Arrays .. INTEGER*2 IG(86),NUM(86) C .. SAVE C .. Data statements .. C C---- Initialize look-up tables of pairs (igroup,abstyp) C DATA IG/517,520,521,524,527,20,21,35,36,37,63,64,65,66,67,68,5,8, + 9,12,15,38,39,40,41,23,24,44,45,46,71,72,73,74,79,80,82,87, + 88,97,98,107,108,109,110,119,120,121,122,139,140,141,142,197, + 199,204,206,211,214,217,220,229,230,146,148,155,160,161,166, + 167,22,42,43,69,70,196,202,203,209,210,216,219,225,226,227, + 228/ DATA NUM/16*1,5*2,4*3,38*4,7*5,16*6/ C .. C C G = IGROUP IF ((IGROUP.GT.4) .AND. (IGROUP.LT.16) .AND. + (CELL(5).NE.90.0)) G = IGROUP + 512 IF ((IGROUP.GT.145) .AND. (IGROUP.LT.168) .AND. + (CELL(1).EQ.CELL(2)) .AND. (CELL(2).EQ.CELL(3)) .AND. + (CELL(4).EQ.CELL(5)) .AND. (CELL(5).EQ.CELL(6))) G = IGROUP + + 512 ABSTYP = 0 C C DO 10 J = 1,86 IF (G.EQ.IG(J)) ABSTYP = NUM(J) 10 CONTINUE C C END C C C SUBROUTINE AUTOIX(IGROUP,CELL,FIXF,Q,ORGX,ORGY,F,ED,AXIS,S0,RCELL, + U,SDU,SDXY,SDPHI,SDCELL,MXBEST,NBEST,ORGXB, + ORGYB,FB,EDB,S0B,RCELLB,UB,ESDU,ESDXY,ESDPHI, + ESDCEL,ESD,RMS,NACCB,MXSPOT,INTY,VC,FHKL,NSPOT, + MXRFL,NRFL,IX,IY,IPHI,IH,IK,IL,V,MXY,TABLE,NVEC, + IER) C ================================================================= C C C---- Autoindexing algorithm (w.kabsch 10-1986) C C REVISION 11-1987 ( DEVICE INDEPENDENT ) C C This routine attempts to find an initial value C for the orientation matrix connecting the crystal C system with the laboratory coordinate system C and refines the relevant parameters. C C C C---- Description of calling parameters **** C C IGROUP - space group number (given) C CELL - unit cell parameter (given) C FIXF - refinement control parameter (given) C -1: refine detector orientation and distance C 0: refine direct beam direction, unit cell C orientation and detector distance C 1: refine direct beam direction, unit cell C orientation and independent cell parameters C Q - length of a pixel in millimeters (given) C ORGX - x- and y-coordinates (pixels) on detector such (updated) C ORGY - that the detector normal would intersect the (updated) C crystal C F - crystal to detector distance (mm). (updated) C ED - real array(3,3) specifying lab coordinates of (updated) C detector orientation C AXIS - array(3) specifying lab coordinates of rotation axis (given) C S0 - coordinates of incident x-ray beam wave vector. (updated) C only the direction of s0 is refined. C length of s is 1.0/lambda. ( reciprocal angstroem ) C RCELL - array(6) of reciprocal unit cell parameters (updated) C in reciprocal angstroem and degrees. C U - array(3,3) containing orientation matrix (result) C SDU - estimated standard deviation (degrees) for (result) C orientation matrix u C SDXY - standard deviation of spot position on detector (result) C (pixels). a positive input value is C used to calculate the weights of the observed C spot positions of the reflections. C a non-positive value indicates that a default C value of 3 pixels should be used. C on return, sdxy contains the estimated standard C deviation of the spot positions on the detector. C SDPHI - standard deviation of angular position of spindle (result) C ( degrees ). a positive input value is C used to calculate the weights of the observed C spindle positions of the reflections. C a non-positive value indicates that a default C value of 0.1 degrees should be used. C on return, sdphi contains the estimated standard C deviation of the spindle angle at which the C reflections were diffracting. C SDCELL - array(6) of standard deviations for reciprocal (result) C unit cell parameters (rec. angstroem and degrees) C MXBEST - maximum number of solutions to be remembered (given) C NBEST - actual number of independent solutions found (result) C ORGXB - array(mxbest) of x-coordinates of origin, (result) C ORGYB - array(mxbest) of y-coordinates of origin (result) C FB - array(mxbest) of crystal to detector distances (result) C EDB - real array(3,3,mxbest) specifying lab coordinates (result) C of detector orientation C S0B - array(3,mxbest) of refined direct beam coordinates (result) C RCELLB - array(6,mxbest) of refined recipr. unit cell param. (result) C UB - array(3,3,mxbest) of refined orientation matricesm. (result) C ESDU - array(mxbest) estimated standard deviations of ub (result) C ESDXY - array(mxbest) estimated standard deviations of spot (result) C positions (pixels). C ESDPHI - array(mxbest) of estimated standard deviations of (result) C spindle positions of spots. C ESDCEL - array(6,mxbest) of estimated standard deviations of (result) C refined reciprocal unit cell parameters. C ESD - array(mxbest) specifying quality of solutions as (result) C obtained from subroutine "orient" C RMS - array(mxbest) specifying rms-discrepancy between (result) C expected and refined unit cell parameters C NACCB - integer array(mxbest) specifying number of indexed (result) C spots for each solution. each number means: C >0: number of indexed spots. no error. C 0: insufficient number of accepted reflections by C subroutines "refine" or "indix" C -1: q is not positive ( returned by "refine" or "indix") C -2: f is zero ( returned by "refine" or "indix") C -3: illegal cell param.( returned by "refine" or "indix") C -4: illegal direct beam( returned by "refine" or "indix") C -5: orientation matrix is singular (returned by "refine") C -6: ncycle not between 1..10 (returned by "refine") C -7: error from "gels" (returned by "relix" ) C -8: error using "gels" (returned by "relix" ) C -9: nv<2 or illegal setting matrix (returned by "relix" ) C -10: error in "u3best" (returned by "relix" ) C MXSPOT - maximum number of difference vector clusters (given) C defines dimension of arrays inty,vc,fhkl C INTY - integer*2 array(mxspot) containing population of (result) C difference vector clusters C VC - real array(3,mxspot) of difference vector coordinates(result) C FHKL - real array(4,mxspot,mxbest) containing (result) C real valued indices determined by the local C indexing method and rms-deviation with respect C to the original reciprocal unit cell geometry. C NSPOT - number of difference vector clusters found (result) C MXRFL - dimension of arrays ix,iy,iphi,ih,ik,il (given) C NRFL - number of reflections to be indexed (given) C IX - integer*2 array specifying x-position on detector (modified) C for each spot in the list ( tenth of a pixel) C IY - integer*2 array specifying y-position on detector (modified) C for each spot in the list ( tenth of a pixel) C IPHI - integer*2 array specifying angular position of (modified) C spindle. the values are given in units of a C hundreth of a degree. iphi(jrfl) refers to the C spot located at ix(jrfl),iy(jrfl). C IH - integer*2 array specifying h-indices for the spots (result) C located at ix,iy,iphi C IK - integer*2 array specifying k-indices for the spots (result) C located at ix,iy,iphi. C IL - integer*2 array specifying l-indices for the spots (result) C located at ix,iy,iphi. C V - real array(3,mxrfl) contains coordinates of recipr. (result) C lattice points corresponding to reflections located C at ix,iy,iphi C MXY - specifies dimensions of array table (given) C TABLE - integer*2 array(mxy,mxy) used for collecting (scratch) C difference vectors C NVEC - integer*4 value specifying number of difference (result) C vectors between reflections found in clusters. C IER - number specifying results obtained by this routine (result) C 1 : no error. solution is unique. C 2 : warning. solution may not be unique. C 3 : error. best solution cannot explain data. C 4 : error. value of nbest is zero. C 5 : error. insufficient number of accepted spots. C 6 : error. dimension of difference vector set <2. C 7 : error. illegal space group number. C 8 : error. illegal cell parameters. C 9 : error. illegal incident beam direction or rotation axis C 10 : error. crystal to detector distance is zero. C 11 : error. undefined detector orientation or size of a C detector pixel <=0. C C---- External subroutines reqired C C ABSHKL,CRYSYS,DIFVEC,INDIX,INVCEL,INVERS,MATMUL, C ORIENT,PXTOMM,RQSORT,RQSORT2,REFINE,RELIX,SETMAT, C SETTING C C C C C---- check precondition C C .. Parameters .. INTEGER NBXY,NBZ,MZ,MZSTEP,NCYCLE,NRFL0*4,NRFLX,NNRFLX,NRFLX3 PARAMETER (NBXY=4,NBZ=3,MZ=8,MZSTEP=2,NCYCLE=4,NRFL0=0,NRFLX=120, + NNRFLX= (NRFLX*NRFLX+NRFLX)/2,NRFLX3=3*NRFLX) REAL RAD,RAD100,SNR,DELTA PARAMETER (RAD=57.29578,RAD100=100.0*RAD,SNR=3.0,DELTA=3.0) C .. C .. Scalar Arguments .. REAL F,ORGX,ORGY,Q,SDPHI,SDU,SDXY INTEGER FIXF,IER,IGROUP,MXBEST,MXRFL,MXSPOT,MXY,NBEST,NRFL,NSPOT INTEGER NVEC C .. C .. Array Arguments .. REAL AXIS(3),CELL(6),ED(3,3),EDB(3,3,MXBEST),ESD(MXBEST), + ESDCEL(6,MXBEST),ESDPHI(MXBEST),ESDU(MXBEST),ESDXY(MXBEST), + FB(MXBEST),FHKL(4,MXSPOT,MXBEST),ORGXB(MXBEST),ORGYB(MXBEST), + RCELL(6),RCELLB(6,MXBEST),RMS(MXBEST),S0(3),S0B(3,MXBEST), + SDCELL(6),U(3,3),UB(3,3,MXBEST),V(3,MXRFL),VC(3,MXSPOT) INTEGER NACCB(MXBEST) INTEGER*2 IH(*),IK(*),IL(*),INTY(MXSPOT),IPHI(*),IX(*),IY(*), + TABLE(MXY,MXY) C .. C .. Local Scalars .. REAL CPHI,GRID,R,RRMIN,RX,RY,SPHI,SS,X,Y,Z INTEGER ABSTYP,IBEST,ICS,ICS0,J,JRFL,K,L,M,MRFL,NACC,NFAR,NOVLAP, + NV C .. C .. Local Arrays .. REAL A(3,3),AA(3,3),ACHSE(3),B(3,3),COEF(NNRFLX),RIGHT(NRFLX3), + T(3,NRFLX),WEIGHT(NRFLX),XCALC(3,NRFLX),XOBS(3,NRFLX) C .. C .. External Subroutines .. EXTERNAL ABSHKL,CRYSYS,DIFVEC,INDIX,INVCEL,INVERS,MATMUL,ORIENT, + PXTOMM,REFINE,RELIX,RQSORT,RQSORT2,SETMAT,SETTING,UNORM C .. C .. Intrinsic Functions .. INTRINSIC COS,MAX,MIN,SIN,SQRT C .. SAVE C C IF ((MXRFL.LT.MXSPOT) .OR. (NRFL.GT.MXRFL) .OR. + (MXBEST.LT.1)) THEN STOP ' Programming error in using "AUTOIX" ' ELSE C C---- get constraint number of general systematic absences C C ************************** CALL ABSHKL(IGROUP,CELL,ABSTYP) C ************************** C C---- get cystal system number C C *********************** CALL CRYSYS(IGROUP,CELL,ICS) C *********************** C IF (ICS.GE.1) THEN C C---- get reciprocal cell C C ******************** CALL INVCEL(CELL,RCELL,X) C ******************** C IF (X.GT.0.0) THEN C C---- check calling parameters C IF (NRFL.GE.25) THEN C C *********** CALL UNORM(ED,J) C *********** C IF ((J.NE.0) .OR. (Q.LE.0.0)) THEN IER = 11 GO TO 340 ELSE IF (F.EQ.0.0) THEN IER = 10 GO TO 340 ELSE SS = S0(1)**2 + S0(2)**2 + S0(3)**2 IF (SS.GT.0.0) THEN R = SQRT(AXIS(1)**2+AXIS(2)**2+AXIS(3)**2) IF (R.GT.0.0) THEN ACHSE(1) = AXIS(1)/R ACHSE(2) = AXIS(2)/R ACHSE(3) = AXIS(3)/R C C---- calculate setting matrix of reciprocal cell C C *************** CALL SETMAT(RCELL,A) C *************** C C---- calculate setting matrix of unit cell C C ************* CALL INVERS(A,B,R) C ************* C IF (R.LE.0.0) THEN GO TO 330 ELSE C C---- set starting values of parameters C DO 40 IBEST = 1,MXBEST NACCB(IBEST) = 0 ESDPHI(IBEST) = SDPHI ESDXY(IBEST) = SDXY ESDU(IBEST) = 0.0 ORGXB(IBEST) = ORGX ORGYB(IBEST) = ORGY FB(IBEST) = F RMS(IBEST) = 1000.0 C C DO 20 J = 1,3 S0B(J,IBEST) = S0(J) C C DO 10 K = 1,3 EDB(J,K,IBEST) = ED(J,K) 10 CONTINUE 20 CONTINUE C C DO 30 J = 1,6 ESDCEL(J,IBEST) = 0.0 RCELLB(J,IBEST) = RCELL(J) 30 CONTINUE 40 CONTINUE C C---- calculate transformation from pixel to MM C C ***************** CALL PXTOMM(Q,F,ED,AA) C ***************** C C---- calculate pixel coordinates of direct beam C X = ED(1,1)*S0(1) + ED(2,1)*S0(2) + ED(3,1)*S0(3) Y = ED(1,2)*S0(1) + ED(2,2)*S0(2) + ED(3,2)*S0(3) Z = ED(1,3)*S0(1) + ED(2,3)*S0(2) + ED(3,3)*S0(3) C C IF (Z.NE.0.0) THEN R = F/ (Q*Z) X = R*X + ORGX Y = R*Y + ORGY C C---- sort spots in increasing resolution C C *************************** CALL RQSORT(NRFL,IX,IY,IPHI,X,Y) C *************************** C C---- calculate laboratory coordinates of observed spots at phi=0 C DO 50 NV = 1,NRFL R = IPHI(NV)/RAD100 CPHI = COS(R) SPHI = SIN(R) RX = IX(NV)/10.0 - ORGX RY = IY(NV)/10.0 - ORGY X = AA(1,1)*RX + AA(1,2)*RY + AA(1,3) Y = AA(2,1)*RX + AA(2,2)*RY + AA(2,3) Z = AA(3,1)*RX + AA(3,2)*RY + AA(3,3) R = SQRT((X*X+Y*Y+Z*Z)/SS) IF (R.GT.0.0) THEN X = X/R - S0(1) Y = Y/R - S0(2) Z = Z/R - S0(3) R = (ACHSE(1)*X+ACHSE(2)*Y+ACHSE(3)*Z)* + (1.0-CPHI) END IF V(1,NV) = ACHSE(1)*R + X*CPHI - + (ACHSE(2)*Z-ACHSE(3)*Y)*SPHI V(2,NV) = ACHSE(2)*R + Y*CPHI - + (ACHSE(3)*X-ACHSE(1)*Z)*SPHI V(3,NV) = ACHSE(3)*R + Z*CPHI - + (ACHSE(1)*Y-ACHSE(2)*X)*SPHI 50 CONTINUE C C MRFL = MIN(NRFL,600) IF (MRFL.LT.25) THEN GO TO 320 ELSE C C---- find most frequent differences of laboratory coordinates C of observed spots. C GRID = MIN(RCELL(1),RCELL(2),RCELL(3))*0.1 C C **************************************** CALL DIFVEC(MRFL,V,DELTA,SNR,GRID,NBXY,NBZ, + MXY,MZ,MZSTEP,MXSPOT,TABLE,IH,IK, + IL,INTY,VC,NSPOT,NVEC) C **************************************** C C---- sort list of difference vectors in decreasing population C C ********************** CALL RQSORT2(NSPOT,VC,INTY) C ********************** C IF (NSPOT.GT.10) NSPOT = 10 C C---- get possible orientations of unit cell C C **************************************** CALL ORIENT(IGROUP,CELL,NSPOT,VC,MXBEST,ESD, + UB,NBEST) C **************************************** C IF (NBEST.LE.0) THEN C C IF ((NBEST+1).LT.0) THEN IER = 6 GO TO 340 ELSE IF ((NBEST+1).EQ.0) THEN GO TO 330 ELSE IF ((NBEST+1).GT.0) THEN IER = 4 GO TO 340 END IF END IF C C---- refine all parameters for each of the nbest solutions C RRMIN = (MAX(RCELL(1),RCELL(2),RCELL(3))* + 2.5)**2 NV = MIN(NRFLX,MRFL) C C DO 170 IBEST = 1,NBEST ORGX = ORGXB(IBEST) ORGY = ORGYB(IBEST) F = FB(IBEST) SDPHI = ESDPHI(IBEST) SDXY = ESDPHI(IBEST) SDU = ESDU(IBEST) C C DO 60 J = 1,6 RCELL(J) = RCELLB(J,IBEST) SDCELL(J) = ESDCEL(J,IBEST) 60 CONTINUE C C DO 80 J = 1,3 S0(J) = S0B(J,IBEST) C C DO 70 K = 1,3 ED(J,K) = EDB(J,K,IBEST) U(J,K) = UB(J,K,IBEST) 70 CONTINUE 80 CONTINUE C C ************************************ CALL RELIX(ABSTYP,RRMIN,A,NSPOT,VC,XOBS, + XCALC,WEIGHT,IH,IK,IL,T,RIGHT, + COEF,U,IER) C ************************************ C IF (IER.GE.0) THEN IF (IER.LE.0) THEN C C ************** CALL MATMUL(U,A,AA) C ************** C DO 100 JRFL = 1,NSPOT FHKL(1,JRFL,IBEST) = RIGHT(JRFL) FHKL(2,JRFL,IBEST) = RIGHT(JRFL+NSPOT) FHKL(3,JRFL,IBEST) = RIGHT(JRFL+NSPOT+ + NSPOT) X = IH(JRFL) Y = IK(JRFL) Z = IL(JRFL) R = 0.0 C C DO 90 K = 1,3 R = (VC(K,JRFL)-AA(K,1)*X-AA(K,2)*Y- + AA(K,3)*Z)**2 + R 90 CONTINUE C C FHKL(4,JRFL,IBEST) = SQRT(R) 100 CONTINUE C C ******************************* CALL RELIX(ABSTYP,RRMIN,A,NV,V,XOBS, + XCALC,WEIGHT,IH,IK,IL,T, + RIGHT,COEF,U,IER) C ******************************* C JRFL = NV IF (IER.LT.0) THEN GO TO 120 ELSE IF (IER.LE.0) THEN C C---- Refine origin and possibly f C IF (FIXF.GT.0) ICS0 = -3 IF (FIXF.LE.0) ICS0 = -2 C C ********************************* CALL REFINE(NCYCLE,JRFL,IH,IK,IL,IX, + IY,IPHI,ACHSE,ICS0,Q,ORGX, + ORGY,F,S0,ED,U,RCELL,SDU, + SDCELL,SDPHI,SDXY,IER) C ********************************* C IF (IER.GT.0) THEN 110 CONTINUE C C---- Refine all parameters and extend indexing C C C IF (FIXF.GT.0) THEN ICS0 = ICS ELSE ICS0 = FIXF END IF C C ******************************* CALL REFINE(NCYCLE,JRFL,IH,IK,IL,IX, + IY,IPHI,ACHSE,ICS0,Q, + ORGX,ORGY,F,S0,ED,U, + RCELL,SDU,SDCELL,SDPHI, + SDXY,IER) C ******************************* C IF (IER.GT.0) THEN J = JRFL JRFL = MIN(JRFL+JRFL,NRFL) C C ***************************** CALL INDIX(ABSTYP,Q,ORGX,ORGY,F, + ED,ACHSE,S0,RCELL,U, + MXRFL,NRFL0,JRFL,IX,IY, + IPHI,IH,IK,IL,NACC, + NOVLAP,NFAR,SDPHI,SDXY) C ***************************** C IF ((NACC.GT.0) .AND. + (J.LT.JRFL)) THEN GO TO 110 ELSE GO TO 130 END IF END IF END IF NACC = IER GO TO 130 END IF END IF NACC = -7 GO TO 130 END IF 120 NACC = IER - 7 130 NACCB(IBEST) = NACC IF (NACC.GT.0) THEN ESDU(IBEST) = SDU ESDXY(IBEST) = SDXY ESDPHI(IBEST) = SDPHI ORGXB(IBEST) = ORGX ORGYB(IBEST) = ORGY FB(IBEST) = F C C **************** CALL SETMAT(RCELL,AA) C **************** C R = 0.0 C C DO 150 J = 1,3 S0B(J,IBEST) = S0(J) C C DO 140 K = 1,3 EDB(J,K,IBEST) = ED(J,K) UB(J,K,IBEST) = U(J,K) R = (A(J,K)-AA(J,K))**2 + R 140 CONTINUE 150 CONTINUE C C RMS(IBEST) = SQRT(R) C C DO 160 J = 1,6 ESDCEL(J,IBEST) = SDCELL(J) RCELLB(J,IBEST) = RCELL(J) 160 CONTINUE END IF 170 CONTINUE C C C C---- Sort solutions in decreasing quality C IF (NBEST.GE.2) THEN NACC = -1 SDXY = 1.0E06 SDPHI = 1.0E06 C C DO 180 J = 1,NBEST IF (NACCB(J).GT.0) THEN IF (NACCB(J).GT.NACC) NACC = NACCB(J) IF ((ESDXY(J).GE.0.0) .AND. + (ESDXY(J).LT.SDXY)) SDXY = ESDXY(J) IF ((ESDPHI(J).GE.0.0) .AND. + (ESDPHI(J).LT.SDPHI)) + SDPHI = ESDPHI(J) END IF 180 CONTINUE C C IF (NACC.GT.0) THEN 190 CONTINUE C C DO 210 J = 1,NBEST - 1 X = ((NACCB(J)-NACC)/ (NACC+1.0))**2 + + ((ESDXY(J)-SDXY)/ (SDXY+0.01))**2 + + ((ESDPHI(J)-SDPHI)/ (SDPHI+0.01))**2 C C DO 200 K = J + 1,NBEST Y = ((NACCB(K)-NACC)/ (NACC+1.0))**2 + + ((ESDXY(K)-SDXY)/ (SDXY+0.01))** + 2 + ((ESDPHI(K)-SDPHI)/ + (SDPHI+0.01))**2 IF (X.GT.Y) GO TO 220 200 CONTINUE C C 210 CONTINUE GO TO 280 220 L = NACCB(J) NACCB(J) = NACCB(K) NACCB(K) = L Z = ESDXY(J) ESDXY(J) = ESDXY(K) ESDXY(K) = Z Z = ESDPHI(J) ESDPHI(J) = ESDPHI(K) ESDPHI(K) = Z Z = ESDU(J) ESDU(J) = ESDU(K) ESDU(K) = Z Z = ORGXB(J) ORGXB(J) = ORGXB(K) ORGXB(K) = Z Z = ORGYB(J) ORGYB(J) = ORGYB(K) ORGYB(K) = Z Z = FB(J) FB(J) = FB(K) FB(K) = Z Z = RMS(J) RMS(J) = RMS(K) RMS(K) = Z C C DO 240 L = 1,3 Z = S0B(L,J) S0B(L,J) = S0B(L,K) S0B(L,K) = Z C C DO 230 M = 1,3 Z = EDB(L,M,J) EDB(L,M,J) = EDB(L,M,K) EDB(L,M,K) = Z Z = UB(L,M,J) UB(L,M,J) = UB(L,M,K) UB(L,M,K) = Z 230 CONTINUE 240 CONTINUE C C DO 250 L = 1,6 Z = ESDCEL(L,J) ESDCEL(L,J) = ESDCEL(L,K) ESDCEL(L,K) = Z Z = RCELLB(L,J) RCELLB(L,J) = RCELLB(L,K) RCELLB(L,K) = Z 250 CONTINUE C C Z = ESD(J) ESD(J) = ESD(K) ESD(K) = Z C C DO 270 JRFL = 1,NSPOT DO 260 L = 1,4 Z = FHKL(L,JRFL,J) FHKL(L,JRFL,J) = FHKL(L,JRFL,K) FHKL(L,JRFL,K) = Z 260 CONTINUE 270 CONTINUE C C GO TO 190 C C---- Assign best values for parameters orgx,orgy,f,rcell,u,sdxy,ETC. C 280 IBEST = 1 ORGX = ORGXB(IBEST) ORGY = ORGYB(IBEST) F = FB(IBEST) SDPHI = ESDPHI(IBEST) SDXY = ESDXY(IBEST) SDU = ESDU(IBEST) C C DO 290 J = 1,6 RCELL(J) = RCELLB(J,IBEST) SDCELL(J) = ESDCEL(J,IBEST) 290 CONTINUE C C DO 310 J = 1,3 S0(J) = S0B(J,IBEST) C C DO 300 K = 1,3 ED(J,K) = EDB(J,K,IBEST) U(J,K) = UB(J,K,IBEST) 300 CONTINUE 310 CONTINUE C C---- Find best setting for the best solution C C ************************************ CALL SETTING(ICS,ABSTYP,A,RCELL,SDCELL,U, + R) C ************************************ C RMS(IBEST) = R END IF END IF C C---- check if solution is unique C IF (NACCB(1).LT. (0.75*NRFL)) THEN IER = 3 ELSE IF (NBEST.NE.1) THEN IF (ESDXY(2).LE. (2.0*ESDXY(1))) THEN IF (ESDPHI(2).LT. (2.0*ESDPHI(1))) THEN IER = 2 GO TO 340 END IF END IF END IF C C---- set result number and exit C IER = 1 END IF GO TO 340 END IF END IF END IF END IF END IF IER = 9 GO TO 340 END IF END IF 320 IER = 5 GO TO 340 END IF END IF 330 IER = 8 340 RETURN END IF END C C C SUBROUTINE BASIDX(ABSTYP,NEQL,MLAUE,RCELL,BASIS,IRANK,MXBEST, + NBEST,ESD,UB) C ====================================================== C C C---- Returns the "nbest" explanations of the observed vector base C in terms of rotations of the reciprocal unit cell. C C W.KABSCH VERSION 3-1987 revised 2-1988 C C---- external subroutines required C C GENABS,INVCEL,SETMAT,U3BEST,ULIST C C---- Description of calling parameters C C ABSTYP - an integer value describing conditions limiting (given) C possible reflections C 0: no conditions on H,K,L C 1: H+K must be even C 2: H+L must be even C 3: H+K+L must be even C 4: -H+K+L must be a multiple of 3 C 5: H+K AND H+L must be even C NEQL - number of laue operators (given) C MLAUE - laue operators (given) C RCELL - reciprocal cell parameters (given) C BASIS - contains observed vector base ( stored columnwise ) (given) C IRANK - number of vectors given ( at present only 2) (given) C MXBEST - maximum number of best unit cell orientations (given) C NBEST - actual number of unit cell orientations obtained (result) C nbest=0 indicates an error in the range of input C parameters C ESD - array(mxbest) describing quality of explanation (result) C UB - array(3,3,mxbest) describing orientations (result) C C C C C C .. Parameters .. INTEGER NMAX PARAMETER (NMAX=50) C .. C .. Scalar Arguments .. INTEGER ABSTYP,IRANK,MXBEST,NBEST,NEQL C .. C .. Array Arguments .. REAL BASIS(3,3),ESD(MXBEST),RCELL(6),UB(3,3,MXBEST) INTEGER MLAUE(216) C .. C .. Local Scalars .. REAL D,D1,D2,R,R1,RMS INTEGER I,IH,IHMAX,II,IK,IKMAX,IL,ILMAX,J,L,M,N1,N2 C .. C .. Local Arrays .. REAL A0(3,3),CELL(6),DELTA(NMAX,2),T(3),U(3,3),V(3,NMAX,2),W(3), + XCALC(3,3),XOBS(3,3) INTEGER N(2) C .. C .. External Subroutines .. EXTERNAL GENABS,INVCEL,SETMAT,U3BEST,ULIST C .. C .. Intrinsic Functions .. INTRINSIC ABS,NINT,SQRT C .. SAVE C C NBEST = 0 IF ((NEQL.GE.1) .AND. (NEQL.LE.24) .AND. (MXBEST.GE.1)) THEN IF (IRANK.EQ.2) THEN C C---- get unit cell parameters from reciprocal cell C C ******************** CALL INVCEL(RCELL,CELL,R) C ******************** C IF (R.GT.0.0) THEN C C---- calculate setting matrix from reciprocal unit cell C C **************** CALL SETMAT(RCELL,A0) C **************** C C---- generate lists of most plausible reciprocal lattice vectors C corresponding to the observed reciprocal vector base C DO 10 I = 1,3 W(I) = 1.0 XOBS(I,1) = BASIS(I,1) XOBS(I,2) = BASIS(I,2) XOBS(I,3) = 0.0 XCALC(I,3) = 0.0 10 CONTINUE C C DO 100 L = 1,IRANK N(L) = 0 C C---- calculate maximum values for h,k,l C R1 = SQRT(BASIS(1,L)**2+BASIS(2,L)**2+BASIS(3,L)**2) IF (R1.LE.0.0) THEN GO TO 170 ELSE IHMAX = NINT(CELL(1)*R1+0.5) IKMAX = NINT(CELL(2)*R1+0.5) ILMAX = NINT(CELL(3)*R1+0.5) C C---- calculate list of low resolution reflections C keep this list sorted in decreasing quality C DO 90 IH = 0,IHMAX DO 80 IK = -IKMAX,IKMAX C C IF ((IH.NE.0) .OR. (IK.GE.0)) THEN C C DO 70 IL = -ILMAX,ILMAX IF ((IH.NE.0) .OR. (IK.NE.0) .OR. + (IL.GT.0)) THEN C C ************************* CALL GENABS(ABSTYP,IH,IK,IL,I) C ************************* C IF (I.EQ.0) THEN R = 0.0 C C DO 20 I = 1,3 T(I) = A0(I,1)*IH + A0(I,2)*IK + + A0(I,3)*IL R = T(I)**2 + R 20 CONTINUE C C R = SQRT(R) D = ABS(R-R1) IF (D.LE. (0.5*R1)) THEN D = D*D M = N(L) + 1 30 CONTINUE IF (M.NE.1) THEN IF (DELTA(M-1,L).GE.D) THEN M = M - 1 GO TO 30 END IF END IF I = N(L) 40 CONTINUE IF (I.GE.M) THEN IF (I.LT.NMAX) THEN DELTA(I+1,L) = DELTA(I,L) C C DO 50 J = 1,3 V(J,I+1,L) = V(J,I,L) 50 CONTINUE END IF C C I = I - 1 GO TO 40 END IF IF (M.LE.NMAX) THEN DELTA(M,L) = D C C DO 60 J = 1,3 V(J,M,L) = T(J) 60 CONTINUE C C IF (N(L).LT.NMAX) N(L) = N(L) + 1 END IF END IF END IF END IF 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE C C IF (N(L).LT.1) GO TO 170 END IF 100 CONTINUE C C---- find interpretation of basis vectors and save orientation C matrices and qualities of the "mxbest" solutions C D = DELTA(N(1),1) + DELTA(N(2),2) C C DO 160 N1 = 1,N(1) D1 = DELTA(N1,1) C C DO 110 I = 1,3 XCALC(I,1) = V(I,N1,1) 110 CONTINUE C C DO 150 N2 = 1,N(2) D2 = DELTA(N2,2) IF ((D1+D2).GT.D) THEN GO TO 160 ELSE C C DO 120 I = 1,3 XCALC(I,2) = V(I,N2,2) 120 CONTINUE C C DO 140 II = 1,2 C C DO 130 I = 1,3 XCALC(I,2) = -XCALC(I,2) 130 CONTINUE C C---- find orientation matrix u and quality "rms" of superposition C C ********************************** CALL U3BEST(W,XCALC,XOBS,3,1,RMS,U,T,I) C ********************************** C IF (I.LT.0) THEN GO TO 150 ELSE C C---- save assignment if among the best "mxbest" independent choices C C ******************************************* CALL ULIST(NEQL,MLAUE,A0,U,RMS,MXBEST,NBEST,UB, + ESD) C ******************************************* C IF (NBEST.EQ.MXBEST) D = ESD(NBEST) END IF 140 CONTINUE END IF 150 CONTINUE 160 CONTINUE END IF END IF END IF C C---- return C 170 RETURN END C C C SUBROUTINE BMATRX(B,R,C,WAVE) C ============================= C C C C---- Build cell orthogonalization matrix B. R C are the recip cell parameters C and C the real cell parameters. C C Ref: W.R.Busing & H.A.Levy, Acta Cryst. (1967) 22, 457-464 C C B = ( a* b* cos gamma* c* cos beta* ) C ( 0 b* sin gamma* -c* sin beta* cos alpha ) C ( 0 0 lambda / c ) C C C C C .. Scalar Arguments .. REAL WAVE C .. C .. Array Arguments .. REAL B(3,3),C(6),R(6) C .. C .. Intrinsic Functions .. INTRINSIC COSD,SIND C .. SAVE C C B(1,1) = R(1) B(1,2) = COSD(R(6)) * R(2) B(1,3) = COSD(R(5)) * R(3) B(2,1) = 0.000000 B(2,2) = SIND(R(6)) * R(2) B(2,3) = -SIND(R(5)) * R(3) * COSD(C(4)) B(3,1) = 0.00000 B(3,2) = 0.00000 B(3,3) = WAVE/C(3) C END C C C SUBROUTINE CHECKU(U) C ==================== C C C---- Check that U is a pure rotation matrix C C C .. Array Arguments .. REAL U(3,3) C .. INTEGER LP,LINOUT LOGICAL ONLINE C .. Local Scalars .. REAL DET INTEGER I,II,J,JJ C .. C .. Local Arrays .. REAL UWINV(3,3),UWTRANS(3,3) C .. C .. External Subroutines .. EXTERNAL INVERS,TRANSP C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C COMMON /RINOUT/ ONLINE,LP,LINOUT SAVE C C C ******************* CALL INVERS(U,UWINV,DET) C ******************* C C---- Check the determinant C IF (ABS(DET-1.000).GT.0.001) THEN IF (ONLINE) WRITE (LINOUT,FMT=6000) DET IF (ONLINE) WRITE (LINOUT,FMT=6002) ((U(II,JJ),JJ=1,3),II=1,3) WRITE (LP,FMT=6000) DET WRITE (LP,FMT=6002) ((U(II,JJ),JJ=1,3),II=1,3) STOP ELSE C C---- Check that the transpose of U is equal to its inverse C C ***************** CALL TRANSP(UWTRANS,U) C ***************** C DO 20 J = 1,3 DO 10 I = 1,3 IF (ABS(UWTRANS(I,J)-UWINV(I,J)).GT.0.001) GO TO 30 10 CONTINUE 20 CONTINUE C C RETURN 30 IF (ONLINE) WRITE (LINOUT,FMT=6004) IF (ONLINE) WRITE (LINOUT,FMT=6002) ((U(II,JJ),JJ=1,3),II=1,3) WRITE (LP,FMT=6004) WRITE (LP,FMT=6002) ((U(II,JJ),JJ=1,3),II=1,3) STOP END IF C C---- Format statements C 6000 FORMAT (/' ** ERROR ** The determinant of U is ',F9.5,' A value ', + 'of 1.0 would correspond to a simple rotation.') 6002 FORMAT (/' The matrix U defining the standard setting is:', + /3 (20X,3F9.5,/)) 6004 FORMAT (/' ** ERROR ** The matrix Umat is not a simple rotation ', + 'matrix') C C END C C C SUBROUTINE CLCALC(CELL,A) C ========================= C C---- Calculate cell (real or reciprocal) from metric tensor A C Angles in degrees C C C C .. Array Arguments .. REAL A(3,3),CELL(6) C .. C .. Local Scalars .. INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC SQRT,ACOSD C .. C SAVE C C DO 10 I = 1,3 CELL(I) = SQRT(A(I,I)) 10 CONTINUE C C CELL(4) = ACOSD( A(2,3) / ( CELL(2)*CELL(3) ) ) CELL(5) = ACOSD( A(1,3) / ( CELL(1)*CELL(3) ) ) CELL(6) = ACOSD( A(1,2) / ( CELL(1)*CELL(2) ) ) C END C C C SUBROUTINE CROSS(A,ASTT,BSTT,CSTT,VOLUME,WLAMDA) C =============================================== C C C---- To calculate C = A cross B C C C C .. Scalar Arguments .. REAL VOLUME,WLAMDA C .. C .. Array Arguments .. REAL A(3,3),ASTT(3),BSTT(3),CSTT(3) C .. SAVE C C ASTT(1) = ( A(2,2)*A(3,3) - A(2,3)*A(3,2) ) *WLAMDA*VOLUME ASTT(2) = ( A(2,3)*A(3,1) - A(2,1)*A(3,3) ) *WLAMDA*VOLUME ASTT(3) = ( A(2,1)*A(3,2) - A(2,2)*A(3,1) ) *WLAMDA*VOLUME BSTT(1) = ( A(3,2)*A(1,3) - A(3,3)*A(1,2) ) *WLAMDA*VOLUME BSTT(2) = ( A(3,3)*A(1,1) - A(3,1)*A(1,3) ) *WLAMDA*VOLUME BSTT(3) = ( A(3,1)*A(1,2) - A(3,2)*A(1,1) ) *WLAMDA*VOLUME CSTT(1) = ( A(1,2)*A(2,3) - A(1,3)*A(2,2) ) *WLAMDA*VOLUME CSTT(2) = ( A(1,3)*A(2,1) - A(1,1)*A(2,3) ) *WLAMDA*VOLUME CSTT(3) = ( A(1,1)*A(2,2) - A(1,2)*A(2,1) ) *WLAMDA*VOLUME C C END C C C SUBROUTINE CRYSYS(IGROUP,CELL,ICS) C ================================== C C C C C---- Determine crystal system number from spacegroup number C and unit cell parameters. w.kabsch 11-1987 C C C IGROUP - space group number of one of the 65 enantiomorphic (given) C groups as obtained from the international tables. C CELL - unit cell parameters in angstroem and degrees (given) C ICS - number specifying crystal system (result) C 0 : unit cell parameters are inconsistent with C space group C 1 : triclinic C 2 : monoclinic first setting C 3 : monoclinic second setting C 4 : orthorhombic C 5 : tetragonal C 6 : trigonal C 7 : hexagonal C 8 : cubic C C C C .. Scalar Arguments .. INTEGER ICS,IGROUP C .. C .. Array Arguments .. REAL CELL(6) C .. C .. Local Scalars .. INTEGER IG,J C .. SAVE C C ICS = 0 IG = IGROUP C C DO 10 J = 1,6 IF (CELL(J).LE.0.0) RETURN 10 CONTINUE C C---- determine crystal system number from spacegroup number C and unit cell parameters C IF (IG.EQ.1) ICS = 1 C C IF ((IG.GE.3) .AND. (IG.LE.5)) THEN IF ((CELL(4).EQ.90.0) .AND. (CELL(5).EQ.90.0)) ICS = 2 IF ((CELL(4).EQ.90.0) .AND. (CELL(6).EQ.90.0)) ICS = 3 IF (ICS.EQ.3) IG = IG + 512 END IF C C IF ((IG.GE.16) .AND. (IG.LE.24) .AND. (CELL(4).EQ.90.0) .AND. + (CELL(5).EQ.90.0) .AND. (CELL(6).EQ.90.0)) ICS = 4 IF ((IG.GE.75) .AND. (IG.LE.98) .AND. (CELL(1).EQ.CELL(2)) .AND. + (CELL(4).EQ.90.0) .AND. (CELL(5).EQ.90.0) .AND. + (CELL(6).EQ.90.0)) ICS = 5 C C IF ((IG.GE.143) .AND. (IG.LE.155) .AND. (CELL(1).EQ.CELL(2))) THEN IF ((CELL(2).EQ.CELL(3)) .AND. (CELL(4).EQ.CELL(5)) .AND. + (CELL(5).EQ.CELL(6))) ICS = 6 IF ((CELL(4).EQ.90.0) .AND. (CELL(5).EQ.90.0) .AND. + (CELL(6).EQ.120.0)) ICS = 7 IF ((ICS.EQ.6) .AND. ((IG.EQ.146).OR. (IG.EQ.155))) IG = IG + + 512 END IF C C IF ((IG.GE.168) .AND. (IG.LE.182) .AND. (CELL(1).EQ.CELL(2)) .AND. + (CELL(4).EQ.90.0) .AND. (CELL(5).EQ.90.0) .AND. + (CELL(6).EQ.120.0)) ICS = 7 IF ((IG.GE.195) .AND. (IG.LE.214) .AND. (CELL(1).EQ.CELL(2)) .AND. + (CELL(2).EQ.CELL(3)) .AND. (CELL(4).EQ.90.0) .AND. + (CELL(5).EQ.90.0) .AND. (CELL(6).EQ.90.0)) ICS = 8 C C END C C C SUBROUTINE DCOSFD(ROT,DIRCOS,KAPPA) C =================================== C C C Given a rotation matrix ROTN, expressing a rotation C through an angle KAPPA right-handedly about an axis C with direction cosines DIRCOS(I), C DCOSFD determines DIRCOS() and KAPPA C KAPPA is returned in degrees in the range 0 to 180 C C Expression for rotation matrix is (see J & J, "Mathl. C Phys.", p. 122):- C C cw + n(1)n(1)(1-cw) n(1)n(2)(1-cw)-n(3)sw n(1)n(3)(1-cw)+n(3)sw C n(1)n(2)(1-cw)+n(3)sw cw + n(2)n(2)(1-cw) n(2)n(3)(1-cw)-n(1)sw C n(3)n(1)(1-cw)-n(2)sw n(3)n(2)(1-cw)+n(1)sw cw + n(3)n(3)(1-cw) C C where cw = cos(KAPPA), sw = sin(KAPPA), C & n(1), n(2), n(3) are the direction cosines. C C C C C C C C .. Scalar Arguments .. REAL KAPPA C .. C .. Array Arguments .. REAL DIRCOS(3),ROT(3,3) C .. INTEGER LP,LINOUT LOGICAL ONLINE C .. Local Scalars .. REAL COSKAP,DIFF,PI,R2,R3,SINKAP,TRACE INTEGER J,JMX,KNTRL C .. C .. Local Arrays .. REAL D(3),P(3),S(3) INTEGER NCS(3),ND(3),NP(3),NS(3) C .. C .. External Subroutines .. EXTERNAL ORDR3 C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,ATAN2,SQRT C .. COMMON /RINOUT/ ONLINE,LP,LINOUT C C SAVE C C DIFF = 1.0E-04 PI = ATAN(1.0)*4.0 KNTRL = 0 C C---- Trace = 1 + 2.cos w C TRACE = ROT(1,1) + ROT(2,2) + ROT(3,3) C C---- S(1)=2.n(1).sin w C S(2)=2.n(2).sin w C S(3)=2.n(3).sin w C S(1) = ROT(3,2) - ROT(2,3) S(2) = ROT(1,3) - ROT(3,1) S(3) = ROT(2,1) - ROT(1,2) C C ***************** CALL ORDR3(KNTRL,S,NS) C ***************** C C---- Is biggest S() zero (i.e. is sin w = 0)? C IF (ABS(S(NS(1))).GT.DIFF) THEN R2 = S(NS(2))/S(NS(1)) R3 = S(NS(3))/S(NS(1)) DIRCOS(NS(1)) = 1.0/SQRT(R2**2+1+R3**2) DIRCOS(NS(2)) = DIRCOS(NS(1))*R2 DIRCOS(NS(3)) = DIRCOS(NS(1))*R3 ELSE C C---- Calculation when sin w = 0 (esp. when w = 180'). C P(1)=n(2).n(3).{0(W=0)or2(w=180)} C P(2)=n(3).n(1).{0(W=0)or2(w=180)} C P(3)=n(1).n(2).{0(W=0)or2(w=180)} C P(1) = ROT(3,2) P(2) = ROT(1,3) P(3) = ROT(2,1) C C ***************** CALL ORDR3(KNTRL,P,NP) C ***************** C C---- Is biggest P() zero (i.e. are all off-diag. terms zero)? C IF (ABS(P(NP(1))).LT.DIFF) THEN IF (TRACE.GT.0.0) THEN C C---- Matrix is unit matrix. C KAPPA = 0.0 DIRCOS(1) = 0.0 DIRCOS(2) = 0.0 DIRCOS(3) = 1.0 ELSE C C---- Trace -ve, so dyad about x,y or z. C KAPPA = PI C C DO 10 J = 1,3 D(J) = ROT(J,J) + 1.0 10 CONTINUE C C ***************** CALL ORDR3(KNTRL,D,ND) C ***************** C DIRCOS(ND(1)) = 1.0 DIRCOS(ND(2)) = 0.0 DIRCOS(ND(3)) = 0.0 END IF RETURN C C ELSE IF (ABS(P(NP(2))).LT.DIFF) THEN C C---- One d.c. is zero. C DIRCOS(NP(1)) = 0.0 DIRCOS(NP(2)) = SQRT((ROT(NP(2),NP(2))+1.0)*0.5) DIRCOS(NP(3)) = SQRT((ROT(NP(3),NP(3))+1.0)*0.5) IF (P(NP(1)).LT.0.0) DIRCOS(NP(2)) = -DIRCOS(NP(2)) ELSE C C---- All d.c's are nonzero. C R2=n(NP(1))/n(NP(2)) C R3=n(NP(1))/n(NP(3)) C R2 = P(NP(2))/P(NP(1)) R3 = P(NP(3))/P(NP(1)) DIRCOS(NP(2)) = 1.0/SQRT(R2**2+1+ (R2/R3)**2) DIRCOS(NP(3)) = 1.0/SQRT(R3**2+1+ (R3/R2)**2) DIRCOS(NP(1)) = SQRT(1.0-DIRCOS(NP(2))**2-DIRCOS(NP(3))**2) C C---- Adjust signs of DIRCOS(). C JMX = 0 C C DO 20 J = 1,3 IF (P(J).GT.0.0) THEN IF (JMX.GT.0) JMX = -1 IF (JMX.EQ.0) JMX = J END IF 20 CONTINUE C C IF (JMX.GT.0) DIRCOS(JMX) = -DIRCOS(JMX) IF (JMX.EQ.0 .AND. ONLINE) WRITE (LINOUT,FMT=6000) IF (JMX.EQ.0) WRITE (LP,FMT=6000) END IF END IF C C---- Given d.c's, calculate KAPPA. C KNTRL = 1 C C *********************** CALL ORDR3(KNTRL,DIRCOS,NCS) C *********************** C C---- Find KAPPA. C SINKAP = ROT(NCS(3),NCS(2)) - ROT(NCS(2),NCS(3)) SINKAP = 0.5*SINKAP/DIRCOS(NCS(1)) COSKAP = (TRACE-1.0)*0.5 KAPPA = ATAN2(SINKAP,COSKAP) KAPPA = KAPPA*180.0/PI C C---- If kappa negative, negate kappa and invert DIRCOS C IF (KAPPA.LT.0.0) THEN KAPPA = -KAPPA DIRCOS(1) = -DIRCOS(1) DIRCOS(2) = -DIRCOS(2) DIRCOS(3) = -DIRCOS(3) END IF C C---- Format statements C 6000 FORMAT (' All P()s. are negative or zero') C C END C C C SUBROUTINE DECOMP(A,U,B) C ======================== C C C C---- Calculate U and B given A. C C C C C---- Orientation matrix given explicitly C C Set UMAT to (UB)transpose C C .. Array Arguments .. REAL A(3,3),B(3,3),U(3,3) C .. C .. Scalars in Common .. REAL CCOM,CCX,CCY,DPHI,F,FIXF,FLAMDA,ORGX,ORGY,Q,THRESH INTEGER IGROUP LOGICAL DCOMP,FILM,LCAMC,TARGET INTEGER LP,LINOUT LOGICAL ONLINE C .. C .. Arrays in Common .. REAL ACHSE,CELL,ED,S0,TARMAT C .. C .. Local Scalars .. REAL D INTEGER I C .. C .. Local Arrays .. REAL CELLC(6),GMAT(3,3),RCELL(6),XMAT(3,3) C .. C .. External Subroutines .. EXTERNAL BMATRX,CHECKU,CLCALC,INVERS,MATMUL,TRANSP C .. C .. Common blocks .. COMMON /REFCOM/IGROUP,FLAMDA,CELL(6),F,THRESH,Q,DPHI,ED(3,3), + ACHSE(3),S0(3),ORGX,ORGY,FIXF,TARMAT(3,3),TARGET,DCOMP, + LCAMC,CCX,CCY,CCOM,FILM COMMON /RINOUT/ ONLINE,LP,LINOUT C .. SAVE C C C *********** CALL TRANSP(U,A) C *********** C C---- Get reciprocal metric tensor G**-1 in GMAT, = (UB)T.(UB) C C **************** CALL MATMUL(U,A,GMAT) C **************** C C---- Get reciprocal cell dimensions RCELL C C ****************** CALL CLCALC(RCELL,GMAT) C ****************** C C---- Invert to get real metric tensor G in UMAT C C **************** CALL INVERS(GMAT,U,D) C **************** C IF (D.EQ.0.0) THEN IF (ONLINE) WRITE (LINOUT,FMT=6000) WRITE (LP,FMT=6000) STOP ELSE C C---- Get real cell dimensions from metric tensor C C *************** CALL CLCALC(CELLC,U) C *************** C DO 10 I = 1,3 CELLC(I) = CELLC(I)*FLAMDA 10 CONTINUE C C---- Rebuild cell orthogonalization matrix B in BMAT C C **************************** CALL BMATRX(B,RCELL,CELLC,FLAMDA) C **************************** C C---- Get matrix U = (UB).(B**-1) to UMAT C C **************** CALL INVERS(B,XMAT,D) CALL MATMUL(A,XMAT,U) CALL CHECKU(U) C ***************** C END IF C C---- Format statements C 6000 FORMAT (//' !!!! Zero reciprocal cell volume from orientation ma', + 'trix !!!!',//) C C END C C C SUBROUTINE DIFVEC(NRFL,V,DELTA,SNR,GRID,NBXY,NBZ,MXY,MZ,MZSTEP, + MXSPOT,TABLE,I2X,I2Y,I2Z,INTY,VC,NSPOT,NVEC) C ================================================================ C C C C C---- Obtain list of difference vector clusters C C W.KABSCH VERSION 3-1987 C C C---- Description of calling parameters C C NRFL - number of reflections in array v (given) C V - real array(3,*) of reflection coordinates (given) C DELTA - estimated positional error (number of grid points) (given) C of difference vector cluster. C SNR - signal to noise ratio passed to subroutine spots (given) C for searching difference vector clusters C GRID - grid spacing such that v(i,j)/grid are grid (given) C coordinates of vector #j C NBXY - defines size of search box when looking for (given) C difference vector clusters. passed to subroutine C "spots". C NBZ - defines width (grid units) along z' as 2*nbz+1 (given) C when searching for difference vector clusters C MXY - specifies dimensions in array table(mxy,mxy) (given) C MZ - defines number of grid steps along direction z' (given) C 0,...,mz when searching for difference vector C clusters C MZSTEP - increment along direction z' (grid units) (given) C MXSPOT - maximum number of difference vector clusters. (given) C it defines dimensions of arrays i2x,i2y,i2z,inty,vc C TABLE - used to accumulate short vector differences (changed) C between reflection coordinates C I2X - grid coordinates (grid units * 10 ) along x', (result) C I2Y - grid coordinates (grid units * 10 ) along y', (result) C I2Z - grid coordinates (grid units * 10 ) along z' and (result) C INTY - population of difference vector clusters. (result) C VC - reciprocal space coordinates of difference vector (result) C clusters. C NSPOT - number of difference vector clusters found (result) C NVEC - number of difference vectors between reflections (result) C lying in clusters C C---- The coordinate system {x',y',z'} is determined by the C principal axes of an ellipsoid fitted to the list of C short difference vectors. the ellipsoid is defined with C decreasing length of the principal axes along x',y',z'. C C C C---- external subroutines required: SPOTS,EIGEN C C C C C C C .. Scalar Arguments .. REAL DELTA,GRID,SNR INTEGER MXSPOT,MXY,MZ,MZSTEP,NBXY,NBZ,NRFL,NSPOT INTEGER NVEC C .. C .. Array Arguments .. REAL V(3,*),VC(3,MXSPOT) INTEGER*2 I2X(MXSPOT),I2Y(MXSPOT),I2Z(MXSPOT),INTY(MXSPOT), + TABLE(MXY,MXY) C .. C .. Local Scalars .. REAL FAC,P,R,RDELTA,SIGMA,X,Y,Z INTEGER I,ISPOT,IZ,J,JRFL,JSPOT,JX,JY,JZ,K,KRFL,N,N1,N1TEN,NSPOT1 INTEGER NSPOT0,NV C .. C .. Local Arrays .. REAL A(6),D(3),E(3,3) C .. C .. External Subroutines .. EXTERNAL EIGEN,SPOTS C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD,NINT C .. SAVE C C NSPOT0 = 0 NSPOT = 0 NVEC = 0 N = (MXY-1)/2 N1 = N + 1 N1TEN = 10*N1 C C IF ((GRID.GT.0.0) .AND. (NRFL.GE.2) .AND. (MXY.GE.3)) THEN IF ((MXSPOT.GE.1) .AND. (MZ.GE.0) .AND. (MZSTEP.GE.1)) THEN IF ((NBXY.GE.1) .AND. (N.GE.NBXY) .AND. (NBZ.GE.0)) THEN RDELTA = N*GRID SIGMA = 10.0*DELTA C C---- fit a flat disc to short difference vectors C P = RDELTA**2 C C DO 10 K = 1,6 A(K) = 0.0 10 CONTINUE C C DO 60 JRFL = 1,NRFL - 1 DO 50 KRFL = JRFL + 1,NRFL R = 0.0 C C DO 20 I = 1,3 D(I) = V(I,JRFL) - V(I,KRFL) R = D(I)**2 + R 20 CONTINUE C C IF (R.LE.P) THEN K = 0 C C DO 40 I = 1,3 DO 30 J = 1,I K = K + 1 A(K) = D(I)*D(J) + A(K) 30 CONTINUE 40 CONTINUE END IF 50 CONTINUE C C 60 CONTINUE C C ************ CALL EIGEN(3,A,E) C ************ C C---- accumulate difference vectors to locate clusters C DO 120 IZ = 0,MZ,MZSTEP C C DO 80 J = 1,MXY DO 70 K = 1,MXY TABLE(J,K) = 0 70 CONTINUE 80 CONTINUE C C NV = 0 C C DO 100 JRFL = 1,NRFL - 1 DO 90 KRFL = JRFL + 1,NRFL X = (V(1,JRFL)-V(1,KRFL))/GRID Y = (V(2,JRFL)-V(2,KRFL))/GRID Z = (V(3,JRFL)-V(3,KRFL))/GRID JX = NINT(E(1,1)*X+E(2,1)*Y+E(3,1)*Z) IF (ABS(JX).LE.N) THEN JY = NINT(E(1,2)*X+E(2,2)*Y+E(3,2)*Z) IF (ABS(JY).LE.N) THEN IF ((JX.NE.0) .OR. (JY.NE.0)) THEN JZ = NINT(E(1,3)*X+E(2,3)*Y+E(3,3)*Z) IF (ABS(IZ-JZ).LE.NBZ) THEN I = TABLE(N1+JX,N1+JY) IF (I.LT.32767) THEN NV = NV + 1 TABLE(N1+JX,N1+JY) = I + 1 END IF END IF IF (ABS(IZ+JZ).LE.NBZ) THEN I = TABLE(N1-JX,N1-JY) IF (I.LT.32767) THEN NV = NV + 1 TABLE(N1-JX,N1-JY) = I + 1 END IF END IF END IF END IF END IF 90 CONTINUE 100 CONTINUE C C IF (NV.GE.2) THEN C C *********************************************** CALL SPOTS(MXY,MXY,NBXY,NBXY,SNR,0.0,1.0,TABLE,TABLE, + MXSPOT,NSPOT0,NSPOT,I2X,I2Y,INTY) C *********************************************** C IF (NSPOT.GE.1) THEN NSPOT1 = NSPOT0 + NSPOT C C DO 110 ISPOT = NSPOT0 + 1,NSPOT1 J = MOD(ISPOT-1,MXSPOT) + 1 I2Z(J) = IZ 110 CONTINUE C C NSPOT0 = NSPOT1 END IF END IF 120 CONTINUE C C---- clean-up list of difference vector clusters C IF (NSPOT0.GT.MXSPOT) NSPOT0 = MXSPOT IF (NSPOT0.GE.2) THEN C C DO 140 ISPOT = 1,NSPOT0 - 1 DO 130 JSPOT = ISPOT + 1,NSPOT0 IF (ABS(I2Y(JSPOT)-I2Y(ISPOT)).LE.SIGMA) THEN IF (ABS(I2X(JSPOT)-I2X(ISPOT)).LE.SIGMA) THEN C C IF (ABS(INTY(ISPOT)).LT.ABS(INTY(JSPOT))) THEN INTY(ISPOT) = -ABS(INTY(ISPOT)) ELSE INTY(JSPOT) = -ABS(INTY(JSPOT)) END IF END IF END IF 130 CONTINUE 140 CONTINUE C C END IF NSPOT = 0 C C DO 150 ISPOT = 1,NSPOT0 IF (INTY(ISPOT).GT.0) THEN NSPOT = NSPOT + 1 I2X(NSPOT) = I2X(ISPOT) - N1TEN I2Y(NSPOT) = I2Y(ISPOT) - N1TEN I2Z(NSPOT) = I2Z(ISPOT)*10 INTY(NSPOT) = INTY(ISPOT) END IF 150 CONTINUE C C IF (NSPOT.GE.2) THEN C C DO 170 ISPOT = 1,NSPOT - 1 DO 160 JSPOT = ISPOT + 1,NSPOT IF (ABS(I2X(JSPOT)+I2X(ISPOT)).LE.SIGMA) THEN IF (ABS(I2Y(JSPOT)+I2Y(ISPOT)).LE.SIGMA) THEN IF (ABS(I2Z(JSPOT)+I2Z(ISPOT)).LE.SIGMA) THEN IF (INTY(ISPOT).GT.0) THEN IF (INTY(JSPOT).GT.0) THEN I2X(ISPOT) = NINT((I2X(ISPOT)-I2X(JSPOT))* + 0.5) I2Y(ISPOT) = NINT((I2Y(ISPOT)-I2Y(JSPOT))* + 0.5) I2Z(ISPOT) = NINT((I2Z(ISPOT)-I2Z(JSPOT))* + 0.5) INTY(JSPOT) = -ABS(INTY(JSPOT)) END IF END IF END IF END IF END IF 160 CONTINUE 170 CONTINUE END IF C C NSPOT1 = NSPOT NSPOT = 0 FAC = GRID/10.0 C C DO 190 ISPOT = 1,NSPOT1 IF (INTY(ISPOT).GT.0) THEN NSPOT = NSPOT + 1 X = I2X(ISPOT)*FAC Y = I2Y(ISPOT)*FAC Z = I2Z(ISPOT)*FAC I2X(NSPOT) = I2X(ISPOT) I2Y(NSPOT) = I2Y(ISPOT) I2Z(NSPOT) = I2Z(ISPOT) INTY(NSPOT) = INTY(ISPOT) NVEC = INTY(NSPOT) + NVEC C C DO 180 J = 1,3 VC(J,NSPOT) = E(J,1)*X + E(J,2)*Y + E(J,3)*Z 180 CONTINUE END IF 190 CONTINUE C C END IF END IF END IF END C C C REAL FUNCTION DOT(A,B) C ====================== C C C---- Dot product of two vectors C C C C .. Array Arguments .. REAL A(3),B(3) C .. SAVE C C DOT = A(1)*B(1) + A(2)*B(2) + A(3)*B(3) C C END C C C SUBROUTINE EIGEN(N,A,E) C ======================= C C C---- Calculates eigenvalues and eigenvectors of a real symmetric C matrix "a" by the jacobi method. C C C N - order of matrices "a" and "e" (n>1 required) (given) C A - upper triangular part of a real symmetric matrix. (modified) C on return the diagonal of "a" contains the C eigenvalues in decreasing order. C dimension of "a" must be at least n*(n+1)/2. C element i,j must have been stored at address C I+(J*J-J)/2 IF I0 : no error. C 0 : insufficient number of accepted reflections C -1 : q is not positive C -2 : f is zero C -3 : illegal cell parameters C -4 : illegal rotation axis or direct beam wave vector C NOVLAP - number of spots which could not be indexed because (result) C more than one predicted position was found to be C within 3*sdxy and 3*sdphi. to indicate this case C reflection indices 0,0,0 are assigned. C NFAR - number of spots which could not be indexed because (result) C of a distance greater than 3*sdxy or 3*sdphi to any C possible predicted position. to indicate this case C reflection indices 0,0,0 are assigned. C SDPHI - standard deviation of angular position of spindle (modified) C ( degrees ). a positive input value is used to C compare differences between observed and calculated C spindle positions of the reflections. C a non-positive value indicates that a default C value of 0.1 degrees should be used. C on return, sdphi contains the estimated standard C deviation of the spindle angle at which the C reflections were diffracting. C SDXY - standard deviation of spot position on detector (modified) C (pixels). a positive input value is used to C compare differences between observed and calculated C spot positions of the reflections. C a non-positive value indicates that a default C value of 3 pixels should be used. C on return, sdxy contains the estimated standard C deviation of the spot positions on the detector. C C C C---- check calling parameters C C .. Parameters .. REAL RAD,RAD100,CUT PARAMETER (RAD=57.29578,RAD100=100.0*RAD,CUT=9.0) C .. C .. Scalar Arguments .. REAL F,ORGX,ORGY,Q,SDPHI,SDXY INTEGER ABSTYP,MXRFL,NACC,NFAR,NOVLAP,NRFL INTEGER NRFL0 C .. C .. Array Arguments .. REAL ACHSE(3),ED(3,3),RCELL(6),S0L(3),U(3,3) INTEGER*2 IH(*),IK(*),IL(*),IPHI(*),IX(*),IY(*) C .. C .. Local Scalars .. REAL CPHI,CPHIC,DELPHI,DELTA,DELTA0,DELXY,ERRPHI,ERRXY,FH,FK,FL,R, + RHOQ,RR,SPHI,SPHIC,SS,WPHI,WXY,X,XOBS,Y,YOBS,Z INTEGER I1,I2,I3,IDX,J,J1,J2,J3,JRFL,K,NIX,NY INTEGER*2 I2PHI C .. C .. Local Arrays .. REAL A(3,3),A0(3,3),B(3,3),EG(3,3),PD(3),S0G(3),SD(3),SL(3), + X0G(3),X0L(3),XG(3),XL(3) INTEGER IHKL(6) C .. C .. External Subroutines .. EXTERNAL GENABS,GONSYS,INVERS,MATMUL,REFLEX,SETMAT C .. C .. Intrinsic Functions .. INTRINSIC COS,MOD,SIGN,SIN,SQRT C .. SAVE C C IF ((MXRFL.LT.1) .OR. (NRFL0.LT.0) .OR. (NRFL.LT.1)) THEN NACC = 0 C C ELSE IF (Q.LE.0.0) THEN NACC = -1 C C ELSE IF (F.EQ.0.0) THEN NACC = -2 ELSE C C---- calculate setting matrix and check for positive volume C C **************** CALL SETMAT(RCELL,A0) CALL MATMUL(U,A0,A) CALL INVERS(A,B,R) C **************** C IF (R.LE.0.0) THEN NACC = -3 ELSE C C---- define goniostat system C C ************************** CALL GONSYS(ACHSE,S0L,S0G,EG,J) C ************************** C SS = S0L(1)**2 + S0L(2)**2 + S0L(3)**2 C C IF ((J.NE.0) .OR. (SS.LE.0.0)) THEN NACC = -4 ELSE C C---- get test values to compare differences between observed C and predicted quantities. C R = 0.1 IF (SDPHI.GT.0.0) R = SDPHI WPHI = (RAD/R)**2 WXY = 1.0/ (3.0*Q)**2 IF (SDXY.GT.0.0) WXY = 1.0/ (SDXY*Q)**2 C C---- indexing loop C NACC = 0 NFAR = 0 NOVLAP = 0 SDXY = 0.0 SDPHI = 0.0 I2PHI = 0 CPHI = 1.0 SPHI = 0.0 C C DO 80 JRFL = 1,NRFL NY = MOD(JRFL-1+NRFL0,MXRFL) + 1 XOBS = (IX(NY)/10.0-ORGX)*Q YOBS = (IY(NY)/10.0-ORGY)*Q C C IF (I2PHI.NE.IPHI(NY)) THEN I2PHI = IPHI(NY) R = I2PHI/RAD100 CPHI = COS(R) SPHI = SIN(R) END IF C C X = ED(1,1)*XOBS + ED(1,2)*YOBS + ED(1,3)*F Y = ED(2,1)*XOBS + ED(2,2)*YOBS + ED(2,3)*F Z = ED(3,1)*XOBS + ED(3,2)*YOBS + ED(3,3)*F R = SQRT((X*X+Y*Y+Z*Z)/SS) X = X/R - S0L(1) Y = Y/R - S0L(2) Z = Z/R - S0L(3) R = (ACHSE(1)*X+ACHSE(2)*Y+ACHSE(3)*Z)* (1.0-CPHI) X0L(1) = ACHSE(1)*R + X*CPHI - + (ACHSE(2)*Z-ACHSE(3)*Y)*SPHI X0L(2) = ACHSE(2)*R + Y*CPHI - + (ACHSE(3)*X-ACHSE(1)*Z)*SPHI X0L(3) = ACHSE(3)*R + Z*CPHI - + (ACHSE(1)*Y-ACHSE(2)*X)*SPHI C C DO 10 J = 1,3 R = B(J,1)*X0L(1) + B(J,2)*X0L(2) + B(J,3)*X0L(3) K = SIGN(0.5,R) + R IHKL(J) = K IF (R.LT.K) K = K - 2 IHKL(J+3) = K + 1 10 CONTINUE DELTA0 = 1.0E06 NIX = 0 C C DO 70 I1 = 1,6,3 J1 = IHKL(I1) FH = J1 C C DO 60 I2 = 2,6,3 J2 = IHKL(I2) FK = J2 C C DO 50 I3 = 3,6,3 J3 = IHKL(I3) FL = J3 C C *************************** CALL GENABS(ABSTYP,J1,J2,J3,IDX) C *************************** C IF (IDX.EQ.0) THEN C C---- get laboratory coordinates of reflection at phi=0 C DO 20 J = 1,3 X0L(J) = A(J,1)*FH + A(J,2)*FK + A(J,3)*FL 20 CONTINUE C C---- calculate diffraction geometry for this reflection C C ******************************************** CALL REFLEX(EG,CPHI,SPHI,S0L,S0G,X0L,X0G,XL,XG,SL, + CPHIC,SPHIC,RR,RHOQ,IDX) C ******************************************** C IF (IDX.EQ.0) THEN C C---- get detector coordinates of diffracted beam wavevector C DO 30 J = 1,3 SD(J) = ED(1,J)*SL(1) + ED(2,J)*SL(2) + + ED(3,J)*SL(3) 30 CONTINUE C C IF (SD(3).NE.0.0) THEN C C---- get detector coordinates of intersection point of diffracted C beam with detector surface C DO 40 J = 1,3 PD(J) = SD(J)*F/SD(3) 40 CONTINUE C C DELXY = ((PD(1)-XOBS)**2+ (PD(2)-YOBS)**2)*WXY IF (DELXY.LE.CUT) THEN DELPHI = ((CPHI-CPHIC)**2+ (SPHI-SPHIC)**2)* + WPHI*RHOQ/RR IF (DELPHI.LE.CUT) THEN DELTA = DELPHI + DELXY C C IF (DELTA.LT.DELTA0) THEN DELTA0 = DELTA ERRXY = DELXY ERRPHI = DELPHI NIX = NIX + 1 IH(NY) = J1 IK(NY) = J2 IL(NY) = J3 END IF END IF END IF END IF END IF END IF 50 CONTINUE C C 60 CONTINUE 70 CONTINUE C C IF ((NIX-1).EQ.0) THEN NACC = NACC + 1 SDXY = SDXY + ERRXY SDPHI = SDPHI + ERRPHI ELSE IF ((NIX-1).GT.0) THEN NOVLAP = NOVLAP + 1 ELSE C C NFAR = NFAR + 1 END IF IH(NY) = 0 IK(NY) = 0 IL(NY) = 0 END IF 80 CONTINUE IF (NACC.GE.1) THEN SDXY = SQRT(SDXY/ (WXY*NACC))/Q SDPHI = SQRT(SDPHI/ (WPHI*NACC))*RAD END IF RETURN END IF END IF END IF END C C C SUBROUTINE INVCEL(CELL,RCELL,V) C =============================== C C C---- Calculate reciprocal cell parameters from unit cell constants. C this routine may also be used to calculate the unit cell C parameters from the reciprocal cell constants. C C CELL - unit cell constants : a,b,c,alfa,beta,gamma (given) C in angstroem and degrees C RCELL - calculated reciprocal cell parameters (result) C in reciprocal angstroem and degrees C V - calculated unit cell volume in angstroem**3 . (result) C v<=0.0 indicates an error situation. C C C C C .. Scalar Arguments .. REAL V C .. C .. Array Arguments .. REAL CELL(6),RCELL(6) C .. C .. Local Scalars .. REAL ARG INTEGER I,J,K C .. C .. Local Arrays .. REAL CA(3),SA(3) C .. C .. Intrinsic Functions .. INTRINSIC ACOS,COS,MOD,SIN,SQRT C .. SAVE C C DO 10 I = 1,3 ARG = CELL(I+3)/57.29578 CA(I) = COS(ARG) SA(I) = SIN(ARG) 10 CONTINUE C C V = CELL(1)*CELL(2)*CELL(3)*SQRT(CA(1)*2.0*CA(2)*CA(3)+1.0- + CA(1)**2-CA(2)**2-CA(3)**2) IF (V.GT.0.0) THEN C C DO 20 I = 1,3 J = MOD(I,3) + 1 K = MOD(I+1,3) + 1 RCELL(I) = CELL(J)*CELL(K)*SA(I)/V ARG = (CA(J)*CA(K)-CA(I))/ (SA(J)*SA(K)) RCELL(I+3) = ACOS(ARG)*57.29578 20 CONTINUE END IF RETURN C C END C C C SUBROUTINE INVERS(A,B,DET) C ========================== C C C A - 3 x 3 matrix (given) C B - 3 x 3 matrix is calculated inverse matrix of a (result) C DET - calculated determinant of given matrix a (result) C C C C C .. Scalar Arguments .. REAL DET C .. C .. Array Arguments .. REAL A(3,3),B(3,3) C .. C .. Local Scalars .. INTEGER I,J,K C .. SAVE C C DO 10 I = 1,3 J = I + 1 IF (J.GT.3) J = J - 3 K = I + 2 IF (K.GT.3) K = K - 3 B(I,1) = A(2,J)*A(3,K) - A(2,K)*A(3,J) B(I,2) = A(3,J)*A(1,K) - A(3,K)*A(1,J) B(I,3) = A(1,J)*A(2,K) - A(1,K)*A(2,J) 10 CONTINUE C C DET = A(1,1)*B(1,1) + A(1,2)*B(2,1) + A(1,3)*B(3,1) IF (DET.NE.0.0) THEN C C DO 30 I = 1,3 DO 20 J = 1,3 B(I,J) = B(I,J)/DET 20 CONTINUE 30 CONTINUE END IF RETURN C C END C C C SUBROUTINE LAUESY(CELL,IGROUP,NEQ,ML) C ===================================== C C C C C---- Obtain equivalent positions of laue-group from space-group number. C The space-group number for your crystal can be looked up in the C international tables for x-ray crystallography vol.i C C W.KABSCH 7-1986 C C C CELL - unit cell parameters in angstroem and degrees (given) C IGROUP - space group number as obtained from the international (given) C tables for x-ray crystallography vol.i C NEQ - number of equivalent positions in laue group (result) C derived from space-group number. a value of zero C indicates an illegal space-group number or illegal C unit cell parameters. otherwise neq assumes values C between 1 and 24. C ML - integer array of length at least neq*9 specified (result) C in the main program, to describe symmetry. all C friedel related positions are omitted. the largest C number of neq is 24, as mentioned. therefore, the C calling program should declare C integer ml(216) to handle all possible cases!!! C the operators returned by this subroutine are C proper rotations. to obtain the complete laue-group C you only have to add the 'neq' friedel mates to C the matrices given by this routine. usually this is C not wanted. C each operator is represented by 9 consecutive C numbers as follows: C if x',y',z' is an equivalent position of x,y,z then C X'=X*ML(1)+Y*ML(2)+Z*ML(3) C Y'=X*ML(4)+Y*ML(5)+Z*ML(6) C Z'=X*ML(7)+Y*ML(8)+Z*ML(9) C The next operator occupies positions ml(10)...ml(18), C and so on... C C .. Parameters .. REAL EPS PARAMETER (EPS=0.001) C .. C .. Scalar Arguments .. INTEGER IGROUP,NEQ C .. C .. Array Arguments .. REAL CELL(6) INTEGER ML(*) C .. C .. Local Scalars .. INTEGER I,IG,J,K,L,M C .. C .. Local Arrays .. INTEGER*2 MAT(288),MATN(56),NA(14),NM(14) C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. SAVE C .. Data statements .. DATA MAT/1,3*0,1,3*0,1,-1,3*0,-1,3*0,1,1,3*0,-1,3*0,-1,-1,3*0,1, + 3*0,-1,0,1,0,-1,4*0,1,0,-1,0,1,4*0,1,0,-1,0,-1,4*0,-1,0,1,0, + 1,4*0,-1,0,0,1,1,3*0,1,0,0,0,1,-1,3*0,-1,0,0,1,3*0,1,1,0,0,0, + 1,3*0,-1,-1,0,0,0,0,-1,1,3*0,-1,0,0,0,-1,-1,3*0,1,0,0,-1,3*0, + 1,-1,0,0,0,-1,3*0,-1,1,0,0,-1,4*0,-1,0,-1,0,0,0,-1,0,-1,0,-1, + 0,0,-1,4*0,1,0,1,0,0,0,-1,0,1,0,1,0,0,1,4*0,-1,0,1,0,0,0,1,0, + -1,0,1,0,0,1,4*0,1,0,-1,0,0,0,1,0,1,0,-1,0,0,0,-1,0,1,-1,3*0, + 1,-1,1,0,-1,4*0,1,-1,0,0,-1,1,3*0,-1,1,-1,0,0,-1,3*0,-1,0,1, + 0,-1,1,3*0,1,1,-1,0,1,4*0,1,1,0,0,1,-1,3*0,-1,-1,1,0,0,1,3*0, + -1/ DATA MATN/1,4,3,2,9,10,11,12,13,14,15,16,5,6,7,8,17,18,19,20,21, + 22,23,24,1,2,5,6,3,4,7,8,1,25,26,8,27,28,1,9,11,7,17,18,1,25, + 26,2,29,30,7,31,32,8,27,28/ DATA NM/1,2,4,12,24,2,4,8,3,6,3,6,6,12/ DATA NA/0,0,0,0,0,24,24,24,32,32,38,38,44,44/ C .. C C---- check unit cell parameters C IG = 0 NEQ = 0 C C DO 10 J = 1,6 IF (CELL(J).LE.0.0) GO TO 40 10 CONTINUE C C---- find laue group from space-group number and unit cell parameters C ..... TRICLINIC AXES C IF ((IGROUP.GT.0) .AND. (IGROUP.LT.3)) IG = 1 C C ..... Monoclinic axes C IF ((IGROUP.GT.2) .AND. (IGROUP.LT.16) .AND. + (ABS(CELL(4)-90.0).LT.EPS)) THEN IF (ABS(CELL(5)-90.0).LT.EPS) IG = 6 IF (ABS(CELL(6)-90.0).LT.EPS) IG = 2 END IF C C ..... Orthorhombic axes C IF ((ABS(CELL(4)-90.0).LT.EPS) .AND. + (ABS(CELL(5)-90.0).LT.EPS) .AND. + (ABS(CELL(6)-90.0).LT.EPS)) THEN IF ((IGROUP.GT.15) .AND. (IGROUP.LT.75)) IG = 3 C C ..... Tetragonal axes C IF (ABS(CELL(1)-CELL(2)).LT.EPS) THEN IF ((IGROUP.GT.74) .AND. (IGROUP.LT.89)) IG = 7 IF ((IGROUP.GT.88) .AND. (IGROUP.LT.143)) IG = 8 C C ..... Cubic axes C IF (ABS(CELL(1)-CELL(3)).LT.EPS) THEN IF ((IGROUP.GT.194) .AND. (IGROUP.LT.207)) IG = 4 IF ((IGROUP.GT.206) .AND. (IGROUP.LT.231)) IG = 5 END IF END IF END IF C C ..... hexagonal axes C IF ((ABS(CELL(1)-CELL(2)).LT.EPS) .AND. + (ABS(CELL(4)-90.0).LT.EPS) .AND. + (ABS(CELL(5)-90.0).LT.EPS) .AND. + (ABS(CELL(6)-120.0).LT.EPS)) THEN IF ((IGROUP.GT.142) .AND. (IGROUP.LT.149)) IG = 9 IF ((IGROUP.GT.148) .AND. (IGROUP.LT.168)) IG = 10 IF ((IGROUP.GT.167) .AND. (IGROUP.LT.177)) IG = 13 IF ((IGROUP.GT.176) .AND. (IGROUP.LT.195)) IG = 14 END IF C C ..... rhombohedral axes C IF ((ABS(CELL(1)-CELL(2)).LT.EPS) .AND. + (ABS(CELL(1)-CELL(3)).LT.EPS) .AND. + (ABS(CELL(4)-CELL(5)).LT.EPS) .AND. + (ABS(CELL(4)-CELL(6)).LT.EPS)) THEN IF ((IGROUP.EQ.146) .OR. (IGROUP.EQ.148)) IG = 11 IF ((IGROUP.EQ.155) .OR. (IGROUP.EQ.160) .OR. + (IGROUP.EQ.161) .OR. (IGROUP.EQ.166) .OR. + (IGROUP.EQ.167)) IG = 12 END IF C C---- get matrices of equivalent positions for laue group C IF (IG.GE.1) THEN NEQ = NM(IG) K = NA(IG) I = 0 C C DO 30 J = 1,NEQ L = (MATN(J+K)-1)*9 DO 20 M = 1,9 I = I + 1 ML(I) = MAT(M+L) 20 CONTINUE 30 CONTINUE END IF 40 RETURN C C END C C C SUBROUTINE MATCOP(A,B) C ====================== C C C C C .. Array Arguments .. REAL A(3,3),B(3,3) C .. C .. Local Scalars .. INTEGER I,J C .. SAVE C C DO 20 I = 1,3 DO 10 J = 1,3 B(I,J) = A(I,J) 10 CONTINUE 20 CONTINUE C C END C C C SUBROUTINE MATMUL(A,B,C) C ======================== C C C C C .. Array Arguments .. REAL A(3,3),B(3,3),C(3,3) C .. C .. Local Scalars .. INTEGER I,J C .. SAVE C C DO 20 I = 1,3 DO 10 J = 1,3 C(I,J) = A(I,1)*B(1,J) + A(I,2)*B(2,J) + A(I,3)*B(3,J) 10 CONTINUE 20 CONTINUE C C END C C C SUBROUTINE ORDR3(KNTRL,Q,NQ) C ============================ C C C C Find NQ(1) so that |Q(NQ(1)| is biggest Q. C If KNTRL=1, make NQ(1),NQ(2),NQ(3) a positive permutn. C If KNTRL= anything else, C find nos. NQ() so that|Q(NQ(1)|>|Q(NQ(2))|>|Q(NQ(3))| C C C C C C---- Find no. NQ(1) of Q() with largest |Q()|. C C .. Scalar Arguments .. INTEGER KNTRL C .. C .. Array Arguments .. REAL Q(3) INTEGER NQ(3) C .. C .. Local Scalars .. INTEGER J C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD C .. SAVE C C NQ(1) = 1 IF (ABS(Q(2)).GT.ABS(Q(1))) NQ(1) = 2 IF (ABS(Q(3)).GT.ABS(Q(NQ(1)))) NQ(1) = 3 C C---- Find nos. for other two Q()'s & order so |Q(NQ(2))|>|Q(NQ(3))| C NQ(2) = MOD(NQ(1),3) + 1 NQ(3) = MOD(NQ(2),3) + 1 IF (KNTRL.NE.1) THEN C C IF (ABS(Q(NQ(3))).GT.ABS(Q(NQ(2)))) THEN J = NQ(3) NQ(3) = NQ(2) NQ(2) = J END IF END IF RETURN C C END C C C SUBROUTINE ORIENT(IGROUP,CELL,NV,V,MXB,ESDB,UB,NB) C ================================================== C C C C C---- Find orientation of crystal from a list of reciprocal lattice C vectors if space group and unit cell dimensions are known C C W.KABSCH VERSION 3-1987 C C C C---- Subroutines required: ABSHKL,BASIDX,INVCEL,LAUESY, C MATMUL,RELIX,SETMAT,ULIST C C C C IGROUP - space group number of crystal (given) C CELL - unit cell parameters (angstroem and degrees) (given) C NV - number of reciprocal lattice vectors (given) C V - real array(3,nv) of vector coordinates (given) C MXB - maximum number of possible orientations (given) C ESDB - real array(mxb) of fit-values for each orientation (result) C solutions are sorted in decreasing quality. C UB - real array(3,3,mxb) of orientation matrices (result) C NB - >0: number of possible orientations (result) C 0: no solution found C -1: illegal unit cell parameter C -2: illegal space group number C -3: number of linear independent lattice vectors C is less than two. C C C C C---- get constraint number of general systematic absences C C .. Parameters .. INTEGER NVX,NNVX,NVX3,NXL,MXBEST PARAMETER (NVX=40,NNVX= (NVX*NVX+NVX)/2,NVX3=NVX*3,NXL=15, + MXBEST=6) C .. C .. Scalar Arguments .. INTEGER IGROUP,MXB,NB,NV C .. C .. Array Arguments .. REAL CELL(6),ESDB(MXB),UB(3,3,MXB),V(3,*) C .. C .. Local Scalars .. REAL CUT,RRMAX,RRMIN,X,Y,Z INTEGER ABSTYP,I,IBEST,IER,IRANK,J,K,L,M,MBEST,N,NBEST,NEQL,NL C .. C .. Local Arrays .. REAL A(3,3),AA(3,3),BASIS(3,3),COEF(NNVX),ESD(MXBEST),R(NVX3), + RCELL(6),T(3,NVX),U(3,3),UBEST(3,3,MXBEST),WEIGHT(NVX), + XCALC(3,NVX),XOBS(3,NVX) INTEGER IP(2,NXL),MLAUE(216) INTEGER*2 IH(NVX),IK(NVX),IL(NVX) C .. C .. External Subroutines .. EXTERNAL ABSHKL,BASIDX,INVCEL,LAUESY,MATMUL,RELIX,SETMAT,ULIST C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN,SQRT C .. SAVE C C C ************************** CALL ABSHKL(IGROUP,CELL,ABSTYP) C ************************** C C---- get reciprocal cell parameters C C ******************** CALL INVCEL(CELL,RCELL,X) C ******************** C NB = -1 IF (X.GT.0.0) THEN C C---- get matrix representation of reciprocal cell C C *************** CALL SETMAT(RCELL,A) C *************** C C---- get operators of laue-group C C ****************************** CALL LAUESY(CELL,IGROUP,NEQL,MLAUE) C ****************************** C NB = -2 IF ((NEQL.GE.1) .AND. (NEQL.LE.24)) THEN C C---- set maximum distance between reciprocal lattice vectors C for use by the local indexing method C RRMAX = (MAX(RCELL(1),RCELL(2),RCELL(3))*2.5)**2 C C---- set minimum length of an acceptable reciprocal lattice vector C RRMIN = (MIN(RCELL(1),RCELL(2),RCELL(3))*0.6)**2 C C---- set maximum allowed rms-error for a good interpretation C of the difference vector set C CUT = MIN(RCELL(1),RCELL(2),RCELL(3))*0.3 C C---- generate a list of linearly independent vector pairs from C input list of reciprocal lattice vectors. C NB = -3 IF (NV.GE.2) THEN NL = 0 C C DO 40 I = 1,NV - 1 X = 0.0 C C DO 10 J = 1,3 X = V(J,I)**2 + X BASIS(J,1) = V(J,I) 10 CONTINUE C C IF (X.GT.RRMIN) THEN C C DO 30 K = I + 1,NV Y = 0.0 Z = 0.0 C C DO 20 J = 1,3 Y = V(J,K)**2 + Y Z = V(J,I)*V(J,K) + Z BASIS(J,2) = V(J,K) 20 CONTINUE C C IF (Y.GT.RRMIN) THEN Z = Z/SQRT(X*Y) IF (ABS(Z).LE.0.7) THEN NL = NL + 1 IP(1,NL) = I IP(2,NL) = K IF (NL.EQ.NXL) GO TO 50 END IF END IF 30 CONTINUE END IF 40 CONTINUE C C 50 IF (NL.GE.1) THEN NB = 0 IRANK = 2 N = MIN(NVX,NV) C C DO 130 L = 1,NL C C---- assign indices to a pair of linear independent vectors C and obtain an initial orientation matrix C DO 70 I = 1,2 K = IP(I,L) C C DO 60 J = 1,3 BASIS(J,I) = V(J,K) 60 CONTINUE 70 CONTINUE C C ************************************************* CALL BASIDX(ABSTYP,NEQL,MLAUE,RCELL,BASIS,IRANK,MXBEST, + NBEST,ESD,UBEST) C ************************************************* C IF (NBEST.NE.0) THEN MBEST = ABS(NBEST) C C---- assign indices to the remaining vectors by the C local indexing method C DO 120 IBEST = 1,MBEST DO 90 J = 1,3 DO 80 K = 1,3 U(J,K) = UBEST(J,K,IBEST) 80 CONTINUE 90 CONTINUE C C ********************************************** CALL RELIX(ABSTYP,RRMAX,A,N,V,XOBS,XCALC,WEIGHT,IH, + IK,IL,T,R,COEF,U,IER) C ********************************************** C IF (IER.EQ.0) THEN C C ************** CALL MATMUL(U,A,AA) C ************** C Z = 0.0 C C DO 110 J = 1,N DO 100 K = 1,3 Z = (V(K,J)-AA(K,1)*IH(J)-AA(K,2)*IK(J)- + AA(K,3)*IL(J))**2 + Z 100 CONTINUE 110 CONTINUE C C Z = SQRT(Z/N) C C---- save assignment if among the best "mxb" independent choices C C ************************************** CALL ULIST(NEQL,MLAUE,A,U,Z,MXB,NB,UB,ESDB) C ************************************** C END IF 120 CONTINUE END IF 130 CONTINUE C C---- find number of solutions which must be considered further C IF (NB.GE.1) THEN M = NB + 1 140 CONTINUE M = M - 1 IF ((M.GT.1) .AND. (ESDB(M).GT. (2.0*ESDB(1)))) NB = M - + 1 IF (M.GT.2) GO TO 140 END IF RETURN END IF END IF END IF END IF END C C C SUBROUTINE POLANG(DC,OMEGA,PHI) C =============================== C C---- Find polar angles OMEGA & PHI (degrees) from C direction cosines DC C OMEGA in range 0 to 180, PHI in range 0 to 360 C C C C C C .. Scalar Arguments .. REAL OMEGA,PHI C .. C .. Array Arguments .. REAL DC(3) C .. C .. Local Scalars .. REAL CONV,SOMEGA C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN,ATAN2,SQRT C .. SAVE C C CONV = 45.0/ATAN(1.0) C C SOMEGA = SQRT(ABS(1.0-DC(3)*DC(3))) OMEGA = ATAN2(SOMEGA,DC(3))*CONV C C IF (DC(1).EQ.0.0 .AND. DC(2).EQ.0.0) THEN PHI = 0.0 ELSE PHI = ATAN2(DC(2),DC(1))*CONV END IF C C END C C C SUBROUTINE PXTOMM(Q,F,ED,QED) C ============================= C C C C C .. Scalar Arguments .. REAL F,Q C .. C .. Array Arguments .. REAL ED(3,3),QED(3,3) C .. C .. Local Scalars .. INTEGER I C .. SAVE C C DO 10 I = 1,3 QED(I,1) = ED(I,1)*Q QED(I,2) = ED(I,2)*Q QED(I,3) = ED(I,3)*F 10 CONTINUE C C END C C C SUBROUTINE RECCEL(R,C,WAVE) C =========================== C C C---- Convert between real and reciprocal C cell parameters. C input, R output C Angles in degrees C C C C .. Scalar Arguments .. REAL WAVE C .. C .. Array Arguments .. REAL C(6),R(6) C .. C .. Local Scalars .. REAL SUM,V,WV INTEGER I C .. C .. External Functions .. REAL ACOSD,COSD,SIND cc EXTERNAL ACOSD,COSD,SIND C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. SAVE C C SUM = 0.0 C C DO 10 I = 4,6 SUM = C(I)*0.5 + SUM 10 CONTINUE C C V = 1.0 C C DO 20 I = 4,6 V = SIND(SUM-C(I))*V 20 CONTINUE C C WV = WAVE/ (C(1)*2.0*C(2)*C(3)*SQRT(SIND(SUM)*V)) R(1) = C(2)*C(3)*SIND(C(4))*WV R(2) = C(1)*C(3)*SIND(C(5))*WV R(3) = C(1)*C(2)*SIND(C(6))*WV C C R(4) = ACOSD((COSD(C(5))*COSD(C(6))-COSD(C(4)))/ + (SIND(C(5))*SIND(C(6)))) R(5) = ACOSD((COSD(C(6))*COSD(C(4))-COSD(C(5)))/ + (SIND(C(6))*SIND(C(4)))) R(6) = ACOSD((COSD(C(4))*COSD(C(5))-COSD(C(6)))/ + (SIND(C(4))*SIND(C(5)))) C C END C C C SUBROUTINE REFINE(NCYCLE,N,IH,IK,IL,IX,IY,IPHI,ACHSE,ICS,Q,ORGX, + ORGY,F,S0L,ED,U,RCELL,SDU,SDCELL,SDPHI,SDXY,IER) C ================================================================ C C C---- Constrained/restrained least-squares refinement C Subroutine for detector and crystal parameters C WOLFGANG KABSCH revised version 11-1987 C C C C C Subroutines required: GONSYS,MATCOP,MATMUL,REFLEX,SOLVER,UNORM C C C C NCYCLE - number of refinement cycles (given) C N - number of reflections (given) C IH - integer*2 array specifying h-indices for a list (given) C of n reflections. C IK - integer*2 array specifying k-indices for a list (given) C of n reflections. C IL - integer*2 array specifying l-indices for a list (given) C of n reflections. C IX - integer*2 array specifying x-position on detector (given) C for each spot in the list ( tenth of a pixel) C IY - integer*2 array specifying y-position on detector (given) C for each spot in the list ( tenth of a pixel) C IPHI - integer*2 array specifying angular position of (given) C spindle where reflection with indices ih,ik,il C was diffracting. the values are given in units C of a hundreth of a degree. C ACHSE - array(3) specifying lab coordinates of rotation axis (given) C ICS - number specifying refinement constraints (given) C -3: refine origin C -2: refine origin and distance between detector and C crystal. C -1: refine distance between detector and crystal C as well as origin and orientation of the detector C 0: refine direct beam direction,unit cell orientation C and distance between detector and crystal C >0: refine direct beam direction,unit cell orientation C and independent reciprocal cell parameters. C these independent parameters are derived from the C exact value of "ics". C 1: triclinic C 2: monoclinic first setting C 3: monoclinic second setting C 4: orthorhombic C 5: tetragonal C 6: trigonal C 7: hexagonal C 8: cubic C Q - length of a pixel in millimeters (given) C ORGX - x- and y-coordinates (pixels) on detector such (updated) C ORGY - that the detector normal would intersect the (updated) C crystal C F - crystal to detector distance (mm). (updated) C S0L - coordinates of incident x-ray beam wave vector. (updated) C only the direction of s0l is refined. C length of s0l is 1.0/lambda. ( reciprocal angstroem ) C ED - real array(3,3) specifying lab coordinates of (updated) C detector axes. C U - array(3,3) containing orientation matrix (updated) C RCELL - array of length 6 of reciprocal unit cell parameters(updated) C in reciprocal angstroem and degrees. C SDU - estimated standard deviation (degrees) (result) C for orientation matrix u C SDCELL - array(6) of standard deviations (result) C ( reciprocal angstroem and degrees ) for C reciprocal unit cell parameters. C SDPHI - standard deviation of angular position of spindle (updated) C ( degrees ). a positive input value is C used to calculate the weights of the observed C spindle positions of the reflections. C a non-positive value indicates that a default C value of 0.1 degrees should be used. C on return, sdphi contains the estimated standard C deviation of the spindle angle at which the C reflections were diffracting. C SDXY - standard deviation of spot position on detector (updated) C (pixels). a positive input value is C used to calculate the weights of the observed C spot positions of the reflections. C a non-positive value indicates that a default C value of 3 pixels should be used. C on return, sdxy contains the estimated standard C deviation of the spot positions on the detector. C IER - error indicator set by subroutine . normally ier is (result) C number of accepted reflections. C >0 : no error. C 0 : insufficient number of accepted reflections C -1 : q is not positive C -2 : f is zero or illegal detector orientation C -3 : illegal cell parameters C -4 : illegal rotation axis or direct beam wave vector C -5 : orientation matrix u is singular C -6 : ncycle must be between 1 and 10 C C C .. Parameters .. REAL RAD,RAD100 PARAMETER (RAD=57.29578,RAD100=100.0*RAD) C .. C .. Scalar Arguments .. REAL F,ORGX,ORGY,Q,SDPHI,SDU,SDXY INTEGER ICS,IER,N,NCYCLE C .. C .. Array Arguments .. REAL ACHSE(3),ED(3,3),RCELL(6),S0L(3),SDCELL(6),U(3,3) INTEGER*2 IH(*),IK(*),IL(*),IPHI(*),IX(*),IY(*) C .. C .. Local Scalars .. REAL CHIQ,CPHI,CPHIC,DELPHI,DELX,DELY,DRHOQ,DRR,ESDPHI,ESDXY,FH, + FK,FL,R,RCPHI,RHO,RHOQ,RR,RSPHI,SIGMA,SPHI,SPHIC,SS,T1,T2,T3, + T4,WPHI,WXY,XOBS,YOBS INTEGER I,ICYCLE,J,M,MFREE,MM,MY1,MY2,NACC,NY C .. C .. Local Arrays .. REAL A(3,3),A0(3,3),B(6),DA(3,3,6),DCPHIC(11),DED(3,3,3), + DPD(3,11),DS0L(3,2),DSD(3,11),DSL(3,11),DSPHIC(11), + DX0G(3,11),DX0L(3,11),DXG(3,11),DXL(3,11),E(121),EG(3,3), + ESD(11),EVAL(11),G(11),H(66),PD(3),RESULT(11),S0G(3),SD(3), + SDAP(11),SL(3),X0G(3),X0L(3),XG(3),XL(3) INTEGER IP(6,8),IP231(3),IP312(3),MCS(8) C .. C .. External Subroutines .. EXTERNAL GONSYS,MATCOP,MATMUL,REFLEX,SOLVER,UNORM C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN,SQRT C .. SAVE C .. Data statements .. DATA IP/1,2,3,4,5,6,1,2,3,0,0,4,1,2,3,0,4,0,1,2,3,0,0,0,1,1,2,0,0, + 0,1,1,1,2,2,2,1,1,2,0,0,0,1,1,1,0,0,0/ DATA IP231,IP312/2,3,1,3,1,2/ DATA MCS/6,4,4,3,2,2,2,1/ C .. C C SDU = 0.0 C C DO 10 MY1 = 1,6 SDCELL(MY1) = 0.0 10 CONTINUE C C---- get number of independent parameters and set C a priori standard deviations of parameters C M = 0 IF (ICS.EQ.0) THEN MFREE = 6 C C DO 20 MY1 = 1,5 SDAP(MY1) = 0.1/RAD 20 CONTINUE C C SDAP(6) = 0.1 ELSE IF (ICS.GT.0) THEN IF (ICS.LT.9) M = MCS(ICS) MFREE = M + 5 C C DO 30 J = 1,5 SDAP(J) = 0.2/RAD 30 CONTINUE C C IER = -3 I = 0 C C DO 40 MY1 = 1,6 R = RCELL(MY1) IF (R.LE.0.0) THEN GO TO 650 ELSE C C IF (MY1.LE.3) THEN RHO = 0.005*R ELSE IF ((R.LT.1.0) .OR. (R.GT.179.0)) THEN GO TO 650 ELSE RHO = 0.2/RAD END IF C C IF (M.GE.1) THEN MY2 = IP(MY1,ICS) IF (MY2.GT.I) THEN I = MY2 SDAP(MY2+5) = RHO END IF END IF END IF 40 CONTINUE ELSE C C C C DO 50 I = 1,3 SDAP(I) = 0.1 SDAP(I+3) = 0.2/RAD 50 CONTINUE C C IF ((ICS+2).EQ.0) THEN MFREE = 3 ELSE IF ((ICS+2).GT.0) THEN MFREE = 6 ELSE C C MFREE = 2 END IF END IF C C CHIQ = MFREE C C---- check for sufficient number of reflections C IER = 0 IF (N.GT.MFREE) THEN C C---- check q,f C IF (Q.LE.0.0) IER = -1 IF (F.EQ.0.0) IER = -2 IF (IER.GE.0) THEN C C---- get weights for observational equations C R = 0.1 IF (SDPHI.GT.0.0) R = SDPHI WPHI = (RAD/R)**2 IF (ICS.LT.0) WPHI = 0.0 WXY = 1.0/ (3.0*Q)**2 IF (SDXY.GT.0.0) WXY = 1.0/ (SDXY*Q)**2 C C---- renormalize orientation matrix to orthogonal form to C prevent accumulation of rounding errors. C C ************ CALL UNORM(U,IER) C ************ C IF (IER.NE.0) THEN IER = -5 ELSE C C---- check number of refine cycles C IER = -6 IF ((NCYCLE.GE.1) .AND. (NCYCLE.LE.10)) THEN C C---- refinemment loop C DO 620 ICYCLE = 1,NCYCLE ESDPHI = 0.0 ESDXY = 0.0 SIGMA = 0.0 C C---- clear normal matrix and right hand side C MM = 0 C C DO 70 MY1 = 1,MFREE G(MY1) = 0.0 DO 60 MY2 = 1,MY1 MM = MM + 1 H(MM) = 0.0 60 CONTINUE 70 CONTINUE C C---- define goniostat system C C **************************** CALL GONSYS(ACHSE,S0L,S0G,EG,IER) C **************************** C IF (IER.EQ.0) THEN C C---- calculate setting matrix C C C DO 100 I = 1,3 R = RCELL(I+3)/RAD B(I) = COS(R) B(I+3) = SIN(R) C C DO 90 J = 1,3 A0(I,J) = 0.0 C C DO 80 MY1 = 1,6 DA(I,J,MY1) = 0.0 80 CONTINUE 90 CONTINUE 100 CONTINUE C C A0(1,1) = RCELL(1)*B(5) A0(1,2) = (B(3)-B(1)*B(2))*RCELL(2)/B(5) A0(2,2) = SQRT((RCELL(2)*B(4))**2-A0(1,2)**2) A0(3,1) = RCELL(1)*B(2) A0(3,2) = RCELL(2)*B(1) A0(3,3) = RCELL(3) C C ************** CALL MATMUL(U,A0,A) C ************** C C---- calculate derivatives of direct beam wavevector C IF (ICS.GE.0) THEN C C DO 110 I = 1,3 DS0L(I,1) = EG(I,1)*S0G(3) DS0L(I,2) = EG(I,3)*S0G(2) - EG(I,2)*S0G(3) 110 CONTINUE C C IF (M.GE.1) THEN C C---- calculate derivatives of setting matrix with respect to C independent cell parameters C DA(1,1,1) = B(5) DA(3,1,1) = B(2) DA(1,2,2) = A0(1,2)/RCELL(2) DA(2,2,2) = A0(2,2)/RCELL(2) DA(3,2,2) = A0(3,2)/RCELL(2) DA(3,3,3) = 1.0 DA(1,2,4) = RCELL(2)*B(4)*B(2)/B(5) DA(2,2,4) = (A0(3,2)*B(4)-A0(1,2)*DA(1,2,4))/ + A0(2,2) DA(3,2,4) = -RCELL(2)*B(4) DA(1,1,5) = A0(3,1) DA(3,1,5) = -A0(1,1) DA(1,2,5) = A0(3,2) - A0(1,2)*B(2)/B(5) DA(2,2,5) = -A0(1,2)*DA(1,2,5)/A0(2,2) DA(1,2,6) = -RCELL(2)*B(6)/B(5) DA(2,2,6) = -A0(1,2)*DA(1,2,6)/A0(2,2) C C DO 160 I = 1,3 DO 150 J = 1,3 C C DO 120 MY1 = 1,M B(MY1) = 0.0 120 CONTINUE C C DO 130 MY1 = 1,6 MY2 = IP(MY1,ICS) IF (MY2.GT.0) B(MY2) = B(MY2) + DA(I,J,MY1) 130 CONTINUE C C DO 140 MY1 = 1,M DA(I,J,MY1) = B(MY1) 140 CONTINUE 150 CONTINUE 160 CONTINUE C C DO 210 MY1 = 1,M DO 180 I = 1,3 DO 170 J = 1,3 A0(I,J) = U(I,1)*DA(1,J,MY1) + + U(I,2)*DA(2,J,MY1) + + U(I,3)*DA(3,J,MY1) 170 CONTINUE 180 CONTINUE C C DO 200 I = 1,3 DO 190 J = 1,3 DA(I,J,MY1) = A0(I,J) 190 CONTINUE 200 CONTINUE 210 CONTINUE END IF END IF C C---- calculate derivative of detector orientation matrix C IF (ICS.LT.0) THEN C C DO 220 J = 1,3 DED(1,J,1) = 0.0 DED(2,J,1) = -ED(3,J) DED(3,J,1) = ED(2,J) DED(1,J,2) = ED(3,J) DED(2,J,2) = 0.0 DED(3,J,2) = -ED(1,J) DED(1,J,3) = -ED(2,J) DED(2,J,3) = ED(1,J) DED(3,J,3) = 0.0 220 CONTINUE END IF C C---- enter contributions of observational equations C NACC = 0 C C DO 510 NY = 1,N FH = IH(NY) FK = IK(NY) FL = IL(NY) XOBS = (IX(NY)/10.0-ORGX)*Q YOBS = (IY(NY)/10.0-ORGY)*Q R = IPHI(NY)/RAD100 CPHI = COS(R) SPHI = SIN(R) C C---- get laboratory coordinates of reflection at phi=0 C DO 230 I = 1,3 X0L(I) = A(I,1)*FH + A(I,2)*FK + A(I,3)*FL 230 CONTINUE C C---- calculate diffraction geometry for this reflection C C ******************************************** CALL REFLEX(EG,CPHI,SPHI,S0L,S0G,X0L,X0G,XL,XG,SL, + CPHIC,SPHIC,RR,RHOQ,I) C ******************************************** C IF (I.EQ.0) THEN C C---- get detector coordinates of diffracted beam wavevector C DO 240 I = 1,3 SD(I) = ED(1,I)*SL(1) + ED(2,I)*SL(2) + + ED(3,I)*SL(3) 240 CONTINUE C C IF (SD(3).NE.0.0) THEN C C---- get detector coordinates of intersection point of diffracted C beam with detector surface C DO 250 I = 1,3 PD(I) = SD(I)*F/SD(3) 250 CONTINUE C C---- calculate all derivatives C NACC = NACC + 1 IF (ICS.EQ.0) THEN C C---- Get laboratory coordinates of derivatives of "x0l" C change of detector distance C C C DO 260 I = 1,3 DX0L(I,6) = 0.0 260 CONTINUE ELSE IF (ICS.LE.0) THEN C C---- derivatives of calculated detector coordinates with respect to C changes of origin, distance and orientation of the detector. C the calculated diffraction pattern is fixed. C C C DO 280 MY1 = 1,3 DO 270 I = 1,3 DPD(I,MY1+3) = DED(1,I,MY1)*SL(1) + + DED(2,I,MY1)*SL(2) + + DED(3,I,MY1)*SL(3) 270 CONTINUE 280 CONTINUE C C---- Change of detector origin C DPD(1,1) = 1.0 DPD(1,2) = 0.0 DPD(2,1) = 0.0 DPD(2,2) = 1.0 C C DO 300 I = 1,2 C C---- change of detector distance C DPD(I,3) = SD(I)/SD(3) C C---- Change of detector orientation C DO 290 MY1 = 4,6 DPD(I,MY1) = (DPD(I,MY1)- + DPD(3,MY1)*SD(I)/SD(3))*F/ + SD(3) 290 CONTINUE 300 CONTINUE C C GO TO 480 END IF C C---- derivatives with respect to eta and kappa are 0 C C C DO 310 I = 1,3 DX0L(I,1) = 0.0 DX0L(I,2) = 0.0 310 CONTINUE C C---- derivatives with respect to unit cell orientation C DX0L(1,3) = 0.0 DX0L(2,3) = -X0L(3) DX0L(3,3) = X0L(2) DX0L(1,4) = X0L(3) DX0L(2,4) = 0.0 DX0L(3,4) = -X0L(1) DX0L(1,5) = -X0L(2) DX0L(2,5) = X0L(1) DX0L(3,5) = 0.0 IF (M.GE.1) THEN C C---- derivatives with respect to independent cell parameters C DO 330 MY1 = 1,M DO 320 I = 1,3 DX0L(I,MY1+5) = DA(I,1,MY1)*FH + + DA(I,2,MY1)*FK + + DA(I,3,MY1)*FL 320 CONTINUE 330 CONTINUE END IF C C---- Get goniostat coordinates of derivatives of "x0g" C derivatives with respect to eta C DX0G(1,1) = -X0G(3) DX0G(2,1) = 0.0 DX0G(3,1) = X0G(1) C C---- all other derivatives of "x0g" C DO 350 J = 2,MFREE DO 340 I = 1,3 DX0G(I,J) = EG(1,I)*DX0L(1,J) + + EG(2,I)*DX0L(2,J) + + EG(3,I)*DX0L(3,J) 340 CONTINUE 350 CONTINUE C C---- Get goniostat coordinates of derivatives of "xg" C derivatives of spindle component C DO 360 J = 1,MFREE DXG(2,J) = EG(1,2)*DX0L(1,J) + + EG(2,2)*DX0L(2,J) + + EG(3,2)*DX0L(3,J) 360 CONTINUE C C---- derivatives with respect to eta are 0 C DXG(1,1) = 0.0 DXG(3,1) = 0.0 C C---- derivatives with respect to kappa C DXG(3,2) = (X0G(2)*S0G(3)-XG(3)*S0G(2))/S0G(3) DXG(1,2) = -XG(3)*DXG(3,2)/XG(1) C C---- all other derivatives C DO 370 J = 3,MFREE DRR = X0L(1)*DX0L(1,J) + X0L(2)*DX0L(2,J) + + X0L(3)*DX0L(3,J) DXG(3,J) = - (S0G(2)*DX0G(2,J)+DRR)/S0G(3) DXG(1,J) = (DRR-XG(2)*DXG(2,J)- + XG(3)*DXG(3,J))/XG(1) 370 CONTINUE C C---- get laboratory coordinates of derivatives of "xl" C DO 390 I = 1,3 C C---- derivatives with respect to eta C DXL(I,1) = EG(I,1)*XG(3) - EG(I,3)*XG(1) C C---- all other derivatives C DO 380 J = 2,MFREE DXL(I,J) = EG(I,1)*DXG(1,J) + + EG(I,2)*DXG(2,J) + + EG(I,3)*DXG(3,J) 380 CONTINUE 390 CONTINUE C C---- calculate derivatives of sin/cos of spindle position C DO 400 J = 1,MFREE DRHOQ = (X0G(1)*DX0G(1,J)+X0G(3)*DX0G(3,J))* + 2.0 DCPHIC(J) = (DXG(1,J)*X0G(1)+XG(1)*DX0G(1,J)+ + DXG(3,J)*X0G(3)+XG(3)*DX0G(3,J)- + CPHIC*DRHOQ)/RHOQ DSPHIC(J) = (DXG(1,J)*X0G(3)+XG(1)*DX0G(3,J)- + DXG(3,J)*X0G(1)-XG(3)*DX0G(1,J)- + SPHIC*DRHOQ)/RHOQ 400 CONTINUE C C---- get laboratory coordinates of derivatives of "sl" C DO 420 I = 1,3 C C---- derivatives with respect to eta C DSL(I,1) = EG(I,1)*S0G(3) + DXL(I,1) C C---- derivatives with respect to kappa C DSL(I,2) = DXL(I,2) - EG(I,2)*S0G(3) + + EG(I,3)*S0G(2) C C---- all other derivatives C DO 410 J = 3,MFREE DSL(I,J) = DXL(I,J) 410 CONTINUE 420 CONTINUE C C---- get detector coordinates of diffracted beam wavevector C derivatives C DO 440 I = 1,3 DO 430 J = 1,MFREE DSD(I,J) = ED(1,I)*DSL(1,J) + + ED(2,I)*DSL(2,J) + + ED(3,I)*DSL(3,J) 430 CONTINUE 440 CONTINUE C C---- get detector coordinates of derivatives of intersection point C of diffracted beam with detector surface C DO 460 I = 1,3 DO 450 J = 1,MFREE DPD(I,J) = (DSD(I,J)*F-PD(I)*DSD(3,J))/SD(3) 450 CONTINUE 460 CONTINUE C C IF (ICS.EQ.0) THEN C C DO 470 I = 1,3 DPD(I,MFREE) = PD(I)/F 470 CONTINUE END IF C C---- calculate residuals C 480 RCPHI = CPHIC - CPHI RSPHI = SPHIC - SPHI DELX = PD(1) - XOBS DELY = PD(2) - YOBS C C---- add contributions to determine root-mean-square of residuals C DELPHI = (RCPHI**2+RSPHI**2)*RHOQ/RR ESDPHI = ESDPHI + DELPHI ESDXY = DELX**2 + ESDXY + DELY**2 SIGMA = WPHI*DELPHI + SIGMA + + (DELX**2+DELY**2)*WXY C C---- add contribution of this reflection to normal equations C MM = 0 C C DO 500 MY1 = 1,MFREE T1 = DSPHIC(MY1)*WPHI*RHOQ/RR T2 = DCPHIC(MY1)*WPHI*RHOQ/RR T3 = DPD(1,MY1)*WXY T4 = DPD(2,MY1)*WXY G(MY1) = G(MY1) - RSPHI*T1 - RCPHI*T2 - + DELX*T3 - DELY*T4 DO 490 MY2 = 1,MY1 MM = MM + 1 H(MM) = DSPHIC(MY2)*T1 + H(MM) + + DCPHIC(MY2)*T2 + DPD(1,MY2)*T3 + + DPD(2,MY2)*T4 490 CONTINUE 500 CONTINUE END IF END IF 510 CONTINUE C C---- get estimated standard deviations of observed reflection C positions and angles C IF (NACC.LE.MFREE) THEN GO TO 650 ELSE SDPHI = SQRT(ESDPHI/ (NACC-MFREE))*RAD SDXY = SQRT(ESDXY/ (NACC-MFREE))/Q SIGMA = SQRT(SIGMA/ (3.0*NACC-MFREE)) C C---- solve normal equations C C ******************************************** CALL SOLVER(H,G,SDAP,E,EVAL,RESULT,MFREE,CHIQ,T1) C ******************************************** C I = 0 C C DO 520 MY1 = 1,MFREE I = I + MY1 ESD(MY1) = H(I) G(MY1) = RESULT(MY1) 520 CONTINUE C C---- apply corrections and calculate estimated standard deviations C for the refined parameters C IF (ICS.GE.0) THEN C C---- incident x-ray beam C SS = 0.0 C C DO 530 I = 1,3 SS = S0L(I)**2 + SS S0L(I) = DS0L(I,1)*G(1) + S0L(I) + + DS0L(I,2)*G(2) 530 CONTINUE C C R = SQRT((S0L(1)**2+S0L(2)**2+S0L(3)**2)/SS) C C DO 540 I = 1,3 S0L(I) = S0L(I)/R 540 CONTINUE C C---- orientation matrix C SDU = 0.0 C C DO 550 I = 1,3 B(I) = G(I+2) SDU = ESD(I+2) + SDU 550 CONTINUE C C SDU = RAD*SIGMA*SQRT(SDU) SPHI = SQRT(B(1)**2+B(2)**2+B(3)**2) IF (SPHI.GT.0.0) THEN C C DO 560 I = 1,3 B(I) = B(I)/SPHI 560 CONTINUE C C IF (SPHI.GT.0.5) SPHI = 0.5 END IF CPHI = SQRT(1.0-SPHI**2) C C DO 580 I = 1,3 DO 570 J = 1,3 R = 0.0 IF (I.EQ.J) R = CPHI A0(I,J) = B(I)*B(J)* (1.0-CPHI) + R 570 CONTINUE 580 CONTINUE C C A0(3,2) = B(1)*SPHI + A0(3,2) A0(2,3) = A0(2,3) - B(1)*SPHI A0(1,3) = B(2)*SPHI + A0(1,3) A0(3,1) = A0(3,1) - B(2)*SPHI A0(2,1) = B(3)*SPHI + A0(2,1) A0(1,2) = A0(1,2) - B(3)*SPHI C C ************** CALL MATMUL(A0,U,A) CALL MATCOP(A,U) C ************** C IF (ICS.EQ.0) F = F + G(6) C C---- reciprocal cell dimensions C IF (M.GE.1) THEN C C DO 590 MY1 = 1,6 MY2 = IP(MY1,ICS) IF (MY2.GT.0) THEN IF (MY1.GT.3) THEN R = G(MY2+5)*RAD IF (R.LT.-3.0) R = -3.0 IF (R.GT.3.0) R = 3.0 R = RCELL(MY1) + R IF (R.LT.1.0) R = 1.0 IF (R.GT.179.0) R = 179.0 RCELL(MY1) = R R = RAD ELSE R = G(MY2+5)/RCELL(MY1) IF (R.LT.-0.03) R = -0.03 IF (R.GT.0.03) R = 0.03 RCELL(MY1) = (1.0+R)*RCELL(MY1) R = 1.0 END IF SDCELL(MY1) = SIGMA*R*SQRT(ESD(MY2+5)) END IF 590 CONTINUE END IF ELSE ORGX = G(1)/Q + ORGX ORGY = G(2)/Q + ORGY IF (ICS.NE.-3) THEN F = G(3) + F IF (ICS.NE.-2) THEN C C DO 610 J = 1,3 DO 600 I = 1,3 ED(I,J) = G(4)*DED(I,J,1) + ED(I,J) + + G(5)*DED(I,J,2) + + G(6)*DED(I,J,3) 600 CONTINUE 610 CONTINUE C C *********** CALL UNORM(ED,I) C *********** C IF (I.NE.0) THEN GO TO 640 ELSE C C END IF END IF END IF END IF END IF ELSE GO TO 630 END IF 620 CONTINUE C C---- end of refinement loop *********** C IER = NACC GO TO 650 630 IER = -4 GO TO 650 640 IER = -2 END IF END IF END IF END IF 650 RETURN END C C C SUBROUTINE REFLEX(EG,CPHI,SPHI,S0L,S0G,X0L,X0G,XL,XG,SL,CPHIC, + SPHIC,RR,RHOQ,IER) C ======================================================== C C C C C---- Calculate diffraction geometry for a reflection C C C C EG - real array(3,3) specifying goniostat system (given) C eg(i,1)=lab coordinates of unit vector along vector C cross product between spindle axis and direct beam. C eg(i,2)=lab coordinates of unit vector along spindle. C eg(i,3)=lab coordinates of eg(i,1) x eg(i,2). C CPHI - cosine of observed spindle position at diffraction (given) C SPHI - sine of observed spindle position at diffraction (given) C S0L - lab coordinates of direct beam wavevector (given) C S0G - coordinates of s0l with respect to eg (given) C X0L - lab coordinates of reflection at phi=0 (given) C X0G - goniostat coordinates of x0l at phi=0 (result) C XL - lab coordinates of reflection at diffraction (result) C XG - goniostat coordinates of reflection at diffraction (result) C SL - lab coordinates of diffracted beam wavevector (result) C CPHIC - cosine of calculated spindle position at diffraction (result) C SPHIC - sine of calculated spindle position at diffraction (result) C RR - square of x0l (result) C RHOQ - square of distance of x0l to spindle axis (result) C IER - result indicator 0: o.k. -1:blind reflection (result) C C C C C C .. Scalar Arguments .. REAL CPHI,CPHIC,RHOQ,RR,SPHI,SPHIC INTEGER IER C .. C .. Array Arguments .. REAL EG(3,3),S0G(3),S0L(3),SL(3),X0G(3),X0L(3),XG(3),XL(3) C .. C .. Local Scalars .. INTEGER I C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. SAVE C C IER = -1 IF (S0G(3).NE.0.0) THEN C C---- Get goniostat coordinates of reflection at phi=0 C DO 10 I = 1,3 X0G(I) = EG(1,I)*X0L(1) + EG(2,I)*X0L(2) + EG(3,I)*X0L(3) 10 CONTINUE C C---- get goniostat coordinates of reflection at diffraction C RR = X0L(1)**2 + X0L(2)**2 + X0L(3)**2 RHOQ = RR - X0G(2)**2 XG(3) = - (X0G(2)*S0G(2)+RR*0.5)/S0G(3) XG(2) = X0G(2) XG(1) = RHOQ - XG(3)**2 IF (XG(1).GT.0.0) THEN XG(1) = SQRT(XG(1)) IF ((X0G(1)*CPHI+X0G(3)*SPHI).LT.0.0) XG(1) = -XG(1) C C---- calculate sin/cos of spindle position at diffraction C CPHIC = (XG(1)*X0G(1)+XG(3)*X0G(3))/RHOQ SPHIC = (XG(1)*X0G(3)-XG(3)*X0G(1))/RHOQ C C DO 20 I = 1,3 C C---- get laboratory coordinates of reflection at diffraction C XL(I) = EG(I,1)*XG(1) + EG(I,2)*XG(2) + EG(I,3)*XG(3) C C---- get laboratory coordinates of diffracted beam wavevector C SL(I) = S0L(I) + XL(I) 20 CONTINUE C C---- return C IER = 0 END IF END IF RETURN END C C C SUBROUTINE RELIX(ABSTYP,RRMIN,A0,NV,V,XOBS,XCALC,WEIGHT,IH,IK,IL, + T,R,C,U,IER) C ================================================================== C C C---- Relative indexing algorithm C C External subroutines required: GELS,GENABS,INVERS,MATMUL,U3BEST C C C C ABSTYP - an integer value describing conditions limiting (given) C possible reflections C RRMIN - maximum value of squared distance between (given) C connected reciprocal lattice points C A0 - setting matrix of reciprocal cell (given) C NV - number of reciprocal lattice points (given) C V - array(3,nv) containing coordinates of reciprocal (given) C lattice points C XOBS, - scratch array(3,nv) (modified) C XCALC C WEIGHT - scratch array(nv) (modified) C IH, - integer*2 array(nv) containing indices (result) C IK, C IL C T - scratch array(3,nv) (modified) C R - scratch array(3*nv) (modified) C C - scratch array((nv*nv+nv)/2) (modified) C U - refined orientation matrix (updated) C IER - >0: warning message from subroutine gels (result) C 0: no error C -1: error in using subroutine gels C -2: nv<2 or illegal setting matrix C -3: error in subroutine u3best C C C C C C .. Parameters .. REAL TWOPI PARAMETER (TWOPI=6.283185) C .. C .. Scalar Arguments .. REAL RRMIN INTEGER ABSTYP,IER,NV C .. C .. Array Arguments .. REAL A0(3,3),C(*),R(*),T(3,*),U(3,3),V(3,*),WEIGHT(*),XCALC(3,*), + XOBS(3,*) INTEGER*2 IH(*),IK(*),IL(*) C .. C .. Local Scalars .. REAL D,RMS,RR,X,Y,Z INTEGER CHANGE,I,I1,I2,I3,J,J1,J2,J3,K,L,M,M1,M2,M3,NV2,NV3,NVV C .. C .. Local Arrays .. REAL A(3,3),B(3,3),FHKL(6),TT(3) INTEGER IHKL(6),MHKL(3) C .. C .. External Subroutines .. EXTERNAL GELS,GENABS,INVERS,MATMUL,U3BEST C .. C .. Intrinsic Functions .. INTRINSIC ABS,ATAN2,COS,NINT,SIGN,SIN C .. C .. Equivalences .. EQUIVALENCE (M1,MHKL(1)), (M2,MHKL(2)), (M3,MHKL(3)) C .. SAVE C C IF (NV.GE.2) THEN C C ************** CALL MATMUL(U,A0,A) CALL INVERS(A,B,D) C ************** C IF (D.GT.0.0) THEN C C---- set-up linear system of equations for the unknown indices C NV2 = NV + NV NV3 = NV + NV2 C C DO 10 J = 1,NV3 R(J) = 0.0 10 CONTINUE C C NVV = (NV*NV+NV)/2 C C DO 20 J = 1,NVV C(J) = 0.0 20 CONTINUE C C---- get all short difference vectors C DO 100 J = 1,NV - 1 DO 90 K = J + 1,NV X = V(1,J) - V(1,K) Y = V(2,J) - V(2,K) Z = V(3,J) - V(3,K) RR = X*X + Y*Y + Z*Z IF (RR.LE.RRMIN) THEN C C---- find indices for this short difference vector C DO 30 L = 1,3 D = B(L,1)*X + B(L,2)*Y + B(L,3)*Z M = SIGN(0.5,D) + D IHKL(L) = M FHKL(L) = ABS(D-M) IF (D.LT.M) M = M - 2 IHKL(L+3) = M + 1 FHKL(L+3) = 1.0 - FHKL(L) 30 CONTINUE C C D = 100.0 M = 0 C C DO 60 I1 = 1,6,3 J1 = IHKL(I1) DO 50 I2 = 2,6,3 J2 = IHKL(I2) DO 40 I3 = 3,6,3 J3 = IHKL(I3) M = M + 1 C C---- check if indices are allowed C C ************************* CALL GENABS(ABSTYP,J1,J2,J3,I) C ************************* C IF (I.EQ.0) THEN Z = FHKL(I1) + FHKL(I2) + FHKL(I3) IF (Z.LE.D) THEN D = Z M1 = J1 M2 = J2 M3 = J3 IF (M.EQ.1) GO TO 70 END IF END IF 40 CONTINUE 50 CONTINUE 60 CONTINUE C C---- enter contribution to linear system of equations C 70 I = 0 C C DO 80 L = 1,3 M = MHKL(L) R(J+I) = R(J+I) + M R(K+I) = R(K+I) - M I = I + NV 80 CONTINUE C C I = (K*K-K)/2 + J C(I) = C(I) - 1.0 I = (K*K+K)/2 C(I) = C(I) + 1.0 I = (J*J+J)/2 C(I) = C(I) + 1.0 END IF 90 CONTINUE 100 CONTINUE C C---- partition reflections into equivalence classes C DO 110 J = 1,NV IK(J) = 0 IL(J) = J 110 CONTINUE 120 CONTINUE C C CHANGE = 0 L = 1 C C DO 140 K = 2,NV DO 130 J = 1,K - 1 IF (C(L+J).EQ.-1.0) THEN C C IF ((IL(J)-IL(K)).NE.0) THEN IF ((IL(J)-IL(K)).GT.0) THEN IL(J) = IL(K) ELSE C C IL(K) = IL(J) END IF CHANGE = 1 END IF END IF 130 CONTINUE L = L + K 140 CONTINUE IF (CHANGE.EQ.1) GO TO 120 C C---- Complete system of equations such that the sum of indices C for each equivalence class is zero. C L = 0 C C DO 160 K = 1,NV J = IL(K) IK(J) = IK(J) + 1 C C DO 150 J = 1,K IF (IL(K).EQ.IL(J)) C(L+J) = C(L+J) + 1.0 150 CONTINUE C C L = L + K 160 CONTINUE C C---- Solve system of equations C C ********************************** CALL GELS(R,C,NV,3,0.000001,IER,WEIGHT) C ********************************** C IF (IER.EQ.0) THEN C C---- Determine centroids for each equivalence class C DO 190 J = 1,NV L = IL(J) IF (L.EQ.J) THEN C C DO 170 K = 1,3 T(K,L) = V(K,J) 170 CONTINUE ELSE C C DO 180 K = 1,3 T(K,L) = T(K,L) + V(K,J) 180 CONTINUE C C END IF 190 CONTINUE C C---- Determine best rotation between crystal and laboratory system C DO 220 J = 1,NV L = IL(J) IF (J.EQ.L) THEN C C DO 200 K = 1,3 T(K,L) = T(K,L)/IK(L) 200 CONTINUE END IF C C WEIGHT(J) = 1.0 C C DO 210 K = 1,3 XOBS(K,J) = V(K,J) - T(K,L) XCALC(K,J) = R(J+NV)*A0(K,2) + A0(K,1)*R(J) + + R(J+NV2)*A0(K,3) 210 CONTINUE 220 CONTINUE C C ******************************************* CALL U3BEST(WEIGHT,XCALC,XOBS,NV,1,RMS,U,TT,IER) C ******************************************* C IF (IER.NE.0) THEN IER = -3 ELSE C C---- Add fractional indices of centroids to reflection indices C C ************** CALL MATMUL(U,A0,A) CALL INVERS(A,B,D) C ************** C DO 260 J = 1,NV L = IL(J) IF (L.EQ.J) THEN C C DO 230 K = 1,3 TT(K) = B(K,1)*T(1,L) + B(K,2)*T(2,L) + + B(K,3)*T(3,L) 230 CONTINUE C C DO 240 K = 1,3 T(K,L) = TT(K) 240 CONTINUE END IF C C I = 0 C C DO 250 K = 1,3 R(J+I) = R(J+I) + T(K,L) I = I + NV 250 CONTINUE 260 CONTINUE C C---- Remove possible offset from real numbered reflection indices C I = 0 C C DO 290 K = 1,3 X = 0.0 Y = 0.0 C C DO 270 J = 1,NV Z = R(J+I)*TWOPI X = COS(Z) + X Y = SIN(Z) + Y 270 CONTINUE C C IF ((X.NE.0.0) .OR. (Y.NE.0.0)) X = ATAN2(Y,X)/TWOPI TT(K) = X C C DO 280 J = 1,NV R(J+I) = R(J+I) - X 280 CONTINUE C C I = I + NV 290 CONTINUE C C I = 0 C C DO 330 L = 1,3 DO 320 J = 1,NV IF (IL(J).EQ.J) THEN X = 0.0 Y = 0.0 C C DO 300 K = J,NV IF (IL(K).EQ.IL(J)) THEN Z = R(K+I)*TWOPI X = COS(Z) + X Y = SIN(Z) + Y END IF 300 CONTINUE C C IF ((X.NE.0.0) .OR. (Y.NE.0.0)) X = ATAN2(Y,X)/TWOPI C C DO 310 K = J,NV IF (IL(K).EQ.IL(J)) R(K+I) = R(K+I) - X 310 CONTINUE END IF 320 CONTINUE C C I = I + NV 330 CONTINUE C C---- Index all reflections C DO 340 J = 1,NV IH(J) = NINT(R(J)) IK(J) = NINT(R(J+NV)) IL(J) = NINT(R(J+NV2)) 340 CONTINUE C C END IF END IF RETURN END IF END IF C C---- Error exit C IER = -2 END C C C SUBROUTINE RNGANG(A) C ==================== C C C---- Put angle A into range 0-360 degrees C C C C C C C .. Scalar Arguments .. REAL A C .. SAVE C C 10 CONTINUE IF (A.LT.0.0) THEN A = A + 360.0 GO TO 10 END IF 20 CONTINUE C C IF (A.GE.360.0) THEN A = A - 360.0 GO TO 20 END IF C C END C C C SUBROUTINE RQSORT(N,IX,IY,IPHI,XORG,YORG) C ======================================== C C C C C---- sorts arrays ix,iy,iphi into increasing distance to film origin C non-recursive version of quicksort algorithm (c.a.r. hoare) C C N - number of items to be sorted (given) C IX - integer*2 array containing x-coordinates (updated) C ( units are "pixel"/10 ) C IY - integer*2 array containing y-coordinates (updated) C ( units are "pixel"/10 ) C IPHI - integer*2 array containing angular positions of (updated) C spindle ( units are degrees/100 ) C XORG - x-coordinate of origin ( units in "pixel" ) (given) C YORG - y-coordinate of origin ( units in "pixel" ) (given) C C C C C .. Scalar Arguments .. REAL XORG,YORG INTEGER N C .. C .. Array Arguments .. INTEGER*2 IPHI(*),IX(*),IY(*) C .. C .. Local Scalars .. REAL ORGX,ORGY,WK,WL,WR INTEGER IL,IR,J,JL,JR,K C .. C .. Local Arrays .. INTEGER ISTACK(32,2) C .. SAVE C C IF (N.GE.2) THEN ORGX = 10.0*XORG ORGY = 10.0*YORG J = 1 ISTACK(1,1) = 1 ISTACK(1,2) = N 10 CONTINUE C C---- take top request from stack C JL = ISTACK(J,1) JR = ISTACK(J,2) J = J - 1 20 CONTINUE C C---- split jl...jr C IL = JL IR = JR K = (JL+JR)/2 WK = (IX(K)-ORGX)**2 + (IY(K)-ORGY)**2 30 CONTINUE WL = (IX(IL)-ORGX)**2 + (IY(IL)-ORGY)**2 IF (WL.LT.WK) THEN IL = IL + 1 IF (IL.LE.JR) THEN GO TO 30 ELSE IL = JR END IF END IF 40 CONTINUE WR = (IX(IR)-ORGX)**2 + (IY(IR)-ORGY)**2 IF (WK.GE.WR) THEN GO TO 50 ELSE IR = IR - 1 IF (IR.GE.JL) GO TO 40 END IF IR = JL 50 IF (IL.LE.IR) THEN K = IX(IL) IX(IL) = IX(IR) IX(IR) = K K = IY(IL) IY(IL) = IY(IR) IY(IR) = K K = IPHI(IL) IPHI(IL) = IPHI(IR) IPHI(IR) = K IL = IL + 1 IR = IR - 1 IF (IL.LE.IR) GO TO 30 END IF IF ((IR-JL).GE. (JR-IL)) THEN C C---- stack request for sorting left partition C IF (JL.LT.IR) THEN J = J + 1 IF (J.GT.32) THEN GO TO 70 ELSE ISTACK(J,1) = JL ISTACK(J,2) = IR END IF END IF C C---- continue sorting right partition C JL = IL ELSE C C---- stack request for sorting right partition C IF (IL.LT.JR) THEN J = J + 1 IF (J.GT.32) THEN GO TO 60 ELSE ISTACK(J,1) = IL ISTACK(J,2) = JR END IF END IF C C---- continue sorting left partition C JR = IR END IF IF (JL.LT.JR) GO TO 20 IF (J.GT.0) THEN GO TO 10 ELSE RETURN END IF 60 STOP ' Stack overflow in "RQSORT" ' 70 STOP ' Stack overflow in "RQSORT" ' END IF END C C C SUBROUTINE RQSORT2(N,V,INTY) C =========================== C C C C C---- Non-recursive version of quicksort algorithm (c.a.r. hoare) C C N - number of items to be sorted (given) C V - real array(3,*) containing cluster coordinates (updated) C indices. arrays v,inty are sorted to decreasing C cluster occupancies. C INTY - integer*2 array of cluster occupancies (updated) C C C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL V(3,*) INTEGER*2 INTY(*) C .. C .. Local Scalars .. REAL X INTEGER IL,IR,J,JJ,JL,JR,K INTEGER*2 W,WK,WL,WR C .. C .. Local Arrays .. INTEGER ISTACK(32,2) C .. SAVE C C IF (N.GE.2) THEN J = 1 ISTACK(1,1) = 1 ISTACK(1,2) = N 10 CONTINUE C C---- take top request from stack C JL = ISTACK(J,1) JR = ISTACK(J,2) J = J - 1 20 CONTINUE C C---- split jl...jr C IL = JL IR = JR K = (JL+JR)/2 WK = -INTY(K) 30 CONTINUE WL = -INTY(IL) IF (WL.LT.WK) THEN IL = IL + 1 IF (IL.LE.JR) THEN GO TO 30 ELSE IL = JR END IF END IF 40 CONTINUE WR = -INTY(IR) IF (WK.GE.WR) THEN GO TO 50 ELSE IR = IR - 1 IF (IR.GE.JL) GO TO 40 END IF IR = JL 50 IF (IL.LE.IR) THEN W = INTY(IL) INTY(IL) = INTY(IR) INTY(IR) = W C C DO 60 JJ = 1,3 X = V(JJ,IL) V(JJ,IL) = V(JJ,IR) V(JJ,IR) = X 60 CONTINUE C C IL = IL + 1 IR = IR - 1 IF (IL.LE.IR) GO TO 30 END IF IF ((IR-JL).GE. (JR-IL)) THEN C C---- stack request for sorting left partition C IF (JL.LT.IR) THEN J = J + 1 IF (J.GT.32) THEN GO TO 80 ELSE ISTACK(J,1) = JL ISTACK(J,2) = IR END IF END IF C C---- continue sorting right partition C JL = IL ELSE C C---- stack request for sorting right partition C IF (IL.LT.JR) THEN J = J + 1 IF (J.GT.32) THEN GO TO 70 ELSE ISTACK(J,1) = IL ISTACK(J,2) = JR END IF END IF C C---- continue sorting left partition C JR = IR END IF IF (JL.LT.JR) GO TO 20 IF (J.GT.0) THEN GO TO 10 ELSE RETURN END IF 70 STOP ' Stack overflow in "RQSORT2" ' 80 STOP ' Stack overflow in "RQSORT2" ' END IF END C C C SUBROUTINE RTOMIS(R,PHI) C ======================== C C C---- Routine returns missetting angles PHI C for a given rotation matrix R C PHI returned in degrees C C [R] = [phiz] . [phiy] . [phix] C C | cz -sz 0 | | cy 0 sy | | 1 0 0 | C [R] = | sz cz 0 | . | 0 1 0 | . | 0 cx -sx | C | 0 0 1 | |-sy 0 cy | | 0 sx cx | C C C | cz*cy cz*sy*sx-sz*cx cz*sy*cx+sz*sx | C [R] = | sz*cy sz*sy*sx+cz*cx sz*sy*cx-cz*sx | C | -sy cy*sx cy*cx | C C C C C .. Array Arguments .. REAL PHI(3),R(3,3) C .. C .. Local Scalars .. REAL ZILCH C .. C .. External Functions .. REAL ASIND,ATAN2D cc EXTERNAL ASIND,ATAN2D C .. C .. Intrinsic Functions .. INTRINSIC ABS,SIGN C .. SAVE C C ZILCH = 0.0000001 C C IF ((1.0-ABS(R(3,1))).LT.ZILCH) THEN C C---- Special case: sin(phiy)= +- 1 and thus cos(phiy) = 0 C Assume phix=0 C PHI(1) = 0.0 PHI(2) = SIGN(90.0,-R(3,1)) PHI(3) = ATAN2D(-R(1,2),R(2,2)) ELSE C C C---- General case: put phiy in range -90 to 90. C Hence cos(phiy) is positive C PHI(2) = ASIND(-R(3,1)) PHI(3) = ATAN2D(R(2,1),R(1,1)) PHI(1) = ATAN2D(R(3,2),R(3,3)) END IF C C END C C C SUBROUTINE RTUMAT(ASTV,BSTV,CSTV) C ================================= C C C Given the setting matrix found by auto-indexing, this routine C operates the point group matrices on it to find the one which C gives the smallest misetting angles which respect to a defined C U matrix. C C The routine DCOSFD is used to work out the unit vector defining C the axis of rotation and the angle. If we find such vectors for C each matrix, the test matrix which gives the largest dot product C with the axis of the target matrix should have the most similar C rotation axis. But there may be an ambiguity about its sense. C This can be resolved by finding the matrix with the most C similar angle of rotation about the axis. C C C NOTES ON ROBUSTNESS: C We often find two solutions with dot product of the axis vectors C which are identical. But because of rounding errors in REAL*4 C representation these errors are compounded over the matrix operations C and conspire to reduce the identity of these products. Here we allow C a tolerance to the fourth decimal place. For added robustness MATMUL C and INVERS should be double precision. At present single precision C seems adequate but in certain unforseen circumstances the wrong C solution may be chosen. Either reduce 1E-4 below or use double C precision routines from the NAG library.PMcL C C C---- External subroutines required C C C DECOMP,INVERS,MATMUL,MATCOP,LAUSEY,RTOMISSET,DCOSFD,UNORM C C C C cc IMPLICIT NONE C .. Array Arguments .. REAL ASTV(3),BSTV(3),CSTV(3) C .. C .. Scalars in Common .. REAL CCOM,CCX,CCY,DPHI,F,FLAMDA,ORGX,ORGY,Q,THRESH INTEGER FIXF,IGROUP LOGICAL DCOMP,FILM,LCAMC,TARGET INTEGER LP,LINOUT LOGICAL ONLINE C .. C .. Arrays in Common .. REAL ACHSE,CELL,ED,S0,TARMAT C .. C .. Local Scalars .. DOUBLE PRECISION BEST,PROD REAL BESTA,D,TEST,TRACE,TRACEU INTEGER IBEST,IER,J,J1,K,K1,K2,NEQ C .. C .. Local Arrays .. REAL AMAT(3,3),ANG(3,23),BINV(3,3),BMAT(3,3),DUM(3,3),EUTAR(3), + EUTEST(3),I(3,3),L(3,3),LBINV(3,3),LINV(3,3),LUINV(3,3), + PHI(3),RMAT(3,3),UINV(3,3),UMAT(3,3),UTAR(3,3),UTEST(3,3) INTEGER MLAUE(432) C .. C .. External Subroutines .. EXTERNAL DCOSFD,DECOMP,INVERS,LAUESY,MATCOP,MATMUL,RTOMIS,UNORM C .. C .. Intrinsic Functions .. INTRINSIC ABS C .. C .. Common blocks .. COMMON /REFCOM/IGROUP,FLAMDA,CELL(6),F,THRESH,Q,DPHI,ED(3,3), + ACHSE(3),S0(3),ORGX,ORGY,FIXF,TARMAT(3,3),TARGET,DCOMP, + LCAMC,CCX,CCY,CCOM,FILM COMMON /RINOUT/ ONLINE,LP,LINOUT C .. SAVE C .. Data statements .. C C DATA I/1,0,0,0,1,0,0,0,1/ DATA MLAUE/432*0/ DATA BEST/1.0E-12/ C .. C C C DO 10 K = 1,3 AMAT(K,1) = ASTV(K) AMAT(K,2) = BSTV(K) AMAT(K,3) = CSTV(K) 10 CONTINUE C C---- DECOMP decomposes AMAT into UMAT*BMAT C C ********************** CALL DECOMP(AMAT,UMAT,BMAT) C ********************** C C---- If the target matrix was given as an A MATRIX C then decompose it to get UTAR C IF (DCOMP) THEN C C *********************** CALL DECOMP(TARMAT,UTAR,DUM) C *********************** C WRITE (LP,FMT= +'('' U-MATRIX DERIVED FROM INPUT A-MATRIX'',/,3(3F10.5 +,/))') ((UTAR(K,J),J=1,3),K=1,3) IF (ONLINE) WRITE (LINOUT,FMT= +'('' U-MATRIX DERIVED FROM INPUT A-MATRIX'',/,3(3F10.5 +,/))') ((UTAR(K,J),J=1,3),K=1,3) ELSE C C ******************* CALL MATCOP(TARMAT,UTAR) C ******************* C END IF C C ******************* CALL MATCOP(UTAR,TARMAT) CALL INVERS(UTAR,UINV,D) CALL INVERS(BMAT,BINV,D) C ******************* C C---- Get unit vector along axis of rotation (EUTAR) and C angle of rotation about that axis (TRACEU) from the C target setting matrix UTAR C C ************************* CALL DCOSFD(UTAR,EUTAR,TRACEU) C ************************* C C---- Read Permutation Matrices for the Point Group C C ***************************** CALL LAUESY(CELL,IGROUP,NEQ,MLAUE) C ****************************** C C---- Main loop over the point group matrices C J1 = 0 IF (ONLINE) WRITE (LINOUT,FMT='(//,1X, + ''LAUE OPERATOR PHIX PHIY PHIZ'')') WRITE (LP,FMT='(//,1X, + ''LAUE OPERATOR PHIX PHIY PHIZ'')') C C DO 50 K = 1,NEQ C C---- Gather point group matrix number K into L(3,3) C DO 30 K1 = 1,3 DO 20 K2 = 1,3 J1 = J1 + 1 L(K2,K1) = MLAUE(J1) 20 CONTINUE 30 CONTINUE C C---- Calculate the UMAT given this point group matrix. C C ******************** CALL INVERS(L,LINV,D) CALL MATMUL(L,UINV,LUINV) CALL MATMUL(UMAT,LUINV,RMAT) CALL RTOMIS(RMAT,PHI) CALL MATMUL(L,BINV,LBINV) CALL MATMUL(AMAT,LBINV,UTEST) C ************************ C C---- Renormalise UTEST to reduce errors built up by matrix C concatenation. C C **************** CALL UNORM(UTEST,IER) C **************** C C---- Calculate how similar this UMAT is to the target matrix. C C Get unit vector along axis of rotation for the candidate C matrix and the angle of rotation about this axis C C ************************* CALL DCOSFD(UTEST,EUTEST,TRACE) C ************************** C TEST = ABS(TRACE-TRACEU) C C PROD = EUTEST(1)*EUTAR(1) + EUTEST(2)*EUTAR(2) + + EUTEST(3)*EUTAR(3) C C---- If the dot product of the axis vectors are larger than the C best so far, make the current the best. If they are similar C to the fourth decimal place, then take the one describing C the smaller angle of rotation C IF ((PROD-BEST).GT.1.0E-4) THEN BESTA = ABS(TEST) BEST = DABS(PROD) IBEST = K ELSE IF (DABS(PROD-BEST).LE.1.0E-4 .AND. + DABS(PROD-BEST).GT.0.0) THEN IF (TEST.LT.BESTA) THEN BESTA = ABS(TEST) BEST = DABS(PROD) IBEST = K END IF END IF C C DO 40 K1 = 1,3 ANG(K1,K) = PHI(K1) 40 CONTINUE C C IF (ONLINE) WRITE (LINOUT,FMT='(10X,I2,2x,3(f8.2),/)') + K,PHI WRITE (LP,FMT='(10X,I2,2x,3(f8.2),/)') K,PHI 50 CONTINUE C C---- End of Main Loop C C---- Now operate the IBEST Laue group matrix on the C input A matrix C J1 = (IBEST-1)*9 C C DO 70 K1 = 1,3 DO 60 K2 = 1,3 J1 = J1 + 1 L(K2,K1) = MLAUE(J1) 60 CONTINUE 70 CONTINUE C C ****************** CALL MATMUL(AMAT,L,DUM) C ****************** C DO 80 K = 1,3 ASTV(K) = DUM(K,1) BSTV(K) = DUM(K,2) CSTV(K) = DUM(K,3) 80 CONTINUE C C END C C C SUBROUTINE SETMAT(RCELL,A) C ========================== C C C---- Calculate setting matrix from reciprocal unit cell. C C RCELL - reciprocal unit cell parameters (given) C in reciprocal angstroem and degrees. C A - setting matrix in standard orientation (result) C C C C C .. Array Arguments .. REAL A(3,3),RCELL(6) C .. C .. Local Scalars .. REAL ARG INTEGER I C .. C .. Local Arrays .. REAL RC(3),RS(3) C .. C .. Intrinsic Functions .. INTRINSIC COS,SIN,SQRT C .. SAVE C C DO 10 I = 1,3 ARG = RCELL(I+3)/57.29578 RC(I) = COS(ARG) RS(I) = SIN(ARG) 10 CONTINUE C C A(1,1) = RCELL(1)*RS(2) A(1,2) = (RC(3)-RC(1)*RC(2))*RCELL(2)/RS(2) A(1,3) = 0.0 A(2,1) = 0.0 A(2,2) = SQRT((RCELL(2)*RS(1))**2-A(1,2)**2) A(2,3) = 0.0 A(3,1) = RCELL(1)*RC(2) A(3,2) = RCELL(2)*RC(1) A(3,3) = RCELL(3) END C C C SUBROUTINE SETTING(ICS,ABSTYP,A0,RCELL,SDCELL,U,RMSMIN) C ======================================================== C C C---- Find setting of unit cell which agrees best with the original C unit cell parameters. W.KABSCH 11-1987 C C C C---- subroutines required: INVERS,MATMUL,SETMAT,UNORM C C C C ICS - number specifying crystal system (given) C 1 : triclinic C 2 : monoclinic first setting C 3 : monoclinic second setting C 4 : orthorhombic C 5 : tetragonal C 6 : trigonal C 7 : hexagonal C 8 : cubic C ABSTYP - constraint number for general conditions limiting (given) C possible reflections C 0: no conditions on h,k,l C 1: h+k must be even C 2: h+l must be even C 3: k+l must be even C 4: h+k+l must be even C 5: -h+k+l must be a multiple of three C 6: h+k and h+l must be even C A0 - setting matrix of original reciprocal cell with (given) C respect to crystal coordinate system. C RCELL - reciprocal cell parameters (rec.angstroem & degrees)(changed) C on return these cell parameters may have been C replaced by a different setting of the unit cell C which agrees better with the original cell parameters. C SDCELL - standard deviations of rec. cell parameters. (changed) C on return these values may have been changed in C accordance with the reciprocal cell parameters. C U - orientation matrix of reciprocal unit cell. (changed) C if a better setting is found the orientation C matrix will be changed accordingly. C RMSMIN - deviation between original rec. cell parameters (result) C and the rec. cell parameters at the best setting. C C C .. Scalar Arguments .. REAL RMSMIN INTEGER ABSTYP,ICS C .. C .. Array Arguments .. REAL A0(3,3),RCELL(6),SDCELL(6),U(3,3) C .. C .. Local Scalars .. REAL RMS,V INTEGER I,I1,I2,I3,IPERM,J,K,M C .. C .. Local Arrays .. REAL A(3,3),A1(3,3),B(3,3),C(3,3),RCELL1(6),RCELL2(6),SDCEL1(6), + SDCEL2(6) INTEGER I123(3),IP(3,6) C .. C .. External Subroutines .. EXTERNAL INVERS,MATMUL,SETMAT,UNORM C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. C .. Equivalences .. EQUIVALENCE (I123(1),I1), (I123(2),I2), (I123(3),I3) C .. SAVE C .. Data statements .. DATA IP/1,2,3,2,3,1,3,1,2,2,1,3,3,2,1,1,3,2/ C .. C C **************** CALL SETMAT(RCELL,A1) CALL MATMUL(U,A1,A) C **************** C DO 10 J = 1,6 RCELL2(J) = RCELL(J) SDCEL2(J) = SDCELL(J) 10 CONTINUE C C RMSMIN = 0.0 C C DO 30 I = 1,3 DO 20 J = 1,3 RMSMIN = (A0(I,J)-A1(I,J))**2 + RMSMIN 20 CONTINUE 30 CONTINUE C C DO 110 IPERM = 1,6 C C IF ((ICS-2).GE.0) THEN IF ((ICS-2).GT.0) THEN C C IF ((ICS-4).EQ.0) THEN C C---- C-axis must be preserved for orthorhombic c-centered cell C IF ((ABSTYP.EQ.1) .AND. (IP(3,IPERM).NE.3)) THEN GO TO 110 C C---- B-axis must be preserved for orthorhombic b-centered cell C ELSE IF ((ABSTYP.EQ.2) .AND. (IP(2,IPERM).NE.2)) THEN GO TO 110 END IF ELSE IF ((ICS-4).GT.0) THEN GO TO 110 C C---- B-axis must be preserved for monoclinic second setting C ELSE IF (IP(2,IPERM).NE.2) THEN GO TO 110 END IF C C---- C-axis must be preserved for monoclinic first setting C ELSE IF (IP(3,IPERM).NE.3) THEN GO TO 110 END IF END IF C C---- apply +/- signs to axes C C C DO 100 I1 = -1,1,2 DO 90 I2 = -1,1,2 DO 80 I3 = -1,1,2 M = I1*I2*I3 C C---- check that unit cell volume would come out positive C IF ((M.GE.0) .OR. (IPERM.GE.4)) THEN IF ((M.LE.0) .OR. (IPERM.LE.3)) THEN C C DO 40 J = 1,3 K = IP(J,IPERM) RCELL1(K) = RCELL(J) SDCEL1(K) = SDCELL(J) SDCEL1(K+3) = SDCELL(J+3) V = RCELL(J+3) IF (M.NE.I123(J)) V = 180.0 - V RCELL1(K+3) = V 40 CONTINUE C C ***************** CALL SETMAT(RCELL1,A1) C ***************** C RMS = 0.0 C C DO 60 I = 1,3 K = IP(I,IPERM) DO 50 J = 1,3 B(J,K) = I123(I)*A(J,I) RMS = (A0(I,J)-A1(I,J))**2 + RMS 50 CONTINUE 60 CONTINUE C C IF (RMS.LT.RMSMIN) THEN RMSMIN = RMS C C ************** CALL INVERS(A1,C,V) C ************** C IF (V.LE.0.0) THEN GO TO 140 ELSE C C ************* CALL MATMUL(B,C,U) CALL UNORM(U,I) C ************* C IF (I.NE.0) THEN GO TO 130 ELSE C C DO 70 J = 1,6 RCELL2(J) = RCELL1(J) SDCEL2(J) = SDCEL1(J) 70 CONTINUE END IF END IF END IF END IF END IF 80 CONTINUE C C 90 CONTINUE 100 CONTINUE 110 CONTINUE C C RMSMIN = SQRT(RMSMIN) C C DO 120 I = 1,6 RCELL(I) = RCELL2(I) SDCELL(I) = SDCEL2(I) 120 CONTINUE C C RETURN 130 STOP ' Program error #2 in subroutine "SETTING" ' 140 STOP ' Program error #1 in subroutine "SETTING" ' END C C C SUBROUTINE SOLVER(A,B,SD,E,EVAL,S,N,CHIQ,Z) C =========================================== C C C---- Solves a system of normal equations using eigenvalue filtering. C C The solution is compatible with a priori assumptions. C The a priori assumption is zero for the mean and sd(i) C for the standard deviation of parameter #i. C C C Reference: DIAMOND,R.(1966)ACTA CRYST.21,253-266. C C C C Subroutines required: EIGEN C C C C A - upper triangular part of a real symmetric matrix. (modified) C on return "a" is replaced by its filtered invers. C dimension of "a" must be at least n*(n+1)/2. C element i,j must have been stored at address C i+(j*j-j)/2 if i1 required) (given) C CHIQ - the solution must satisfy sum((s(j)/sd(j))**2)<=chiq (given) C usually we choose chiq=n. C Z - actual value of sum((s(j)/sd(j))**2) (result) C C C C C C .. Parameters .. REAL ZERO,TWO,TOL PARAMETER (ZERO=0.0,TWO=2.0,TOL=5.0E-8) C .. C .. Scalar Arguments .. REAL CHIQ,Z INTEGER N C .. C .. Array Arguments .. REAL A(*),B(*),E(*),EVAL(*),S(*),SD(*) C .. C .. Local Scalars .. REAL G,GMAX,GMIN,X INTEGER I,IA,J,L,LL,NMAX C .. C .. External Subroutines .. EXTERNAL EIGEN C .. SAVE C C IF (N.GE.2) THEN C C---- rescale input matrix C LL = 0 C C DO 20 I = 1,N IF (SD(I).LE.ZERO) THEN GO TO 160 ELSE C C DO 10 J = 1,I LL = LL + 1 A(LL) = A(LL)*SD(I)*SD(J) 10 CONTINUE END IF 20 CONTINUE C C---- Get eigenvalues and eigenvectors of "a". C C ************ CALL EIGEN(N,A,E) C ************ C C---- Save eigenvalues C LL = 0 NMAX = N C C DO 30 L = 1,N LL = LL + L X = A(LL) C C IF (X.LE.ZERO) THEN NMAX = NMAX - 1 X = ZERO END IF C C EVAL(L) = X 30 CONTINUE C C---- Find coordinates of right-hand side with respect to eigenvectors C LL = 0 C C DO 50 I = 1,N X = ZERO IF (I.LE.NMAX) THEN C C DO 40 J = 1,N X = E(J+LL)*B(J)*SD(J) + X 40 CONTINUE END IF C C A(I) = X LL = LL + N 50 CONTINUE C C---- Find lagrange multiplier C G = EVAL(1)*0.001 Z = ZERO IF (NMAX.GE.1) THEN C C DO 60 I = 1,NMAX Z = (A(I)/ (EVAL(I)+G))**2 + Z 60 CONTINUE C C IF ((Z-CHIQ).LT.0.0) THEN GMAX = G G = ZERO Z = ZERO C C DO 70 I = 1,NMAX Z = (A(I)/ (EVAL(I)+G))**2 + Z 70 CONTINUE C C IF ((Z-CHIQ).GT.0.0) THEN GMIN = G C C DO 90 L = 1,20 G = (GMIN+GMAX)/TWO IF ((GMAX-GMIN).LE. (TOL*EVAL(1))) THEN GO TO 100 ELSE Z = ZERO C C DO 80 I = 1,NMAX Z = (A(I)/ (EVAL(I)+G))**2 + Z 80 CONTINUE C C IF ((Z-CHIQ).EQ.0.0) THEN GO TO 100 ELSE IF ((Z-CHIQ).GT.0.0) THEN GMIN = G ELSE C C GMAX = G END IF END IF 90 CONTINUE END IF END IF END IF C C---- Calculate filtered inverse C 100 IA = 0 C C DO 130 I = 1,N DO 120 J = 1,I IA = IA + 1 X = ZERO C C IF (NMAX.GE.1) THEN LL = 0 C C DO 110 L = 1,NMAX X = E(I+LL)*E(J+LL)/ (EVAL(L)+G) + X LL = LL + N 110 CONTINUE C C END IF C C A(IA) = SD(I)*X*SD(J) 120 CONTINUE 130 CONTINUE C C---- Solve the normal equations by applying the filtered invers C to right hand side vector. C DO 150 I = 1,N X = ZERO C C DO 140 J = 1,N IF (I.LE.J) IA = I + (J*J-J)/2 IF (J.LT.I) IA = J + (I*I-I)/2 X = A(IA)*B(J) + X 140 CONTINUE C C S(I) = X 150 CONTINUE C C RETURN 160 STOP ' !!! ERROR !!! Illegal a priori estimate in "SOLVER".' END IF END C C C SUBROUTINE SPOTS(NX,NY,NBX,NBY,SNR,CB,CS,IBKG,IFRAME,MXRFL,NRFL0, + NRFL,I2X,I2Y,INTY) C ======================================================= C C C---- spot finding routine C C NX - number of pixels along x in array iframe (given) C NY - number of pixels along y in array iframe (given) C NBX, - width of a spot area is (2*nbx+1)*(2*nby+1). (given) C NBY both numbers must be between 1 and 25. (given) C SNR - minimum value of signal to noise ratio for a spot (given) C to become accepted. C CB - average counts/pixel in background frame "ibkg". (given) C when called with cb<=0.0 then no use is made of "ibkg". C "ibkg" and "iframe" may then refer to the same array. C CS - average counts/pixel in data frame "iframe" (given) C IBKG - integer*2 array(nx,ny) of background counts (given) C IFRAME - integer*2 array(nx,ny) of x-ray data counts (given) C MXRFL - dimension of circular arrays i2x,i2y,inty (given) C NRFL0 - number of spots already found in earlier frames. (given) C NRFL - number of spots found on this frame. (result) C a spot is found if its centroid coincides with the C center of the spot box. the background is estimated C as mean value of the counts at the box boundary. C nrfl=0 indicates calling parameter error or no peaks. C nrfl>mxrfl indicates that too many peaks have been C found. the x,y pixel coordinates and C intensities of the last 'mxrfl' peaks C are saved in arrays i2x,i2y,inty. C I2X - integer*2 array(mxrfl) to keep x-coordinates, (modified) C I2Y - integer*2 array(mxrfl) to keep y-coordinates, (modified) C of frame position (pixels) of spots for this frame. C the x,y- coordinates are multiplied by 10. C INTY - integer*2 array(mxrfl) to keep spot intensities (modified) C C C C .. Parameters .. INTEGER MXMBX,MXMBY PARAMETER (MXMBX=51,MXMBY=51) C .. C .. Scalar Arguments .. REAL CB,CS,SNR INTEGER MXRFL,NBX,NBY,NRFL,NX,NY INTEGER NRFL0 C .. C .. Array Arguments .. INTEGER*2 I2X(MXRFL),I2Y(MXRFL),IBKG(NX,NY),IFRAME(NX,NY), + INTY(MXRFL) C .. C .. Local Scalars .. REAL BKG,FAC,FMBXY,PEAK,SIGNAL,SNRQ,VAR INTEGER I,IX,IX1,IX2,IXC0,IY,IYC0,J,JRFL,JRFL0,JRFL1,MBX,MBY,NBX1, + NBY1 INTEGER I32,IXC,IYC,J32,K32,SUM,SUMB C .. C .. Local Arrays .. INTEGER SUMCOL(MXMBX),SUMROW(MXMBY) C .. C .. Intrinsic Functions .. INTRINSIC ABS,MOD,NINT,REAL C .. SAVE C C NBX1 = NBX + 1 NBY1 = NBY + 1 MBX = NBX + NBX1 MBY = NBY + NBY1 FMBXY = MBX*MBY NRFL = 0 C C---- Check calling parameters C IF ((MXRFL.GE.1) .AND. (NRFL0.GE.0) .AND. (NBX.GE.1) .AND. + (NBY.GE.1) .AND. (MBX.LE.MXMBX) .AND. (MBY.LE.MXMBY) .AND. + (NX.GE.MBX) .AND. (NY.GE.MBY) .AND. (SNR.GT.0.0) .AND. + (CS.GT.0.0)) THEN C C---- Peak finding loop C FAC = REAL(MBX*MBY)/REAL((MBX+MBY)*2) SNRQ = SNR*SNR JRFL0 = 0 IXC0 = -10*MBX IYC0 = -10*MBY C C DO 130 IY = NBY1,NY - NBY JRFL1 = NRFL C C DO 20 IX = 1,MBX J32 = 0 C C DO 10 J = -NBY,NBY I32 = IFRAME(IX,IY+J) IF (I32.LT.0) I32 = I32 + 65536 J32 = J32 + I32 10 CONTINUE C C SUMCOL(IX) = J32 20 CONTINUE SUM = 0 C C DO 40 J = 1,MBY J32 = 0 C C DO 30 IX = 1,MBX I32 = IFRAME(IX,IY-NBY1+J) IF (I32.LT.0) I32 = I32 + 65536 J32 = J32 + I32 30 CONTINUE C C SUMROW(J) = J32 SUM = SUM + J32 40 CONTINUE C C DO 120 IX = NBX1,NX - NBX IX1 = MOD(IX-NBX1,MBX) + 1 IX2 = MOD(IX+NBX-1,MBX) + 1 C C---- Estimate background and check signal to noise ratio C BKG = (SUMROW(1)+SUMROW(MBY)+SUMCOL(IX1)+SUMCOL(IX2))*FAC SIGNAL = SUM - BKG IF (SIGNAL.GT.0.0) THEN VAR = FAC*BKG + SUM IF ((SIGNAL**2).GT. (SNRQ*VAR)) THEN C C---- Check if centroid of distribution coincides with center of box C IYC = 0 C C DO 50 J = 1,NBY IYC = (SUMROW(NBY1+J)-SUMROW(NBY1-J))*J + IYC 50 CONTINUE C C IYC = NINT(10.0*IYC/SIGNAL) IF (ABS(IYC).LE.5) THEN IXC = 0 C C DO 60 I = 1,NBX IX1 = MOD(IX-I-1,MBX) + 1 IX2 = MOD(IX+I-1,MBX) + 1 IXC = (SUMCOL(IX2)-SUMCOL(IX1))*I + IXC 60 CONTINUE C C IXC = NINT(10.0*IXC/SIGNAL) IF (ABS(IXC).LE.5) THEN C C---- Check if this spot has already been recorded C IXC = 10*IX + IXC IYC = 10*IY + IYC IF ((ABS(IXC-IXC0).GE.10) .OR. + (ABS(IYC-IYC0).GE.10)) THEN C C DO 70 JRFL = JRFL0,JRFL1 IF (JRFL.GE.1) THEN I = MOD(NRFL0+JRFL-1,MXRFL) + 1 IF ((ABS(I2X(I)-IXC).LT.10) .AND. + (ABS(I2Y(I)-IYC).LT.10)) GO TO 100 END IF 70 CONTINUE C C---- Use a better background for estimating significance of spots C IF (CB.GT.0.0) THEN SUMB = 0 VAR = 0.0 C C DO 90 J = -NBY,NBY DO 80 I = -NBX,NBX I32 = IBKG(IX+I,IY+J) IF (I32.LE.0) THEN GO TO 100 ELSE VAR = I32*I32 + VAR SUMB = SUMB + I32 END IF 80 CONTINUE 90 CONTINUE C C BKG = SUMB/CB PEAK = SUM/CS SIGNAL = PEAK - BKG IF (SIGNAL.LE.0.0) THEN GO TO 100 ELSE VAR = VAR/ (CB*CB) + PEAK/CS - BKG*BKG/FMBXY IF (VAR.LE.0.0) THEN GO TO 100 ELSE IF ((SIGNAL**2).LE. (VAR*SNRQ)) THEN GO TO 100 END IF END IF END IF C C---- Record a new spot C IXC0 = IXC IYC0 = IYC I = MOD(NRFL0+NRFL,MXRFL) + 1 NRFL = NRFL + 1 I2X(I) = IXC I2Y(I) = IYC J32 = NINT(SIGNAL) IF (J32.GT.32767) J32 = 32767 INTY(I) = J32 END IF END IF END IF END IF END IF 100 IF ((IX+NBX).NE.NX) THEN C C---- Update loop invariants C IX2 = MOD(IX+NBX,MBX) + 1 J32 = 0 C C DO 110 J = -NBY,NBY I32 = IFRAME(IX+NBX1,IY+J) IF (I32.LT.0) I32 = I32 + 65536 K32 = IFRAME(IX-NBX,IY+J) IF (K32.LT.0) K32 = K32 + 65536 J32 = J32 + I32 SUMROW(J+NBY1) = SUMROW(J+NBY1) + I32 - K32 110 CONTINUE C C SUM = (J32-SUMCOL(IX2)) + SUM SUMCOL(IX2) = J32 END IF 120 CONTINUE JRFL0 = JRFL1 130 CONTINUE END IF RETURN END C C C SUBROUTINE TRANSP(A,B) C ====================== C C C C---- Transpose a 3x3 matrix C C A = BT C C C C .. Array Arguments .. REAL A(3,3),B(3,3) C .. C .. Local Scalars .. INTEGER I,J C .. SAVE C C DO 20 I = 1,3 DO 10 J = 1,3 A(I,J) = B(J,I) 10 CONTINUE 20 CONTINUE C C END C C C SUBROUTINE U3BEST(W,X,Y,N,MODE,RMS,U,T,IER) C =========================================== C C C---- Calculates a best rotation & translation between two vector sets C such that u*x+t is the closest approximation to y. C C The calculated best superposition may not be unique as indicated C by a result value ier=-1. however it is guarantied that within C numerical tolerances no other superposition exists giving a C smaller value for rms. C C This version of the algorithm is optimized for three-dimensional C real vector space. C C Use of this routine is restricted to non-profit academic C APPLICATIONS. C PLEASE REPORT ERRORS TO C PROGRAMMER: W.KABSCH MAX-PLANCK-INSTITUTE FOR MEDICAL RESEARCH C JAHNSTRASSE 29, 6900 HEIDELBERG, FRG. C REFERENCES: W.KABSCH ACTA CRYST.(1978).A34,827-828 C W.KABSCH ACTA CRYST.(1976).A32,922-923 C C C C W - w(m) is weight for atom pair # m (given) C X - x(i,m) are coordinates of atom # m in set x (given) C Y - y(i,m) are coordinates of atom # m in set y (given) C N - n is number of atom pairs (given) C MODE - 0:calculate rms only (given) C 1:calculate rms,u,t (takes longer) C RMS - sum of w*(ux+t-y)**2 over all atom pairs (result) C U - u(i,j) is rotation matrix for best superposition (result) C T - t(i) is translation vector for best superposition (result) C IER - 0: a unique optimal superposition has been determined(result) C -1: superposition is not unique but optimal C -2: no result obtained because of negative weights w C or all weights equal to zero. C C .. Scalar Arguments .. REAL RMS INTEGER IER,MODE,N C .. C .. Array Arguments .. REAL T(3),U(3,3),W(1),X(3,1),Y(3,1) C .. C .. Local Scalars .. REAL*8 COF,CTH,D,DET,E0,E1,E2,E3,G,H,ONE,P,RR1,RR2,RR3,RR4,RR5, + RR6,SPUR,SQRT3,SQRTH,SS1,SS2,SS3,SS4,SS5,SS6,STH,THREE,TOL, + TWO,WC,ZERO REAL SIGMA INTEGER I,J,K,L,M,M1 C .. C .. Local Arrays .. REAL*8 A(3,3),B(3,3),E(3),R(3,3),RR(6),SS(6),XC(3),YC(3) INTEGER IP(9),IP2312(4) C .. C .. Intrinsic Functions .. INTRINSIC DABS,DATAN2,DCOS,DSIN,DSQRT C .. C .. Equivalences .. EQUIVALENCE (RR1,RR(1)), (RR2,RR(2)), (RR3,RR(3)), (RR4,RR(4)), + (RR5,RR(5)), (RR6,RR(6)), (SS1,SS(1)), (SS2,SS(2)), + (SS3,SS(3)), (SS4,SS(4)), (SS5,SS(5)), (SS6,SS(6)), + (E1,E(1)), (E2,E(2)), (E3,E(3)) C .. SAVE C .. Data statements .. DATA SQRT3,TOL/1.73205080756888D00,1.0D-2/ DATA ZERO,ONE,TWO,THREE/0.0D00,1.0D00,2.0D00,3.0D00/ DATA IP/1,2,4,2,3,5,4,5,6/ DATA IP2312/2,3,1,2/ C .. C C C WC = ZERO RMS = 0.0 E0 = ZERO C C DO 20 I = 1,3 E(I) = ONE XC(I) = ZERO YC(I) = ZERO T(I) = 0.0 C C DO 10 J = 1,3 D = ZERO IF (I.EQ.J) D = ONE U(I,J) = D A(I,J) = D R(I,J) = ZERO 10 CONTINUE 20 CONTINUE C C IER = -1 IF (N.GE.1) THEN C C---- Determine centroids of both vector sets x and y C IER = -2 C C DO 40 M = 1,N IF (W(M).LT.0.0) THEN RETURN ELSE WC = W(M) + WC DO 30 I = 1,3 XC(I) = W(M)*X(I,M) + XC(I) YC(I) = W(M)*Y(I,M) + YC(I) 30 CONTINUE END IF 40 CONTINUE C C IF (WC.GT.ZERO) THEN C C DO 50 I = 1,3 XC(I) = XC(I)/WC YC(I) = YC(I)/WC 50 CONTINUE C C---- Determine correlation matrix r between vector sets y and x C DO 80 M = 1,N DO 70 I = 1,3 E0 = ((X(I,M)-XC(I))**2+ (Y(I,M)-YC(I))**2)*W(M) + E0 D = (Y(I,M)-YC(I))*W(M) DO 60 J = 1,3 R(I,J) = (X(J,M)-XC(J))*D + R(I,J) 60 CONTINUE 70 CONTINUE 80 CONTINUE C C---- Calculate determinant of r(i,j) C SIGMA = (R(2,2)*R(3,3)-R(2,3)*R(3,2))*R(1,1) - + (R(2,1)*R(3,3)-R(2,3)*R(3,1))*R(1,2) + + (R(2,1)*R(3,2)-R(2,2)*R(3,1))*R(1,3) C C---- Form upper triangle of transposed(r)*r C M = 0 C C DO 100 J = 1,3 DO 90 I = 1,J M = M + 1 RR(M) = R(1,I)*R(1,J) + R(2,I)*R(2,J) + R(3,I)*R(3,J) 90 CONTINUE 100 CONTINUE C C---- Rescale matrix such that the sum of the new main diagonal C elements becomes 3. C SPUR = (RR1+RR3+RR6)/THREE IF (SPUR.GT.ZERO) THEN C C DO 110 M = 1,6 RR(M) = RR(M)/SPUR 110 CONTINUE C C---- Eigenvalues C C form characteristic cubic X**3-3*X**2+3*COF*X-DET=0 C COF = (RR3*RR6-RR5*RR5+RR1*RR6-RR4*RR4+RR1*RR3-RR2*RR2)/ + THREE DET = (RR3*RR6-RR5*RR5)*RR1 - (RR2*RR6-RR5*RR4)*RR2 + + (RR2*RR5-RR3*RR4)*RR4 C C---- Reduce cubic to standard form y**3-3hy+2g=0 by putting x=y+one C H = ONE - COF G = (COF-DET)/TWO - H C C---- Solve cubic. roots are e1,e2,e3 in decreasing order C IF (H.GT.ZERO) THEN SQRTH = DSQRT(H) D = G/ (H*SQRTH) P = (ONE-D)* (ONE+D) C C IF (P.LT.ZERO) THEN P = ZERO ELSE P = DSQRT(P) END IF C C D = DATAN2(P,-D)/THREE CTH = DCOS(D)*SQRTH STH = SQRTH*SQRT3*DSIN(D) E1 = ONE + CTH + CTH E2 = ONE - CTH + STH E3 = ONE - CTH - STH IF (MODE.NE.0) THEN C C---- Eigenvectors C C C DO 150 L = 1,3,2 D = E(L) SS1 = (D-RR3)* (D-RR6) - RR5*RR5 SS2 = (D-RR6)*RR2 + RR4*RR5 SS3 = (D-RR1)* (D-RR6) - RR4*RR4 SS4 = (D-RR3)*RR4 + RR2*RR5 SS5 = (D-RR1)*RR5 + RR2*RR4 SS6 = (D-RR1)* (D-RR3) - RR2*RR2 J = 1 IF (DABS(SS1).LT.DABS(SS3)) THEN J = 2 IF (DABS(SS3).GE.DABS(SS6)) GO TO 120 ELSE IF (DABS(SS1).GE.DABS(SS6)) THEN GO TO 120 END IF J = 3 120 D = ZERO J = (J-1)*3 DO 130 I = 1,3 K = IP(I+J) A(I,L) = SS(K) D = SS(K)*SS(K) + D 130 CONTINUE IF (D.GT.ZERO) D = ONE/DSQRT(D) C C DO 140 I = 1,3 A(I,L) = A(I,L)*D 140 CONTINUE 150 CONTINUE C C D = A(1,1)*A(1,3) + A(2,1)*A(2,3) + A(3,1)*A(3,3) M1 = 3 M = 1 IF ((E1-E2).LE. (E2-E3)) THEN M1 = 1 M = 3 END IF P = ZERO C C DO 160 I = 1,3 A(I,M1) = A(I,M1) - A(I,M)*D P = A(I,M1)**2 + P 160 CONTINUE C C IF (P.LE.TOL) THEN P = ONE C C DO 170 I = 1,3 IF (P.GE.DABS(A(I,M))) THEN P = DABS(A(I,M)) J = I END IF 170 CONTINUE C C K = IP2312(J) L = IP2312(J+1) P = DSQRT(A(K,M)**2+A(L,M)**2) IF (P.LE.TOL) THEN GO TO 270 ELSE A(J,M1) = ZERO A(K,M1) = -A(L,M)/P A(L,M1) = A(K,M)/P END IF ELSE P = ONE/DSQRT(P) C C DO 180 I = 1,3 A(I,M1) = A(I,M1)*P 180 CONTINUE C C END IF A(1,2) = A(2,3)*A(3,1) - A(2,1)*A(3,3) A(2,2) = A(3,3)*A(1,1) - A(3,1)*A(1,3) A(3,2) = A(1,3)*A(2,1) - A(1,1)*A(2,3) ELSE GO TO 290 END IF C C---- Handle special case of 3 identical roots C ELSE IF (MODE.EQ.0) THEN GO TO 290 END IF C C---- Rotation matrix C C C DO 210 L = 1,2 D = ZERO C C DO 190 I = 1,3 B(I,L) = R(I,1)*A(1,L) + R(I,2)*A(2,L) + R(I,3)*A(3,L) D = B(I,L)**2 + D 190 CONTINUE C C IF (D.GT.ZERO) D = ONE/DSQRT(D) DO 200 I = 1,3 B(I,L) = B(I,L)*D 200 CONTINUE 210 CONTINUE C C D = B(1,1)*B(1,2) + B(2,1)*B(2,2) + B(3,1)*B(3,2) P = ZERO C C DO 220 I = 1,3 B(I,2) = B(I,2) - B(I,1)*D P = B(I,2)**2 + P 220 CONTINUE C C IF (P.LE.TOL) THEN P = ONE C C DO 230 I = 1,3 IF (P.GE.DABS(B(I,1))) THEN P = DABS(B(I,1)) J = I END IF 230 CONTINUE C C K = IP2312(J) L = IP2312(J+1) P = DSQRT(B(K,1)**2+B(L,1)**2) IF (P.LE.TOL) THEN GO TO 270 ELSE B(J,2) = ZERO B(K,2) = -B(L,1)/P B(L,2) = B(K,1)/P END IF ELSE P = ONE/DSQRT(P) C C DO 240 I = 1,3 B(I,2) = B(I,2)*P 240 CONTINUE C C END IF B(1,3) = B(2,1)*B(3,2) - B(2,2)*B(3,1) B(2,3) = B(3,1)*B(1,2) - B(3,2)*B(1,1) B(3,3) = B(1,1)*B(2,2) - B(1,2)*B(2,1) C C DO 260 I = 1,3 DO 250 J = 1,3 U(I,J) = B(I,1)*A(J,1) + B(I,2)*A(J,2) + B(I,3)*A(J,3) 250 CONTINUE 260 CONTINUE END IF C C---- Translation vector C 270 CONTINUE C C DO 280 I = 1,3 T(I) = YC(I) - U(I,1)*XC(1) - U(I,2)*XC(2) - U(I,3)*XC(3) 280 CONTINUE C C---- rms error C 290 CONTINUE C C DO 300 I = 1,3 D = E(I)*SPUR IF (D.LT.ZERO) D = ZERO E(I) = DSQRT(D) 300 CONTINUE C C IER = 0 IF (E2.LE. (E1*1.0D-05)) IER = -1 D = E3 IF (SIGMA.LT.0.0) THEN D = -D IF ((E2-E3).LE. (E1*1.0D-05)) IER = -1 END IF D = D + E2 + E1 RMS = E0 - D - D IF (RMS.LT.0.0) RMS = 0.0 END IF END IF END C C C SUBROUTINE ULIST(NEQL,MLAUE,A0,U,RMS,MXBEST,NBEST,UB,ESD) C ========================================================= C C C C C---- Checks a given orientation matrix of the reciprocal cell for C symmetry equivalence with an already existing list of possible C orientation matrices. if the new orientation matrix is not C related to any other and of sufficient quality it is entered C into the list. the list is kept sorted in decreasing quality C of the orientation matrices. C C VERSION 5-1987 C C---- External subroutines required; C C MATMUL,INVERS C C----Describtion of calling parameters C C C C NEQL - number of laue operators (given) C MLAUE - laue operators (given) C A0 - setting matrix of reciprocal unit cell (given) C U - orientation matrix to be checked for symmetry (given) C independence and quality. if the matrix is new C and of sufficient quality it will be saved in the C list of possible orientations ub. C RMS - quality measure for u (given) C MXBEST - maximum number of best unit cell orientations (given) C NBEST - actual number of unit cell orientations in the list (updated) C UB - array(3,3,mxbest) describing orientations (updated) C ESD - array(mxbest) describing quality of orientations (updated) C C C---- Save assignment if among the "mxbest" choices C save orientation matrix if list of possible C solutions is empty C C .. Scalar Arguments .. REAL RMS INTEGER MXBEST,NBEST,NEQL C .. C .. Array Arguments .. REAL A0(3,3),ESD(MXBEST),U(3,3),UB(3,3,MXBEST) INTEGER MLAUE(216) C .. C .. Local Scalars .. REAL MINUS,PLUS,R INTEGER I,J,K,L,M C .. C .. Local Arrays .. REAL A1(3,3),B2(3,3),ML(9) C .. C .. External Subroutines .. EXTERNAL INVERS,MATMUL C .. C .. Intrinsic Functions .. INTRINSIC ABS,MAX,MIN C .. SAVE C C M = 1 IF (NBEST.LT.1) GO TO 160 C C---- check orientation matrix for being kept in the list C IF ((NBEST.EQ.MXBEST) .AND. (RMS.GE.ESD(NBEST))) RETURN C C---- Check if a symmetry equivalent orientation is already present C C *************** CALL MATMUL(U,A0,A1) CALL INVERS(A1,B2,R) C *************** C IF (R.LE.0.0) RETURN C C DO 80 I = 1,NBEST DO 20 J = 1,3 DO 10 K = 1,3 A1(J,K) = UB(J,1,I)*A0(1,K) + UB(J,2,I)*A0(2,K) + + UB(J,3,I)*A0(3,K) 10 CONTINUE 20 CONTINUE C C L = 0 C C DO 40 K = 1,3 DO 30 J = 1,3 L = L + 1 ML(L) = B2(J,1)*A1(1,K) + B2(J,2)*A1(2,K) + B2(J,3)*A1(3,K) 30 CONTINUE 40 CONTINUE C C L = 0 C C DO 70 J = 1,NEQL MINUS = 0.0 PLUS = 0.0 C C DO 50 K = 1,9 L = L + 1 R = ABS(ML(K)-MLAUE(L)) PLUS = MAX(PLUS,R) R = ABS(ML(K)+MLAUE(L)) MINUS = MAX(MINUS,R) 50 CONTINUE C C R = MIN(MINUS,PLUS) IF (R.GT.0.2) GO TO 60 C C---- A symmetry equivalent orientation matrix has been found C is the new orientation matrix better ? C C---- NO C IF (RMS.GE.ESD(I)) RETURN C C---- YES. replace old by new matrix C M = I GO TO 200 60 CONTINUE 70 CONTINUE 80 CONTINUE C C---- Determine lowest index of entries which have to be shifted C 90 M = NBEST + 1 100 IF (M.EQ.1) GO TO 110 C C IF (ESD(M-1).LT.RMS) GO TO 110 M = M - 1 GO TO 100 C C---- Shift entries of lower quality than the new one C 110 I = NBEST 120 IF (I.LT.M) GO TO 160 IF (I.GE.MXBEST) GO TO 150 ESD(I+1) = ESD(I) C C DO 140 J = 1,3 DO 130 K = 1,3 UB(J,K,I+1) = UB(J,K,I) 130 CONTINUE 140 CONTINUE C C 150 I = I - 1 GO TO 120 C C---- Insert new entry C 160 IF (M.GT.MXBEST) GO TO 190 ESD(M) = RMS C C DO 180 J = 1,3 DO 170 K = 1,3 UB(J,K,M) = U(J,K) 170 CONTINUE 180 CONTINUE C C IF (NBEST.LT.MXBEST) NBEST = NBEST + 1 190 CONTINUE C C---- return C RETURN C C 200 IF (M.EQ.NBEST) GO TO 240 C C DO 230 I = M + 1,NBEST ESD(I-1) = ESD(I) DO 220 J = 1,3 DO 210 K = 1,3 UB(J,K,I-1) = UB(J,K,I) 210 CONTINUE 220 CONTINUE 230 CONTINUE C C 240 NBEST = NBEST - 1 GO TO 90 END C C C SUBROUTINE UNORM(U,IER) C ======================= C C C---- Normalize matrix u(3,3) to orthogonal form C C C C C .. Scalar Arguments .. INTEGER IER C .. C .. Array Arguments .. REAL U(3,3) C .. C .. Local Scalars .. REAL P,Q INTEGER J C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. SAVE C C IER = -1 P = U(1,2)**2 + U(2,2)**2 + U(3,2)**2 Q = U(1,2)*U(1,3) + U(2,2)*U(2,3) + U(3,2)*U(3,3) C C DO 10 J = 1,3 U(J,3) = U(J,3)*P - U(J,2)*Q 10 CONTINUE C C Q = SQRT(U(1,3)**2+U(2,3)**2+U(3,3)**2) IF (Q.GE.0.000001) THEN P = SQRT(P) C C DO 20 J = 1,3 U(J,2) = U(J,2)/P U(J,3) = U(J,3)/Q 20 CONTINUE C C U(1,1) = U(2,2)*U(3,3) - U(2,3)*U(3,2) U(2,1) = U(3,2)*U(1,3) - U(3,3)*U(1,2) U(3,1) = U(1,2)*U(2,3) - U(1,3)*U(2,2) IER = 0 END IF RETURN END