Mix colors from different objects at collision (Rigid Bodies and Particles)
Using Animation Nodes
Updated for Animation Nodes 2.0!
Method with the Molecular Modifier and CubeSurfer below!!!
Another way to do this is using @Jacques Lucke's Animation Nodes addon. Here is the final result:

And here is a picture of the node setup:

This is done with a Rigid Body Simulation and Vertex Colors. After installing Animation Nodes (linked above), we can get to work.
1. Setting Up the Scene
First, we have to set up our Rigid Bodies. Add a ground plane, and under the Physics tab of the Properties Editor, click on Rigid Body and uncheck Active. For every object that we want to actively collide, click on Rigid Body, but leave Active checked.

In addition, under the Mesh tab of the Properties Editor, add a vertex color by clicking the + button under the Vertex Colors panel. Since we have to do this for every object, we can save ourselves some typing by leaving the default name, Col.

The last thing we have is give it a material. Call it whatever you want, and then create the following (Cycles) NodeTree:

This can be modified as needed, as long as the vertex color is still a part of the nodetree.
Repeat the above for all of the objects you want to mix color, giving the the same material, or just duplicate the one object. Finally, select all of the objects you want to mix colors and group them together (Ctrl + G). Don't forget to name it, either in the T panel or by hitting F6.

The scene is set up now. Let's dive into the Animation Nodes!
2. Creating the Nodetree
In the node window, click on the Animation Nodes. Add an Objects From Group node (to add a node, either press Shift + A and find the node, or press Ctrl + A to search for it). Select the Group.

The first thing we want to do is to assign random colors to each object on Frame 1. Let's create a subprogram to do that. Add a Loop Input node and name it SetRandCol. Press the + button next to New Iterator and type &quot;Object List&quot;. Set up the following nodes:

To activate the loop, add an Invoke Subprogram node, select SetRandCol, and plug in the Object list. It should take effect immediately.

We don't want this to happen every frame, though. Just the first. Add these nodes to make it only work on Frame 1:

Cool. The next step is to calculate all of the different pairs of objects we will have to check collision against. We could do this with nodes, but it is faster to just use a short script. First, add a Get List Length node and plug in the Object list. Next, add a Script node and name it FindPairs. Press the + button to add a new script. Use the buttons to add 1 input (an integer) and 2 outputs (integer lists). Name them as follows:

In a Text Editor window, open the FindPairs script and paste the following:

obj1, obj2 = [],[]

for i in range(0,objcount):
    for j in range(i,objcount):
        if i != j:
            obj1.append(i)
            obj2.append(j)

To use the script, add an Invoke Subprogram node and select the &quot;FindPairs&quot; option. Plug the Get List Length node into the script node. obj1 and obj2 are all the pairs of objects (i.e. index 0 goes with index 0, etc).

The next thing to do is to create a loop that checks every pair for collision. We will do this with the Is Inside Volume node. Create a loop the same way we did the first time. Name it BVHCollide. Give it two &quot;Integer List&quot; iterators (named obj1 and obj2), an &quot;Object List&quot; parameter, and a &quot;Float&quot; parameter named Epsilon.

Add these nodes:

Next, we need to check for collision between the two objects. We'll use the Is Inside Volume node. After we check for collision, we'll combine our data to get a single boolean value for whether or not the two objects intersect:

The Epsilon parameter determines how far outside the loop the BVH node will look. Usually a super-low, non-zero value will work, like 0.0002. The last thing to do is to change the color of the objects. Add one more loop:

In the Advanced Node Settings dropdown of the N panel for the Expression nodes, be sure to uncheck &quot;Correct Type&quot;:

Finally, invoke this last loop, and your simulation should work!


Blend File (requires Animation Nodes):


If you want to set the color manually, you can create this loop. Just provide it a list of colors with corresponding indices to the list of objects:


Incorporating the Molecular Modifier
Incorporating the Molecular Modifier has a similar workflow to the Rigid Body Method. Begin by setting up a Molecular simulation. Here is mine:

Here is the final result:

Node setup:

After setting up the molecular simulation, we can work on the coloring. To color the particles, we are exploiting a feature of the Molecular Modifier: since the entire simulation is baked with its own physics, the particles never rotate and therefore never use the Angular Velocity attribute. We will store our colors there. Set all of the particles to render as an Icosphere object and give the Icosphere this material:

Now lets get into Animation Nodes. Unfortunately, in Animation Nodes 2.0, the Particle datatype has been removed. The most efficient way to get around this is to create our own in Python. First, let's create an Extract Particles node. Add a Script node, and give it Particle System input and a Generic List output. Name the List ParticleList. Create a new script file and paste in this code:
ParticleList = []

for p in ParticleSystem.particles:
    ParticleList.append(p)


Next, we need a Filter Particles node. Add another Script node, and give it a Generic List input, four Boolean inputs, and a Generic List output. Name them as follows:

Here is the script:
Filtered = []

