Copyright (c) 2002 Peter Guntert. All rights reserved.
c     ==================================================================
      program mapper
c
      parameter (maxaa=50,maxs=6,maxr=300,maxf=maxr,maxa=4,maxfr=50,
     *           undef=0.0,big=1.0E20)
      common /mapdat/csh(maxa,maxfr,maxf),cs(maxs,maxa,maxaa),
     *               offset(maxa,maxr),devmax,offpro(maxa),
     *               mode(maxa),iraa(maxr),nfr(maxf)
c
      parameter (maxfrn=2*maxr,maxbst=50,maxar=maxa*maxr)
      dimension ir(maxr),indx(maxr,maxf),score(maxr,maxf),
     *          last(maxf),icount(maxf),p(maxr,maxf),f(maxf+1),
     *          fbest(maxbst),ibest(maxf,maxbst),icorr(maxf),
     *          ifrnum(maxfr,maxf),ifruse(maxr,maxf),tol(maxa),
     *          ifnum(maxf),akey(maxf),iord(maxf),prob(maxf),
     *          irank(maxf),jrank(maxf),scale(maxa),off(maxa),
     *          iovrlp(2,maxbst),ncsh(maxf),nused(maxf),snorm(maxf)
      logical pfr(maxr,maxf),selr(maxr),ovrlap(maxfr,maxf),found,
     *        compat(-1:maxfrn,maxfrn),used(maxr,maxf),best,incomp
      character aanam(maxaa)*5,rnam(maxr)*5,anam(maxa)*5,esc,
     *          seqfil*100,libfil*100,line*500,str*200,fmt*20,ch,rcode
      data nfr,icorr,ifnum,offset,scale
     *     /maxf*0,maxf*0,maxf*0,maxar*0.0,maxa*1.0/
c
      esc=char(int(92))
c     ----------------------------------------------- read fragment file
      read '(A)',libfil,seqfil
      i=index(libfil,'#')
      if (i.gt.0) libfil(i:)=' '
      i=index(seqfil,'#')
      if (i.gt.0) seqfil(i:)=' '
c     .................................................... read sequence
      open (unit=1,file=seqfil,status='OLD')
      call getseq (1,rnam,ir,nr)
      close (1)
c     ....................................................... atom types
      read *,maxmap,dmax,mfmax,cut,cutn,cutg,nbest
      devmax=0.5*dmax**2
   10 read '(A)',line
      i=index(line,'#')
      if (i.gt.0) line(i:)=' '
      if (line.eq.' ') go to 10
      read (line,'(10A8)') (anam(i),i=1,maxa)
      do i=1,maxa
        call lefstr (anam(i))
        if (anam(i).ne.' ') na=i
      end do
      read *,(tol(i),i=1,na)
c      read *,(scale(i),i=1,na)
      read *,(offpro(i),i=1,na)
      read '(A)',line
      do i=1,na
        mode(i)=1
        if (index(line(8*(i-1)+1:8*i),'rel').gt.0) mode(i)=2
      end do
c     .................................................. chemical shifts
      write (fmt,'(''('',I2,''F8.2,A)'')') na
      nf=1
   20 read (*,'(A)',end=30) line
        if (line.eq.' ') then
          nf=max(1,nf)
          if (nfr(nf).gt.0) then
            if (ifnum(nf).ge.0) then
              nf=nf+1
            else
              nfr(nf)=0
            end if
          end if
        end if
        i=index(line,'#')
        if (i.gt.0) line(i:)=' '
        if (line.eq.' ') go to 20
        if (line(1:7).eq.'correct') then
          read (line(8:),'(I10)') k
          do l=1,nr
            if (ir(l).eq.k) icorr(nf)=l
          end do
          go to 20
        else if (line(1:8).eq.'fragment') then
          read (line(9:),'(I10)') ifnum(nf)
          go to 20
        else if (line(1:7).eq.'offsets') then
          nf=0
          go to 20
        end if
c       ....................................................... fragment
        if (nf.ge.1) then
          nfr(nf)=nfr(nf)+1
          k=nfr(nf)
          read (line,fmt) (csh(i,k,nf),i=1,na),str
          i=index(str,'overlap')
          ovrlap(k,nf)=i.gt.0
          if (i.gt.0) str(i:i+6)=' '
          call select (str,rnam,ir,nr,selr)
c          print '(A)',str(1:lenstr(str))
c          do j=1,nr
c            if (selr(j)) print '(A,I5)',rnam(j),j
c          end do
          if (k.eq.1) then
            do j=1,nr
              pfr(j,nf)=selr(j)
            end do
          else
            do j=1,nr
              l=j+k-1
              pfr(j,nf)=pfr(j,nf) .and. selr(min(l,nr)) .and. l.le.nr
            end do
          end if
c       ......................................... chemical shift offsets
        else
          read (line,fmt) (off(i),i=1,na),str
          call select (str,rnam,ir,nr,selr)
          do i=1,na
            do j=1,nr
              if (selr(j)) offset(i,j)=off(i)
            end do
          end do
        end if
      go to 20
   30 if (nfr(nf).eq.0) nf=nf-1
      if (ifnum(nf).lt.0) nf=nf-1
      if (ifnum(1).eq.0) ifnum(1)=1
      do i=2,nf
        if (ifnum(i).eq.0) ifnum(i)=ifnum(i-1)+1
      end do
      do j=1,nr-1
        do i=1,na
          if (rnam(j+1).eq.'PRO') offset(i,j)=offset(i,j)+offpro(i)
        end do
      end do
