anchor.diff:    end do
     spare2 = 0.0d0
     do i = 1, N_latt_moves_eq
  -    call random_number(spare)
  +    spare = univar_3Nplus10001(3*N + i)
       if (spare > 0.5d0) then ! stretch move
  +      u_rnd6 = univar_60k(6*(i-1)+1:6*(i-1)+6)
         call roam_latticeshape_stretch(h0, h0_temp, s, 1.0d0, 0.35d0, min_height, &
  -                             & 1.0d0, 1.d0, spare_int, spare, spare2, skip_enth_calc=.true. )
  +                             & 1.d0, u_rnd6, spare_int, spare, spare2, skip_ener_calc=.true. )
         h0=h0_temp
       else ! shear move
  +      u_rnd3 = univar_30k(2*(i-1)+1:2*(i-1)+3)
  +      g_rnd6 = nvar_60k(6*(i-1)+1:6*(i-1)+6)
         call roam_latticeshape_shear(h0, h0_temp, s, 1.0d0, 0.5d0, min_height, &
  -                             & 1.0d0, 1.d0, spare_int, spare, spare2, skip_enth_calc=.true. )
  +                             & 1.d0, u_rnd3, g_rnd6, spare_int, spare, spare2, &
  +                             & skip_ener_calc=.true. )
         h0=h0_temp
       end if
     end do  
   
  -  ! initialise volume with prior V^N over (0,Vmax)
  -  call polynomial_rand(real(N,dp),Vmax,single=V)
  -
  -  enth = enthalpy( s, h0, V, pressure, cutoff )
  +  if (flat_v_prior) then
  +    V = univar_3Nplus10001(3*N + 10001)
  +    V = V*Vmax
  +  else
  +    V = univar_3Nplus10001(3*N + 10001)
  +    V = V**(1.0d0/real(N+1,dp))
  +    V = V*Vmax
  +  end if
  +
  +  ener = energy( s, h0, V, cutoff )
   
   end subroutine initialise_config
   
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  -!! ROUTINES FOR CALCULATING THE ENTHALPY OF A LENNARD JONES SYSTEM !!!!!!!!!!!!!!!!!!!!!!!!!!
  +!! ROUTINES FOR CALCULATING THE ENERGY OF A LENNARD JONES SYSTEM !!!!!!!!!!!!!!!!!!!!!!!!!!
   !! EPSILON = 1.0 , SIGMA = 1.0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !! STANDARD MEAN FIELD CORRECTION FOR r > cutoff !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
     !***********************************************************************************
  -  ! function to calculate the enthalpy PE + PV + longrange (mean field) correction
  +  ! function to calculate the energy PE + longrange (mean field) correction
     ! if atom is specified, PE only contains the contribution due to the interactions of atom "atom"
     !***********************************************************************************
  -  function enthalpy( s, h0, V, press, cutoff, atom, n_call_enthalpy )
  -
  -    real(dp) :: enthalpy
  +  function energy( s, h0, V, cutoff, atom, n_call_energy )
  +
  +    real(dp) :: energy
   
       real(dp), dimension(:,:), intent(in) :: s ! fractional coordinates
       real(dp), dimension(3,3), intent(in) :: h0 ! normalised cell matrix
       real(dp), intent(in) :: V ! volume
  -    real(dp), intent(in) :: press ! units of eps/(sigma**3)
       real(dp), intent(in) :: cutoff ! units of sigma
       integer, optional, intent(in)    :: atom
  -    real(dp), optional, intent(inout)    :: n_call_enthalpy
  +    real(dp), optional, intent(inout)    :: n_call_energy
   
       real(dp), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_dp
       real(dp) :: mfc
   
       if ( present(atom) )  then
  -       enthalpy = lj_energy_correct_longrange( &
  +       energy = lj_energy_correct_longrange( &
            & s,h0,V,cutoff,iat=atom)
       else ! atom not specified
  -       enthalpy = lj_energy_correct_longrange( &
  +       energy = lj_energy_correct_longrange( &
            & s,h0,V,cutoff)
       end if
   
       ! add mean field correction beyond cutoff
       mfc = (2.0_dp/3.0_dp)*pi*4.0d0*dble((size(s)/3)**2)
       mfc = mfc*((((1.0d0/cutoff)**9)/3.0_dp)-((1.0d0/cutoff)**3))/V
  -    enthalpy = enthalpy + mfc
  -
  -    enthalpy = enthalpy + ( press * V )
  -
  -    if(present(n_call_enthalpy)) n_call_enthalpy = n_call_enthalpy + 1.0
  -
  -  end function enthalpy
  +    energy = energy + mfc
  +
  +    if(present(n_call_energy)) n_call_energy = n_call_energy + 1.0
  +
  +  end function energy
   
     function lj_energy_correct_longrange(s,h0,V,cutoff,iat)
       ! Calculate explicit interaction part of LJ potential energy 
  
positive.code: ! 1. MC MOVES
  ! 2. ROUTINE FOR CREATING AN INITIAL CONFIGURATION
  ! It should not be neccessary to read the last two sections
  ! 3. ROUTINES FOR CALCULATING THE ENTHALPY OF A LENNARD JONES SYSTEM
  ! 4. VARIOUS NUMERICAL FUNCTIONS AND ROUTINES REQUIRED BY THE ABOVE
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
negative.code:   !   ************************************************************
    !   * Single atom MC move                                      *
    !   ************************************************************
    subroutine roam_oneatom_onestep(s_in, s_out, h0, V, delta_ener, step_1a1s, iat, cutoff, &
                                  & univar_rnd4, n_call_ener, skip_ener_calc)
  
      real(dp), dimension(:,:), intent(in)   :: s_in  ! fractional coords in
      real(dp), dimension(:,:), intent(out)  :: s_out ! fractional coords out
      real(dp), dimension(3,3), intent(in)   :: h0    ! normalised cell matrix
      real(dp), intent(in)                   :: V     ! volume
      real(dp), intent(out)                  :: delta_ener ! energy change due to single atom step
      real(dp), intent(in)      :: step_1a1s          ! step length
                                                      ! defined as a fraction of unit_cell_volume**(1/3)
      integer, intent(out)      :: iat                ! index of random atom that was chosen moved 
                                                      ! element of (1,2,...,N)
      real(dp), intent(in)      :: cutoff             ! radial cutoff of potential, defined in units of sigma
      real(dp), dimension(4), intent(in) :: univar_rnd4      ! an array of 4 iid univariate random numbers 
                                                      ! on (0,1)
      real(dp), intent(inout)   :: n_call_ener        ! number of entahlpy calls (tally) - updated by +1 
      logical, optional, intent(in) :: skip_ener_calc ! OPTIONAL. if .true. do not calculate delta_energy 
      ! internal arguments
      real(dp), dimension(3) :: translation_roos
      real(dp), dimension(4) :: univar_rnd4mod
      integer :: i, N
  
      N = size(s_in)/3
  
      iat = int(univar_rnd4(1)*real(N)) + 1
      ! choose random particle between 1 and N for taking the step
  
      univar_rnd4mod = univar_rnd4 * 2.0_dp * step_1a1s - step_1a1s  
  
      s_out = s_in
      do i = 2,4
         ! do move in fractional coordinates
         s_out(i-1,iat) = s_out(i-1,iat) + univar_rnd4mod(i)
      end do
  
      ! map back into cell if atom has stepped out of periodic boundaries
      translation_roos = nint(s_out(:,iat) - 0.5_dp, dp)
      s_out(:,iat) = s_out(:,iat) - translation_roos
  
      if (.not. (present(skip_ener_calc).and.skip_ener_calc)) then 
        delta_ener = energy( s_out, h0, V, cutoff, atom=iat, n_call_energy=n_call_ener )
        delta_ener = delta_ener - energy( s_in, h0, V, cutoff, atom=iat ) 
        ! only count one energy call for this move
      end if
  
    end subroutine roam_oneatom_onestep
  