subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max,
     &   ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2,
     &   nc_frac,nc_req_set,istr,llocal,ib,iprot,lprn)
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont),
     &   ishift,ishif2,nc_match
      integer i,is,js,ishif1,jfrag,nc_match1_max,n_shif1,n_shif2,
     &  nc_req_set,istr,nc_match_max,nc_req,nc_match1,ib,iprot
      double precision nc_frac
      logical llocal,lprn
      nc_match_max=0
      do i=1,ncont_ref
        nc_match_max=nc_match_max+
     &   min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
      enddo
      if (istr.eq.3) then
        nc_req=0
      else if (nc_req_set.eq.0) then
        nc_req=nc_match_max*nc_frac
      else
        nc_req = dmin1(nc_match_max*nc_frac+0.5d0,
     &    dfloat(nc_req_set)+1.0d-7)
      endif
c      write (iout,*) "match_contact: nc_req:",nc_req
c      write (iout,*) "nc_match_max",nc_match_max
c      write (iout,*) "jfrag",jfrag," n_shif1",n_shif1,
c     &   " n_shif2",n_shif2
C Match current contact map against reference contact map; exit, if at least
C half of the contacts match
      call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,
     &    ncont,icont,jfrag,llocal,ib,iprot,lprn)
      nc_match1_max=nc_match1
      if (lprn .and. nc_match.gt.0) write (iout,*) 
     &  "Shift:",0,0," nc_match1",nc_match1,
     &  " nc_match=",nc_match," req'd",nc_req
      if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
     &    nc_req.eq.0 .and. nc_match.eq.1) then
         ishif1=0
         ishif2=0
         return
      endif
C If sufficient matches are not found, try to shift contact maps up to three
C positions.
      if (n_shif1.gt.0) then
      do is=1,n_shif1
C The following four tries help to find shifted beta-sheet patterns
C Shift "left" strand backward
        call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,
     &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
        if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
        if (lprn .and. nc_match.gt.0) write (iout,*) 
     &    "Shift:",-is,0," nc_match1",nc_match1,
     &    " nc_match=",nc_match," req'd",nc_req
        if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
     &     nc_req.eq.0 .and. nc_match.eq.1) then
          ishif1=-is
          ishif2=0
          return
        endif
C Shift "left" strand forward
        call ncont_match(nc_match,nc_match1,is,0,ncont_ref,
     &      icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
        if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
        if (lprn .and. nc_match.gt.0) write (iout,*) 
     &   "Shift:",is,0," nc_match1",nc_match1,
     &   " nc_match=",nc_match," req'd",nc_req
        if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
     &     nc_req.eq.0 .and. nc_match.eq.1) then
          ishif1=is
          ishif2=0
          return
        endif
      enddo
      if (nc_req.eq.0) return
C Shift "right" strand backward
      do is=1,n_shif1
        call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,
     &     icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
        if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
        if (lprn .and. nc_match.gt.0) write (iout,*) 
     &    "Shift:",0,-is," nc_match1",nc_match1,
     &    " nc_match=",nc_match," req'd",nc_req
        if (nc_match.ge.nc_req) then
          ishif1=0
          ishif2=-is
          return
        endif
