anchor.diff:  
        call simple_mpcd_step(solvent, solvent_cells, mt)
   
  -     kin_co = (mass(1)*sum(colloids% vel(:,1)**2)+mass(2)*sum(colloids% vel(:,2)**2))/2
  +     kin_co = (colloids% mass(1)*sum(colloids% vel(:,1)**2)+ colloids% mass(2)*sum(colloids% vel(:,2)**2))/2
        write(*,'(1i16,6f16.3,1e16.8)') i, e1, e2, &
             kin_co, sum(solvent% vel**2)/2, &
             e1+e2+kin_co+sum(solvent% vel**2)/2, &
  
positive.code:         colloids% pos_rattle = colloids% pos
          do k=1, colloids% Nmax
             colloids% pos(:,k) = colloids% pos(:,k) + dt * colloids% vel(:,k) + &
                  dt**2 * colloids% force(:,k) / (2 * mass(colloids%species(k)))
          end do
  
          call rattle_dimer_pos(colloids, d, dt, solvent_cells% edges)
  
negative.code:   sigma(2,2) = 2*sigma_N
    sigma_cut = sigma*2**(1.d0/6.d0)
  
    call colloid_lj% init(epsilon, sigma, sigma_cut)
  
    mass(1) = rho * sigma_C**3 * 4 * 3.14159265/3
  