anchor.diff:  #endif
         write(coll_scatterUnit,'(1x,i2,2x,i4,2x,i5,2x,i1,3(2x,e18.10))') icoll,iturn,ipart,nabs_tmp,(p-pBef)/pBef,xp-xpBef,zp-zpBef
   #ifdef HDF5
  -      endif
  +      end if
   #endif
       end if
       return
     end if
   
  -!++  Otherwise do multi-coulomb scattering.
  -!++  REGULAR STEP IN ITERATION LOOP
  +  ! Otherwise do multi-coulomb scattering.
  +  ! REGULAR STEP IN ITERATION LOOP
     call k2coll_mcs(s)
   
  -!++  Check if particle is outside of collimator (X.LT.0) after
  -!++  MCS. If yes, calculate output longitudinal position (s),
  -!++  reduce momentum (output as dpop) and return.
  -!++  PARTICLE LEFT COLLIMATOR BEFORE ITS END.
  -  if(x.le.zero) then
  -    s=(zlm-rlen)+s
  -
  -#ifndef MERLINSCATTER
  -    call k2coll_calcIonLoss(mat,p,rlen,m_dpodx)
  -    p=p-m_dpodx*s
  -#endif
  +  ! Check if particle is outside of collimator (X.LT.0) after
  +  ! MCS. If yes, calculate output longitudinal position (s),
  +  ! reduce momentum (output as dpop) and return.
  +  ! PARTICLE LEFT COLLIMATOR BEFORE ITS END.
  +  if(x <= zero) then
  +    s = (zlm-rlen)+s
  +
   #ifdef MERLINSCATTER
       call merlinscatter_calc_ion_loss(p,edens(mat),pleng(mat),exenergy(mat),s,m_dpodx)
  -    p=p-m_dpodx
  -#endif
  -    dpop=(p-p0)/p0
  +    p = p-m_dpodx
  +#else
  +    call k2coll_calcIonLoss(mat,p,rlen,m_dpodx)
  +    p = p-m_dpodx*s
  +#endif
  +    dpop = (p-p0)/p0
   
       if(dowrite_impact) then
         ! write coll_scatter.dat for complete scattering histories
  
positive.code: 
      ! Now add Coulomb
      csect(0,ma) = csect(0,ma) + csect(5,ma)
    
      ! Interaction length in meter
      xintl(ma) = (c1m2*anuc(ma))/(((fnavo * rho(ma))*csect(0,ma))*1e-24_fPrec)
  
  
negative.code:     return
    end if
  
    ! Initialize the interaction length to input interaction length
    rlen = zlm
  
    ! Do a step for a point-like interaction.
    ! Get monte-carlo interaction length.
  10 continue
    zlm1     = (-one*xintl(mat))*log_mb(real(rndm4(),fPrec))
    nabs_tmp = 0  ! type of interaction reset before following scattering process
    xpBef    = xp ! save angles and momentum before scattering
    zpBef    = zp
    pBef     = p
  
    ! If the monte-carlo interaction length is longer than the
    ! remaining collimator length, then put it to the remaining
    ! length, do multiple coulomb scattering and return.
    ! LAST STEP IN ITERATION LOOP
    if(zlm1 > rlen) then
      zlm1 = rlen
      call k2coll_mcs(s)
      s = (zlm-rlen)+s
  #ifdef MERLINSCATTER
      call merlinscatter_calc_ion_loss(p,edens(mat), pleng(mat),exenergy(mat),s,m_dpodx)
      p = p-m_dpodx
  #else
      call k2coll_calcIonLoss(mat,p,rlen,m_dpodx)  ! DM routine to include tail
      p = p-m_dpodx*s
  #endif
  
      dpop = (p-p0)/p0
      if(dowrite_impact) then
        ! write coll_scatter.dat for complete scattering histories
  #ifdef HDF5
  