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   7 Sep 1992  Phil Evans     MRC LMB, Cambridge
C      Added CHAIN option to match SFALL
C      Some minor tidying
C
C===================================================================
C     PROGRAM TO CALCULATE THE OVERLAP OF TWO MAPS
C     the two maps could be on different scale and averages !
C     this is to calculate the CORRELATION COEFFICIENT of two maps
C     ******************************************************************
C     **                        ave(X*Y) - ave(X)*ave(Y)
C     * corr. coeff.= -------------------------------------------------
C     **             SQRT(ave(X**2)-ave(X)**2)*SQRT(ave(Y**2)-ave(Y)**2)
C     ******************************************************************
C     ******************************************************************
C     ******************************************************************
C  Most correlation between an "observed" map (isomorphous or
C  solventflattened or 2Fo -Fc   and the "Fcalc" maps generated from
C   known  coordinates.
C     ******************************************************************
C     ******************************************************************
C  January 1988  ejd
C  EXTRA MOD TO Average maps
C         or to exclude points from output map
C  when input fc map exists - ie to leave solvent density only
C         or to only include points in output map
C  when input fc map exists - ie to include the density fitted already.
C
C  15-Sep-93  IJT (at York)
C  Replace MCLOSE with MWCLOS and remove variables
C  Calc correlation coefficient using X' = X - Xo, Y' = Y - Yo where
C  Xo and Yo are working origins (uses mean values read from map files).
C  This is to avoid rounding errors when map mean is large compared
C  with map sigma.
C
C     ******************************************************************
C
      program overlapmap
C     ==================
C
C      implicit none
C
      integer ngrlim
      parameter(ngrlim=250000)
      real rin1(ngrlim)
      real rin2(ngrlim)
      real rin3(ngrlim)
      logical sigm
      real correl(10001),xave(10001),yave(10001),xxave(10001),
     1 yyave(10001),xyave(10001),rfac(10001)
      integer iave(10001)
C Chains
      integer maxchn
      parameter (maxchn=20)
      common /chain/ nchain, iresn(maxchn), iresc(maxchn), ires0(maxchn)
      integer nchain, iresn, iresc, ires0
      common /cchain/ chnlab(maxchn)
      character chnlab*1
C
C Things for PARSER ----
      integer maxtok
      parameter (maxtok=100)
      character key*4,cvalue(maxtok)*4,line*80
      integer ibeg(maxtok),iend(maxtok),ityp(maxtok),
     .     idec(maxtok),ntok
      real    fvalue(maxtok)
      logical lend
C
      character*4 word2
      character*80 tit1,tit2,tit3
      real cell3(6)
      integer nxyz3(3),iuvw3(3),nsec3
      integer jsec3,jz13,jz23,jx13,jx23,lspgr3
      real rhmen3, rhmin3, rhmax3
      common/inptt1/tit1
      common /input1/ nxyz1(3),ix11,ix21,iy11,iy21,iz11,iz21,
     &  jsec1,jx11,jx21,jz11,jz21,iuvw1(3),nsec1,kp1(3)
     & ,cell1(6),lspgr1,rhmin1,rhmax1,rhmen1,rhmnsq1
      integer nxyz1,ix11,ix21,iy11,iy21,iz11,iz21,
     &  jsec1,jx11,jx21,jz11,jz21,iuvw1,nsec1,kp1,lspgr1
      real cell1,rhmin1,rhmax1,rhmen1,rhmnsq1
      common/inptt2/tit2
      common /input2/ nxyz2(3),ix12,ix22,iy12,iy22,iz12,iz22,
     &  jsec2,jx12,jx22,jz12,jz22,iuvw2(3),nsec2,kp2(3)
     & ,cell2(6),lspgr2,rhmin2,rhmax2,rhmen2,scale1,scale2,
     &  rhmnsq2,sigm,sigmin,sigmax
      integer nxyz2,ix12,ix22,iy12,iy22,iz12,iz22,
     &  jsec2,jx12,jx22,jz12,jz22,iuvw2,nsec2,kp2,lspgr2
      real cell2,rhmin2,rhmax2,rhmen2,rhmnsq2,sigmin,sigmax
C
      integer nhpts,nvpts,jmax,imax,iun1,iun2,iun3,iun11,mode
C
C
C Locals
      integer lmode1, lmode2, lmode3, irssum, nchk
      real rhmnsq, scchk
C
C%
C
C READ HEADER TO GET INPUT MAP LIMITS
C
      call ccpfyp
      CALL CCPRCS(6,'OVERLAPMAP', '$Date: 2002/08/06 12:12:44 $')
      iun1=3
      iun2=4
      iun3=2
      lmode1=2
      lmode2=2
      sigm = .false.
C
      nchain = 0
      irssum = 0
C  Map 1 
      call mrdhdr(iun1,'MAPIN1',tit1,nsec1,iuvw1,nxyz1,jsec1,jz11,jz21,
     & jx11,jx21,cell1,lspgr1,lmode1,rhmin1,rhmax1,rhmen1,rhmnsq1)
C  Map 2 
      call mrdhdr(iun2,'MAPIN2',tit2,nsec2,iuvw2,nxyz2,jsec2,jz12,jz22,
     & jx12,jx22,cell2,lspgr2,lmode2,rhmin2,rhmax2,rhmen2,rhmnsq2)
C
C Check maps 1 and 2
      if(iuvw1(1) .ne. iuvw2(1) .or. iuvw1(2) .ne. iuvw2(2)
     & .or. iuvw1(3) .ne. iuvw2(3)) goto 400
      if(nxyz1(1) .ne. nxyz2(1) .or. nxyz1(2) .ne. nxyz2(2)
     & .or. nxyz1(3) .ne. nxyz2(3)) goto 500
      if(jsec1 .ne. jsec2) goto 600
      if(nsec1 .ne. nsec2) goto 700
      if(jz11  .ne.  jz12) goto 800
      if(jz21  .ne.  jz22) goto 900
      if(jx11  .ne.  jx12) goto 1000
      if(jx21  .ne.  jx22) goto 1100
C
      imax=1
      jmax=1
      nhpts=jz21-jz11+1
      nvpts=jx21-jx11+1
C  Check NHPTS*NVPTS LT NGRLIM
       nchk = nhpts*nvpts
       if(nchk.gt.ngrlim) call ccperr
     1  (1,'Reset NGRLIM dimension: currently NGRLIM=250000')
      write(6,'(//,a)') ' Input key words now'
      write(6,'(//,a,/,a)') '  EITHER',
     1' CORRELATE SECTION - to calculate correlation section by section'
      write(6,'(/,a,/,a)') '  OR',
     1' CORRELATE RESIDUE - to calculate correlation residue by residue'
      write(6,'(/,a,/,a)') '  OR',
     1 ' CORRELATE ATOM    - to calculate correlation atom by atom'
      write(6,'(/,a,/,a,/,a,/,a,/,a)') '  OR',
     1   ' REAL space R factor - calculated  residue by residue',
     1 ' This requires Map1 = FC  ',
     1 '               Map2 = Fobs ( scaled to FC)',
     1 '               Map3 = "Atom map" with residue flag'
      write(6,'(/,a,/,a)') '  OR',
     1 ' MAP AVERAGE - average two maps and output new map'
      write(6,'(//,a,/,a,a,/,a,a)') '  OR',
     1 ' MAP EXCLUDE - output map2 with points excluded if they ',
     1 ' exist in map1 ',
     1 ' eg output solvent density after masking with FC from ',
     1 ' protein coordinates.'
      write(6,'(//,a,/,a,a,/,a,a//)') '  OR',
     1 ' MAP INCLUDE - output map2 with points included only  ',
     1 'if they exist in map1 ',
     1 ' eg output protein density after masking with FC ',
     1 'from protein coordinates.'
C
      mode=-1000
C
C---- read command input
 10   line=' '
      ntok=maxtok
      call parser(
     .  key,line,ibeg,iend,ityp,fvalue,cvalue,idec,ntok,lend,.true.)
      if(lend) go to 200
      if(ntok.eq.0) go to 10
      call ccpupc(key)
C
      if(key .eq.'CORR' )then
         word2 = cvalue(2)
         call ccpupc(word2)
C
         if(word2.eq.'SECT')then
            write(6,*) ' Correlation section by section'
            mode=2
            imax=1
            jmax=500
C
         elseif(word2.eq.'RESI')then
            mode=3
            imax=2
            jmax=4999
            write(6,*) ' Correlation residue by residue'
C
         elseif(word2.eq.'ATOM')then
            mode=4
            imax=1
            jmax=5000
            write(6,*) ' Correlation atom by atom '
         endif
C
C  correlation on sigma level only compare patterson points when
C  either of patterson pixels are in range sigmin < RHO/sig < sigmax
C
      elseif (key .eq. 'PATT')then
         sigmin = 3.0
         sigmax = 12.0
         CALL GTPREA(2,sigmin,NTOK,ITYP,FVALUE)
         CALL GTPREA(3,sigmax,NTOK,ITYP,FVALUE)
         temp = sigmin
         if(sigmax.lt.sigmin) then
          sigmin = sigmax
          sigmax = temp
         end if
         write(6,*) ' Correlation section by section'
         mode=2
         imax=1
         jmax=500
         Sigm = .true.
      elseif (key .eq. 'REAL')then
         mode = 5
         imax=2
         jmax=4999
         write(6,*) ' Real space R factor  residue by residue'
C     Check "Fobs" and "Fc" map are on approximately the same scale.
         scchk = abs(rhmen2/rhmen1 -1.0)
         if( scchk.ge.0.3) write(6,'(a,/,a,/,2(a,F14.5))')
     1    ' Check Fobs and Fc map are on approximately the same scale.',
     1       ' Are you sure this is sensible? Mean densities different',
     1        ' RHO_mean map1 ',rhmen1,' RHO_mean map2',rhmen2
C
      elseif (key .eq. 'MAP')then
         word2 = cvalue(2)
         call ccpupc(word2)
C
         if(word2.eq.'AVER')then
            mode=0
            write(6,*) ' average two maps and output new map'
C
         elseif(word2.eq.'ADD ')then
            mode=-2
            scale1 = 1.0
            scale2 = 1.0
            if(ntok.gt.2) call gtprea(3,SCALE1,NTOK,ITYP,fvalue)
            if(ntok.gt.3) call gtprea(4,SCALE2,NTOK,ITYP,fvalue)
            write(6,*) ' add two maps and output new map'
            write(6,'(a,f5.1,a,f5.1)') ' Scale for mapin1',scale1,
     +                                 ' scale for mapin2',scale2
C
         elseif(word2.eq.'EXCL')then
            mode=1
            write(6,*)
     +           'exclude points from output map if they exist in map1'
C
C
         elseif(word2.eq.'INCL')then
            mode=-1
            write(6,*)
     +           'exclude points from output map if they exist in map1'
         endif
C
      elseif (key .eq. 'CHAI') then
         nchain = nchain + 1
         if (nchain .gt. maxchn) then
            write (6,'(A,I5)')
     $           ' **** Too many chains, maximum = ',maxchn
            call ccperr(1,' ** Too many chains **')
         endif
         chnlab(nchain) = line(ibeg(2):ibeg(2))
         call gtpint(3,iresn(nchain),ntok,ityp,fvalue)
         call gtpint(4,iresc(nchain),ntok,ityp,fvalue)
         ires0(nchain) = irssum
         irssum = irssum + iresc(nchain) - iresn(nchain) + 1
      elseif (key .eq. 'END') then
         goto 200
      endif
C
      go to 10
C
C----------- All command input read
C
 200  if(mode.eq.-1000) then
         call ccperr(1,' *****  Wrong key words given ! **')
      endif
C
      write(6,'(/,a,i5,///)') ' Check print - Mode = ',mode
C
      if(mode.ge.3) then
C   We need to open Map 3
         call mrdhdr(
     $       iun3,'MAPIN3',tit3,nsec3,iuvw3,nxyz3,jsec3,jz13,jz23,
     &       jx13,jx23,cell3,lspgr3,lmode3,rhmin3,rhmax3,rhmen3,rhmnsq)
C
C Check maps 1 and 3
         if(iuvw1(1) .ne. iuvw3(1) .or. iuvw1(2) .ne. iuvw3(2)
     &        .or. iuvw1(3) .ne. iuvw3(3)) goto 400
         if(nxyz1(1) .ne. nxyz3(1) .or. nxyz1(2) .ne. nxyz3(2)
     &        .or. nxyz1(3) .ne. nxyz3(3)) goto 500
         if(jsec1 .ne. jsec3) goto 600
         if(nsec1 .ne. nsec3) goto 700
         if(jz11  .ne.  jz13) goto 800
         if(jz21  .ne.  jz23) goto 900
         if(jx11  .ne.  jx13) goto 1000
         if(jx21  .ne.  jx23) goto 1100
      end if
C
      if(mode.le.1)then
         iun11=11
C      title=tit2(1:50)//' atom density set to zero'
         call mwrhdr(iun11,tit2,nsec1,iuvw1,nxyz1,jsec1,jz11,jz21,jx11,
     .        jx21,cell1,lspgr1,2)
C
C
C Copy symmetry to output file
         call msycpy(iun2,iun11)
      endif
C
C
      call mapovlp(rin1,rin2,rin3,nhpts,nvpts,
     1 jmax,imax,iun1,iun2,iun3,iun11,mode,
     1 correl,xave,yave,xxave,yyave,xyave,iave,
     1 rfac)
C
      call ccperr(0,' Normal termination ')
C
C
C
400   write(6,410)
410   format(5x,'***Fast,medium,slow axis do not match.***')
      call ccperr(1,' Fatal Error')
500   write(6,510)
510   format(5x,'***Grid sampling does not match.***')
      call ccperr(1,' Fatal Error')
600   write(6,610)
610   format(5x,'***Start of section does not match.***')
      call ccperr(1,' Fatal Error')
700   write(6,710)
710   format(5x,'***End of section does not match.***')
      call ccperr(1,' Fatal Error')
800   write(6,810)
810   format(5x,'***First point of column does not match.***')
      call ccperr(1,' Fatal Error')
900   write(6,910)
910   format(5x,'***Last point of column does not match.***')
      call ccperr(1,' Fatal Error')
1000  write(6,1010)
1010  format(5x,'***First point of row does not match.***')
      call ccperr(1,' Fatal Error')
1100  write(6,1110)
1110  format(5x,'***Last point of row does not match.***')
      call ccperr(1,' Fatal Error')
      end
C
C
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      subroutine mapovlp(rin1,rin2,rin3,nhpts,nvpts,
     1     jmax,imax,iun1,iun2,iun3,iun11,mode,
     1     correl,xave,yave,xxave,yyave,xyave,iave,
     1     rfac)
C
C      implicit none
C
      integer imax, jmax, nhpts,nvpts
      real rin1(nhpts,nvpts)
      real rin2(nhpts,nvpts)
      real rin3(nhpts,nvpts)
      real correl(0:jmax,imax)
C
      real xave(0:jmax,imax),yave(0:jmax,imax),xxave(0:jmax,imax),
     1     yyave(0:jmax,imax),xyave(0:jmax,imax),rfac(0:jmax,imax),
     1     xavesj(2),yavesj(2),xxavesj(2),
     1     yyavesj(2),xyavesj(2)
      integer iave(0:jmax,imax),iavesj(2)
      real    dtop(500),dbot(500)
      integer ndd(500)
      real xavet,yavet,xxavet,yyavet,xyavet
      real aa,bb,cc,aat,bbt,cct
      real rtot1,rtot2,ctot1,ctot2
C
C     Chains
      integer maxchn
      parameter (maxchn=20)
      common /chain/ nchain, iresn(maxchn), iresc(maxchn), ires0(maxchn)
      integer nchain, iresn, iresc, ires0
      common /cchain/ chnlab(maxchn)
      character chnlab*1
C
      logical sigm
      character*80 tit1,tit2
      common/inptt1/tit1
      common /input1/ nxyz1(3),ix11,ix21,iy11,iy21,iz11,iz21,
     &  jsec1,jx11,jx21,jz11,jz21,iuvw1(3),nsec1,kp1(3)
     & ,cell1(6),lspgr1,rhmin1,rhmax1,rhmen1,rhmnsq1
      integer nxyz1,ix11,ix21,iy11,iy21,iz11,iz21,
     &  jsec1,jx11,jx21,jz11,jz21,iuvw1,nsec1,kp1,lspgr1
      real cell1,rhmin1,rhmax1,rhmen1,rhmnsq1
      common/inptt2/tit2
      common /input2/ nxyz2(3),ix12,ix22,iy12,iy22,iz12,iz22,
     &  jsec2,jx12,jx22,jz12,jz22,iuvw2(3),nsec2,kp2(3)
     & ,cell2(6),lspgr2,rhmin2,rhmax2,rhmen2,scale1,scale2,
     &  rhmnsq2,sigm,sigmin,sigmax
      integer nxyz2,ix12,ix22,iy12,iy22,iz12,iz22,
     &  jsec2,jx12,jx22,jz12,jz22,iuvw2,nsec2,kp2,lspgr2
      real cell2,rhmin2,rhmax2,rhmen2,rhmnsq2,sigmin,sigmax
C
C Locals
      integer mnsc, mxsc, mnrs, mxrs, mnat, mxat
      integer i,j,k,n, lsec, iun1, iun2, iun3, iun11, ier, mode,
     $      ii, jj, mchflg, ires, jatm, navetx, navety, iavet,
     $     navetnegx, ichain
      real r1, yy, yyavetpos,  yyavetposx,  yyavetnegx, ttcorrel,
     $     bbcc, xwork, ywork
C
C do you want to look at map points sigma only
      if(sigm) then
        sig1min = sigmin*rhmnsq1
        sig1max = sigmax*rhmnsq1
        sig2min = sigmin*rhmnsq2
        sig2max = sigmax*rhmnsq2
        write(6,*) ' sigma cutoffs map1  = ',sig1min,sig1max
        write(6,*) ' sigma cutoffs map2  = ',sig2min,sig2max
      endif
      mnsc=1000000
      mxsc=-1000000
      mnrs=1000000
      mxrs=-1000000
      mnat=1000000
      mxat=-1000000
      aat = 0.0
      bbt = 0.0
      cct = 0.0
      iavet = 0
      navetx = 0
      navety = 0
      yyavetpos = 0.0
      yyavetposx = 0.0
      xyavet=0.0
      xavet=0.0
      yavet=0.0
      xxavet=0.0
      yyavet=0.0
      dtopt =0.0
      dbott =0.0
C     SET ZERO
      do 1 i=1,imax
         iavesj(i)=0
         yavesj(i)=0.0
         yyavesj(i)=0.0
         xyavesj(i)=0.0
         xavesj(i)=0.0
         xxavesj(i)=0.0
         do 2 j=0,jmax
            iave(j,i)=0
            yave(j,i)=0.0
            yyave(j,i)=0.0
            xyave(j,i)=0.0
            xave(j,i)=0.0
            xxave(j,i)=0.0
 2       continue
 1    continue
      do 555 l=1,500
        ndd(l) = 0
        dtop(l)=0.0
        dbot(l)=0.0
555   continue
C
      do 210 lsec=jsec1,jsec1+nsec1-1
C     WRITE(6,*)LSEC
C......
C     POSITION TO SECTION
         call mposn(iun1,lsec)
C
C     CORRECT SECTION, GENERATE OUTPUT SECTION
         call mgulpr(iun1,rin1,ier)
         if(ier.ne.0) go to 1200
         call mposn(iun2,lsec)
         call mgulpr(iun2,rin2,ier)
         if(ier.ne.0) go to 1200
         if(mode.ge.3) then
            call mposn(iun3,lsec)
            call mgulpr(iun3,rin3,ier)
            if(ier.ne.0) go to 1200
         end if
C......
C
C
         do 102 j=1,nvpts
            do 100 k=1,nhpts
C
C     WRITE(6,*)RIN1(K,J),RIN2(K,J)
C
C...............addtwomaps
               if(mode.eq.-2)then
                  rin1(k,j)=(scale1*rin1(k,j) + scale2*rin2(k,j))
               endif
C
C...............averagetwomaps
               if(mode.eq.0)then
                  rin1(k,j)=0.5*(rin1(k,j) + rin2(k,j))
               endif
C...............excludepointsfrom output map if they
C.....existinmap1
               if(mode.eq.1)then
                  r1=rin1(k,j)
                  rin1(k,j)=0.
                  if(abs(r1) .lt. 1.0E-06)rin1(k,j)=rin2(k,j)
               endif
C...............includepointsin output map if they
C.....existinmap1
               if(mode.eq.-1)then
                  r1=rin1(k,j)
                  rin1(k,j)=0.
                  if(r1.gt.0.)rin1(k,j)=rin2(k,j)
               endif
C......
C     WRITE(6,*)RIN1(K,J),RIN2(K,J)
               if(mode.le.1) go to 100
C
C
C     IF MODE EQ 2 DO CORRELATION section by section
               if (mode.eq.2)then
                  ii=1
                  jj=lsec
                  if (jj.lt.mnsc) mnsc = jj
                  if (jj.gt.mxsc) mxsc = jj
               endif
C
C
C     IF MODE EQ 3 DO CORRELATION RESIDUE BY RESIDUE
C     IF MODE EQ 5 DO Real space Rfactor RESIDUE BY RESIDUE
C     In either case we must extract the residue flag from RIN1.
               if (mode.eq.3 .or. mode.eq.5)then
                  ii=0
                  jj=-999
                  if(rin3(k,j).ne.0.)then
                     mchflg=nint(rin3(k,j)/1000000.)
                     rin3(k,j)=rin3(k,j)-mchflg*1000000.
                     ires=nint(rin3(k,j)/100.)
                     rin3(k,j)=rin3(k,j)-ires*100.
C     OK - now we have used RIN1 to get flag - if MODE = 5 overwrite
C     RIN1 with RIN3 - save rewriting correlation stuff..
                     ii=mchflg
C     mchflg set to 0 or 2 for side chain
C     mchflg set to 1 for main chain
                     jj=ires
                     if (jj.lt.mnrs) mnrs = jj
                     if (jj.gt.mxrs) mxrs = jj
                     if(ii.eq.0.and.jj.ne.0)ii=2
C     write(6,1101)ii,jj,rin1(k,j),k,j,lsec,ires,ii
 1101                format('  funny residue number at grid point',
     1                    /,2i8,f15.6,3i4,5x,2i8)
                  endif
               endif
C
C
C     IF MODE EQ 4 DO CORRELATION Atom by Atom
               if (mode.eq.4)then
                  ii=0
                  jj=-999
                  if(rin3(k,j).ne.0.)then
                     jatm=nint(rin3(k,j)/100.)
                     rin3(k,j)=rin3(k,j)-jatm*100.
C     DO NOT USE POINTS WHERE THERE IS VERY LITTLE "FC" DENSITY FOR ATOM
C     IMPORTANT FOR NEUTRON DATA
                     if(abs(rin3(k,j)).gt.0.0000001)then
                        ii=1
                        jj=jatm
                        if (jj.lt.mnat) mnat = jj
                        if (jj.gt.mxat) mxat = jj
                     endif
C     write(6,*)k,j,ii,jj,rin1(k,j)
                  endif
C     no correlation
               endif
               if(sigm) then
                if((rin1(k,j).ge.sig1min.and.rin1(k,j).le.sig1max).or.
     *           (rin2(k,j).ge.sig2min.and.rin2(k,j).le.sig2max)) then
                 xwork=rin1(k,j)-rhmen1
                 ywork=rin2(k,j)-rhmen2
                 if(jj.lt.0) go to 101
                 dtop(jj+1) = rin1(k,j)*rin2(k,j) + dtop(jj+1)
                 dbot(jj+1) = abs(rin1(k,j)*rin2(k,j)) + dbot(jj+1)
                 ndd(jj+1) = ndd(jj+1) + 1
                 dtopt = rin1(k,j)*rin2(k,j) + dtopt
                 dbott = abs(rin1(k,j)*rin2(k,j)) + dbott
                 xyave(jj,ii)=xyave(jj,ii)+xwork*ywork
                 xave(jj,ii)=xave(jj,ii)+xwork
                 yave(jj,ii)=yave(jj,ii)+ywork
                 iave(jj,ii)=iave(jj,ii)+1
                 xxave(jj,ii)=xxave(jj,ii)+xwork**2
                 yyave(jj,ii)=yyave(jj,ii)+ywork**2
C   Get sums for ii= 1 and ii=2
                 xyavesj(ii)=xyavesj(ii)+xwork*ywork
                 xavesj(ii)=xavesj(ii)+xwork
                 yavesj(ii)=yavesj(ii)+ywork
                 iavesj(ii)=iavesj(ii)+1
                 xxavesj(ii)=xxavesj(ii)+xwork**2
                 yyavesj(ii)=yyavesj(ii)+ywork**2
 101             continue
C
C     point inside fc boundary
                 xyavet=xyavet+xwork*ywork
                 xavet=xavet+xwork
C     count positive points in map 2
                 if(rin2(k,j).gt.0.0)then
                  navety=navety+1
                  yyavetpos=yyavetpos+ywork**2
C     count positive points in map 2 AND map1
                  if(rin1(k,j).gt.0)then
                     navetx=navetx+1
                     yyavetposx=yyavetposx+ywork**2
                  endif
                 endif
                 yavet=yavet+ywork
                 iavet=iavet + 1
                 xxavet=xxavet+xwork**2
                 yyavet=yyavet+ywork**2
                end if
               else
                xwork=rin1(k,j)-rhmen1
                ywork=rin2(k,j)-rhmen2
                if(jj.lt.0) go to 103
                xyave(jj,ii)=xyave(jj,ii)+xwork*ywork
                xave(jj,ii)=xave(jj,ii)+xwork
                yave(jj,ii)=yave(jj,ii)+ywork
                iave(jj,ii)=iave(jj,ii)+1
                xxave(jj,ii)=xxave(jj,ii)+xwork**2
                yyave(jj,ii)=yyave(jj,ii)+ywork**2
C   Get sums for ii= 1 and ii=2
                xyavesj(ii)=xyavesj(ii)+xwork*ywork
                xavesj(ii)=xavesj(ii)+xwork
                yavesj(ii)=yavesj(ii)+ywork
                iavesj(ii)=iavesj(ii)+1
                xxavesj(ii)=xxavesj(ii)+xwork**2
                yyavesj(ii)=yyavesj(ii)+ywork**2
 103            continue
C
C     point inside fc boundary
               xyavet=xyavet+xwork*ywork
               xavet=xavet+xwork
C     count positive points in map 2
               if(rin2(k,j).gt.0.0)then
                  navety=navety+1
                  yyavetpos=yyavetpos+ywork**2
C     count positive points in map 2 AND map1
                  if(rin1(k,j).gt.0)then
                     navetx=navetx+1
                     yyavetposx=yyavetposx+ywork**2
                  endif
                endif
                yavet=yavet+ywork
                iavet=iavet + 1
                xxavet=xxavet+xwork**2
                yyavet=yyavet+ywork**2
               end if
C
C...............
 100        continue
 102     continue
C...............
         if(mode.le.1 )then
            call mspew(iun11,rin1)
         endif
C...............
C     go back and test next section
 210  continue
CWRITE(6,*)LSEC
C
C     Close map
      if(mode.le.1 )then
         call mwclose(11)
         call ccperr(0,' Map written ')
      endif
C
C.............
      xyavet=xyavet/iavet
      xavet=xavet/iavet
      yavet=yavet/iavet
      xxavet=xxavet/iavet
      yyavet=yyavet/iavet
C.............
C
C     rms value of  points in map 2 which match map1
C     and points in map 2 which do not match map1
      navetnegx=navety- navetx
      yyavetnegx=yyavetpos - yyavetposx
      yy=0.
      if(navety.gt.0)yy=sqrt(yyavetpos/navety)
      yyavetpos=yy
      yy=0.
      if(navetx.gt.0)yy=sqrt(yyavetposx/navetx)
      yyavetposx=yy
      yy=0.
      if(navetnegx.gt.0)yy=sqrt(yyavetnegx/navetnegx)
      yyavetnegx=yy
C..........
      aat=xyavet-xavet*yavet
      bbt=sqrt(xxavet-xavet**2)
      cct=sqrt(yyavet-yavet**2)
      ttcorrel=aat/(bbt*cct)
C..................
      write(6,220)ttcorrel,yyavetpos,yyavetposx,yyavetnegx,
     1     navety,navetx,navetnegx
 220  format(6x,'Total correlation coefficient is : ',F9.4,
     1 /,' RMS positive density in  1.map2 (all positive  pts)',F12.3,
     1 /,' RMS positive density in  1.map2 (matching map1 pts)',F12.3,
     1 /,' RMS positive density in  1.map2 (nonmatching map1 pts)',
     +     F12.3,
     1 /,' Numbers of positive pts  in  above ranges',3I10)
C......
C     ******************************************************************
C     **                        ave(X*Y) - ave(X)*ave(Y)
C     * corr. coeff.= -------------------------------------------------
C     **             SQRT(ave(X**2)-ave(X)**2)*SQRT(ave(Y**2)-ave(Y)**2)
C     ******************************************************************
C     ******************************************************************
C......
C
C
      if(mode.eq.2) write(6,'(a,/,a,/a,/,a,/,a, i3,a,i3,a/,a,/,a)')
     $ ' ********************************************************',
     $ ' ********************************************************',
     $ ' $TABLE:',
     $ '  Correlation section by section :',
     $ ' $GRAPHS: Section correlation :',mnsc,'|',mxsc,'x-1|1:1,2: $$',
     1 'Section_number  Correlation_coefficient $$',
     1 'Section_number  Correlation_coefficient $$'
C
      if(mode.eq.3 .and. nchain .eq. 0)
     $  write(6,6003) ' ',' ',mnrs, mxrs
 6003 format(
     $ ' ********************************************************'/
     $ ' ********************************************************'/
     $ ' $TABLE:',
     $ '  Correlation Residue by Residue ',a,' :'/
     $ ' $GRAPHS: Residue correlation ',A,':',
     $     i3,'|',i3,'x-1|1:1,2,3: $$'/
     1 'Residue_number  Correlation_coeff_Main_chain ',/,
     1 ' Correlation_coeff_Side_chain $$',/,
     1 'ResNo      Main-corr-coeff Side-Corr-coeff $$')
C
      if(mode.eq.4) write(6,'(a,/,a,/,a,/,a,/,a,i3,a,i3,a,/,a,/,a)')
     $ ' ********************************************************',
     $ ' ********************************************************',
     $ ' $TABLE:',
     1 '  Atom number   Correlation coefficient  No of grid pts :',
     $ ' $GRAPHS: Atom correlation :',mnat,'|',mxat,'x-1|1:1,2: $$',
     1 'Atom_number  Correlation_coefficient Number_of_grid_points $$',
     1 'Atom_number  Correlation_coefficient No_of_grid_points $$'
C
      if(mode.eq.5 .and. nchain .eq. 0)
     $  write(6,6005) ' ',mnrs, mxrs, mnrs, mxrs, mnrs, mxrs
 6005 format(
     $ ' ********************************************************'/
     $ ' ********************************************************'/
     $ ' $TABLE:',
     $ '  Correlation Residue by Residue  , Real space R factor ',
     $ A,' :',/
     $ ' $GRAPHS: Residue correlation :',i3,'|',i3,'x-1|1:1,2,3:',/,
     $ ' : Real space R factor -main chain:',i3,'|',i3,'x-1|1:1,4:',/,
     $ ' : Real space R factor -side chain:',i3,'|',i3,'x-1|1:1,5:',/,
     1 '$$ Residue_number  Correlation_coeff_Main_chain ',
     1 ' Correlation_coeff_Side_chain ',
     1 ' R_factor_main_chain  ',
     1 ' R_factor_side_chain  ',
     1 ' Fobs(mc) Fcalc(mc) Fobs(sc) Fcalc(sc) $$ ',/,
     1 '  ResNo    Main-Corr.coeff  Side-corr-coeff ',
     1 ' RealR_main_chain  RealR_side_chain ',
     1 ' Fobs(mc) Fcalc(mc) Fobs(sc) Fcalc(sc) $$')
C
      ichain = -1000000
C
C  If imax eq 2 get main chain andside chain correlation.
       if(imax.gt.1) then
         do 215 ii=1,imax
            correl(1,ii) = 0.0
            rfac(1,ii)=0.0
            if(iavesj(ii).ne.0) then
               xyavesj(ii)=xyavesj(ii)/iavesj(ii)
               xavesj(ii)=xavesj(ii)/iavesj(ii)
               yavesj(ii)=yavesj(ii)/iavesj(ii)
               xxavesj(ii)=xxavesj(ii)/iavesj(ii)
               yyavesj(ii)=yyavesj(ii)/iavesj(ii)
               aa=xyavesj(ii)-xavesj(ii)*yavesj(ii)
               bb=sqrt(max(xxavesj(ii)-xavesj(ii)**2,0.))
               cc=sqrt(max(yyavesj(ii)-yavesj(ii)**2,0.))
               bbcc=bb*cc
               if(bbcc.gt.0.0)then
                  correl(1,ii)=aa/bbcc
                  if(mode.eq.5)rfac(1,ii)=(xavesj(ii)- yavesj(ii))/
     +                 (xavesj(ii)+ yavesj(ii))
               endif
            end if
 215     continue
            ctot1=correl(1,1)
            ctot2=correl(1,2)
C            write(6,'(a,10x,2f9.4)')' Total: ',correl(1,1),correl(1,2)
            if(mode.eq.5) then
              rtot1=rfac(1,1)
              rtot2=rfac(1,2)
C     +      write(6,'(a,10x,2f9.4)')' Total: ',rfac(1,1),rfac(1,2)
            end if
       end if
c
      do 221 jj=0,jmax
         do 222 ii=1,imax
            correl(jj,ii) = 0.0
            rfac(jj,ii)=0.0
            if(iave(jj,ii).ne.0) then
               xyave(jj,ii)=xyave(jj,ii)/iave(jj,ii)
               xave(jj,ii)=xave(jj,ii)/iave(jj,ii)
               yave(jj,ii)=yave(jj,ii)/iave(jj,ii)
               xxave(jj,ii)=xxave(jj,ii)/iave(jj,ii)
               yyave(jj,ii)=yyave(jj,ii)/iave(jj,ii)
               aa=xyave(jj,ii)-xave(jj,ii)*yave(jj,ii)
               bb=sqrt(max(xxave(jj,ii)-xave(jj,ii)**2,0.))
               cc=sqrt(max(yyave(jj,ii)-yave(jj,ii)**2,0.))
               bbcc=bb*cc
               if(bbcc.gt.0.0)then
                  correl(jj,ii)=aa/bbcc
                  if(mode.eq.5)rfac(jj,ii)=(xave(jj,ii)- yave(jj,ii))/
     +                 (xave(jj,ii)+ yave(jj,ii))
               endif
            end if
 222     continue
         if(mode.eq.2 .and. jj.lt.lsec)then
c         if(mode.eq.2 .and. correl(jj,1).ne.0.)then
            if(ndd(jj+1).eq.0.and.sigm) then
              cor=-999.
            elseif (sigm) then
              cor = dtop(jj+1)/dbot(jj+1)
            end if
C
            if(sigm) then
              write(6,253)jj,correl(jj,1),cor,ndd(jj+1)
            else
              write(6,251)jj,correl(jj,1)
            end if
C
         endif
 251     format(8x,i5,13x,2f9.4,5x,2f9.3)
 253     format(8x,'section cor1 cor2 num ',i5,13x,2f9.4,5x,i8)
 305     format(3x,i5,10x,2f9.4,5x,2f9.3,3x,4f12.3)
C
         if((mode.eq.3 ).and.
     1        (iave(jj,1) .ne. 0 .or. iave(jj,2) .ne. 0)) then
            j = jj
            if (nchain .gt. 0) then
               do 230, n = nchain, 1, -1
                  j = jj - ires0(n) + iresn(n) - 1
                  if (j .ge. iresn(n)) go to 235
 230           continue
               call ccperr(1,'** Impossible residue number **')
C
 235           if (n .ne. ichain) then
                  if (ichain .gt. -999999) write (6, '(1X,A)') '$$'
C  Write header for new chain
                  write (6, 6003)
     $                 '(chain '//chnlab(n)//')',
     $                 '(chain '//chnlab(n)//')',
     $                 iresn(n), iresc(n)
                  ichain = n
               endif
            endif
C
            write(6,305)j,correl(jj,1),correl(jj,2)
         endif
         if(( mode.eq.5).and.
     1        (iave(jj,1) .ne. 0 .or. iave(jj,2) .ne. 0)) then
            j = jj
            if (nchain .gt. 0) then
               do 240, n = nchain, 1, -1
                  j = jj - ires0(n) + iresn(n) - 1
                  if (j .ge. iresn(n)) go to 245
 240           continue
               call ccperr(1,'** Impossible residue number **')
C
 245           if (n .ne. ichain) then
                  if (ichain .gt. -999999) write (6, '(1X,A)') '$$'
C  Write header for new chain
                  write (6, 6005) '(chain '//chnlab(n)//')',
     $                 iresn(n), iresc(n), iresn(n), iresc(n),
     $                 iresn(n), iresc(n)
                  ichain = n
               endif
            endif
C
            write(6,305)
     $           j,correl(jj,1),correl(jj,2),rfac(jj,1),rfac(jj,2),
     1           xave(jj,1),yave(jj,1),xave(jj,2),yave(jj,2)
         endif
C
         if(mode.eq.4 .and. correl(jj,1).ne.0.)then
C     write(22,252)jj,iave(jj,1),XYAVE(JJ,1),XAVE(JJ,1),
C     1   YAVE(JJ,1),XXAVE(JJ,1),YYAVE(JJ,1),CORREL(JJ,1)
            write(6,252)jj,correl(jj,1),iave(jj,1)
         endif
 252     format(9x,i5,15x, f9.4,8x,i5)
C
C
 221  continue
      write(6,*)' $$'
      cort=-999.
      if(dbott.gt.0.0) then
         cort=dtopt/dbott
         write(6,254) cort
      end if
254   format(/,6x,'  cor2=Sum(map1*,map2)/Sum(abs(map1*map2)) ',f10.4)
C
            write(6,'(a,10x,2f9.4)')' Total: ',ctot1,ctot2
            if(mode.eq.5)
     +      write(6,'(a,10x,2f9.4)')' Total: ',rtot1,rtot2
      return
 1200 write(6,1210)
 1210 format(5x,'***Error in map input***')
      call ccperr(1,' Fatal error')
      end