anchor.diff:  
     call colloid_lj% init(epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:))
   
  -  epsilon(1,1) = 1.d0 
  -  epsilon(1,2) = 1.d0 
  -  epsilon(2,1) = 1.d0 
  -  epsilon(2,2) = 1.d0 
  +  epsilon(1,1) = 1.d0
  +  epsilon(1,2) = 1.d0
  +  epsilon(2,1) = 1.d0
  +  epsilon(2,2) = 1.d0
     sigma(1,:) = [sigma_C, sigma_N]
     sigma(2,:) = [sigma_C, sigma_N]
     sigma_cut = sigma*2**(1.d0/6.d0)
  
positive.code: 
    call solvent_colloid_lj% init(epsilon, sigma, sigma_cut)
  
    epsilon(1,1) = 1.d0 !PTread_d(config, 'epsilon_C_C')
    epsilon(1,2) = 1.d0 !PTread_d(config, 'epsilon_N_C')
    epsilon(2,1) = 1.d0 !PTread_d(config, 'epsilon_N_C')
    epsilon(2,2) = 1.d0 !PTread_d(config, 'epsilon_N_N')
  
    sigma(1,1) = 2*sigma_C
    sigma(1,2) = sigma_C + sigma_N
  
negative.code:   double precision :: conc_z(400)
    double precision :: colloid_pos(3,2)
  
    type(PTo) :: config
    integer(c_int64_t) :: seed
    integer :: i, L(3),  n_threads
    integer :: j, k, m
  