implicit double precision(a-h,o-z)
      real*8 v(100000,3),fa(100000,3),fb(100000,3),x(100000,3)
      real*8 av(3),sv(3),dx(3)

      open( 1,file='xe_ini_1.d',status='unknown')
      open( 3,file='xe_ini_2.d',status='unknown')
      open(10,file='xe_moni_1.d',status='unknown')

      wkb = 1.38065812d-23
      avo = 6.02213673d+23
c read initial data
      read (1,*) npart,nstep,sigma,e,wm,cell,tset,dt
      do i=1,npart
        read(1,*) x(i,1),x(i,2),x(i,3)
        read(1,*) v(i,1),v(i,2),v(i,3)
      enddo
      read(1,*) time0
c remove momentum of the cell
      sv=0.d0
      do i=1,npart
        sv(:) = sv(:) + v(i,:)
      enddo
      av(:) = sv(:)/npart
      do i=1,npart
        v(i,:) = v(i,:) - av(:)
      enddo
c main calculation
      cellh = cell * 0.5d0
      call force(fb,cell,cellh,x,sigma,e,u,npart)
c main loop
      do k=1,nstep
        call cpu_time(t1)
        time0 = dt + time0
        do i=1,npart
          v(i,:)=v(i,:)+dt*fb(i,:)/(2.0d0*wm)
          x(i,:)=x(i,:)+dt*v(i,:)
        enddo
c periodic boundary condition
        do i=1,npart
          do l=1,3
            dx(l)=x(i,l)-cellh
            if(dabs(dx(l)).gt.cellh) x(i,l)=x(i,l)-dsign(cell,dx(l))
          enddo
        enddo
        call force(fa,cell,cellh,x,sigma,e,u,npart)
        do i=1,npart
          v(i,:)=v(i,:)+dt*fa(i,:)/(2.0d0*wm)
        enddo
        fb=fa
c summation of the velocities
        sv2=0.d0
        do i=1,npart
          sv2=sv2+v(i,1)**2+v(i,2)**2+v(i,3)**2
        enddo
        av2=sv2/npart
        t=wm*av2/(3.0d0*wkb)
c calculation of the total energy
        epot = u * avo / npart
        ekin = av2 * wm / 2.0d0 * avo
        etot = epot + ekin
c monitor
        if(mod(k,1).eq.0)then
          write(10,200) time0,ekin,epot,etot,t
          write(6,*) 'Ek ,Ep, E(',k,')=',ekin,epot,etot
          write(6,*) 'time=',time0,'  T=',t
        end if
c velocity scaling
        call cpu_time(t2)
        if(mod(k,1).eq.0) write(*,*)'ElapsedCPUtime= ',
     $   t2-t1
      enddo
c output
      write(3,400) npart,nstep,sigma,e,wm,cell,tset,dt
      do i=1,npart
        write(3,200) x(i,1),x(i,2),x(i,3)
        write(3,200) v(i,1),v(i,2),v(i,3)
      enddo
      write(3,200) time0
  400 format(2i6,6d23.15)
  200 format(5d23.15)
      end

c subroutine force
      subroutine force(f,cell,cellh,x,sigma,e,u,npart)
      implicit double precision(a-h,o-z)
      real*8 f(100000,3),x(100000,3),dx(3),fdx(3)
      rc=3.0d0*sigma
      r2=sigma*sigma
      rc2=rc*rc

      u = 0.0d0
      f = 0.0d0

      do i=1,npart-1
        do j=i+1,npart
          do l=1,3
            dx(l)=x(i,l)-x(j,l)
            if(dabs(dx(l)).gt.cellh) dx(l)=dx(l)-dsign(cell,dx(l))
          enddo 
          d2=dx(1)*dx(1)+dx(2)*dx(2)+dx(3)*dx(3)
          if(d2.lt.rc2) then
            rd2=r2/d2
            rd6=rd2*rd2*rd2
            rd12=rd6*rd6
            fk=48*e*(rd12-rd6*0.5d0)
            fdx(:)=fk*dx(:)/d2
            f(i,:)=f(i,:)+fdx(:)
            f(j,:)=f(j,:)-fdx(:)
            u=u+4*e*(rd12-rd6)
          end if
        enddo
      enddo
      return
      end