c      print '(I5,2F10.2)',(j,(offset(i,j),i=1,na),j=1,nr)
c     ...................................... read chemical shift library
      if (libfil.eq.' ') libfil='shifts.lib'
      open (unit=1,file=libfil,status='OLD')
      call getlib (1,anam,na,aanam,naa,cs,maxs,maxa)
      close (1)
      do i=1,nr
        do j=1,naa
          if (rnam(i)(1:3).eq.aanam(j)) iraa(i)=j
        end do
      end do
      do j=1,naa
        do k=1,na
          cs(2,k,j)=cs(2,k,j)*scale(k)
        end do
      end do
c     ----------------------------------------------------- count strips
      n=0
      do i1=1,nf
        do j1=1,nfr(i1)
          n=n+1
          if (n.gt.maxfrn) then
            print '(A)','ERROR: Too many chemical shifts.'
            stop
          end if
          ifrnum(j1,i1)=n
          compat(0,n)=.true.
          compat(-1,n)=.false.
        end do
      end do
c     -------------------------------------------- set up overlap matrix
      do i1=1,nf
        do j1=1,nfr(i1)
          n1=ifrnum(j1,i1)
          do i2=i1+1,nf
            do j2=1,nfr(i2)
              n2=ifrnum(j2,i2)
              compat(n1,n2)=ovrlap(j1,i1) .or. ovrlap(j2,i2)
              if (compat(n1,n2)) then
                compat(n1,n2)=.false.
                do k=1,na
                  c1=csh(k,j1,i1)
                  c2=csh(k,j2,i2)
                  if (c1.ne.undef .and. c2.ne.undef .and.
     *                abs(c1-c2).gt.tol(k)) go to 40
                end do
                compat(n1,n2)=.true.
   40           continue
              end if
              compat(n2,n1)=compat(n1,n2)
            end do
          end do
        end do
      end do
c     ------------------------------------------------- calculate scores
      ntot=0
      do i=1,nf
        ncsh(i)=0
        snorm(i)=0.0
        do jr=1,nr
          if (pfr(jr,i)) then
            n=0
            t=0.0
            do j=1,nfr(i)
              ifr=ifrnum(j,i)
              k=iraa(jr+j-1)
              do l=1,na
                d=rscore(l,j,i,jr)
                if (d.eq.big) then
                  pfr(jr,i)=.false.
                  score(jr,i)=big
                  go to 50
                else if (d.ge.0.0) then
                  n=n+1
                  t=t+d
                end if
              end do
            end do
            if (ncsh(i).eq.0) ncsh(i)=n
c            if (ncsh(i).ne.n) print *,i,jr,ncsh(i),n
            if (ncsh(i).gt.0) then
              g=gammq(0.5*ncsh(i),t)
            else
              g=1.0
            end if
            if (g.lt.cut) then
              pfr(jr,i)=.false.
              score(jr,i)=big
            else
              score(jr,i)=t
              if (t.lt.50.0) snorm(i)=snorm(i)+exp(-t)
            end if
          else
            score(jr,i)=big
          end if
   50     continue
c          print *,jr,rnam(jr),score(jr,i)
        end do
        ntot=ntot+ncsh(i)
c       ........................... apply cutoff on relative probability
        if (snorm(i).gt.0.0) snorm(i)=-log(snorm(i))
        do jr=1,nr
          if (pfr(jr,i)) then
            if (exp(-score(jr,i)+snorm(i)).lt.cutn) then
              pfr(jr,i)=.false.
              score(jr,i)=big
            end if
          end if
        end do
      end do
c     --------------------------------------- map fragments individually
      print '(A)','___________________________________________________',
     *            ' ',
     *            'MAPPER 2.0',
     *            ' ',
     *            'Copyright (c) 2002 Peter Guntert',
     *            'All rights reserved.',
     *            '___________________________________________________'
      do i=1,nf
        m=0
        do j=1,nr
          if (pfr(j,i)) then
            m=m+1
            p(m,i)=score(j,i)
            indx(m,i)=j
            used(m,i)=.false.
          end if
        end do
        last(i)=min(m,maxmap)
        print '(//A,I4,A,I4,A)',
     *        'Fragment',ifnum(i),' (acceptable mappings:',last(i),'):'
        if (last(i).lt.1) then
          print '(A)','No acceptable mapping!'
          stop
        end if
        call rsort (p(1,i),indx(1,i),last(i),1)
        found=.false.
        do ii=1,last(i)
          l=indx(ii,i)
          found=found .or. l.eq.icorr(i)
          ch=' '
          if (l.eq.icorr(i)) ch='*'
          write (line,'(I2,F8.2,A,50I6)')
     *          ii,p(ii,i),ch,(ir(j),j=l,l+nfr(i)-1)
          write (line(lenstr(line)+1:),'(F8.2,A)')
     *          100.0*gammq(0.5*ncsh(i),p(ii,i)),'%'
          print '(/A)',line(1:lenstr(line))
          write (line,'(11X,50(3X,A3))') (rnam(j)(1:3),j=l,l+nfr(i)-1)
          write (line(lenstr(line)+1:),'(F8.2,A)')
     *          100.0*exp(-p(ii,i)+snorm(i)),'%'
          print '(A)',line(1:lenstr(line))
          do k=1,na
            sum=0.0
            do j=1,nfr(i)
              t=rscore(k,j,i,l)
              if (t.ge.0.0 .and. t.lt.big) then
                s=1.0
                if (mode(k).eq.1)
     *            s=csh(k,j,i)-(cs(1,k,iraa(l+j-1))-offset(k,l))
                write (line(6*(j-1)+1:),'(F6.2)') sign(t,s)
                sum=sum+t
              else
                line(6*(j-1)+1:)='     -'
              end if
            end do
            if (na.gt.1) write (line(lenstr(line)+1:),'(F8.2)') sum
            print '(A10,1X,A)',
     *            anam(k)(1:lenstr(anam(k))),line(1:lenstr(line))
          end do
        end do
        if (icorr(i).gt.0 .and. .not.found)
     *    print '(A)','Correct fragment not found!'
      end do
