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