double precision function rmscalc_frag(ishif,i,j,jcon,kkk,
     &  lprn)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.COMPAR'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.CONTROL'
      double precision przes(3),obrot(3,3)
      double precision creff(3,maxres2),cc(3,maxres2)
      logical iadded(maxres)
      integer inumber(2,maxres)
      common /ccc/ creff,cc,iadded,inumber
      logical lprn
      logical non_conv
      integer ishif,i,j
      integer kkk
      if (lprn) then
        write (iout,*) "i",i," j",j," jcont",jcon," ishif",ishif
        write (iout,*) "npiece",npiece(j,i)
        write (iout,*)  "kkk",kkk
        call flush(iout)
      endif
c      write (iout,*) "symetr",symetr
c      call flush(iout)
c      nperm=1
c      do idup=1,symetr
c      nperm=nperm*idup
c      enddo
c      write (iout,*) "nperm",nperm
c      call flush(iout)
c      do kkk=1,nperm
      idup=0
      do l=1,nres
        iadded(l)=.false.
      enddo
c      write (iout,*) "kkk",kkk
c      call flush(iout)
      do k=1,npiece(j,i)
        if (i.eq.1) then
          if (lprn) then
            write (iout,*) "Level 1: j=",j,"k=",k," adding fragment",
     &         ifrag(1,k,j),ifrag(2,k,j)
            call flush(iout)
          endif
          call cprep(ifrag(1,k,j),ifrag(2,k,j),ishif,idup,kkk)
c          write (iout,*) "Exit cprep"
c          call flush(iout)
c          write (iout,*) "ii=",ii
        else
          kk = ipiece(k,j,i)
c          write (iout,*) "kk",kk," npiece",npiece(kk,1)
          do l=1,npiece(kk,1)
            if (lprn) then
              write (iout,*) "Level",i,": j=",j,"k=",k," kk=",kk,
     &          " l=",l," adding fragment",
     &          ifrag(1,l,kk),ifrag(2,l,kk)
              call flush(iout)
            endif
            call cprep(ifrag(1,l,kk),ifrag(2,l,kk),ishif,idup,kkk)
c            write (iout,*) "After cprep"
c            call flush(iout)
          enddo 
        endif
      enddo
      if (lprn) then
        write (iout,*) "tuszukaj"
c        do kkk=1,nperm
          do k=1,idup
            write(iout,'(5i4,2(3f10.5,5x))') i,j,k,inumber(1,k),
     &        inumber(2,k),(creff(l,k),l=1,3),(cc(l,k),l=1,3)
          enddo
c        enddo
        call flush(iout)
      endif
c      rminrms=1.0d10
c      do kkk=1,nperm
      call fitsq(rms,cc(1,1),creff(1,1),idup,przes,obrot,non_conv)
      if (non_conv) then
        print *,'Error: FITSQ non-convergent, jcon',jcon,i
        rms = 1.0d10
      else if (rms.lt.-1.0d-6) then 
        print *,'Error: rms^2 = ',rms,jcon,i
        rms = 1.0d10
      else if (rms.ge.1.0d-6 .and. rms.lt.0) then
        rms = 0.0d0
      endif
c      write (iout,*) "rmsmin", rminrms, "rms", rms
c      if (rms.le.rminrms) rminrms=rms
c      enddo
      rmscalc_frag = dsqrt(rms)
c      write (iout, *) "analysys", rmscalc,anatemp
      return
      end
c-------------------------------------------------------------------------
      subroutine cprep(if1,if2,ishif,idup,kwa)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.CONTROL'
      include 'COMMON.IOUNITS'
      include 'COMMON.COMPAR'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      double precision przes(3),obrot(3,3)
      double precision creff(3,maxres2),cc(3,maxres2)
      logical iadded(maxres)
      integer inumber(2,maxres),iistrart,kwa,blar
      common /ccc/ creff,cc,iadded,inumber
      integer ll,iperm
c      write (iout,*) "Calling cprep if1",if1," if2",if2," ishif",ishif,
c     &    " kwa",kwa
c      nperm=1
c      do blar=1,symetr
c      nperm=nperm*blar
c      enddo
c      write (iout,*) "nperm",nperm
c      kkk=kwa
c      ii=0
      do l=if1,if2
c        write (iout,*) "l",l," iadded",iadded(l)," ireschain",
c     &     ireschain(l),ireschain(l+ishif)
c        call flush(iout)
        if (l+ishif.gt.1 .and. l+ishif.le.nres .and. .not.iadded(l) 
     &   .and. ireschain(l+ishif).gt.0 .and. ireschain(l).gt.0 .and.
     &   ireschain(l).eq.ireschain(l+ishif)) then
          idup=idup+1
          iadded(l)=.true.
          inumber(1,idup)=l
          inumber(2,idup)=l+ishif
          ll=iperm(l+ishif,kwa)
          do m=1,3
            creff(m,idup)=cref(m,l)
            cc(m,idup)=c(m,ll)
          enddo
c          write (iout,'(2i5,3f10.5,5x,3f10.5)') l,ll,
c     &     (creff(m,idup),m=1,3),(cc(m,idup),m=1,3)
        endif
      enddo
c      write (iout,*) "idup",idup
      return
      end
c-------------------------------------------------------------------------
      double precision function rmsnat(jcon,ipermmin)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'DIMENSIONS.ZSCOPT'
      include 'DIMENSIONS.COMPAR'
      include 'COMMON.IOUNITS'
      include 'COMMON.COMPAR'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      include 'COMMON.CONTROL'
      integer ipermmin
      rmsnat = rmscalc(c(1,1),cref(1,1),ipermmin)
      return
      end