c     ------------------------------------------------ check consistency
      do mf=2,mfmax
        incomp=.false.
        do i=1,mf
          iord(i)=nf
        end do
        iord(mf)=0
  100   do imf=1,mf
          if (iord(imf).lt.nf) then
            do jmf=imf,1,-1
              if (jmf.eq.imf) then
                iord(jmf)=iord(imf)+1
              else
                iord(jmf)=iord(jmf+1)+1
                if (iord(jmf).gt.nf) go to 100
              end if
            end do
c            print '(10I4)',(iord(k),k=1,mf)
            do j=1,nr
              ifruse(j,mf+1)=0
            end do
            do i=1,mf-1
              icount(i)=last(iord(i))
            end do
            icount(mf)=0
c           ..................... loop over all possible global mappings
  120       do i=1,mf
              if (icount(i).lt.last(iord(i))) then
                do j=i,1,-1
                  jj=iord(j)
                  if (j.eq.i) then
                    k=icount(i)+1
                  else
                    k=1
                  end if
                  icount(j)=k
                  j1=j+1
c                 ........................................ overlap check
                  ibeg=indx(k,jj)
                  iend=ibeg+nfr(jj)-1
                  i1=ibeg-1
                  do l=ibeg,iend
                    ll=ifrnum(l-i1,jj)
                    lr=ifruse(l,j1)
                    if (.not.compat(lr,ll)) go to 120
                    if (lr.eq.0) then
                      ifruse(l,j)=ll
                    else
                      ifruse(l,j)=-1
                    end if
                  end do
                  do l=1,ibeg-1
                    ifruse(l,j)=ifruse(l,j1)
                  end do
                  do l=iend+1,nr
                    ifruse(l,j)=ifruse(l,j1)
                  end do
                end do
c               ................................... global mapping found
                go to 100
              end if
            end do
            if (.not.incomp) print '(/)'
            print '(A,10I4)','Incompatible fragments:',
     *            (ifnum(iord(k)),k=mf,1,-1)
            incomp=.true.
            go to 100
          end if
        end do
  150   if (incomp) stop
      end do
      if (mfmax.ge.2) print '(//A,I3,A)',
     *  'All subsets of up to',mfmax,' fragments are consistent.'
c     --------------------------------------------------- sort fragments
      do i=1,nf
        akey(i)=-last(i)+0.001*nfr(i)
      end do
      call rsort (akey,iord,nf,0)
c     ------------------------------------------- map fragments globally
      do i=1,nbest
        fbest(i)=big
      end do
c
      f(nf+1)=0.0
      do j=1,nr
        ifruse(j,nf+1)=0
      end do
      do i=1,nf-1
        icount(i)=last(iord(i))
      end do
      icount(nf)=0
      m=0
      n=0
c     --------------------------- loop over all possible global mappings
  200 do i=1,nf
        if (icount(i).lt.last(iord(i))) then
          do j=i,1,-1
            jj=iord(j)
            if (j.eq.i) then
              k=icount(i)+1
            else
              k=1
            end if
            icount(j)=k
            j1=j+1
            f(j)=f(j1)+p(k,jj)
c           .............................................. overlap check
            ibeg=indx(k,jj)
            iend=ibeg+nfr(jj)-1
            i1=ibeg-1
            do l=ibeg,iend
              ll=ifrnum(l-i1,jj)
              lr=ifruse(l,j1)
              if (.not.compat(lr,ll)) go to 200
              if (lr.eq.0) then
                ifruse(l,j)=ll
              else
                ifruse(l,j)=-1
              end if
            end do
            do l=1,ibeg-1
              ifruse(l,j)=ifruse(l,j1)
            end do
            do l=iend+1,nr
              ifruse(l,j)=ifruse(l,j1)
            end do
          end do
c         ......................................... global mapping found
          m=m+1
c          print '(I8,F8.2,50I4)',
c     *          m,f(1),(icount(j),j=1,nf),(indx(icount(j),j),j=1,nf)
          j=lookr(f(1),fbest,nbest)+1
          if (j.le.nbest) then
            do k=nbest,j+1,-1
              k1=k-1
              fbest(k)=fbest(k1)
              do l=1,nf
                ibest(l,k)=ibest(l,k1)
              end do
            end do
            fbest(j)=f(1)
            do l=1,nf
              ll=iord(l)
              ibest(ll,k)=indx(icount(l),ll)
            end do
          end if
          if (j.le.nbest .or. gammq(0.5*ntot,f(1)).ge.cutg) then
            n=n+1
            do l=1,nf
              ll=iord(l)
              used(icount(l),ll)=.true.
            end do
          end if
          go to 200
        end if
      end do
