Module mdmod
   implicit none

    type atomtype
       integer :: many,Z
       real(8) :: mass
       real(8), pointer :: r(:,:),v(:,:),F(:,:)
    end type

       real(8)  :: latticevec(3,3)

    integer :: totalparticles, numtypes, timesteps
    real(8) :: temperature, dt
    type (atomtype), pointer :: At(:)

contains

  function ranx1()
     real(8):: ranx1
     integer :: m=100001,konst=125

     m=m*konst
     m=m-2796203*(m/2796203)
     ranx1=m/2796203.d0
     
  end function ranx1

  subroutine initialyze
     integer :: i,j,k,l,m,n,many
     real(8) :: x
     real(8), parameter ::   Tcon=9118.37304d0


     write(6,*) 'Input number of types of particles'
     read(5,*) numtypes

     allocate(At(numtypes))
 
     totalparticles=0
       do i=1, numtypes
          write(6,*) 'For type ', i, &
            '  enter Z,  mass (amu),  No. particles of this type'
          read(5,*) At(i)%Z,At(i)%mass,At(i)%many
          many=At(i)%many
          allocate(At(i)%r(3,many),At(i)%v(3,many),At(i)%F(3,many))
          totalparticles=totalparticles+At(i)%many
       enddo

     write(6,*) 'Pgm initialyzed with total particle number = ', totalparticles

     write(6,*) 'Enter the 3 lattice vectors for this simulations in Angsroms'
     read(5,*) latticevec(:,1)
     read(5,*) latticevec(:,2)
     read(5,*) latticevec(:,3)
     write(6,*) 'Enter initial temperature for simulation'
     read(5,*) temperature
     write(6,*) 'Pgm initialyzed with temperature = ', temperature

     ! set initial positions and velocities
      do i=1,numtypes
         !  prepare scale factor for velocity
         x=Tcon*sqrt(temperature/At(i)%mass)
         do j=1,At(i)%many
            At(i)%r(:,j)=0
            do k=1,3
               At(i)%r(:,j)=At(i)%r(:,j)+latticevec(:,k)*(ranx1()-0.5d0)
               At(i)%v(k,j)=x*(ranx1()-0.5d0)
            enddo
         enddo
       enddo

     write(6,*) 'Input the total number of timesteps and dt (in ps)'
     read(5,*) timesteps,dt
      
  end subroutine initialyze

  subroutine writecoords(icount,iunit)
     integer, intent(IN) :: icount,iunit

     integer:: i,j

     write(iunit,*) 'PRIMCOORD', icount
     write(iunit,*)  totalparticles, '  1'

     do i=1,numtypes     
        do j=1,At(i)%many
           write(iunit,'(i10, 3f15.6)') At(i)%Z,At(i)%r(:,j)
        enddo
     enddo
  end subroutine

  subroutine forcefunc(r,FF)
     real(8), intent(IN) :: r(3)
     real(8), intent(OUT) :: FF(3)
     real(8), parameter :: eps=1.d0, sig=2.0

     real(8) :: x,rr
  
     rr=sqrt(r(1)**2+r(2)**2+r(3)**2)
     x=sig/rr
     FF=0
     if (x < 0.02d0) return
     x=-6*eps*(2*(x**13)-(x**7))*x
     FF=r*x
   end  subroutine forcefunc

   subroutine forceonparticle(it,ip,F)
      integer, intent(IN) :: ip,it
      real(8), intent(OUT) :: F(3)

      integer :: i,j,k,n1,n2,n3
      real(8) :: F1(3),r(3)

      F=0

      do i=1,numtypes
         do j=1,At(i)%many
               do n1=-3,3
                  do n2=-3,3
                     do n3=-3,3
                        if (.not.(i==it.and.j==ip.and.n1==0.and.n2==0.and.n3==0)) then
                        r=n1*latticevec(:,1)+n2*latticevec(:,2) &
                           +  n3*latticevec(:,3)+At(i)%r(:,j)-At(it)%r(:,ip)
                        call forcefunc(r,F1)
                        F=F+F1
                        endif
                      enddo
                   enddo
                enddo
           enddo
       enddo        

  end subroutine forceonparticle

   subroutine moveparticles
      integer :: i,j,k
      real(8) :: F(3)

      do i=1,numtypes
         do j=1,At(i)%many
            call forceonparticle(i,j,F)
          ! update positions
            At(i)%r(:,j)=At(i)%r(:,j)+At(i)%v(:,j)*dt+(F*(dt**2)/2)/At(i)%mass
          ! update velocities
            At(i)%v(:,j)=At(i)%v(:,j)+(F*dt/2)/At(i)%mass
         enddo
      enddo
   end subroutine moveparticles
      


end module