! PROGRAM TO CONVERT DCD TRAJECTORY FILE FROM NAMD TO A SERIES OF TINKER XYZ FILES



! MAIN PROGRAM
! WRITTEN BY JOSHUA D. MOORE, jdmoore@ncsu.edu
! NORTH CAROLINA STATE UNIVERSITY

      program dcdtoxyz

      implicit none

      include 'system.inc'

	integer i
	integer tinker_stride

      tinkerframe = 0

	write(*,*) "Enter name of dcd file you want to read"
	read(*,*) dcdfilename


! Initialize dcd file (read header information)
      call read_dcd(0)


      write(*,*) "Enter number of frames you want to jump in dcd file"
	read(*,*) nframeskip

	write(*,*) "The number of dcd frames remaining for use is",
     &            nframes - nframeskip

	write(*,*) "Enter the total number of frames the in dcd file you", 
     &           " want to use to write TINKER files"
	read(*,*) nframes_actual

	write(*,*) "How many dcd frames do you want between",
     &	           " writing Tinker files?"
	read(*,*) tinker_stride

      write(*,*) "Enter number of carbon atoms listed before fluid", 
     &            " atoms in dcd file"
	read(*,*) ncarbon
     
      write(*,*) "Enter number of fluid atoms"
      read(*,*) nummol

	write(*,*) "Enter name of Atom (e.g., Ar)"
	read(*,*) atomname

	write(*,*) "Enter type of Atom (e.g., 186)"
	read(*,*) atomtype

	write(*,*) "Do you want to wrap the coordinates before writing", 
     &            " (type yes or no)"

      read(*,*) wrapcoordinates


      if(wrapcoordinates.eq."yes") then

	   write(*,*) "Enter X,Y,Z Box Lengths"
	   read(*,*) boxx,boxy,boxz

	endif



! Skip over the number of frames at the beginning of file you dont want to write out

      call read_dcd(2)

! Read through the dcd file to write the Tinker xyz files for the number of frames you want

      do i = 1, nframes_actual


	   call read_dcd(1)

	   if(mod(frame,tinker_stride).eq.0) then

         call write_tinkerxyz

	   endif

	enddo 

      print*, "DONE WRITING TINKER FILES"

      end program





!_______________________________________________________________________
!                   SUBROUTINE TO OPEN AND READ DCD FILE
!_______________________________________________________________________

! READ IN INITIAL INFO
      subroutine read_dcd(switch)

      implicit none
      include 'system.inc'

      integer mxmols
      parameter(mxmols=99999)
      integer switch
      character*4 HDR 
      integer ICNTRL(20)
      integer NTITL
      character*80 TITLE(50)
      integer FREEAT(100)
      integer NFREAT
      real*8 XTLABC(6)
      real*4 x(mxmols),y(mxmols),z(mxmols)
      integer j,k,i
      integer nframecount

!_______________________________________________________________________
!OPEN DCD FILE AND READ HEADER
!_______________________________________________________________________

      if (switch.eq.0) then
        open(unit=8, file=dcdfilename, form = 'unformatted', 
     &  status= 'old', access= 'sequential')
        read(8) HDR, (ICNTRL(j),j=1,20)
        read(8) 
        read(8) totatoms
        NFREAT = totatoms - ICNTRL(9)
        print*, 'TOTAL ATOMS = ', totatoms
        print*, 'NUMBER OF MOLECULES = ', nummol
         nframes = ICNTRL(1)
         print*, "There are", nframes,"in dcd file"
!         print*, 'NUMBER OF FRAMES BETWEEN TIME ORIGINS = ', it0
!         print*, 'NUMBER OF TIME ORIGINS =', nframes/it0
         print*, 'DONE OPENING DCD'
!_______________________________________________________________________
!READ DCD FILE EACH STEP
!_______________________________________________________________________
      
          else if (switch.eq.1) then
          frame = frame + 1
          read(8) (XTLABC(j),j=1,6)
          read(8) (x(j),j=1,totatoms)
          read(8) (y(j),j=1,totatoms)
          read(8) (z(j),j=1,totatoms)
          do j = ncarbon+1, totatoms
             k = j - ncarbon
             xcenter(k) = (x(j))
             ycenter(k) = (y(j))
             zcenter(k) = (z(j))

          enddo


!SKIP A NUMBER OF SPECIFIED FRAMES IN DCD FILE (EQUILIBRATION TIME ETC)
!_____________________________________________________________________

      else if(switch.eq.2) then

      nframecount = 0

      do i = 1,nframeskip 
     
      nframecount = nframecount + 1

      read(8) (XTLABC(j),j=1,6)
      read(8) (x(j),j=1,totatoms)
      read(8) (y(j),j=1,totatoms)
      read(8) (z(j),j=1,totatoms)

      enddo

      print*, 'I AM SKIPPING ',nframecount,' FRAMES'

      endif

      return
      endsubroutine

!_______________________________________________________________________
!                     END OF SUBROUTINE TO READ DCD FILE
!_______________________________________________________________________


! SUBROUTINE TO WRITE TINKER XYZ FORMAT


      subroutine write_tinkerxyz


      implicit none

	include 'system.inc'

	integer i
	character*60 tinkerfilename


      tinkerframe = tinkerframe + 1

      if(wrapcoordinates.eq."yes") then

	do i = 1, nummol

         xcenter(i) = xcenter(i) - boxx*anint(xcenter(i)/boxx)
         ycenter(i) = ycenter(i) - boxy*anint(xcenter(i)/boxy)
	   zcenter(i) = zcenter(i) - boxz*anint(xcenter(i)/boxz)

      enddo

	endif


      if (tinkerframe.lt.10) then
      write(tinkerfilename,'(A9,I1)')  'TINKER.00',tinkerframe

      elseif (tinkerframe.lt.100.and.tinkerframe.gt.9) then
      write(tinkerfilename,'(A8,I2)')  'TINKER.0',tinkerframe

      elseif (tinkerframe.lt.1000.and.tinkerframe.gt.99) then
      write(tinkerfilename,'(A7,I3)')  'TINKER.',tinkerframe

      elseif (tinkerframe.lt.10000.and.tinkerframe.gt.999) then
      write(tinkerfilename,'(A7,I4)')  'TINKER.',tinkerframe

      elseif (tinkerframe.lt.100000.and.tinkerframe.gt.9999) then
      write(tinkerfilename,'(A7,I5)')  'TINKER.',tinkerframe

      elseif (tinkerframe.lt.1000000.and.tinkerframe.gt.99999) then
      write(tinkerfilename,'(A7,I6)')  'TINKER.',tinkerframe

	endif

      open(unit=21,file=tinkerfilename,status="new")

	write(21,'(I6)') nummol
      
      do i = 1, nummol

	write(21,'(I6,2X,A2, F13.6,F12.6,F12.6,I6)') i,atomname,
     &                                               xcenter(i),
     &                                               ycenter(i),
     &                                               zcenter(i),
     &                                               atomtype



      enddo

	close(21)

     
      return

	endsubroutine