implicit real*8 (a-h,o-z)
      include 'DIMENSIONS'
      include 'COMMON.CHAIN'
      include 'COMMON.INTERACT'
      include 'COMMON.SBRIDGE'
      real*4 coord(3,2*maxres)
      real*4 prec,time,potE,uconst,t_bath,qfrag(100)
      real*8 etot
      real*8 cm(3,maxchain),shift_coord(3,maxchain),ccm(3),ccm_prev(3)
      character*80 arg,seqfile,pdbfile
      character*3 sequenc(maxres)
      character*50 tytul
      character*8 onethree,cfreq
      character*8 ucase
      external ucase
      logical oneletter
      integer rescode
      external rescode
      logical iblnk
      external iblnk
      logical previous,change
      integer lmin(maxchain),lshift(maxchain)
      
      previous=.false.
      ifreq=1
      is=1
      ie=1000000000
      if (iargc().lt.6) then
       print '(2a)',
     & "Usage: xdrf2pdb one/three seqfile boxx boxy boxz cxfile [freq]",
     &" [start] [end] [pdbfile]"
        stop
      endif
      call getarg(1,onethree)
      onethree = ucase(onethree)
      if (onethree.eq.'ONE') then
        oneletter = .true.
      else if (onethree.eq.'THREE') then
        oneletter = .false.
      else
        print *,"ONE or THREE must be specified"
      endif
      call getarg(2,seqfile)
      open (1,file=seqfile,status='old')
      if (oneletter) then
        read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
   10   continue
        nres=i
        i=0
        do while (.not.iblnk(sequenc(i+1)(1:1)))
          i=i+1
        enddo 
        nres=i
        do i=1,nres
          itype(i)=rescode(i,sequenc(i),1)
        enddo
      else
        read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
   11   continue
        nres=i
        i=0
        do while (.not.iblnk(sequenc(i+1)(1:1)))
          i=i+1
        enddo 
        nres=i
        do i=1,nres
          itype(i)=rescode(i,sequenc(i),0)
        enddo
      endif
      call seq2chains(nres,itype,nchain,chain_length,chain_border,
     &  ireschain)
      print *,"nchain",nchain
      print *,"chain_border, chain_length"
      nrestot=0
      do i=1,nchain
        nrestot=nrestot+chain_length(i)
        print *,chain_border(1,i),chain_border(2,i),chain_length(i)
      enddo
      print *,"nrestot",nrestot
      call getarg(3,arg)
      read(arg,*) boxsize(1)
      call getarg(4,arg)
      read(arg,*) boxsize(2)
      call getarg(5,arg)
      read(arg,*) boxsize(3)
      call getarg(6,arg)
      iext = index(arg,'.cx') - 1
      if (iext.lt.0) then
        print *,"Error - not a cx file"
        stop
      endif
      if (iargc().gt.6) then
        call getarg(7,cfreq)
        read (cfreq,*) ifreq
      endif
      if (iargc().gt.7) then
        call getarg(8,cfreq)
        read (cfreq,*) is
      endif
      if (iargc().gt.8) then
        call getarg(9,cfreq)
        read (cfreq,*) ie
      endif
      if (iargc().gt.9) then
        call getarg(10,pdbfile)
      else
        pdbfile=arg(:iext)//'.pdb'
      endif
      open(9,file=pdbfile)
      nnt = 1
      if (itype(1).eq.ntyp1) nnt = 2
      nct=nres
      if (itype(nres).eq.ntyp1) nct = nres-1
      print *,"nnt",nnt," nct",nct
      print *,"file",arg
      call xdrfopen(ixdrf,arg, "r", iret)
      if (iret.eq.0) stop
      print *,"iret",iret
      print *,"is",is," ie",ie
      kk = 0
      do while(is.eq.0 .or. kk.lt.ie) 
       call xdrffloat(ixdrf, time, iret)
       print *,"time",time," iret",iret
       if(iret.eq.0) exit
       kk = kk + 1
       call xdrffloat(ixdrf, potE, iret)
       call xdrffloat(ixdrf, uconst, iret)
       print *,"potE",potE," uconst",uconst