for p in Particles:
    if Unborn:
        if p.alive_state == 'UNBORN':
            Filtered.append(p)
    if Alive:
        if p.alive_state == 'ALIVE':
            Filtered.append(p)
    if Dying:
        if p.alive_state == 'DYING':
            Filtered.append(p)
    if Dead:
        if p.alive_state == 'DEAD':
            Filtered.append(p)

Last, we need a Particle Data node, to get information from the particles. Add another Script node, and give it a Generic List input, a Vector List output, a Float List output, and a Color List output. Name them as follows:

Here is the script:
Locations, Sizes, AngularVelocities = [],[],[]

for p in ParticleList:
    Locations.append(p.location)
    Sizes.append(p.size)
    AngularVelocities.append(p.angular_velocity)

For some reason, Animation Nodes doesn't like us using a Color List to store vector data, so disable Correct Output Types in the Advanced Node Settings panel of the N sidebar for the Script Node:

Now that we have that set up, we can get on to coloring the particles. For each of our particle emitters, add a Particle Systems from Object node. Add one of our ExtractParticles nodes for each Particle System:

We need a way to color the particles. Create a new Script subprogram named Color Particles. Give it two Color inputs and one Generic List input:

Code:
for p in particles:
    if color1 != color2:
        p.angular_velocity = [(color1[0]+color2[0])/2,(color1[1]+color2[1])/2,(color1[2]+color2[2])/2]

To set the initial colors for each particles system, use a FilterParticles node to get unborn particles, then the Color Particles node to set the color. Note that the two color inputs have to be different for the node to take effect:

Now, join all the particles together with a Combine List node and use a FilterParticles node to get a list of all Alive particles:

Next, we need to calculate all the different pairs of particles to check collision for. We can just use the same node we created above:

At this point, the numbers start to get really big, and unless you have an über-fast computer, you're going to start to lag. Even then, this NodeTree works better if it's only executed once a frame. If you haven't already, you should turn off Always in the Auto Execution panel of the T sidebar:

Ok. Next, let's check collision. Since everything is a perfect sphere, the math is pretty simple. The Loop node will need 2 Integer List iterators (obj1 and obj2), a Vector List parameter, a Float List parameter, a Float parameter, and a Boolean List output. Here is the NodeTree:

We can invoke it like so:

We now have a problem. We have an enormous amount of data. For my little 8000 particle simulation, there are over 2 million particle pairs that have to be checked for collision. If we also have to check if they need to have their colors modified, it would take forever. So, we should optimize the lists to only have pairs for which collision is true. Create a Script node with an  Integer List and a Boolean List input and an Integer List output:

Here is the script:
n = [i[a] for a in range(len(i)) if b[a]]

Use this Prune node on both of the obj lists:

Lastly, we need to change the color of all the colliding pairs. Since we already have a script to calculate and change the color of a pair of particles, we don't have to do much. Create this Loop:

Make sure to change the Slice Type from Index to Length. To finish, put in the Update Particles node:

Done!!!
Notes:
This will not run in real-time. For an 8000 particle simulation, it takes my (about average) computer 6-7 seconds to calculate each frame. Realize that the number of pairs increases quadratically, not linearly, so a 10000 particle simulation will take roughly four times as long as a 5000 particle simulation, depending on how many collisions there are.
Be sure to save before rendering. With both the Molecular script and Animation Nodes running, my computer froze once while rendering (I only have 8GB of RAM, though).


Blend File (requires Animation Nodes and Molecular):


Using CubeSurfer
We can even mesh our simulation and transfer the colors from the particles to the mesh. I'm using CubeSurfer to mesh my particles. It's pretty easy to pick up (controls are at the bottom of the Object panel) and there are tutorials online. Doing this will slow down the simulation a little more, but the NodeTree is pretty simple:

The two Reroutes are bringing in the list of active particles, and the locations of those particles, as labeled, from their respective nodes in the Main frame of the previous NodeTree. The object in the Object Input node is the Isosurface object created by CubeSurfer. You'll need to enable visibility for the Polygon Centers output of the Object Mesh Data node in the N panel. Here is the SetVertexCol script:
Object.data.vertex_colors.new() #Add a new vertex color, since it gets deleted at every frame change.
v_colors = Object.data.vertex_colors[&quot;Col&quot;].data #Get the vertex color data
colors = [ParticleList[i].angular_velocity for i in Indices] #Get the colors for each face from the particles
i,j = 0,0
for f in Object.data.polygons:        #Iterate through each face
    for l in f.loop_indices:          #Iterate through each vertex
        v_colors[i].color = colors[j] #Paint vertex with color
        i += 1                        #Increment vertex
    j += 1                            #Increment face

The above code works because Vertex Color data is stored as a list of colors sorted by face, then vertex, which happens to be in the same order the Animation Nodes gets the faces.
Final Result:

It looks pretty awesome. Unfortunately, converting it to .gif form took down the quality. You can watch the full clip on YouTube here. I recommend watching in slow motion so you can see it frame by frame.

Blend File (Requires Animation Nodes, Molecular, and CubeSurfer):