c
c ----------------------------------------------------------------------
c
c The routines CCP4MAPHEAD, CCP4MAPIN and CCP4MAPOUT were written
c by Kevin Cowtan. As of version 3.5 of CCP4 these routines will
c be placed in the library maplib.f, however they are also examples
c of using the maplib routines.
c
c This program will read a complete map into an array, automatically
c sorted as X,Y,Z. The map array is then manipulated in the subroutine
c 'mapmanipulate' (in this example the values in the map are inverted,
c but users can easily alter/add to the code to perform more complicated
c manipulations). The modified map is then written out to a second file.
c
c It should be run as:
c
c map_rw MAPIN foo1.map MAPOUT foo2.map
c
c Nb there are no keywords; see the parsing templates for examples
c of adding keywords.
c
c (Thanks to Ho-Leung Ng for fixing some problems with this template!)
c
program map_rw
c
implicit none
c
c ..Parameters ..
INTEGER MAXMAP
PARAMETER (MAXMAP=1000000)
c
c ..Variables to store the map info
INTEGER nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2,
+ mapsize
REAL cell(6)
REAL my_map(maxmap)
CHARACTER*120 title
c
c CCP4 initialisation
call ccpfyp
c
c Get header info for map
c
write(6,FMT=9100)
write(6,fmt=9000)
write(6,FMT=9100)
9000 format('***** RETRIEVING HEADER INFORMATION... *****')
c
call ccp4maphead('MAPIN',nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,
+ nv2,nw2)
c
c Check that the map will fit into the array
c
mapsize = (nu2-nu1)*(nv2-nv1)*(nw2-nw1)
c
if (mapsize.gt.maxmap) then
write(6,fmt='(a,i10,a)') 'Warning: increase parameter MAXMAP to'
+ ,mapsize,' and recompile'
call ccperr (1, 'Map array too small')
end if
c
c Write some information from the header to the standard output
c
write(6,FMT=9100)
write(6,FMT=9010)nspgrp
write(6,FMT=9020)cell
write(6,FMT=9030)nu1,nu2,nv1,nv2,nw1,nw2
c
9010 format(/,'Spacegroup number ',I4)
9020 format(/,'Cell dimensions:',3F8.3,' A',/,
+ 'Cell angles :',3F8.1,' deg.')
9030 format(/,'Size of map: nu1 nu2 nv1 nv2 nw1 nw2',
+ /,15X,6I5,/)
c
c Read in the map
c
write(6,FMT=9100)
write(6,FMT=9040)
write(6,FMT=9100)
9040 format('***** READING IN THE MAP ... *****')
c
call ccp4mapin('MAPIN',title,my_map,nu1,nv1,nw1,nu2,nv2,nw2)
c
c Write more information to standard output
c
write(6,FMT=9100)
write(6,FMT=9050)title
9050 format(/,'Map title:',//,A120)
c
c ********************************************************
c **** Put code in MAPMANIPULATE to use/alter the map ****
c ********************************************************
c
call mapmanipulate(my_map,nu1,nv1,nw1,nu2,nv2,nw2)
c
c Write out the modified map
c
write(6,FMT=9100)
write(6,FMT=9060)
write(6,FMT=9100)
9060 format('***** WRITING OUT THE NEW MAP ... *****')
9100 format('*************************************************')
c
call ccp4mapout('MAPOUT',title,my_map,nspgrp,cell,
+ nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2)
c
c End the program cleanly
c
call ccperr(0,'Normal termination')
end
c
c ----------------------------------------------------------------------
c
c_BEGIN_MAPMANIPULATE
c
subroutine mapmanipulate(map,nu1,nv1,nw1,nu2,nv2,nw2)
c =====================================================
c
c---- The map held in map(nu1:nu2,nv1:nv2,nw1:nw2) can be manipulated
c within this subroutine. In this example the values in the map
c are simply scaled by a factor of -1.
c
c_END_MAPMANIPULATE
c
implicit none
c
c ..Subroutine arguments
integer nu1,nv1,nw1,nu2,nv2,nw2
real map(nu1:nu2,nv1:nv2,nw1:nw2)
c
c ..Local variables
integer ic,jc,kc
c
c *********************************************
c **** Put code here to manipulate the map ****
c *********************************************
c
c In this example the map is inverted, i.e. scaled by a factor of -1.
c
do 10 ic = nu1,nu2
do 20 jc = nv1,nv2
do 30 kc = nw1,nw2
c
map(ic,jc,kc) = -map(ic,jc,kc)
c
30 continue
20 continue
10 continue
c
return
c
end
c
c
c ----------------------------------------------------------------------
c
c_BEGIN_CCP4MAPHEAD
c
subroutine ccp4maphead(name,nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,
+ nv2,nw2)
c ================================================================
c
c CCP4MAPhead - get limits of mask from file
c
c---- Function: read header information from a map file. It is used to
c get the map limits before calling ccp4mapin.
c
c Call: CALL ccp4maphead(name,nspgrp,cell,nu,nv,nw,nu1,nv1,nw1,nu2,
c + nv2,nw2)
c
c---- Note that This differs from MRDHDS (which it calls) in that the
c file is not left open but is closed (by a call to MRCLOS) before
c the subroutine terminates. Note also that the map limits are
c returned in x,y,z order rather than in fast, medium, slow order.
c
c
c---- Arguments:
c =========
c
c name (I) Logical file name (type character*8) e.g.'MAPIN'
c nspgrp (O) Space group number (integer)
c cell (O) 6 word array for cell dimensions in Angstroms and
c degrees (real)
c nu (O) Sampling interval along whole cell on X (integer)
c nv (O) Sampling interval along whole cell on Y (integer)
c nw (O) Sampling interval along whole cell on Z (integer)
c nu1 (O) Start of map on X axis, in grid units (integer)
c nv1 (O) Start of map on Y axis, in grid units (integer)
c nw1 (O) Start of map on Z axis, in grid units (integer)
c nu2 (O) End of map on X axis (integer)
c nv2 (O) End of map on Y axis (integer)
c nw2 (O) End of map on Z axis (integer)
c
c_END_CCP4MAPHEAD
c
implicit none
c
character name*(*)
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(7,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(7)
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
c_BEGIN_CCP4MAPIN
c
subroutine ccp4mapin(name,title,map,nu1,nv1,nw1,nu2,nv2,nw2)
c ============================================================
c
c CCP4MAPIN - read a ccp4 logical map and store in xyz order
c
c---- Function: read whole map into an array and store in x,y,z order.
c The map limits are required as input to dimension the array holding
c the map, and can be obtained with a call to the subroutine
c ccp4maphead.
c
c Call: CALL ccp4mapin (name,title,map,nu1,nv1,nw1,nu2,nv2,nw2)
c
c---- ccp4mapin is a "wrapper" subroutine which utilises calls to the
c following maplib routines: MRDHDS, MGULPR, MRCLOS.
c
c
c---- Arguments:
c ==========
c
c name (O) logical file name (type character) e.g. 'MAPIN'
c title (O) Map title (type character)
c map (O) Real array of dimension (nu1:nu2,nv1:nv2,nw1:nw2)
c which stores the map which is read in
c nu1 (I) Start of map on X axis, in grid units (integer)
c nv1 (I) Start of map on Y axis, in grid units (integer)
c nw1 (I) Start of map on Z axis, in grid units (integer)
c nu2 (I) End of map on X axis (integer)
c nv2 (I) End of map on Y axis (integer)
c nw2 (I) End of map on Z axis (integer)
c
c_END_CCP4MAPIN
c
implicit none
c
integer nu1,nv1,nw1,nu2,nv2,nw2
real map(nu1:nu2,nv1:nv2,nw1:nw2)
character name*(*),title*(*)
c
real lsec(0:100000)
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,iu,iv,iw
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(7,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,'ccp4mapin - mask 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
c
if (lfast*lmedm.gt.100000)
+ call ccperr(1,' ccp4mapin - 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(7,lsec,ierr)
if (ierr.ne.0) call ccperr(1,' ccp4mapin - ccp4 read error')
c now sort into uvw map
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)
iu=ifms(juvw(1))
iv=ifms(juvw(2))
iw=ifms(juvw(3))
map(iu,iv,iw)=lsec(ifast+lfast*imedm)
190 continue
200 continue
c close the map file
call mrclos(7)
c
write (*,910)
910 format (/' MAP/MASK READ SUCCESSFUL'//)
c
return
end
c
c ----------------------------------------------------------------------
c
c_BEGIN_CCP4MAPOUT
c
subroutine ccp4mapout(name,title,map,nspgrp,cell,nu,nv,nw,nu1,
+ nv1,nw1,nu2,nv2,nw2)
c ==============================================================
c
c CCP4MAPOUT - write a ccp4 logical map in xyz order
c
c---- Function: Write out a whole map in x,y,z order.
c
c Call: CALL ccpmap4out(name,title,map,nspgrp,cell,nu,nv,nw,nu1,
c + nv1,nw1,nu2,nv2,nw2)
c
c---- ccpmap4out is a "wrapper" routine which utilises calls to the
c following maplib subroutines: MWRHDL, MSYWRT, MSPEW, MWCLOSE. There
c is also a call to the symlib routine MSYMLB.
c
c
c---- Arguments:
c ==========
c
c name (I) Logical file name (type character) e.g.'MAPIN'
c title (I) Map title (type character)
c map (I) Real array of dimension (nu1:nu2,nv1:nv2,nw1:nw2)
c which stores the map being written out
c nspgrp (I) Space group number (integer)
c cell (I) 6 word array for cell dimensions in Angstroms and degrees (real)
c nu (I) Sampling interval along whole cell on X (integer)
c nv (I) Sampling interval along whole cell on Y (integer)
c nw (I) Sampling interval along whole cell on Z (integer)
c nu1 (I) Start of map on X axis, in grid units (integer)
c nv1 (I) Start of map on Y axis, in grid units (integer)
c nw1 (I) Start of map on Z axis, in grid units (integer)
c nu2 (I) End of map on X axis (integer)
c nv2 (I) End of map on Y axis (integer)
c nw2 (I) End of map on Z axis (integer)
c
c_END_CCP4MAPOUT
c
implicit none
c
integer nspgrp,nu,nv,nw,nu1,nv1,nw1,nu2,nv2,nw2
real cell(6),map(nu1:nu2,nv1:nv2,nw1:nw2)
character name*(*),title*(*)
c
real lsec(0:100000)
c
integer ifast,imedm,islow,lfast,lmedm,lslow
integer m1(3),m2(3),ifms(3),jfms(3),juvw(3),mxyz(3)
integer i,iu,iv,iw
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
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=2
c
call msymlb(7,nspgrp,namspg,nampg,nsym,nsymp,rsym)
c
call mwrhdl(7,name,title,nsec,jfms,mxyz,
+ m1(3),m1(1),m2(1),m1(2),m2(2),cell,nspgrp,mode)
call msywrt(7,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.100000)
+ call ccperr(1,' ccp4mapout - 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+m1(3)
do 190 imedm=0,lmedm-1
ifms(2)=imedm+m1(2)
do 190 ifast=0,lfast-1
ifms(1)=ifast+m1(1)
iu=ifms(juvw(1))
iv=ifms(juvw(2))
iw=ifms(juvw(3))
lsec(ifast+lfast*imedm)=map(iu,iv,iw)
190 continue
c write the section
call mspew(7,lsec)
200 continue
c close the map file
call mwclose(7)
c
return
end
c
c ----------------------------------------------------------------------
c