#ifdef NEWUNRES
       call xdrffloat(ixdrf, uconst_back, iret)
#endif
       print *,"uconst_back",uconst_back
       call xdrffloat(ixdrf, t_bath, iret)
       print *,"t_bath",t_bath
       call xdrfint(ixdrf, nss, iret) 
       do j=1,nss
        call xdrfint(ixdrf, ihpb(j), iret)
        call xdrfint(ixdrf, jhpb(j), iret)
       enddo
       print *,"nss",nss
       call xdrfint(ixdrf, nfrag, iret)
       do i=1,nfrag
        call xdrffloat(ixdrf, qfrag(i), iret)
       enddo
       print *,"nfrag",nfrag
       prec=10000.0

       isize=0
       call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
c       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
c       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
c       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
       if (kk.ge.is .and. mod(kk,ifreq).eq.0) then
         if (isize .ne. nres+nct-nnt+1) then
           print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
         endif
         do i=1,nres
           do j=1,3
             c(j,i)=coord(j,i)
           enddo
         enddo
         ii = 0
         do i=nnt,nct
           ii = ii + 1 
           do j=1,3
             c(j,i+nres)=coord(j,ii+nres)
           enddo
         enddo
         do ichain=1,nchain
c           print *,"ichain",ichain
           do j=1,3
             cm(j,ichain)=0.0d0
           enddo
           do i=chain_border(1,ichain),chain_border(2,ichain)
c             print *,i,c(:,i)
             do j=1,3
               cm(j,ichain)=cm(j,ichain)+c(j,i)
             enddo
           enddo
           do j=1,3
             cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
           enddo
         enddo
c         write (*,*) "CM"
c         do i=1,nchain
c           write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
c         enddo
         do j=1,3
           lmin=0
           sumodl_min=1.0d10
           do i=0,3**nchain-1
             ii=i
             do ichain=1,nchain
               lshift(ichain)=mod(ii,3)-1
               ii=ii/3
             enddo
c             print '(10i3)',(lshift(ichain),ichain=1,nchain)
             sumodl=0.0d0
             do ichain=1,nchain
               do jchain=1,ichain-1
               sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
     &             -cm(j,jchain)-lshift(jchain)*boxsize(j))
               enddo
             enddo
c             print *,"sumodl",sumodl," sumodl_min",sumodl_min
             if (sumodl.lt.sumodl_min) then
               sumodl_min=sumodl
               lmin=lshift
             endif
           enddo
           do ichain=1,nchain
             cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
             shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
           enddo
         enddo
         do ichain=1,nchain
           do j=1,3
             ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
           enddo
         enddo
         do j=1,3
           ccm(j)=ccm(j)/nrestot
         enddo
c         write (*,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
c         write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
c         write (*,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
c         write (*,*) "shift_coord"
c         do i=1,nchain
c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
c         enddo
#ifdef PREVIOUS
         if (previous) then
           do j=1,3
             change=.true.
             do while (change)
             if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
               do ichain=1,nchain
                 shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
               enddo
               ccm(j)=ccm(j)-boxsize(j)
             else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
               do ichain=1,nchain
                 shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
               enddo
               ccm(j)=ccm(j)+boxsize(j)
             else
               change=.false.
             endif
           enddo
           enddo
           do j=1,3
             ccm_prev(j)=ccm(j)
           enddo
         else
           previous=.true.
           do j=1,3
             ccm_prev(j)=ccm(j)
           enddo
         endif
#else
         do ichain=1,nchain
           do j=1,3
             shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
           enddo
         enddo 
#endif
c         write (*,*) "shift_coord"
c         do i=1,nchain
c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
c         enddo
         do ichain=1,nchain
           do i=chain_border(1,ichain),chain_border(2,ichain)
             do j=1,3
               c(j,i)=c(j,i)+shift_coord(j,ichain)
               c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
             enddo
           enddo
         enddo 
         etot=potE
         write (tytul,'(a,i6)') "Structure",kk
         call pdbout(etot,tytul,9)
       endif
      enddo
     
      end