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:     colloids%force(:,1) = colloids%force(:,1) - f
      colloids%force(:,2) = colloids%force(:,2) + f
  
      ! U = k/2 * (dist - d)**2
      ! f = k * (dist-d) * d dist / d x
  
  
negative.code: 
      double precision :: f(3), x12(3), dist
  
      x12 = rel_pos(colloids%pos(:,1), colloids%pos(:,2), solvent_cells%edges)
      dist = norm2(x12)
  
  