subroutine SetTemperature(v,N,m,T)
  USE GLOBAL
  implicit none

  !
  ! Subroutine sets the temperature of the atoms in the system
  ! to T assuming a Maxwell-Boltzmann velocity distribution
  !
  type(vector) :: v(*)
  integer :: N
  real(double) :: T,m

  ! stuff removed

end subroutine SetTemperature


subroutine GetTemperature(v,N,m,T)
  USE global
  implicit none

  !
  ! Subroutine gets the temperature of the atoms in the system
  !
  type(vector) :: v(*)
  integer :: N
  real(double) :: T,m

  real(double) :: Ekinsum,gaussianrand
  integer :: i

  ! Get sum of kinetic energies
  Ekinsum=0
  do i=1,N

     Ekinsum=Ekinsum+v(i)%x*v(i)%x+v(i)%y*v(i)%y+v(i)%z*v(i)%z

  enddo
  ! Do necessary unit transformations to get E in J
  Ekinsum=0.5*m*u*Ekinsum*vunit*vunit
  
  ! and get temperature using E=3/2 N k T
  T = Ekinsum/(1.5*N*kB)
  
  return

end subroutine GetTemperature

subroutine GetEnergies(v,N,m,Ekinsum,Epotsum,Etot,Ekin,Epot)
  USE GLOBAL
  implicit none

  type(vector) :: v(*)
  real(double) :: Ekinsum,Epotsum,Etot,m,Ekin(*),Epot(*)
  integer :: N
  
  real(double) :: help1
  integer :: i
  
  Ekinsum=0.0; Epotsum=0.0;
  help1=0.5*m*u*vunit*vunit/e;
  do i=1,N
     Ekin(i)=help1*(v(i)%x*v(i)%x+v(i)%y*v(i)%y+v(i)%z*v(i)%z)
     Ekinsum=Ekinsum+Ekin(i)
     Epotsum=Epotsum+Epot(i)
  enddo
  Etot=Ekinsum+Epotsum

  return

end subroutine GetEnergies
  

!---------------------------------------------------------------------------

double precision function uniformrand(seed)
  implicit none

  ! -----------------------------------------------------------
  ! Park-Miller "minimal" Random number generator
  ! uniform distribution ]0,1[ . See Numerical Recipes ch. 7.0
  ! -----------------------------------------------------------

  integer :: seed,IA,IM,IQ,IR,MASK
  double precision :: AM
  integer :: k
  
  parameter (IA=16807,IM=2147483647,AM=1.0/IM,IQ=127773,IR=2836,MASK=123459876)

  seed=ieor(seed,MASK)
  k=seed/IQ
  seed=IA*(seed-k*IQ)-IR*k
  if (seed < 0) seed=seed+IM
  uniformrand=AM*seed
  seed=ieor(seed,MASK)
  
  return
  
end function uniformrand



double precision function gaussianrand(seed)
  implicit none

  ! ---------------------------------------------------------
  ! Random numbers with normal (Gaussian) distribution.
  ! Mean is 0 and standard deviation is 1
  ! ---------------------------------------------------------

  integer :: seed

  ! Stuff removed
  
end function gaussianrand