PROGRAM DATA_QUALITY c note: as of Dec 1998 the XSCALE program distributed by Wolfgang Kabsch as c part of the XDS suite has Rmrgd-F output. Rmeas is in both XDS and XSCALE. c c the CCP4 Program SCALA has Rmeas output. c c version Wed Jul 15 11:41:41 MDT 1998 KD some minor problems like c uninitialized avef, divide by zero in empty resol range and an err= clause c instead of end= for direct access XDS.HKL c also don't print out negative intensities implicit none C C Program for calculating data quality indicators vs resolution C from files generated by SCALEPACK or XDS C C modified to read SCALEPACK output produced using 'no merge original index' C C K. Diederichs & A. Karplus (1997) Nature Structural Biology 4, 269 C C Program version 1.0 - November 1996 c 1.1 - January 1997 weighting of I's for Rmrgd calculations c 1.1.1 Feb 2001 print out bad input lines in SCALEPACK files C C output quality indicators C C resolution is given both in Angstrom and as sin(theta)/lamda C ntotal - # collected reflections (number .le.0 in parentheses) C nuniq - # unique reflections (number .le.0 in parentheses) C n.ge.2 - # unique reflections with nobs > 1 C nsingl - # unique reflections with nobs = 1 C C Rmeas - Rmeas as described in the paper C Rsym - Rsym given in parentheses only for reference C purposes. Rmeas should always be reported. C PCV - The pooled coefficient of variation (see paper). C If the measured sigma values are accurate, the PCV C should be equal to 1/ of the raw data C I-over-sigma C measured data - based on all ummerged measured reflections C In parentheses is the same quantity based only on C reflections observed > 1 time (so as to compare with C the PCV) C merged data - for the final merged intensities C C fract - The distribution of how reflections contribute to the C overall Rmeas. C multi - The average multiplicity C C Imrgd1 - The Rmrgd-I value expected for the merged data, C assuming no systematic errors between the P and Q subsets, C and P and Q of same size (even for odd number of refl). C *** note for comparison that the Rmrgd-I value which C could result in the case of multiplicity 2 is exactly C Rmeas*sqrt(2).*** C Imrgd2 - The Rmrgd-I observed for random P and Q partitioning. C If 'Random P,Q subsets(1)' is chosen, Imrgd2 should be C almost equal to Imrgd1 (it will generally be a bit higher C because for an odd multiplicity refl P and Q differ by 1). C If 'based on frame ranges(2)' is chosen, then the P and Q C sizes match those of Imrgd3 to allow direct comparison, C and Imrgd2 is no longer directly comparable with Imrgd1. C Imrgd3 - The Rmrgd-I observed given the chosen non-random P and Q C partitioning. The Imrgd3 value will equal Riso-on-I between C data resulting from independent reduction of the P and Q C data subsets. Imrgd3=0 if random partitioning is chosen. C C Fmrgd1,2, and 3 - the same as above but for Rmrgd-F. Fmrgd2 based on C 'Random P,Q subsets(1)' is what we suggest should C be used as an estimate of the final reduced data C quality, for instance for comparisons with Rcryst C and Riso. As noted by Diederichs and Karplus, C the Fmrgd2 value in theory is underestimates the C data quality by up to a factor of sqrt(2). C - The average structure factor amplitudes output C CHARACTER NAME*80 Integer nbin,nref(101),ntot(101),nunique(101) integer ih,ik,il,iend, nin,itype,store,ifirst2,ilast2 integer nobstot,ifirst,ilast,i,nobs,nlores,nhires integer singlet(101),ipq,nuse(101),nneg1(101),nneg2(101),nneg, & npick real fpsq,fqsq,frac,resmid,xobs real reslim1,reslim2,temp,s2,fsq,sigfsq,delff,delf,fsqp, & fsqq,delf1,delf2 real stol3max,stolmax3,stolmin3,stolinc,stol2min,stol3mid real SUMIOS(101),iosmeas(101), $AVEF(101),fobs,ios,sumios2(101) real STOL2MAX(100),STOL2MID(100),RESOL(100),sig real sumrsym(101),sumrmeas(101), & sumpcv(101), sumdenom(101),fak,iosmeas2(101) real Imrgd1(101), Imrgd2(101), Imrgd3(101), sumI(101), sumnI(101), & Fmrgd1(101), Fmrgd2(101), Fmrgd3(101), sumF(101), sumnF(101), & sumpqi(101), sumpqf(101) DATA nhires/0/,nlores/0/,nref/101*0/,ntot/101*0/, & nobstot/0/,nneg1/101*0/,nneg2/101*0/,nneg/0/ & sumios/101*0./,sumrsym/101*0./, & sumrmeas/101*0./, & nunique/101*0/, & sumpcv/101*0./,sumdenom/101*0./ & ,singlet /101*0/,nuse/101*0/, iosmeas/101*0./, & iosmeas2/101*0./, sumios2 /101*0./ DATA Imrgd1/101*0./, Imrgd2/101*0./, Imrgd3/101*0./, & sumI/101*0./, sumnI/101*0./ , & Fmrgd1/101*0./, Fmrgd2/101*0./, Fmrgd3/101*0./, & sumF/101*0./, sumnF/101*0./, & sumpqi/101*0./,sumpqf/101*0./, delf1/0./, avef/101*0./ c store=1 on VAX and SGI, 4 otherwise e.g. Linux (see XDS docu) data store/1/ C C Open input and output files and get cell constants C write(*,*) 'SCALEPACK (1) or XDS (2)?' read(*,*) itype if (itype.eq.1) then WRITE(*,*) 'File name for SCALEPACK no merge original ... output' READ(*,1000) NAME OPEN(1,FILE=NAME,STATUS='OLD') elseif (itype.eq.2) then write(*,*) 'XDS unmerged output file (normally XDS.HKL)' read(*,1000) name open(1,file=name,access='direct', & recl=(store*34)/2,form='unformatted',status='old') endif C C Set up the equal volume resolution bins C WRITE(*,*) 'Input upper and lower resolution limits (Angstroms)' READ(*,*) RESLIM1,RESLIM2 IF (RESLIM1.LT.RESLIM2) THEN TEMP=RESLIM1 RESLIM1=RESLIM2 RESLIM2=TEMP ENDIF WRITE(*,*) 'Number of sin(theta)/lamda bins (up to 100 allowed)?' READ(*,*) NBIN STOLMAX3=(1./(2.*RESLIM2))**3 STOLMIN3=(1./(2.*RESLIM1))**3 STOLINC=(STOLMAX3-STOLMIN3)/NBIN STOL2MIN=EXP(ALOG(STOLMIN3)*2./3.) DO 10 I=1,NBIN STOL3MAX=STOLINC*I+STOLMIN3 STOL2MAX(I)=EXP(ALOG(STOL3MAX)*2./3.) STOL3MID=STOL3MAX-STOLINC/2. STOL2MID(I)=EXP(ALOG(STOL3MID)*2./3.) 10 CONTINUE C C get fram ranges for analysis C WRITE(*,*)'Random P,Q subsets(1) or based on frame ranges(2)?' read(*,*) ipq if (ipq.eq.1) then write(*,*) 'first, last frame number for analysis?' read(*,*) ifirst,ilast ifirst2 = -1 ilast2 = -1 else write(*,*) 'first, last frame number to be included in setP?' read(*,*) ifirst,ilast if (ifirst.ge.ilast) stop 'ifirst must be less than ilast' write(*,*) 'first, last frame number to be included in setQ?' read(*,*) ifirst2,ilast2 if (ifirst2.ge.ilast2) stop 'ifirst must be less than ilast' if (ifirst2.le.ilast) & stop 'first setQ frame must be > last setP frame' endif C C**** LOOP OVER ALL REFLECTIONS C CALL GETRES(IH,IK,IL,S2) nin = 0 C C read in the data for a given unique reflection and apply the resolution C cutoffs. Change the subroutine for new data formats. C iend=0 20 call gethkl(ih,ik,il,fpsq,fqsq,iend,nobs, & ifirst,ilast,fsq,sigfsq,delff,delf,itype,ifirst2, & ilast2,sig,ios,nneg,npick,fsqp,fsqq,delf1,delf2) if (iend.eq.1) goto 99 CALL GETRES(IH,IK,IL,S2) IF (MOD(NIN,1000).EQ.0) print *,'REF #', $ NIN,IH,Ik,IL,Fpsq,Fqsq,fsq,S2,iend IF (S2.LT.STOL2MIN) THEN NLORES= NLORES + 1 GOTO 20 ELSEif (S2.GT.STOL2MAX(nbin)) THEN NHIRES= NHIRES + 1 GOTO 20 ENDIF C C Count overall number of reflections read C NIN = NIN +1 nobstot = nobstot + nobs C C Increment the sums for the appropriate resolution bin C DO 70 I=1,NBIN IF (S2.GT.STOL2MAX(I)) GOTO 70 Nunique(I) = Nunique(I) + 1 ntot(i)=ntot(i)+nobs nneg1(i) = nneg1(i) + nneg iosmeas(i) = iosmeas(i) + ios sumios(i) = sumios(i) + fsq/sigfsq if (nobs.eq.1) singlet(i) = singlet(i) + 1 if (nobs.gt.1) then xobs = float(nobs) nref(i) = nref(i) + 1 nuse(i)=nuse(i)+nobs if (fsq.le.0.) nneg2(i) = nneg2(i) + 1 iosmeas2(i) = iosmeas2(i) + ios sumios2(i) = sumios2(i) + fsq/sigfsq C fak = sqrt(real(nobs)/real(nobs-1)) sumrsym(i)=sumrsym(i)+delff sumrmeas(i)=sumrmeas(i)+fak*delff sumpcv(i) = sumpcv(i) + sig sumdenom(i) = sumdenom(i) + (nobs-1) C fak = fak/sqrt(2.) sumnI(i)= sumnI(i)+ fsq*nobs sumI(i) = sumI(i) + fsq Imrgd1(i)=Imrgd1(I) + fak*2.*delff*sqrt(2./xobs) if (npick.gt.0) then Imrgd2(i)=Imrgd2(I) + abs(fpsq - fqsq) Imrgd3(i)=Imrgd3(I) + abs(fsqp - fsqq) sumpqi(i)=sumpqi(i) + fsq endif C fobs = sign(sqrt(abs(fsq)),fsq) sumnF(i)= sumnF(i)+ fobs*nobs sumF(i) = sumF(i) + fobs Fmrgd1(i)=Fmrgd1(I) + fak*2.*delf*sqrt(2./xobs) if (npick.gt.0) then Fmrgd2(i)=Fmrgd2(I) + delf2 Fmrgd3(i)=Fmrgd3(I) + delf1 sumpqf(i)=sumpqf(i) + fobs endif C endif GOTO 20 70 CONTINUE GOTO 20 C***********************************************************************!PAK 99 WRITE(*,*)'NUMBER UNIQUE REFLECTIONS IN RESOL RANGE',NIN WRITE(*,*)'NUMBER OBSERVED REFLECTIONS IN RESOL RANGE',nobstot WRITE(*,*)'NUMBER UNIQUE WITH TOO LOW RES',NLORES WRITE(*,*)'NUMBER UNIQUE WITH TOO HIGH RES',NHIRES C***********************************************************************!PAK C*****CALCULATION AND OUTPUT OF STATISTICS !PAK C***********************************************************************!PAK 999 DO 550 I=1,NBIN+1 if (i.eq.nbin+1) goto 545 C NREF(NBIN+1)=NREF(NBIN+1)+NREF(I) ntot(nbin+1)=ntot(nbin+1)+ntot(i) Nunique(NBIN+1)=Nunique(NBIN+1)+Nunique(I) Nuse(NBIN+1)=Nuse(NBIN+1)+Nuse(I) nneg1(nbin+1)=nneg1(nbin+1)+nneg1(i) nneg2(nbin+1)=nneg2(nbin+1)+nneg2(i) singlet(nbin+1)=singlet(nbin+1)+singlet(i) RESOL(I) = SQRT(1./STOL2Max(I))/2. C sumI(nbin+1) = sumI(nbin+1) + sumI(i) sumnI(nbin+1) = sumnI(nbin+1) + sumnI(i) sumpqi(nbin+1) = sumpqi(nbin+1) + sumpqi(i) Imrgd1(nbin+1) = Imrgd1(nbin+1) + Imrgd1(i) Imrgd2(nbin+1) = Imrgd2(nbin+1) + Imrgd2(i) Imrgd3(nbin+1) = Imrgd3(nbin+1) + Imrgd3(i) sumF(nbin+1) = sumF(nbin+1) + sumF(i) sumnF(nbin+1) = sumnF(nbin+1) + sumnF(i) sumpqf(nbin+1) = sumpqf(nbin+1) + sumpqf(i) Fmrgd1(nbin+1) = Fmrgd1(nbin+1) + Fmrgd1(i) Fmrgd2(nbin+1) = Fmrgd2(nbin+1) + Fmrgd2(i) Fmrgd3(nbin+1) = Fmrgd3(nbin+1) + Fmrgd3(i) C sumrsym(nbin+1)=sumrsym(nbin+1)+sumrsym(i) sumrmeas(nbin+1)=sumrmeas(nbin+1)+sumrmeas(i) sumpcv(nbin+1)=sumpcv(nbin+1)+sumpcv(i) sumdenom(nbin+1)=sumdenom(nbin+1) + sumdenom(i) iosmeas(nbin+1) = iosmeas(nbin+1) + iosmeas(i) iosmeas2(nbin+1) = iosmeas2(nbin+1) + iosmeas2(i) sumios(nbin+1) = sumios(nbin+1) + sumios(i) sumios2(nbin+1) = sumios2(nbin+1) + sumios2(i) 545 continue if (sumI(i).gt.0.) & sumpcv(i)=100.*nref(i)*sqrt(sumpcv(i)/sumdenom(i))/sumI(i) if (sumnI(i).gt.0.) then sumrsym(i) =100.*sumrsym(i) /sumnI(i) sumrmeas(i)=100.*sumrmeas(i)/sumnI(i) Imrgd1(i)=100.*Imrgd1(i)/sumnI(i) endif if (sumnI(i).gt.0.) then Imrgd2(i)=100.*Imrgd2(i)/sumpqi(i) Imrgd3(i)=100.*Imrgd3(i)/sumpqi(i) endif if (sumnF(i).gt.0.) Fmrgd1(i)=100.*Fmrgd1(i)/sumnF(i) if (sumpqf(i).gt.0.) then Fmrgd2(i)=100.*Fmrgd2(i)/sumpqf(i) Fmrgd3(i)=100.*Fmrgd3(i)/sumpqf(i) endif IF (NREF(I).NE.0) AVEF(I)= sumF(I)/NREF(I) if (nunique(i).ne.0) sumios(i) = sumios(I)/nunique(i) if (nref(i).ne.0) sumios2(i) = sumios2(I)/nref(i) if (ntot(i).ne.0) iosmeas(i) = iosmeas(I)/ntot(i) if (nuse(i).ne.0) iosmeas2(i)= iosmeas2(I)/nuse(i) 550 CONTINUE C WRITE (8,1500) DO 150 I = 1 , NBIN resmid = SQRT(1./STOL2MID(I))/2. WRITE (8,1520) I, $ RESOL(I),1/(2.*RESOL(I)),resmid,1/(2.*resmid), $ ntot(i),nneg1(i), nunique(i),nneg2(i),nref(i), $ singlet(i) 150 CONTINUE I = NBIN + 1 WRITE (8,1519) $ ntot(i), nneg1(i),nunique(i),nneg2(i), $ nref(i), singlet(i) C WRITE (8,1501) WRITE (8,1503) DO 160 I = 1 , NBIN resmid = SQRT(1./STOL2MID(I))/2. frac = 100.*float(nuse(i))/float(nuse(nbin+1)) WRITE (8,1525) I,RESmid,sumrmeas(i),sumrsym(i),sumpcv(i), $ iosmeas(i),iosmeas2(i), sumios(i),sumios2(i), frac, $ real(ntot(i))/max(1,nunique(i)) c the max() is only there to avoid divide by zero 160 CONTINUE I = NBIN + 1 frac = 100.*float(nuse(i))/float(nuse(nbin+1)) WRITE (8,1524) sumrmeas(i),sumrsym(i),sumpcv(i), $ iosmeas(i),iosmeas2(i), sumios(i),sumios2(i), frac, $ real(ntot(i))/max(1,nunique(i)) c the max() is only there to avoid divide by zero C WRITE (8,1502) WRITE (8,1504) DO 170 I = 1 , NBIN resmid = SQRT(1./STOL2MID(I))/2. WRITE (8,1530) I,RESmid,Imrgd1(i),imrgd2(i),imrgd3(i), $ Fmrgd1(i),Fmrgd2(i),Fmrgd3(i),avef(i) 170 CONTINUE I = NBIN + 1 WRITE (8,1529) Imrgd1(i),imrgd2(i),imrgd3(i), $ Fmrgd1(i),Fmrgd2(i),Fmrgd3(i),avef(i) if (ifirst2.ne.-1) then WRITE(8,1505) min(ifirst,ifirst2),max(ilast2,ilast),NAME else WRITE(8,1505) ifirst,ilast,NAME endif STOP 1000 FORMAT(A) 1500 FORMAT('Bin Resol limit & midpoint ', $ 'Ntotal (#<=0) Nuniq (#<=0) N.ge.2 Nsingl') 1501 FORMAT(/,'Bin Res ', $ ' Rmeas ( Rsym) PCV ', $ 'fract multi') 1503 FORMAT(' (A) ', $ ' (%) (%) (%) measured data merged data', $ ' (%)') 1504 FORMAT(' (A) (%) (%) (%) ', $ ' (%) (%) (%)') 1502 FORMAT(/,'Bin Res ', $ 'Imrgd1 Imrgd2 Imrgd3 Fmrgd1 Fmrgd2 Fmrgd3 ') 1505 FORMAT(/,' Frames ',I8,' to ',I8,' from the file:',/,A80) 1519 FORMAT(' Totals ',19x,I9,'(',I5,')',I7,' (',I4,')',2I7,I9) 1520 FORMAT(I3,2(f6.2,f6.3),I9,' (',I4,')',I7,' (',I4,')',2I7,I9) 1524 FORMAT(' Totals ',f8.1,' (',f5.1,')',F8.1, $ 2(f7.1,' (',F4.1,')'),2f6.1) 1525 FORMAT(I3,f6.2,f7.1,' (',f5.1,')',F8.1, $ 2(f7.1,' (',F4.1,')'),2f6.1) 1529 FORMAT(' Totals ',3F8.1,2x,3F8.1,F8.2) 1530 FORMAT(I3,f6.2,f7.1,2F8.1,2x,3F8.1,F8.2) END C*********************************************************************** SUBROUTINE GETRES(IH,IK,IL,S2) integer icall,ih,ik,il real as,bs,cs,cosgs,cosbs,cosas,cosa,cosb,cosg, & sina,sinb,sing real s2,a,b,c,alf,bet,gam,vol common /cell/ as,bs,cs,cosgs,cosbs,cosas save icall DATA ICALL /0/ IF (ICALL.EQ.1) THEN S2=((IH*AS)**2+(IK*BS)**2+(IL*CS)**2+2*IH*IK*AS*BS*COSGS $+2*IH*IL*AS*CS*COSBS+2*IK*IL*BS*CS*COSAS)/4. c STOL = SQRT(S2) c RESOL=1/(2.*STOL) ELSE WRITE(*,*) 'Unit cell parameters (a,b,c,alf,bet,gam)' READ(*,*) A,B,C,ALF,BET,GAM ALF=ALF/57.2958 BET=BET/57.2958 GAM=GAM/57.2958 COSA=COS(ALF) COSB=COS(BET) COSG=COS(GAM) SINA=SIN(ALF) SINB=SIN(BET) SING=SIN(GAM) COSAS=(COSB*COSG-COSA)/(SINB*SING) COSBS=(COSA*COSG-COSB)/(SINA*SING) COSGS=(COSA*COSB-COSG)/(SINA*SINB) VOL=A*B*C*SQRT(1-COSA**2-COSB**2-COSG**2+2*COSA*COSB*COSG) AS=B*C*SINA/VOL BS=A*C*SINB/VOL CS=A*B*SING/VOL ICALL = 1 ENDIF RETURN END C************************************************************************ SUBROUTINE gethkl(IHold,IKold,ILold,FPSQ,FMSQ,IEND, & nobs,ifirst,ilast,fsq,sigfsq,delff,delf,itype,ifirst2, & ilast2,sig,ios,nneg,npick,fsqp,fsqq,delf1,delf2) c new 26.1.97 Kay: average I's considering their weights c 23.2.01 Kay: treat bad input lines implicit none integer ih,ik,il,iend,itype,irec,nneg integer icall,ifrm,nskip,ihold,ikold,ilold,ifirst,nq, & ilast,i,nobs,minfrm,maxfrm,ifirst2,ilast2,npick real fpsq,fmsq,tmp_fsq,tmp_sigfsq,sig,ios,delf1,delf2 real fsq,sigfsq,ff(100),sumfsq,sumsigfsq,delff,delf real ftmp,fitmp,gg(100),fsqp,fsqq,sd(100),te(100) real sigfmsq,sigfpsq,sigfsqp,sigfsqq integer*2 ihkl(3),idum6(6),idum5(5),ifrm2,idum1 logical new,bytswp character*38 buf,swap character*80 line c the following variables should be saved across calls: save icall,minfrm,maxfrm,irec common idum1,ihkl,idum6,tmp_fsq,tmp_sigfsq,idum5,ifrm2 equivalence (swap,ihkl) data icall/0/,minfrm/100000/,maxfrm/-100000/, & irec/0/,bytswp/.false./ C 1234 format(12x,3i4,i6,7x,2f8.1) new=.true. if (iend.eq.2) then print*,'first and last frame in file is',minfrm,maxfrm iend=1 return endif c read leading lines for option 'no merge original indices' if (itype.eq.1.and.icall.eq.0) then icall=1 read(1,'(i5)') nskip do 100 i=1,nskip*2 100 read(1,'(x)') endif nobs=0 nneg = 0 sumfsq=0. sumsigfsq=0. fmsq=0. sigfmsq=0. fpsq=0. sigfpsq=0. fsqp=0. sigfsqp=0. fsqq=0. sigfsqq=0. ios = 0. nq=0 10 if (itype.eq.1) then read(1,1234,err=12,end=99)ih,ik,il,ifrm,fsq,sigfsq else irec=irec+1 read(1,rec=irec,err=99) buf if (bytswp) then call swab(swap,buf) else swap=buf endif ih=ihkl(1) if (ih.eq.10000) goto 99 ik=ihkl(2) il=ihkl(3) ifrm=ifrm2 fsq=tmp_fsq sigfsq=tmp_sigfsq c if (fsq.lt.0.) print*,ih,ik,il,ifrm,fsq,sigfsq endif minfrm=min(minfrm,ifrm) maxfrm=max(maxfrm,ifrm) if (ifrm.ge.ifirst.and.ifrm.le.ilast) goto 11 if (ifrm.ge.ifirst2.and.ifrm.le.ilast2) goto 11 goto 10 12 print*,'a read error occured in this line:' backspace(1) read(1,'(a80)') line print*,line stop 'file read problem' 11 if (new.or.(ih.eq.ihold.and.ik.eq.ikold.and.il.eq.ilold)) then if (new) then ihold=ih ikold=ik ilold=il new=.false. endif nobs=nobs+1 if (nobs.gt.100) stop 'nobs > 100' sumfsq=sumfsq+fsq sumsigfsq=sumsigfsq+1./sigfsq**2 ff(nobs)=fsq sd(nobs)=sigfsq if (fsq.le.0.) nneg = nneg + 1 ios = ios + fsq/sigfsq c systematic error estimation: assign non-random P and Q based on ifrm if (ifirst2.ne.-1) then if (ifrm.ge.ifirst2.and.ifrm.le.ilast2) then nq=nq+1 fsqq=fsqq+fsq/sigfsq**2 sigfsqq=sigfsqq+1./sigfsq**2 else fsqp=fsqp+fsq/sigfsq**2 sigfsqp=sigfsqp+1./sigfsq**2 endif endif c end of ifrm-based assignment to P and Q goto 10 else if (itype.eq.1) then backspace(1) else irec=irec-1 endif endif goto 101 99 IEND = 2 101 continue c the following deals with Rsym, Rmeas if (nobs.gt.0) then fsq=sumfsq/nobs sigfsq=sqrt(1./sumsigfsq) delff=0. delf=0. sig=0. do 200 i=1,nobs ftmp = sign(sqrt(abs(fsq)),fsq) fitmp = sign(sqrt(abs(ff(i))),ff(i)) delf=delf+abs(fitmp-ftmp) delff=delff+abs(ff(i)-fsq) sig = sig + (ff(i)-fsq)**2 200 continue endif if (nobs.lt.2) goto 999 c the following deals with Rmrgd if (ifirst2.eq.-1) then c P and Q of "equal" size, for Rmrgd to be compared w/ Rcryst, Rfree npick=nobs/2 elseif (nq.gt.0.and.nobs.gt.nq) then c P and Q of possibly unequal size (comparison w/ Riso) npick=nq fsqp=fsqp/sigfsqp fsqq=fsqq/sigfsqq c don't know if sigfsqp,sigfsqq will be used later but I compute them anyway sigfsqp=1./sqrt(sigfsqp) sigfsqq=1./sqrt(sigfsqq) ftmp = sign(sqrt(abs(fsqp)),fsqp) fitmp = sign(sqrt(abs(fsqq)),fsqq) delf1=abs(fitmp-ftmp) else c doesn't contribute to syst. err. estim. as either np=nobs-nq or nq is 0 npick=0 endif if (npick.gt.0) then call pickmofn(ff,sd,nobs,gg,te,npick) do 300 i=1,npick fpsq=fpsq+gg(i)/te(i)**2 sigfpsq=sigfpsq+1./te(i)**2 300 continue do 400 i=npick+1,nobs fmsq=fmsq+gg(i)/te(i)**2 sigfmsq=sigfmsq+1./te(i)**2 400 continue fpsq=fpsq/sigfpsq fmsq=fmsq/sigfmsq c don't know if sigfpsq,sigfmsq will be used later but I compute them anyway sigfpsq=1./sqrt(sigfpsq) sigfmsq=1./sqrt(sigfmsq) ftmp = sign(sqrt(abs(fpsq)),fpsq) fitmp = sign(sqrt(abs(fmsq)),fmsq) delf2=abs(fitmp-ftmp) endif 999 RETURN END c function ran(i) c primitive random number generator implicit none integer i real ran i=mod(151*i+1,20011) ran=i/20011. return end c subroutine swab(swap,buf) c byte swapping, for integer*2 that's easy but for reals also exchange c the first with the second half! tested with SGI-xds.hkl on Linux PC integer i character*38 swap,buf character*2 c2 do 100 i=1,37,2 100 swap(i:i)=buf(i+1:i+1) do 200 i=2,38,2 200 swap(i:i)=buf(i-1:i-1) c2=swap(19:20) swap(19:20)=swap(21:22) swap(21:22)=c2 c2=swap(23:24) swap(23:24)=swap(25:26) swap(25:26)=c2 return end c subroutine pickmofn(arr,sd,n,brr,te,m) c purpose: randomly pick m out of n entries in arr and put them into the c beginning of brr. The remaining entries of arr go to the end of brr. KD 11/96 c 26.1.97 KD: sd is passively assigned to te the same way as arr to brr implicit none integer n,m,i,iseed,ind real arr(n),brr(n),sd(n),te(n),ran,save data iseed/12345/ do 100 i=1,n te(i)=sd(i) 100 brr(i)=arr(i) do 200 i=1,m c choose at random an index ind such that i<=ind<=n . requires 0<=ran<1 ind=ran(iseed)*(n-i+1)+i if (ind.lt.i.or.ind.gt.n) stop 'check ran() generator' c exchange brr(i) and brr(ind) save=brr(i) brr(i)=brr(ind) brr(ind)=save save=te(i) te(i)=te(ind) te(ind)=save 200 continue return end