! This progrum is for reduction of dynamical matrix.
! In slab model, bottom layers are fixted at bulk position.
! So, negative frequencies would appear.
! Due to this progrum, calculation range is limited in only surface.
!
! set natm
! set reduced_natm
!
      program reduce_dynamical_matrix
      implicit none

      integer, parameter :: natm=16
      integer, parameter :: reduced_natm= 6
      integer :: nimage, i_devi_atom(natm*6), i_atom_image(natm,natm*6)
      integer :: itmp(2)
      real(8) :: devi_atom(3,natm*6), force_atom(3,natm,natm*6)
      integer :: iimage, j, iat


      open(20,file='force.data')
      open(21,file='force.data_new')

      read(20,*) nimage, itmp(1), itmp(2)
      if ( nimage .ne. natm*6 ) then
        write(*,*) 'Stop !!'
        write(*,*) 'number of atoms is strange'
        write(*,*) 'number of atoms = ', natm
        write(*,*) 'number of images =', nimage
        stop
      endif

      do iimage=1,nimage
        read(20,*) i_devi_atom(iimage), (devi_atom(j,iimage),j=1,3)
        do iat=1,natm
          read(20,*) i_atom_image(iat,iimage),
     &              (force_atom(j,iat,iimage),j=1,3)
        end do
      end do

  777 format(i4,3e26.12)
      write(21,'(3i8)') reduced_natm*6, itmp(1), itmp(2)
      do iimage=1,reduced_natm*6
        write(21,777) i_devi_atom(iimage), (devi_atom(j,iimage),j=1,3)
        do iat=1,reduced_natm
          write(21,777) i_atom_image(iat,iimage),
     &                 (force_atom(j,iat,iimage),j=1,3)
        end do
      end do

      end