c     ---------------------------------------- report used possibilities
      print '(//A,I8,A)',
     *      'Possibilities used in the',n,' best global mappings:'
      if (n.eq.0) stop
      do l=1,nf
        nused(l)=0
        do i=1,last(l)
          if (used(i,l)) nused(l)=nused(l)+1
        end do
      end do
      print '(/A,100I4)','Fragment ',(ifnum(l),l=1,nf)
      print '(A,100I4)', 'Possib.  ',(nused(l),l=1,nf)
      print '(A)','Residue'
      iprev=0
      do j=1,nr
        line=' '
        do l=1,nf
          do k=1,last(l)
            if (used(k,l) .and.
     *          j.ge.indx(k,l) .and. j.lt.indx(k,l)+nfr(l)) then
              if (j.eq.indx(k,l)) then
                write (line(4*(l-1)+1:4*l),'(I4)') ifnum(l)
              else if (line(4*(l-1)+1:4*l).eq.' ') then
                if (indx(k,l).eq.icorr(l)) then
                  line(4*(l-1)+1:4*l)='   *'
                else
                  line(4*(l-1)+1:4*l)='   .'
                end if
              end if
            end if
          end do
        end do
        if (line.ne.' ') then
          if (iprev.eq.j-2) print '(I3,1X,A5,A)',ir(j-1),rnam(j-1)
          print '(I3,1X,A5,A)',ir(j),rnam(j),line(1:lenstr(line))
          iprev=j
        else if (iprev.eq.j-2) then
          print '(4X,A)','...'
        end if
      end do
c     -------------------------------------- report best global mappings
      nbest=min(m,nbest)
      print '(//A,I8,A)',
     *      'Best global mappings (possible global mappings:',m,'):'
      if (nbest.eq.0) stop
      print '(/A,50I4)','Residue ',(k,k=1,nbest)
      iprev=0
      do j=1,nr
        line=' '
        do k=1,nbest
          m=0
          do l=1,nf
            if (j.ge.ibest(l,k) .and. j.lt.ibest(l,k)+nfr(l)) then
              m=m+1
              iovrlp(m,k)=l
              if (line(4*(k-1)+1:4*k).eq.' ' .or.
     *            .not.ovrlap(j-ibest(l,k)+1,l)) then
                if (k.gt.1 .and. ibest(l,k).eq.ibest(l,1)) then
                  line(4*(k-1)+1:4*k)='  . '
                else
                  ch=' '
                  if (ibest(l,k).eq.icorr(l)) ch='*'
                  write (line(4*(k-1)+1:4*k),'(I3,A)') ifnum(l),ch
                end if
              end if
            end if
          end do
          if (m.gt.1 .and. line(4*k:4*k).eq.' ' .and.
     *        (k.eq.1 .or. iovrlp(1,k).ne.iovrlp(1,1) .or.
     *        iovrlp(2,k).ne.iovrlp(2,1))) line(4*k:4*k)='|'
        end do
        if (line.ne.' ') then
          if (iprev.eq.j-2) print '(I3,1X,A5,A)',ir(j-1),rnam(j-1)
          print '(I3,1X,A5,A)',ir(j),rnam(j),line(1:lenstr(line))
          iprev=j
        else if (iprev.eq.j-2) then
          print '(4X,A)','...'
        end if
      end do
c     --------------------------------- overview of best global mappings
      print '(//A4,A8,1X,100I4)','Rank','Score',(ifnum(l),l=1,nf)
      best=.false.
      do k=1,nbest+1
        ch='*'
        line=' '
        do l=1,nf
          if (k.le.nbest) then
            j=ibest(l,k)
          else
            j=indx(1,l)
          end if
          if (k.gt.1 .and. j.eq.ibest(l,1)) then
            line(4*(l-1)+1:4*l+1)='   .'
            if (j.eq.icorr(l)) line(4*(l-1)+1:4*l+1)='   *'
          else
            write (line(4*(l-1)+1:4*l+1),'(I4)') ir(j)
          end if
          if (j.ne.icorr(l)) ch=' '
        end do
        if (k.le.nbest) then
          best=best .or. ch.eq.'*'
          print '(I4,F8.2,2A)',k,fbest(k),ch,line(1:lenstr(line))
        else
          t=0.0
          do i=1,nf
            t=t+p(1,i)
          end do
          print '(/A4,F8.2,2A)','min',t,ch,line(1:lenstr(line))
        end if
      end do
      n=0
      do l=1,nf
        do i=1,last(l)
          do k=1,nbest
            if (ibest(l,k).eq.indx(i,l)) then
              if (k.eq.1) irank(l)=i
              jrank(l)=i
            end if
          end do
        end do
      end do
c
      do i=1,nf
        do j=1,nfr(i)
          do l=1,na
            t=rscore(l,j,i,ibest(i,1))
          end do
        end do
        prob(i)=gammq(0.5*ncsh(i),score(ibest(i,1),i))
      end do
      print '(/A,100I4)','Length       ',(nfr(l),l=1,nf)
      print '(A,100I4)', 'Possibilities',(last(l),l=1,nf)
      print '(A,100I4)', 'Used possib. ',(nused(l),l=1,nf)
      print '(A,100I4)', 'Highest rank ',(jrank(l),l=1,nf)
      print '(/A)','Best mapping:'
      print '(A,100I4)', 'Rank         ',(irank(l),l=1,nf)
      print '(A,100I4)', 'Score %      ',
     *      (nint(100.0*score(ibest(l,1),l)/fbest(1)),l=1,nf)
      print '(A,100I4)', 'Probability %',
     *      (nint(100.0*prob(l)),l=1,nf)
      print '(A,F7.2,A)','Total probability:',
     *      100.0*gammq(0.5*ntot,fbest(1)),' %'
