!LIGGGHTSの計算結果の経時変化をgnuplot用に項目別にまとめる
      ! t1 atom1 atom2 atom3 ……
      ! t2 atom1 atom2 atom3 ……
      ! t3 atom1 atom2 atom3 ……
      !  :   :     :     :   \
      !といった具合にpost/に項目ごとにファイルを作ってまとめる

!出力するデータとかファイルは適宜追加してください
!配列は(粒子数,出力ステップ数)なので、適宜配列を書き換えたり、
!出力ステップ数＝run/dumpなので、入力ファイルのdumpの値を大きくしてください。
!ここではy,z,xyz,vvの4ファイルのみ出力

!2345678901234567890
      implicit none
      real ::  timestep,dump,run
      integer :: i,j,k,l,o
      character :: filename*15,filenamename*30,m*10
      !配列は(粒子数,出力ステップ数)なので、適宜増やしてください
      integer :: atomnum,id(10000,1000),atomtype(10000,1000)
      real :: x(10000,1000),y(10000,1000)
      real :: z(10000,1000),ix(10000,1000)
      real :: iy(10000,1000),iz(10000,1000)
      real :: vx(10000,1000),vy(10000,1000)
      real :: vz(10000,1000),fx(10000,1000)
      real :: fy(10000,1000),fz(10000,1000)
      real :: omegax(10000,1000),omegay(10000,1000)
      real :: omegaz(10000,1000),radius(10000,1000)
      real :: t(1000),vv(10000,1000),xyz(10000,1000)

!LIGGGHTSの入力ファイルに従い、書き換える!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      timestep=0.00001 !1ステップで計算する時間（秒）
      dump=100 !出力ファイルに記録するステップの間隔
      run=100000 !全部で計算するステップ数
      filename="saka" !LIGGGHTSで計算させた結果のファイル名
!post/以下を分かりやすい名前に書き換える、書き込み専用!!!!!!!!!!!!!!!!!!!!!!!!!!!
      open(2,file=('post/tz.dat'),status='replace') !グラフt-zの名前
      open(3,file=('post/tvv.dat'),status='replace') !グラフt-vvの名前
      open(4,file=('post/ty.dat'),status='replace') !グラフt-yの名前
      open(5,file=('post/txyz.dat'),status='replace') !グラフt-xyzの名前
!出力ファイルを増やしたかったら、上に「open(n,~」と一番下に「close(n)」を書き足す

      l=run/dump
      do i=1,l
      j=i*dump
      write(m,"(i10)") j
      filenamename="post/dump"//trim(adjustl(m))//"."//filename
      open(1,file=(filenamename),status='old') !読み込み専用

      do o =1,3
      read(1,'()')
      end do
      read(1,*) atomnum
      do o=1,5
      read(1,'()')
      end do

      t(i)=timestep*j !実際の経過時間(sec)の計算
      
      do k=1,atomnum
      read(1,*) id(k,i),atomtype(id(k,i),i)&
          &,x(id(k,i),i),y(id(k,i),i),z(id(k,i),i)&
          &,ix(id(k,i),i),iy(id(k,i),i),iz(id(k,i),i)&
          &,vx(id(k,i),i),vy(id(k,i),i),vz(id(k,i),i)&
          &,fx(id(k,i),i),fy(id(k,i),i),fz(id(k,i),i)&
          &,omegax(id(k,i),i),omegay(id(k,i),i),omegaz(id(k,i),i)&
          &,radius(id(k,i),i)
      
      vv(id(k,i),i)=sqrt((vx(id(k,i),i))**2+(vy(id(k,i),i))**2&
          &+(vz(id(k,i),i))**2) !速度の合成
      xyz(id(k,i),i)=sqrt((x(id(k,i),i)-x(id(k,i),1))**2&
          &+(y(id(k,i),i)-y(id(k,i),1))**2&
          &+(z(id(k,i),i)-z(id(k,i),1))**2) !座標の合成(初期位置からどれだけ動いたか)
      end do
      
      write(2,'(1pe13.5,a,$)') t(i)," "
      write(3,'(1pe13.5,a,$)') t(i)," "
      write(4,'(1pe13.5,a,$)') t(i)," "
      write(5,'(1pe13.5,a,$)') t(i)," "
      do k=1,atomnum-1
!出力したいデータに書き換える!!!!!!!!!!!!!!!!
      write(2,'(1pe13.5,a,$)') z(k,i)," "
      write(3,'(1pe13.5,a,$)') vv(k,i)," "
      write(4,'(1pe13.5,a,$)') y(k,i)," "
      write(5,'(1pe13.5,a,$)') xyz(k,i)," "
      end do
      write(2,'(1pe13.5)') z(atomnum,i)
      write(3,'(1pe13.5)') vv(atomnum,i)
      write(4,'(1pe13.5)') y(atomnum,i)
      write(5,'(1pe13.5)') xyz(atomnum,i)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      end do

      close(1)
      close(2)
      close(3)
      close(4)
      close(5)
      end
