c c ---------------------------------------------------------------------- c subroutine domnin(dommap,filnam) c DOMNIN - read domain mask include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' c real dommap(mu1:mu2,mv1:mv2,mw1:mw2) character filnam*(*),title*80 c integer iu,iv,iw,ju,jv,jw c c first read in the mask call dommskin(dommap,filnam) c c Now expand the mask one pixel c c Note that we need to average to just outside the mask so that values c just outside the mask will be available for averaging back into the c unit cells. c c mark inside and outside as +/-1000 do 100 iw=mw1,mw2 do 100 iv=mv1,mv2 do 100 iu=mu1,mu2 if (dommap(iu,iv,iw).gt.0.5) then dommap(iu,iv,iw)=1000.0 else dommap(iu,iv,iw)=-1000.0 endif 100 continue c now set all points within 1 pixel of an inside point to 0 do 200 iw=mw1+1,mw2-1 do 200 iv=mv1+1,mv2-1 do 200 iu=mu1+1,mu2-1 if (dommap(iu,iv,iw).gt.500.0) then do 150 jw=iw-1,iw+1 do 150 jv=iv-1,iv+1 do 150 ju=iu-1,iu+1 dommap(ju,jv,jw)=max(dommap(ju,jv,jw),0.0) 150 continue endif 200 continue c and set the inside to 0 do 300 iw=mw1,mw2 do 300 iv=mv1,mv2 do 300 iu=mu1,mu2 if (dommap(iu,iv,iw).gt.500.0) dommap(iu,iv,iw)=0.0 300 continue c return end c c ---------------------------------------------------------------------- c subroutine dommskin(dommap,filnam) c DOMMSKIN - read domain mask include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' c real dommap(mu1:mu2,mv1:mv2,mw1:mw2) character filnam*(*),title*80 c integer iu,iv,iw,iu1,iv1,iw1,iu2,iv2,iw2,lfast,lmedm,lslow integer jfms(3),juvw(3),mxyz(3),m1(3),m2(3),msec,nspgrp,mode integer msku,mskv,mskw,nn,mm integer i,ifast,imedm,islow,ifms(3),ifail,iprint,ierr real cell(6),rmin,rmax,rmean,rrms integer lsec(maxsec) c integer lint external lint c c clear mask mm=0 nn=0 do 100 iw=mw1,mw2 do 100 iv=mv1,mv2 do 100 iu=mu1,mu2 dommap(iu,iv,iw)=0.0 100 continue c c read header ifail=0 iprint=1 call mrdhds(lmsks,filnam,title,msec,jfms,mxyz,m1(3),m1(1),m2(1), + m1(2),m2(2),cell,nspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint) m2(3)=m1(3)+msec-1 do 110 i=1,3 juvw(jfms(i))=i 110 continue c find out the mask grid dimensions lfast=m2(1)-m1(1)+1 lmedm=m2(2)-m1(2)+1 lslow=m2(3)-m1(3)+1 msku=mxyz(1) mskv=mxyz(2) mskw=mxyz(3) c check format if (lfast*lmedm.gt.4*maxsec) + call ccperr(1,' domnin - Mask section > lsec: recompile') if (mode.ne.0) + call ccperr(1,' domnin - Mask file is not logical*1') c c now read the map in the order it is on the file do 200 islow=0,lslow-1 c get a section ierr=0 call mgulp(lmsks,lsec,ierr) if (ierr.ne.0) call ccperr(1,' domnin - ccp4 read error') c now sort into uvw mask and extrapolate onto program grid ifms(3)=islow+m1(3) do 190 imedm=0,lmedm-1 ifms(2)=imedm+m1(2) do 190 ifast=0,lfast-1 ifms(1)=ifast+m1(1) call ccpgtb(i,lsec,ifast+lfast*imedm+1) if (i.ne.0) mm=mm+1 iu=ifms(juvw(1)) iv=ifms(juvw(2)) iw=ifms(juvw(3)) iu1=lint((real(iu)-0.50)*mu/msku) iu2=lint((real(iu)+0.50)*mu/msku) iv1=lint((real(iv)-0.50)*mv/mskv) iv2=lint((real(iv)+0.50)*mv/mskv) iw1=lint((real(iw)-0.50)*mw/mskw) iw2=lint((real(iw)+0.50)*mw/mskw) c we must fill every program grid point covered by this mask point do 180 iw=iw1+1,iw2 do 180 iv=iv1+1,iv2 do 180 iu=iu1+1,iu2 if (iu.ge.mu1.and.iv.ge.mv1.and.iw.ge.mw1.and. + iu.le.mu2.and.iv.le.mv2.and.iw.le.mw2) then if (i.ne.0) dommap(iu,iv,iw)=1.0 endif 180 continue 190 continue 200 continue c close the map file call mrclos(lmsks) c c print some stats do 250 iw=mw1,mw2 do 250 iv=mv1,mv2 do 250 iu=mu1,mu2 if (dommap(iu,iv,iw).gt.0.5) nn=nn+1 250 continue if (iprint.ne.0) write (lpt,260) + mm,real(mm)/real(msku*mskv*mskw),nn,real(nn)/real(mu*mv*mw) 260 format( + ' Mask grid : Grid points set- ',i7,' Cell fraction- ',f6.3/ + ' Program grid: Grid points set- ',i7,' Cell fraction- ',f6.3/) c c print a mask sec do 350 iv=mv1,mv2 do 350 iu=mu1,mu2 lsec((iu-mu1)+(mu2-mu1+1)*(iv-mv1)+1)= + nint(dommap(iu,iv,(mw1+mw2)/2)) 350 continue call masksec(lsec,mu2-mu1+1,mv2-mv1+1) c return end c c ---------------------------------------------------------------------- c subroutine domzero(dommap) c DOMZERO - clear a domain include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' c real dommap(mu1:mu2,mv1:mv2,mw1:mw2) c integer iu,iv,iw c c clear the domain map leaving only the flags -1000,0,1000 do 100 iw=mw1,mw2 do 100 iv=mv1,mv2 do 100 iu=mu1,mu2 dommap(iu,iv,iw)=real(1000*nint(dommap(iu,iv,iw)/1000.0)) 100 continue c return end c c ---------------------------------------------------------------------- c c Routines for averaging from a unit cell into a dommain store... c c ---------------------------------------------------------------------- c subroutine averagefromcell(hkls,map,dommul,domoff,xtlwt) c AVERAGEFROMCELL - average from unit cells into domain masks include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' include 'xavginfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c real hkls(nhklc,0:*),map(0:nu-1,0:nv-1,0:nw-1),xtlwt real dommul(*) integer domoff(*) c integer xncs,domn real dtox(4,4),xtod(4,4),temp(4,4) c c c loop through each ncs in turn ... c do 500 xncs=1,nncs c . average to which domain? domn=srcdom(xncs) call getdomn(domn) c . mat = (xfrac)(rot)(dorth) c . xtod goes from cell to domain c . dtox goes from domain to cell call symmatmul(rsymncs(1,1,xncs),dq,temp) call symmatmul(qt,temp,dtox) call invsym(dtox,xtod) c . do the average (this in subroutine so the domain map can be indexed) call doaveragefrom(map,dommul(domoff(domn)),dtox,xtod,xtlwt) 500 continue c return end c c ---------------------------------------------------------------------- c subroutine doaveragefrom(map,dommap,dtox,xtod,wt) c DOAVERAGEFROM - single step average from unit cell into domain mask include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c real map(0:nu-1,0:nv-1,0:nw-1),dommap(mu1:mu2,mv1:mv2,mw1:mw2) real dtox(4,4),xtod(4,4),wt c integer iu,iv,iw real du,dv,dw,xu,xv,xw c real bsplinterp,gridinterp external bsplinterp,gridinterp c do 500 iw=mw1,mw2 do 500 iv=mv1,mv2 do 500 iu=mu1,mu2 if (dommap(iu,iv,iw).gt.-500.0) then du=real(iu)/mu dv=real(iv)/mv dw=real(iw)/mw call symuvw(dtox,du,dv,dw,xu,xv,xw) dommap(iu,iv,iw)=dommap(iu,iv,iw)+ + wt*bsplinterp(map,xu*nu,xv*nv,xw*nw,nu,nv,nw,q) endif 500 continue c return end c c + wt*gridinterp(map,xu*nu,xv*nv,xw*nw,0,0,0,nu-1,nv-1,nw-1,q) c c ---------------------------------------------------------------------- c c Routines for expanding from a domain store into a unit cell... c c ---------------------------------------------------------------------- c subroutine averagetocell(hkls,map,dommul,domoff,reseff) c AVERAGETOCELL - average from unit cells into domain masks include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' include 'xavginfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c real hkls(nhklc,0:*),map(0:nu-1,0:nv-1,0:nw-1),reseff real dommul(*) integer domoff(*) c integer xncs,domn real dtox(4,4),xtod(4,4),temp(4,4),wt,correl c real stdres external stdres c c loop through each ncs in turn ... c do 500 xncs=1,nncs c . average to which domain? domn=srcdom(xncs) call getdomn(domn) c . get normalisation for this dommain wt=stdres(1.0/reseff**2) c . mat = (xfrac)(rot)(dorth) c . xtod goes from cell to domain c . dtox goes from domain to cell call symmatmul(rsymncs(1,1,xncs),dq,temp) call symmatmul(qt,temp,dtox) call invsym(dtox,xtod) c . do the average (this in subroutine so the domain map can be indexed) call doaverageto(map,dommul(domoff(domn)),dtox,xtod,wt,correl) write (lpt,490)xncs,domn,correl 490 format (' Averaging operation: ',i2,' Averaging from domain:', + i2,' Correlation: ',f6.3) 500 continue c return end c c ---------------------------------------------------------------------- c subroutine doaverageto(map,dommap,dtox,xtod,stdprt,correl) c DOAVERAGETO - single step expand from domain mask into unit cell include 'dmmulti.fh' include 'xdominfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c real map(0:nu-1,0:nv-1,0:nw-1),dommap(mu1:mu2,mv1:mv2,mw1:mw2) real dtox(4,4),xtod(4,4),stdprt,correl c integer iu,iv,iw,ju,jv,jw,ku,kv,kw,isym integer u1,v1,w1,u2,v2,w2 real du,dv,dw,xu,xv,xw,xu1,xv1,xw1 real rx,ry,sn,sx,sy,sx2,sy2,sxy,domave,domstd,mapave,mapstd c integer lmod real gridinterp logical gridinside external lmod,gridinterp,gridinside c c first make a box in the xtal coords covering the domain mask c and do stats on domain density c sn =0.0 sx =0.0 sx2=0.0 u1=100000 v1=100000 w1=100000 u2=-100000 v2=-100000 w2=-100000 c do 200 iw=mw1,mw2 dw=real(iw)/mw do 200 iv=mv1,mv2 dv=real(iv)/mv do 200 iu=mu1,mu2 du=real(iu)/mu c? if (gridinside(dommap, c? + real(iu),real(iv),real(iw),mu1,mv1,mw1,mu2,mv2,mw2)) then if (dommap(iu,iv,iw).gt.-500.0) then call symuvw(dtox,du,dv,dw,xu,xv,xw) u1=min(u1,nint(xu*nu-0.5)) u2=max(u2,nint(xu*nu+0.5)) v1=min(v1,nint(xv*nv-0.5)) v2=max(v2,nint(xv*nv+0.5)) w1=min(w1,nint(xw*nw-0.5)) w2=max(w2,nint(xw*nw+0.5)) sn =sn +1.0 sx =sx +dommap(iu,iv,iw) sx2=sx2+dommap(iu,iv,iw)**2 endif 200 continue c mapave=sx/sn mapstd=sqrt(sn*sx2-sx**2)/sn domave=aveprt domstd=stdprt c c now build the map density c sn=0.0 sx=0.0 sy=0.0 sx2=0.0 sy2=0.0 sxy=0.0 c do 600 iw=w1,w2 jw=lmod(iw,nw) xw=real(iw)/nw do 600 iv=v1,v2 jv=lmod(iv,nv) xv=real(iv)/nv do 600 iu=u1,u2 ju=lmod(iu,nu) xu=real(iu)/nu c . already checked if this bit is in mask c . get interpolant and do correlations c . note: we must do mask check again because a point in the map may c . correspond to two points in the mask for wrapped masks call symuvw(xtod,xu,xv,xw,du,dv,dw) if (gridinside + (dommap,du*mu,dv*mv,dw*mw,mu1,mv1,mw1,mu2,mv2,mw2)) then rx=map(ju,jv,jw) ry=domave+(domstd/mapstd)*(gridinterp + (dommap,du*mu,dv*mv,dw*mw,mu1,mv1,mw1,mu2,mv2,mw2,dq)-mapave) sn=sn+1.0 sx=sx+rx sy=sy+ry sx2=sx2+rx**2 sy2=sy2+ry**2 sxy=sxy+rx*ry c . and expand by crystallographic symmetry into map do 550 isym=1,nsym call symuvw(rsym(1,1,isym),xu,xv,xw,xu1,xv1,xw1) ku=lmod(nint(xu1*nu),nu) kv=lmod(nint(xv1*nv),nv) kw=lmod(nint(xw1*nw),nw) map(ku,kv,kw)=ry 550 continue endif 600 continue c correl=(sn*sxy-sx*sy)/ + sqrt(max((sn*sx2-sx**2)*(sn*sy2-sy**2),1.0e-12)) c return end c c ---------------------------------------------------------------------- c c Routines for refining averaging operators from domain back to map c c ---------------------------------------------------------------------- c subroutine refinefromcell(map,hklmul,dommul,xtlwt,lfa,lfb) c REFINEFROMCELL - refine from unit cells into domain masks include 'dmmulti.fh' include 'io.fh' include 'xtalinfo.fh' include 'xdominfo.fh' include 'xavginfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c real map(*) real hklmul(*),dommul(*),xtlwt(*) integer lfa,lfb c integer xtal,domn,xncs,srcxtl,srcncs real temp(4,4),xtod(4,4),dtox(4,4) c c srcxtl=1 srcncs=1 c c first make copy the best copy of the domain into the domain map c call getxtal(srcxtl) cv loop over all domains do 100 domn=1,ndomn call getdomn(domn) c . mat = (xfrac)(rot)(dorth) c . xtod goes from cell to domain c . dtox goes from domain to cell call symmatmul(rsymncs(1,1,srcncs),dq,temp) call symmatmul(qt,temp,dtox) call invsym(dtox,xtod) c . do the average (this in subroutine so the domain map can be indexed) call fftlistmap(map,hklmul(hkloff(srcxtl)),lfa,lfb) call doaveragefrom(map,dommul(domoff(domn)),dtox,xtod, + xtlwt(srcxtl)) 100 continue cv c c loop through each ncs in turn ... c do 500 xtal=1,nxtal c . get the map for this xtal call getxtal(xtal) if (refdr.gt.0.0.and.xtlwt(xtal).gt.0.05) then call fftlistmap(map,hklmul(hkloff(xtal)),lfa,lfb) do 490 xncs=1,nncs c .. refine with which domain? domn=srcdom(xncs) call getdomn(domn) c .. do the refine (this in subroutine so the domain map can be indexed) if (xtal.ne.srcxtl.or.xncs.ne.srcncs) then write (lpt,460)xncs,domn 460 format (' Averaging operation number ',i3,' from domain ',i1) call dorefinefrom(map,dommul(domoff(domn)),rsymncs(1,1,xncs)) else write (lpt,470)xncs,domn 470 format (' Source averaging operator ',i3,' from domain ',i1) endif 490 continue c . and store the refined matrices call ccpmvr(xtaldat(lsrcncs+16,xtal),rsymncs,16*nncs) endif 500 continue c return end c c ---------------------------------------------------------------------- c subroutine dorefinefrom(map,dommap,rmat) c DOREFINEFROM - single step refine from unit cell into domain mask include 'dmmulti.fh' include 'io.fh' include 'xdominfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c real map(0:nu-1,0:nv-1,0:nw-1),dommap(mu1:mu2,mv1:mv2,mw1:mw2) real rmat(4,4) c integer i,j,k,iu,iv,iw,iparam real nr,a,b,x,y,y2,z0,z(2),step,dp(6),dp1(6),dy(6) real temp1(4,4),temp2(4,4),dtox(4,4),rncs(4,4,2),drdx(4,4,6) real du,dv,dw,xd,yd,zd,xc,yc,zc c c data ((temp2(i,j),j=1,4),i=1,4)/ c + 0.99947,-0.02798,-0.01660,0.5, c + 0.02777,0.99953,-0.01287,0.0, c + 0.01695,0.01241,0.99978,-0.3, c + 0.0,0.0,0.0,1.0/ c call symcpy(rmat,temp1) c call symmatmul(temp2,temp1,rmat) c first get molecule centre in orthogonal coordinates du=0.0 dv=0.0 dw=0.0 nr=0.0 do 50 iw=mw1,mw2 do 50 iv=mv1,mv2 do 50 iu=mu1,mu2 if (dommap(iu,iv,iw).gt.-50.0) then du=du+real(iu)/mu dv=dv+real(iv)/mv dw=dw+real(iw)/mw nr=nr+1.0 endif 50 continue call symuvw(dq,du/nr,dv/nr,dw/nr,xd,yd,zd) call symuvw(rmat,xd,yd,zd,xc,yc,zc) c make partial derivative matrices (jacobian) call makemat(temp1,0,0.0,xc,yc,zc) do 70 k=1,6 call makemat(temp2,k,0.0001,xc,yc,zc) do 70 j=1,4 do 70 i=1,4 drdx(i,j,k)=(temp2(i,j)-temp1(i,j))/0.0001 70 continue c make coordinate shifts do 100 i=1,3 dp(i)=0.05 dp(i+3)=0.5 100 continue c step=1.0 i=0 c 200 i=i+1 c FIRST WE WORK OUT THE STEP DIRECTION c calculate the correlation, and its gradient with each of the 6 params call ncscorgrd(map,dommap,rmat,drdx,y,dy) y2=sqrt(dy(1)**2+dy(2)**2+dy(3)**2+dy(4)**2+dy(5)**2+dy(6)**2) c???????????????????????????????????????????????????????????????????????? c write (*,*)'molc',xc,yc,zc c write (*,301)'rmat',((rmat(k,j),j=1,4),k=1,4) c write (*,*)y c do 300 iparam=1,6 c a=0.005 c if (iparam.ge.4) a=0.01 c call makemat(temp1,iparam,a,xc,yc,zc) c call symmatmul(temp1,rmat,rncs(1,1,1)) c call ncscorgrd(map,dommap,rncs(1,1,1),drdx,z0,temp2) c dp1(iparam)=(y-z0)/a c write (*,*)z0 c 300 continue c 301 format (1x,a4,1x,4f12.6/3(6x,4f12.6/)) c write (*,'(6f12.6)'),dy,dp1 c???????????????????????????????????????????????????????????????????????? c c NOW WE USE QUADRATIC TRIAL TO GET STEP LENGTH c make matrices for shifts of lengths 1 and 2 call symcpy(rmat,rncs(1,1,1)) call symcpy(rmat,rncs(1,1,2)) do 340 iparam=1,6 dp1(iparam)=-step*dy(iparam)*dp(iparam)/y2 c 1/2 x input step call makemat(temp1,iparam,0.5*dp1(iparam),xc,yc,zc) call symcpy(rncs(1,1,1),temp2) call symmatmul(temp1,temp2,rncs(1,1,1)) c 1 x input step call makemat(temp1,iparam,dp1(iparam),xc,yc,zc) call symcpy(rncs(1,1,2),temp2) call symmatmul(temp1,temp2,rncs(1,1,2)) 340 continue c now generate correlations z0=y call ncscorgrd(map,dommap,rncs(1,1,1),drdx,z(1),dy) call ncscorgrd(map,dommap,rncs(1,1,2),drdx,z(2),dy) c b=4.0*z(1)-3.0*z0-z(2) a=2.0*z(2)+2.0*z0-4.0*z(1) x=0.0 if (a.lt.0.0) then x=min(max(-0.5*b/a,0.0),1.0) else if (z(2).gt.z0) x=1.0 endif c NOW DO THE SHIFT do 360 iparam=1,6 call makemat(temp1,iparam,x*dp1(iparam),xc,yc,zc) call symcpy(rmat,temp2) call symmatmul(temp1,temp2,rmat) 360 continue write (lpt,810)z0,(x*dp1(j)/dtor,j=1,3),(x*dp1(j),j=4,6) c guess a step size for next step step=step*sqrt(x) c next refinement cycle if (step.gt.0.25.and.i.lt.8) goto 200 c 810 format (' Correlation: ',f6.3,' Step: ',6f8.2) c return end c c + wt*gridinterp(map,xu*nu,xv*nv,xw*nw,0,0,0,nu-1,nv-1,nw-1,q) c c ---------------------------------------------------------------------- c subroutine ncscorgrd(map,dommap,rmat,rderiv,correl,dcorrel) include 'dmmulti.fh' include 'xdominfo.fh' include 'xtlparam.fh' include 'xcrystal.fh' c ncscorgrd - perform a ncs correlation using the given matrix c and calculate correlation gradiants along directions given by c partial derivative tensors rderiv. c real map(0:nu-1,0:nv-1,0:nw-1),dommap(mu1:mu2,mv1:mv2,mw1:mw2) real rmat(4,4),rderiv(4,4,6),correl,dcorrel(6) c integer i,iu,iv,iw real u,v,w,u1,v1,w1,x1,y1,z1,dx1,dy1,dz1,g(3),dtox(4,4),temp(4,4) real rx,ry,sn,sx,sy,sx2,sy2,sxy,dr(6),gy(6),gy2(6),gxy(6) real bsplinterp c c c mat = (xfrac)(rot)(dorth) c dtox goes from domain to cell call symmatmul(rmat,dq,temp) call symmatmul(qt,temp,dtox) c sn=0.0 sx=0.0 sy=0.0 sx2=0.0 sy2=0.0 sxy=0.0 do 50 i=1,6 gy(i)=0.0 gy2(i)=0.0 gxy(i)=0.0 50 continue c c accumulate correlation statistics c note: rmat maps fractional domain to fractional crystal do 100 iw=mw1,mw2 w=real(iw)/mw do 100 iv=mv1,mv2 v=real(iv)/mv do 100 iu=mu1,mu2 u=real(iu)/mu if (dommap(iu,iv,iw).gt.-500.0) then call symuvw(dtox,u,v,w,u1,v1,w1) rx=dommap(iu,iv,iw) call bsplintgrd(ry,g,map,nu*u1,nv*v1,nw*w1,nu,nv,nw,q) sn=sn+1.0 sx=sx+rx sx2=sx2+rx**2 sy=sy+ry sy2=sy2+ry**2 sxy=sxy+rx*ry c calculate gradient along parameter directions c note: derivs and grads are on orthogonal coordinates call symuvw(q,u1,v1,w1,x1,y1,z1) do 90 i=1,6 call symuvw(rderiv(1,1,i),x1,y1,z1,dx1,dy1,dz1) dr(i)=dx1*g(1)+dy1*g(2)+dz1*g(3) gy(i)=gy(i)+dr(i) gy2(i)=gy2(i)+2.0*ry*dr(i) gxy(i)=gxy(i)+rx*dr(i) 90 continue endif 100 continue c c and calculate correlation and gradients correl=(sn*sxy-sx*sy)/sqrt((sn*sx2-sx**2)*(sn*sy2-sy**2)) do 150 i=1,6 dcorrel(i)=(sn*gxy(i)-sx*gy(i))/ + sqrt((sn*sx2-sx**2)*(sn*sy2-sy**2)) c FUDGE: assumes variance is constant (tests OK). correct code is: c + ( (sn*gxy(i)-sx*gy(i))*(sn*sx2-sx**2)*(sn*sy2-sy**2) - c + 0.5*(sn*sxy-sx*sy)*(sn*sx2-sx**2)*(sn*gy2(i)-2.0*sy*gy(i)) )/ c + ((sn*sx2-sx**2)*(sn*sy2-sy**2))**1.5 150 continue c c write (*,*)correl return end c c this code checks the gradients but not rderiv:- c real x2,y2,z2,u2,v2,w2,dr1(6) c x2=x1+0.001* c +(rderiv(1,1,i)*x1+rderiv(1,2,i)*y1+rderiv(1,3,i)*z1+rderiv(1,4,i)) c y2=y1+0.001* c +(rderiv(2,1,i)*x1+rderiv(2,2,i)*y1+rderiv(2,3,i)*z1+rderiv(2,4,i)) c z2=z1+0.001* c +(rderiv(3,1,i)*x1+rderiv(3,2,i)*y1+rderiv(3,3,i)*z1+rderiv(3,4,i)) c call symuvw(qt,x2,y2,z2,u2,v2,w2) c dr1(i)=1000*(bsplinterp(map,nu*u2,nv*v2,nw*w2,nu,nv,nw,q)-ry) c write (*,'(6f12.6)')dr c write (*,'(6f12.6)')dr1 c c ---------------------------------------------------------------------- c subroutine makemat(mat,iparam,dp,cx,cy,cz) include 'dmmulti.fh' c makemat - make a matrix to shift parameter i c integer iparam real mat(4,4),cx,cy,cz,dp c real temp(4,4),dphi integer i,j c dphi=dp c do 110 j=1,4 do 100 i=1,4 mat(i,j)=0.0 100 continue mat(j,j)=1.0 110 continue c if (iparam.le.0) return c if (iparam.le.3) then i=mod(iparam,3)+1 j=mod(i,3)+1 c rotation about 0 0 0 mat(i,i)=cos(dphi) mat(j,j)=cos(dphi) mat(i,j)=sin(dphi) mat(j,i)=-sin(dphi) c adjust for rotation about cx cy cz mat(1,4)=cx-(mat(1,1)*cx+mat(1,2)*cy+mat(1,3)*cz) mat(2,4)=cy-(mat(2,1)*cx+mat(2,2)*cy+mat(2,3)*cz) mat(3,4)=cz-(mat(3,1)*cx+mat(3,2)*cy+mat(3,3)*cz) else c translation along x,y,z mat(iparam-3,4)=dp endif c return end c c ---------------------------------------------------------------------- c c Interpolation routines of different sorts... c c ---------------------------------------------------------------------- c subroutine bsplinewt(hkls,lfa,lfb,lfwa,lfwb) include 'dmmulti.fh' include 'xcrystal.fh' include 'xhklinfo.fh' include 'xtlparam.fh' c BSPLINEWT - apply reciprocal space weights to correct for B-spline integer lfa,lfb,lfwa,lfwb real hkls(nhklc,0:*) c integer i,ih,ik,il real rh,rk,rl,x1,x2,x3,wt complex phase external phase c hkls(lfwa,0)=hkls(lfa,0) hkls(lfwb,0)=hkls(lfb,0) do 100 i=1,nrflo ih=(nint(hkls(1,i))) ik=(nint(hkls(2,i))) il=(nint(hkls(3,i))) rh=ih+(qt(1,2)*ik+qt(1,3)*il)/qt(1,1) rk=ik+qt(2,3)*il/qt(2,2) rl=il x1=max(pi*abs(rh)/real(nu),0.001) x2=max(pi*abs(rk)/real(nv),0.001) x3=max(pi*abs(rl)/real(nw),0.001) wt=(sin(x1)*sin(x2)*sin(x3)/(x1*x2*x3))**3 hkls(lfwa,i)=hkls(lfa,i)/wt hkls(lfwb,i)=hkls(lfb,i)/wt 100 continue c return end c c ---------------------------------------------------------------------- c c c ---------------------------------------------------------------------- c real function bsplinterp(rho,u,v,w,nu,nv,nw,q) include 'dmmulti.fh' c bsplinterp - Interpolate a density point by bspline convolution c integer nu,nv,nw real rho(0:nu-1,0:nv-1,0:nw-1),q(4,4),u,v,w c integer i,j,k,l,m,i1,j1,k1,iu(3),iv(3),iw(3),lmod real vw,uvw,bw,bvw,rw,sw,swrho,du(3),dv(3),dw(3),bsplconv external lmod,bsplconv c i=nint(u) j=nint(v) k=nint(w) c c set up tables of map indices and coordinate differences do 100 l=1,3 m=l-2 i1=i+m iu(l)=lmod(i1,nu) du(l)=i1-u j1=j+m iv(l)=lmod(j1,nv) dv(l)=j1-v k1=k+m iw(l)=lmod(k1,nw) dw(l)=k1-w 100 continue c c now perform summation - approx 300 flops sw=0.0 swrho=0.0 do 400 k1=1,3 k=iw(k1) w=dw(k1) bw=bsplconv(w) vw=q(2,3)*w/q(2,2) do 400 j1=1,3 j=iv(j1) v=dv(j1) bvw=bsplconv(v+vw)*bw uvw=(q(1,2)*v+q(1,3)*w)/q(1,1) do 400 i1=1,3 i=iu(i1) u=du(i1) rw=bsplconv(u+uvw)*bvw c sw=sw+rw swrho=swrho+rw*rho(i,j,k) 400 continue c bsplinterp=swrho/sw return end c c ---------------------------------------------------------------------- c subroutine bsplintgrd(b,g,rho,u,v,w,nu,nv,nw,q) include 'dmmulti.fh' c bsplintgrd - Interpolate a density point and gradient by bspline convolution c integer nu,nv,nw real b,g(3),rho(0:nu-1,0:nv-1,0:nw-1),q(4,4),u,v,w c integer i,j,k,l,m,i1,j1,k1,iu(3),iv(3),iw(3),lmod real vw,uvw,bu,bv,bw,gu,gv,gw,rw,rx,ry,rz real sw,swrho,sxrho,syrho,szrho,du(3),dv(3),dw(3) external lmod c i=nint(u) j=nint(v) k=nint(w) c c set up tables of map indices and coordinate differences do 100 l=1,3 g(l)=0.0 m=l-2 i1=i+m iu(l)=lmod(i1,nu) du(l)=i1-u j1=j+m iv(l)=lmod(j1,nv) dv(l)=j1-v k1=k+m iw(l)=lmod(k1,nw) dw(l)=k1-w 100 continue c c now perform summation - approx 300 flops sw=0.0 swrho=0.0 sxrho=0.0 syrho=0.0 szrho=0.0 do 400 k1=1,3 k=iw(k1) w=dw(k1) call bsplcgrd(w,bw,gw) vw=q(2,3)*w/q(2,2) do 400 j1=1,3 j=iv(j1) v=dv(j1) call bsplcgrd(v+vw,bv,gv) uvw=(q(1,2)*v+q(1,3)*w)/q(1,1) do 400 i1=1,3 i=iu(i1) u=du(i1) call bsplcgrd(u+uvw,bu,gu) c rw=bu*bv*bw rx=gu*bv*bw ry=bu*gv*bw rz=bu*bv*gw c sw=sw+rw swrho=swrho+rw*rho(i,j,k) sxrho=sxrho+rx*rho(i,j,k) syrho=syrho+ry*rho(i,j,k) szrho=szrho+rz*rho(i,j,k) 400 continue c g(1)=(sxrho/sw)*(nu/q(1,1)) g(2)=(syrho/sw)*(nv/q(2,2)) g(3)=(szrho/sw)*(nw/q(3,3)) b=swrho/sw return end c c ---------------------------------------------------------------------- c real function bsplconv(dr) include 'dmmulti.fh' c bsplconv - quandratic b-spline convolution function real dr,x x=abs(dr) if (x.ge.1.5) then bsplconv=0.0 else if (x.ge.0.5) then bsplconv=0.5*(x-1.5)**2 else bsplconv=-x**2+0.75 endif return end c c ---------------------------------------------------------------------- c subroutine bsplcgrd(dr,b,grad) include 'dmmulti.fh' c bsplcgrd - gradient of quandratic b-spline convolution function real dr,b,grad if (dr.ge.1.5) then b=0.0 grad=0.0 else if (dr.ge.0.5) then grad=dr-1.5 b=0.5*grad**2 else if (dr.gt.-0.5) then b=-dr**2+0.75 grad=-2.0*dr else if (dr.gt.-1.5) then grad=dr+1.5 b=0.5*grad**2 else b=0.0 grad=0.0 endif return end c c ---------------------------------------------------------------------- c logical function + gridinside(r,u,v,w,nu1,nv1,nw1,nu2,nv2,nw2) c GRIDINSIDE - check if all interpolation source points are inside mask include 'dmmulti.fh' c integer nu1,nv1,nw1,nu2,nv2,nw2 real r(nu1:nu2,nv1:nv2,nw1:nw2),u,v,w c integer iu0,iv0,iw0,iu1,iv1,iw1,iu2,iv2,iw2,iu3,iv3,iw3 real sum c integer lint external lint c iu0=lint(u)-1 iv0=lint(v)-1 iw0=lint(w)-1 iu1=iu0+1 iv1=iv0+1 iw1=iw0+1 iu2=iu1+1 iv2=iv1+1 iw2=iw1+1 iu3=iu2+1 iv3=iv2+1 iw3=iw2+1 c if (iu0.ge.nu1.and.iv0.ge.nv1.and.iw0.ge.nw1.and. + iu3.le.nu2.and.iv3.le.nv2.and.iw3.le.nw2) then c middle box + u squares + v squares + w squares sum= + r(iu1,iv1,iw1)+r(iu1,iv1,iw2)+r(iu1,iv2,iw1)+r(iu1,iv2,iw2)+ + r(iu2,iv1,iw1)+r(iu2,iv1,iw2)+r(iu2,iv2,iw1)+r(iu2,iv2,iw2)+ + r(iu0,iv1,iw1)+r(iu0,iv1,iw2)+r(iu0,iv2,iw1)+r(iu0,iv2,iw2)+ + r(iu3,iv1,iw1)+r(iu3,iv1,iw2)+r(iu3,iv2,iw1)+r(iu3,iv2,iw2)+ + r(iu1,iv0,iw1)+r(iu1,iv0,iw2)+r(iu2,iv0,iw1)+r(iu2,iv0,iw2)+ + r(iu1,iv3,iw1)+r(iu1,iv3,iw2)+r(iu2,iv3,iw1)+r(iu2,iv3,iw2)+ + r(iu1,iv1,iw0)+r(iu1,iv2,iw0)+r(iu2,iv1,iw0)+r(iu2,iv2,iw0)+ + r(iu1,iv1,iw3)+r(iu1,iv2,iw3)+r(iu2,iv1,iw3)+r(iu2,iv2,iw3) else sum=-1000.0 endif c gridinside=(sum.gt.-500.0) c return end c c ---------------------------------------------------------------------- c real function + gridinterp(rho,x,y,z,nu1,nv1,nw1,nu2,nv2,nw2,q) c GRIDINTERP - non-orthogonal grid interpolation include 'dmmulti.fh' c integer nu1,nv1,nw1,nu2,nv2,nw2 real rho(0:nu2-nu1,0:nv2-nv1,0:nw2-nw1),q(4,4),x,y,z c integer nu,nv,nw,iu,iv,iw,iu1,iv1,iw1,iu2,iv2,iw2 real u,v,w,du1,dv1,dw1,du2,dv2,dw2 real r111,rx111,ry111,rz111,r112,rx112,ry112,rz112 real r121,rx121,ry121,rz121,r122,rx122,ry122,rz122 real r211,rx211,ry211,rz211,r212,rx212,ry212,rz212 real r221,rx221,ry221,rz221,r222,rx222,ry222,rz222 real ra11,ra12,ra21,ra22,rb1,rb2,rzb1,rzb2 real rya11,rya12,rya21,rya22,rza11,rza12,rza21,rza22 c integer lmod,lint real rcubiccub,rcubiclin external rcubiccub,rcubiclin,lmod,lint c nu=nu2-nu1+1 nv=nv2-nv1+1 nw=nw2-nw1+1 c u=x-nu1 iu=lint(u) iu1=lmod(iu,nu) iu2=lmod(iu+1,nu) du1=u-iu v=y-nv1 iv=lint(v) iv1=lmod(iv,nv) iv2=lmod(iv+1,nv) dv1=v-iv w=z-nw1 iw=lint(w) iw1=lmod(iw,nw) iw2=lmod(iw+1,nw) dw1=w-iw c call rcubicgrad(rho,iu1,iv1,iw1,nu,nv,nw,r111,rx111,ry111,rz111) call rcubicgrad(rho,iu1,iv1,iw2,nu,nv,nw,r112,rx112,ry112,rz112) call rcubicgrad(rho,iu1,iv2,iw1,nu,nv,nw,r121,rx121,ry121,rz121) call rcubicgrad(rho,iu1,iv2,iw2,nu,nv,nw,r122,rx122,ry122,rz122) call rcubicgrad(rho,iu2,iv1,iw1,nu,nv,nw,r211,rx211,ry211,rz211) call rcubicgrad(rho,iu2,iv1,iw2,nu,nv,nw,r212,rx212,ry212,rz212) call rcubicgrad(rho,iu2,iv2,iw1,nu,nv,nw,r221,rx221,ry221,rz221) call rcubicgrad(rho,iu2,iv2,iw2,nu,nv,nw,r222,rx222,ry222,rz222) c ra11=rcubiccub(du1,r111,r211,rx111,rx211) ra12=rcubiccub(du1,r112,r212,rx112,rx212) ra21=rcubiccub(du1,r121,r221,rx121,rx221) ra22=rcubiccub(du1,r122,r222,rx122,rx222) rya11=rcubiclin(du1,ry111,ry211) rya12=rcubiclin(du1,ry112,ry212) rya21=rcubiclin(du1,ry121,ry221) rya22=rcubiclin(du1,ry122,ry222) rza11=rcubiclin(du1,rz111,rz211) rza12=rcubiclin(du1,rz112,rz212) rza21=rcubiclin(du1,rz121,rz221) rza22=rcubiclin(du1,rz122,rz222) c rb1=rcubiccub(dv1,ra11,ra21,rya11,rya21) rb2=rcubiccub(dv1,ra12,ra22,rya12,rya22) rzb1=rcubiclin(dv1,rza11,rza21) rzb2=rcubiclin(dv1,rza12,rza22) c gridinterp=rcubiccub(dw1,rb1,rb2,rzb1,rzb2) return end c subroutine rcubicgrad(rho,iu,iv,iw,nu,nv,nw,r,rx,ry,rz) include 'dmmulti.fh' integer iu,iv,iw,nu,nv,nw,lmod real rho(0:nu-1,0:nv-1,0:nw-1),rx,ry,rz,r external lmod r=rho(iu,iv,iw) rx=0.5*(rho(lmod(iu+1,nu),iv,iw)-rho(lmod(iu-1,nu),iv,iw)) ry=0.5*(rho(iu,lmod(iv+1,nv),iw)-rho(iu,lmod(iv-1,nv),iw)) rz=0.5*(rho(iu,iv,lmod(iw+1,nw))-rho(iu,iv,lmod(iw-1,nw))) return end c real function rcubiccub(x,y1,y2,dy1,dy2) include 'dmmulti.fh' real x,y1,y2,dy1,dy2,a,b,c,d a=(dy1+dy2)+2.0*(y1-y2) b=3.0*(y2-y1)-2.0*dy1-dy2 c=dy1 d=y1 rcubiccub=d+x*(c+x*(b+x*a)) return end c real function rcubiclin(x,y1,y2) include 'dmmulti.fh' real x,y1,y2 rcubiclin=(1.0-x)*y1+x*y2 return end c c ---------------------------------------------------------------------- c real function stdres(s) c STDRES - return expected protein standard deviation as a function of res. include 'dmmulti.fh' c real s c integer ndatares,i real res(50),std(50),s0,s1,w0,w1 data ndatares/20/ data res/ 1.0e6,10.00, 8.00, 7.00, 6.00, 5.00, 4.50, 4.00, 3.60, + 3.30, 3.00, 2.80, 2.60, 2.40, 2.20, 2.00, 1.80, 1.60, 1.50, + 1.0,30*0.0/ data std/ 0.000,0.076,0.107,0.121,0.143,0.200,0.248,0.316,0.375, + 0.427,0.476,0.512,0.552,0.599,0.654,0.718,0.783,0.850,0.890, + 1.0,30*0.0/ c do 100 i=2,ndatares s1=1.0/res(i)**2 if (s.lt.s1) goto 110 100 continue 110 s0=1.0/res(i-1)**2 w0=(s1-s)/(s1-s0) w1=1.0-w0 stdres=w0*std(i-1)+w1*std(i) c return end c c ---------------------------------------------------------------------- c