anchor.diff:  
     call solvent_colloid_lj% init(epsilon, sigma, sigma_cut)
   
  -  epsilon(1,1) = PTread_d(config, 'epsilon_C_C')
  -  epsilon(1,2) = PTread_d(config, 'epsilon_N_C')
  -  epsilon(2,1) = PTread_d(config, 'epsilon_N_C')
  -  epsilon(2,2) = PTread_d(config, 'epsilon_N_N')
  +  epsilon(1,1) = PTread_d(config, 'epsilon_C_C', loc=params_group)
  +  epsilon(1,2) = PTread_d(config, 'epsilon_N_C', loc=params_group)
  +  epsilon(2,1) = epsilon(1,2)
  +  epsilon(2,2) = PTread_d(config, 'epsilon_N_N', loc=params_group)
   
     sigma(1,1) = 2*sigma_C
     sigma(1,2) = sigma_C + sigma_N
  
positive.code:   call time_change%init('change')
  
    call h5open_f(error)
  
    prob = PTread_d(config,'probability')
    bulk_rmpcd = PTread_l(config, 'bulk_rmpcd')
    bulk_rate = PTread_d(config, 'bulk_rate')
  
    L = PTread_ivec(config, 'L', 3)
    rho = PTread_i(config, 'rho')
    N = rho *L(1)*L(2)*L(3)
  
    T = PTread_d(config, 'T')
    d = PTread_d(config, 'd')
    
    tau = PTread_d(config, 'tau')
    N_MD_steps = PTread_i(config, 'N_MD')
    dt = tau / N_MD_steps
    N_loop = PTread_i(config, 'N_loop')
  
    sigma_C = PTread_d(config, 'sigma_C')
    sigma_N = PTread_d(config, 'sigma_N')
  
    ! solvent index first, colloid index second, in solvent_colloid_lj
    epsilon(:,1) = PTread_dvec(config, 'epsilon_C', 2)
    epsilon(:,2) = PTread_dvec(config, 'epsilon_N', 2)
  
    sigma(:,1) = sigma_C
    sigma(:,2) = sigma_N
  
negative.code:   solvent% force_old = solvent% force
    colloids% force_old = colloids% force
  
    call h5fflush_f(hfile%id, H5F_SCOPE_GLOBAL_F, error)
    write(*,*) 'Running for', N_loop, 'loops'
    write(*,*) 'mass', mass 
    call main%tic()
  