C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C c ...................................................................... program mama2ccp4 c mama2ccp4: convert a mama old or new format mask to ccp4 c array is dimensional to 20Mb, and elements are stored as c bytes using ccp4 byte handling conventions. If you can't fit c your mask in here, then you need a new program. c change as required. c c going the other way however elements are stored as reals. c implicit none integer maxsize,maxsec parameter (maxsize=5000000,maxsec=100000) integer mask(maxsize) c integer i,j,k,n,nn,ifail,lspgrp,nsymp,nsym,mamatype integer imask,icheck,imodew,imodec,nchitm,itest integer nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2,lu,lv,lw,iold,inew real cell(6),rot(4,4,96) character namspg*40,nampg*40,title*80,line*132,byteline*3 c integer ibeg(40),iend(40),ityp(40),idec(40) real fvalue(40) character*4 cvalue(40) c call ccpfyp CALL CCPRCS(6,'MAMA2CCP4','$Date: 2002/11/04 16:41:24 $') c c first fetch the spacegroup namspg=' ' call ugtenv('SPGRP',namspg) title = ' ' C now we try to work out what the MASKIN file is.... C The 53rd word of a CCP4 mask file is the string 'MAP ', so we C try to find it. imask = 8 imodew = 2 nchitm = 0 imodec = 0 byteline=' ' icheck = 0 call qopen(imask,'MASKIN','READONLY') call qmode(imask,imodew,nchitm) call qseek(imask,1,53,0) call qmode(imask,imodec,nchitm) call qreadc(imask,byteline,icheck) call qclose(imask) if (byteline.eq.'MAP') goto 1000 c c MAMA -> CCP4 --------------------------------------------------------- c write (*,80) 80 format (/,' CONVERTING MAMA -> CCP4',/) call ccpupc(namspg) if (index('PABCFIR',namspg(1:1)).eq.0) namspg='P1' lspgrp=0 call msymlb(8,lspgrp,namspg,nampg,nsymp,nsym,rot) write (*,90)lspgrp,namspg 90 format (' Space group no. ',i4,' Name ',a) c c READ IN A MAMA MASK -------------------------------------------------- c i=0 ifail=0 call ccpdpn(7,'MASKIN','READONLY','F',i,ifail) c c ***HEADER*** c c check for mask format read (7,91)line if (line(1:8).eq.'COMPRESS') then c PARSE COMPRESSED FORMAT MAMA MASK HEADER mamatype=3 itest = 0 read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 3) then call gtpint(1,nu1,n,ityp,fvalue) call gtpint(2,nv1,n,ityp,fvalue) call gtpint(3,nw1,n,ityp,fvalue) if (n .ge. 4) title=line(ibeg(4):len(line)) itest = itest + 1 endif read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 3) then call gtpint(1,lu,n,ityp,fvalue) call gtpint(2,lv,n,ityp,fvalue) call gtpint(3,lw,n,ityp,fvalue) itest = itest + 1 endif read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 3) then call gtpint(1,nu,n,ityp,fvalue) call gtpint(2,nv,n,ityp,fvalue) call gtpint(3,nw,n,ityp,fvalue) itest = itest + 1 endif read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 6) then call gtnrea(1,6,cell,n,ityp,fvalue) itest = itest + 1 endif if (itest .ne. 4) then write(6,92) 92 format(/,' *** Mask header corrupted or new O-mask format',/) call ccperr(1,' *** ERROR: program terminated ***') endif else if (line(1:8).eq.'NEW_MASK') then c PARSE NEW FORMAT MAMA MASK HEADER mamatype=2 itest=0 read (7,91)title 91 format (a) 1110 read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (cvalue(1).eq.'ORIG') then if (n.ge.4) then call gtpint(2,nu1,n,ityp,fvalue) call gtpint(3,nv1,n,ityp,fvalue) call gtpint(4,nw1,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'EXTE') then if (n.ge.4) then call gtpint(2,lu,n,ityp,fvalue) call gtpint(3,lv,n,ityp,fvalue) call gtpint(4,lw,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'GRID') then if (n.ge.4) then call gtpint(2,nu,n,ityp,fvalue) call gtpint(3,nv,n,ityp,fvalue) call gtpint(4,nw,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'CELL') then if (n.ge.6) then call gtnrea(2,6,cell,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).ne.'MAP') goto 1110 if (itest .ne. 4) then write(6,92) call ccperr(1,' *** ERROR: program terminated ***') endif c REALLY NEW MASK FORMAT CAN BE COMPRESSED OR UNCOMPRESSED else if (line(1:11) .eq. '.MASK_INPUT') then itest = 0 mamatype = 0 c read (7,91)title title = ' ' 1210 read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (cvalue(1).eq.'ORIG') then if (n.ge.4) then call gtpint(2,nu1,n,ityp,fvalue) call gtpint(3,nv1,n,ityp,fvalue) call gtpint(4,nw1,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'EXTE') then if (n.ge.4) then call gtpint(2,lu,n,ityp,fvalue) call gtpint(3,lv,n,ityp,fvalue) call gtpint(4,lw,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'GRID') then if (n.ge.4) then call gtpint(2,nu,n,ityp,fvalue) call gtpint(3,nv,n,ityp,fvalue) call gtpint(4,nw,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'CELL') then if (n.ge.6) then call gtnrea(2,6,cell,n,ityp,fvalue) itest = itest + 1 endif endif if (cvalue(1).eq.'COMP') then mamatype = 3 else if (cvalue(1).eq.'EXPL') then mamatype = 1 else goto 1210 endif if (itest.ne.4 .or. mamatype.eq.0) then write(6,92) call ccperr(1,' *** ERROR: program terminated ***') endif else c PARSE OLD FORMAT MAMA MASK HEADER mamatype=1 itest = 0 n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 3) then call gtpint(1,nu1,n,ityp,fvalue) call gtpint(2,nv1,n,ityp,fvalue) call gtpint(3,nw1,n,ityp,fvalue) title=line(ibeg(4):len(line)) itest = itest + 1 endif read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 3) then call gtpint(1,lu,n,ityp,fvalue) call gtpint(2,lv,n,ityp,fvalue) call gtpint(3,lw,n,ityp,fvalue) itest = itest + 1 endif read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 3) then call gtpint(1,nu,n,ityp,fvalue) call gtpint(1,nv,n,ityp,fvalue) call gtpint(1,nw,n,ityp,fvalue) itest = itest + 1 endif read (7,91)line n=-40 call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,n) if (n .ge. 6) then call gtnrea(1,6,cell,n,ityp,fvalue) itest = itest + 1 endif if (itest .ne. 4) then write(6,92) call ccperr(1,' *** ERROR: program terminated ***') endif endif c c done header, now read the rest in any format c ***MASK*** c if (mamatype.eq.3) then c COMPRESSED FORMAT do 140 i=1,lu*lv*lw call ccpstb(0,mask,i) 140 continue n=0 150 nn=-40 read (7,110,end=180,err=180)line call parse(line,ibeg,iend,ityp,fvalue,cvalue,idec,nn) call gtpint(1,j,nn,ityp,fvalue) call gtpint(2,k,nn,ityp,fvalue) do 160 i=j,k call ccpstb(1,mask,i) n=n+1 160 continue goto 150 c 180 continue else c OLD OR NEW FORMAT j=0 n=0 c get a line 100 line=' ' read (7,110,end=200,err=200)line 110 format (a) c process the line do 120 i=1,len(line) k=index('01',line(i:i))-1 if (k.ge.0) then j=j+1 call ccpstb(k,mask,j) n=n+k endif 120 continue c goto 100 c 200 continue if (j.ne.lu*lv*lw) call ccperr(1,' error reading mama mask') c endif c c done mask read c close (unit=7,status='KEEP') c write (*,205)real(n)/real(nu*nv*nw) 205 format (' Fraction of cell covered by mask = ',f8.4) c c WRITE OUT A CCP4 MASK ------------------------------------------------ c nu2=nu1+lu-1 nv2=nv1+lv-1 nw2=nw1+lw-1 write (*,10)nu1,nu2,nv1,nv2,nw1,nw2 10 format (' nu1 : nu2, nv1 : nv2, nw1 : nw2 '/6i6) call ccpmskout('MASKOUT',title,mask,lspgrp,cell, + nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2) CALL CCPERR(0,'Normal termination') c c MAMA -> CCP4 --------------------------------------------------------- c 1000 write (*,1080) 1080 format (/,' CONVERTING CCP4 -> MAMA',/) c c read a ccp4 mask c call ccpmaphead + ('MASKIN ',lspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2) call ccpmskin + ('MASKIN ',title,mask,nu1,nv1,nw1,nu2,nv2,nw2) c c write a mama mask c lu=nu2-nu1+1 lv=nv2-nv1+1 lw=nw2-nw1+1 c i=0 ifail=0 call ccpdpn(7,'MASKOUT','NEW','F',i,ifail) write (7,2110) write (7,2120)nu1,nv1,nw1,title write (7,2130)lu,lv,lw write (7,2140)nu,nv,nw write (7,2150)cell iold=0 do 1100 i=1,lu*lv*lw call ccpgtb(inew,mask,i) if (iold.eq.0.and.inew.eq.1) j=i if (iold.eq.1.and.inew.eq.0) write (7,2160)j,i-1 iold=inew 1100 continue if (iold.eq.1) write (7,2160)j,i-1 close (unit=7,status='KEEP') 2110 format('COMPRESSED_MASK') 2120 format(3i5,1x,a) 2130 format(3i5) 2140 format(3i5) 2150 format(6f10.4) 2160 format(2i9) c CALL CCPERR(0,'Normal termination') end c c ---------------------------------------------------------------------- c c c ---------------------------------------------------------------------- c subroutine ccpmaphead + (name,nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2) c CCPMAPhead - get limits of mask from file implicit none c character name*8 integer nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2 real cell(6) c integer jfms(3),juvw(3),mxyz(3),m1(3),m2(3),nsec,mode integer i,ifail,iprint real rmin,rmax,rmean,rrms character title*80 c c read the file headers ifail=0 iprint=1 call mrdhds(8,name,title,nsec,jfms,mxyz,m1(3),m1(1),m2(1), + m1(2),m2(2),cell,nspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint) call mrclos(8) m2(3)=m1(3)+nsec-1 do 110 i=1,3 juvw(jfms(i))=i 110 continue nu1=m1(juvw(1)) nu2=m2(juvw(1)) nv1=m1(juvw(2)) nv2=m2(juvw(2)) nw1=m1(juvw(3)) nw2=m2(juvw(3)) nu=mxyz(1) nv=mxyz(2) nw=mxyz(3) c return end c c ---------------------------------------------------------------------- c subroutine ccpmskin + (name,title,map,nu1,nv1,nw1,nu2,nv2,nw2) c CCPMSKIN - read a ccp4 logical map and store in xyz order implicit none c integer nu1,nv1,nw1,nu2,nv2,nw2 integer map(*) character name*(*),title*(*) c integer maxsec parameter (maxsec=50000) real lsec(0:maxsec) c integer ifast,imedm,islow,lfast,lmedm,lslow,ierr,ifail,iprint integer m1(3),m2(3),ifms(3),jfms(3),juvw(3),mxyz(3) integer i,k,iu,iv,iw,lu,lv,lw integer nsec,mode,mspgrp real cell(6),rmin,rmax,rmean,rrms c c c Note: juvw convert from fast/med/slow to u/v/w c jfms convert from u/v/w to fast/med/slow c c now open map header and read map, re-ordering it as necessary c ifail=0 iprint=0 call mrdhds(8,name,title,nsec,jfms,mxyz,m1(3),m1(1),m2(1), + m1(2),m2(2),cell,mspgrp,mode,rmin,rmax,rmean,rrms,ifail,iprint) m2(3)=m1(3)+nsec-1 do 110 i=1,3 juvw(jfms(i))=i 110 continue c check we got the right header info: if (nu1.ne.m1(juvw(1)).or.nu2.ne.m2(juvw(1)).or. + nv1.ne.m1(juvw(2)).or.nv2.ne.m2(juvw(2)).or. + nw1.ne.m1(juvw(3)).or.nw2.ne.m2(juvw(3))) + call ccperr(1,'ccpmskin - map grid sizes are inconsistent') c find out the map grid dimensions lfast=m2(1)-m1(1)+1 lmedm=m2(2)-m1(2)+1 lslow=m2(3)-m1(3)+1 lu=nu2-nu1+1 lv=nv2-nv1+1 lw=nw2-nw1+1 c if (lfast*lmedm.gt.4*maxsec) + call ccperr(1,' ccpmskin - Map section > lsec: recompile') c now read the map in the order it is on file do 200 islow=0,lslow-1 c get a section ierr=0 call mgulpr(8,lsec,ierr) if (ierr.ne.0) call ccperr(1,' ccpmskin - ccp4 read error') c now sort into uvw map ifms(3)=islow do 190 imedm=0,lmedm-1 ifms(2)=imedm do 190 ifast=0,lfast-1 ifms(1)=ifast iu=ifms(juvw(1)) iv=ifms(juvw(2)) iw=ifms(juvw(3)) k=nint(lsec(ifast+lfast*imedm)) call ccpstb(k,map,iu+lu*(iv+lv*iw)+1) 190 continue 200 continue c close the map file call mrclos(8) c write (6,910) 910 format (/' MAP/MASK READ SUCCESSFUL'//) c return end c c ---------------------------------------------------------------------- c subroutine ccpmskout + (name,title,map,nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2) c CCPMSKOUT implicit none c integer nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2,map(*) real cell(6) character name*(*),title*(*) c integer maxsec parameter (maxsec=50000) integer lsec(0:maxsec) c integer ifast,imedm,islow,lfast,lmedm,lslow integer m1(3),m2(3),ifms(3),jfms(3),juvw(3),mxyz(3) integer i,k,iu,iv,iw,lu,lv,lw integer nsec,mode,nsym,nsymp real rsym(4,4,192) character*10 namspg,nampg c c c spacegroups with yxz=1 and zxy=2 axis ordering integer axis(230) data axis/2,2,2,2,1,1,1,1,1,2,1,1,1,1,1,2,2,2,1,2,2,1,2,207*1/ c c Note: juvw convert from fast/med/slow to u/v/w c jfms convert from u/v/w to fast/med/slow c c now open map header and read map, re-ordering it as necessary c if (axis(mod(nspgrp,1000)).eq.1) then jfms(1)=2 jfms(2)=1 jfms(3)=3 else jfms(1)=3 jfms(2)=1 jfms(3)=2 endif c do 110 i=1,3 juvw(jfms(i))=i 110 continue c lu=nu2-nu1+1 lv=nv2-nv1+1 lw=nw2-nw1+1 mxyz(1)=nu mxyz(2)=nv mxyz(3)=nw m1(juvw(1))=nu1 m2(juvw(1))=nu2 m1(juvw(2))=nv1 m2(juvw(2))=nv2 m1(juvw(3))=nw1 m2(juvw(3))=nw2 nsec=m2(3)-m1(3)+1 mode=0 c call msymlb(8,nspgrp,namspg,nampg,nsym,nsymp,rsym) c call mwrhdl(8,name,title,nsec,jfms,mxyz, + m1(3),m1(1),m2(1),m1(2),m2(2),cell,nspgrp,mode) call msywrt(8,nsym,rsym) c 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 c if (lfast*lmedm.gt.4*maxsec) + call ccperr(1,' ccpmapin - Mask section > lsec: recompile') c now write the map in the new order do 200 islow=0,lslow-1 c sort onto section ifms(3)=islow do 190 imedm=0,lmedm-1 ifms(2)=imedm do 190 ifast=0,lfast-1 ifms(1)=ifast iu=ifms(juvw(1)) iv=ifms(juvw(2)) iw=ifms(juvw(3)) call ccpgtb(k,map,iu+lu*(iv+lv*iw)+1) call ccpstb(k,lsec,ifast+lfast*imedm+1) 190 continue c write the section call mspew(8,lsec) 200 continue c close the map file call mwclose(8) c return end c c ---------------------------------------------------------------------- c