anchor.diff:  real(c_double), parameter :: mass = 1.0
   real(c_double), parameter :: C6 = 4.0
   real(c_double), parameter :: C12 = 4.0
  +real(c_double), parameter :: cutoff = 3.0
  +real(c_double), parameter :: cutoffsq = cutoff*cutoff
   
   integer :: timesteps_read = 0
   
  +!interparticle distance
  +real(c_double), dimension(3) :: del
  +real(c_double) :: rsq
  +
   !intermediate variables
  -real(c_double) :: tcon_numerator = 0.0
  -real(c_double) :: tcon_denominator = 0.0
  -real(c_double) :: tconF = 0.0
  +real(c_double) :: tconF_numerator = 0.0
  +real(c_double) :: tconF_denominator = 0.0
   real(c_double) :: tcon1 = 0.0
   real(c_double) :: tkin = 0.0
  +real(c_double) :: tconF = 0.0
   
   !accumulators for averages over entire simulation
   real(c_double) :: tkin_ave = 0.0
  
positive.code: real(c_double) :: tconF_denominator_sq_ave = 0.0
  real(c_double) :: tcon1sq_ave = 0.0
  
  call get_command_argument(1, filename, filename_len)
  
  file_opened = OpenFile(filename, filename_len)
  
negative.code:         print*, "Couldn't open file ", filename
  else
          open(unit=output_file, file=output_filename, action="write", status="replace")
          write(output_file,*) "# Timestep   Tkin   Tcon"
          do
                  loop_continue = readframe(timestep, N, box_lo, box_hi, x, v, f)
                  if (.not. loop_continue) exit
                  print*, timestep
  
                  ! calculate the periodic replica distance
                  ! needed for the minimum_image subroutune
                  prd = box_hi - box_lo
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
                  !!! START MODIFYING HERE !!!
  
                  ! the following have been read from the trajectory for you:
                  ! integer timestep
                  ! integer N (number of atoms)
                  ! real box_lo(3) minimum coordinates in all three directions
                  ! real box_hi(3) maximum coordinates in all three directions
                  ! real prd(3) periodic replica distance in all three directions
                  ! real x(N,3) atomic positons
                  ! real v(N,3) atomic velocities
                  ! real f(N,3) total force on each atom F_i
  
                  ! the following subprogrammes are provided:
                  ! - minimum_image(vector, prd) : apply the minimum image convention
                  ! - veclengthsq(vector) : return the squared length of the vector x**2 + y**2 + z**2
  
                  ! You need to give values to the following variables:
                  ! - tkin
                  ! - tconF_denominator
                  ! - tconF_numerator
                  ! the existing code will take care of the averages and standard deviations
  
                  ! two loop counters have been predefined for you: i, j
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
                  !!! STOP MODIFYING HERE !!!
                  ! (Unless you wish to add a function/subroutine) !
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
                  tkin = tkin / N
                  tkin_ave = tkin_ave + tkin
                  tkinsq_ave = tkinsq_ave + tkin*tkin
  
                  tcon1 = tconF_denominator / tconF_numerator
                  tcon1_ave = tcon1_ave + tcon1
                  tcon1sq_ave = tcon1sq_ave + tcon1*tcon1
  
                  tconF_numerator_ave = tconF_numerator_ave + tconF_numerator
                  tconF_numerator_sq_ave = tconF_numerator_sq_ave + tconF_numerator*tconF_numerator
                  tconF_denominator_ave = tconF_denominator_ave + tconF_denominator
                  tconF_denominator_sq_ave = tconF_denominator_sq_ave + tconF_denominator*tconF_denominator
                  !write the current temperatures to the output file
                  write(output_file,*) timestep, tkin, tcon1
  
                  timesteps_read = timesteps_read + 1
  
                  !reset the instantaneous variables for next time
                  tconF_numerator = 0.0
                  tconF_denominator = 0.0
                  tcon1 = 0.0
                  tkin = 0.0
                  tconF = 0.0
          end do
          close(output_file)
  
          tkin = tkin_ave / timesteps_read
          tkin_err = stderr(tkin, tkinsq_ave / timesteps_read, timesteps_read)
          tcon1 = tcon1_ave / timesteps_read
          tcon1_err = stderr(tcon1, tcon1sq_ave / timesteps_read, timesteps_read)
          tconF_numerator = tconF_numerator_ave / timesteps_read
          tconF_numerator_err = stderr(tconF_numerator, tconF_numerator_sq_ave / timesteps_read, timesteps_read)
          tconF_denominator = tconF_denominator_ave / timesteps_read
          tconF_denominator_err = stderr(tconF_denominator, tconF_denominator_sq_ave / timesteps_read, timesteps_read)
  
          tconF = tconF_denominator / tconF_numerator
          tconF_err = tconF**2 * ((tconF_numerator_err/tconF_numerator)**2 + (tconF_denominator_err/tconF_denominator)**2)
  
          print*, "Read", timesteps_read, "timesteps."
          print*, "Temperatures"
          print*, "Kinetic: ", tkin, "+/-", tkin_err 
          print*, "Con1", tcon1, "+/-", tcon1_err
          print*, "ConF", tconF, "+/-", tconF_err
  end if
  contains
  real(c_double) function stderr(X, X2, N)
          implicit none
          real(c_double), intent(in) :: X, X2
          integer(c_int), intent(in) :: N
  
          stderr = X2 - X**2
          stderr = sqrt(stderr/N)
  end function stderr
  real(c_double) function veclengthsq(vec)
          implicit none
          real(c_double), dimension(3), intent(in) :: vec
  
          veclengthsq = vec(1)*vec(1) + vec(2)*vec(2) + vec(3)*vec(3)
  end function veclengthsq
  subroutine minimum_image(vec, prd)
          real(c_double), dimension(3), intent(inout) :: vec
          real(c_double), dimension(3), intent(in) :: prd
          integer i
  
          do i=1,3
                  if (vec(i) < -prd(i)/2.0) then
                          vec(i) = vec(i) + prd(i)
                  else if (vec(i) > prd(i)/2.0) then
                          vec(i) = vec(i) - prd(i)
                  end if
          end do
  end subroutine minimum_image
  end program analyse
  