anchor.diff:  
   open(11, file='energies_temp.dat', status='unknown')
   write(11, *) '# step   ek   eu    etot   temp'
  +open(12, file='energies_temp.dat', status='unknown')
  +write(12, *) '# step   ek   eu    etot   temp'
   
   !!!!!!!!!!!!!!!!!!  Termalization !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   do k = 1, 1000
  
positive.code: end do
  
  !!!!!!!!!!!!!!!!!  Integration of motion   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  do k =1, nstep
      t = t0 + dt * k
      call integrate(f, eu, ek, part, vel, npart, a, r_cut2, e_cut, dt, e_tot, temp_k)
      part = part - l * nint(part/l)
  
negative.code: !!!!!!!!!!!!!!!!!!!!  Formats for printing !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  34 format(I12, 2X, 3(ES14.6e2, 2X))
  36 format(F12.6, 2X, 4(ES14.6e2, 2X))
  !38 format(F12.6, 2X, 6(ES14.6e2, 2X))
  
  
  !!!!!!!!!!!!!!!!!!!!  Initial conditions and settings!!!!!!!!!!!!!!!!!!!!!!!!!!
  