C Shift "right" strand upward
        call ncont_match(nc_match,nc_match1,0,is,ncont_ref,
     &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
        if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
        if (lprn .and. nc_match.gt.0) write (iout,*) 
     &    "Shift:",0,is," nc_match1",nc_match1,
     &    " nc_match=",nc_match," req'd",nc_req
        if (nc_match.ge.nc_req) then
          ishif1=0
          ishif2=is
          return
        endif
      enddo ! is
C Now try to shift both residues in contacts.
      do is=1,n_shif1
        do js=1,is
          if (js.ne.is) then
            call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,
     &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
            if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
            if (lprn .and. nc_match.gt.0) write (iout,*) 
     &         "Shift:",-is,-js," nc_match1",nc_match1,
     &         " nc_match=",nc_match," req'd",nc_req
            if (nc_match.ge.nc_req) then
              ishif1=-is
              ishif2=-js
              return
            endif
            call ncont_match(nc_match,nc_match1,is,js,ncont_ref,
     &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
            if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
            if (lprn .and. nc_match.gt.0) write (iout,*) 
     &        "Shift:",is,js," nc_match1",nc_match1,
     &        " nc_match=",nc_match," req'd",nc_req
            if (nc_match.ge.nc_req) then
              ishif1=is
              ishif2=js
              return
            endif
c
            call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,
     &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
            if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
            if (lprn .and. nc_match.gt.0) write (iout,*) 
     &        "Shift:",-js,-is," nc_match1",nc_match1,
     &        " nc_match=",nc_match," req'd",nc_req
            if (nc_match.ge.nc_req) then
              ishif1=-js
              ishif2=-is
              return
            endif
c
            call ncont_match(nc_match,nc_match1,js,is,ncont_ref,
     &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
            if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
            if (lprn .and. nc_match.gt.0) write (iout,*) 
     &        "Shift:",js,is," nc_match1",nc_match1,
     &        " nc_match=",nc_match," req'd",nc_req
            if (nc_match.ge.nc_req) then
              ishif1=js
              ishif2=is
              return
            endif
          endif
c
          if (is+js.le.n_shif1) then
          call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,
     &      icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
          if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
          if (lprn .and. nc_match.gt.0) write (iout,*) 
     &     "Shift:",-is,js," nc_match1",nc_match1,
     &     " nc_match=",nc_match," req'd",nc_req
          if (nc_match.ge.nc_req) then
            ishif1=-is
            ishif2=js
            return
          endif
c
          call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,
     &      icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
          if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
          if (lprn .and. nc_match.gt.0) write (iout,*) 
     &     "Shift:",js,-is," nc_match1",nc_match1,
     &     " nc_match=",nc_match," req'd",nc_req
          if (nc_match.ge.nc_req) then
            ishif1=js
            ishif2=-is
            return
          endif
          endif
c
        enddo !js
      enddo !is
      endif

      if (n_shif2.gt.0) then
      do is=1,n_shif2
        call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,
     &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
        if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
        if (lprn .and. nc_match.gt.0) write (iout,*) 
     &     "Shift:",-is,-is," nc_match1",nc_match1,
     &     " nc_match=",nc_match," req'd",nc_req
        if (nc_match.ge.nc_req) then
          ishif1=-is
          ishif2=-is
          return
        endif
        call ncont_match(nc_match,nc_match1,is,is,ncont_ref,
     &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
        if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
        if (lprn .and. nc_match.gt.0) write (iout,*) 
     &    "Shift:",is,is," nc_match1",nc_match1,
     &    " nc_match=",nc_match," req'd",nc_req
        if (nc_match.ge.nc_req) then
          ishif1=is
          ishif2=is
          return
        endif
      enddo
      endif
C If this point is reached, the contact maps are different. 
      nc_match=0
      ishif1=0
      ishif2=0
      return
      end
c-------------------------------------------------------------------------
      subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,
     &   ncont_ref,icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.INTERACT'
      include 'COMMON.GEO'
      include 'COMMON.COMPAR'
      integer nc_match,nc_match1,ishif1,jfrag,i,j,ic1,ic2,ib,iprot
      double precision diffang,fract
      logical llocal,lprn
      integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont),
     &   icont_match(2,maxcont),ishift,ishif2,nang_pair,
     &   iang_pair(2,maxres)
C Compare the contact map against the reference contact map; they're stored
C in ICONT and ICONT_REF, respectively. The current contact map can be shifted.
      if (lprn) write (iout,'(80(1h*))')
      nc_match=0
      nc_match1=0
c Check the local structure by comparing dihedral angles.
      if (llocal .and. ncont_ref.eq.0) then
c If there are no contacts just compare the dihedral angles and exit.
        call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag,ib,iprot),
     &    diffang,fract,ib,iprot,lprn)
        if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
     &   " ang_cut:",ang_cut(jfrag,ib,iprot)*rad2deg," fract",fract
        if (diffang.le.ang_cut(jfrag,ib,iprot) .and. 
     &   fract.ge.frac_min(jfrag,ib,iprot)) 
     &  then
          nc_match=1
        else
          nc_match=0
        endif
        return
      endif
      nang_pair=0
      do i=1,ncont
        ic1=icont(1,i)+ishif1
        ic2=icont(2,i)+ishif2
c        write (iout,*) "i",i," ic1",ic1," ic2",ic2
        if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
        do j=1,ncont_ref
          if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
            nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
            nc_match1=nc_match1+1
            icont_match(1,nc_match1)=ic1
            icont_match(2,nc_match1)=ic2
c            call add_angpair(icont(1,i),icont_ref(1,j),
c     &         nang_pair,iang_pair)
c            call add_angpair(icont(2,i),icont_ref(2,j),
c     &         nang_pair,iang_pair) 
            if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),
     &       " match",icont_ref(1,j),icont_ref(2,j),
     &       " shifts",ishif1,ishif2
            goto 10
          endif
        enddo 
   10   continue
      enddo
      if (lprn) then
        write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
        write (iout,*) "icont_match"
        do i=1,nc_match1
          write (iout,*) icont_match(1,i),icont_match(2,i)
        enddo
      endif
      if (llocal .and. nc_match.gt.0) then
        call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,
     &    ang_cut1(jfrag,ib,iprot),diffang,fract,ib,iprot)
        if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
     &   " ang_cut:",ang_cut(jfrag,ib,iprot)*rad2deg,
     &   " ang_cut1",ang_cut1(jfrag,ib,iprot)
        if (diffang.gt.ang_cut(jfrag,ib,iprot) 
     &    .or. fract.lt.frac_min(jfrag,ib,iprot)) nc_match=0
      endif
