c c ---------------------------------------------------------------------- c subroutine phextsetup(hkls) include 'dmmulti.fh' include 'io.fh' include 'xtlparam.fh' include 'xhklinfo.fh' include 'cyclsync.fh' c PHEXTSETUP - set up phase extension. c We must first select a starting set of data. c real hkls(nhklc,0:*) c integer nresfine,nres,nhist parameter (nresfine=60,nres=30,nhist=100) real refls(nresfine),sfom(nresfine),sigf2(nres),hist(nhist) real hfom(0:100) c integer i,j,ires,ifom,imin,imax,idmax,ip,nrefli real s,ds,scut,fom,hsum,fdrop,fdmax,sfomlo,sfomhi,refllo,reflhi real fo,eo,eps,pmag,pres,pfom,pcomb,rflcut,rfltot,fracinc,flag real pcut(maxcyc) logical centre c integer lpcomb parameter (lpcomb=lcycin) c real cumprc,cumprnc external cumprc,cumprnc c do 50 i=1,nresfine refls(i)=0.0 sfom(i)=0.0 50 continue do 55 i=0,100 hfom(i)=0.0 55 continue c ds=smax/(nresfine-1) c first make a histogram of FOM's to find out where the good data goes out to do 100 i=1,nrflo c get data from array s=hkls(4,i) fom=hkls(lfomo,i) ires=int(s/ds)+1 refls(ires)=refls(ires)+1.0 sfom(ires)=sfom(ires)+fom ifom=nint(100.0*fom) hfom(ifom)=hfom(ifom)+1.0 100 continue c c now search for the biggest drop in mean fom either side of an arbitrary c resolution cutoff (but inc 10% to 80% of reflections) imin=nresfine*(0.10**0.666667) imax=nresfine*(0.70**0.666667) idmax=-1 fdmax=0.0 do 200 i=imin,imax sfomlo=0.0 refllo=0.0 sfomhi=0.0 reflhi=0.0 do 180 j=1,i-1 sfomlo=sfomlo+sfom(j) refllo=refllo+refls(j) 180 continue do 190 j=i,nresfine sfomhi=sfomhi+sfom(j) reflhi=reflhi+refls(j) 190 continue fdrop=sfomlo/max(refllo,1.0)-sfomhi/max(reflhi,1.0) if (fdrop.gt.fdmax) then idmax=i fdmax=fdrop endif 200 continue c c Identify the nature of the input data: if (wmag+wfom+wres.ne.0.0) then c c first deal with a user extension scheme if (fracinit.eq.0.0) fracinit=0.3 nrefli=nrflo*fracinit write (lpt,250) 250 format (/' USER EXTENSION SCHEME SELECTED ') if (wmag+wfom.eq.0) write(lpt,260) if (wmag+wres.eq.0) write(lpt,270) if (wres+wfom.eq.0) write(lpt,220) 260 format (/' Extension is in steps of resolution only.') 270 format (/' Extension is on the basis of FOMs only.') 220 format (/' Extension is on the basis of magnitudes only.') write (lpt,290)nrefli,nrflo 290 format ( + /' Approximate number of starting reflections ',i7,' of ',i7) c else if (idmax.ge.0.and.fdmax.gt.0.25) then c c deal with phase extension-like cases: c now count number of reflections to include scut=idmax*ds nrefli=0 do 300 i=1,nrflo s=hkls(4,i) fom=hkls(lfomo,i) if (s.lt.scut.and.fom.ne.0.0) nrefli=nrefli+1 300 continue c if (fracinit.eq.0.0) fracinit=(scut/smax)**1.5 nrefli=nrflo*fracinit wmag=1.0 wfom=1.0 wres=2.0 c write (lpt,310)sqrt(1.0/scut),nrefli,nrflo 310 format ( + /' INPUT DATA IDENTIFIED: Data appears to have good FOMs at low' + /' resolution and poor FOMs at high resolution. Calculation' + /' will be phase extension by resolution.' + /' Approximate resolution cutoff for starting data: ',f7.3,' A' + /' Approximate number of starting reflections ',i7,' of ',i7) c else c c now deal with phase refinement type cases: if (fracinit.eq.0.0) fracinit=0.3 nrefli=nrflo*fracinit wmag=1.0 wfom=1.0 wres=1.0 c write (lpt,320)nrefli,nrflo 320 format ( + /' INPUT DATA IDENTIFIED: Data appears to have approximate' + /' phasing at all resolutions. Calculation will be phase' + /' refinement by magnitude and FOM.' + /' Approximate number of starting reflections ',i7,' of ',i7) endif c c Now write out a graph of FOMs c write (lpt,510) 510 format (/ +' ------------------ xloggraph ------------------'/ +' \$TABLE: Analysis of input figures of merit :'/ +' \$GRAPHS: Mean FOM with resolution :N:2,4:'/ +' \$\$'/ +' res 4sin2(th) number \$\$'/ +' (lim) /lambda2 (reflns) \$\$') do 520 i=1,nresfine-1 write (lpt,530) + 1.0/sqrt(i*ds),i*ds,refls(i),sfom(i)/max(refls(i),1.0) 520 continue 530 format(1x,f5.2,2x,f6.3,4x,f7.0,4x,f6.3) write (lpt,540) 540 format (' \$\$'/' ------------------ xloggraph ------------------') c c write (lpt,*)' wmag,wfom,wres ',wmag,wfom,wres c c Now we pick the number of reflections to use at each resolution c do 650 i=1,nres refls(i)=0.0 sigf2(i)=0.0 650 continue do 660 i=1,nhist hist(i)=0.0 660 continue c ds=smax/(nres-1) c get stuff to calc E's, number of triplet terms do 700 i=1,nrflo c get data from array s=hkls(4,i) eps=abs(hkls(5,i)) centre=(hkls(5,i).lt.0.0) fo=hkls(lfo,i) ires=int(s/ds)+1 refls(ires)=refls(ires)+1.0 sigf2(ires)=sigf2(ires)+fo**2/eps 700 continue c do 750 i=1,nres sigf2(i)=sigf2(i)/max(refls(i),1.0) 750 continue c c get histogram of reflection pcomb's do 800 i=1,nrflo c get data from array s=hkls(4,i) eps=abs(hkls(5,i)) centre=(hkls(5,i).lt.0.0) fo=hkls(lfo,i) fom=hkls(lfomo,i) ires=int(s/ds)+1 c form the reflection quality parameters eo=fo/max(sqrt(sigf2(ires)*eps),1.0) if (centre) then pmag=cumprc(eo) else pmag=cumprnc(eo) endif pfom=fom pres=1.0-(s/smax)**1.5 pcomb=(pmag*wmag+pfom*wfom+pres*wres)/(wmag+wfom+wres) hkls(lpcomb,i)=pcomb ip=int(pcomb*(nhist-1))+1 hist(ip)=hist(ip)+1.0 800 continue c c find pcut to include the correct number of reflections do 830 j=1,ncycl-1 fracinc=fracinit*(1.0/fracinit)**(real(j-1)/real(ncycl-1)) rflcut=(1.0-fracinc)*nrflo rfltot=0.0 do 810 i=1,nhist-1 rfltot=rfltot+hist(i) if (rfltot.gt.rflcut) goto 820 810 continue 820 pcut(j)=(i-1-(rfltot-rflcut)/hist(i))/(nhist-1) 830 continue pcut(ncycl)=0.0 c c now label the reflections ip=0 do 900 i=1,nrflo c set flag for incl/excl, or count free-r reflections pcomb=hkls(lpcomb,i) do 850 j=1,ncycl-1 if (pcomb.ge.pcut(j)) goto 860 850 continue 860 hkls(lcycin,i)=real(j) if (j.eq.1) ip=ip+1 900 continue c return end c c ---------------------------------------------------------------------- c subroutine phextstep(hkls,print,reseff) include 'dmmulti.fh' include 'io.fh' include 'xtlparam.fh' include 'xhklinfo.fh' include 'cyclsync.fh' c PHEXTSTEP - take a phase extension step. c - work out some sort of intelligent criterion for including c refllections. This includes some weighted combination of: c 1. The magnitude of the structure factor (converted to a probability) c 2. The estimated FOM attached to its phase c 3. A poor estimate of the number of triplet relns for that refln c All of the above are normalised to a scale of 0...1 c c Included reflections are flagged with +1 c c Maximum Entropy results suggest that it would be good to include S.F.s c which are poorly extrapolated by density modification. This applies for c us if the extrapolated F is large and its obs. value is small, or if the c MIR phase is good and phase or magnitude are significantly different. c real hkls(nhklc,0:*) real reseff logical print c integer i,ninc,nexc real s,fo real sfoinc,sfoexc,flag,fninc,fsinc c c c now include reflections ninc=0 sfoinc=0.0 nexc=0 sfoexc=0.0 do 300 i=1,nrflo c get data from array s=hkls(4,i) fo=hkls(lfo,i) flag=hkls(lflag,i) c set flag for incl/excl, or count free-r reflections if (icycl.ge.nint(hkls(lcycin,i))) then ninc=ninc+1 sfoinc=sfoinc+fo hkls(lflag,i)=1.0 else nexc=nexc+1 sfoexc=sfoexc+fo hkls(lflag,i)=0.0 endif 300 continue c print=(icycl.eq.1) fninc=real(ninc)/real(ninc+nexc) fsinc=sfoinc/(sfoinc+sfoexc) reseff=rmax/fninc**0.3333 c write (lpt,220)icycl,ncycl,reseff,ninc,nexc 220 format (/' PHASE EXTENSION CYCLE ',i3,' OF ',i3 +/' Effective resolution: ',f6.3,' A', +/' Number of reflections included: ',i7, +/' Number of reflections excluded: ',i7) return end c c ---------------------------------------------------------------------- c logical function converged() include 'dmmulti.fh' include 'cyclsync.fh' c CONVERGED - convergence check for phase extension c icycl=icycl+1 converged=icycl.gt.ncycl return end c c ---------------------------------------------------------------------- c subroutine setxfr(hkls) include 'dmmulti.fh' include 'xhklinfo.fh' c SETXFR - Set up cross-validated free-r weighting real hkls(nhklc,0:*) c integer i c do 100 i=1,nrflo hkls(lffree,i)=-1.0 100 continue return end c c ---------------------------------------------------------------------- c subroutine accxfr(hkls,lftmp,lphitmp) include 'dmmulti.fh' include 'xhklinfo.fh' c ACCXFR - Accumulate F's for cross-validated free-r weighting c real hkls(nhklc,0:*) integer lftmp,lphitmp c integer i c do 100 i=1,nrflo if (hkls(lflag,i).lt.0.0) hkls(lffree,i)=hkls(lftmp,i) 100 continue return end c c ---------------------------------------------------------------------- c subroutine combxfr(hkls,lftmp,lphitmp) include 'dmmulti.fh' include 'io.fh' include 'cyclsync.fh' include 'xhklinfo.fh' include 'output.fh' c COMBXFR - phase combination by Sim, Nixon & North (197?), c Free-R multi stage cross-validating version. c real hkls(nhklc,0:*) integer lftmp,lphitmp c real rres(maxbin),refls(maxbin),reflsi(maxbin),swcen(maxbin), + sigmaq(maxbin),sumfq(maxbin),sfofc(maxbin), + sfomo(maxbin),sfomc(maxbin),sfomn(maxbin), + sphioc(maxbin),sphion(maxbin),sphicn(maxbin) logical centre c integer i,nrbin,nres,nrefli,ires,iresbin,iref real s,eps,fo,fc,wcen,x,phio,phic,fomo,fomc,hla,hlb,hlc,hld real sim,hlamod,hlbmod,phin,fomn,rdenom,rmod external iresbin,sim,rmod c c c INITIALISATION c ============== c initialise sin and cos tables call probsetup() c initialise arrays do 10 i=1,maxbin refls(i)=0.0 reflsi(i)=0.0 swcen(i)=0.0 sumfq(i)=0.0 sfomo(i)=0.0 sfomc(i)=0.0 sfomn(i)=0.0 sphioc(i)=0.0 sphion(i)=0.0 sphicn(i)=0.0 10 continue c set up resolution bins nrbin=25*max(ncfr,2) nrbin=min(nrbin,nrflo/50) call resbins(hkls,lffree,nrbin,maxbin,rres,nres) c c free-r stats and data scaling c call freerstat1(hkls,rres,nres) call scalefofc(hkls,rres,nres,lfo,lftmp,lffree,sfofc) c c clear new cols and count number of reflections included and in bins nrefli=0 do 20 i=1,nrflo if (hkls(lffree,i).lt.0.0) goto 20 c get data from array s=hkls(4,i) eps=abs(hkls(5,i)) centre=(hkls(5,i).lt.0.0) fo=hkls(lfo,i) fc=hkls(lffree,i) wcen=2.0 if (centre) wcen=1.0 c pick shell and include reflection ires=iresbin(s,rres,nres) nrefli=nrefli+1 refls(ires)=refls(ires)+1.0 swcen(ires)=swcen(ires)+wcen sumfq(ires)=sumfq(ires)+(fo**2+fc**2)*wcen/eps 20 continue c c CALCULATE SIGMAQ c ================ c iterative refinment by Nixon & North, Garib do 160 iref=1,10 c c first calculate sigmaq do 100 i=1,nres sigmaq(i)=sumfq(i)/max(swcen(i),1.0) sumfq(i)=0.0 100 continue c c now calc Sim weight: c w = tanh(X/2) ... centric c w = sim(X) ... non-centric c where X = 2 Fo Fc / eps * SIGMAq do 140 i=1,nrflo if (hkls(lffree,i).lt.0.0) goto 140 c get data from array s=hkls(4,i) eps=abs(hkls(5,i)) centre=(hkls(5,i).lt.0.0) fo=hkls(lfo,i) fc=hkls(lffree,i) wcen=2.0 if (centre) wcen=1.0 c pick shell and estimate fomc ires=iresbin(s,rres,nres) x=wcen*fo*fc/(eps*sigmaq(ires)) if (centre) then fomc=tanh(x) else fomc=sim(x) endif c update sigmaq estimate sumfq(ires)=sumfq(ires)+(fo**2+fc**2-2.0*fomc*fo*fc)*wcen/eps c next reflection 140 continue c next cycle of refinement 160 continue c c NOW DO THE PHASE RECOMBINATION c ============================== do 200 i=1,nrflo c NOTE: DO NOT INCLUDE FREE-R SET IN RECOMBINATION - c Only the F is available for free-r calc if (hkls(lflag,i).gt.0.0) then c get data from the array s=hkls(4,i) eps=abs(hkls(5,i)) centre=(hkls(5,i).lt.0.0) c obs data fo=hkls(lfo,i) phio=hkls(lphio,i) fomo=hkls(lfomo,i) hla=hkls(lhla,i) hlb=hkls(lhlb,i) hlc=hkls(lhlc,i) hld=hkls(lhld,i) c calc data fc=hkls(lftmp,i) phic=hkls(lphitmp,i) c weight for centrics wcen=2.0 if (centre) wcen=1.0 c res range ires=iresbin(s,rres,nres) reflsi(ires)=reflsi(ires)+1.0 c calc x x=wcen*fo*fc/(eps*sigmaq(ires)) if (centre) then fomc=tanh(x) else fomc=sim(x) endif c new h-l coeffs hlamod=x*cos(phic) hlbmod=x*sin(phic) hla=hla+hlamod hlb=hlb+hlbmod c calculate new phase and fom (suffix n=new for combined phase) if (centre) then phin=phio call probcen(hla,hlb,phin,fomn) else call probhl(hla,hlb,hlc,hld,phin,fomn) endif c store the new phase and fom hkls(lfdm,i)=fo hkls(lphidm,i)=phin hkls(lfomdm,i)=fomn c accumulate statistics sfomo(ires)=sfomo(ires)+fomo sfomc(ires)=sfomc(ires)+fomc sfomn(ires)=sfomn(ires)+fomn sphioc(ires)=sphioc(ires)+abs(rmod(phio-phic+pi,2.0*pi)-pi) sphion(ires)=sphion(ires)+abs(rmod(phio-phin+pi,2.0*pi)-pi) sphicn(ires)=sphicn(ires)+abs(rmod(phic-phin+pi,2.0*pi)-pi) c else zero excluded reflections else if (hkls(lflag,i).eq.0.0) then hkls(lfdm,i)=fo hkls(lphidm,i)=0.0 hkls(lfomdm,i)=0.0 c else leave free-r reflections alone endif c next reflection 200 continue c c now compile the statistics do 300 i=1,nres rdenom=max(reflsi(i),1.0) sfomo(i)=sfomo(i)/rdenom sfomc(i)=sfomc(i)/rdenom sfomn(i)=sfomn(i)/rdenom sphioc(i)=(sphioc(i)/dtor)/rdenom sphion(i)=(sphion(i)/dtor)/rdenom sphicn(i)=(sphicn(i)/dtor)/rdenom phcmbst(i,1)=1.0/rres(i)**2 phcmbst(i,2)=sphion(i) phcmbst(i,3)=sfomo(i) phcmbst(i,4)=sfomn(i) 300 continue nresio=nres c and display write (lpt,400)nrefli 400 format (' Phase combination statistics:'/ +' Number of reflections included in phase recombination =',i8// +' RMIN - RMAX S^2 NREFLS SIGMAQ SCALE ', +'MEAN FOM MEAN FOM MEAN FOM MEAN DPHI MEAN DPHI MEAN DPHI'/ +' (Fo/Fc) ', +' (obs.) (calc) (comb) |obs-cal| |obs-com| |cal-com|') write (lpt,410) + 100.0,rres(1),1.0/rres(1)**2,refls(1),sigmaq(1),sfofc(1), + sfomo(1),sfomc(1),sfomn(1),sphioc(1),sphion(1),sphicn(1) do 405 i=2,nres write (lpt,410) + rres(i-1),rres(i),1.0/rres(i)**2,refls(i),sigmaq(i),sfofc(i), + sfomo(i),sfomc(i),sfomn(i),sphioc(i),sphion(i),sphicn(i) 405 continue 410 format (1x,2f6.2,2x,f5.3,1x,f6.0,1x,f9.0,1x,f6.2,2x, + 3(f8.3,2x),3(f9.1,2x)) c return end c c ---------------------------------------------------------------------- c subroutine resbins(hkls,lflagx,nrange,maxranges,rres,nres) include 'dmmulti.fh' include 'xhklinfo.fh' c RESBINS - set up resolution bins subject to a minimum c number of reflections and minimum shell thickness c counting only reflections with the desired flag value c lflagxx = flag column, sflag = desired flag, c nrange = number in a range, maxranges = max. no. of ranges c real hkls(nhklc,0:*) real rres(*) integer lflagx,nrange,maxranges,nres c integer nhist parameter (nhist=1000) integer hist(nhist) c integer i,ires,ibin,nbin real s,ds,sbin,dsmin c c set up histogram ds=smax/(nhist-1) dsmin=smax/maxranges do 10 i=1,nhist hist(i)=0 10 continue c accumulate resolution histograms do 100 i=1,nrflo if (hkls(lflagx,i).le.0.0) goto 100 ires=int(hkls(4,i)/ds)+1 hist(ires)=hist(ires)+1 100 continue c now set resolution limits to accomodate minimum s and no. of refls ibin=0 nbin=0 sbin=0 do 200 i=1,nhist nbin=nbin+hist(i) sbin=sbin+ds if (nbin.gt.nrange.and.sbin.gt.dsmin) then ibin=ibin+1 rres(ibin)=1.0/sqrt(i*ds) nbin=0 sbin=0 endif 200 continue c tidy up last range nres=ibin rres(nres)=1.0/sqrt(smax) c return end c c ---------------------------------------------------------------------- c integer function iresbin(s,rres,nres) include 'dmmulti.fh' c IRESBIN - find the resolution bin for this reflection real s,rres(*) integer nres c integer i c do 10 i=1,nres-1 if (s.lt.1.0/rres(i)**2) goto 20 10 continue 20 iresbin=i return end c c ---------------------------------------------------------------------- c subroutine setresstat(s,rres,nres) include 'dmmulti.fh' c SETRESSTAT - set up weights for resolution averaged bins real s,rres(*) integer nres common /resdat/ires0,ires1,w0,w1 save /resdat/ c integer ires0,ires1,iresbin real s0,s1,w0,w1 external iresbin c ires1=iresbin(s,rres,nres) ires0=ires1-1 s0=0.0 if (ires0.gt.0) s0=1.0/rres(ires0)**2 s1=1.0/rres(ires1)**2 w0=(s1-s)/(s1-s0) w1=(s-s0)/(s1-s0) c return end c c ---------------------------------------------------------------------- c subroutine addresstat(stat,val) include 'dmmulti.fh' c ADDRESSTAT - accumulate resolution averaged statistic real stat(0:*),val common /resdat/ires0,ires1,w0,w1 save /resdat/ c integer ires0,ires1 real w0,w1 c stat(ires0)=stat(ires0)+w0*val stat(ires1)=stat(ires1)+w1*val c return end c c ---------------------------------------------------------------------- c real function resstat(stat) include 'dmmulti.fh' c GETRESSTAT - return a resolution averaged statistic real stat(0:*) common /resdat/ires0,ires1,w0,w1 save /resdat/ c integer ires0,ires1 real w0,w1 c resstat=w0*stat(ires0)+w1*stat(ires1) c return end c c ---------------------------------------------------------------------- c subroutine probsetup() include 'dmmulti.fh' c PROBSETUP - set up sin/cos lookup for phase integration common /trigwork/sintab(144),costab(144) save /trigwork/ integer i real x,sintab,costab c do 10 i=1,144 x=(i-1)*5*dtor sintab(i)=sin(x) costab(i)=cos(x) 10 continue return end c c ---------------------------------------------------------------------- c subroutine probhl(a,b,c,d,phi,fom) include 'dmmulti.fh' c PROBHL - integrate round phase circle using Henrickson-Lattman coeffs c taken from CCP4 code real a,b,c,d,phi,fom common /trigwork/sintab(144),costab(144) save /trigwork/ real prod(72),sintab,costab c integer i real explim,qmin,qmax,qscale,q,q0,pq,sind,cosd,sigd,fomc,foms,sim external sim parameter (explim=20.0) c if (c.eq.0.0.and.d.eq.0.0) then c a,b coefficients only fom=sim(sqrt(a**2+b**2)) phi=0.0 if (fom.gt.0.001) phi=atan2(b,a) c else c a,b,c,d coefficients qmin=1.e30 qmax=-1.e30 c loop round phase circle at 5 degree intervals, finding range of exponentials do 100 i=1,72 c---- Get log of calculated phase probability q=a*costab(i)+b*sintab(i)+c*costab(2*i)+d*sintab(2*i) qmin=min(qmin,q) qmax=max(qmax,q) prod(i)=q 100 continue c If the Q's go outside the range -EXPLIM to EXPLIM, find shift and c scale factors to get them in range qscale=1.0 if (qmax-qmin.gt.2.*explim) qscale=2.*explim/(qmax-qmin) q0=(qmax+qmin)/2 c Now run round phase circle again integrating probabilities cosd=0.0 sind=0.0 sigd=0.0 do 120 i=1,72 q=qscale*(prod(i)-q0) pq=exp(q) cosd=cosd+pq*costab(i) sind=sind+pq*sintab(i) sigd=sigd+pq 120 continue c Get new figure of merit, depending on whether centric or not fomc=cosd/sigd foms=sind/sigd fom=sqrt(fomc**2+foms**2) c and new phase phi=0.0 if (fom.gt.0.001) phi=atan2(foms,fomc) c endif c return end c c ---------------------------------------------------------------------- c subroutine probcen(a,b,phi,fom) include 'dmmulti.fh' c PROBCEN - Centric phase integration c IN: PHI = POSSIBLE VALUE FOR CENTRIC PHASE c taken from CCP4 code real a,b,phi,fom c common /trigwork/sintab(144),costab(144) save /trigwork/ real prod(72),sintab,costab c integer i,lmod real q,explim,eq1,eq2,cosd,sind,sigd,fomc,foms external lmod parameter (explim=20.0) c c extract initial phase i=nint(phi/dtor) i=lmod(i,180)/5+1 c c Get log of calculated phase probability q=a*costab(i)+b*sintab(i) q=min(q,explim) q=max(q,-explim) c c integrate probabilities eq1=exp(q) eq2=1.0/eq1 cosd=(eq1-eq2)*costab(i) sind=(eq1-eq2)*sintab(i) sigd=eq1+eq2 c c get new figure of merit fomc=cosd/sigd foms=sind/sigd fom=sqrt(fomc**2+foms**2) c and new phase phi=0.0 if (fom.gt.0.001) phi=atan2(foms,fomc) c return end c c ---------------------------------------------------------------------- c subroutine scalefofc(hkls,rres,nres,lf0,lf1,lf2,sfofc) include 'dmmulti.fh' include 'io.fh' include 'xhklinfo.fh' c SCALEFOFC - scale columns f1 and f2 to fit f1 to f0 c scale fc data to fo using weighted bins and interpolation c real hkls(nhklc,0:*) real rres(*),sfofc(*) integer nres,lf0,lf1,lf2 c real sfo(0:maxbin),sfc(0:maxbin),sw(0:maxbin) integer i real s,fo,fc,forms,fcrms,sc,resstat external resstat c do 90 i=0,nres sw(i)=0.0 sfo(i)=0.0 sfc(i)=0.0 90 continue c do 100 i=1,nrflo c included data only in calculating scale factors if (hkls(lflag,i).le.0.0) goto 100 s=hkls(4,i) fo=hkls(lf0,i) fc=hkls(lf1,i) call setresstat(s,rres,nres) call addresstat(sw,1.0) call addresstat(sfo,fo**2) call addresstat(sfc,fc**2) 100 continue c do 190 i=0,nres sfo(i)=sqrt(sfo(i)/max(sw(i),0.01)) sfc(i)=sqrt(sfc(i)/max(sw(i),0.01)) 190 continue c do 200 i=1,nrflo c scale all data by overall scale if (hkls(lflag,i).eq.0.0) goto 200 s=hkls(4,i) fo=hkls(lf0,i) fc=hkls(lf1,i) call setresstat(s,rres,nres) forms=resstat(sfo) fcrms=resstat(sfc) sc=forms/max(fcrms,0.1) c?? hkls(lf1,i) =sc*hkls(lf1,i) c?? hkls(lf2,i) =sc*hkls(lf2,i) 200 continue c c calculate scale for o/p statistics do 250 i=1,nres sfofc(i)=sfo(i)/max(sfc(i),0.01) 250 continue c return end c c ---------------------------------------------------------------------- c subroutine freerstat1(hkls,rres,nres) include 'dmmulti.fh' include 'io.fh' include 'xhklinfo.fh' include 'cyclsync.fh' include 'output.fh' c scale fc data to fo using weighted bins and interpolation c real hkls(nhklc,0:*) real rres(*) integer nres c real sfo(0:maxbin),sfc(0:maxbin),sw(0:maxbin) integer i real s,fo,fc,forms,fcrms,freefd,freefo,resstat external resstat c freefd=0.0 freefo=0.0 do 90 i=0,nres sw(i)=0.0 sfo(i)=0.0 sfc(i)=0.0 90 continue c do 100 i=1,nrflo c scale by free r set alone if (hkls(lffree,i).lt.0.0) goto 100 s=hkls(4,i) fo=hkls(lfo,i) fc=hkls(lffree,i) call setresstat(s,rres,nres) call addresstat(sw,1.0) call addresstat(sfo,fo**2) call addresstat(sfc,fc**2) 100 continue c do 190 i=0,nres sfo(i)=sqrt(sfo(i)/max(sw(i),0.01)) sfc(i)=sqrt(sfc(i)/max(sw(i),0.01)) 190 continue c do 200 i=1,nrflo c scale free r set alone if (hkls(lffree,i).lt.0.0) goto 200 s=hkls(4,i) fo=hkls(lfo,i) fc=hkls(lffree,i) call setresstat(s,rres,nres) forms=resstat(sfo) fcrms=resstat(sfc) c?? fc=fc*forms/max(fcrms,0.1) freefd=freefd+abs(fo-fc) freefo=freefo+fo 200 continue c freer(icycl)=1.0 if (freefo.gt.0.0) freer(icycl)=freefd/freefo write (lpt,250)thisxtal,freer(icycl) 250 format (/ + ' *************************'/ + ' Free density modification residual - XTAL: ',I2,f10.6/ + ' *************************'/) c return end c c ---------------------------------------------------------------------- c c c ---------------------------------------------------------------------- c subroutine addfpar(hkls,lf0,lf1,lphi1,lf2,lphi2) include 'dmmulti.fh' include 'io.fh' include 'xhklinfo.fh' c scale and add Fdm and Fpart for best agreement with Fo c real hkls(nhklc,0:*) integer lf0,lf1,lphi1,lf2,lphi2 c integer i,j,k,iter,niter real f1,f2,phi1,phi2,r,s,dxa,dxb,sca,scb,a1(maxbin),a2(maxbin) real phaseof,schess,resstat real x(2),dx(2),dr(2),h(2,2) real x01(2),x10(2),x11(2),dr1(2),h1(2,2),d,g1,g2 integer ires,nres,iresbin real rres(maxbin) real avres(maxbin),oldres,w0,w1,s0,s1,s2 complex fc,fc1,fc2,phase external iresbin,phase,phaseof,schess,resstat parameter (niter=10) c c call resbins(hkls,lflag,100,10,rres,nres) c x(1)=0.8 x(2)=1.2 d=0.05 do i=1,2 do j=1,2 do k=1,2 x01(k)=x(k) x10(k)=x(k) x11(k)=x(k) enddo x01(i)=x01(i)+d x10(j)=x10(j)+d x11(i)=x11(i)+d x11(j)=x11(j)+d g1=(schess(hkls,lf0,lf1,lphi1,lf2,lphi2,x01,dr,h,1,rres,nres)- + schess(hkls,lf0,lf1,lphi1,lf2,lphi2,x ,dr,h,1,rres,nres))/d g2=(schess(hkls,lf0,lf1,lphi1,lf2,lphi2,x11,dr,h,1,rres,nres)- + schess(hkls,lf0,lf1,lphi1,lf2,lphi2,x10,dr,h,1,rres,nres))/d dr1(i)=g1 h1(i,j)=(g2-g1)/d enddo enddo write (lpt,20)dr,dr1 write (lpt,20)h write (lpt,20)h1 20 format (2(2g16.5/)) c c now do refinement c oldres=1.0e6 do 110 ires=1,nres x(1)=1.0 x(2)=1.0 c do 100 iter=1,niter r=schess(hkls,lf0,lf1,lphi1,lf2,lphi2,x,dr,h,ires,rres,nres) call gausselim(h,dr,dx,2) dxa=sqrt(dx(1)**2+dx(2)**2) sca=min(0.25,dxa)/dxa x(1)=x(1)-sca*dx(1) x(2)=x(2)-sca*dx(2) c write (lpt,*)x 100 continue a1(ires)=x(1) a2(ires)=x(2) write (lpt,*)' ires,rres,a1,a2 ',ires,rres(ires),x c avres(ires)=1.0/sqrt(1.0/rres(ires)**2+1.0/oldres**2) oldres=rres(ires) 110 continue c do 300 i=1,nrflo if (hkls(lflag,i).le.0.0) goto 300 s=hkls(4,i) do 310 ires=2,nres if (s.lt.1.0/avres(ires)**2) goto 320 310 continue 320 s0=1.0/avres(ires-1)**2 s1=1.0/avres(ires)**2 w0=(s1-s)/(s1-s0) w1=(s-s0)/(s1-s0) s1=w0*a1(ires-1)+w1*a1(ires) s2=w0*a2(ires-1)+w1*a2(ires) f1=s1*hkls(lf1,i) f2=s2*hkls(lf2,i) phi1=hkls(lphi1,i) phi2=hkls(lphi2,i) fc1=f1*phase(phi1) fc2=f2*phase(phi2) fc=fc1+fc2 hkls(lf1,i)=abs(fc) hkls(lphi1,i)=phaseof(fc) 300 continue c return end c c ---------------------------------------------------------------------- c subroutine gausselim(a,b,x,n) include 'dmmulti.fh' c gaussian elimination for solution of ax=b c integer i,j,k,n real a(n,n),b(n),x(n) real s,pivot c do 100 i=1,n pivot=a(i,i) do 90 j=1,n if (j.ne.i) then s=a(j,i)/pivot do 80 k=i+1,n a(j,k)=a(j,k)-s*a(i,k) 80 continue b(j)=b(j)-s*b(i) endif 90 continue 100 continue c do 200 i=1,n x(i)=b(i)/a(i,i) 200 continue c return end c c ---------------------------------------------------------------------- c real function schess(hkls,lf0,lf1,lphi1,lf2,lphi2,x,dr,h, + ires,rres,nres) include 'dmmulti.fh' include 'xhklinfo.fh' c hessian of non-linear scaling system c real hkls(nhklc,0:*) integer lf0,lf1,lphi1,lf2,lphi2,ires,nres real x(2),dr(2),h(2,2),rres(*) c integer i,j,iresbin real a1,a2,s,fo,f1,f2,phi1,phi2,fc,phic real df,r,dfda1,dfda2,dfda11,dfda22,dfda12,cos12,phaseof complex cfc,phase external iresbin,phase,phaseof c a1=x(1) a2=x(2) c r=0.0 do 50 i=1,2 dr(i)=0.0 do 50 j=1,2 h(i,j)=0.0 50 continue c do 100 i=1,nrflo s=hkls(4,i) if (hkls(lflag,i).le.0.0.or.ires.ne.iresbin(s,rres,nres)) + goto 100 fo=hkls(lf0,i) f1=hkls(lf1,i) f2=hkls(lf2,i) phi1=hkls(lphi1,i) phi2=hkls(lphi2,i) cfc=a1*f1*phase(phi1)+a2*f2*phase(phi2) fc=abs(cfc) phic=phaseof(cfc) df=fo**2-fc**2 cos12=cos(phi1-phi2) c dfda1=2.0*a1*f1**2+2.0*a2*f1*f2*cos12 dfda2=2.0*a2*f2**2+2.0*a1*f1*f2*cos12 dfda11=2.0*f1**2 dfda22=2.0*f2**2 dfda12=2.0*f1*f2*cos12 c r=r+df**2 dr(1)=dr(1)-2.0*df*dfda1 dr(2)=dr(2)-2.0*df*dfda2 h(1,1)=h(1,1)-2.0*(df*dfda11-dfda1**2) h(2,2)=h(2,2)-2.0*(df*dfda22-dfda2**2) h(1,2)=h(1,2)-2.0*(df*dfda12-dfda1*dfda2) 100 continue c h(2,1)=h(1,2) schess=r c return end c c ---------------------------------------------------------------------- c