c     -------------------------------------------------- correct mapping
c      if (.not.best) then
c        t=0.0
c        do l=1,nf
c          t=t+score(icorr(l),l)
c        end do
c        print '(A)','Correct global mapping not found:'
c        print '(F12.2,A,100I4)',t,'*',(ir(icorr(l)),l=1,nf)
c      end if
c     ------------------------------------------------------------- plot
      open (1,file='mapper.grf',status='unknown')
      write (1,'(A)') '# Output from mapper'
      write (1,'(A)') '#RECTANGLE y'
      write (1,'(A)') '#ROTATE'
      do k=1,na
        if (k.eq.1) then
          write (1,'(A,4F10.2)')
     *          '#VIEWPORT',-1.0,1.0,0.7-0.5*min(0.65,1.3/na),0.7
        else
          write (1,'(A)') '#NEXTPLOT'
          write (1,'(A)') '#VIEWPORT v'
        end if
c       ................. plot chemical shifts from sequence and library
        write (1,'(A)') '#RIGHTSTRING -0.05 0.5 1.0 '//esc//'Symbol w'
        write (1,'(2A)') '#CENTERSTRING 0.5 0.9 1.2 ',anam(k)
        tmin=1.0e20
        tmax=-1.0e20
        smax=0.0
        do j=1,nr
          t=cs(1,k,iraa(j))
          if (t.ne.0.0) then
            s=cs(2,k,iraa(j))
            tmin=min(tmin,t-3.0*s)
            tmax=max(tmax,t+3.0*s)
            smax=max(smax,s)
          end if
        end do
        smax=3.0*smax
        write (1,'(A,4F10.2)') '#MINMAX',0.5,real(nr)+0.5,tmin,tmax
        write (1,'(A)') '#XTICSONLY'
        write (1,'(A)') '#REPRES p'
        write (1,'(A)') '#CONNECT t'
        write (1,'(A)') '#LINES 0'
        do j=1,nr
          t=cs(1,k,iraa(j))
          if (t.ne.0.0) then
            t=t+sign(0.06*(tmax-tmin),t-0.5*(tmin+tmax))
            write (1,'(''#CENTERTEXT'',I4,2F10.2,2X,A)')
     *            j,t,0.4,rcode(rnam(j))
          end if
        end do
        do j=1,nr
          t=cs(1,k,iraa(j))
          if (t.eq.0.0) t=999999.0
          write (1,'(I4,F10.2)') j,t
        end do
c       .................. plot deviations for fragments in best mapping
        write (1,'(A)') '#NEXTPLOT'
        write (1,'(A)') '#VIEWPORT v'
        write (1,'(A)') '#RIGHTSTRING -0.05 0.5 1.0 '//esc//'Symbol Dw'
        if (k.lt.na) write (1,'(A)') '#XTICSONLY'
        write (1,'(A,4F10.2)')
     *        '#MINMAX',0.5,real(nr)+0.5,-smax,smax
        if (k.eq.na) write (1,'(A)') '#XTEXT Sequence'
        write (1,'(''#CENTERTEXT'',3F10.2,I5)')
     *        (ibest(l,1)+0.5*(nfr(l)-1),0.9*smax,0.5,ifnum(l),l=1,nf)
        do l=1,nf
          if (l.ne.1) write (1,'(A)') '#NEXTCURVE'
          write (1,'(A)') '#REPRES -'
          write (1,'(A)') '#CONNECT t'
          write (1,'(A)') '#LINES 2'
          write (1,'(2F10.2)') ibest(l,1)-0.5,0.0,ibest(l,1)-0.5,smax
          write (1,'(A)') '#NEXTCURVE'
          write (1,'(2F10.2)')
     *          ibest(l,1)+nfr(l)-0.5,0.0,ibest(l,1)+nfr(l)-0.5,smax
          write (1,'(A)') '#NEXTCURVE'
          write (1,'(A)') '#CONNECT -'
          do j=2,nbest
            if (ibest(l,j).ne.ibest(l,1)) go to 410
          end do
  410     if (j.le.nbest) then
            write (1,'(A)') '#REPRES h'
          else
            write (1,'(A)') '#REPRES b'
          end if
          write (1,'(I4,F10.2)') ibest(l,1)-1,0.0
          do j=1,nfr(l)
            jj=ibest(l,1)+j-1
            ss=cs(1,k,iraa(jj))-offset(k,jj)
            t=csh(k,j,l)-ss
            if (csh(k,j,l).eq.0.0) t=0.0
            t=min(0.9999*smax,max(-0.9999*smax,t))
            write (1,'(I4,F10.2)') jj,t
          end do
          write (1,'(I4,F10.2)') ibest(l,1)+nfr(l),0.0
        end do
      end do
      end
c     ==================================================================
c     l  = 1,...,maxa     atom
c     j  = 1,...,nfr(i)   residue position in fragment
c     i  = 1,...,nf       fragment
c     jr = 1,...,nr       position of first residue in sequence
c
      function rscore(l,j,i,jr)
