anchor.diff:    ! init Janus particle
     colloids%species = 1
     colloids%species(1:n_colloids/2) = 1
  -  colloids%species(n_colloids/2+1:n_colloids) = 1
  +  colloids%species(n_colloids/2+1:n_colloids) = 2
     colloids%vel = 0
     colloids%force = 0
     colloids%pos = reshape([13.3745, 8.98251, 4.69069, 8.5896, 7.72357, 4.80796, 9.47026, &
  
positive.code:           end do
         end if
      end do
  
    end subroutine flag_particles_nl
  
    subroutine change_species
      double precision :: dist_to_C_sq
  
negative.code:                 x = rel_pos(colloids% pos(:,i),solvent% pos(:,r),solvent_cells% edges)
                  dist_to_C_sq = dot_product(x, x)
                  if (dist_to_C_sq < solvent_colloid_lj%cut_sq(1,1)) then
                     if (threefry_double(state(thread_id)) <= prob) then
                        solvent% flag(r) = 1
                     end if
                  end if
  