anchor.diff:  
           call md_vel(solvent, solvent_cells% edges, dt)
           
  -        colloids% vel = colloids% vel + &
  -             dt * ( colloids% force + colloids% force_old ) / (2 * mass)
  +        do k=1, colloids% Nmax
  +           colloids% vel(:,k) = colloids% vel(:,k) + &
  +             dt * ( colloids% force(:,k) + colloids% force_old(:,k) ) / (2 * mass(colloids%species(k)))
  +        end do
   
           call rattle_dimer_vel(colloids, d, dt)
   
  
positive.code:   use mt19937ar_module
    use mpcd
    use md
    use iso_c_binding
    implicit none
  
  
negative.code:   write(*,*) '    i           |    e co so     |   e co co     |   kin co      |   kin so      |   total       |   temp        |'
    write(*,*) ''
  
    do i = 1, N_loop
       md: do j = 1, N_MD_steps
          call md_pos(solvent, dt)
  
          ! Extra copy for rattle
          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)
  
  