anchor.diff:  
     e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
     e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
  -  e3 = compute_harmonic()
  +  if (do_harmonic) then
  +     e3 = compute_harmonic()
  +  else
  +     e3 = 0
  +  end if
  +
     if (do_zwall) then
        e_wall = lj93_zwall(colloids, solvent_cells% edges, walls_colloid_lj)
     else
  
positive.code:   mass(1) = rho * sigma(1,1)**3 * 4 * 3.14159265/3
    mass(2) = rho * sigma(1,2)**3 * 4 * 3.14159265/3
  
    call solvent% init(N,2, system_name='solvent') !there will be 2 species of solvent particles
  
    call colloids% init(2,2, mass, system_name='colloids') !there will be 2 species of colloids
  
    call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)
  
  
negative.code: 
      double precision :: f(3), x12(3), dist
  
      x12 = rel_pos(colloids%pos(:,1), colloids%pos(:,2), solvent_cells%edges)
      dist = norm2(x12)
  
  