PROGRAM LEVY C C ROBERT H. BLESSING C MEDICAL FOUNDATION OF BUFFALO C 73 HIGH STREET C BUFFALO, NEW YORK 14203, USA C TELEPHONE: (716) 856-9600, EXTENSION 335 C ELECTRONIC MAIL: Blessing@HWI.Buffalo.Edu C C FEBRUARY 1995,..., JANUARY 2001 C C A PRIORI SCALE FACTOR AND OVERALL ANISOTROPIC MEAN-SQUARE ATOMIC C DISPLACEMENT PARAMETERS BY AN EXTENSION OF THE METHOD OF C H.A. LEVY, W.E. THIESSEN, AND (IN PART) G.M. BROWN (1970). A LEAST- C SQUARES METHOD FOR CONVERTING OBSERVED INTENSITIES TO NORMALIZED C STRUCTURE FACTORS. AM. CRYST. ASSOC. MEETING, TULANE UNIV., NEW C ORLEANS, MARCH 1970. ABSRACT NO. B3. C C LEAST-SQUARES FIT OF WILSON EXPECTATION VALUES TO MEASURED INTENSITIES C C CHISQ = SUM W*DELTA**2 = SUM (DELTA/SIGMA)**2 C DELTA = YOBS - YCALC C YOBS = FSQ(MEAS)/[EPSILON*SUM F(J)**2] C YCAL = (SCALEK**-2)*EXP[-2*{H(J)*B(I,J)*H(I) - [H(J)*C(I,J)*H(I)]**2}] C W = M/[SIGMA(YOBS)**2 + YCALC**2*VAR] C M = MULTIPLICITY OF REFLECTION C VAR = 1 OR 2 FOR, RESPECTIVELY, ACENTRIC OR CENTRIC WILSON C DISTRIBUTIONS OF NORMALIZED INTENSITIES. C C YOBS IS A MEASURED SQUARED STRUCTURE FACTOR MAGNITUDE NORMALIZED TO C THE ABSOLUTE SCALE OF SCATTERING BY ATOMS AT REST. C C YCALC IS A WILSON EXPECTATION VALUE ON THE RELATIVE EXPERIMENTAL C SCALE FOR VIBRATING OR DISORDERED ATOMS. C C SCALEK AND B(I,J) ARE DEFINED ACCORDING TO: C C F(ABSOLUTE) = SCALEK*F(RELATIVE), C C f(T) = f(0)*EXP[-SUM(I) SUM(J) H(I)*H(J)*B(I,J)]. C C THE C(I,J) ARE DISPERSION PARAMETERS OF THE DISTRIBUTION OVER THE C UNIT CELL OF THE ANISOTROPIC ATOMIC DISPLACEMENT PARAMETERS B(I,J). C C IN THE ISOTROPIC APPROXIMATION, C C YCALC = (SCALEK**-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2}, C C WHERE S = SIN(THETA)/LAMBDA, AND BISO AND CISO ARE, RESPECTIVELY, THE C ESTIMATED MEAN AND STANDARD DEVIATION OF THE B-DISTRIBUTION. C CHARACTER ATIME*8,ADATE*9,TITLE*80,FILE*80,XRAY*2,SYST*12,LATT*1, & PTGP*5,LAUE*5,EL*2,ELANOM*2 INTEGER ZCELL DIMENSION EL(20),XI(20),ZI(20),RM(20),A1(20),B1(20),A2(20),B2(20), & A3(20),B3(20),A4(20),B4(20),C0(20),FP(20),FPP(20) DIMENSION CELL(6),GA(3,3),GB(3,3),B(3,3),C(3,3),U(3,3),V(3,3), & COV(13,13),BIBJ(6),VAR(13),TEMP(109) DATA B,C,COV,TEMP /9*0,9*0,169*0,109*0/ C C FORTRAN-90 TIME AND DATE TRANSLATIONS C CALL VTIME (ATIME) CALL VDATE (ADATE) c ATIME='hh:mm:ss' c ADATE='dd-mmm-yy' IO1=10 IO2=20 IO3=30 ILP=60 C C READ CONTROL DATA. C OPEN (UNIT=IO1,FILE='levy.dat',STATUS='OLD') READ (IO1,*) TITLE READ (IO1,*) ITYPE READ (IO1,*) FILE READ (IO1,*) LATT READ (IO1,*) PTGP READ (IO1,*) CELL CALL DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE) CALL METRIC (CELL,VCELL,GA,GB) ZCELL=0 XRAY=' ' NEL=0 READ (IO1,*,END=4) ZCELL READ (IO1,*,END=4) XRAY READ (IO1,*,END=4) NEL 4 CONTINUE IF (ZCELL.EQ.0) ZCELL=1 IF (XRAY.EQ.' ') XRAY='X ' IF (NEL.GT.0) THEN DO I=1,NEL READ (IO1,*) EL(I),XI(I) END DO ICELL=1 ELSE ICELL=0 C C DEFAULT EMPIRICAL FORMULA FOR AN AVERAGE PROTEIN, C4 H6 N O. C EL(1)='C ' EL(2)='H ' EL(3)='N ' EL(4)='O ' XI(1)=4 XI(2)=6 XI(3)=1 XI(4)=1 NEL=4 END IF F000=0 RMCELL=0 DO I=1,NEL CALL FTABLE (EL(I),ZI(I),RM(I),A1(I),B1(I),A2(I),B2(I),A3(I), & B3(I),A4(I),B4(I),C0(I),XRAY,FP(I),FPP(I)) XI(I)=ZCELL*XI(I) F000=F000+XI(I)*ZI(I) RMCELL=RMCELL+XI(I)*RM(I) END DO DX=(RMCELL/VCELL)/0.60221 C C READ ANOMALOUS DISPERSION CORRECTIONS FOR WAVELENGTHS OTHER THAN CU OR C MO K-ALPHA. C NANOM=0 READ (IO1,*,END=5) NANOM IF (NANOM.GT.0) THEN DO I=1,NEL FP(I)=0 FPP(I)=0 END DO DO I=1,NANOM READ (IO1,*) ELANOM,FPI,FPPI DO J=1,NEL IF (EL(J).EQ.ELANOM) THEN FP(J)=FPI FPP(J)=FPPI END IF END DO END DO END IF 5 CONTINUE SMIN=0 SMAX=0 ISOLV=0 ASOLV=0 BSOLV=0 READ (IO1,*,END=6) SMIN,SMAX READ (IO1,*,END=6) ISOLV,ASOLV,BSOLV 6 CONTINUE CLOSE (UNIT=IO1,STATUS='KEEP') IF (ICELL.EQ.0) THEN C C ALONG WITH THE DEFAULT EMPIRICAL FORMULA, C4 H6 N O, FOR AN AVERAGE C PROTEIN, ASSUME PROTEIN DENSITY 1.35 MG MM**-3 AND 50% SOLVENT VOLUME C FRACTION. C RMTEMP=1.35*0.5*VCELL*0.60221 T=RMTEMP/RMCELL RMCELL=RMTEMP F000=F000*NINT(T) DO I=1,4 XI(I)=XI(I)*NINT(T) END DO END IF C C CALCULATE MATTHEWS COEFFICIENT, SOLVENT VOLUME FRACTION, AND NUMBER OF C WATER MOLECULES AT 29.9 A**3 PER H2O MOLECULE TO FILL THE SOLVENT c VOLUME IN PROTEIN CRYSTALS. C VM=VCELL/RMCELL VS=1-1.23*RMCELL/VCELL NWATER=NINT(VS*VCELL/29.9) IF (VS.GT.0.1.AND.NWATER.GT.0) THEN C C CALCULATE BULK-SOLVENT-CORRECTED UNIT-CELL SCATTERING POWER AND MASS C DENSITY. C F000=F000+NWATER*10 DX=((RMCELL+NWATER*18.0153)/VCELL)/0.60221 END IF C C ECHO CONTROL DATA. C OPEN (UNIT=ILP,FILE='levy.lp',STATUS='NEW') WRITE (ILP,600) ATIME,ADATE,TITLE WRITE (ILP,610) ITYPE WRITE (ILP,601) FILE WRITE (ILP,602) SYST,LATT,PTGP,LAUE IF (PTGP.EQ.LAUE) THEN WRITE (ILP,604) ELSE WRITE (ILP,603) END IF WRITE (ILP,606) CELL,VCELL IF (ICELL.EQ.0) WRITE (ILP,670) WRITE (ILP,607) ZCELL,(EL(I),XI(I)/ZCELL,I=1,NEL) WRITE (ILP,608) NINT(F000),RMCELL/ZCELL,DX,VM,100*VS,NWATER/ZCELL WRITE (ILP,609) XRAY,(EL(I),FP(I),FPP(I),I=1,NEL) IF (SMIN.NE.0) WRITE (ILP,611) SMIN,1/(2*SMIN) IF (SMAX.NE.0) WRITE (ILP,612) SMAX,1/(2*SMAX) IF (ISOLV.NE.0) THEN IF (ASOLV.EQ.0) ASOLV=0.75 IF (BSOLV.EQ.0) BSOLV=250 WRITE (ILP,613) ASOLV,BSOLV END IF 600 FORMAT (/1H1,130('-')/1X,'PROGRAM LEVY. ',A,', ',A,'. ',A) 610 FORMAT (/1X,'INPUT REFLECTION FILE TYPE: ITYPE = ',I2) 601 FORMAT (/1X,'INPUT REFLECTION FILE NAME: ',A) 602 FORMAT (/1X,A//1X,A,'-LATTICE'//1X, &'CRYSTAL CLASS POINT GROUP ',A//1X, &'LAUE GROUP ',A) 603 FORMAT (/1X,'NONCENTROSYMMETRIC STRUCTURE') 604 FORMAT (/1X,'CENTROSYMMETRIC STRUCTURE') 606 FORMAT (/1X,9X,'A',9X,'B',9X,'C',5X,'ALPHA',6X,'BETA',5X,'GAMMA', &14X,'V'/1X,3F10.4,3F10.3,F15.2) 670 FORMAT (/1X, &'DEFAULT CRYSTAL CHEMICAL COMPOISITION FOR AN AVERAGE PROTEIN:'/ &1X, &'-------------------------------------------------------------'/ &1X, &' EMPIRICAL FORMULA C4 H6 N O'/1X, &' SPECIFIC PARTIAL VOLUME 0.74 MM**3 MG-1'/1X, &' MOLECULAR MASS DENSITY 1.35 MG MM**-3'/1X, &' SOLVENT VOLUME FRACTION 50% H2O') 607 FORMAT (/1X, &'UNIT CELL CONTENTS:'/1X, &'-------------------'//1X, &' Z = ',I2,' CRYSTAL CHEMICAL UNITS PER UNIT CELL'//1X, &' CRYSTAL CHEMICAL UNIT:'//(' ',A2,F10.2)) 608 FORMAT (/1X, &'F000 = ',I10,' E (ELECTRONS)'/1X, &'MOL. WT. = ',F10.2,' DA (DALTONS)'/1X, &'D(X-RAY) = ',E10.3,' MG MM**-3 (MILLIGRAMS PER CUBIC MILLIMETER)' &//1X, &'MATTHEWS COEFFICIENT AND SOLVENT VOLUME FRACTION:'/1X, &'-------------------------------------------------'/1X, &' VM = ',E10.3,' A**3 DA**-1 (CUBIC ANGSTROMS PER DALTON)'/ &1X, &' VS = ',E10.3,' %'/1X, &' N(H2O) = ',I10,' WATER MOLECULES PER CRYSTAL CHEMICAL UNIT NE', &'EDED TO FILL THE SOLVENT VOLUME') 609 FORMAT (/1X, &'RADIATION: ',A2//1X, &'ANOMALOUS DISPERSION CORRECTION TERMS:'/1X, &'--------------------------------------'/1X, &'EL FP FPP'/1X, &'-- -- ---'/ &(1X,A2,2F10.3)) 611 FORMAT (/1X, &'MINIMUM SIN(THETA)/LAMBDA FOR INCLUDING REFLECTION IN LEAST-SQU', &'ARES FITTING'/1X, &' SMIN = ',F7.4,' RECIPROCAL ANGSTROMS'/1X, &' DMAX = ',F5.2,' ANGSTROMS') 612 FORMAT (/1X, &'MAXIMUM SIN(THETA)/LAMBDA FOR INCLUDING REFLECTION IN LEAST-SQU', &'ARES FITTING'/1X, &' SMAX = ',F7.4,' RECIPROCAL ANGSTROMS'/1X, &' DMIN = ',F5.2,' ANGSTROMS') 613 FORMAT (/1X, &'BULK SOLVENT COMPENSATION PARAMETERS:'/1X, &'-------------------------------------'/1X, &' F(PROT) = F(MEAS)/(1 - ASOLV*EXP{-BSOLV*[SIN(THETA)/LAMBDA]**', &'2})'/1X, &' ASOLV = ',E10.3/1X, &' BSOLV = ',E10.3' A**2 (SQUARED ANGSTROMS)') C C BUILD WORKING REFLECTION DATA FILES. C CALL DATAIN (FILE,ITYPE,IO1,IO2,IO3,GB,SMIN,SMAX,SLOW,SHIGH, & ASOLV,BSOLV,MLATT,IPTGP,ILAUE,NEL,XI,A1,B1,A2,B2,A3,B3,A4,B4,C0, & FP,FPP,MDATA,NDATA) C C PRINT SOME DATA SET STATISTICS. C WRITE (ILP,600) ATIME,ADATE,TITLE WRITE (ILP,614) MDATA,NDATA,SMIN,SMAX,1/(2*SMIN),1/(2*SMAX), & SLOW,SHIGH,1/(2*SLOW),1/(2*SHIGH) 614 FORMAT (/1X, &'NUMBER OF REFLECTIONS IN THE INPUT DATA FILE:'/1X, &' MDATA = ',I6/1X, &'NUMBER OF REFLECTIONS USED IN THE LEAST-SQUARES FITTING:'/1X, &' NDATA = ',I6//1X, &'RESOLUTION LIMITS OF THE INPUT DATA SET:'/1X, &' S(MIN) = SIN(THETA(MIN))/LAMBDA = ',F6.3,' RECIPROCAL ANGSTRO', &'MS'/1X, &' S(MAX) = ',F6.3/1X, &' D(MIN) = 1/(2*S(MAX)) = ',F5.2,' ANGSTROMS'/1X, &' D(MAX) = 1/(2*S(MIN)) = ',F5.2//1X, &'RESOLUTION LIMITS OF THE LEAST-SQUARES DATA SET:'/1X, &' S(MIN) = SIN(THETA(MIN))/LAMBDA = ',F6.3,' RECIPROCAL ANGSTRO', &'MS'/1X, &' S(MAX) = ',F6.3/1X, &' D(MIN) = 1/(2*S(MAX)) = ',F5.2,' ANGSTROMS'/1X, &' D(MAX) = 1/(2*S(MIN)) = ',F5.2) C C GET STARTING SCALEK, BISO, AND CISO. C CALL WILSON (IO1,NDATA,SCALEK,BISO,CISO,COV) C C REFINE SCALEK, BISO, AND CISO. C CALL KBCISO (IO1,NDATA,SCALEK,BISO,CISO,NCYCLE,R,Z,NFREE,NPAR,COV) C C STORE REFINED SCALEK, BISO, AND CISO. C B(1,1)=BISO C(1,1)=CISO CALL STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,TEMP) C C IF CISO = 0, TRY REFINING AGAIN FROM PERTURBED BISO AND CISO. C IF (CISO.EQ.0) THEN CISO=BISO/2 BISO=2*BISO T=R NPAR=3 CALL KBCISO (IO1,NDATA,SCALEK,BISO,CISO,NCYCLE,R,Z,NFREE,NPAR, & COV) C C TEST FOR AN (R .LT. T) IMPROVED FIT. C IF ((R-T)/T.LT.-1E-6) THEN C C STORE RE-REFINED SCALEK AND BISO AND CISO. C B(1,1)=BISO C(1,1)=CISO CALL STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,TEMP) ELSE C C RECALL SCALEK AND BISO WITH CISO = 0. C CALL RECALL (TEMP,NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV) BISO=B(1,1) CISO=0 END IF END IF C C PRINT RESULTS OF ISOTROPIC FITTING. C WRITE (ILP,600) ATIME,ADATE,TITLE WRITE (ILP,616) NCYCLE,R,Z,NFREE,NPAR WRITE (ILP,617) SCALEK,BISO,CISO DO I=1,3 VAR(I)=COV(I,I) IF (VAR(I).EQ.0) VAR(I)=1 END DO WRITE (ILP,618) (SQRT(COV(I,I)),I=1,3),(I,(COV(I,J)/SQRT(VAR(I)* & VAR(J)),J=1,3),I=1,3) 616 FORMAT (/1X, &'SCALE AND ISOTROPIC ATOMIC DISPLACEMENT AND DISTRIBUTION DISPER', &'SION'/1X, &'PARAMETERS:'//1X, &' YOBS = FSQ(MEAS)/[EPSILON*SUM F(J)**2]'/1X, &' YCALC = (SCALEK**-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2)}'//1X, &' WHERE'/1X, &' EXP{-2*[BISO - (CISO*S)**2]*S**2} = '/1X, &' IS THE STATISTICAL EXPECTATION VALUE OF THE DEBYE-WALLER FACT', &'OR FOR'/1X, &' AN ASSUMED-TO-BE-NORMAL DISTRIBUTION OF ATOMIC B-VALUES WITH ', &'MEAN'/1X, &' BISO = '/1X, &' AND STANDARD DEVIATION'/1X, &' CISO = <(B - )**2>**(1/2).'//1X, &'STATISTICS OF FIT:'//1X, &'NCYCLE = ',I11/1X, &'R = ',E11.4/1X, &'Z = ',E11.4/1X, &'NFREE = ',I11/1X, &'NPAR = ',I11) 617 FORMAT (/1X, &'RESULTS:'//1X,' SCALEK BISO CISO'/1X,3E12.4) 618 FORMAT (/1X, &'PARAMETER E.S.D.''S'//1X, &' SIGMA(K) SIGMA(B) SIGMA(C)'/1X,3E12.2//1X, &'CORRELATION MATRIX'//1X, &' K BISO CISO'//1X, &' 1 2 3'//(1X,I2,3F6.2)) C C PLOT RESULTS OF ISOTROPIC FITTING. C CALL WPLOT (IO3,ILP,ATIME,ADATE,TITLE,SCALEK,BISO,CISO,MDATA) C C STORE SCALEK AND ISOTROPIC B AND C. C B(1,1)=BISO C(1,1)=CISO CALL STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,TEMP) C C EVALUATE EQUIVALENT ANISOTROPIC PARAMETERS. C CALL EXPAND (GB,BISO,B) CALL EXPAND (GB,CISO,C) CALL RESET (SYST,B) CALL RESET (SYST,C) C C CALCULATE U'S FROM B'S. C PI=3.14159 UISO=BISO/(8*PI**2) VISO=CISO/(8*PI**2) DO I=1,3 DO J=I,3 U(I,J)=B(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J))) U(J,I)=U(I,J) V(I,J)=C(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J))) V(J,I)=V(I,J) END DO END DO WRITE (ILP,600) ATIME,ADATE,TITLE WRITE (ILP,619) BISO,CISO,UISO,VISO, & B(1,1),B(2,2),B(3,3),B(1,2),B(1,3),B(2,3), & C(1,1),C(2,2),C(3,3),C(1,2),C(1,3),C(2,3), & U(1,1),U(2,2),U(3,3),U(1,2),U(1,3),U(2,3), & V(1,1),V(2,2),V(3,3),V(1,2),V(1,3),V(2,3) 619 FORMAT (/1X, &'BISO = ',E12.4,' CISO = ',E12.4/1X, &'UISO = ',E12.4,' VISO = ',E12.4//1X, &'EQUIVALENT ANISOTROPIC ATOMIC DISPLACEMENT PARAMETERS'/1X, &'AND DISTRIBUTION DISPERSION PARAMETERS:'//1X, &' B(1,1) B(2,2) B(3,3) B(1,2) B(1,3) ', &' B(2,3)'/1X,6E12.4/1X, &' C(1,1) C(2,2) C(3,3) C(1,2) C(1,3) ', &' C(2,3)'/1X,6E12.4//1X, &' U(1,1) U(2,2) U(3,3) U(1,2) U(1,3) ', &' U(2,3)'/1X,6E12.4/1X, &' V(1,1) V(2,2) V(3,3) V(1,2) V(1,3) ', &' V(2,3)'/1X,6E12.4) C C REFINE SCALEK, B(I,J), AND C(I,J). C T=R IF (CISO.EQ.0) THEN NPAR=7 ELSE NPAR=13 END IF CALL KBCIJ (IO1,NDATA,GA,GB,SCALEK,B,C,SYST,NCYCLE,R,Z,NFREE,NPAR, & COV) C C CALCULATE EQUIVALENT ISOTROPIC B AND C. C CALL CONTRACT (GA,B,BISO) CALL CONTRACT (GA,C,CISO) C C CALCULATE U'S FROM B'S. C UISO=BISO/(8*PI**2) VISO=CISO/(8*PI**2) DO I=1,3 DO J=I,3 U(I,J)=B(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J))) U(J,I)=U(I,J) V(I,J)=C(I,J)/(2*PI**2*SQRT(GB(I,I)*GB(J,J))) V(J,I)=V(I,J) END DO END DO WRITE (ILP,600) ATIME,ADATE,TITLE WRITE (ILP,624) NCYCLE,R,Z,NFREE,NPAR WRITE (ILP,625) SCALEK, & B(1,1),B(2,2),B(3,3),B(1,2),B(1,3),B(2,3), & C(1,1),C(2,2),C(3,3),C(1,2),C(1,3),C(2,3), & U(1,1),U(2,2),U(3,3),U(1,2),U(1,3),U(2,3), & V(1,1),V(2,2),V(3,3),V(1,2),V(1,3),V(2,3) WRITE (ILP,626) BISO,CISO,UISO,VISO WRITE (ILP,600) ATIME,ADATE,TITLE BIBJ(1)=GB(1,1) BIBJ(2)=SQRT(GB(1,1)*GB(2,2)) BIBJ(3)=SQRT(GB(1,1)*GB(3,3)) BIBJ(4)=GB(2,2) BIBJ(5)=SQRT(GB(2,2)*GB(3,3)) BIBJ(6)=GB(3,3) DO I=1,13 VAR(I)=COV(I,I) IF (VAR(I).EQ.0) VAR(I)=1 END DO WRITE (ILP,627) (SQRT(COV(1+I,1+I))/(2*PI**2*BIBJ(I)),I=1,6), & (SQRT(COV(7+I,7+I))/(2*PI**2*BIBJ(I)),I=1,6), & (SQRT(COV(I,I)),I=1,13), & (I,(NINT(100*COV(I,J)/SQRT(VAR(I)*VAR(J))),J=1,13),I=1,13) 624 FORMAT (/1X, &'SCALE AND ANISOTROPIC ATOMIC DISPLACEMENT AND DISTRIBUTION DISP', &'ERSION'/1X, &'PARAMETERS:'//1X, &' YOBS = FSQ(MEAS)/(EPSILON*SUM F(J)**2)'/1X, &' YCALC = (SCALEK**-2)*EXP(-2*H(J)*B(I,J)*H(I) + 2*(H(J)*C(I,J)', &'*H(I))**2)'//1X, &'STATISTICS OF FIT:'//1X, &'NCYCLE = ',I11/1X, &'R = ',E11.4/1X, &'Z = ',E11.4/1X, &'NFREE = ',I11/1X, &'NPAR = ',I11) 625 FORMAT (/1X, &'RESULTS:'//1X, &'SCALEK, B(I,J), C(I,J), U(I,J), AND B(I,J) ARE DEFINED ACCORDIN', &'G TO:'//1X, &' F(ABSOLUTE) = SCALEK*F(RELATIVE)'/1X, &' f(T) = f(0)*EXP(-SUM(I)SUM(J) H(I)*H(J)*B(I,J))'/1X, &' f(T) = f(0)*EXP(-2*PI**2*SUM(I)SUM(J) H(I)*H(J)*ASTAR(I)*ASTA', &'R(J)*U(I,J))'/1X, &' = EXP(-2*H(J)**H(I) + 2*(H', &'(J)**H(I))**2)'//1X, &'SCALEK'/1X,E12.4//1X, &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENT PARAMETERS:'/1X, &' B(1,1) B(2,2) B(3,3) B(1,2) B(1,3) ', &' B(2,3)'/1X,6E12.4,' (DIMENSIONLESS)'//1X, &'B-DISTRIBUTION DISPERSION PARAMETERS:'/1X, &' C(1,1) C(2,2) C(3,3) C(1,2) C(1,3) ', &' C(2,3)'/1X,6E12.4,' (DIMENSIONLESS)'//1X, &'OVERALL MEAN-SQUARE ATOMIC DISPLACEMENTS:'/1X, &' U(1,1) U(2,2) U(3,3) U(1,2) U(1,3) ', &' U(2,3)'/1X,6E12.4,' SQUARED ANGSTROMS'//1X, &'U-DISTRIBUTION DISPERSION PARAMETERS:'/1X, &' V(1,1) V(2,2) V(3,3) V(1,2) V(1,3) ', &' V(2,3)'/1X,6E12.4,' SQUARED ANGSTROMS') 626 FORMAT (/1X, &'EQUIVALENT ISOTROPIC VALUES (ALL IN SQUARED ANGSTROM):'//1X, &'BISO = ',E12.4,' = '/1X, &'CISO = ',E12.4,' = <(B - )**2>**(1/2)'//1X, &'UISO = ',E12.4,' = '/1X, &'VISO = ',E12.4,' = <(u**2 - )**2>**(1/2)'//1X, &'WHERE BISO = 8*PI**2*UISO') 627 FORMAT (/1X, &'PARAMETER E.S.D.''S'//1X, &'NOTE THAT THE ORDER OF THE MEAN-SQUARE DISPLACEMENT PARAMETER E'', &RRORS IS'/1X, &'11, 12, 13, 22, 23, 33 NOT 11, 22, 33, 12, 13, 23.'/1X, &' ---'//!X, &' SIGMA(U11) SIGMA(U12) SIGMA(U13) SIGMA(U22) SIGMA(U23) S', &'IGMA(U33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X, &' SIGMA(V11) SIGMA(V12) SIGMA(V13) SIGMA(V22) SIGMA(V23) S', &'IGMA(V33)'/1X,6E12.2,' SQUARED ANGSTROMS'//1X, &' SIGMA(K)'/1X,E12.2//1X, &' SIGMA(B11) SIGMA(B12) SIGMA(B13) SIGMA(B22) SIGMA(B13) S', &'IGMA(B33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X, &' SIGMA(C11) SIGMA(C12) SIGMA(C13) SIGMA(C12) SIGMA(C23) S', &'IGMA(C33)'/1X,6E12.2,' (DIMENSIONLESS)'//1X, &'CORRELATION MATRIX 100*CORR(I,J):'//1X, &' K B11 B12 B13 B22 B23 B33 C11 C12 C13 C22 C23 C33'//1X, &' 1 2 3 4 5 6 7 8 9 10 11 12 13'/ &/(1X,I2,2X,13I4)) C C IF THE ANISOTROPIC MODEL DID NOT IMPROVE THE FIT, RECALL THE ISOTROPIC C MODEL C IF ((R-T)/T.LT.-1E-6) THEN CONTINUE ELSE CALL RECALL (TEMP,NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV) END IF C C WRITE OUTPUT FILE. C CLOSE (UNIT=IO1,STATUS='DELETE') OPEN (UNIT=IO1,FILE='eval.dat',STATUS='NEW') WRITE (IO1,700) ATIME WRITE (IO1,700) ADATE WRITE (IO1,700) TITLE WRITE (IO1,702) ITYPE WRITE (IO1,700) FILE WRITE (IO1,700) LATT WRITE (IO1,700) PTGP WRITE (IO1,701) CELL WRITE (IO1,700) XRAY WRITE (IO1,702) NEL DO I=1,NEL WRITE (IO1,703) EL(I),XI(I),ZI(I),A1(I),B1(I),A2(I),B2(I) WRITE (IO1,704) A3(I),B3(I),A4(I),B4(I),C0(I),FP(I),FPP(I) END DO WRITE (IO1,701) SMIN,SMAX WRITE (IO1,706) ASOLV,BSOLV WRITE (IO1,705) SCALEK WRITE (IO1,705) B(1,1),B(1,2),B(1,3),B(2,2),B(2,3),B(3,3) WRITE (IO1,705) C(1,1),C(1,2),C(1,3),C(2,2),C(2,3),C(3,3) WRITE (IO1,706) ((COV(I,J),J=1,13),I=1,13) 700 FORMAT (1X,A) 701 FORMAT (1X,3F10.4,3F10.3) 702 FORMAT (1X,I6) 703 FORMAT (1X,A2,F10.2,F6.0,4(1X,F10.6)) 704 FORMAT (1X,5(1X,F10.6),2X,2(1X,F10.6)) 705 FORMAT (1X,6E12.5) 706 FORMAT (1X,7E10.3) CLOSE (UNIT=IO1,STATUS='KEEP') STOP 'PROGRAM LEVY FINIS!' END C----------------------------------------------------------------------- SUBROUTINE AGREE (IFILE,NDATA,SCALEK,B,C,NPAR,NFREE,Z,R) C C AGREEMENT STATISTICS C DIMENSION H(3),B(3,3),C(3,3) IF (NPAR.LT.7) THEN BISO=B(1,1) CISO=C(1,1) END IF NOBS=0 CHISQ=0 SUMWY2=0 REWIND IFILE DO I=1,NDATA READ (IFILE) J,K,L,SSQ,Y,SIGY,MULT,VAR IF (NPAR.LT.7) THEN HTBH=BISO*SSQ HTCH=CISO*SSQ ELSE H(1)=J H(2)=K H(3)=L HTBH=VTMV(H,B) HTCH=VTMV(H,C) END IF YCALC=(1/SCALEK**2)*EXP(-2*(HTBH-HTCH**2)) W=MULT/(SIGY**2+YCALC**2*VAR) CHISQ=CHISQ+W*(Y-YCALC)**2 SUMWY2=SUMWY2+W*Y**2 NOBS=NOBS+MULT END DO NFREE=NOBS-NPAR Z=SQRT(CHISQ/NFREE) R=SQRT(CHISQ/SUMWY2) RETURN END C----------------------------------------------------------------------- SUBROUTINE CONTRACT (G,B,BISO) C C CONTRACT TENSOR BIJ'S TO EQUIVALENT SCALAR BISO. C DIMENSION G(3,3),B(3,3) BISO=0 DO I=1,3 DO J=1,3 BISO=BISO+G(I,J)*B(I,J) END DO END DO BISO=BISO*4/3 RETURN END C----------------------------------------------------------------------- SUBROUTINE CSYM (PTGP,H,K,L,C) C C TABLE OF CENTROSYMMETRIC ROWS AND ZONES IN NON-CENTROSYMMETRIC C RECIPROCAL LATTICE POINT GROUPS C C TRICL. MONOCL. ORTHO. TETRAG. TRIG. RHOMB. CUB. C C (1) 1 (3) 2 (6) 222 (9) 4 (17) 3 (33) R3 (38) 23 C (2) -1 (4) M (7) MM2 (10) -4 (18) -3 (34) R-3 (39) M3 C (5) 2/M (8) MMM (11) 4/M C (19) 312 (35) R32 (40) 432 C (12) 422 (20) 321 (36) R3M (41) -43M C (13) 4MM (21) 31M (37) R-3M (42) M3M C (14) -42M (22) 3M1 C (15) -4M2 (23) -31M C (16) 4/MMM (24) -3M1 C C HEXAG. C C (25) 6 C (26) -6 C (27) 6/M C C (28) 622 C (29) 6MM C (30) -6M2 C (31) -62M C (32) 6/MMM C C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS: C C -42M -4M2 C 3 R3 C -3 R-3 C 312 321 R32 C 31M 3M1 R3M C -31M -3M1 R-3M C -62M -6M2 C C C = 0 FOR ACENTRIC DISTRIBUTIONS C C = 1 FOR CENTRIC DISTRIBUTIONS C INTEGER PTGP,H,C C=0 GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22, & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) PTGP 1 RETURN 2 C=1 RETURN 3 IF (K.EQ.0) C=1 RETURN 4 IF (H.EQ.0.AND.L.EQ.0) C=1 RETURN 5 C=1 RETURN 6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 RETURN 7 IF (L.EQ.0) C=1 RETURN 8 C=1 RETURN 9 IF (L.EQ.0) C=1 RETURN 10 IF (L.EQ.0) C=1 IF (H.EQ.0.AND.K.EQ.0) C=1 RETURN 11 C=1 RETURN 12 IF (L.EQ.0.OR.K.EQ.0.OR.H.EQ.0) C=1 IF (ABS(K).EQ.ABS(H)) C=1 RETURN 13 IF (L.EQ.0) C=1 RETURN 14 IF (L.EQ.0.OR.K.EQ.0.OR.H.EQ.0) C=1 RETURN 15 IF (L.EQ.0) C=1 IF (ABS(K).EQ.ABS(H)) C=1 RETURN 16 C=1 RETURN 17 RETURN 18 C=1 RETURN 19 IF (K.EQ.H.OR.K.EQ.-2*H) C=1 RETURN 20 IF (K.EQ.0.OR.H.EQ.0) C=1 RETURN 21 IF ((K.EQ.0.AND.L.EQ.0).OR.(H.EQ.0.AND.L.EQ.0)) C=1 RETURN 22 RETURN 23 C=1 RETURN 24 C=1 RETURN 25 IF (L.EQ.0) C=1 RETURN 26 IF (H.EQ.0.AND.K.EQ.0) C=1 RETURN 27 C=1 RETURN 28 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 IF (K.EQ.H.OR.K.EQ.-2*H) C=1 RETURN 29 IF (L.EQ.0) C=1 RETURN 30 IF (K.EQ.H.OR.K.EQ.-2*H) C=1 RETURN 31 IF (H.EQ.0.OR.K.EQ.0) C=1 RETURN 32 C=1 RETURN 33 RETURN 34 C=1 RETURN 35 IF (H.EQ.K.OR.K.EQ.L.OR.L.EQ.H) C=1 RETURN 36 IF (L.EQ.0.AND.K.EQ.-H) C=1 IF (K.EQ.0.AND.L.EQ.-H) C=1 IF (H.EQ.0.AND.L.EQ.-K) C=1 RETURN 37 C=1 RETURN 38 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 RETURN 39 C=1 RETURN 40 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 IF (ABS(H).EQ.ABS(K).OR.ABS(K).EQ.ABS(L).OR.ABS(L).EQ.ABS(H)) C=1 RETURN 41 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) C=1 RETURN 42 C=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE DATAIN (FILE,ITYPE,IO1,IO2,IO3,GINV,S1,S2,SMIN,SMAX, & ASOLV,BSOLV,MLATT,IPTGP,ILAUE,NEL,X,A1,B1,A2,B2,A3,B3,A4,B4,C0, & FP,FPP,MDATA,NDATA) C C CREATE WORKING DATA FILES. C CHARACTER FILE*80,FORM*11 DIMENSION GINV(3,3),X(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1), & A4(1),B4(1),C0(1),FP(1),FPP(1) C C FORTRAN-90 DYNAMIC STORAGE ALLOCATION C REAL, ALLOCATABLE::DATA(:) C C COPY DATA FROM THE INPUT REFLECTION DATA FILE TO A TEMPORARY FILE: C IF (ITYPE.LT.3) THEN FORM='FORMATTED' ELSE FORM='UNFORMATTED' END IF OPEN (UNIT=IO1,FILE=FILE,STATUS='OLD',FORM=FORM) OPEN (UNIT=IO2,STATUS='SCRATCH',FORM='UNFORMATTED') IF (S1.EQ.0) S1=-9 IF (S2.EQ.0) S2=+9 SMIN=+9 SMAX=-9 IF (IPTGP.LT.ILAUE) THEN JMIN=+999 JMAX=-999 KMIN=+999 KMAX=-999 LMIN=+999 LMAX=-999 END IF MDATA=0 NDATA=0 IEND=0 1 CALL READ1 (ITYPE,IO1,IEND,IH,IK,IL,FSQ,SIGFSQ) IF (IEND.NE.0) GO TO 9 S=SINTHL(IH,IK,IL,GINV) SMIN=MIN(SMIN,S) SMAX=MAX(SMAX,S) IF (ASOLV.NE.0) THEN C C COMPUTE BULK SOLVENT COMPENSATION FACTOR. C Q=1/(1-ASOLV*EXP(-BSOLV*S**2)) ELSE Q=1 END IF WRITE (IO2) IH,IK,IL,Q**2*FSQ,Q**2*SIGFSQ MDATA=MDATA+1 IF (IPTGP.LT.ILAUE) THEN C C STORE THE LAUE-UNIQUE INDEX LIMITS. C J=IH K=IK L=IL CALL EQUIV (ILAUE,J,K,L) JMIN=MIN(JMIN,J) JMAX=MAX(JMAX,J) KMIN=MIN(KMIN,K) KMAX=MAX(KMAX,K) LMIN=MIN(LMIN,L) LMAX=MAX(LMAX,L) END IF GO TO 1 9 CLOSE (UNIT=IO1,STATUS='KEEP') ENDFILE IO2 IF (IPTGP.LT.ILAUE) THEN C C STORE A TEMPORARY ARRAY OF THE BIJVOET- OR FRIEDEL-PAIR POPULATIONS C (0, 1, OR 2) INDEXED BY THE PACKED LAUE-UNIQUE INDICES. C NJ=JMAX-JMIN+1 NK=KMAX-KMIN+1 NL=LMAX-LMIN+1 NMAX=NJ*NK*NL+NK*NL+NL C C ALLOCATE STORAGE FOR DATA ARRAY. C ALLOCATE (DATA(NMAX)) DO I=1,NMAX DATA(I)=0 END DO REWIND IO2 DO I=1,MDATA READ (IO2) J,K,L CALL EQUIV (ILAUE,J,K,L) JKL=(J-JMIN)*NK*NL+(K-KMIN)*NL+(L-LMIN) DATA(JKL)=DATA(JKL)+1 END DO END IF C C EVALUATE ADDITIONAL REFLECTION INFORMATION AND STORE IT IN THE WORKING C DATA FILES. C REWIND IO2 OPEN (UNIT=IO1,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN (UNIT=IO3,STATUS='SCRATCH',FORM='UNFORMATTED', & ACCESS='DIRECT',RECL=12) DO I=1,MDATA READ (IO2) IH,IK,IL,FSQ,SIGFSQ C C GET ACENTRIC OR CENTRIC REFLECTION FLAG. C CALL CSYM (IPTGP,IH,IK,IL,IC) C C GET POINT GROUP DEGENERACY AND MULTIPLICITY. C CALL ESYM (IPTGP,IH,IK,IL,IE,M) C C ADJUST DEGENERACY FOR LATTICE MULTIPLICITY. C IE=MLATT*IE C C DOUBLE THE POINT GROUP MULTIPLICITY FOR ACENTRIC REFLECTIONS THAT ARE C REPRESENTED BY ONLY ONE MEMBER OF THE BIJVOET OR FRIEDEL PAIR. C IF (IC.EQ.0) THEN J=IH K=IK L=IL CALL EQUIV (ILAUE,J,K,L) JKL=(J-JMIN)*NK*NL+(K-KMIN)*NL+(L-LMIN) IF (DATA(JKL).EQ.1) M=2*M END IF C C GET SIN(THETA)/LAMBDA. C S=SINTHL(IH,IK,IL,GINV) C C CALCULATE Y = FSQ(MEAS)/(EPSILON*SUM F(J)**2) AND SIGMA(Y). C CALL FCALC (NEL,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S,SUMFSQ) XI=S**2 YI=FSQ/(IE*SUMFSQ) SIGYI=SIGFSQ/(IE*SUMFSQ) IF (YI.LT.-3*SIGYI) YI=-3*SIGYI C C WILSON DISTRIBUTION VARIANCE = 1 OR 2 FOR ACENTRIC OR CENTRIC C DISTRIBUTIONS, RESPECTIVELY. C IF (IC.EQ.0) V=1 IF (IC.EQ.1) V=2 C C WRITE A DATA FILE FOR SORTING AND PLOTTING. C WRITE (IO3,REC=I) XI,YI,M C C WRITE THE LEAST-SQUARES DATA FILE. C IF (S1.LE.S.AND.S.LE.S2) THEN WRITE (IO1) IH,IK,IL,XI,YI,SIGYI,M,V NDATA=NDATA+1 END IF END DO S1=MAX(S1,SMIN) S2=MIN(S2,SMAX) CLOSE (UNIT=IO2,STATUS='DELETE') DEALLOCATE (DATA) RETURN END C----------------------------------------------------------------------- SUBROUTINE DECODE (PTGP,LAUE,SYST,LATT,MLATT,IPTGP,ILAUE) C C CRYSTAL SYSTEM, LATTICE TYPE AND MULTIPLICITY, AND LAUE GROUP AND C POINT GROUP SYMMETRIES C C THERE ARE TEN CASES OF ALTERNATIVE AXES FOR SEVEN OF THE 32 C CRYSTALLOGRAPHIC POINT GROUPS (SEE SUBROUTINE ESYM), THUS A TOTAL OF C 42 CASES IS TABULATED. C C SIMILARLY, THERE ARE THREE CASES OF ALTERNATIVE AXES FOR TWO OF THE 11 C LAUE GROUPS, SO 14 CASES ARE TABULATED. C CHARACTER SYST*12,LATT*1,PTGP*5,LAUE*5,SYMBOL(42)*5 DATA SYMBOL &/'1 ','-1 ', & '2 ','M ','2/M ', & '222 ','MM2 ','MMM ', & '4 ','-4 ','4/M ', & '422 ','4MM ','-42M ','-4M2 ','4/MMM', & '3 ','-3 ', & '312 ','321 ','31M ','3M1 ','-31M ','-3M1 ', & '6 ','-6 ','6/M ', & '622 ','6MM ','-6M2 ','-62M ','6/MMM', & 'R3 ','R-3 ', & 'R32 ','R3M ','R-3M ', & '23 ','M3 ', & '432 ','-43M ','M3M '/ C C POINT GROUP INDEX NUMBERS: C C TRICLINIC 1, 2; C MONOCLINIC 3, 4, 5; C ORTHORHOMBIC 6, 7, 8; C TETRAGONAL 9, 10, 11, C 12, 13, 14, 15, 16; C TRIGONAL 17, 18, C 19, 20, 21, 22, 23, 24; C HEXAGONAL 25, 26, 27, C 28, 29, 30, 31, 32; C RHOMBOHEDRAL 33, 34, C 35, 36, 37; C CUBIC 38, 39, C 40, 41, 42. C DO I=1,42 IF (PTGP.EQ.SYMBOL(I)) GO TO 1 END DO STOP 'LEVY/DECODE: ERROR IN POINT GROUP SYMBOL' 1 IPTGP=I C C LAUE GROUP C IF (I.GE. 1) LAUE='-1 ' IF (I.GE. 3) LAUE='2/M ' IF (I.GE. 6) LAUE='MMM ' IF (I.GE. 9) LAUE='4/M ' IF (I.GE.12) LAUE='4/MMM' IF (I.GE.17) LAUE='-3 ' IF (I.GE.19) LAUE='-31M ' IF (I.EQ.20.OR.I.EQ.22.OR.I.EQ.24) LAUE='-3M1 ' IF (I.GE.25) LAUE='6/M ' IF (I.GE.28) LAUE='6/MMM' IF (I.GE.33) LAUE='R-3 ' IF (I.GE.35) LAUE='R-3M ' IF (I.GE.38) LAUE='M3 ' IF (I.GE.40) LAUE='M3M ' DO I=1,42 IF (LAUE.EQ.SYMBOL(I)) GO TO 2 END DO 2 ILAUE=I C C CRYSTAL SYSTEM C IF (I.GE. 1) SYST='TRICLINIC' IF (I.GE. 3) SYST='MONOCLINIC' IF (I.GE. 6) SYST='ORTHORHOMBIC' IF (I.GE. 9) SYST='TETRAGONAL' IF (I.GE.17) SYST='TRIGONAL' IF (I.GE.25) SYST='HEXAGONAL' IF (I.GE.33) SYST='RHOMBOHEDRAL' IF (I.GE.38) SYST='CUBIC' C C LATTICE CENTERING MULTIPLICITY C IF (LATT.EQ.'P') MLATT=1 IF (LATT.EQ.'A') MLATT=2 IF (LATT.EQ.'B') MLATT=2 IF (LATT.EQ.'C') MLATT=2 IF (LATT.EQ.'F') MLATT=4 IF (LATT.EQ.'I') MLATT=2 IF (LATT.EQ.'R') MLATT=1 IF (LATT.EQ.'R'.AND.SYST.NE.'RHOMBOHEDRAL') MLATT=3 RETURN END C----------------------------------------------------------------------- FUNCTION DELTA(I,J) C C KRONECKER DELTA C IF (I.EQ.J) THEN DELTA=1 ELSE DELTA=0 END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE EQUIV (PTGP,H,K,L) C C EQUIVALENT, UNIQUE REFLECTION INDICES FOR THE CRYSTALLOGRAPHIC POINT C GROUPS C C TRICL. MONOCL. ORTHO. TETRAG. TRIG. RHOMB. CUB. C C (1) 1 (3) 2 (6) 222 (9) 4 (17) 3 (33) R3 (38) 23 C (2) -1 (4) M (7) MM2 (10) -4 (18) -3 (34) R-3 (39) M3 C (5) 2/M (8) MMM (11) 4/M C (19) 312 (35) R32 (40) 432 C (12) 422 (20) 321 (36) R3M (41) -43M C (13) 4MM (21) 31M (37) R-3M (42) M3M C (14) -42M (22) 3M1 C (15) -4M2 (23) -31M C (16) 4/MMM (24) -3M1 C C HEXAG. C C (25) 6 C (26) -6 C (27) 6/M C C (28) 622 C (29) 6MM C (30) -6M2 C (31) -62M C (32) 6/MMM C C SEE INTERNATIONAL TABLES FOR CRYSTALLOGRAPHY (1983). VOLUME A, PP. C 750-770. C C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS: C C -42M -4M2 C 3 R3 C -3 R-3 C 312 321 R32 C 31M 3M1 R3M C -31M -3M1 R-3M C -62M -6M2 C INTEGER PTGP,H,X,Y,Z GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22, & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) PTGP C C 1 C 1 RETURN C C -1 C 2 IF (L.GT.0) RETURN H=-H K=-K L=-L IF (L.GT.0) RETURN IF (K.GT.0) RETURN H=-H K=-K IF (K.GT.0) RETURN IF (H.GE.0) RETURN H=-H RETURN C C 2 C 3 IF (L.GT.0) RETURN H=-H L=-L IF (L.GT.0) RETURN IF (H.GE.0) RETURN H=-H RETURN C C M C 4 K=ABS(K) RETURN C C 2/M C 5 K=ABS(K) GO TO 3 C C 222 C 6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) THEN H=ABS(H) K=ABS(K) L=ABS(L) RETURN END IF IF (H.GE.0.AND.K.GE.0) RETURN X=H Y=K Z=L H=-X K=-Y L=Z IF (H.GE.0.AND.K.GE.0) RETURN H=-X K=Y L=-Z IF (H.GE.0.AND.K.GE.0) RETURN H=X K=-Y L=-Z RETURN C C MM2 C 7 H=ABS(H) K=ABS(K) RETURN C C MMM C 8 L=ABS(L) GO TO 7 C C 4 C 9 IF (H.EQ.0.AND.K.EQ.0) RETURN IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K H=-Y K=X IF (H.GT.0.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=-X RETURN C C -4 C 10 IF (H.EQ.0.AND.K.EQ.0) THEN L=ABS(L) RETURN END IF IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K H=-X K=-Y IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=-X L=-L IF (H.GT.0.AND.K.GE.0) RETURN H=-Y K=X RETURN C C 4/M C 11 L=ABS(L) GO TO 9 C C 422 C 12 IF (H.EQ.0.OR.K.EQ.0.OR.ABS(H).EQ.ABS(K)) L=ABS(L) 512 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=X IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=-X IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=Y L=-L GO TO 512 C C 4MM C 13 H=ABS(H) K=ABS(K) IF (H.GE.K) RETURN X=H H=K K=X RETURN C C -42M C 14 IF (H.EQ.0.OR.K.EQ.0) L=ABS(L) 514 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=-X L=-L IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=X IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=Y GO TO 514 C C -4M2 C 15 H=ABS(H) K=ABS(K) IF (H.EQ.K) L=ABS(L) IF (H.GE.K) RETURN X=H H=K K=X L=-L RETURN C C 4/MMM C 16 L=ABS(L) GO TO 13 C C 3 C 17 IF (H.EQ.0.AND.K.EQ.0) RETURN IF (H.LT.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.LT.0.AND.K.GE.0) RETURN H=Y K=I RETURN C C -3 C 18 IF (H.EQ.0.AND.K.EQ.0) THEN L=ABS(L) RETURN END IF 518 IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN X=H Y=K I=-H-K H=I K=X IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN H=Y K=I IF (H.LT.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.LT.-K))) RETURN H=-X K=-Y L=-L GO TO 518 C C 312 C 19 IF (H.EQ.0.OR.K.EQ.0.OR.-H.EQ.K) L=ABS(L) 519 IF (H.GE.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.0.AND.K.GE.0) RETURN H=-Y K=-X L=-L GO TO 519 C C 321 C 20 IF (H.EQ.K.OR.-2*H.EQ.K.OR.-H.EQ.2*K) L=ABS(L) 520 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=X L=-L GO TO 520 C C 31M C 21 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=X GO TO 21 C C 3M1 C 22 IF (H.GE.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.0.AND.K.GE.0) RETURN H=-Y K=-X GO TO 22 C C -31M C 23 IF (H.EQ.0.OR.K.EQ.0.OR.-H.EQ.K) L=ABS(L) 523 IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=Y K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=Y K=X IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=X K=I IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=I K=Y IF (H.GE.0.AND.-H/2.LE.K.AND.K.LE.H.AND.(L.GT.0.OR.(L.EQ.0.AND. & K.GE.0))) RETURN H=-X K=-Y L=-L GO TO 523 C C -3M1 C 24 IF (H.EQ.K.OR.-2*H.EQ.K.OR.-H.EQ.2*K) L=ABS(L) 524 IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=Y K=I IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=Y K=X L=-L IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=X K=I IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=I K=Y IF (H.GE.0.AND.K.GE.0.AND.(L.GT.0.OR.(L.EQ.0.AND.H.GE.K))) RETURN H=-X K=-Y GO TO 524 C C 6 C 25 IF (H.EQ.0.AND.K.EQ.0) RETURN IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GT.0.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GT.0.AND.K.GE.0) RETURN H=-I K=-X IF (H.GT.0.AND.K.GE.0) RETURN H=-Y K=-I RETURN C C -6 C 26 L=ABS(L) GO TO 17 C C 6/M C 27 L=ABS(L) GO TO 25 C C 622 C 28 IF (H.EQ.0.OR.K.EQ.0.OR.ABS(H).EQ.ABS(K).OR.ABS(2*H).EQ.ABS(K).OR. & ABS(H).EQ.ABS(2*K)) L=ABS(L) 528 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=-I K=-X IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=-I IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=X L=-L GO TO 528 C C 6MM C 29 IF (H.GE.K.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=I IF (H.GE.K.AND.K.GE.0) RETURN H=-X K=-Y IF (H.GE.K.AND.K.GE.0) RETURN H=-I K=-X IF (H.GE.K.AND.K.GE.0) RETURN H=-Y K=-I IF (H.GE.K.AND.K.GE.0) RETURN H=Y K=X GO TO 29 C C -6M2 C 30 L=ABS(L) IF (H.EQ.0.AND.K.EQ.0) RETURN 530 IF (H.GT.0.AND.K.GE.0) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GT.0.AND.K.GE.0) RETURN H=Y K=I IF (H.GT.0.AND.K.GE.0) RETURN H=-Y K=-X GO TO 530 C C -62M C 31 L=ABS(L) IF (H.EQ.0.AND.K.EQ.0) RETURN 531 IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN X=H Y=K I=-H-K H=I K=X IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=I IF (H.GT.0.AND.-H/2.LE.K.AND.K.LE.H) RETURN H=Y K=X GO TO 531 C C 6/MMM C 32 L=ABS(L) GO TO 29 C C R3 C 33 STOP 'LEVY/EQUIV: RE-INDEX ON HEXAGONAL AXES' C C R-3 C 34 GO TO 33 C C R32 C 35 GO TO 33 C C R3M C 36 GO TO 33 C C R-3M C 37 GO TO 33 C C 23 C 38 IF (H.EQ.K.AND.K.EQ.ABS(L)) RETURN 538 IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN X=H Y=K Z=L H=-X K=-Y L=Z IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN H=-X K=Y L=-Z IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN H=X K=-Y L=-Z IF (H.GT.ABS(L).AND.K.GE.ABS(L)) RETURN H=Z K=X L=Y GO TO 538 C C M3 C 39 H=ABS(H) K=ABS(K) L=ABS(L) IF (H.EQ.K.AND.K.EQ.L) RETURN 539 IF (H.GT.L.AND.K.GE.L) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GT.L.AND.K.GE.L) RETURN H=Y K=Z L=X GO TO 539 C C 432 C 40 IF (ABS(H).EQ.ABS(K).AND.ABS(K).EQ.ABS(L)) THEN H=ABS(H) K=ABS(K) L=ABS(L) RETURN END IF 540 IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Y K=Z L=X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-X K=-Y L=Z IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Z K=-X L=-Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Y K=Z L=-X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-X K=Y L=-Z IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Z K=-X L=Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Y K=-Z L=-X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=X K=-Y L=-Z IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Z K=X L=-Y IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=-Y K=-Z L=X IF (H.GT.L.AND.K.GE.L.AND.L.GE.0) RETURN H=Y K=X L=-Z GO TO 540 C C -43M C 41 IF (H.GE.K.AND.K.GE.ABS(L)) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Y K=Z L=X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-X K=-Y L=Z IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Z K=-X L=-Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Y K=Z L=-X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-X K=Y L=-Z IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Z K=-X L=Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Y K=-Z L=-X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=X K=-Y L=-Z IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Z K=X L=-Y IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=-Y K=-Z L=X IF (H.GE.K.AND.K.GE.ABS(L)) RETURN H=Y K=X L=Z GO TO 41 C C M3M C 42 H=ABS(H) K=ABS(K) L=ABS(L) IF (H.GE.K.AND.K.GE.L) RETURN X=H Y=K Z=L H=Z K=X L=Y IF (H.GE.K.AND.K.GE.L) RETURN H=Y K=Z L=X IF (H.GE.K.AND.K.GE.L) RETURN H=Z K=Y L=X IF (H.GE.K.AND.K.GE.L) RETURN H=X K=Z L=Y IF (H.GE.K.AND.K.GE.L) RETURN H=Y K=X L=Z RETURN END C------------------------------------------------------------------------ SUBROUTINE ESYM (PTGP,H,K,L,E,M) C C RECIPROCAL LATTICE POINT DEGENERACIES, E = EPSILON(HKL), AND C MULTIPLICITIES, M = M(HKL), FOR THE CRYSTALLOGRAPHIC POINT GROUPS C C AS TABULATED BY: C C HITOSHI IWASAKI AND TETSUZO ITO (1977). VALUES OF EPSILON FOR C OBTAINING NORMALIZED STRUCTURE FACTORS. ACTA CRYST. A33, 227-229. C C TRICL. MONOCL. ORTHO. TETRAG. TRIG. RHOMB. CUB. C C (1) 1 (3) 2 (6) 222 (9) 4 (17) 3 (33) R3 (38) 23 C (2) -1 (4) M (7) MM2 (10) -4 (18) -3 (34) R-3 (39) M3 C (5) 2/M (8) MMM (11) 4/M C (19) 312 (35) R32 (40) 432 C (12) 422 (20) 321 (36) R3M (41) -43M C (13) 4MM (21) 31M (37) R-3M (42) M3M C (14) -42M (22) 3M1 C (15) -4M2 (23) -31M C (16) 4/MMM (24) -3M1 C C HEXAG. C C (25) 6 C (26) -6 C (27) 6/M C C (28) 622 C (29) 6MM C (30) -6M2 C (31) -62M C (32) 6/MMM C C A TOTAL OF 42 CASES IS TABULATED BECAUSE THERE ARE TEN CASES OF C ALTERNATIVE AXES FOR SEVEN OF THE 32 CRYSTALLOGRAPHIC POINT GROUPS: C C -42M -4M2 C 3 R3 C -3 R-3 C 312 321 R32 C 31M 3M1 R3M C -31M -3M1 R-3M C -62M -6M2 C INTEGER PTGP,H,E E=1 GO TO (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22, & 23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42) PTGP 1 M=1 GO TO 99 2 M=2 GO TO 99 3 M=2 IF (H.EQ.0.AND.L.EQ.0) E=2 GO TO 99 4 M=2 IF (K.EQ.0) E=2 GO TO 99 5 M=4 IF (K.EQ.0) E=2 IF (H.EQ.0.AND.L.EQ.0) E=2 GO TO 99 6 M=4 IF (K.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=2 GO TO 99 7 M=4 IF (H.EQ.0) E=2 IF (K.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 8 M=8 IF (H.EQ.0) E=2 IF (K.EQ.0) E=2 IF (L.EQ.0) E=2 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 9 M=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 10 M=4 IF (H.EQ.0.AND.K.EQ.0) E=2 GO TO 99 11 M=8 IF (L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 12 M=8 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 13 M=8 IF (K.EQ.0.OR.H.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=8 GO TO 99 14 M=8 IF (K.EQ.H.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 15 M=8 IF (K.EQ.0.OR.H.EQ.0) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 16 M=16 IF (K.EQ.0.OR.H.EQ.0.OR.L.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-H)) E=4 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=8 GO TO 99 17 M=3 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 18 M=6 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 19 M=6 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 20 M=6 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 21 M=6 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 22 M=6 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 23 M=12 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 24 M=12 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 25 M=6 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 26 M=6 IF (L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=3 GO TO 99 27 M=12 IF (L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 28 M=12 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=2 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 29 M=12 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.0) E=12 GO TO 99 30 M=12 IF (L.EQ.0) E=2 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 31 M=12 IF (L.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=6 GO TO 99 32 M=24 IF (L.EQ.0) E=2 IF (K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K) E=2 IF (L.EQ.0.AND.(K.EQ.H.OR.K.EQ.-2*H.OR.H.EQ.-2*K)) E=4 IF (K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H) E=2 IF (L.EQ.0.AND.(K.EQ.0.OR.H.EQ.0.OR.K.EQ.-H)) E=4 IF (H.EQ.0.AND.K.EQ.0) E=12 GO TO 99 33 M=3 IF (K.EQ.H.AND.L.EQ.H) E=3 GO TO 99 34 M=6 IF (K.EQ.H.AND.L.EQ.H) E=3 GO TO 99 35 M=6 IF (L.EQ.0.AND.K.EQ.-H) E=2 IF (K.EQ.0.AND.L.EQ.-H) E=2 IF (H.EQ.0.AND.K.EQ.-L) E=2 IF (K.EQ.H.AND.L.EQ.H) E=3 GO TO 99 36 M=6 IF (K.EQ.H.OR.L.EQ.H.OR.L.EQ.K) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=2 IF (K.EQ.H.AND.L.EQ.H) E=6 GO TO 99 37 M=12 IF (K.EQ.H.OR.L.EQ.H.OR.K.EQ.L) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=2 IF (K.EQ.-H.AND.L.EQ.0) E=2 IF (L.EQ.-H.AND.K.EQ.0) E=2 IF (L.EQ.-K.AND.H.EQ.0) E=2 IF (K.EQ.H.AND.L.EQ.H) E=6 GO TO 99 38 M=12 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3 IF (K.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.L.EQ.0) E=2 IF (H.EQ.0.AND.K.EQ.0) E=2 GO TO 99 39 M=24 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) E=2 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 40 M=24 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=3 IF (L.EQ.0.AND.ABS(K).EQ.ABS(H)) E=2 IF (K.EQ.0.AND.ABS(L).EQ.ABS(H)) E=2 IF (H.EQ.0.AND.ABS(K).EQ.ABS(L)) E=2 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 41 M=24 IF (ABS(K).EQ.ABS(H).OR.ABS(L).EQ.ABS(H).OR.ABS(K).EQ.ABS(L)) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=6 IF (K.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.L.EQ.0) E=4 IF (H.EQ.0.AND.K.EQ.0) E=4 GO TO 99 42 M=48 IF (ABS(K).EQ.ABS(H).OR.ABS(L).EQ.ABS(H).OR.ABS(K).EQ.ABS(L)) E=2 IF (ABS(K).EQ.ABS(H).AND.ABS(L).EQ.ABS(H)) E=6 IF (H.EQ.0.OR.K.EQ.0.OR.L.EQ.0) E=2 IF (L.EQ.0.AND.ABS(K).EQ.ABS(H)) E=4 IF (K.EQ.0.AND.ABS(L).EQ.ABS(H)) E=4 IF (H.EQ.0.AND.ABS(K).EQ.ABS(L)) E=4 IF (K.EQ.0.AND.L.EQ.0) E=8 IF (H.EQ.0.AND.L.EQ.0) E=8 IF (H.EQ.0.AND.K.EQ.0) E=8 99 M=M/E RETURN END C----------------------------------------------------------------------- SUBROUTINE EXPAND (GINV,BISO,B) C C EXPAND SCALAR BISO TO EQUIVALENT TENSOR BIJ'S. C DIMENSION GINV(3,3),B(3,3) DO I=1,3 DO J=I,3 B(I,J)=0.25*BISO*GINV(I,J) B(J,I)=B(I,J) END DO END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE FCALC (N,X,A1,B1,A2,B2,A3,B3,A4,B4,C0,FP,FPP,S, & SUMFSQ) DIMENSION X(1),A1(1),B1(1),A2(1),B2(1),A3(1),B3(1),A4(1),B4(1), & C0(1),FP(1),FPP(1) SUMFSQ=0 DO I=1,N SUMFSQ=SUMFSQ+X(I)*((A1(I)*EXP(-B1(I)*S**2) & +A2(I)*EXP(-B2(I)*S**2) & +A3(I)*EXP(-B3(I)*S**2) & +A4(I)*EXP(-B4(I)*S**2)+C0(I)+FP(I))**2 & +FPP(I)**2) END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE FTABLE (EL,Z,RM,A1,B1,A2,B2,A3,B3,A4,B4,C,XRAY,FP,FPP) C C NEUTRAL ATOM SCATTERING FACTOR COEFFICIENTS FROM CROMER & WABER C (INTERNATIONAL TABLES, VOL. IV, 1974) C C F(S) = A1*EXP(-B1*S**2) + A2*EXP(-B2*S**2) + A3*EXP(-B3*S**2) C + A4*EXP(-B4*S**2) + C C CHARACTER EL*2,SYMBOL*2,XRAY*2 DIMENSION SYMBOL(2*98),RMASS(98),AA1(98),BB1(98),AA2(98),BB2(98), & AA3(98),BB3(98),AA4(98),BB4(98),CC(98),FPMO(98),FPPMO(98), & FPCU(98),FPPCU(98),FPKISS(98) C C ELEMENT SYMBOLS C DATA SYMBOL / & 'H ', 'He', & 'Li','Be','B ','C ','N ','O ','F ','Ne', & 'Na','Mg','Al','Si','P ','S ','Cl','Ar', & 'K ','Ca', & 'Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn', & 'Ga','Ge','As','Se','Br','Kr', & 'Rb','Sr', & 'Y ','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd', & 'In','Sn','Sb','Te','I ','Xe', & 'Cs','Ba', & 'La', & 'Ce','Pr','Nd','Pm','Sm','Eu','Gd', & 'Tb','Dy','Ho','Er','Tm','Yb','Lu', & 'Hf','Ta','W ','Re','Os','Ir','Pt','Au','Hg', & 'Tl','Pb','Bi','Po','At','Rn', & 'Fr','Ra', & 'Ac', & 'Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf', C & 'H ', 'HE', & 'LI','BE','B ','C ','N ','O ','F ','NE', & 'NA','MG','AL','SI','P ','S ','CL','AR', & 'K ','CA', & 'SC','TI','V ','CR','MN','FE','CO','NI','CU','ZN', & 'GA','GE','AS','SE','BR','KR', & 'RB','SR', & 'Y ','ZR','NB','MO','TC','RU','RH','PD','AG','CD', & 'IN','SN','SB','TE','I ','XE', & 'CS','BA', & 'LA', & 'CE','PR','ND','PM','SM','EU','GD', & 'TB','DY','HO','ER','TM','YB','LU', & 'HF','TA','W ','RE','OS','IR','PT','AU','HG', & 'TL','PB','BI','PO','AT','RN', & 'FR','RA', & 'AC', & 'TH','PA','U ','NP','PU','AM','CM','BK','CF'/ C C RELATIVE ATOMIC MASSES (NATURAL ISOTOPIC ABUNDANCES) C DATA RMASS / 1.00797, 4.0026, 6.939, 9.0122, 10.811, 12.01115, & 14.0067, 15.9994, 18.9984, 20.183, 22.9898, 24.312, 26.9815, & 28.086, 30.9738, 32.064, 35.453, 39.948, 39.102, 40.08, 44.956, & 47.90, 50.942, 51.996, 54.938, 55.847, 58.933, 58.71, 63.54, & 65.37, 69.72, 72.59, 74.922, 78.96, 79.909, 83.80, 85.47, 87.62, & 88.905, 91.22, 92.906, 95.94, 98, 101.07, 102.905, 106.4, & 107.870, 112.40, 114.82, 118.69, 121.75, 127.60, 126.904, 131.30, & 132.905, 137.34, 138.91, 140.12, 140.907, 144.24, 147, 150.35, & 151.96, 157.25, 158.924, 162.50, 164.930, 167.26, 168.934, & 173.04, 174.97, 178.49, 180.948, 183.85, 186.2, 190.2, 192.2, & 195.09, 196.967, 200.59, 204.37, 207.19, 208.980, 210, 210, 222, & 223, 226, 227, 232.038, 231, 238.03, 237, 242, 243, 247, 247, & 249/ C C SCATTERING FACTOR COEFFICIENTS C DATA AA1 / 0.493002, 0.873400, 1.128200, 1.591900, 2.054500, & 2.310000, 12.212590, 3.048500, 3.539200, 3.955300, 4.762600, & 5.420400, 6.420200, 6.291500, 6.434500, 6.905300, 11.460390, & 7.484500, 8.218600, 8.626600, 9.189000, 9.759500, 10.297100, & 10.640600, 11.281900, 11.769490, 12.284090, 12.837590, 13.338000, & 14.074290, 15.235400, 16.081600, 16.672300, 17.000590, 17.178890, & 17.355490, 17.178400, 17.566290, 17.776000, 17.876490, 17.614190, & 3.702500, 19.130090, 19.267390, 19.295700, 19.331890, 19.280800, & 19.221400, 19.162390, 19.188900, 19.641790, 19.964400, 20.147200, & 20.293300, 20.389200, 20.336100, 20.578000, 21.167090, 22.044000, & 22.684490, 23.340490, 24.004190, 24.627390, 25.070900, 25.897590, & 26.507000, 26.904900, 27.656290, 28.181900, 28.664090, 28.947600, & 29.143990, 29.202390, 29.081800, 28.762100, 28.189400, 27.304900, & 27.005900, 16.881890, 20.680890, 27.544600, 31.061700, 33.368890, & 34.672600, 35.316290, 35.563090, 35.929900, 35.763000, 35.659690, & 35.564490, 35.884700, 36.022790, 36.187390, 36.525400, 36.670590, & 36.648800, 36.788100, 36.918500/ DATA BB1 / 10.510890, 9.103700, 3.954600, 43.642700, 23.218500, & 20.843900, 0.005700, 13.277090, 10.282500, 8.404200, 3.285000, & 2.827500, 3.038700, 2.438600, 1.906700, 1.467900, 0.010400, & 0.907200, 12.794890, 10.442090, 9.021300, 7.850800, 6.865700, & 6.103800, 5.340900, 4.761100, 4.279100, 3.878500, 3.582800, & 3.265500, 3.066900, 2.850900, 2.634500, 2.409800, 2.172300, & 1.938400, 1.788800, 1.556400, 1.402900, 1.276180, 1.188650, & 0.277200, 0.864132, 0.808520, 0.751536, 0.698655, 0.644600, & 0.594600, 0.547600, 5.830300, 5.303400, 4.817420, 4.347000, & 3.928200, 3.569000, 3.216000, 2.948170, 2.812190, 2.773930, & 2.662480, 2.562700, 2.472740, 2.387900, 2.253410, 2.242560, & 2.180200, 2.070510, 2.073560, 2.028590, 1.988900, 1.901820, & 1.832620, 1.773330, 1.720290, 1.671910, 1.629030, 1.592790, & 1.512930, 0.461100, 0.545000, 0.655150, 0.690200, 0.704000, & 0.700999, 0.685870, 0.663100, 0.646453, 0.616341, 0.589092, & 0.563359, 0.547751, 0.529300, 0.511929, 0.499384, 0.483629, & 0.465154, 0.451018, 0.437533/ DATA AA2 / 0.322912, 0.630900, 0.750800, 1.127800, 1.332600, & 1.020000, 3.132200, 2.286800, 2.641200, 3.112500, 3.173600, & 2.173500, 1.900200, 3.035300, 4.179100, 5.203400, 7.196400, & 6.772300, 7.439800, 7.387300, 7.367900, 7.355800, 7.351100, & 7.353700, 7.357300, 7.357300, 7.340900, 7.292000, 7.167600, & 7.031800, 6.700600, 6.374700, 6.070100, 5.819600, 5.235800, & 6.728600, 9.643500, 9.818400, 10.294590, 10.947990, 12.014390, & 17.235590, 11.094790, 12.918190, 14.350090, 15.501700, 16.688500, & 17.644390, 18.559600, 19.100490, 19.045500, 19.013790, 18.994900, & 19.029800, 19.106200, 19.296990, 19.598990, 19.769500, 19.669690, & 19.684690, 19.609490, 19.425790, 19.088590, 19.079800, 18.218500, & 17.638300, 17.294000, 16.428490, 15.885100, 15.434490, 15.220800, & 15.172590, 15.229290, 15.430000, 15.718890, 16.154990, 16.729590, & 17.763900, 18.591290, 19.041700, 19.158400, 13.063690, 12.951000, & 15.473290, 19.021100, 21.281600, 23.054700, 22.906400, 23.103190, & 23.421900, 23.294790, 23.412790, 23.596400, 23.808300, 24.099190, & 24.409600, 24.773600, 25.199490/ DATA BB2 / 26.125700, 3.356800, 1.052400, 1.862300, 1.021000, & 10.207500, 9.893300, 5.701100, 4.294400, 3.426200, 8.842200, & 79.261090, 0.742600, 32.333690, 27.156990, 22.215100, 1.166200, & 14.840700, 0.774800, 0.659900, 0.572900, 0.500000, 0.438500, & 0.392000, 0.343200, 0.307200, 0.278400, 0.256500, 0.247000, & 0.233300, 0.241200, 0.251600, 0.264700, 0.272600, 16.579600, & 16.562300, 17.315090, 14.098790, 12.800600, 11.916000, 11.765990, & 1.095800, 8.144870, 8.434670, 8.217580, 7.989290, 7.472600, & 6.908900, 6.377600, 0.503100, 0.460700, 0.420885, 0.381400, & 0.344000, 0.310700, 0.275600, 0.244475, 0.226836, 0.222087, & 0.210628, 0.202088, 0.196451, 0.194200, 0.181951, 0.196143, & 0.202172, 0.197940, 0.223545, 0.238849, 0.257119, 9.985190, & 9.599900, 9.370460, 9.225900, 9.092270, 8.979480, 8.865530, & 8.811740, 8.621600, 8.448400, 8.707510, 2.357600, 2.923800, & 3.550780, 3.974580, 4.069100, 4.176190, 3.871350, 3.651550, & 3.462040, 3.415190, 3.325300, 3.253960, 3.263710, 3.206470, & 3.089970, 3.046190, 3.007750/ DATA AA3 / 0.140191, 0.311200, 0.617500, 0.539100, 1.097900, & 1.588600, 2.012500, 1.546300, 1.517000, 1.454600, 1.267400, & 1.226900, 1.593600, 1.989100, 1.780000, 1.437900, 6.255600, & 0.653900, 1.051900, 1.589900, 1.640900, 1.699100, 2.070300, & 3.324000, 3.019300, 3.522200, 4.003400, 4.443800, 5.615800, & 5.165200, 4.359100, 3.706800, 3.431300, 3.973100, 5.637700, & 5.549300, 5.139900, 5.422000, 5.726290, 5.417320, 4.041830, & 12.887590, 4.649010, 4.863370, 4.734250, 5.295370, 4.804500, & 4.461000, 4.294800, 4.458500, 5.037100, 6.144870, 7.513800, & 8.976700, 10.661990, 10.887990, 11.372690, 11.851300, 12.385600, & 12.774000, 13.123490, 13.439590, 13.760290, 13.851790, 14.316690, & 14.559590, 14.558300, 14.977890, 15.154190, 15.308690, 15.100000, & 14.758600, 14.513500, 14.432700, 14.556400, 14.930500, 15.611490, & 15.713100, 25.558190, 21.657500, 15.538000, 18.442000, 16.587700, & 13.113800, 9.498870, 8.003700, 12.143890, 12.473890, 12.597700, & 12.747300, 14.189100, 14.949090, 15.640190, 16.770700, 17.341500, & 17.399000, 17.891900, 18.331690/ DATA BB3 / 3.142360, 22.927590, 85.390500,103.483000, 60.349790, & 0.568700, 28.997490, 0.323900, 0.261500, 0.230600, 0.313600, & 0.380800, 31.547190, 0.678500, 0.526000, 0.253600, 18.519390, & 43.898300,213.186900, 85.748390,136.108000, 35.633800, 26.893790, & 20.262600, 17.867400, 15.353500, 13.535900, 12.176300, 11.396590, & 10.316300, 10.780500, 11.446800, 12.947890, 15.237190, 0.260900, & 0.226100, 0.274800, 0.166400, 0.125599, 0.117622, 0.204785, & 11.003990, 21.570690, 24.799690, 25.874890, 25.205200, 24.660500, & 24.700800, 25.849890, 26.890890, 27.907390, 28.528390, 27.766000, & 26.465890, 24.387890, 20.207300, 18.772590, 17.608300, 16.766900, & 15.885000, 15.100890, 14.399600, 13.754590, 12.933090, 12.664790, & 12.189900, 11.440690, 11.360400, 10.997500, 10.664690, 0.261033, & 0.275116, 0.295977, 0.321703, 0.350500, 0.382661, 0.417916, & 0.424593, 1.482600, 1.572900, 1.963470, 8.618000, 8.793700, & 9.556420, 11.382390, 14.042200, 23.105190, 19.988690, 18.598990, & 17.830900, 16.923490, 16.092690, 15.362190, 14.945500, 14.313590, & 13.434590, 12.894590, 12.404390/ DATA AA4 / 0.040810, 0.178000, 0.465300, 0.702900, 0.706800, & 0.865000, 1.166300, 0.867000, 1.024300, 1.125100, 1.112800, & 2.307300, 1.964600, 1.541000, 1.490800, 1.586300, 1.645500, & 1.644200, 0.865900, 1.021100, 1.468000, 1.902100, 2.057100, & 1.492200, 2.244100, 2.304500, 2.348800, 2.380000, 1.673500, & 2.410000, 2.962300, 3.683000, 4.277900, 4.354300, 3.985100, & 3.537500, 1.529200, 2.669400, 3.265880, 3.657210, 3.533460, & 3.742900, 2.712630, 1.567560, 1.289180, 0.605844, 1.046300, & 1.602900, 2.039600, 2.466300, 2.682700, 2.523900, 2.273500, & 1.990000, 1.495300, 2.695900, 3.287190, 3.330490, 2.824280, & 2.851370, 2.875160, 2.896040, 2.922700, 3.545450, 2.953540, & 2.965770, 3.638370, 2.982330, 2.987060, 2.989630, 3.716010, & 4.300130, 4.764920, 5.119820, 5.441740, 5.675890, 5.833770, & 5.783700, 5.860000, 5.967600, 5.525930, 5.969600, 6.469200, & 7.025880, 7.425180, 7.443300, 2.112530, 3.210970, 4.086550, & 4.807030, 4.172870, 4.188000, 4.185500, 3.479470, 3.493310, & 4.216650, 4.232840, 4.243910/ DATA BB4 / 57.799690, 0.982100,168.261000, 0.542000, 0.140300, & 51.651190, 0.582600, 32.908900, 26.147590, 21.718390,129.423900, & 7.193700, 85.088590, 81.693690, 68.164500, 56.171990, 47.778390, & 33.392890, 41.684090,178.436900, 51.353100,116.104900,102.477900, & 98.739890, 83.754300, 76.880490, 71.169200, 66.342100, 64.812600, & 58.709700, 61.413490, 54.762490, 47.797190, 43.816290, 41.432800, & 39.397200,164.934000,132.376000,104.354000, 87.662700, 69.795700, & 61.658400, 86.847190, 94.292800, 98.606200, 76.898600, 99.815590, & 87.482490, 92.802900, 83.957100, 75.282500, 70.840300, 66.877590, & 64.265790,213.904000,167.201900,133.123900,127.113000,143.643900, &137.902900,132.720900,128.007000,123.173900,101.397900,115.361900, &111.873900, 92.656600,105.703000,102.960900,100.417000, 84.329800, & 72.029000, 63.364390, 57.055990, 52.086100, 48.164700, 45.001090, & 38.610300, 36.395590, 38.324600, 45.814890, 47.257900, 48.009290, & 47.004500, 45.471490, 44.247290,150.645000,142.324900,117.020000, & 99.172190,105.251000,100.613000, 97.490790,105.979900,102.272900, & 88.483390, 86.003000, 83.788100/ DATA CC / 0.003038, 0.006400, 0.037700, 0.038500, -0.193200, & 0.215600,-11.529000, 0.250800, 0.277600, 0.351500, 0.676000, & 0.858400, 1.115100, 1.140700, 1.114900, 0.866900, -9.557400, & 1.444500, 1.422800, 1.375100, 1.332900, 1.280700, 1.219900, & 1.183200, 1.089600, 1.036900, 1.011800, 1.034100, 1.191000, & 1.304100, 1.718900, 2.131300, 2.531000, 2.840900, 2.955700, & 2.825000, 3.487300, 2.506400, 1.912130, 2.069290, 3.755910, & 4.387500, 5.404280, 5.378740, 5.328000, 5.265930, 5.179000, & 5.069400, 4.939100, 4.782100, 4.590900, 4.352000, 4.071200, & 3.711800, 3.335200, 2.773100, 2.146780, 1.862640, 2.058300, & 1.984860, 2.028760, 2.209630, 2.574500, 2.419600, 3.583240, & 4.297280, 4.567960, 5.920460, 6.756210, 7.566720, 7.976280, & 8.581540, 9.243540, 9.887500, 10.472000, 11.000490, 11.472200, & 11.688300, 12.065790, 12.608900, 13.174590, 13.411800, 13.578200, & 13.677000, 13.710800, 13.690500, 13.724690, 13.621100, 13.526590, & 13.431400, 13.428700, 13.396590, 13.357290, 13.381190, 13.359190, & 13.288700, 13.275400, 13.267390/ C C ANOMALOUS DISPERSION CORRECTIONS FOR MO AND CU X-RAYS FROM CROMER & C LIBERMAN (INTERNATIONAL TABLES, VOL. IV, 1974) C DATA FPMO / & 0.000, 0.000, 0.000, 0.000, 0.000, 0.002, 0.004, 0.008, & 0.014, 0.021, 0.030, 0.042, 0.056, 0.072, 0.090, 0.110, & 0.132, 0.155, 0.179, 0.203, 0.226, 0.248, 0.267, 0.284, & 0.295, 0.301, 0.299, 0.285, 0.263, 0.222, 0.163, 0.081, & -0.030, -0.178, -0.374, -0.652, -1.044, -1.657, -2.951, -2.965, & -2.197, -1.825, -1.590, -1.420, -1.287, -1.177, -1.085, -1.005, & -0.936, -0.873, -0.816, -0.772, -0.726, -0.684, -0.644, -0.613, & -0.588, -0.564, -0.530, -0.535, -0.530, -0.533, -0.542, -0.564, & -0.591, -0.619, -0.666, -0.723, -0.795, -0.884, -0.988, -1.118, & -1.258, -1.421, -1.598, -1.816, -2.066, -2.352, -2.688, -3.084, & -3.556, -4.133, -4.861, -5.924, -7.444, -8.862, -7.912, -7.620, & -7.725, -8.127, -8.960,-10.673,-11.158, -9.725, -8.926, -8.416, & -7.990, -7.683/ DATA FPPMO / & 0.000, 0.000, 0.000, 0.000, 0.001, 0.002, 0.003, 0.006, & 0.010, 0.016, 0.025, 0.036, 0.052, 0.071, 0.095, 0.124, & 0.159, 0.201, 0.250, 0.306, 0.372, 0.446, 0.530, 0.624, & 0.729, 0.845, 0.973, 1.113, 1.266, 1.431, 1.609, 1.801, & 2.007, 2.223, 2.456, 2.713, 2.973, 3.264, 3.542, 0.560, & 0.621, 0.688, 0.759, 0.836, 0.919, 1.007, 1.101, 1.202, & 1.310, 1.424, 1.546, 1.675, 1.812, 1.958, 2.119, 2.282, & 2.452, 2.632, 2.845, 3.018, 3.225, 3.442, 3.669, 3.904, & 4.151, 4.410, 4.678, 4.958, 5.248, 5.548, 5.858, 6.185, & 6.523, 6.872, 7.232, 7.605, 7.990, 8.388, 8.798, 9.223, & 9.659, 10.102, 10.559, 11.042, 9.961, 10.403, 7.754, 8.105, & 8.472, 8.870, 9.284, 9.654, 4.148, 4.330, 4.511, 4.697, & 4.908, 5.107/ DATA FPCU / & 0.000, 0.000, 0.001, 0.003, 0.008, 0.017, 0.029, 0.047, & 0.069, 0.097, 0.129, 0.165, 0.204, 0.244, 0.283, 0.319, & 0.348, 0.366, 0.365, 0.341, 0.285, 0.189, 0.035, -0.198, & -0.568, -1.179, -2.464, -2.956, -2.019, -1.612, -1.354, -1.163, & -1.011, -0.879, -0.767, -0.665, -0.574, -0.465, -0.386, -0.314, & -0.248, -0.191, -0.145, -0.105, -0.077, -0.059, -0.060, -0.079, & -0.126, -0.194, -0.287, -0.418, -0.579, -0.783, -1.022, -1.334, & -1.716, -2.170, -2.939, -3.431, -4.357, -5.696, -7.718, -9.242, & -9.498,-10.423,-12.255, -9.733, -8.488, -7.701, -7.713, -6.715, & -6.351, -6.048, -5.790, -5.581, -5.391, -5.233, -5.096, -4.990, & -4.883, -4.818, -4.776, -4.756, -4.772, -4.787, -4.833, -4.898, & -4.994, -5.091, -5.216, -5.359, -5.529, -5.712, -5.930, -6.176, & -6.498, -6.798/ DATA FPPCU / & 0.000, 0.000, 0.000, 0.001, 0.004, 0.009, 0.018, 0.032, & 0.053, 0.083, 0.124, 0.177, 0.246, 0.330, 0.434, 0.557, & 0.702, 0.872, 1.066, 1.286, 1.533, 1.807, 2.110, 2.443, & 2.808, 3.204, 3.608, 0.509, 0.589, 0.678, 0.777, 0.886, & 1.006, 1.139, 1.283, 1.439, 1.608, 1.820, 2.025, 2.245, & 2.482, 2.735, 3.005, 3.296, 3.605, 3.934, 4.282, 4.653, & 5.045, 5.459, 5.894, 6.352, 6.835, 7.348, 7.904, 8.460, & 9.036, 9.648, 10.535, 10.933, 11.614, 12.320, 11.276, 11.946, & 9.242, 9.748, 3.704, 3.937, 4.181, 4.432, 4.693, 4.977, & 5.271, 5.577, 5.891, 6.221, 6.566, 6.925, 7.297, 7.686, & 8.089, 8.505, 8.930, 9.383, 9.843, 10.317, 10.803, 11.296, & 11.799, 12.330, 12.868, 13.409, 13.967, 14.536, 15.087, 15.634, & 16.317, 16.930/ C C WAVELENGTH-INDEPENDENT KISSEL AND PRATT CORRECTIONS TO THE CROMER AND C LIBERMAN FP VALUES [L. KISSEL AND R.H. PRATT (1990). ACTA CRYST. A46, C 170-175]. C DATA FPKISS / & 0.000, 0.000, 0.001, 0.000, 0.001, 0.001, 0.002, 0.003, & 0.004, 0.004, 0.006, 0.008, 0.008, 0.011, 0.012, 0.014, & 0.017, 0.020, 0.022, 0.025, 0.028, 0.031, 0.035, 0.039, & 0.042, 0.048, 0.052, 0.057, 0.061, 0.067, 0.073, 0.079, & 0.085, 0.092, 0.099, 0.106, 0.114, 0.122, 0.130, 0.138, & 0.147, 0.156, 0.166, 0.175, 0.186, 0.196, 0.207, 0.219, & 0.230, 0.242, 0.255, 0.267, 0.281, 0.294, 0.308, 0.323, & 0.338, 0.354, 0.369, 0.386, 0.402, 0.419, 0.437, 0.455, & 0.474, 0.493, 0.512, 0.532, 0.553, 0.574, 0.596, 0.617, & 0.640, 0.663, 0.687, 0.711, 0.736, 0.762, 0.788, 0.814, & 0.842, 0.870, 0.899, 0.928, 0.957, 0.988, 1.018, 1.050, & 1.083, 1.115, 1.149, 1.184, 1.219, 1.255, 1.292, 1.330, & 1.368, 1.407/ IF (EL.EQ.'J ') EL='I ' DO I=1,2*98 IF (EL.EQ.SYMBOL(I)) GO TO 1 END DO STOP 'ELEMENT SYMBOL IS NOT IN THE TABLE (Z .LE. 98)' 1 CONTINUE IF (I.GT.98) I=I-98 Z=I RM=RMASS(I) A1=AA1(I) B1=BB1(I) A2=AA2(I) B2=BB2(I) A3=AA3(I) B3=BB3(I) A4=AA4(I) B4=BB4(I) C=CC(I) IF (XRAY.EQ.'MO') THEN FP=FPMO(I)+FPKISS(I) FPP=FPPMO(I) ELSE IF (XRAY.EQ.'CU') THEN FP=FPCU(I)+FPKISS(I) FPP=FPPCU(I) ELSE FP=0 FPP=0 END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE KBCISO (IFILE,NDATA,SCALEK,BISO,CISO,NCYCLE,R,Z,NFREE, & NPAR,COV) C C NON-LINEAR LEAST-SQUARES REFINEMENT OF SCALE AND ISOTROPIC C DISPLACEMENT AND DISTRIBUTION DISPERSION PARAMETERS C C YOBS = FSQ(MEAS)/[EPSILON*SUM F(J)**2] C YCALC = SCALEK**(-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2} C DIMENSION AA(3,3),XX(3),BB(3),COV(13,13) DOUBLE PRECISION AA MPAR=3 NPAR=3 C C CALCULATE INITIAL AGREEMENT STATISTICS. C CALL AGREE (IFILE,NDATA,SCALEK,BISO,CISO,NPAR,NFREE,Z,R) C C ITERATE THROUGH LEAST-SQUARES LOOP. C NCYCLE=0 1 NCYCLE=NCYCLE+1 C C STORE THE AGREEMENT INDICES FROM THE PRECEDING CYCLE. C ROLD=R ZOLD=Z C C REFINE DISPERSION PARAMETER ONLY IF IT IS NONZERO. C IF (CISO.EQ.0) NPAR=2 C C EVALUATE THE DERIVATIVES, AND BUILD THE NORMAL MATRIX AND VECTOR. C DO J=1,NPAR BB(J)=0 DO K=1,NPAR AA(J,K)=0 END DO END DO REWIND IFILE DO I=1,NDATA READ (IFILE) IH,IK,IL,SSQ,YI,SIGYI,MULT,VAR Y=(1/SCALEK**2)*EXP(-2*(BISO*SSQ-(CISO*SSQ)**2)) DYDK=-2*Y/SCALEK DYDB=-2*Y*SSQ DYDS=+4*Y*SSQ**2*CISO XX(1)=DYDK XX(2)=DYDB XX(3)=DYDS WI=MULT/(SIGYI**2+Y**2*VAR) DO J=1,NPAR DO K=J,NPAR AA(J,K)=AA(J,K)+WI*XX(J)*XX(K) AA(K,J)=AA(J,K) END DO BB(J)=BB(J)+WI*XX(J)*(YI-Y) END DO END DO C C INVERT THE NORMAL MATRIX. C CALL MATINV (NPAR,MPAR,AA,DET) IF (DET.EQ.0) STOP 'LEVY/KBCISO: SINGULAR NORMAL MATRIX' C C CALCULATE THE FULL PARAMETER SHIFTS. C DO I=1,NPAR XX(I)=0 DO J=1,NPAR XX(I)=XX(I)+AA(I,J)*BB(J) END DO END DO IF (NPAR.LT.MPAR) THEN DO I=NPAR+1,MPAR XX(I)=0 END DO END IF C C INITIALIZE THE SHIFT-DAMPING FACTOR TO UNITY. C F=2 2 F=0.5*F IF (F.LT.1E-6) THEN NCYCLE=NCYCLE-1 GO TO 9 END IF C C APPLY THE DAMPED PARAMETER SHIFTS. C TK=SCALEK+F*XX(1) TB=BISO +F*XX(2) IF (NPAR.LE.2) THEN TC=0 ELSE TC=CISO+F*XX(3) END IF C C REQUIRE POSITIVE SCALEK AND NON-NEGATIVE BISO AND CISO. C IF (TK.LE.0.OR.TB.LT.0.OR.TC.LT.0) GO TO 2 C C FIT A PARABOLA TO THE AGREEMENT INDICES FOR ZERO, HALF, AND FULL C SHIFTS TO ESTIMATE AN OPTIMUM SHIFT-DAMPING FACTOR. C C SHIFT-DAMPING FACTORS, 0 .LT. F .LT. 1, MIGHT BE NECESSARY BECAUSE THE C WEIGHTS, AS WELL AS THE RESIDUALS, DEPEND ON THE FITTED PARAMETERS. C THIS COMPOUNDS THE NON-LINEARITY OF THE LEAST-SQUARES PROBLEM, AND C OVER-SHIFTS CAN LEAD TO A DIVERGENT REFINEMENT. C CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R) F0=0 R0=ROLD F1=F R1=R X=0.5*F N=0 31 N=N+1 IF (N.LT.10) THEN F=0.5*(F0+F1) ELSE IF (F.LT.X) F=F0 IF (F.GT.X) F=F1 GO TO 39 END IF TK=SCALEK+F*XX(1) TB=BISO +F*XX(2) IF (NPAR.GT.2) TC=CISO+F*XX(3) CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R) RMIN=MIN(R0,R,R1) IF (RMIN.EQ.R0) THEN R1=R F1=F GO TO 31 ELSE IF (RMIN.EQ.R1) THEN R0=R F0=F GO TO 31 ELSE CALL PMIN (F0,R0,F,R,F1,R1,FMIN,RMIN) F=FMIN END IF 39 TK=SCALEK+F*XX(1) TB=BISO +F*XX(2) IF (NPAR.GT.2) TC=CISO+F*XX(3) CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R) C C TEST FOR CONVERGENCE. C IF (R.LT.ROLD) THEN C C CONVERGING OR CONVERGED. STORE SHIFTED PARAMETERS. C SCALEK=TK BISO=TB IF (NPAR.GT.2) CISO=TC C C EVALUATE THE VARIANCE-COVARIANCE MATRIX. C ZSQ=Z**2 DO I=1,NPAR DO J=1,I COV(I,J)=ZSQ*AA(I,J) COV(J,I)=COV(I,J) END DO END DO C C TEST MAXIMUM SHIFT-TO-ERROR RATIO. C X=0 DO I=1,NPAR X=MAX(X,F*ABS(XX(I))/SQRT(COV(I,I))) END DO IF (X.LT.1E-6) THEN GO TO 9 ELSE IF (NCYCLE.LE.20) THEN GO TO 1 END IF ELSE C C CONVERGED OR DIVERGING. C NCYCLE=NCYCLE-1 R=ROLD Z=ZOLD END IF 9 CONTINUE DO I=NPAR+1,13 DO J=1,I COV(I,J)=0 COV(J,I)=0 END DO END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE KBCIJ (IFILE,NDATA,G,GINV,SCALEK,B,C,SYSTEM,NCYCLE,R,Z, & NFREE,NPAR,COV) C C REFINE THE FIT OF THE SCALE FACTOR, OVERALL ANISOTROPIC ATOMIC C DISPLACEMENT PARAMETERS, AND DISTRIBUTION DISPERSION PARAMETERS. C C YOBS = FSQ(MEAS)/[EPSILON*SUM F(J)**2] C YCAL = SCALEK**(-2)*EXP(-2*(H(J)*B(I,J)*H(I) - (H(J)*C(I,J)*H(I))**2)) C CHARACTER SYSTEM*12 DIMENSION G(3,3),GINV(3,3),H(3),B(3,3),C(3,3),AA(13,13), & XX(13),BB(13),COV(13,13),TB(3,3),TC(3,3) DOUBLE PRECISION AA MPAR=13 NPAR=13 C C CALCULATE INITIAL AGREEMENT STATISTICS. C CALL AGREE (IFILE,NDATA,SCALEK,B,C,NPAR,NFREE,Z,R) C C ITERATE THROUGH LEAST-SQUARES LOOP. C NCYCLE=0 1 NCYCLE=NCYCLE+1 C C STORE THE AGREEMENT INDICES FROM THE PRECEDING CYCLE. C ROLD=R ZOLD=Z C C REFINE DISPERSION PARAMETERS ONLY IF THE C-TENSOR IS NONZERO. C IF (C(1,1).EQ.0) NPAR=7 C C EVALUATE THE DERIVATIVES, AND BUILD THE NORMAL MATRIX AND VECTOR. C DO J=1,NPAR BB(J)=0 DO K=1,NPAR AA(J,K)=0 END DO END DO REWIND IFILE DO I=1,NDATA READ (IFILE) IH,IK,IL,SSQ,YI,SIGYI,MULT,VAR H(1)=IH H(2)=IK H(3)=IL HTBH=VTMV(H,B) HTCH=VTMV(H,C) Y=(1/SCALEK**2)*EXP(-2*(HTBH-HTCH**2)) DYDK=-2*Y/SCALEK XX(1)=DYDK L=1 DO J=1,3 DO K=J,3 DJK=(2-DELTA(J,K))*H(J)*H(K) DYDB=-2*Y*DJK DYDC=+4*Y*DJK*HTCH L=L+1 XX(L)=DYDB XX(L+6)=DYDC END DO END DO WI=MULT/(SIGYI**2+Y**2*VAR) DO J=1,NPAR DO K=J,NPAR AA(J,K)=AA(J,K)+WI*XX(J)*XX(K) AA(K,J)=AA(J,K) END DO BB(J)=BB(J)+WI*XX(J)*(YI-Y) END DO END DO C C INVERT THE NORMAL MATRIX. C CALL MATINV (NPAR,MPAR,AA,DET) IF (DET.EQ.0) STOP 'LEVY/KBCIJ: SINGULAR NORMAL MATRIX' C C CALCULATE THE FULL PARAMETER SHIFTS. C DO I=1,NPAR XX(I)=0 DO J=1,NPAR XX(I)=XX(I)+AA(I,J)*BB(J) END DO END DO IF (NPAR.LT.MPAR) THEN DO I=NPAR+1,MPAR XX(I)=0 END DO END IF C C INITIALIZE THE SHIFT-DAMPING FACTOR TO UNITY. C F=2 2 F=0.5*F IF (F.LT.1E-6) THEN NCYCLE=NCYCLE-1 GO TO 9 END IF C C APPLY THE DAMPED PARAMETER SHIFTS. C TK=SCALEK+F*XX(1) C C REQUIRE POSITIVE SCALEK. C IF (TK.LE.0) GO TO 2 L=1 DO J=1,3 DO K=J,3 L=L+1 TB(J,K)=B(J,K)+F*XX(L) TC(J,K)=0 IF (NPAR.GT.7) TC(J,K)=C(J,K)+F*XX(L+6) END DO END DO CALL RESET (SYSTEM,TB) CALL RESET (SYSTEM,TC) C C REQUIRE POSITIVE DEFINITE B- AND C-TENSORS. C CALL POSDEF (TB,G,GINV,SYSTEM) IF (TB(1,1).LE.0) GO TO 2 IF (NPAR.GT.7) THEN CALL POSDEF (TC,G,GINV,SYSTEM) IF (TC(1,1).LE.0) GO TO 2 END IF C C FIT A PARABOLA TO THE AGREEMENT INDICES FOR ZERO, HALF, AND FULL C SHIFTS TO ESTIMATE AN OPTIMUM SHIFT-DAMPING FACTOR. C C SHIFT-DAMPING FACTORS, 0 .LT. F .LT. 1, MIGHT BE NECESSARY BECAUSE THE C WEIGHTS, AS WELL AS THE RESIDUALS, DEPEND ON THE FITTED PARAMETERS. C THIS COMPOUNDS THE NON-LINEARITY OF THE LEAST-SQUARES PROBLEM, AND C OVER-SHIFTS CAN LEAD TO A DIVERGENT REFINEMENT. C CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R) F0=0 R0=ROLD F1=F R1=R X=0.5*F N=0 31 N=N+1 IF (N.LT.10) THEN F=0.5*(F0+F1) ELSE IF (F.LT.X) F=F0 IF (F.GT.X) F=F1 GO TO 39 END IF TK=SCALEK+F*XX(1) L=1 DO J=1,3 DO K=J,3 L=L+1 TB(J,K)=B(J,K)+F*XX(L) IF (NPAR.GT.7) TC(J,K)=C(J,K)+F*XX(L+6) END DO END DO CALL RESET (SYSTEM,TB) CALL RESET (SYSTEM,TC) CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R) RMIN=MIN(R0,R,R1) IF (RMIN.EQ.R0) THEN R1=R F1=F GO TO 31 ELSE IF (RMIN.EQ.R1) THEN R0=R F0=F GO TO 31 ELSE CALL PMIN (F0,R0,F,R,F1,R1,FMIN,RMIN) F=FMIN END IF 39 TK=SCALEK+F*XX(1) L=1 DO J=1,3 DO K=J,3 L=L+1 TB(J,K)=B(J,K)+F*XX(L) IF (NPAR.GT.7) TC(J,K)=C(J,K)+F*XX(L+6) END DO END DO CALL RESET (SYSTEM,TB) CALL RESET (SYSTEM,TC) CALL AGREE (IFILE,NDATA,TK,TB,TC,NPAR,NFREE,Z,R) C C TEST FOR CONVERGENCE. C IF (R.LT.ROLD) THEN C C CONVERGING OR CONVERGED. STORE SHIFTED PARAMETERS. C SCALEK=TK DO J=1,3 DO K=1,3 B(J,K)=TB(J,K) IF (NPAR.GT.7) C(J,K)=TC(J,K) END DO END DO C C EVALUATE THE VARIANCE-COVARIANCE MATRIX. C DO I=1,NPAR DO J=1,I COV(I,J)=Z**2*AA(I,J) COV(J,I)=COV(I,J) END DO END DO C C TEST MAXIMUM SHIFT-TO-ERROR RATIO. C X=0 DO I=1,NPAR X=MAX(X,F*ABS(XX(I))/SQRT(COV(I,I))) END DO IF (X.LT.1E-6) THEN GO TO 9 ELSE IF (NCYCLE.LE.20) THEN GO TO 1 END IF ELSE C C CONVERGED OR DIVERGING. C NCYCLE=NCYCLE-1 R=ROLD Z=ZOLD END IF 9 CONTINUE IF (NPAR.LT.13) THEN DO I=NPAR+1,13 DO J=1,I COV(I,J)=0 COV(J,I)=0 END DO END DO END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE MATINV (N,M,A,D) C C INVERT A SQUARE MATRIX BY GAUSS-JORDAN ELIMINATION, AND CALCULATE ITS C DETERMINANT. C C A - INPUT MATRIX, WHICH IS REPLACED BY ITS INVERSE C N - ORDER OF MATRIX C - LOGICAL SIZE OF ARRAY A C M - PHYSICAL SIZE OF ARRAY A, M .GE. N C D - DETERMINANT OF THE INPUT MATRIX C C RETURNS D = 0 TO FLAG A SINGULAR MATRIX. C C MATRIX IS SCALED TO AVOID ARITHMETIC OVERFLOW OR UNDERFLOW DURING C INVERSION. C C MATRIX SCALING USES A DIAGONAL MATRIX Q WITH C C Q(I,I) = 1/SQRT(ABS(A(I,I))). C C FORM B = Q A Q, C C B(I,J) = Q(I,I)*A(I,J)*Q(J,J), C C AND INVERT B TO OBTAIN B**-1. THEN, C C B**-1 = (Q A Q)**-1 C = (A Q)**-1 Q**-1 C B**-1 Q = Q**-1 A**-1 C Q B**-1 Q = A**-1 C C IN ADDITION, C C DET(A) = DET(B)/PROD(I=1,N) Q(I,I)**2 C C IF N EXCEEDS MAXN, THE WORK VECTORS Q, IK, AND JK MUST BE RE- C DIMENSIONED TO A LENGTH OF N. C PARAMETER (MAXN=25) DIMENSION A(M,M),Q(MAXN),IK(MAXN),JK(MAXN) DOUBLE PRECISION A,Q,AMAX,T C C STORE RECIPROCAL SQUARE-ROOT DIAGONAL MAGNITUDES FOR SCALING. C DO 10 I=1,N T=SQRT(ABS(A(I,I))) IF (T.EQ.0) THEN D=0 RETURN ELSE Q(I)=1/T END IF 10 CONTINUE C C SCALE MATRIX. C DO 11 I=1,N DO 11 J=1,N A(I,J)=A(I,J)*Q(I)*Q(J) 11 CONTINUE C C INVERT SCALED MATRIX AND CALCULATE ITS DETERMINANT. C D=1 DO 100 K=1,N C C FIND LARGEST ELEMENT A(I,J) IN REST OF MATRIX. C AMAX=0 DO 30 I=K,N DO 30 J=K,N IF (ABS(A(I,J)).GT.ABS(AMAX)) THEN AMAX=A(I,J) IK(K)=I JK(K)=J END IF 30 CONTINUE D=D*AMAX IF (D.EQ.0) RETURN C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN A(K,K). C I=IK(K) IF (I.GT.K) THEN DO 50 J=1,N T=A(K,J) A(K,J)=A(I,J) A(I,J)=-T 50 CONTINUE END IF J=JK(K) IF (J.GT.K) THEN DO 60 I=1,N T=A(I,K) A(I,K)=A(I,J) A(I,J)=-T 60 CONTINUE END IF C C ACCUMULATE ELEMENTS OF INVERSE MATRIX. C DO 70 I=1,N IF (I.NE.K) A(I,K)=-A(I,K)/AMAX 70 CONTINUE DO 80 I=1,N DO 80 J=1,N IF (I.NE.K.AND.J.NE.K) A(I,J)=A(I,J)+A(I,K)*A(K,J) 80 CONTINUE DO 90 J=1,N IF (J.NE.K) A(K,J)=+A(K,J)/AMAX 90 CONTINUE A(K,K)=1/AMAX 100 CONTINUE C C RESTORE ORDERING OF MATRIX. C DO 130 K=N,1,-1 J=IK(K) IF (J.GT.K) THEN DO 110 I=1,N T=A(I,K) A(I,K)=-A(I,J) A(I,J)=T 110 CONTINUE END IF I=JK(K) IF (I.GT.K) THEN DO 120 J=1,N T=A(K,J) A(K,J)=-A(I,J) A(I,J)=T 120 CONTINUE END IF 130 CONTINUE C C SCALE INVERSE MATRIX. C DO 140 I=1,N DO 140 J=1,N A(I,J)=A(I,J)*Q(I)*Q(J) 140 CONTINUE C C SCALE DETERMINANT, IF POSSIBLE. C T=0 DO 150 I=1,N T=T+LOG(Q(I)) 150 CONTINUE T=LOG(ABS(D))-2*T C C TEST AGAINST RANGE OF MACHINE ALLOWED MAGNITUDES. C C EIGHT-BIT EXPONENT HAS MAXIMUM MAGNITUDE 2**7 = 128. C XMIN = 2**(-128) = 0.29E-38 LOG(XMIN) = -88.7 C XMAX = (2**(+128))/2 = 1.70E+38 LOG(XMAX) = +88.0 C IF (-87.LT.T.AND.T.LT.+87) D=(D/ABS(D))*EXP(T) RETURN END C----------------------------------------------------------------------- SUBROUTINE MATINV3 (A,B,D) C C INVERT SYMMETRICAL (3 X 3) MATRIX. C C A IS THE MATRIX; B IS ITS INVERSE; D IS ITS DETERMINANT. C DIMENSION A(3,3),B(3,3) DOUBLE PRECISION A,B,D DO I=1,3 K=1+MOD(I+1,3) M=1+MOD(I+3,3) DO J=1,3 L=1+MOD(J+1,3) N=1+MOD(J+3,3) B(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K) END DO END DO D=0 DO I=1,3 D=D+A(I,1)*B(1,I) END DO DO I=1,3 DO J=1,3 B(I,J)=B(I,J)/D END DO END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE METRIC (A,V,GA,GB) C C CALCULATE RECIPROCAL LATTICE PARAMETERS AND DIRECT AND RECIPROCAL C LATTICE METRIC TENSORS FROM DIRECT LATTICE PARAMETERS. C DIMENSION A(6),B(6),GA(3,3),GB(3,3) PI=ACOS(-1.0) RAD=PI/180 DEG=180/PI CA=COS(A(4)*RAD) CB=COS(A(5)*RAD) CG=COS(A(6)*RAD) SA=SIN(A(4)*RAD) SB=SIN(A(5)*RAD) SG=SIN(A(6)*RAD) V=A(1)*A(2)*A(3)*SQRT(1-CA**2-CB**2-CG**2+2*CA*CB*CG) B(1)=A(2)*A(3)*SA/V B(2)=A(1)*A(3)*SB/V B(3)=A(1)*A(2)*SG/V B(4)=(CB*CG-CA)/(SB*SG) B(5)=(CA*CG-CB)/(SA*SG) B(6)=(CA*CB-CG)/(SA*SB) GA(1,1)=A(1)*A(1) GA(1,2)=A(1)*A(2)*CG GA(1,3)=A(1)*A(3)*CB GA(2,1)=GA(1,2) GA(2,2)=A(2)*A(2) GA(2,3)=A(2)*A(3)*CA GA(3,1)=GA(1,3) GA(3,2)=GA(2,3) GA(3,3)=A(3)*A(3) GB(1,1)=B(1)*B(1) GB(1,2)=B(1)*B(2)*B(6) GB(1,3)=B(1)*B(3)*B(5) GB(2,1)=GB(1,2) GB(2,2)=B(2)*B(2) GB(2,3)=B(2)*B(3)*B(4) GB(3,1)=GB(1,3) GB(3,2)=GB(2,3) GB(3,3)=B(3)*B(3) B(4)=ACOS(B(4))*DEG B(5)=ACOS(B(5))*DEG B(6)=ACOS(B(6))*DEG RETURN END C----------------------------------------------------------------------- SUBROUTINE PMIN (X1,Y1,X2,Y2,X3,Y3,XMIN,YMIN) C C CALCULATE THE COEFFICIENTS OF THE PARABOLA, C Y = C0 + C1*X + C2*X**2, C THROUGH THREE POINTS, C (X1,Y1), (X2,Y2), AND (X3,Y3) WITH Y2 .LT. Y1, Y3, C AND FIND THE MINIMUM OF THE CURVE. C DIMENSION A(3,3),AINV(3,3),B(3),C(3) DOUBLE PRECISION A,AINV,DETA A(1,1)=1 A(1,2)=X1 A(1,3)=X1*X1 A(2,1)=1 A(2,2)=X2 A(2,3)=X2*X2 A(3,1)=1 A(3,2)=X3 A(3,3)=X3*X3 B(1)=Y1 B(2)=Y2 B(3)=Y3 CALL MATINV3 (A,AINV,DETA) DO I=1,3 C(I)=0 DO J=1,3 C(I)=C(I)+AINV(I,J)*B(J) END DO END DO XMIN=X2 IF (C(3).GT.0) THEN XMIN=-0.5*C(2)/C(3) IF (XMIN.LT.X1) XMIN=X1 IF (XMIN.GT.X3) XMIN=X3 YMIN=C(1)+C(2)*XMIN+C(3)*XMIN*XMIN ELSE XMIN=X2 YMIN=Y2 END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE POSDEF (B,G,GINV,SYSTEM) C C TEST NECESSARY AND SUFFICIENT CONDITION FOR POSITIVE-DEFINITE, C SYMMETRIC B-TENSOR. C CHARACTER SYSTEM*12 DIMENSION B(3,3),G(3,3),GINV(3,3) T1=B(1,1) T2=B(1,1)*B(2,2)-B(2,1)*B(1,2) T3=B(1,1)*B(2,2)*B(3,3)+B(1,2)*B(2,3)*B(3,1)+B(1,3)*B(2,1)*B(3,2) & -B(3,1)*B(2,2)*B(1,3)-B(3,2)*B(2,3)*B(1,1)-B(3,3)*B(2,1)*B(1,2) IF (MIN(T1,T2,T3).LE.0) THEN CALL CONTRACT (G,B,BISO) IF (BISO.LT.0) BISO=0 CALL EXPAND (GINV,BISO,B) CALL RESET (SYSTEM,B) END IF RETURN END C----------------------------------------------------------------------- SUBROUTINE READ1 (ITYPE,IUNIT,IEND,IH,IK,IL,FSQ,SIGFSQ) 1 IF (ITYPE.EQ.0) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF IF (ITYPE.EQ.1) READ (IUNIT,*,END=9) IH,IK,IL,FSQ,SIGFSQ IF (ITYPE.EQ.2) READ (IUNIT,*,END=9) IH,IK,IL,F ,SIGF IF (ITYPE.EQ.3) READ (IUNIT, END=9) IH,IK,IL,FSQ,SIGFSQ,F,SIGF IF (ITYPE.EQ.4) READ (IUNIT, END=9) IH,IK,IL,FSQ,SIGFSQ IF (ITYPE.EQ.5) READ (IUNIT, END=9) IH,IK,IL,F ,SIGF IF (ITYPE.EQ.2.OR.ITYPE.EQ.5) THEN IF (F.LT.0) F=0 IF (SIGF.LT.0) SIGF=0 FSQ=F**2 SIGFSQ=MAX(2*F*SIGF,4*SIGF**2) END IF IF (SIGFSQ.LE.0) GO TO 1 RETURN 9 IEND=1 RETURN END C----------------------------------------------------------------------- SUBROUTINE RECALL (T,NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV) DIMENSION T(109),B(3,3),C(3,3),COV(13,13) NCYCLE=T(1) R =T(2) Z =T(3) NFREE =T(4) NPAR =T(5) SCALEK=T(6) I=6 DO J=1,3 DO K=J,3 I=I+1 B(J,K)=T(I) END DO END DO DO J=1,3 DO K=J,3 I=I+1 C(J,K)=T(I) END DO END DO DO J=1,13 DO K=J,13 I=I+1 COV(J,K)=T(I) END DO END DO DO I=1,3 DO J=I,3 B(J,I)=B(I,J) C(J,I)=C(I,J) END DO END DO DO I=1,13 DO J=1,13 COV(J,I)=COV(I,J) END DO END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE RESET (SYSTEM,B) CHARACTER SYSTEM*12 DIMENSION B(3,3) IF (SYSTEM.EQ.'MONOCLINIC') THEN B(1,2)=0 B(2,3)=0 ELSE IF (SYSTEM.EQ.'ORTHORHOMBIC') THEN B(1,2)=0 B(1,3)=0 B(2,3)=0 ELSE IF (SYSTEM.EQ.'TETRAGONAL') THEN Q=(B(1,1)+B(2,2))/2 B(1,1)=Q B(2,2)=Q B(1,2)=0 B(1,3)=0 B(2,3)=0 ELSE IF (SYSTEM.EQ.'TRIGONAL'.OR.SYSTEM.EQ.'HEXAGONAL') THEN Q=(B(1,1)+B(2,2))/2 B(1,1)=Q B(2,2)=Q B(1,2)=Q/2 B(1,3)=0 B(2,3)=0 ELSE IF (SYSTEM.EQ.'RHOMBOHEDRAL') THEN Q=(B(1,1)+B(2,2)+B(3,3))/3 B(1,1)=Q B(2,2)=Q B(3,3)=Q Q=(B(1,2)+B(1,3)+B(2,3))/3 B(1,2)=Q B(1,3)=Q B(2,3)=Q ELSE IF (SYSTEM.EQ.'CUBIC') THEN Q=(B(1,1)+B(2,2)+B(3,3))/3 B(1,1)=Q B(2,2)=Q B(3,3)=Q B(1,2)=0 B(1,3)=0 B(2,3)=0 END IF B(2,1)=B(1,2) B(3,1)=B(1,3) B(3,2)=B(2,3) RETURN END C----------------------------------------------------------------------- SUBROUTINE SORT (N,DATA,INDEX) C C INDEXES THE ARRAY DATA(N) AND RETURNS THE ARRAY INDEX(N) SORTED SUCH C THAT THE VALUES DATA(INDEX(I)) ARE IN ASCENDING ORDER FOR I = 1, 2, C ..., N. THE INPUT VARIABLES N AND DATA(N) ARE NOT CHANGED. C C EMPLOYS THE "HEAPSORT" ALGORITHM. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C ART OF SCIENTIFIC COMPUTING, PP. 229-233. CAMBRIDGE, ENGLAND: C CAMBRIDGE UNIVERSITY PRESS. C DIMENSION DATA(N),INDEX(N) DO 10 I=1,N INDEX(I)=I 10 CONTINUE I=N/2+1 J=N 1 CONTINUE IF (I.GT.1) THEN I=I-1 INDEXT=INDEX(I) T=DATA(INDEXT) ELSE INDEXT=INDEX(J) T=DATA(INDEXT) INDEX(J)=INDEX(1) J=J-1 IF (J.EQ.1) THEN INDEX(1)=INDEXT RETURN END IF END IF K=I L=I+I 2 IF (L.LE.J) THEN IF (L.LT.J) THEN IF (DATA(INDEX(L)).LT.DATA(INDEX(L+1))) L=L+1 END IF IF (T.LT.DATA(INDEX(L))) THEN INDEX(K)=INDEX(L) K=L L=L+L ELSE L=J+1 END IF GO TO 2 END IF INDEX(K)=INDEXT GO TO 1 END C----------------------------------------------------------------------- SUBROUTINE SPLINE (N,X,Y,YPP) C C NATURAL CUBIC SPLINE FIT TO A TABULATED FUNCTION Y(I) = Y(X(I)), I = C 1, 2,..., N, WITH X(1) .LT. X(2) .LT. ... .LT. X(N). RETURNS THE C TABULATED SECOND DERIVATIVES y'' = d2y/dx2 = YPP(I), I = 1, 2,..., N, C WITH YPP(1) = YPP(N) = 0. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C ART OF SCIENTIFIC COMPUTING, PP. 86-89. CAMBRIDGE, ENGLAND: C CAMBRIDGE UNIVERSITY PRESS. C C MUST HAVE NMAX .GE. N. C PARAMETER (NMAX=1000) DIMENSION X(N),Y(N),YPP(N),U(NMAX) YPP(1)=0 YPP(N)=0 U(1)=0 U(N)=0 DO I=2,N-1 P=(X(I)-X(I-1))/(X(I+1)-X(I-1)) Q=P*YPP(I-1)+2 YPP(I)=(P-1)/Q U(I)=(6*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))/ & (X(I)-X(I-1)))/(X(I+1)-X(I-1))-P*U(I-1))/Q END DO DO I=N-1,1,-1 YPP(I)=YPP(I)*YPP(I+1)+U(I) END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE SPLINT (N,X,Y,YPP,XI,YI) C C CUBIC SPLINE INTERPOLATION OF A TABULATED FUNCTION Y(I) = Y(X(I)), I = C 1, 2,..., N, USING THE TABULATED SECOND DERIVATIVES YPP(I) FROM C SUBROUTINE SPLINE. C C FORTRAN CODE QUOTED FROM WILLIAM H. PRESS, BRIAN P. FLANNERY, SAUL A. C TEUKOLSKY, AND WILLIAM T. VETTERLING (1986). NUMERICAL RECIPIES: THE C ART OF SCIENTIFIC COMPUTING, PP. 86-89. CAMBRIDGE, ENGLAND: C CAMBRIDGE UNIVERSITY PRESS. C DIMENSION X(N),Y(N),YPP(N) IF (XI.LE.X(1)) THEN YI=Y(1) RETURN END IF IF (XI.GE.X(N)) THEN YI=Y(N) RETURN END IF I1=1 I2=N 1 IF (I2-I1.GT.1) THEN I=(I2+I1)/2 IF (X(I).GT.XI) THEN I2=I ELSE I1=I END IF GO TO 1 END IF H=X(I2)-X(I1) A=(X(I2)-XI)/H B=(XI-X(I1))/H YI=A*Y(I1)+B*Y(I2)+((A**3-A)*YPP(I1)+(B**3-B)*YPP(I2))*H**2/6 RETURN END C----------------------------------------------------------------------- SUBROUTINE STORE (NCYCLE,R,Z,NFREE,NPAR,SCALEK,B,C,COV,T) DIMENSION B(3,3),C(3,3),COV(13,13),T(109) T(1)=NCYCLE T(2)=R T(3)=Z T(4)=NFREE T(5)=NPAR T(6)=SCALEK I=6 DO J=1,3 DO K=J,3 I=I+1 T(I)=B(J,K) END DO END DO DO J=1,3 DO K=J,3 I=I+1 T(I)=C(J,K) END DO END DO DO J=1,13 DO K=J,13 I=I+1 T(I)=COV(J,K) END DO END DO RETURN END C----------------------------------------------------------------------- FUNCTION SINTHL(IH,IK,IL,GINV) DIMENSION H(3),GINV(3,3) H(1)=IH H(2)=IK H(3)=IL SINTHL=0.5*SQRT(VTMV(H,GINV)) RETURN END C----------------------------------------------------------------------- FUNCTION VTMV(A,B) DIMENSION A(3),B(3,3) VTMV=0 DO I=1,3 DO J=1,3 VTMV=VTMV+A(J)*B(I,J)*A(I) END DO END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE WILSON (IFILE,NDATA,SCALEK,BISO,CISO,COV) C C LOGARITHMICALLY LINEARIZED, WILSON-PLOT FIT TO UNAVERAGED INTENSITIES C TO MINIMIZE C C FSQ(MEAS) C CHISQ = SUM W*[------------------- - Q**2]**2 , C EPSILON*SUM F(J)**2 C C WHERE Q = Q0*EXP(Q1*S**2 + Q2*S**4) , C Q0 = SCALEK**2 , C Q1 = -2* , C Q2 = 2* = 2*<(BISO - )**2> , C AND S = SIN(THETA)/LAMBDA . C C THE LOGARITHMIC FIT MINIMIZES C C FSQ(MEAS) C CHISQ' = SUM W'*{LN ------------------- - [2*LN(SCALEK) C EPSILON*SUM F(J)**2 C - 2*BISO*S**2 C + 2*CISO**2*S**4]}. C C INTENSITIES ARE NOT AVERAGED IN S-SHELLS BEFORE TAKING LOGARITHMS. C LARGE INTENSITIES ARE EFFECTIVELY DOWN-WEIGHTED ON THE LOGARITHMIC C SCALE, HENCE SCALEK IS BIASED TO A VALUE THAT IS TOO LARGE. AND, C SINCE THE LARGE INTENSITIES ARE AT LOW S-VALUES, BISO IS BIASED TO C A VALUE THAT IS TOO SMALL AND CISO IS BIASED TO A VALUE THAT IS TOO C LARGE. C DIMENSION COV(13,13),A(3,3),AINV(3,3),B(3),C(3) REAL*8 SUMW,SUMX,SUMY,SUMX2,SUMX3,SUMX4,SUMY2,SUMXY,SUMX2Y N=0 SUMW=0 SUMX=0 SUMY=0 SUMX2=0 SUMX3=0 SUMX4=0 SUMY2=0 SUMXY=0 SUMX2Y=0 REWIND IFILE DO I=1,NDATA READ (IFILE) J,K,L,SSQ,YOBS,SIGOBS,MULT,VAR IF (YOBS.LT.SIGOBS) THEN C C THRESHOLD INTENSITY FSQMIN/2 FOR (VAR = 1) ACENTRIC REFLECTIONS, AND C FSQMIN/3 FOR (VAR = 2) CENTRIC REFLECTIONS. C IF (VAR.EQ.1) YOBS=SIGOBS/2 IF (VAR.EQ.2) YOBS=SIGOBS/3 END IF C C MINIMIZE CHISQ = SUM W*DELTA**2 = SUM (DELTA/SIGMA(DELTA))**2. C C ON THE DATA SCALE: C C DELTA = YOBS - YCALC C = YOBS - (SCALEK**-2)*EXP{-2*[BISO - (CISO*S)**2]*S**2} C C SIGMA = SQRT[SIGMA(YOBS)**2 + SIGMA(YCALC)**2] C = SQRT[SIGMA(YOBS)**2 + YCALC**2*VAR] C C ON THE LOGARITHMIC FITTING SCALE: C C DELTA = LN(YOBS) - LN(YCALC) C C SIGMA = SQRT{[SIGMA(YOBS)/YOBS]**2 + [SIGMA(YCALC)/YCALC]**2} C = SQRT{[SIGMA(YOBS)/YOBS]**2 + VAR} C X=SSQ Y=LOG(YOBS) W=MULT/((SIGOBS/YOBS)**2+VAR) C C ACCUMULATE LEAST-SQUARES SUMS. C N=N+MULT SUMW=SUMW+W SUMX=SUMX+W*X SUMY=SUMY+W*Y SUMX2=SUMX2+W*X**2 SUMX3=SUMX3+W*X**3 SUMX4=SUMX4+W*X**4 SUMY2=SUMY2+W*Y**2 SUMXY=SUMXY+W*X*Y SUMX2Y=SUMX2Y+W*X**2*Y END DO C C STORE AND THEN INVERT THE NORMAL MATRIX C A(1,1)=SUMW A(1,2)=SUMX A(1,3)=SUMX2 A(2,1)=A(1,2) A(2,2)=SUMX2 A(2,3)=SUMX3 A(3,1)=A(1,3) A(3,2)=A(2,3) A(3,3)=SUMX4 CALL MINV3 (A,AINV,DET) C C STORE THE NORMAL VECTOR AND THEN EVALUATE THE PARAMETERS AND THE C VAR-COV MATRIX. C B(1)=SUMY B(2)=SUMXY B(3)=SUMX2Y CALL MV (AINV,B,C) Q0=C(1) Q1=C(2) Q2=C(3) C C CHISQ = SUM W*(Y - Q)**2 C = SUM W*Y**2 + SUM W*Q**2 - 2*SUM W*Q*Y C C Q = Q0 + Q1*X + Q2*X**2 C C Q**2 = Q0**2 + (Q1*X + Q2*X**2)**2 + 2*Q0*(Q1*X + Q2*X**2) C = Q0**2 + Q1**2*X**2 + Q2**2*X**4 + 2*Q1*Q2*X**3 C + 2*Q0*Q1*X + 2*Q0*Q2*X**2 C CHISQ=SUMY2+Q0**2*SUMW+Q1**2*SUMX2+Q2**2*SUMX4 & +2*(Q0*Q1*SUMX+Q0*Q2*SUMX2+Q1*Q2*SUMX3) & -2*(Q0*SUMY+Q1*SUMXY+Q2*SUMX2Y) NPAR=3 NFREE=N-3 ZSQ=CHISQ/NFREE V0=ZSQ*AINV(1,1) V1=ZSQ*AINV(2,2) V2=ZSQ*AINV(3,3) COV01=ZSQ*AINV(1,2) COV02=ZSQ*AINV(1,3) COV12=ZSQ*AINV(2,3) C C IF THE Q2 = 2*CISO**2 COEFFICIENT OF THE TERM QUADRATIC IN S**2 IS NOT C STATISTICALLY SIGNIFICANT, RESORT TO A SCALING FUNCTION LINEAR IN C S**2. C IF (V2.LE.0.OR.Q2.LT.SQRT(ABS(V2))) THEN DET=SUMW*SUMX2-SUMX*SUMX Q0=(SUMX2*SUMY-SUMX*SUMXY)/DET Q1=(-SUMX*SUMY+SUMW*SUMXY)/DET CHISQ=SUMY2+Q0**2*SUMW+Q1**2*SUMX2 & +2*Q0*Q1*SUMX & -2*(Q0*SUMY+Q1*SUMXY) NPAR=NPAR-1 NFREE=NFREE+1 ZSQ=CHISQ/NFREE V0=ZSQ*SUMX2/DET V1=ZSQ*SUMW/DET COV01=ZSQ*(-SUMX/DET) Q2=0 V2=0 COV02=0 COV12=0 C C IF THE Q1 = -2*BISO COEFFICIENT OF THE TERM LINEAR IN S**2 IS NOT C STATISTICALLY SIGNIFICANT, RESORT TO A CONSTANT SCALING FUNCTION. C IF (V1.LE.0.OR.-Q1.LT.SQRT(ABS(V1))) THEN Q0=SUMY/SUMW CHISQ=SUMY2+Q0**2*SUMW-2*Q0*SUMY NPAR=NPAR-1 NFREE=NFREE+1 ZSQ=CHISQ/NFREE V0=ZSQ/SUMW Q1=0 V1=0 COV01=0 END IF END IF C C TRANSFORM PARAMETERS AND VAR-COV MATRIX FROM THE LOGARITHMIC FITTING C SCALE TO THE DATA SCALE, AND CALCULATE AGREEMENT STATISTICS. C SCALEK=EXP(-0.5*Q0) BISO=-0.5*Q1 CISO=SQRT(0.5*Q2) DKD0=-0.5*SCALEK DBD1=-0.5 IF (CISO.GT.0) THEN DCD2=1/CISO ELSE DCD2=0 END IF COV(1,1)=DKD0**2*V0 COV(2,2)=DBD1**2*V1 COV(3,3)=DCD2**2*V2 COV(1,2)=DKD0*DBD1*COV01 COV(1,3)=DKD0*DCD2*COV02 COV(2,3)=DBD1*DCD2*COV12 COV(2,1)=COV(1,2) COV(3,1)=COV(1,3) COV(3,2)=COV(2,3) RETURN END C----------------------------------------------------------------------- SUBROUTINE MINV3 (A,B,D) C C INVERSE OF A 3 X 3 MATRIX C DIMENSION A(3,3),B(3,3) C C DETERMINANT C D=A(1,1)*A(2,2)*A(3,3)+A(1,2)*A(2,3)*A(3,1)+A(1,3)*A(2,1)*A(3,2) & -A(3,1)*A(2,2)*A(1,3)-A(3,2)*A(2,3)*A(1,1)-A(3,3)*A(2,1)*A(1,2) IF (D.EQ.0) RETURN C C MATRIX OF MINORS C B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3) B(1,2)=A(2,1)*A(3,3)-A(3,1)*A(2,3) B(1,3)=A(2,1)*A(3,2)-A(3,1)*A(2,2) B(2,1)=A(1,2)*A(3,3)-A(3,2)*A(1,3) B(2,2)=A(1,1)*A(3,3)-A(3,1)*A(1,3) B(2,3)=A(1,1)*A(3,2)-A(3,1)*A(1,2) B(3,1)=A(1,2)*A(2,3)-A(2,2)*A(1,3) B(3,2)=A(1,1)*A(2,3)-A(2,1)*A(1,3) B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2) C C MATRIX OF CO-FACTORS (SIGNED MINORS) C DO I=1,3 DO J=1,3 B(I,J)=(-1)**(I+J)*B(I,J) END DO END DO C C ADJOINT MATRIX (TRANSPOSED MATRIX OF CO-FACTORS) C DO I=1,3-1 DO J=I+1,3 T=B(I,J) B(I,J)=B(J,I) B(J,I)=T END DO END DO C C INVERSE MATRIX C DO I=1,3 DO J=1,3 B(I,J)=B(I,J)/D END DO END DO RETURN END C----------------------------------------------------------------------- SUBROUTINE MV (A,B,C) C C SQUARE MATRIX-COLUMN VECTOR MULTIPLICATION C C A(I,J)*B(J) = C(I) C DIMENSION A(3,3),B(3),C(3) DO I=1,3 C(I)=0 DO J=1,3 C(I)=C(I)+A(I,J)*B(J) END DO END DO RETURN END C----------------------------------------------------------------------- C===================================================================== C Machine-dependent system-service routines (such as TIME and DATE) C===================================================================== SUBROUTINE VFDATE(DAYTIME) C C This routine returns the date/time as a 20 charater string. C Truncation occurs if the length of the input string is less than 20. C The format of the returned string is "dd mmm yyyy hh:mm:ss". C Note: this routine is Y2K compliant. C C +++++++++++++++++++++++++++++++++++++++ C Machine dependent intel/Win32 version. C +++++++++++++++++++++++++++++++++++++++ C IMPLICIT NONE C input/output CHARACTER*(*) DAYTIME C local CHARACTER DATE*8, TIME*10, ZONE*5, MONTH(12)*3, TEMP*20 INTEGER VALUES(8) DATA MONTH/'Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', & 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'/ C CALL DATE_AND_TIME(DATE, TIME, ZONE, VALUES) TEMP = DATE(7:8) // ' ' // MONTH(VALUES(2)) // ' ' // & DATE(1:4) // ' ' // TIME(1:2) // ':' // TIME(3:4) // & ':' // TIME(5:6) DAYTIME = TEMP RETURN END C===================================================================== SUBROUTINE VDATE(ADATE) C C This routine returns the date as a 9 charater string. Truncation C occurs if the length of the input string is less than 9. The format C of the returned string is "dd-mmm-yy". C C Note: this routine is not Y2K compliant. C C +++++++++++++++++++++++++++++++++++++++ C Machine dependent intel/Win32 version. C +++++++++++++++++++++++++++++++++++++++ C IMPLICIT NONE C input/output CHARACTER*(*) ADATE C local CHARACTER TEMP*20 C CALL VFDATE(TEMP) ADATE = TEMP(1:2) // '-' // TEMP(4:6) // '-' // TEMP(10:11) RETURN END C===================================================================== SUBROUTINE VTIME(ATIME) C C This routine returns the time of day as an 8 character string. C Truncation occurs if the length of the input string is less than 8. C The format of the returned string is "hh:mm:ss". C C +++++++++++++++++++++++++++++++++++++++ C Machine dependent intel/Win32 version. C +++++++++++++++++++++++++++++++++++++++ C IMPLICIT NONE C input/output CHARACTER*(*) ATIME C local CHARACTER TEMP*20 C CALL VFDATE(TEMP) ATIME = TEMP(13:20) RETURN END C===================================================================== C----------------------------------------------------------------------- SUBROUTINE WPLOT (IFILE,ILP,ATIME,ADATE,TITLE,SCALEK,BISO,CISO, & NDATA) C C WILSON-PLOT DISPLAYS OF DATA AND FITTED CURVES. C C X = [SIN(THETA)/LAMBDA]**2 C Y = FSQ(MEAS)/[EPSILON*SUM(J) F(J)**2] C M = RECIPROCAL LATTICE POINT GROUP MULTIPLICITY C CHARACTER ATIME*8,ADATE*9,TITLE*80,FILEO*80 DIMENSION X(100),Y(100),YPP(100) C C FORTRAN-90 DYNAMIC STORAGE ALLOCATION C REAL, ALLOCATABLE::DATA(:) INTEGER, ALLOCATABLE::INDEX(:) ALLOCATE (DATA(NDATA),INDEX(NDATA)) C C SORT ON [SIN(THETA)/LAMBDA]**2. C DO I=1,NDATA READ (IFILE,REC=I) XI,YI,M DATA(I)=XI END DO CALL SORT (NDATA,DATA,INDEX) C C EVALUATE LOCAL AVERAGE INTENSITIES IN SHELLS OF SIN(THETA)/LAMBDA. C C OVERLAPPING MOVING WINDOW AVERAGES C NSHELL=100 NMOVE=NDATA/NSHELL NMOVE=MAX(NMOVE,1) NLOCAL=100 NLOCAL=MIN(NLOCAL,NDATA) NLOCAL=MAX(NLOCAL,NMOVE) DO I=1,NSHELL X(I)=0 Y(I)=0 SUMM=0 SUMX=0 SUMY=0 IREC=(I-1)*NMOVE DO N=1,NLOCAL IREC=IREC+1 IF (IREC.GT.NDATA) GO TO 5 READ (IFILE,REC=INDEX(IREC)) XI,YI,M SUMM=SUMM+M SUMX=SUMX+M*XI SUMY=SUMY+M*YI END DO 5 CONTINUE IF (SUMM.GT.0) THEN X(I)=SUMX/SUMM Y(I)=SUMY/SUMM END IF END DO CLOSE (UNIT=IFILE,STATUS='DELETE') DEALLOCATE (DATA,INDEX) C C ELIMINATE ANY EMPTY SHELLS. C N=0 DO I=1,NSHELL IF (X(I).NE.0.OR.Y(I).NE.0) THEN N=N+1 X(N)=X(I) Y(N)=Y(I) END IF END DO NSHELL=N C C ENSURE INCREASING X VALUES FOR SPLINE INTERPOLATION. C N=1 DO I=2,NSHELL IF (X(I).GT.X(I-1)) THEN N=N+1 X(N)=X(I) Y(N)=Y(I) END IF END DO NSHELL=N CALL SPLINE (NSHELL,X,Y,YPP) WRITE (ILP,6000) ATIME,ADATE,TITLE CALL XYPLOT (NSHELL,X,Y,YPP,SCALEK,BISO,CISO,ILP,0) WRITE (ILP,6001) SCALEK,BISO,CISO WRITE (ILP,6000) ATIME,ADATE,TITLE CALL XYPLOT (NSHELL,X,Y,YPP,SCALEK,BISO,CISO,ILP,1) WRITE (ILP,6002) SCALEK,BISO,CISO 6000 FORMAT (/1H1,130('-')/1X,'PROGRAM LEVY. ',A,', ',A,'. ',A/) 6001 FORMAT (/1X, &'Y = VERSUS X = = <(SI', &'N(THETA)/LAMBDA)**2>'//1X, &'"*" FITTED FUNCTION: Y = K**(-2)*EXP[-2*(BISO - CISO**2*S**2)*', &'S**2], K = ',E10.3,', BISO = ',E10.3,', CISO = ',E10.3/1X, &'"O" SPLINE-INTERPOLATED DATA CURVE') 6002 FORMAT (/1X, &'LOG(Y) = LOG() VERSUS X = = <(SIN(THETA)/LAMBDA)**2>'//1X, &'"*" FITTED FUNCTION Y = K**(-2)*EXP[-2*(BISO - CISO**2S**2)**2', &'] K = ',E10.3,' BISO = ',E10.3,' CISO = ',E10.3/1X, &'"O" SPLINE-INTERPOLATED DATA CURVE') RETURN END C----------------------------------------------------------------------- SUBROUTINE XYPLOT (N,X,Y,YPP,A0,A1,A2,ILP,ILOG) PARAMETER (NX=100,NY=50) CHARACTER AA*1 DIMENSION AA(0:NX,0:NY),X(N),Y(N),YPP(N) C C BLANK OUT THE PLOT ARRAY. C DO IX=0,NX DO IY=0,NY AA(IX,IY)=' ' END DO END DO C C FILL IN THE GRID MARKS. C DO IX=0,NX AA(IX, 0)='.' AA(IX,NY)='.' END DO DO IX=0,NX,10 AA(IX, 0)='+' AA(IX,NY)='+' END DO DO IY=0,NY AA( 0,IY)='.' AA(NX,IY)='.' END DO DO IY=0,NY,10 AA( 0,IY)='+' AA(NX,IY)='+' END DO C C SCALE PLOT AXES. C XMIN=0 XMAX=X(N) YMIN=+1E9 YMAX=-1E9 DO I=1,N YMIN=MIN(YMIN,Y(I)) YMAX=MAX(YMAX,Y(I)) END DO IF (A0.GT.0) YMAX=MAX(YMAX,A0**(-2)) IF (ILOG.EQ.1) THEN YMIN=LOG(YMIN) YMAX=LOG(YMAX) ELSE YMIN=0 END IF C C FILL IN FITTED FUNCTION. C NI=4*NX DX=(XMAX-XMIN)/NI DO I=0,NI XI=XMIN+I*DX IX=NINT(((XI-XMIN)/(XMAX-XMIN))*NX) IF (A0.GT.0.AND.A1.GE.0.AND.A2.GE.0) THEN YI=A0**(-2)*EXP(-2*A1*XI+2*(A2*XI)**2) IF (ILOG.EQ.1) YI=LOG(YI) IY=NINT(((YMAX-YI)/(YMAX-YMIN))*NY) IF (IX.GE.0.AND.IX.LE.NX.AND.IY.GE.0.AND.IY.LE.NY) & AA(IX,IY)='*' END IF END DO C C FILL IN SPLINE-INTERPOLATED DATA CURVE. C DX=(X(N)-X(1))/NI DO I=0,NI XI=X(1)+I*DX IX=NINT(((XI-XMIN)/(XMAX-XMIN))*NX) CALL SPLINT (N,X,Y,YPP,XI,YI) IF (ILOG.EQ.1) THEN IF (YI.LE.0) THEN YI=YMIN ELSE YI=LOG(YI) END IF END IF IY=NINT(((YMAX-YI)/(YMAX-YMIN))*NY) IF (IX.GE.0.AND.IX.LE.NX.AND.IY.GE.0.AND.IY.LE.NY) & AA(IX,IY)='O' END DO C C PRINT THE PLOT ARRAY. C DY=(YMAX-YMIN)/NY DO IY=0,NY IF (MOD(IY,10).EQ.0) THEN WRITE (ILP,100) YMAX-IY*DY,(AA(IX,IY),IX=0,NX) ELSE WRITE (ILP,101) (AA(IX,IY),IX=0,NX) END IF END DO DX=(XMAX-XMIN)/10 WRITE (ILP,102) (XMIN+IX*DX,IX=0,10) 100 FORMAT (1X,E10.3,101A1) 101 FORMAT (1X,10X, 101A1) 102 FORMAT (1X,1X,11F10.4) RETURN END C-----------------------------------------------------------------------