anchor.diff:          ! Extra copy for rattle
           if (do_rattle) &
                colloids% pos_rattle = colloids% pos
  -        if (do_harmonic) &
  -             e3 = compute_harmonic()
           do k=1, colloids% Nmax
              colloids% pos(:,k) = colloids% pos(:,k) &
                   + dt * colloids% vel(:,k) + &
                   dt**2 * colloids% force(:,k) / (2 * colloids% mass(k))
           end do
  +
           if (do_rattle) &
                call rattle_dimer_pos(colloids, d, dt, solvent_cells% edges)
   
  
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)
  
  