anchor.diff:  
     double precision :: e1, e2
     double precision :: tau, dt , T
  -  double precision :: d
  +  double precision :: d,prob
     double precision :: skin, co_max, so_max
     integer :: N_MD_steps, N_loop
     integer :: n_extra_sorting
  
positive.code:   call init_genrand(mt, int(clock, c_long))
  
    call h5open_f(error)
  
    L = PTread_ivec(config, 'L', 3)
    rho = PTread_i(config, 'rho')
  
negative.code:     double precision :: x(3)
      integer :: n
  
      far = (L(1)/2)**2
  
  
      do n = 1,solvent% Nmax
         if (solvent% species(n) == 2) then
  