program xdrfpdb_m
!      implicit real*8 (a-h,o-z)
!      include 'DIMENSIONS'
!      include 'COMMON.CHAIN'
!      include 'COMMON.INTERACT'
!      include 'COMMON.SBRIDGE'
      use geometry_data, only: nres,c,nfrag,boxxsize,boxysize,boxzsize
      use energy_data, only: itype,nnt,nct,nss,ihpb,jhpb,uconst_back,molnum
      use io_base!, only: maxres,ucase,iblnk,rescode,pdbout
!      use geometry, only: returnbox
      implicit none
      real(kind=4),allocatable,dimension(:,:) :: coord !(3,2*maxres)
      real(kind=4) :: prec,time,potE,uconst,t_bath,qfrag(100)
      real(kind=8) :: etot
      character(len=80) :: arg,seqfile,pdbfileX,boxtraj
      character(len=3),allocatable,dimension(:) :: sequenc !(maxres)
      character(len=50) :: tytul
      character(len=8) :: onethree,cfreq,cntraj,citraj
      character(len=3) :: licz
!      external ucase
      logical oneletter !,iblnk
!      integer rescode
!      external rescode
      integer :: i,ii,j,k,kk,is,ie,ifreq,mnum,molec,iext,isize
      integer :: iret,ixdrf,itraj,ntraj
      print *,"begin program"
!      if(.not.allocated(molnum)) allocate(molnum(maxres))
      ifreq=1
      nres=0
      mnum=1
      molec=1
      allocate(sequenc(maxres))
      allocate(molnum(maxres))
      allocate(itype(maxres,5))
      if (iargc().lt.6) then
        print '(2a)',&
         "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",&
         "boxsize [pdbfile] [freq]"
        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,mnum)=rescode(i,sequenc(i),1,molec)
        enddo
      else
        read(1,'(20(a3,i1))',end=11,err=11) (sequenc(i),molnum(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,molnum(i))=rescode(i,sequenc(i),0,molnum(i))
        enddo
        print *,nres
        print '(20(a3,1x))',(sequenc(i),i=1,nres)
      endif
      call getarg(3,arg)
      iext = index(arg,'.cx') - 1
      if (iext.lt.0) then
        print *,"Error - not a cx file"
        stop
      endif
      call getarg(4,cntraj)
      read (cntraj,*) ntraj
      call getarg(5,citraj)
      read (citraj,*) itraj
      call getarg(6,boxtraj)
      read (boxtraj,*) boxxsize
      boxysize=boxxsize
      boxzsize=boxxsize
      if (iargc().gt.6) then
        call getarg(7,pdbfileX)
      else
        write(licz,'(bz,i3.3)') itraj
        pdbfileX=arg(:iext)//'_'//licz//'.pdb'
      endif
      if (iargc().gt.7) then
        call getarg(8,cfreq)
        read (cfreq,*) ifreq
      endif
!      print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
      open(9,file=pdbfileX)
      nnt = 1
      mnum=molnum(1)
      if (itype(1,mnum).eq.ntyp1) nnt = 2
      nct=nres
      mnum=molnum(nct)
      if (itype(nres,mnum).eq.ntyp1) nct = nres-1
      print *,"nnt",nnt," nct",nct
      call xdrfopen(ixdrf,arg, "r", iret)
      kk = 0
      allocate(coord(3,2*nres))
      allocate(c(3,2*nres))
      do while(.true.) 
       call xdrffloat(ixdrf, time, iret)
       if(iret.eq.0) exit
       kk = kk + 1
       call xdrffloat(ixdrf, potE, iret)
       call xdrffloat(ixdrf, uconst, iret)
!#ifdef NEWUNRES
!       call xdrffloat(ixdrf, uconst_back, iret)
!#endif
       call xdrffloat(ixdrf, t_bath, iret)
#ifdef NEWUNRES
       call xdrffloat(ixdrf, uconst_back, iret)
#endif

       call xdrfint(ixdrf, nss, iret) 
       do j=1,nss
        call xdrfint(ixdrf, ihpb(j), iret)
        call xdrfint(ixdrf, jhpb(j), iret)
       enddo
       call xdrfint(ixdrf, nfrag, iret)
       do i=1,nfrag
        call xdrffloat(ixdrf, qfrag(i), iret)
       enddo
       prec=10000.0

       isize=0
       call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)

       print *,"search here"
       write (*,'(e15.8,2e15.5,2f12.5)') time,potE,uconst,t_bath,&
        uconst_back
!       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
!       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
       do j=1,isize
       write (*,'(3f10.5,i3)') (coord(k,j),k=1,3),j
       enddo
       if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) 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
         call returnbox
         etot=potE
!         write (tytul,'(a,i6)') "Structure",kk
         write (tytul,'(a,i6,a,f8.3)') "Structure",kk,"Temp",t_bath
         call pdbout(etot,tytul,9)
       endif
      enddo
     
      end
      subroutine returnbox
      use geometry_data, only: nres,c,nfrag,boxxsize,boxysize,boxzsize
      use energy_data, only: itype,nnt,nct,nss,ihpb,jhpb,uconst_back,molnum
      use io_base!, only: maxres,ucase,iblnk,rescode,pdbout
      integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
      integer :: chain_end,ireturnval
      real*8 :: difference,xi,boxsize,x,xtemp,box2shift

       do i=1,nres
         if (molnum(i).eq.5) then
          c(1,i)=dmod(c(1,i),boxxsize)
!          if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
          x=c(1,i)-c(1,2)
          print *,"return box,before wrap",c(1,i)
         boxsize=boxxsize
           xtemp=dmod(x,boxsize)
      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
        box2shift=xtemp-boxsize
      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
        box2shift=xtemp+boxsize
      else
        box2shift=xtemp
      endif
          xi=box2shift
!          print *,c(1,2),xi,xtemp,box2shift
          c(1,i)=c(1,2)+xi
          print *,"after",c(1,2),xi,xtemp,box2shift,c(1,i)

          c(2,i)=dmod(c(2,i),boxysize)
!          if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize
          x=c(2,i)-c(2,2)
         boxsize=boxysize
           xtemp=dmod(x,boxsize)
      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
        box2shift=xtemp-boxsize
      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
        box2shift=xtemp+boxsize
      else
        box2shift=xtemp
      endif
          xi=box2shift
          c(2,i)=c(2,2)+xi
!         c(3,i)=dmod(c(3,i),boxzsize)
          if (c(3,i).lt.0) c(3,i)=c(3,i)+boxzsize
          x=c(3,i)-c(3,2)
         boxsize=boxzsize
           xtemp=dmod(x,boxsize)
      if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
        box2shift=xtemp-boxsize
      else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
        box2shift=xtemp+boxsize
      else
        box2shift=xtemp
      endif
          xi=box2shift
          c(3,i)=c(3,2)+xi

          c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
          c(2,i+nres)=dmod(c(2,i+nres),boxysize)
          c(3,i+nres)=dmod(c(3,i+nres),boxzsize)
         endif
        enddo
        return
        end       subroutine returnbox