c      if (nc_match.gt.0) then
c        diffang = angnorm1(nang_pair,iang_pair,lprn)
c        if (diffang.gt.ang_cut(jfrag,ib,iprot)) nc_match=0
c      endif
      if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,
     &   " diffang",rad2deg*diffang," nc_match",nc_match
      return
      end
c------------------------------------------------------------------------------
      subroutine match_secondary(jfrag,isecstr,nsec_match,nsec_nomatch,
     &  avesig,ib,iprot,lprn)
c This subroutine compares the secondary structure (isecstr) of fragment jfrag 
c conformation considered to that of the reference conformation.
c Returns the number of equivalent residues (nsec_match).
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.PEPTCONT'
      include 'COMMON.COMPAR'
      logical lprn
      integer isecstr(maxres)
      integer i,j,j1,k,jc,jc1,jfrag,nsec_match,nsec_nomatch,ib,iprot,
     &  npart,nsig
      double precision sumsig,sumsigs1,sumsigs2,avesig,signature,sig,
     &  strand_signature
      external signature,strand_signature
      npart = npiece(jfrag,1,ib,iprot)
      nsec_match=0
      nsec_nomatch=0
      if (lprn) then
        write (iout,*) "match_secondary jfrag",jfrag," ifrag",
     &        (ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot),
     &        i=1,npart)
        write (iout,'(80i1)') (isec_ref(j,iprot),j=1,nres)
        write (iout,'(80i1)') (isecstr(j),j=1,nres)
      endif
      do i=1,npart
        do j=ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot)
c The residue has equivalent conformational state to that of the reference
c structure, if:
c  a) the conformational states are equal or
c  b) the reference state is a coil and that of the conformation considered 
c     is a strand or
c  c) the conformational state of the conformation considered is a strand
c     and that of the reference conformation is a coil.
c 10/28/02 - case (b) deleted.
          if (isecstr(j).eq.isec_ref(j,iprot) .or. 
c     &        isecstr(j).eq.0 .and. isec_ref(j,iprot).eq.1 .or.
     &        isec_ref(j,iprot).eq.0 .and. isecstr(j).eq.1) 
     &      nsec_match=nsec_match+1 
c 8/12/05 Unmatched secondary structure
          if (isecstr(j).eq.1 .and. isec_ref(j,iprot).eq.2 .or.
     &        isecstr(j).eq.2 .and. 
     &        (isec_ref(j,iprot).eq.0 .or. isec_ref(j,iprot).eq.1)
     &        .or. isecstr(j).eq.3 .and. isec_ref(j,iprot).ne.3)
     &       nsec_nomatch = nsec_nomatch+1
        enddo
      enddo
      avesig=strand_signature(jfrag,ib,iprot,lprn)
      return
      end
c------------------------------------------------------------------------------
      double precision function strand_signature(jfrag,ib,iprot,lprn)
c This subroutine compares the secondary structure (isecstr) of fragment jfrag 
c conformation considered to that of the reference conformation.
c Returns the number of equivalent residues (nsec_match).
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.PEPTCONT'
      include 'COMMON.COMPAR'
      logical lprn
      integer isecstr(maxres)
      integer i,j,j1,k,jc,jc1,jfrag,nsec_match,nsec_nomatch,ib,iprot,
     &  npart,nsig
      double precision sumsig,sumsigs1,sumsigs2,avesig,signature,sig
      external signature
