anchor.diff:          
           pos_rattle = colloids% pos
           colloids% pos_old = colloids% pos + dt * colloids% vel + dt**2 * colloids% force / (2 * mass)
  -        
  +
           call rattle_pos
           
           do k = 1, colloids% Nmax
  
positive.code:         
          call rattle_vel
          
          call flag_particles
          
          call change_species
  
       end do md
  
       write(15,*) colloids% pos + colloids% image * spread(solvent_cells% edges, dim=2, ncopies=colloids% Nmax), &
                   colloids% vel, e1+e2+mass*sum(colloids% vel**2)/2+sum(solvent% vel**2)/2
  
       call solvent% sort(solvent_cells)
       call neigh% update_list(colloids, solvent, 1.5d0*sigma, solvent_cells)
  
       call simple_mpcd_step(solvent, solvent_cells, mt)
  
       write(*,'(7f16.3)') real(i),e1, e2, mass*sum(colloids% vel**2)/2, sum(solvent% vel**2)/2, &
           e1+e2+mass*sum(colloids% vel**2)/2+sum(solvent% vel**2)/2, &
           compute_temperature(solvent, solvent_cells)
  
    end do
    
  
negative.code:   colloids% force_old = colloids% force
  
    write(*,*) ''
    write(*,*) '    e co so     |   e co co     |   kin co      |   kin so      |   total       |   temp        |'
    write(*,*) ''
   
    write(*,'(6f16.3)')e1, e2, mass*sum(colloids% vel**2)/2, sum(solvent% vel**2)/2, &
           e1+e2+mass*sum(colloids% vel**2)/2+sum(solvent% vel**2)/2, &
           compute_temperature(solvent, solvent_cells)
   
    do i = 1, 20
       so_max = 0
       co_max = 0
       md: do j = 1, N_MD_steps
  