dimension s1(3,10000),s2(3,10000),ca1(3,10000),ca2(3,10000)
      dimension xyz1(3,5000),xyz2(3,5000),waby(5000),ires(5000),
     & atnam(5000),rms(500,500),imin(5000),imax(5000),isel(5000),
     & isel_t(5000),issel(5000),diff(5000),iperm(5000)
      character*4 atnam
      dimension t(3),b(3,3),z(3,5000)
      integer resnum,resnum1,resnum2,resnum3,resnum4,frag(5000)
      character*80 karta
      character*80 fnam1,fnam2,fnam3
      logical insel
      ifrag=1
      nnn = iargc()
      if (nnn.lt.3) then
        stop 
     & 'Usage: suppdb pdb_file_1 pdb_file_2 out_file [resnum1 resnum2]'
      endif
      call getarg(1,fnam1)
      call getarg(2,fnam2)
      call getarg(3,fnam3)
      open (2,file=fnam3,status='unknown')
      if (nnn.lt.5) then
        print *,'resnum1,resnum2'
        read (*,*) resnum1,resnum2
      else
        call getarg(4,karta)
        read (karta,*) resnum1
        call getarg(5,karta)
        read (karta,*) resnum2
      endif
      if (nnn.gt.5) then
        call getarg(6,karta)
        read (karta,*) resnum3
      else
        resnum3=resnum1
      endif
      resnum4=resnum3+resnum2-resnum1
      print *,'resnum1',resnum1,' resnum2',resnum2
      print *,'resnum3',resnum3,' resnum4',resnum4
      print *,'opening file ',fnam1
      open (1,file=fnam1,status='old',err=111)
      iat=0
      isup=0
      write (2,'(/2a)') 'Coordinates of the first set from ',
     &  fnam1(:len_trim(fnam1))
      do irec=1,10000
      read (1,'(a80)',end=10) karta
      if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and.
     & karta(14:15).eq.'CA') then
      read (karta(23:26),*) resnum
      if (resnum.ge.resnum1 .and. resnum.le.resnum2) then
      iat=iat+1
      read (karta(14:16),'(a3)') atnam(iat)
      read (karta(31:54),'(3f8.3)') xyz1(1,iat),xyz1(2,iat),
     1                              xyz1(3,iat)
      read (karta(23:26),'(i4)') ires(iat)
      read (karta(18:20),'(a3)') waby(iat)
      if (atnam(iat).eq.'CA ') then
        write (2,'(a)') karta
        isup=isup+1
        do i=1,3
          ca1(i,isup)=xyz1(i,iat)
          s1(i,isup)=xyz1(i,iat)
        enddo
      endif
      endif
      endif
      enddo
   10 nsup=isup
      nat=iat
      print *,'from first file: nat=',nat
      close(1)
      print *,'opening file ',fnam2
      open (1,file=fnam2,status='old',err=112)
      iat=0
      isup=0
      write (2,'(/2a)') 'Coordinates of the first set from ',
     &  fnam2(:len_trim(fnam2))
      do irec=1,10000
      read (1,'(a80)',end=20) karta
      if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and.
     & karta(14:15).eq.'CA') then
      read (karta(23:26),*) resnum
      if (resnum.ge.resnum3 .and. resnum.le.resnum4) then
      iat=iat+1
      read (karta(31:54),'(3f8.3)') xyz2(1,iat),xyz2(2,iat),
     1                              xyz2(3,iat)
      if (atnam(iat).eq.'CA ') then
        write (2,'(a)') karta
        isup=isup+1
        do i=1,3
          ca2(i,isup)=xyz2(i,iat)
          s2(i,isup)=xyz2(i,iat)
        enddo
      endif
      endif
      endif
      enddo
   20 continue 
      print *,'from second file: nat=',iat
!      write (2,'(//4a8)') '# length','  start','   rms'


      do i=1,nsup
        issel(i)=i
      enddo

      nsup_t = nsup

      nsup_old=0

      DO WHILE (NSUP>4 .and. nsup.ne.nsup_old) 

      nsup_old = nsup

      do is1=4,nsup
!        write (2,'(a)')
        do is2=1,nsup-is1+1
          call fitsq(rms(is1,is2),s2(1,is2),s1(1,is2),is1,t,b)
