anchor.diff:  
     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')
  +  epsilon = 1.d0
   
     sigma(1,1) = 2*sigma_C
     sigma(1,2) = sigma_C + sigma_N
  
positive.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
  
negative.code: 
    call colloids% init(2,2, mass) !there will be 2 species of colloids
  
    call PTkill(config)
  
    open(17,file ='dimerdata_FullExp_1.txt')
    open(18,file ='dimerdata_FullExp_2.txt')
    open(19,file ='dimerdata_vx_flow_wall.txt')  
  