#!/bin/bash

###################################################################
#                                                                 #
#   Parse a DCD Trajectory                                        #
#   Requires:  1) dcd   2) psf					  #
#                                                                 #
###################################################################


PS3="Your choice: "
QUIT="QUIT THIS PROGRAM!"

echo " "
echo "Which DCD Trajectory?"
echo " "
touch "$QUIT"
select dcd in *.dcd QUIT*;
do
  case $dcd in
        "$QUIT")
          rm QUIT*
          echo "Exiting."
          exit
          ;;
        *)
          rm QUIT*
          echo "You picked $dcd ($REPLY)"
          ;;
  esac
  break
done



echo " "
echo "Which PSF?"
echo " "
touch "$QUIT"
select psf in *.psf QUIT*;
do
  case $psf in
        "$QUIT")
          rm QUIT*
          echo "Exiting."
          exit
          ;;
        *)
          rm QUIT*
          echo "You picked $psf ($REPLY)"
          ;;
  esac
  break
done

natom=`cat $psf | grep NATOM | awk '{print $1}'`

cat $psf | sed -n '/NATOM/,/NBOND/p' | grep -v "NATOM\|NBOND" | grep -v '^$' | awk '{print $3, $4, $7, $8}' > info

ln -s $dcd tmp.dcd


cat <<HERE > inp1.f
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program readdcd

      character*4 header
      integer control(20)
      integer i,j,ntitle,natom
      character*80 title(4)

c     open input trajectory file
      open(unit=56,file="tmp.dcd",status='old',form='unformatted')

c     READ IN HEADER AND CONTROL INFO
      read(56) header,control
      write (*,*) " "
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) "TRAJECTORY INFORMATION"
      write (*,*) " "
      write (*,*) "Number of Frames:                 ",control(1)
      write (*,*) "Number of Previous Dynamic Steps: ",control(2)
      write (*,*) "Skip Frequency:                   ",control(3)
      write (*,*) "Number of Dynamics Steps:         ",control(4)

c     read in the number of title lines
      read(56) ntitle,(title(j),j=1,ntitle)

c     read in the number of atoms
      read(56) natom
      write (*,*) "Number of Atoms:                  ",natom
      write (*,*) " "
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     
      write (*,*) " "     

      close(56)
   
      end program
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
HERE
g77 inp1.f -o a1.out
./a1.out 


echo "Enter First Frame"
echo " "
read bframe
echo "Enter Last Frame"
echo " "
read eframe
echo "Enter Skip"
echo " "
read skip




cat <<HERE > inp2.f
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      program readdcd

      character*4 header
      character*4 resname($natom)
      integer control(20),resnum($natom)
      integer i,j,ntitle,natom
      character*80 title(4)
      real*4 x($natom)
      real*4 y($natom)
      real*4 z($natom)
      double precision xyz($natom,3),xyz_avg(3)
      double precision m($natom),ch($natom)
      double precision com(3),moi(3),rgy
      
c     OPEN INFO FILE
      open(unit=55,file='info',status='old')
      do i=1,$natom
        read(55,*) resnum(i),resname(i),ch(i),m(i)
      enddo
      close(55)
c     OPEN TRAJECTORY FILE
      open(unit=56,file="tmp.dcd",status='old',form='unformatted')
c     READ IN HEADER AND CONTROL INFO
      read(56) header,control
c     READ IN THE TITLE LINES
      read(56) ntitle,(title(j),j=1,ntitle)
c     READ IN NUMBER OF ATOMS
      read(56) natom
c     READ IN COORDINATES AND CONVERT TO DOUBLE PRECISION
      do i=1,$eframe
        read(56) x
        read(56) y
        read(56) z
        if (mod(i,$skip).eq.0) then
          do j=1,natom
            xyz(j,1)=x(j)
            xyz(j,2)=x(j)
            xyz(j,3)=x(j)
          enddo
c         CALCULATE PHYSICAL QUANTITIES
          call mean_position(natom,xyz,xyz_avg)
          write (*,*) "    Mean Position:  ",xyz_avg
          call center_of_mass(natom,m,xyz,com)
          write (*,*) "   Center of Mass:  ",com
          call moment_of_inertia(natom,m,xyz,moi)
          write (*,*) "Moment of Inertia:  ",moi
          call radius_of_gyration(natom,xyz,xyz_avg,rgy)
          write (*,*) "Radius of Gyration: ",rgy
          write (*,*) " "
        endif
      enddo
      close(56)
      end program

 

      subroutine mean_position(natom,xyz,xyz_avg)
      integer i,natom
      double precision xyz(natom,3),xyz_avg(3)
        xyz_avg(1)=0
        xyz_avg(2)=0
        xyz_avg(3)=0
        do i=1,natom
          xyz_avg(1)=xyz_avg(1)+xyz(i,1)    
          xyz_avg(2)=xyz_avg(1)+xyz(i,2)    
          xyz_avg(3)=xyz_avg(1)+xyz(i,3)    
        enddo
        xyz_avg(1)=xyz_avg(1)/natom
        xyz_avg(2)=xyz_avg(2)/natom
        xyz_avg(3)=xyz_avg(3)/natom
      end subroutine



      subroutine center_of_mass(natom,m,xyz,com)
      integer i,natom
      double precision m(natom),com(3),mtot
      double precision xyz(natom,3)
      mtot=0
      com(1)=0
      com(2)=0
      com(3)=0
      do i=1,natom
        mtot=mtot+m(i)
        com(1)=com(1)+m(i)*xyz(i,1)
        com(2)=com(2)+m(i)*xyz(i,2)
        com(3)=com(3)+m(i)*xyz(i,3)
      enddo
        com(1)=com(1)/mtot
        com(2)=com(2)/mtot
        com(3)=com(3)/mtot
      end subroutine



      subroutine moment_of_inertia(natom,m,xyz,moi)
        integer i
        double precision m(natom),moi(3)
        double precision xyz(natom,3)
        moi(1)=0
        moi(2)=0
        moi(3)=0
        do i=1,natom
          moi(1) = moi(1) + m(i) * (xyz(i,2)**2+xyz(i,3)**2)
          moi(2) = moi(2) + m(i) * (xyz(i,1)**2+xyz(i,3)**2)
          moi(3) = moi(3) + m(i) * (xyz(i,1)**2+xyz(i,2)**2)
        enddo
      end subroutine



      subroutine radius_of_gyration(natom,xyz,xyz_avg,rgy)
        integer i,natom
        double precision xyz(natom,3),xyz_avg(3)
        double precision r,rmean,rgy,rgy2
        rgy=0
        rgy2=0
        rmean=(xyz_avg(1)**2+xyz_avg(2)**2+xyz_avg(3)**2)**0.5
        do i=1,natom
          r=(xyz(i,1)**2+xyz(i,2)**2+xyz(i,3)**2)**0.5 
          rgy2=rgy2+(r-rmean)**2
        enddo
        rgy2=rgy2/natom
        rgy=rgy2**0.5
      end subroutine

c
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
HERE

g77 inp2.f -o a2.out
./a2.out 



rm inp1.f
rm inp2.f
rm a1.out
rm a2.out
rm info
rm tmp.dcd

exit