program locscl c c Robert H. Blessing c Hauptman-Woodward Institute c 73 High St. c Buffalo, NY 14203, USA c Tel. 716-856-9600, ext. 335 c blessing@hwi.buffalo.edu c c Dec. 1995,..., Jan. 2001 c c Scaling by least-squares fitted local scale factors c c q = c c to improve the accuracy of structure factor magnitude differences c c Delta(F) = F1 - q*F2 c c for SIR F1(native)/F2(derivative) or SAS F1(+h+k+l)/F2(-h-k-l) pairs. c c The local scale factors are evaluated in approximately equidimensional c local blocks of c c n = {[2*m(h) + 1] X [2*m(k) + 1] X [2*m(l) + 1]} - 1 c c reciprocal lattice points surrounding, but excluding, each point hkl. c c B.W. Matthews and E.W. Czerwinski (1975). Local Scaling: a Method to c Reduce Systematic Errors in Isomorphous Replacement and Anomalous c Scattering Measurements. Acta Cryst. A31, 480-497. c c W.A. Hendrickson (1988). Structural Analysis by the Method of c Multiwavelength Anomalous Diffraction. In "Crystallographic c Computing 4", edited by N.W. Isaacs and M.R. Taylor, pp. 97-108. c Oxford: Oxford Univ. Press. And references cited therein. c character atime*8,adate*9,title*72,type*3,ptgp*5, & file1*40,file2*40,form1*11,form2*11 c c Fortran-90 dynamic storage allocation c real, allocatable::data(:,:) dimension cell(6),ginv(3,3) data ptgp,itype1,itype2,cutoff,cell,mh,mk,ml,i012,zlimit & /' ',14*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 file "locscl.dat". c open (unit=io1,file='locscl.dat',status='old') read (io1,'(a)') title read (io1,'(a)') type if (type.eq.'sir') type='SIR' if (type.eq.'sas') type='SAS' read (io1,'(a)') file1 if (type.eq.'SIR') read (io1,'(a)') file2 read (io1,*,end=1) itype1 if (type.eq.'SIR') read (io1,*,end=1) itype2 read (io1,*) cutoff read (io1,'(a)') ptgp read (io1,*) (cell(i),i=1,6) read (io1,*,end=1) mh,mk,ml read (io1,*,end=1) i012 read (io1,*,end=1) zlimit 1 close (unit=io1,status='keep') c c Maximum local block dimensions are c (2*mh + 1)*(2*mk + 1)*(2*ml + 1) = 25*25*25. c if (max(mh,mk,ml).gt.12) then write (6,'(/1x, &''locscl/main: Choose smaller mh, mk, and/or ml .le. 12.'')') stop end if c----------------------------------------------------------------------- c c Control data: c ------------ c c title = Job title/crystal name c c type = SAS or SIR = Type of data set pair c c file1 = Name of SAS reflections file OR SIR Native reflections file c (file2 = Name of SIR Derivative reflections file) c c itype1 = File type (see below) of file1 c (itype2 = File type of file2) c c Input reflection data files: c --------------------------- c c For SAS data, one file: c ------------------ c file1 containing both F(+h+k+l) and F(-h-k-l) in separate reflection c records, in any order. c c For SIR data, two files: c ------------------- c file1 for the native data, and c file2 for the derivative data. In both files the reflection records c may be in any order, and the two files need not be of the same type. c c c The files are read in subroutine "read1" according to the values of c itype1 and itype2: c ------ ------ c c itype = 0 ASCII, formatted h,k,l,Fsq,sigFsq,F,sigF,E,sigE c 1 " " h,k,l,Fsq,sigFsq,F,sigF c 2 " " h,k,l,Fsq,sigFsq c 3 " " h,k,l, F,sigF,E,sigE c 4 " " h,k,l, F,sigF c 5 Binary, unformatted h,k,l,Fsq,sigFsq,F,sigF,E,sigE c 6 " " h,k,l,Fsq,sigFsq,F,sigF c 7 " " h,k,l,Fsq,sigFsq c 8 " " h,k,l, F,sigF,E,sigE c 9 " " h,k,l, F,sigF c c c cutoff = minimum F/sigma(F) threshold for input data to be included in c the scaling calculations and the locally scaled output files. c c Default is cutoff = 3. To apply cutoff = 0, enter c cutoff .lt. 0. c c ptgp = Point group symbol c c Must be one of the following case-sensitive character values. c c Tricl. Monocl. Ortho. Tetrag. Trig. Rhomb. Cub. c ------ ------- ------ ------- ----- ------ ---- c 1 2 222 4 3 R3 23 c -1 m mm2 -4 -3 R-3 m3 c 2/m mmm 4/m c 312 R32 432 c 422 321 R3m -43m c 4mm 31m R-3m m3m c -42m 3m1 c -4m2 -31m c 4/mmm -3m1 c c Hexag. c ------ c 6 c -6 c 6/m c c 622 c 6mm c -6m2 c -62m c 6/mmm c c cell = lattice parameters (a, b, c, alpha, beta, gamma) c c mh,mk,ml = semidimensions m(h), m(k), m(l) if fixed semidimensions c are desired for the moving local block of c c n = [2*m(h) + 1]*[2*m(k) + 1]*[2*m(l) + 1] - 1 c c reciprocal lattice points surrounding, but excluding, each c point hkl. c c The default (mh = mk = ml = 0) is to automatically vary and c statistically optimize the semidimensions. c c To bypass the local scale factor fitting and apply constant c local scale factors of unity, supply mh = mk = ml = -1. c c If fixed semidimensions are chosen, choose c mh, mk, ml .le. 12. c c i012 = 0 Keep constant the mean of F1 and F2. c 1 " " " value of F1. c 2 " " " " " F2. c c The default i012 = 0 is recommended. c c zlimit = limiting value of c c abs(Delta) abs(Delta) c ---------- = ----------------------- c sigma 1.25*median[abs(Delta)] c c abs[F1 - F2 - median(F1 - F2)] c = ------------------------------------------- , c 1.25*median{abs[F1 - F2 - median(F1 - F2)]} c c for reflection pairs to be accepted for the output file. c The variable zlimit is used to exclude data with unreliably c large values of abs(F1 - F2) in the tails of the data-set c Delta(F) distribution. c c The zlimit test assumes that the data-set c z = Delta/sigma distribution should approximate a c zero-mean, unit-variance normal distribution, for which c values of z .lt. -zlimit or z .gt. +zlimit are c extremely improbable. c c With the default zlimit = 0 , the zlimit exclusion is not c applied. A recommended trial value is zlimit = 6. c c----------------------------------------------------------------------- open (unit=ilp,file='locscl.lp',status='new') write (ilp,'(/1x,''Program locscl. '',a,'', '',a//1x, &''Least-squares local scale factors''/1x, &'' q = ''/1x, &''to improve the accuracy of structure factor magnitude differe'', &''nces''/1x, &'' Delta(F) = F1 - q*F2''/1x, &''for SIR [F1(native)/F2(derivative)] or SAS [F1(+h+k+l)/F2(-h-'', &''k-l)] pairs.''//1x, &''The local scale factors are evaluated in approximately equidi'', &''mensional''/1x,''local blocks of''/1x, &'' n = {[2*delta(h) + 1] X [2*delta(k) + 1] X [2*delta(l) + 1]'', &''} - 1''/1x, &''reciprocal lattice points surrounding, but excluding, each po'', &''int hkl.'')') atime,adate write (ilp,'(/1x,''Job title:''/1x,a)') title write (ilp,'(/1x,''Lattice parameters:''/1x, &'' a b c alpha beta gamma''/ &1x,3f10.4,3f10.3)') cell call metric (cell,ginv) if (cutoff.eq.0) cutoff=3 if (cutoff.lt.0) cutoff=0 c c For SIR data, open native and derivative reflection data files. c if (type.eq.'SIR') then write (ilp,'(/1x, & ''SIR File-1 [F1(Native)]:''/1x,a//1x, & ''SIR File-2 [F2(Derivative)]:''/1x,a//1x, & ''File-1 type: itype1 = '',i2/1x, & ''File-2 type: itype2 = '',i2)') file1,file2,itype1,itype2 form1='formatted' form2='formatted' if (itype1.ge.5) form1='unformatted' if (itype2.ge.5) form2='unformatted' open (unit=io1,file=file1,form=form1,status='old') open (unit=io2,file=file2,form=form2,status='old') c c For SAS data, prepare separate +h+k+l and -h-k-l reflection data files c from one input file. c else if (type.eq.'SAS') then write (ilp,'(/1x, & ''SAS Input File [F1(+h+k+l) and F2(-h-k-l)]:''/1x,a//1x, & ''File type: itype1 = '',i2)') file1,itype1 form1='formatted' if (itype1.ge.5) form1='unformatted' open (unit=io3,file=file1,status='old',form=form1) open (unit=io2,file='data.neg',status='new',form='unformatted') open (unit=io1,file='data.pos',status='new',form='unformatted') call bijvoet (io1,io2,io3,ilp,itype1,ptgp) close (unit=io3,status='keep') rewind io2 rewind io1 itype1=8 itype2=8 end if c c Count reflection data pairs for scaling. c call count (io1,io2,ilp,itype1,itype2,cutoff,m1,n1,m2,n2, & ihmin,ikmin,ilmin,nh,nk,nl,imax) c c Allocate data storage. c allocate (data(10,0:imax)) c c Store reflection data pairs in memory. c call datain (io1,io2,ilp,itype1,itype2,cutoff,m1,n1,m2,n2, & ihmin,ikmin,ilmin,nh,nk,nl,ginv,imax,data,jmin,jmax,type,ptgp) if (type.eq.'SIR') then close (unit=io1,status='keep') close (unit=io2,status='keep') else if (type.eq.'SAS') then close (unit=io1,status='delete') close (unit=io2,status='delete') end if c c Global pre-scaling of SIR pairs c if (type.eq.'SIR') call pscale (imax,data,jmin,jmax,cutoff,ginv, & ptgp,ihmin,ikmin,ilmin,nk,nl,ilp) c c Local scaling c write (ilp,'(//1x,''Local Scaling''/1x,''-------------'')') if (mh.lt.0.or.mk.lt.0.or.ml.lt.0) then c c Apply fixed local scale factors q = 1. c write (ilp,'(/1x,''All local scale factors fixed at q = 1.'')') go to 15 else if (mh.eq.0.and.mk.eq.0.and.ml.eq.0) then c c Optimize the local block semidimensions mh, mk, and ml to maximize c delta = sqrt(rmsd**2 - esd**2) - esd , c where c rmsd = <(q - )**2>**(1/2) , c esd = **(1/2) , c and c sigma(q) = sqrt{sum w*[(F1/F2) - q]**2/[(n - 1)*sum w]} . c c Both rmsd and esd tend to decrease with increasing local block size. c Rmsd**2 estimates the total, systematic-plus-random, variance. c Esd is the rms error of fit, and esd**2 estimates the random variance. c Thus delta estimates the difference between the systematic and random c local variation. c write (ilp,'(1x, &''Optimize the local block semidimensions mh, mk, and ml to max'', &''imize''/1x, &'' delta = sqrt(rmsd**2 - esd**2) - esd ,''/1x,''where''/1x, &'' rmsd = <(q - )**2>**(1/2) ,''/1x, &'' esd = **(1/2) ,''/1x,''and''/1x, &'' sigma(q) = sqrt{sum w*[(F1/F2) - q]**2/[(n - 1)*sum w]} .'')') write (ilp,'(/1x, &''Both rmsd and esd tend to decrease with increasing local bloc'', &''k size.''/1x, &''Rmsd**2 estimates the total, systematic-plus-random, variance.'' &/1x, &''Esd is the rms error of fit, and esd**2 estimates the random '', &''variance.''/1x, &''Thus delta estimates the difference between the systematic an'', &''d random''/1x,''local variation.'')') write (ilp,'(/1x, &'' m mh mk ml nh nk nl n rmsd esd delta delta(m'', &'')-''/1x, &'' % % % delta('', &''m-1)''/1x, &'' - -- -- -- -- -- -- - --- ---- --- ----- -------'', &''----'')') a=cell(1) b=cell(2) c=cell(3) amax=max(a,b,c) a=a/amax b=b/amax c=c/amax delta2=0 delta1=0 delta=0 m=0 do i=1,12 delta2=delta1 delta1=delta m=m+1 mh=max(nint(m*a),1) mk=max(nint(m*b),1) ml=max(nint(m*c),1) call qfit (imax,data,jmin,jmax,ihmin,ikmin,ilmin,nk,nl, & mh,mk,ml,nrej,ndata,qmin,sigmin,qmax,sigmax,qmean,rmsd,esd, & wqmean,wrmsd,wesd) delta=(sqrt(max((wrmsd**2-wesd**2),0.0))-wesd) write (ilp,'(1x,i2,1x,3i3,1x,3i3,i6,f7.4,3f6.2,e12.3)') & m,mh,mk,ml, & (2*mh+1),(2*mk+1),(2*ml+1),((2*mh+1)*(2*mk+1)*(2*ml+1))-1, & wqmean,100*wrmsd/wqmean,100*wesd/wqmean,100*delta/wqmean, & delta-delta1 if (m.gt.2.and.delta1.gt.max(delta2,delta)) then m=m-1 mh=max(nint(m*a),1) mk=max(nint(m*b),1) ml=max(nint(m*c),1) go to 5 end if end do go to 10 else if (mh.gt.0.or.mk.gt.0.or.ml.gt.0) then c c Use fixed, input local block semidimensions. c write (ilp,'(1x, & ''Fixed, input local block semidimensions:''//1x, & '' mh mk ml nh nk nl n''/1x, & '' -- -- -- -- -- -- -''/1x,2x,3i3,1x,3i3,i6)') & mh,mk,ml, & (2*mh+1),(2*mk+1),(2*ml+1), & ((2*mh+1)*(2*mk+1)*(2*ml+1))-1 go to 5 end if 5 continue c c Use either input or optimized local block dimensions in a final call c to subroutine qfit. c call qfit (imax,data,jmin,jmax,ihmin,ikmin,ilmin,nk,nl,mh,mk,ml, & nrej,ndata,qmin,sigmin,qmax,sigmax,qmean,rmsd,esd,wqmean,wrmsd, & wesd) 10 continue write (ilp,'(/1x, &''Local Scaling Results''/1x,''---------------------''/1x, &''Local blocks of up to''/1x, &'' n = {[2*m(h) + 1] X [2*m(k) + 1] X [2*m(l) + 1]} - 1''/1x, &''reciprocal lattice points surrounding, but excluding, each po'', &''int hkl''/1x, &'' m(h) = '',i3,'' n(h) = '',i3/1x, &'' m(k) = '',i3,'' n(k) = '',i3/1x, &'' m(l) = '',i3,'' n(l) = '',i3,'' n = '',i6)') & mh,(2*mh+1), & mk,(2*mk+1), & ml,(2*ml+1),((2*mh+1)*(2*mk+1)*(2*ml+1))-1 write (ilp,'(/1x, &''Nrej = '',i8,'' reflection pairs rejected due to local scale'', &'' factor''/1x, &'' '',8x,'' error-of-fit uncertainty''//1x, &''N = '',i8,'' locally scaled F1(hkl) and F2(hkl) reflectio'', &''n pairs'')') nrej,ndata write (ilp,'(/1x, &''qmin = '',f8.4,'' sigma(qmin) = '',f8.4/1x, &''qmax = '',f8.4,'' sigma(qmax) = '',f8.4//1x, &'' Unit-Weighted Reciprocal-Variance-'', &''Weighted''/1x, &'' ------------- --------------------'', &''--------''/1x, &''mean = = '',2f8.4/1x, &''rmsd = <(q - )**2>**(1/2) = '',2f8.4/1x, &''esd = **(1/2) = '',2f8.4//1x, &''delta = 100%*[sqrt(wrmsd**2 - wesd**2) - wesd]/wavg = '',f6.2, &''%'')') & qmin,sigmin,qmax,sigmax,qmean,wqmean,rmsd,wrmsd,esd,wesd, & 100*(sqrt(wrmsd**2-wesd**2)-wesd)/wqmean 15 continue c c Apply local scale factors. c call qscale (imax,data,jmin,jmax,ihmin,ikmin,ilmin,nk,nl,io1,io2, & i012,type,ptgp,ndata,nc,ilp) c c Compile locally scaled Delta(F) distribution statistics and histogram, c and write output reflection files. c call output (io1,io2,ilp,type,ptgp,ndata,nc,zlimit,atime,adate, & title) c c Compile subset agreement statistics. c call agree (io1,ilp,ginv) c c Describe the output files. c write (ilp,'(//1x,''Output Files''/1x,''------------'')') if (i012.eq.0) write (ilp,'(1x, &''Output data files have F1 and locally scaled F2 re-scaled to '', &''give''/1x,''the same mean as the unscaled F1 and F2.'')') if (i012.eq.1) write (ilp,'(1x,''Output data files have unchang'', &''ed F1 and locally scaled F2.'')') if (i012.eq.2) write (ilp,'(1x,''Output data files have locally'', &'' scaled F1 and unchanged F2.'')') write (ilp,'(/1x, &''The "data.locscl" output file is sorted on the Laue-group uni'', &''que indices.''/1x, &'' -------------''/1x, &''The file contains header records that label the data records:''/ &1x, &'' h,k,l,F1,sigma(F1),F2,sigma(F2),E1,sigma(E1),E2,sigma(E2),''/ &1x,''where''/1x, &'' 1 designates the Native or the +h+k+l data set, and''/1x, &'' 2 " " Derivative " " -h-k-l " " .'')') write (ilp,'(/1x, &''The "sort.locscl" output file is sorted on''/1x, &'' -------------''//1x, &'' abs[Delta(E)] abs(E1 - E2)''/1x, &'' abs[z(E)] = ------------- = -------------------------------'', &''--''/1x, &'' sigma(Delta) sqrt[sigma(E1)**2 + sigma(E2)**'', &''2]''//1x, &''or, in the absence of E values, the analogous z(F).''/1x, &''The file contains header records that label the data records:''/ &1x, &'' h,k,l,F1,sigma(F1),F2,sigma(F2),E1,sigma(E1),E2,sigma(E2),z'', &''(E or F).'')') if (type.eq.'SAS') write (ilp,'(/1x, &''The "hphn.locscl" output file has adjacent separate records f'', &''or the SAS''/1x, &'' -------------''/1x,''Bijvoet pairs.''/1x, &''The file contains header records that label the data records:''/ &1x, &'' +h,+k,+l,F1,sigma(F1),E1,sigma(E1);''/1x, &'' -h,-k,-l,F2,sigma(F2),E2,sigma(E2).'')') write (ilp,'(/)') stop ' Program locscl finis' end c----------------------------------------------------------------------- subroutine agree (io1,ilp,ginv) dimension ginv(3,3) dimension d(0:20),n(20),sumdf(20),sumf1(20),sumf2(20),sumde(20), & sumae(20),zmin(10) real*8 sumdf,sumf1,sumf2,sumde,sumae write (ilp,'(//1x, &''Agreement Statistics''/1x,''--------------------''/1x, &''Q(F) = /''/1x, &''R(F) = /<0.5*(F1 + F2)>''/1x, &''R(E) = /<0.5*(E1 + E2)>''//1x, &''Statistics in equal-volume resolution shells with''/1x, &''dmax .gt. d .ge. dmin:''//1x, &'' dmax dmin N Q(F) R(F) R(E)''/1x, &'' ---- ---- - ---- ---- ----'')') open (unit=io1,file='data.locscl',status='old') c c Read dummy header records. c do i=1,4 read (io1,'(a1)') dummy end do dlimit=1e9 ndata=0 do i=1,1e6 read (io1,*,end=1) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 dlimit=min(dlimit,1/(2*sinthl(ih,ik,il,ginv))) ndata=ndata+1 end do 1 continue do i=20,1,-1 d(i)=1/((float(i)/20)**(1/3.0)*(1/dlimit)) end do do i=1,20 n(i)=0 sumdf(i)=0 sumf1(i)=0 sumf2(i)=0 sumde(i)=0 sumae(i)=0 end do rewind io1 do i=1,4 read (io1,'(a1)') dummy end do do j=1,ndata read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 di=1/(2*sinthl(ih,ik,il,ginv)) do i=1,20 if (di.ge.d(i)) then n(i)=n(i)+1 sumdf(i)=sumdf(i)+abs(f1-f2) sumf1(i)=sumf1(i)+f1 sumf2(i)=sumf2(i)+f2 sumde(i)=sumde(i)+abs(e1-e2) sumae(i)=sumae(i)+e1+e2 go to 2 end if end do 2 continue end do d(0)=999.99 do i=1,20 qf=0 rf=0 re=0 if (sumf2(i).ne.0) qf=sumf1(i)/sumf2(i) if (sumf1(i)+sumf2(i).ne.0) & rf=sumdf(i)/(0.5*(sumf1(i)+sumf2(i))) if (sumae(i).ne.0) re=sumde(i)/(0.5*sumae(i)) write (ilp,'(1x,2f7.2,i7,3f7.3)') d(i-1),d(i),n(i),qf,rf,re end do write (ilp,'(/1x, &''Statistics for subsets with Y .ge. Ymin, where''/1x, &''Y(E) = min{[E1/sigma(E1)], [E2/sigma(E2)]},''/1x, &''or, in the absence of E-values, the analogous Y(F):''//1x, &'' Ymin N Q(F) R(F) R(E)''/1x, &'' ---- - ---- ---- ----'')') data zmin /0, 1, 2, 3, 4, 6, 8, 10, 30, 100/ do i=1,10 n(i)=0 sumdf(i)=0 sumf1(i)=0 sumf2(i)=0 sumde(i)=0 sumae(i)=0 end do rewind io1 do i=1,4 read (io1,'(a1)') dummy end do read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 if (e1.gt.0.and.e2.gt.0) then ie=1 else ie=0 end if rewind io1 do i=1,4 read (io1,'(a1)') dummy end do do j=1,ndata read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 if (ie.eq.1) then zi=min(e1/t1,e2/t2) else zi=min(f1/s1,f2/s2) end if do i=1,10 if (zi.ge.zmin(i)) then n(i)=n(i)+1 sumdf(i)=sumdf(i)+abs(f1-f2) sumf1(i)=sumf1(i)+f1 sumf2(i)=sumf2(i)+f2 sumde(i)=sumde(i)+abs(e1-e2) sumae(i)=sumae(i)+e1+e2 end if end do end do do i=1,10 qf=0 rf=0 re=0 if (sumf2(i).ne.0) qf=sumf1(i)/sumf2(i) if (sumf1(i)+sumf2(i).ne.0) & rf=sumdf(i)/(0.5*(sumf1(i)+sumf2(i))) if (sumae(i).ne.0) re=sumde(i)/(0.5*sumae(i)) write (ilp,'(1x,2i7,3f7.3)') int(zmin(i)),n(i),qf,rf,re end do write (ilp,'(/1x, &''Statistics for subsets with abs(Z) .ge. Zmin, where''/1x, &''Z(E) = (E1 - E2)/sqrt[sigma(E1)**2 + sigma(E2)**2],''/1x, &''or, in the absence of E-values, the analogous Z(F):''//1x, &'' Zmin N Q(F) R(F) R(E)''/1x, &'' ---- - ---- ---- ----'')') do i=1,10 n(i)=0 sumdf(i)=0 sumf1(i)=0 sumf2(i)=0 sumde(i)=0 sumae(i)=0 end do rewind io1 do i=1,4 read (io1,'(a1)') dummy end do do j=1,ndata read (io1,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 if (ie.eq.1) then zi=(e1-e2)/sqrt(t1**2+t2**2) else zi=(f1-f2)/sqrt(s1**2+s2**2) end if do i=1,10 if (abs(zi).ge.zmin(i)) then n(i)=n(i)+1 sumdf(i)=sumdf(i)+abs(f1-f2) sumf1(i)=sumf1(i)+f1 sumf2(i)=sumf2(i)+f2 sumde(i)=sumde(i)+abs(e1-e2) sumae(i)=sumae(i)+e1+e2 end if end do end do do i=1,10 qf=0 rf=0 re=0 if (sumf2(i).ne.0) qf=sumf1(i)/sumf2(i) if (sumf1(i)+sumf2(i).ne.0) & rf=sumdf(i)/(0.5*(sumf1(i)+sumf2(i))) if (sumae(i).ne.0) re=sumde(i)/(0.5*sumae(i)) write (ilp,'(1x,2i7,3f7.3)') int(zmin(i)),n(i),qf,rf,re end do close (unit=io1,status='keep') return end c----------------------------------------------------------------------- subroutine bijvoet (io1,io2,io3,ilp,itype,ptgp) c c Prepare separate files of +h+k+l reflection data and Bijvoet-related c -h-k-l anti-reflection data from a file of (acentric) point-group c unique data. c character ptgp*5,laue*5 integer h,hp,hn call decode (ptgp,laue,iptgp,ilaue) write (ilp,'(/1x,''Point group: '',a/1x,''Laue group: '',a)') & ptgp,laue nc=0 n1=0 n2=0 iend=0 do i=1,1e9 call read1 (itype,io3,iend,h,k,l,f,sigf,e,sige) if (iend.ne.0) go to 9 c c Transform to point-group unique indicies. c call equiv (iptgp,h,k,l) c c Acentric (ic = 0) or centric (ic = 1) reflection? c call csym (ptgp,h,k,l,ic) if (ic.eq.1) then c c Write centric reflections to both the +h+k+l reflection file (1) and c the -h-k-l anti-reflection file (2). c nc=nc+1 write (io1) h,k,l,f,sigf,e,sige n1=n1+1 write (io2) h,k,l,f,sigf,e,sige n2=n2+1 else c c Test acentric reflections by comparing point-group and Laue-group c unique indicies. c lh=h lk=k ll=l call equiv (ilaue,lh,lk,ll) hp=+lh kp=+lk lp=+ll hn=-lh kn=-lk ln=-ll call equiv (iptgp,hp,kp,lp) call equiv (iptgp,hn,kn,ln) if (h.eq.hp.and.k.eq.kp.and.l.eq.lp) then c c Write to +h+k+l reflection file. c write (io1) hp,kp,lp,f,sigf,e,sige n1=n1+1 else if (h.eq.hn.and.k.eq.kn.and.l.eq.ln) then c c Write to -h-k-l anti-reflection file, but with indices +h+k+l. c write (io2) hp,kp,lp,f,sigf,e,sige n2=n2+1 end if end if end do 9 continue write (ilp,'(/1x, &''n1 = '',i6,'' F1(+h+k+l) SAS data''/1x, &''n2 = '',i6,'' F2(-h-k-l) SAS data'')') n1,n2 if (nc.gt.0) write (ilp,'(/1x, &''nc = '',i6,'' centric reflections stored in''/1x, &'' '',6x,'' both the F1(+h+k+l) and F2(-h-k-l) subsets.''/1x, &'' '',6x,'' ---- ---'')') nc return end c----------------------------------------------------------------------- subroutine count (iunit1,iunit2,ilp,itype1,itype2,cutoff, & m1,n1,m2,n2,hmin,kmin,lmin,nh,nk,nl,imax) c c Count data and find index limits. c integer h,hmin,hmax,hkl write (ilp,'(//1x, &''Data Count for Local Scaling''/1x, &''----------------------------''/1x, &''Minimum F/sigma(F) cutoff:''/1x, &''cutoff = '',f5.2)') cutoff hmin=+1e6 hmax=-1e6 kmin=+1e6 kmax=-1e6 lmin=+1e6 lmax=-1e6 m1=0 n1=0 iend=0 do i=1,1e6 call read1 (itype1,iunit1,iend,h,k,l,f,sigmaf,e,sigmae) if (iend.ne.0) go to 1 m1=m1+1 if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then n1=n1+1 hmin=min(hmin,h) hmax=max(hmax,h) kmin=min(kmin,k) kmax=max(kmax,k) lmin=min(lmin,l) lmax=max(lmax,l) end if end do 1 continue m2=0 n2=0 iend=0 do i=1,1e6 call read1 (itype2,iunit2,iend,h,k,l,f,sigmaf,e,sigmae) if (iend.ne.0) go to 2 m2=m2+1 if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then n2=n2+1 hmin=min(hmin,h) hmax=max(hmax,h) kmin=min(kmin,k) kmax=max(kmax,k) lmin=min(lmin,l) lmax=max(lmax,l) end if end do 2 continue nl=lmax-lmin+1 nk=kmax-kmin+1 nh=hmax-hmin+1 imax=nh*nk*nl+nk*nl+nl write (ilp,'(/1x, &''m1 = '',i6,'' total F1 data''/1x, &''m2 = '',i6,'' total F2 data''//1x, &''n1 = '',i6,'' F1 data with F1 .ge. cutoff*sigma(F1)''/1x, &''n2 = '',i6,'' F2 data with F2 .ge. cutoff*sigma(F2)'')') & m1,m2,n1,n2 return end c----------------------------------------------------------------------- subroutine csym (g,h,k,l,c) c c Table of centrosymmetric rows and zones in non-centrosymmetric c reciprocal lattice point groups c c c = 0 for acentric distributions c c = 1 for centric distributions 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 character g*5 integer h,c c=0 c c Triclinic c if (g.eq.'1 ') then c=0 else if (g.eq.'-1 ') then c=1 c c Monoclinic (b-axis unique) c else if (g.eq.'2 ') then if (k.eq.0) c=1 else if (g.eq.'m ') then if (h.eq.0.and.l.eq.0) c=1 else if (g.eq.'2/m ') then c=1 c c Orthorhombic c else if (g.eq.'222 ') then if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1 else if (g.eq.'mm2 ') then if (l.eq.0) c=1 else if (g.eq.'mmm ') then c=1 c c Tetragonal c else if (g.eq.'4 ') then if (l.eq.0) c=1 else if (g.eq.'-4 ') then if (l.eq.0) c=1 if (h.eq.0.and.k.eq.0) c=1 else if (g.eq.'4/m ') then c=1 else if (g.eq.'422 ') then if (l.eq.0.or.k.eq.0.or.h.eq.0) c=1 if (abs(k).eq.abs(h)) c=1 else if (g.eq.'4mm ') then if (l.eq.0) c=1 else if (g.eq.'-42m ') then if (l.eq.0.or.k.eq.0.or.h.eq.0) c=1 else if (g.eq.'-4m2 ') then if (l.eq.0) c=1 if (abs(k).eq.abs(h)) c=1 else if (g.eq.'4/mmm') then c=1 c c Trigonal, hexagonal axes c else if (g.eq.'3 ') then c=0 else if (g.eq.'-3 ') then c=1 else if (g.eq.'312 ') then if (k.eq.h.or.k.eq.-2*h) c=1 else if (g.eq.'321 ') then if (k.eq.0.or.h.eq.0) c=1 else if (g.eq.'31m ') then if ((k.eq.0.and.l.eq.0).or.(h.eq.0.and.l.eq.0)) c=1 else if (g.eq.'3m1 ') then c=0 else if (g.eq.'-31m ') then c=1 else if (g.eq.'-3m1 ') then c=1 c c Hexagonal c else if (g.eq.'6 ') then if (l.eq.0) c=1 else if (g.eq.'-6 ') then if (h.eq.0.and.k.eq.0) c=1 else if (g.eq.'6/m ') then c=1 else if (g.eq.'622 ') then 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 else if (g.eq.'6mm ') then if (l.eq.0) c=1 else if (g.eq.'-6m2 ') then if (k.eq.h.or.k.eq.-2*h) c=1 else if (g.eq.'-62m ') then if (h.eq.0.or.k.eq.0) c=1 else if (g.eq.'6/mmm') then c=1 c c Trigonal, rhombohedral axes c else if (g.eq.'R3 ') then c=0 else if (g.eq.'R-3 ') then c=1 else if (g.eq.'R32 ') then if (h.eq.k.or.k.eq.l.or.l.eq.h) c=1 else if (g.eq.'R3m ') then 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 else if (g.eq.'R-3m ') then c=1 c c Cubic c else if (g.eq.'23 ') then if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1 else if (g.eq.'m3 ') then c=1 else if (g.eq.'432 ') then 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 else if (g.eq.'-43m ') then if (h.eq.0.or.k.eq.0.or.l.eq.0) c=1 else if (g.eq.'m3m ') then c=1 else stop ' locscl/csym: Invalid point group symbol.' end if return end c----------------------------------------------------------------------- subroutine datain (iunit1,iunit2,ilp,itype1,itype2,cutoff, & m1,n1,m2,n2,hmin,kmin,lmin,nh,nk,nl,ginv,imax,data,jmin,jmax, & type,ptgp) c c Store reflection data in memory. c integer h,hmin,hmax,hkl character type*3,ptgp*5 dimension data(10,0:imax),ginv(3,3) real*8 sumdf,sumf1,sumf2,sumde,sume,sum1sq,sum2sq write (ilp,'(//1x, &''Data Selection for Local Scaling''/1x, &''--------------------------------''/1x, &''Minimum F/sigma(F) cutoff:''/1x, &''cutoff = '',f5.2)') cutoff c c Zero the data array. c do i=0,imax do j=1,10 data(j,i)=0 end do end do c c Store data by packed index. c jmin=imax jmax=0 rewind iunit1 do i=1,m1 call read1 (itype1,iunit1,iend,h,k,l,f,sigmaf,e,sigmae) if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then hkl=(h-hmin)*nk*nl+(k-kmin)*nl+(l-lmin) if (hkl.gt.imax) then write (ilp,1000) imax else jmin=min(jmin,hkl) jmax=max(jmax,hkl) data(1,hkl)=f data(2,hkl)=sigmaf data(5,hkl)=e data(6,hkl)=max(sigmae,e*sigmaf/f) end if end if end do rewind iunit2 do i=1,m2 call read1 (itype2,iunit2,iend,h,k,l,f,sigmaf,e,sigmae) if (sigmaf.gt.0.and.f.ge.cutoff*sigmaf) then hkl=(h-hmin)*nk*nl+(k-kmin)*nl+(l-lmin) if (hkl.gt.imax) then write (ilp,1000) imax else jmin=min(jmin,hkl) jmax=max(jmax,hkl) data(3,hkl)=f data(4,hkl)=sigmaf data(7,hkl)=e data(8,hkl)=max(sigmae,e*sigmaf/f) end if end if end do 1000 format (/1x, &'Warning: Packed index exceeds imax = ',i6/1x, &'Not all reflections can be stored in the data array.') c c Eliminate single, unpaired reflections from the data array. c j1=0 j2=0 do j=jmin,jmax f1=data(1,j) f2=data(3,j) if (f1.le.0.xor.f2.le.0) then if (f1.gt.0) j1=j1+1 if (f2.gt.0) j2=j2+1 do i=1,8 data(i,j)=0 end do end if end do write (ilp,'(/1x, &''j1 = '',i6,'' single, unpaired F1 data eliminated from the da'', &''ta set''/1x, &''j2 = '',i6,'' single, unpaired F2 data eliminated from the da'', &''ta set'')') j1,j2 c c Compile data statistics before scaling. c n=0 sumdf=0 sumf1=0 sumf2=0 sumde=0 sume=0 sum1sq=0 sum2sq=0 smin=9 smax=0 nc=0 do j=jmin,jmax f1=data(1,j) f2=data(3,j) if (f1.gt.0.and.f2.gt.0) then c c Accumulate agreement statistics. c n=n+1 sumdf=sumdf+abs(f1-f2) sumf1=sumf1+f1 sumf2=sumf2+f2 e1=data(5,j) e2=data(7,j) if (e1.gt.0.and.e2.gt.0) then sumde=sumde+abs(e1-e2) sume=sume+e1+e2 sum1sq=sum1sq+e1**2 sum2sq=sum2sq+e2**2 end if c c Unpack indices, hkl = (h - hmin)*nk*nl + (k - kmin)*nl + (l - lmin). c hkl=j h=hkl/(nk*nl) k=(hkl-h*nk*nl)/nl l=hkl-h*nk*nl-k*nl h=h+hmin k=k+kmin l=l+lmin c c Calculate sin(theta)/lambda. c s=sinthl(h,k,l,ginv) smin=min(smin,s) smax=max(smax,s) c c Count centric reflections in SAS data sets. c if (type.eq.'SAS') then call csym (ptgp,h,k,l,ic) if (ic.eq.1) nc=nc+1 end if c c Store an initial local scale factor of unity. c data(9,j)=1 end if end do q=sumf1/sumf2 rf=sumdf/(0.5*(sumf1+sumf2)) if (sume.gt.0) then re=sumde/(0.5*sume) e1sq=sum1sq/n e2sq=sum2sq/n end if write (ilp,'(/1x, &''N = '',i6,'' reflection pairs accepted for scaling'')') n if (nc.gt.0) write (ilp,'(/1x, &''nc = '',i6,'' centric reflections stored in''/1x, &'' '',6x,'' both the F1(+h+k+l) and F2(-h-k-l) subsets.''/1x, &'' '',6x,'' ---- ---'')') nc write (ilp,'(/1x, &''s = sin(theta)/lambda d = lambda/[2*sin(thet'', &''a)]''/1x, &''smin = '',f7.4,'' reciprocal Angstroms dmax = '',f5.1, &'' Angstroms''/1x, &''smax = '',f7.4,'' dmin = '',f6.2)') & smin,1/(2*smin),smax,1/(2*smax) write (ilp,'(/1x, &''Before scaling:''/1x, &''Q(F) = / = '',e10.3/1x, &''R(F) = /<0.5*(F1 + F2)> = '',e10.3/1x, &''R(E) = /<0.5*(E1 + E2)> = '',e10.3/1x, &'' = '',f6.3/1x,'' = '',f6.3)') & q,rf,re,e1sq,e2sq if (cutoff.gt.0) write (ilp,'(/1x, &''Values of .gt. 1 are possible since minimum F/sigma(F)'', &'' = '',f5.2)') cutoff return end c----------------------------------------------------------------------- subroutine decode (ptgp,laue,iptgp,ilaue) c c Point group and Laue group index numbers c c There are ten cases of alternative axes for seven of the 32 c crystallographic point groups, thus a total of 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 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 ' locscl/decode: Invalid input point group symbol.' 1 iptgp=i c c Laue group c if (i.ge. 1) ilaue=2 if (i.ge. 3) ilaue=5 if (i.ge. 6) ilaue=8 if (i.ge. 9) ilaue=11 if (i.ge.12) ilaue=16 if (i.ge.17) ilaue=18 if (i.eq.19.or.i.eq.21.or.i.eq.23) ilaue=23 if (i.eq.20.or.i.eq.22.or.i.eq.24) ilaue=24 if (i.ge.25) ilaue=27 if (i.ge.28) ilaue=32 if (i.ge.33) ilaue=34 if (i.ge.35) ilaue=37 if (i.ge.38) ilaue=39 if (i.ge.40) ilaue=42 laue=symbol(ilaue) 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 if (h.eq.0.and.k.eq.0.and.l.eq.0) return 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.ge.k.and.h.ge.-2*k) return x=h y=k i=-h-k h=i k=x if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return h=y k=i if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return h=y k=x l=-l go to 520 c c 31m c 21 if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return x=h y=k i=-h-k h=i k=x if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) return h=y k=i if (h.ge.0.and.h.ge.k.and.h.ge.-2*k) 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.ge.k.and.h.ge.-2*k.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.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0.and. & k.ge.0))) return h=y k=i if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0.and. & k.ge.0))) return h=y k=x if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0.and. & k.ge.0))) return h=x k=i if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.and.(l.gt.0.or.(l.eq.0.and. & k.ge.0))) return h=i k=y if (h.ge.0.and.h.ge.k.and.h.ge.-2*k.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.ge.k.and.h.gt.-2*k) return x=h y=k i=-h-k h=i k=x if (h.gt.0.and.h.ge.k.and.h.gt.-2*k) return h=y k=i if (h.gt.0.and.h.ge.k.and.h.gt.-2*k) 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 '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 (abs(h).eq.abs(k).and.abs(k).eq.abs(l)) then h=abs(h) k=abs(k) l=abs(l)*(h*k*l)/abs(h*k*l) return end if 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 (g,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 character g*5 integer h e=1 if (g.eq.'1 ') then m=1 else if (g.eq.'-1 ') then m=2 else if (g.eq.'2 ') then m=2 if (h.eq.0.and.l.eq.0) e=2 else if (g.eq.'m ') then m=2 if (k.eq.0) e=2 else if (g.eq.'2/m ') then m=4 if (k.eq.0) e=2 if (h.eq.0.and.l.eq.0) e=2 else if (g.eq.'222 ') then 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 else if (g.eq.'mm2 ') then m=4 if (h.eq.0) e=2 if (k.eq.0) e=2 if (h.eq.0.and.k.eq.0) e=4 else if (g.eq.'mmm ') then 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 else if (g.eq.'4 ') then m=4 if (h.eq.0.and.k.eq.0) e=4 else if (g.eq.'-4 ') then m=4 if (h.eq.0.and.k.eq.0) e=2 else if (g.eq.'4/m ') then m=8 if (l.eq.0) e=2 if (h.eq.0.and.k.eq.0) e=4 else if (g.eq.'422 ') then 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 else if (g.eq.'4mm ') then 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 else if (g.eq.'-42m ') then 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 else if (g.eq.'-4m2 ') then 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 else if (g.eq.'4/mmm') then 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 else if (g.eq.'3 ') then m=3 if (h.eq.0.and.k.eq.0) e=3 else if (g.eq.'-3 ') then m=6 if (h.eq.0.and.k.eq.0) e=3 else if (g.eq.'312 ') then 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 else if (g.eq.'321 ') then 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 else if (g.eq.'31m ') then 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 else if (g.eq.'3m1 ') then 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 else if (g.eq.'-31m ') then 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 else if (g.eq.'-3m1 ') then 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 else if (g.eq.'6 ') then m=6 if (h.eq.0.and.k.eq.0) e=6 else if (g.eq.'-6 ') then m=6 if (l.eq.0) e=2 if (h.eq.0.and.k.eq.0) e=3 else if (g.eq.'6/m ') then m=12 if (l.eq.0) e=2 if (h.eq.0.and.k.eq.0) e=6 else if (g.eq.'622 ') then 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 else if (g.eq.'6mm ') then 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 else if (g.eq.'-6m2 ') then 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 else if (g.eq.'-62m ') then 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 else if (g.eq.'6/mmm') then 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 else if (g.eq.'r3 ') then m=3 if (k.eq.h.and.l.eq.h) e=3 else if (g.eq.'r-3 ') then m=6 if (k.eq.h.and.l.eq.h) e=3 else if (g.eq.'r32 ') then 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 else if (g.eq.'r3m ') then 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 else if (g.eq.'r-3m ') then 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 else if (g.eq.'23 ') then 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 else if (g.eq.'m3 ') then 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 else if (g.eq.'432 ') then 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 else if (g.eq.'-43m ') then 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 else if (g.eq.'m3m ') then 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 end if m=m/e return end c----------------------------------------------------------------------- subroutine fsqtof (fsq,sigfsq,f,sigf) if (fsq.lt.0) then f=0 else f=sqrt(fsq) end if if (fsq.le.sigfsq) then sigf=0.5*sqrt(sigfsq) else sigf=0.5*sigfsq/f end if return end c----------------------------------------------------------------------- subroutine metric (a,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) deg=180/pi rad=pi/180 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 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----------------------------------------------------------------------- subroutine output (ifile,jfile,ilp,type,ptgp,ndata,nc,zlimit, & atime,adate,title) c c Compile Delta(F) distribution statistics and histogram, c optionally trim data from the distribution tails, and c write the hkl-sorted and Delta/sigma(Delta)-sorted output files. c parameter (jmax=25) dimension xh(jmax),yh(jmax) character atime*8,adate*9,title*72,type*3,ptgp*5 c c Fortran-90 dynamic storage allocation c real, allocatable::z(:) integer, allocatable::index(:) allocate (z(ndata),index(ndata)) c c Find median(F1 - F2). c rewind ifile 1000 format (1x,3i5,2(e12.5,e9.2),2(2f7.3)) do i=1,ndata read (ifile,1000) ih,ik,il,f1,s1,f2 z(i)=f1-f2 end do call sort (ndata,z,index) zmedian=z(index(nint(0.5*ndata))) c c Compile Delta(F) = F1 - F2 - median(F1 - F2) histogram. c zmean=0 do i=1,ndata d=z(i)-zmedian zmean=zmean+d z(i)=d end do zmin=z(index(1)) zmax=z(index(ndata)) zmean=zmean/ndata do j=1,jmax xh(j)=0 yh(j)=0 end do j=0 zj=zmin d=(zmax-zmin)/jmax do i=1,ndata zi=z(index(i)) do while (zi.ge.zj) j=j+1 zj=zj+d end do if (j.le.jmax) then xh(j)=xh(j)+zi yh(j)=yh(j)+1 end if end do xh(jmax)=xh(jmax)+zi yh(jmax)=yh(jmax)+1 c c In SAS cases, reduce the count in the histogram bin that contains the c centric-pair zero-differences. c if (type.eq.'SAS') then j=0 zj=zmin do while (zj.lt.0) j=j+1 zj=zj+d end do yh(j)=yh(j)-nc end if yhmax=0 do j=1,jmax if (yh(j).gt.0) then xh(j)=xh(j)/yh(j) yhmax=max(yhmax,yh(j)) end if end do c c Print statistics summary. c write (ilp,'(//1x, &''Delta(F) = F1 - F2 - median(F1 - F2) Distribution Statistics''/ &1x, &''------------------------------------------------------------''/ &1x, &''Before trimming distribution tails from the data set:''/1x, &''------''/1x, &'' N = '',i10,'' Delta(F) values''/1x, &'' median (F1 - F2) = '',e10.3/1x, &'' minimum Delta(F) = '',e10.3/1x, &'' maximum = '',e10.3/1x, &'' mean = '',e10.3)') & ndata,zmedian,zmin,zmax,zmean c c Store median(F1 - F2). c f0=zmedian c c Find the median and mean of the nonzero abs[Delta(F)] values. c nz=0 zmean=0 do i=1,ndata d=abs(z(i)) if (d.eq.0) nz=nz+1 zmean=zmean+d z(i)=d end do call sort (ndata,z,index) zmedian=z(index(nint(0.5*(ndata+nz)))) zmean=zmean/(ndata-nz) write (ilp,'(1x, &'' median abs[Delta(F)] = '',e10.3/1x, &'' mean = '',e10.3)') & zmedian,zmean c c Print Delta(F) histogram. c write (ilp,'(/1x, &''Delta(F) distribution before trimming distribution tails''/1x, &''--------------------------------------------------------''/1x, &'' n ''/1x,'' - ----------'')') do j=1,jmax n=nint(50*yh(j)/yhmax) write (ilp,'(1x,i6,e12.3,1x,50a1)') nint(yh(j)),xh(j), & ('X',i=1,n) end do c c Write the hkl-sorted, locally scaled output file "data.locscl" and c optionally reject data pairs with abnormal Delta(F) values in the c extreme tails of the Delta(F) distribution. c sigma=1.25*zmedian dlimit=zlimit*sigma nrej=0 nrejn=0 nrejp=0 rewind ifile open (unit=jfile,file='data.locscl',status='new') write (jfile,'(1x, &''Program locscl. '',a,'', '',a,''. '',a//1x, &'' h k l F1 sigmaF1 F2 sigmaF2'', &'' E1 sigmaE1 E2 sigmaE2''/1x, &'' - - - -- ------- -- -------'', &'' -- ------- -- -------'')') atime,adate,title j=0 zmean=0 do i=1,ndata read (ifile,1000) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 d=f1-f2-f0 if (dlimit.gt.0.and.abs(d).gt.dlimit) then nrej=nrej+1 if (d.lt.0) nrejn=nrejn+1 if (d.gt.0) nrejp=nrejp+1 if (nrej.eq.1) write (ilp,'(/1x, &''Reflection pairs trimmed from the distribution tails:''/1x, &''-----------------------------------------------------''/1x, &'' h k l F1 sigmaF1 F2 sigmaF2'', &'' E1 sigmaE1 E2 sigmaE2 (F1-F2)/sqrt(sigF1**2+sigF2**2)''/ &1x, &'' - - - -- ------- -- -------'', &'' -- ------- -- ------- ------------------------------'')') if (nrej.le.100) write (ilp, & '(1x,3i5,2(e12.5,e9.2),2(2f7.3),e10.2)') & ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2,(f1-f2)/sqrt(s1**2+s2**2) else write (jfile,1000) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 j=j+1 z(j)=d zmean=zmean+d end if end do close (unit=ifile,status='delete') endfile jfile i=ifile ifile=jfile jfile=i c c Recompile Delta(F) distribution statistics and histogram. c if (nrej.eq.0) then write (ilp,'(/1x, &''No data trimmed from distribution tails with zlimit = '',e10.3) &') zlimit else ndata=ndata-nrej call sort (ndata,z,index) zmin=z(index(1)) zmax=z(index(ndata)) zmedian=z(index(nint(0.5*ndata))) zmean=zmean/ndata do j=1,jmax xh(j)=0 yh(j)=0 end do j=0 zj=zmin d=(zmax-zmin)/jmax do i=1,ndata zi=z(index(i)) do while (zi.ge.zj) j=j+1 zj=zj+d end do if (j.le.jmax) then xh(j)=xh(j)+zi yh(j)=yh(j)+1 end if end do xh(jmax)=xh(jmax)+zi yh(jmax)=yh(jmax)+1 if (type.eq.'SAS') then j=0 zj=zmin do while (zj.lt.0) j=j+1 zj=zj+d end do yh(j)=yh(j)-nc end if yhmax=0 do j=1,jmax if (yh(j).gt.0) then xh(j)=xh(j)/yh(j) yhmax=max(yhmax,yh(j)) end if end do write (ilp,'(/1x, &''Delta(F) Distribution Statistics''/1x, &''--------------------------------''/1x, &''After trimming distribution tails:''/1x, &''-----''//1x, &'' maximum abs[Delta(F)]''/1x &'' zlimit = ------------------------- = '',e10.3/1x &'' 1.25*median abs[Delta(F)]''//1x, &'' Nrej = '',i10,'' rejected data pairs''/1x, &'' Nrejn = '',i10,'' negative Delta(F)''/1x, &'' Nrejp = '',i10,'' positive''//1x, &'' N = '',i10,'' remaining data pairs''/1x, &'' minimum Delta(F) = '',e10.3/1x, &'' maximum = '',e10.3/1x, &'' mean = '',e10.3)') & zlimit,nrej,nrejn,nrejp,ndata,zmin,zmax,zmean zmean=0 do i=1,ndata d=abs(z(i)) zmean=zmean+d z(i)=d end do call sort (ndata,z,index) zmedian=z(index(nint(0.5*(ndata+nz)))) zmean=zmean/(ndata-nz) write (ilp,'(1x, & '' median abs[Delta(F)] = '',e10.3/1x, & '' mean = '',e10.3)') & zmedian,zmean write (ilp,'(/1x, & ''Delta(F) distribution after trimming distribution tails''/1x, & ''-------------------------------------------------------''/1x, & '' n ''/1x,'' - ----------'')') do j=1,jmax n=nint(50*yh(j)/yhmax) write (ilp,'(1x,i6,e12.3,1x,50a1)') nint(yh(j)),xh(j), & ('X',i=1,n) end do end if c c Write an output file sorted on c c abs[Delta(E)] abs(E1 - E2) c abs[z(E)] = ------------- = --------------------------------- c sigma(Delta) sqrt[sigma(E1)**2 + sigma(E2)**2] c c or, in the absence of E values, the analogous z(F). c rewind ifile do i=1,4 read (ifile,'(a1)') dummy end do read (ifile,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 if (e1.gt.0.and.e2.gt.0) then ie=1 else ie=0 end if rewind ifile do i=1,4 read (ifile,'(a1)') dummy end do open (unit=jfile,status='scratch',form='unformatted', & access='direct',recl=44) do i=1,ndata read (ifile,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 if (ie.eq.1) then z(i)=abs(e1-e2)/sqrt(t1**2+t2**2) else z(i)=abs(f1-f2)/sqrt(s1**2+s2**2) end if write (jfile,rec=i) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 end do call sort (ndata,z,index) close (unit=ifile,status='keep') open (unit=ifile,file='sort.locscl',status='new') write (ifile,'(1x, &''Program locscl. '',a,'', '',a,''. '',a//1x, &'' h k l F1 sigmaF1 F2 sigmaF2'', &'' E1 sigmaE1 E2 sigmaE2 Delta/sigma(Delta)''/1x, &'' - - - -- ------- -- -------'', &'' -- ------- -- ------- ------------------'')') &atime,adate,title do i=ndata,1,-1 j=index(i) read (jfile,rec=j) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 if (ie.eq.1) then z(j)=sign(z(j),(e1-e2)) else z(j)=sign(z(j),(f1-f2)) end if write (ifile,'(1x,3i5,2(e12.5,e9.2),2(2f7.3),f8.2)') & ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2,z(j) end do close (unit=ifile,status='keep') close (jfile,status='delete') deallocate (z,index) c c Write an SAS output file with separate records for +h+k+l and -h-k-l. c if (type.eq.'SAS') then open (unit=ifile,file='data.locscl',status='old') open (unit=jfile,file='hphn.locscl',status='new') write (jfile,'(1x, & ''Program locscl. '',a,'', '',a,''. '',a//1x, & '' h k l F sigmaF E sigmaE''/1x, & '' - - - - ------ - ------'')') & atime,adate,title do i=1,4 read (ifile,'(a1)') dummy end do do i=1,ndata read (ifile,*) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 write (jfile,'(1x,3i5,e12.5,e9.2,2f7.3)') & +ih,+ik,+il,f1,s1,e1,t1 call csym (ptgp,ih,ik,il,ic) if (ic.eq.0) write (jfile,'(1x,3i5,e12.5,e9.2,2f7.3)') & -ih,-ik,-il,f2,s2,e2,t2 end do close (jfile,status='keep') end if close (ifile,status='keep') return end c----------------------------------------------------------------------- subroutine pscale (imax,data,minhkl,maxhkl,cutoff,ginv,ptgp, & ihmin,ikmin,ilmin,nk,nl,ilp) c c Global pre-scaling via a logarithmically linearized least-squares fit c of empirical scaling parameters p0, p1, and p2 to minimize c chi**2 = sum w*[(F1/F2) - p]**2, c where c p = p0*exp(p1*s**2 + p2*s**4) and s = sin(theta)/lambda. c c p1 approximates - , and c p2 approximates <(Biso1 - )**2> - <(Biso2 - )**2>. c character ptgp*5 dimension data(10,0:imax),ginv(3,3),a(3,3),ainv(3,3),b(3),c(3) real*8 sumw,sumx,sumy,sumx2,sumx3,sumx4,sumy2,sumxy,sumx2y, & sump,sumv,sumpsq n=0 sumw=0 sumx=0 sumy=0 sumx2=0 sumx3=0 sumx4=0 sumy2=0 sumxy=0 sumx2y=0 do ihkl=minhkl,maxhkl c c Get F and sigma(F) values. c f1=data(1,ihkl) s1=data(2,ihkl) f2=data(3,ihkl) s2=data(4,ihkl) c c Are both reflections significant? c if (s1.gt.0.and.s2.gt.0.and.f1.gt.cutoff*s1.and.f2.gt.cutoff*s2) & then c c Unpack indices, c c ihkl = (ih - ihmin)*nk*nl + (ik - ikmin)*nl + (il - ilmin). c ih=ihkl/(nk*nl) ik=(ihkl-ih*nk*nl)/nl il=ihkl-ih*nk*nl-ik*nl ih=ih+ihmin ik=ik+ikmin il=il+ilmin c c Calculate sin(theta)/lambda. c s=sinthl(ih,ik,il,ginv) c c Get point group multiplicity. c call esym (ptgp,ih,ik,il,epsilon,m) c c Acentric (ic = 0) or centric (ic = 1) reflection? c call csym (ptgp,ih,ik,il,ic) c c Double the multiplicity for acentric reflections. c if (ic.eq.0) m=2*m c c Store data for logarithmic fit. c x=s**2 y=log(f1/f2) w=m/((s1/f1)**2+(s2/f2)**2) c c Accumulate least-squares sums. c n=n+m 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 if 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 their c esd's. c b(1)=sumy b(2)=sumxy b(3)=sumx2y call mv (ainv,b,c) p0=c(1) p1=c(2) p2=c(3) c c chisq = sum w*(y - p)**2 c = sum w*y**2 + sum w*p**2 - 2*sum w*p*y c c p = p0 + p1*x + p2*x**2 c c p**2 = p0**2 + (p1*x + p2*x**2)**2 + 2*p0*(p1*x + p2*x**2) c = p0**2 + p1**2*x**2 + p2**2*x**4 + 2*p1*p2*x**3 c + 2*p0*p1*x + 2*p0*p2*x**2 c chisq=sumy2+p0**2*sumw+p1**2*sumx2+p2**2*sumx4 & +2*(p0*p1*sumx+p0*p2*sumx2+p1*p2*sumx3) & -2*(p0*sumy+p1*sumxy+p2*sumx2y) zsq=chisq/(n-3) 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 p2 coefficient of the term quadratic in s**2 is not c statistically significant, resort to a scaling function linear in s**2. c if (abs(p2).lt.sqrt(v2)) then det=sumw*sumx2-sumx*sumx p0=(sumx2*sumy-sumx*sumxy)/det p1=(-sumx*sumy+sumw*sumxy)/det chisq=sumy2+p0**2*sumw+p1**2*sumx2 & +2*p0*p1*sumx & -2*(p0*sumy+p1*sumxy) zsq=chisq/(n-2) v0=zsq*sumx2/det v1=zsq*sumw/det cov01=zsq*(-sumx/det) p2=0 v2=0 cov02=0 cov12=0 c c If the p1 coefficient of the term linear in s**2 is not statistically c significant, resort to a constant scaling function. c if (abs(p1).lt.sqrt(v1)) then p0=sumy/sumw p1=0 chisq=sumy2+p0**2*sumw-2*p0*sumy zsq=chisq/(n-1) v0=zsq/sumw v1=0 cov01=0 end if end if c c Scale the F2(derivative) values, keeping constant the F1(native) c values. c sumx=0 sumy=0 pmin=+1e9 pmax=-1e9 sumn=0 sump=0 sumv=0 sumpsq=0 do ihkl=minhkl,maxhkl f1=data(1,ihkl) s1=data(2,ihkl) f2=data(3,ihkl) s2=data(4,ihkl) if (s1.gt.0.and.s2.gt.0.and.f1.ge.cutoff*s1.and.f2.ge.cutoff*s2) & then ih=ihkl/(nk*nl) ik=(ihkl-ih*nk*nl)/nl il=ihkl-ih*nk*nl-ik*nl ih=ih+ihmin ik=ik+ikmin il=il+ilmin s=sinthl(ih,ik,il,ginv) p=p0+p1*s**2+p2*s**4 v=v0+v1*s**4+v2*s**8+2*(cov01*s**2+cov02*s**4+cov12*s**6) c c Convert from logarithmic fitting scale to antilogarithmic data scale. c p=exp(p) v=p**2*v c c Apply scaling function, and compile scaling statistics. c s2=sqrt((p*s2)**2+f2**2*v) f2=p*f2 data(3,ihkl)=f2 data(4,ihkl)=s2 sumx=sumx+f1 sumy=sumy+f2 if (p.lt.pmin) then pmin=p sigmin=sqrt(v) end if if (p.gt.pmax) then pmax=p sigmax=sqrt(v) end if sumn=sumn+1 sump=sump+p sumv=sumv+v sumpsq=sumpsq+p**2 end if end do c c Re-scale the re-scaled data to adjust from the reciprocal-variance- c weighted, logarithmic fitting scale to the unit-weighted data scale. c q=sumx/sumy pmin=pmin*q pmax=pmax*q sigmin=sigmin*q sigmax=sigmax*q sump=sump*q sumv=sumv*q**2 sumpsq=sumpsq*q**2 sumx=0 sumy=0 sumw=0 do ihkl=minhkl,maxhkl f1=data(1,ihkl) s1=data(2,ihkl) f2=data(3,ihkl) s2=data(4,ihkl) if (s1.gt.0.and.s2.gt.0.and.f1.ge.cutoff*s1.and.f2.ge.cutoff*s2) & then f2=f2*q s2=s2*q data(3,ihkl)=f2 data(4,ihkl)=s2 sumw=sumw+abs(f1-f2) sumx=sumx+f1 sumy=sumy+f2 end if end do write (ilp,'(//1x, &''Global Pre-scaling''/1x,''------------------''/1x, &''Empirical scaling parameters p0, p1, and p2 obtained by a''/1x, &''logarithmically linearized least-squares fit to minimize''/1x, &'' chi**2 = sum w*[(F1/F2) - p]**2 ,''/1x,''where''/1x, &'' p = p0*exp(p1*s**2 + p2*s**4) and s = sin(theta)/lamb'', &''da .''//1x, &''p1 approximates - , and''/1x, &''p2 approximates <(Biso(1) - )**2> - <(Biso(2) - )**2>.''//1x, &''p0 = '',e10.3/1x, &''p1 = '',e10.3/1x, &''p2 = '',e10.3,'' sqrt(abs(p2)) = '',e10.3)') & q*exp(p0),p1,p2,sqrt(abs(p2)) write (ilp,'(/1x, &''Pre-scaling Statistics:''/1x, &''-----------------------''/1x, &''pmin = '',e10.3,'' sigma(pmin) = ''e10.3/1x, &''pmax = '',e10.3,'' sigma(pmax) = ''e10.3//1x, &''

= '',e10.3/1x, &''**(1/2) = '',e10.3/1x, &''<(p -

)**2>**(1/2) = '',e10.3//1x, &''Q = / = '',e10.3/1x, &''R = /<0.5*(F1 + F2*p)> = '',e10.3)') & pmin,sigmin,pmax,sigmax, & sump/sumn,sqrt(sumv/sumn),sqrt((sumpsq/sumn)-(sump/sumn)**2), & sumx/sumy,sumw/(0.5*(sumx+sumy)) return end c----------------------------------------------------------------------- subroutine qfit (imax,data,minhkl,maxhkl,ihmin,ikmin,ilmin, & nk,nl,mh,mk,ml,nrej,ndata,qmin,sigmin,qmax,sigmax,qmean,rmsd,esd, & wqmean,wrmsd,wesd) c c Fit least-squares local scale factors in blocks of c n = [2*m(h) + 1]*[2*m(k) + 1]*[2*m(l) + 1] - 1 c reciprocal lattice points surrounding, but excluding, each point hkl. c dimension data(10,0:imax) parameter (jmax=25) dimension t(4,-jmax:+jmax,-jmax:+jmax,-jmax:+jmax) real*8 sumw,sumx,sumxsq,sumq,sumqsq,sumv,wsum,wsumq,wsumqsq c c Loop through the data set. c nrej=0 ndata=0 qmin=+1e9 qmax=-1e9 sumq=0 sumqsq=0 sumv=0 wsum=0 wsumq=0 wsumqsq=0 do ihkl=minhkl,maxhkl c c Get F and sigma(F) values. c f1=data(1,ihkl) s1=data(2,ihkl) f2=data(3,ihkl) s2=data(4,ihkl) if (f1.gt.0.and.s1.gt.0.and.f2.gt.0.and.s2.gt.0) then c c Fetch local block c c -mh .le. jh .le. +mh c -mk .le. jk .le. +mk c -ml .le. jj .le. +ml c c surrounding ih,ik,il stored by packed indices: c c ihkl = (ih - ihmin)*nk*nl + (ik - ikmin)*nl + (il - ilmin) c c jhkl = (ih - ihmin + jh)*nk*nl c + (ik - ikmin + jk)*nl c + (il - ilmin + jl) c c jhkl = ihkl + jh*nk*nl + jk*nl + jl c do jl=-ml,+ml do jk=-mk,+mk do jh=-mh,+mh do j=1,4 t(j,jh,jk,jl)=0 end do jhkl=ihkl+jh*nk*nl+jk*nl+jl if (jhkl.ne.ihkl.and. & minhkl.le.jhkl.and.jhkl.le.maxhkl) then do j=1,4 t(j,jh,jk,jl)=data(j,jhkl) end do end if end do end do end do c c Accumulate local least-squares sums for reflections surrounding, but c excluding, the given reflection. c n=0 sumw=0 sumx=0 sumxsq=0 do jl=-ml,+ml do jk=-mk,+mk do jh=-mh,+mh if (jh.ne.0.or.jk.ne.0.or.jl.ne.0) then f1=t(1,jh,jk,jl) s1=t(2,jh,jk,jl) f2=t(3,jh,jk,jl) s2=t(4,jh,jk,jl) if (f1.gt.0.and.s1.gt.0.and.f2.gt.0.and.s2.gt.0) then x=f1/f2 w=1/(x**2*((s1/f1)**2+(s2/f2)**2)) n=n+1 sumw=sumw+w sumx=sumx+w*x sumxsq=sumxsq+w*x**2 end if end if end do end do end do c c Test for ill-determined local scaling data. c if (n.lt.2.or.sumw.le.0.or.sumx.le.0) then q=0 chisq=0 else c c Evaluate the local scale factor, q = , and the local mean- c square error of fit. c q=sumx/sumw chisq=sumxsq-2*q*sumx+q**2*sumw end if c c Test for ill-determined error of fit. c if (chisq.le.0) then q=0 sigmaq=0 nrej=nrej+1 else c c Evaluate the local error-of-fit uncertainty, sigma(q). c sigmaq=sqrt(chisq/((n-1)*sumw)) c c Accumulate the q distribution statistics. c ndata=ndata+1 sumq=sumq+q sumqsq=sumqsq+q**2 sumv=sumv+sigmaq**2 w=1/sigmaq**2 wsum=wsum+w wsumq=wsumq+w*q wsumqsq=wsumqsq+w*q**2 if (q.lt.qmin) then qmin=q sigmin=sigmaq end if if (q.gt.qmax) then qmax=q sigmax=sigmaq end if end if c c Store the local scale factor. c data( 9,ihkl)=q data(10,ihkl)=sigmaq end if end do c c Evaluate the q-distribution statistics. c qmean=sumq/ndata rmsd=sqrt((sumqsq/ndata)-qmean**2) esd=sqrt(sumv/ndata) wqmean=wsumq/wsum wrmsd=sqrt((wsumqsq/wsum)-wqmean**2) wesd=sqrt(ndata/wsum) return end c----------------------------------------------------------------------- subroutine qscale (imax,data,minhkl,maxhkl,ihmin,ikmin,ilmin, & nk,nl,ifile,jfile,i012,type,ptgp,ndata,nc,ilp) c c Apply the local scale factors. c dimension data(10,0:imax) character type*3,ptgp*5 real*8 sumdf,sumde,sumf1,sumf2,sume,sum1sq,sum2sq,zmean c c Loop through the data set. c sumdf=0 sumde=0 sumf1=0 sumf2=0 sume=0 sum1sq=0 sum2sq=0 ndata=0 nc=0 open (unit=ifile,file='temp1.locscl',status='new') do ihkl=minhkl,maxhkl c c Get F and sigma(F), E and sigma(E), and q and sigma(q) values. c f1=data(1,ihkl) s1=data(2,ihkl) f2=data(3,ihkl) s2=data(4,ihkl) e1=data(5,ihkl) t1=data(6,ihkl) e2=data(7,ihkl) t2=data(8,ihkl) q=data( 9,ihkl) sigmaq=data(10,ihkl) c c Are both reflections present? c if (f1.gt.0.and.s1.gt.0.and.f2.gt.0.and.s2.gt.0.and.q.gt.0) then c c Unpack indices, hkl = (h - hmin)*nk*nl + (k - kmin)*nl + (l - lmin). c ih=ihkl/(nk*nl) ik=(ihkl-ih*nk*nl)/nl il=ihkl-ih*nk*nl-ik*nl ih=ih+ihmin ik=ik+ikmin il=il+ilmin c c Evaluate and write out the locally scaled data. c if (i012.eq.0) then c c Keep constant the mean of F1 and F2. c p1=f1+f2 p2=f1+f2*q v1=s1**2+s2**2 v2=s1**2+(s2*q)**2+(f2*sigmaq)**2 p=p1/p2 v=(v1+v2*(p1/p2)**2)/p2**2 c c The factor p = (F1 + F2)/(F1 + q*F2) preserves the mean of F1 and F2, c i.e., F1 + F2 = p*(F1 + q*F2). c s1=sqrt((p*s1)**2+f1**2*v) s2=sqrt((p*q*s2)**2+(f2*q)**2*v+(f2*p*sigmaq)**2) f1=f1*p f2=f2*p*q if (e1.gt.0.and.e2.gt.0) then t1=sqrt((p*t1)**2+e1**2*v) t2=sqrt((p*q*t2)**2+(e2*q)**2*v+(e2*p*sigmaq)**2) e1=e1*p e2=e2*p*q end if else if (i012.eq.1) then c c Keep constant the value of F1. c s2=sqrt((q*s2)**2+(f2*sigmaq)**2) f2=q*f2 if (e2.gt.0) then t2=sqrt((q*t2)**2+(e2*sigmaq)**2) e2=q*e2 end if else c c Keep constant the value of F2. c s1=sqrt((s1**2+(sigmaq*f1/q)**2)/q**2) f1=f1/q if (e1.gt.0) then t1=sqrt((t1**2+(sigmaq*e1/q)**2)/q**2) e1=e1/q end if end if c c For SAS data, average centric reflection pairs. c if (type.eq.'SAS') then call csym (ptgp,ih,ik,il,ic) if (ic.eq.1) then nc=nc+1 f1=0.5*(f1+f2) e1=0.5*(e1+e2) s1=0.5*sqrt(s1**2+s2**2) t1=0.5*sqrt(t1**2+t2**2) f2=f1 e2=e1 s2=s1 t2=t1 end if end if c c Accumulate agreement statistics. c sumdf=sumdf+abs(f1-f2) sumde=sumde+abs(e1-e2) sumf1=sumf1+f1 sumf2=sumf2+f2 sume=sume+e1+e2 sum1sq=sum1sq+e1**2 sum2sq=sum2sq+e2**2 c c Write a locally scaled output file. c ndata=ndata+1 write (ifile,1000) ih,ik,il,f1,s1,f2,s2,e1,t1,e2,t2 end if end do 1000 format (1x,3i5,2(e12.5,e9.2),2(2f7.3)) endfile ifile c c Evaluate and print the agreement statistics. c rf=sumdf/(0.5*(sumf1+sumf2)) if (sume.gt.0) then re=sumde/(0.5*sume) e1sq=sum1sq/ndata e2sq=sum2sq/ndata end if write (ilp,'(/1x,''After scaling:''/1x, &''Q(F) = / = '',f7.4/1x, &''R(F) = /<0.5*(F1 + F2)> = '',f7.4/1x, &''R(E) = /<0.5*(E1 + E2)> = '',f7.4/1x, &'' = '',f6.3/1x,'' = '',f6.3)') & sumf1/sumf2,rf,re,e1sq,e2sq return end c----------------------------------------------------------------------- subroutine read1 (itype,iunit,iend,h,k,l,f,sigf,e,sige) c c itype = 0 ASCII, formatted h,k,l,Fsq,sigFsq,F,sigF,E,sigE c 1 " " h,k,l,Fsq,sigFsq,F,sigF c 2 " " h,k,l,Fsq,sigFsq c 3 " " h,k,l, F,sigF,E,sigE c 4 " " h,k,l, F,sigF c 5 Binary, unformatted h,k,l,Fsq,sigFsq,F,sigF,E,sigE c 6 " " h,k,l,Fsq,sigFsq,F,sigF c 7 " " h,k,l,Fsq,sigFsq c 8 " " h,k,l, F,sigF,E,sigE c 9 " " h,k,l, F,sigF c integer h e=0 sige=0 1 if (itype.eq.0) read (iunit,*,end=9) h,k,l,y,sigy,f,sigf,e,sige if (itype.eq.1) read (iunit,*,end=9) h,k,l,y,sigy,f,sigf if (itype.eq.2) read (iunit,*,end=9) h,k,l,y,sigy if (itype.eq.3) read (iunit,*,end=9) h,k,l, f,sigf,e,sige if (itype.eq.4) read (iunit,*,end=9) h,k,l, f,sigf if (itype.eq.5) read (iunit, end=9) h,k,l,y,sigy,f,sigf,e,sige if (itype.eq.6) read (iunit, end=9) h,k,l,y,sigy,f,sigf if (itype.eq.7) read (iunit, end=9) h,k,l,y,sigy if (itype.eq.8) read (iunit, end=9) h,k,l, f,sigf,e,sige if (itype.eq.9) read (iunit, end=9) h,k,l, f,sigf if (itype.eq.2.or.itype.eq.7) then if (sigy.le.0) go to 1 call fsqtof (y,sigy,f,sigf) end if if (sigf.le.0) go to 1 return 9 iend=1 return end c----------------------------------------------------------------------- function sinthl(ih,ik,il,ginv) dimension h(3),ginv(3,3) h(1)=ih h(2)=ik h(3)=il q=0 do i=1,3 do j=1,3 q=q+h(j)*ginv(i,j)*h(i) end do end do sinthl=0.5*sqrt(q) 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 --------- -------- --- c Art of Scientific Computing, pp. 229-233. Cambridge, England: c --- -- ---------- --------- c Cambridge University Press. c dimension data(n),index(n) do i=1,n index(i)=i end do i=n/2+1 j=n 1 continue if (i.gt.1) then i=i-1 inext=index(i) t=data(inext) else inext=index(j) t=data(inext) index(j)=index(1) j=j-1 if (j.eq.1) then index(1)=inext return end if end if k=i l=i+i do while (l.le.j) 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 end do index(k)=inext go to 1 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=====================================================================