c
      parameter (maxaa=50,maxs=6,maxr=300,maxf=maxr,maxa=4,maxfr=50,
     *           undef=0.0,big=1.0E20)
      common /mapdat/csh(maxa,maxfr,maxf),cs(maxs,maxa,maxaa),
     *               offset(maxa,maxr),devmax,offpro(maxa),
     *               mode(maxa),iraa(maxr),nfr(maxf)
c
      parameter (relfac=0.75)
c
      dev(c,ref,sd)=0.5*((c-ref)/sd)**2
c
      rscore=-1.0
      if (mode(l).eq.1) then
        if (csh(l,j,i).eq.undef) return
        rscore=big
        k=iraa(jr+j-1)
        if (cs(1,l,k).eq.undef) return
        ss=cs(1,l,k)-offset(l,jr+j-1)
        d=dev(csh(l,j,i),ss,cs(2,l,k))
        if (d.gt.devmax) return
        rscore=d
      else if (mode(l).eq.2) then
        if (j.eq.nfr(i)) return
        if (csh(l,j,i).eq.undef .or. csh(l,j+1,i).eq.undef) return
        rscore=big
        k=iraa(jr+j-1)
        kk=iraa(jr+j)
        if (cs(1,l,k).eq.undef .or. cs(1,l,kk).eq.undef) return
        s=cs(1,l,k)-offset(l,jr+j-1)
        ss=cs(1,l,kk)-offset(l,jr+j)
        sd=sqrt(cs(2,l,k)**2+cs(2,l,kk)**2)*relfac
        d=dev(csh(l,j+1,i)-csh(l,j,i),ss-s,sd)
        if (d.gt.devmax) return
        rscore=d
      end if
      end
c     ==================================================================
      function lookr (r,tab,n)
      dimension tab(n)
      lookr=0
      ju=n+1
   10 if (ju-lookr.gt.1) then
        m=(ju+lookr)/2
        if (r.ge.tab(m)) then
          lookr=m
        else
          ju=m
        end if
      go to 10
      end if
      end
c     ==================================================================
      subroutine select (str,rnam,ir,nr,selr)
c
      parameter (maxr=300)
      character*(*) str,rnam(nr),s*200,p*200,op*1
      dimension ir(nr)
      logical selr(nr),rflag(maxr),allr,start,match,illexp
c
      if (str.eq.' ') then
        do j=1,nr
          selr(j)=.true.
        end do
        return
      end if
      op='='
      allr=.true.
      do j=1,nr
        rflag(j)=.false.
        selr(j)=.false.
      end do
c     --------------------------------------------- loop over parameters
      s=str
      call lefstr (s)
      start=.true.
   10 continue
        l=lenstr(s)
        if (l.gt.0) then
          do i=1,l
            if (s(i:i).le.' ') go to 21
          end do
   21     p=s(1:i-1)
        else
          p=' '
        end if
c        print *,'s=',s(:50)
c        print *,'p=',p(:50)
c       .............................................. combine selection
        if (p.eq.'=' .or. p.eq.'+' .or. p.eq.'-' .or.
     *      p.eq.'/' .or. p.eq.'!' .or. p.eq.' ') then
          if (.not.start) then
            do j=1,nr
              rflag(j)=rflag(j) .or. allr
            end do
            if (op.eq.'=') then
              do j=1,nr
                selr(j)=rflag(j)
              end do
            else if (op.eq.'+') then
              do j=1,nr
                selr(j)=selr(j) .or. rflag(j)
              end do
            else if (op.eq.'-') then
              do j=1,nr
                selr(j)=selr(j) .and. .not.rflag(j)
              end do
            else if (op.eq.'/') then
              do j=1,nr
                selr(j)=selr(j) .and. rflag(j)
              end do
            else if (op.eq.'!') then
              do j=1,nr
                selr(j)=.not.rflag(j)
              end do
            end if
            if (p.eq.' ') return
            allr=.true.
            do j=1,nr
              rflag(j)=.false.
            end do
          end if
          op=p(1:1)
c       ................................................... residue name
        else if ((p(1:1).ge.'A' .and. p.le.'Z') .or.
     *           (p(1:1).ge.'a' .and. p.le.'z') .or.
     *           p(1:1).eq.'*' .or. p(1:1).eq.'?') then
          allr=.false.
          if (p.eq.'FIRST') then
            rflag(1)=.true.
          else if (p.eq.'first') then
            rflag(1)=.true.
            do j=2,nr
              rflag(j)=rflag(j) .or. ir(j).ne.ir(j-1)+1
            end do
          else if (p.eq.'LAST') then
            rflag(nr)=.true.
          else if (p.eq.'last') then
            rflag(nr)=.true.
            do j=1,nr-1
              rflag(j)=rflag(j) .or. ir(j).ne.ir(j+1)-1
            end do
          else
            do j=1,nr
              rflag(j)=rflag(j) .or. match(rnam(j),p)
            end do
          end if
c       ........................................... residue number range
        else
          allr=.false.
          ibeg=ir(1)
          iend=ir(nr)
          call irange (p,ibeg,iend,illexp)
c          print *,p(1:10),ibeg,iend
          if (illexp) then
            print '(3A)','Illegal residue range "',p(1:lenstr(p)),'".'
            stop
          end if
          do j=1,nr
            rflag(j)=rflag(j) .or. (ir(j).ge.ibeg .and. ir(j).le.iend)
          end do
        end if
        start=.false.
        p=s(i:)
        s=p
        call lefstr (s)
      go to 10
      end