c 12/15/07 Compute signature for beta structure
      nsig=0
      sumsig=0.0d0
      sumsigs1=0.0d0
      sumsigs2=0.0d0
      npart = npiece(jfrag,1,ib,iprot)
      do i=1,npart
        do j=ifrag(1,i,jfrag,ib,iprot)+1,ifrag(2,i,jfrag,ib,iprot)-1
          if (isec_ref(j,iprot).eq.1) then
            do k=1,ncont_pept_near(j,ib,iprot)
              jc=icont_pept_near(k,j,ib,iprot)
              call find_dir(j,jc,j1,jc1,ib,iprot)
              sig=signature(j,jc,j1,jc1,.false.)
              sumsig=sumsig+sig*sig_ref(k,j,ib,iprot)
              sumsigs1=sumsigs1+sig_ref(k,j,ib,iprot)**2
              sumsigs2=sumsigs2+sig*sig
              nsig=nsig+1
              if (lprn) then
                write (iout,*) "j",j," jc",jc," j1",j1," jc1",jc1,
     &           " isec_ref",isec_ref(j,iprot)," sig",sig,
     &           " sig_ref",sig_ref(k,j,ib,iprot)," sumsig",sumsig,
     &           " nsig",nsig
              endif
            enddo
          endif
        enddo
      enddo
      if (nsig.gt.0) then
        avesig=sumsig/dsqrt(sumsigs1*sumsigs2)
      else
        avesig=1.0d0
      endif
      if (lprn) write (iout,*) "avesig",avesig
      strand_signature=avesig
      return
      end
c------------------------------------------------------------------------------
      double precision function signature(j,jc,j1,jc1,lprn)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.PEPTCONT'
      include 'COMMON.COMPAR'
      integer j,j1,jc,jc1,k
      double precision u(3),v(3),w(3),z(3),unorm,vnorm,znorm
      double precision scalar
      external scalar
      logical lprn
c
c This function computes the signature of a residue involved in a beta-sheet;
c the signature is defined as the scalar product of a vector normal to the
c beta-sheet surface at contact beteen residues j and jc and the unit vector
c of the bisector of the CA(j-1)-CA(j)-CA(j+1) angle. 
c
c A.Liwo 12/15/07
c
c Compute the unit vectors CA(j)-CA(jc) (u) and CA(j-1)-CA(j+1) (v)
c
      do k=1,3
        u(k)=c(k,j1)+c(k,j+1)-c(k,jc1)-c(k,jc+1)
c        u(k)=c(k,j)+c(k,j+1)-c(k,jc)-c(k,jc+1)
        v(k)=c(k,j+1)-c(k,j-1)
      enddo
      unorm=dsqrt(scalar(u,u))
      vnorm=dsqrt(scalar(v,v))
      do k=1,3
        u(k)=u(k)/unorm
        v(k)=v(k)/vnorm
      enddo
      if (lprn) then
        write (iout,*) "vector u",(u(k),k=1,3)
        write (iout,*) "vector v",(v(k),k=1,3)
      endif
c
c Compute the normal vector to the beta-sheet at residue j (defined as the
c vector product of u and v).
c
      call vecpr(u,v,w)
      if (lprn) write (iout,*) "vector w",(w(k),k=1,3)
c
c Compute the unit vector of the CA(j-1)-CA(j)-CA(j+1) bisector (z).
c
      do k=1,3
        u(k)=c(k,j-1)-c(k,j)
        v(k)=c(k,j+1)-c(k,j)
      enddo
      unorm=dsqrt(scalar(u,u))
      vnorm=scalar(v,v)
      do k=1,3
        u(k)=u(k)/unorm
        v(k)=v(k)/vnorm
        z(k)=u(k)+v(k)
      enddo
      znorm=dsqrt(scalar(z,z))
      do k=1,3
        z(k)=z(k)/znorm
      enddo
      if (lprn) write (iout,*) "vector z",(z(k),k=1,3)
c
c Compute the signature
c
      signature=scalar(w,z)
      return
      end
c-------------------------------------------------------------------------------
      subroutine find_dir(j,jc,j1,jc1,ib,iprot)
      implicit none
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN'
      include 'COMMON.PEPTCONT'
      include 'COMMON.COMPAR'
      integer j,j1,jc,jc1,l,ib,iprot
      double precision u(3),v(3),w(3),z(3),unorm,vnorm,znorm
      double precision scalar
      external scalar
      logical lprn
      j1=j
      jc1=jc 
      do l=1,ncont_pept_near(j-1,ib,iprot)
        if (icont_pept_near(l,j-1,ib,iprot).eq.jc-1) then
          j1=j-1
          jc1=jc-1
          return
        else if (icont_pept_near(l,j-1,ib,iprot).eq.jc+1) then
          j1=j-1
          jc1=jc+2
          return
        endif
      enddo
      return
      end