!          write (2,'(2(i4,4x),2f8.2)') is1,is2+resnum1-1,rms(is1,is2)
        enddo
      enddo
      write (2,'(//2a8,3a8)') '# length',' start','rmsmin',
     &  'rmsmax','rmsave'
      do i=4,nsup
        imin(i)=1
        imax(i)=1
        rmin=rms(i,1)
        rmax=rms(i,1)
        rmsave=rms(i,1)**2
        do j=2,nsup-i+1
          if ( rms(i,j).lt.rmin) then
            imin(i)=j
            rmin=rms(i,j)
          else if ( rms(i,j).gt.rmax) then
            imax(i)=j
            rmax=rms(i,j)
          endif
          rmsave=rmsave+rms(i,j)*rms(i,j)
        enddo
        rmsave=sqrt(rmsave/(nsup-i+1))
        write (2,'(2(i4,4x),3(f8.2))') i,imin(i)+resnum1-1,
     &   rms(i,imin(i)), ! rms(i,imin(i))/sqrt(float(i)),
     &   rms(i,imax(i)), ! rms(i,imax(i))/sqrt(float(i)),
     &   rmsave ! ,rmsave/sqrt(float(i))
      enddo
      rmscut=7.0d0
      rmscut1=1.0d0*rmscut
      do i4=nsup,4,-1
        if (rms(i4,imin(i4)).le.rmscut) exit
      enddo    
      call fitsq(rrms,s2(1,imin(i4)),s1(1,imin(i4)),i4,t,b)
      nsup_t=nsup
      do i=1,nsup
        isel(i)=issel(i)
      enddo
      rrms=1000.0d0
c      do while (rrms.gt.rmscut)
      do while (rrms.gt.rmscut)
        nnsup_t=0
        std=0.0d0
        diffmax=-10000.0d0
        do i=1,nsup_t
          ii=isel(i)
          call mvvad(b,ca2(1,ii),t,z(1,i))
          diffi=0.0d0
          do j=1,3
            diffi=diffi+(z(j,i)-ca1(j,ii))**2
          enddo
          std=std+diffi
          diffi=sqrt(diffi)
!          write (2,*) i,ii,diffi
          diff(i)=diffi
          if (diffi.gt.diffmax) diffmax=diffi
        enddo
        do i=1,nsup_t
          iperm(i)=i
        enddo
        write(2,*) 'diff'
        write(2,*) (diff(i),i=1,nsup_t)
        call mysort(nsup_t,diff,iperm)
        write(2,*) 'diff'
        write(2,*) (diff(i),i=1,nsup_t)
        std=sqrt(std/nsup_t)
        write (2,*) 'std=',std
        write (*,*) 'std=',std
        sstd=dmin1(1.5d0*std,0.9999d0*diffmax)
!        do i=1,nsup_t
!          ii=isel(i)
!          if (diff(i).lt.sstd .or. diff(i).lt.rmscut) then
!            nnsup_t=nnsup_t+1
!            isel(nnsup_t)=ii
!          endif
!        enddo
        sumdiff=0.0d0
        do i=1,nsup_t
          ii=isel(iperm(i))
c          sumdiff=sumdiff+diff(i)**2
c          if (dsqrt(sumdiff/dfloat(i)).lt.rmscut1) then
          if (diff(i).lt.rmscut1) then
            nnsup_t=nnsup_t+1
            isel_t(nnsup_t)=ii
          endif
        enddo
        nsup_t=nnsup_t
        do i=1,nsup_t
          iperm(i)=i
          isel(i)=isel_t(i)
        enddo
        if (nsup_t.eq.0 .or. nnsup_t.eq.0) exit
        write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
        write (2,'(20i4)') (isel(i),i=1,nnsup_t)
        call mysort(nsup_t,isel,iperm)
        print *,'nsup_t',nsup_t,' nnsup_t',nnsup_t
        print '(20i4)',(isel(i),i=1,nnsup_t)
        write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
        write (2,'(20i4)') (isel(i),i=1,nnsup_t)
        len_min=4
        i=1
        len=1
        do i=1,nsup_t 
          ii=isel(i)
          do j=1,3
            s1(j,i)=ca1(j,ii)
            s2(j,i)=ca2(j,ii)
          enddo
        enddo
        do i=1,nsup_t
          ii=isel(i)
          do j=1,3
            s1(j,i)=ca1(j,ii)
            s2(j,i)=ca2(j,ii) 
          enddo
        enddo
        if (nsup_t.le.2) exit
        print *,"before fitsq: nsup_r",nsup_t
        call fitsq(rrms,s2(1,1),s1(1,1),nsup_t,t,b)
        write (2,*) 'nsup_t=',nsup_t,' rrms=',rrms
        write (*,*) 'nsup_t=',nsup_t,' rrms=',rrms
      enddo 
      i=1
      do while (i.le.nnsup_t)
        write (*,*) 'i=',i,' isel',isel(i),' isel1',isel(i+1)
        if (i.eq.nnsup_t .or. isel(i+1)-isel(i).gt.1) then
          if (len.lt.len_min) then
            write (*,'(20i4)') (isel(k),k=1,nnsup_t)
            write (*,*) 'len=',len,' len_min=',len_min
            do j=i+1,nnsup_t
              isel(j-len)=isel(j)
            enddo
            nnsup_t=nnsup_t-len