c-----------------------------------------------------------------------------
      double precision function rmscalc(ccc,cccref,ipermmin)
      implicit none
      include 'DIMENSIONS'
      include 'COMMON.IOUNITS'
      include 'COMMON.CHAIN' 
      include 'COMMON.INTERACT'
      include 'COMMON.VAR'
      double precision cccref(3,maxres2),creff(3,maxres2),
     &  ccc(3,maxres2),cc(3,maxres2)
      double precision przes(3),obrot(3,3)
      logical non_conv
      integer i,ii,j,ib,ichain,indchain,ichain1,ichain2,
     &  iperm,ipermmin
      double precision rms,rmsmin
C Loop over chain permutations
c      write (iout,*) "iz_sc",iz_sc
      rmsmin=1.0d10
      DO IPERM=1,NPERMCHAIN
c      write (iout,*) "iperm",iperm
      ii=0
      if (iz_sc.lt.2) then
        do ichain=1,nchain
          indchain=tabpermchain(ichain,iperm)
#ifdef DEBUG
          write (iout,*) "ichain",ichain," indchain",indchain
          write (iout,*) "chain_border",chain_border(1,ichain),
     &      chain_border(2,ichain)
#endif
          do i=1,chain_length(ichain)
c          do i=nstart_sup(ichain),nend_sup(ichain)
            ichain1=chain_border(1,ichain)+i-1
            ichain2=chain_border(1,indchain)+i-1
            if (ichain1.lt.nstart_sup .or. ichain1.gt.nend_sup .or.
     &          ichain2.lt.nstart_sup .or. ichain2.gt.nend_sup) cycle
            ii=ii+1
#ifdef DEBUG
            write (iout,*) "back",ii," ichain1",ichain1,
     &       " ichain2",ichain2," i",i,chain_border(1,ichain)+i-1
#endif
            do j=1,3
              cc(j,ii)=ccc(j,ichain2)
              creff(j,ii)=cccref(j,ichain1)
            enddo
#ifdef DEBUG
            write (iout,'(2i5,3f10.5,5x,3f10.5)')
     &       ichain1,ii,(cc(j,ii),j=1,3),(creff(j,ii),j=1,3)
#endif
          enddo
        enddo
      endif
      if (iz_sc.gt.0) then
        do ichain=1,nchain
          indchain=tabpermchain(ichain,iperm)
          do i=1,chain_length(ichain)
c          do i=nstart_sup(ichain),nend_sup(ichain)
            ichain1=chain_border(1,ichain)+i-1
            ichain2=chain_border(1,indchain)+i-1
            if (ichain1.lt.nstart_sup .or. ichain1.gt.nend_sup .or.
     &          ichain2.lt.nstart_sup .or. ichain2.gt.nend_sup) cycle
            if (itype(ichain1).ne.10) then
              ii=ii+1
#ifdef DEBUG
              write (iout,*) "side",ii," ichain1",ichain1,
     &         " ichain2",ichain2
#endif
              do j=1,3
                cc(j,ii)=ccc(j,ichain2+nres)
                creff(j,ii)=cccref(j,ichain1+nres)
              enddo
#ifdef DEBUG
              write (iout,'(2i5,3f10.5,5x,3f10.5)') 
     &        ichain1+nres,ii,(cc(j,ii),j=1,3),(creff(j,ii),j=1,3)
#endif
            endif
          enddo
        enddo
      endif
c      write (iout,*) "rmscalc: iprot",iprot," nsup",nsup(iprot)," ii",ii
      call fitsq(rms,cc(1,1),creff(1,1),ii,przes,obrot,non_conv)
      if (non_conv) then
        write (iout,*) 'Error: FITSQ non-convergent'
        rms=1.0d2
      else if (rms.lt.-1.0d-6) then 
        print *,'Error: rms^2 = ',rms
        rms = 1.0d2
      else if (rms.ge.1.0d-6 .and. rms.lt.0) then
        rmscalc=0.0d0
      else 
        rms = dsqrt(rms)
      endif
      if (rms.lt.rmsmin) then
        rmsmin=rms
        ipermmin=iperm
      endif
#ifdef DEBUG
      write (iout,*) "iperm",iperm," rms",rms
#endif
      ENDDO
      rmscalc=rmsmin
#ifdef DEBUG
      write (iout,*) "ipermmin",ipermmin," rmsmin",rmsmin
#endif
      return
      end
c-----------------------------------------------------------------------------
      double precision function gyrate(jcon)
      implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.INTERACT'
      include 'COMMON.CHAIN'
      include 'COMMON.IOUNITS'
      double precision cen(3),rg

      do j=1,3
       cen(j)=0.0d0
      enddo

      ii=0
      do i=nnt,nct
        if (itype(i).eq.ntyp1) cycle
        ii=ii+1
        do j=1,3
          cen(j)=cen(j)+c(j,i)
        enddo
      enddo
      do j=1,3
        cen(j)=cen(j)/dble(ii)
      enddo
      rg = 0.0d0
      do i = nnt, nct
        if (itype(i).eq.ntyp1) cycle
        do j=1,3
         rg = rg + (c(j,i)-cen(j))**2
        enddo
      end do
      gyrate = dsqrt(rg/dble(ii))
      return
      end