c     ==================================================================
      subroutine getlib (iunit,anam,na,aanam,naa,cs,ns,maxa)
c
      parameter (undef=0.0)
      dimension cs(ns,maxa,*)
      character*(*) anam(na),aanam(*),rn*5,line*80
c
      rn=' '
      naa=0
   10 read (iunit,'(A)',end=50) line
        i=index(line,'#')
        if (i.gt.0) line(i:)=' '
        if (line.eq.' ') go to 10
        if (line(1:3).ne.rn) then
          naa=naa+1
          rn=line(1:3)
          aanam(naa)=rn
          do j=1,na
            cs(1,j,naa)=undef
          end do
        end if
        do j=1,na
          if (line(5:8).eq.anam(j))
     *      read (line,'(21X,6F8.2)') (cs(i,j,naa),i=1,ns)
        end do
      go to 10
   50 continue
      end
c     ==================================================================
      subroutine getseq (iunit,rnam,ir,nr)
c
      dimension ir(*)
      character*(*) rnam(*),line*81,el*5
c
      nr=0
   10 continue
        read (iunit,'(A)',end=50) line
        line(81:81)=' '
        do i=1,81
          if (line(i:i).eq.'#') go to 30
        end do
   30   k=0
        el=' '
        do j=1,i-1
          if (line(j:j).le.' ') then
            if (k.gt.0) then
              if ((el(1:1).ge.'0' .and. el(1:1).le.'9') .or.
     *            el(1:1).eq.'-' .or. el(1:1).eq.'+') then
                read (el,'(BN,I5)') ir(nr)
              else
                nr=nr+1
                rnam(nr)=el
                if (nr.eq.1) then
                  ir(nr)=1
                else
                  ir(nr)=ir(nr-1)+1
                end if
              end if
              k=0
              el=' '
            end if
          else
            k=k+1
            el(k:k)=line(j:j)
          end if
        end do
      go to 10
   50 return
      end
c     ==================================================================
      subroutine rsort (akey,indx,n,mode)
c
      parameter (nstack=30)
      dimension akey(n),indx(n),isl(nstack),isr(nstack)
c
      if (mode.eq.0) then
        do i=1,n
          indx(i)=i
        end do
      end if
      js=1
      isl(1)=1
      isr(1)=n
   50 continue
        jl=isl(js)
        jr=isr(js)
        js=js-1
   52   continue
          i=jl
          j=jr
          xx=akey((jl+jr)/2)
   53     continue
   54       if (akey(i).lt.xx) then
              i=i+1
            go to 54
            end if
   55       if (akey(j).gt.xx) then
              j=j-1
            go to 55
            end if
            if (i.gt.j) go to 58
            t=akey(i)
            akey(i)=akey(j)
            akey(j)=t
            k=indx(i)
            indx(i)=indx(j)
            indx(j)=k
            i=i+1
            j=j-1
          if (i.le.j) go to 53
   58     if (j-jl.lt.jr-i) then
            if (i.lt.jr) then
              js=js+1
              isl(js)=i
              isr(js)=jr
            end if
            jr=j
          else
            if (jl.lt.j) then
              js=js+1
              isl(js)=jl
              isr(js)=j
            end if
            jl=i
          end if
        if (jl.lt.jr) go to 52
      if (js.gt.0) go to 50
      end
c     ==================================================================
      function lenstr(s)
      character*(*) s
      do l=len(s),1,-1
        if (s(l:l).gt.' ') go to 20
      end do
   20 lenstr=l
      end
c     ==================================================================
      subroutine lefstr (s)
c
      character*(*) s
      logical flag
c
      if (s(1:1).gt.' ' .or. s.eq.' ') return
      i=0
      flag=.false.
      do l=1,len(s)
        flag=flag .or. s(l:l).gt.' '
        if (flag) then
          i=i+1
          s(i:i)=s(l:l)
        end if
      end do
      s(i+1:)=' '
      end
c     ==================================================================
c     MATCH:    Does the string S match the string O which may
c               contain wildcards?
c
c               The two possible wildcards are:
c                 "?" ... matches exactly one character
c                 "*" ... matches any number (including zero)
c                         of characters
c               The two strings S and O will not be changed.
c
c               Written in standard FORTRAN-77.
c               Peter Guentert, 28-11-1988
c     ------------------------------------------------------------------
      logical function match(s,o)
c
      character*(*) s,o
      logical fix,toend
c
      ls=lenstr(s)
      lo=lenstr(o)
      match=.false.
      fix=.true.
      jo=1
      js=1
   10 if (jo.le.lo) then
        if (o(jo:jo).eq.'*') then
          fix=.false.
          jo=jo+1
        else
          j=index(o(jo:lo),'*')-1
          toend=j.lt.0
          if (toend) j=lo-jo+1
          if (fix) then
            imax=js
          else
            imax=ls-j+1
            if (toend) js=imax
          end if
          do i=js,imax
            do l=jo,jo+j-1
              k=i+l-jo
              if (k.gt.ls) return
              if (o(l:l).ne.'?') then
                if (o(l:l).ne.s(k:k)) go to 30
              end if
            end do
            jo=jo+j
            js=i+j
            fix=.true.
            go to 10
   30       continue
          end do
          return
        end if
      go to 10
      end if
      match=js.gt.ls .or. .not. fix
      end