!            write (2,'(20i4)') (isel(k),k=1,nnsup_t)
            i=i-len+1
            len=1
          else
            len=1
            i=i+1
          endif
        else
          len=len+1
          i=i+1
        endif
!        write (2,*) 'i=',i,' len=',len
      enddo
      print *,'final nsup_t',nsup_t
      print *,"final isel"
      print '(20i4)',(isel(i),i=1,nnsup_t)
      write (2,*) 'final nsup_t',nsup_t
      write (2,*) "final isel"
      write (2,'(20i4)') (isel(i),i=1,nnsup_t)
      if (nnsup_t.gt.10) then
       do i=1,nnsup_t
        frag(isel(i))=ifrag
       enddo
       ifrag=ifrag+1
      endif
      print *,"----- end of iteration"
      write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
      nsup_t=nnsup_t
      write (2,'(/4a)') 'REMARK superimposed coordinates from files ',
     &fnam1(:len_trim(fnam1)),' and ',fnam2(:len_trim(fnam2))
      write (2,'(a,f6.2,a)') 
     & 'REMARK Calpha rms deviation:',rrms,
     & ' angstroms'
      write (2,'(a,I5)') 'REMARK number of superposed residues',nsup_t
      i1=isel(1)
      do i=2,nnsup_t
        if (isel(i)-isel(i-1).gt.1) then
          write (2,'(a,i4,a,i4,a,f5.1,a)') 
     &     'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
           i1=isel(i)
        endif
      enddo
      write (2,'(a,i4,a,i4,a,f5.1,a)') 
     &     'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
      do i=1,nsup_t
        ii=isel(i)
        write (2,100) ii+resnum1-1,atnam(ii),waby(ii),ires(ii),
     &                (xyz1(j,ii),j=1,3)
      enddo
      write (2,'(a)') 'TER'
      do i=1,nnsup_t
        ii=isel(i)
        call mvvad(b,xyz2(1,ii),t,z(1,ii))
        write (2,100) ii+resnum3-1,atnam(ii),waby(ii),ires(ii),
     &                (z(j,ii),j=1,3)
      enddo 
      write (2,'(a)') 'TER'

      ind=0
      do i=1,nsup
        insel=.false.
        do j=1,nsup_t
          if (issel(i).eq.isel(j)) then
            insel=.true.
            exit
          endif
        enddo
        if (.not.insel) then
          ind=ind+1 
          issel(ind)=issel(i)
        endif
      enddo
      nsup=nsup-nsup_t
      print *,"nsup",nsup
      if (nsup.gt.0) print *,"issel",(issel(i),i=1,nsup)
      do i=1,nsup
        do j=1,3
          s1(j,i)=ca1(j,issel(i))
          s2(j,i)=ca2(j,issel(i))
        enddo
      enddo

      ENDDO 

      write (2,'(5000i2)') (frag(i),i=1,resnum2-resnum1+1)
      stop 777
  111 print '(2a)','Error opening input file ',fnam1(:len_trim(fnam1))
      stop 111
  112 print '(2a)','Error opening input file ',fnam2(:len_trim(fnam2))
      stop 112
  100 format ('ATOM',i7,2x,a2,2x,a3,i6,4x,3f8.3)
      end
c----------------------------------------------------------------------
      subroutine distout(ncon,rms)
      parameter(ncol=10)
      dimension rms(200,200)
      dimension b(ncol)
      iout=2
c      icant(i,j)=((ncon+ncon-j)*(j-1))/2+i-j
      write (iout,'(a)') 'The distance matrix'
      do 1 i=1,ncon,ncol
      nlim=min0(i+ncol-1,ncon)
      write (iout,1000) (k,k=i,nlim)
      write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 1000 format (/8x,10(i4,3x))
 1020 format (/1x,80(1h-)/)
      do 2 j=i,ncon
      jlim=min0(j,nlim)
      if (jlim.eq.j) then
        b(jlim-i+1)=0.0d0
        jlim1=jlim-1
      else
        jlim1=jlim
      endif
      do 3 k=i,jlim1
    3 b(k-i+1)=rms(k,j)
      write (iout,1010) j,(b(k),k=1,jlim-i+1)
    2 continue
    1 continue
 1010 format (i5,3x,10(f6.2,1x))
      return
      end
!--------------------------------------------------
      subroutine mysort(n, x, ipermut)
      implicit none
      real x(n),xsort(n)
      integer ipermut(n)
      integer i,j,imax,ipm,n
      real xtemp
      do i=1,n
        xtemp=x(i)
        imax=i
        do j=i+1,n
          if (x(j).lt.xtemp) then
            imax=j
            xtemp=x(j)
          endif
        enddo
        x(imax)=x(i)
        x(i)=xtemp
        ipm=ipermut(imax)
        ipermut(imax)=ipermut(i)
        ipermut(i)=ipm
      enddo
      return
      end