c     ==================================================================
      subroutine irange (str,ibeg,iend,illexp)
c
      logical illexp
      character*(*) str,s*20
c
      illexp=.false.
      if (str.eq.'*') return
      jbeg=ibeg
      jend=iend
      i=index(str,'..')
      if (i.eq.0) then
        s=str
        read (s,'(I20)') ibeg
        iend=ibeg
      else
        if (i.gt.1) then
          s=str(1:i-1)
          read (s,'(I20)') ibeg
        end if
        if (i.lt.lenstr(str)-1) then
          s=str(i+2:)
          read (s,'(I20)') iend
        end if
      end if
      illexp=illexp .or. ibeg.lt.jbeg .or. ibeg.gt.jend .or.
     *                   iend.lt.ibeg .or. iend.gt.jend
      end
c     ==================================================================
      character function rcode(rnam)
      character*(*) rnam
      if (rnam(1:3).eq.'ALA') then
        rcode='A'
      else if (rnam(1:3).eq.'ARG') then
        rcode='R'
      else if (rnam(1:3).eq.'ASN') then
        rcode='N'
      else if (rnam(1:3).eq.'ASP') then
        rcode='D'
      else if (rnam(1:3).eq.'CYS') then
        rcode='C'
      else if (rnam(1:3).eq.'GLN') then
        rcode='Q'
      else if (rnam(1:3).eq.'GLU') then
        rcode='E'
      else if (rnam(1:3).eq.'GLY') then
        rcode='G'
      else if (rnam(1:3).eq.'HIS') then
        rcode='H'
      else if (rnam(1:3).eq.'ILE') then
        rcode='I'
      else if (rnam(1:3).eq.'LEU') then
        rcode='L'
      else if (rnam(1:3).eq.'LYS') then
        rcode='K'
      else if (rnam(1:3).eq.'MET') then
        rcode='M'
      else if (rnam(1:3).eq.'PHE') then
        rcode='F'
      else if (rnam(1:3).eq.'PRO') then
        rcode='P'
      else if (rnam(1:3).eq.'SER') then
        rcode='S'
      else if (rnam(1:3).eq.'THR') then
        rcode='T'
      else if (rnam(1:3).eq.'TRP') then
        rcode='W'
      else if (rnam(1:3).eq.'TYR') then
        rcode='Y'
      else if (rnam(1:3).eq.'VAL') then
        rcode='V'
      else
        rcode='?'
      end if
      end
c     ==================================================================
      function gammq (a,x)
c
      if (x.lt.0.0 .or. a.le.0.0) then
c        print '(A)','GAMMQ: Illegal argument.'
c        stop
        gammq=1.0
        return
      end if
      if (x.lt.a+1.0) then
        call gser (gamser,a,x,gln)
        gammq=1.0-gamser
      else
        call gcf (gammcf,a,x,gln)
        gammq=gammcf
      end if
      end
c     ==================================================================
      subroutine gser (gamser,a,x,gln)
c
      parameter (itmax=100,eps=3.0e-7)
c
      gln=gammln(a)
      if (x.le.0.0) then
        if (x.lt.0.0) then
          print '(A)','GSER: Illegal argument.'
          stop
        end if
        gamser=0.0
        return
      end if
      ap=a
      sum=1.0/a
      del=sum
      do n=1,itmax
        ap=ap+1.0
        del=del*x/ap
        sum=sum+del
        if (abs(del).lt.abs(sum)*eps) go to 1
      end do
      print '(A)','GSER: A too large, ITMAX too small.'
      stop
    1 gamser=sum*exp(-x+a*log(x)-gln)
c      print *,'gamser: ',x,a,gln,-x+a*log(x)-gln,sum
      end
c     ==================================================================
      subroutine gcf (gammcf,a,x,gln)
c
      parameter (itmax=100,eps=3.0e-7)
c
      gln=gammln(a)
      gold=0.0
      a0=1.0
      a1=x
      b0=0.0
      b1=1.0
      fac=1.0
      do n=1,itmax
        an=real(n)
        ana=an-a
        a0=(a1+a0*ana)*fac
        b0=(b1+b0*ana)*fac
        anf=an*fac
        a1=x*a0+anf*a1
        b1=x*b0+anf*b1
        if (a1.ne.0.0) then
          fac=1.0/a1
          g=b1*fac
          if (abs((g-gold)/g).lt.eps) go to 1
          gold=g
        end if
      end do
      print '(A)','GCF: A too large, ITMAX too small.'
      stop
    1 t=-x+a*log(x)-gln
      if (t.lt.-50.0) then
        gammcf=0.0
      else
        gammcf=exp(t)*g
      end if
c      print *,'gcf: ',x,a,gln,-x+a*log(x)-gln,g
      end
c     ==================================================================
      function gammln(xx)
c
      double precision cof(6),stp,half,one,fpf,x,tmp,ser
      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
     *             -1.231739516d0,0.120858003d-2,-0.536382d-5,
     *             2.50662827465d0/
      data half,one,fpf/0.5d0,1.0d0,5.5d0/
c
      x=xx-one
      tmp=x+fpf
      tmp=(x+half)*log(tmp)-tmp
      ser=one
      do j=1,6
        x=x+one
        ser=ser+cof(j)/x
      end do
      gammln=tmp+log(stp*ser)
      end
