Shader: Rounded edge between the two intersecting faces? Is there way to do "Rounded Edge Width" in Blender? It's familiar for me from Modo. Anyone know the trick? 

Like Modo's Rounded Edge Width. Couple links for reference:
From the modo documentation:


  The 'Rounded Edge Width' option adjusts the rendered result of surface
  normals at a polygons edge to blend them with adjoining polygon
  normals, giving the impression of a small rounded edge between the two
  intersecting faces. To explain it more simply it shades polygons edges
  at render time to look as if a small rounded bevel has been applied,
  to what normally would be a sharp intersection. Users can define an
  Edge Width to determine the amount or 'roundness' but keep in mind it
  is merely a shading trick and wont round the edges of the actual
  geometry, nor change the objects silhouette in the rendered image, so
  it is best to keep set the length to just a few pixels in width in the
  final rendered image.
  
  emphasis mine


Video demonstrating the effect.
There's an OSL bevel Shader for Cycles that I believe was originally written by the person who built Cycles.  Several people have been making adjustments to it and have been posting the results to the following thread on BlenderArtists.org  http://www.blenderartists.org/forum/showthread.php?329295-Bevel-shader
OSL Code
Paste the following code into a new Text Datablock in Blender, and use it in the script node of an OSL node tree:
void rng_seed(output int rng, int seed)
{
  int chash = seed;
  if (chash == 0) chash = 1;
  rng = chash * 30391861;
}

float rng_uniform(output int rng)
{
  float res = rng / float(2137483647) * 0.5 + 0.5;
  rng *= 30391861;
  return res;
}

void to_unit_disk(float x, float y, output float x_out, output float y_out)
{
  float r, phi;
  float a = 2.0 * x - 1.0;
  float b = 2.0 * y - 1.0;
    
  if(a &gt; -b) 
  { if(a &gt; b) 
    { r = a;
      phi = M_PI_4 *(b/a);
    }
    else 
    { r = b;
      phi = M_PI_4 *(2.0 - a/b);
  } }
  else 
  { if(a &lt; b) 
    { r = -a;
      phi = M_PI_4 *(4.0 + b/a);
    }
    else 
    { r = -b;
      if(b != 0.0) phi = M_PI_4 *(6.0 - a/b);
      else phi = 0.0;
  } }
  x_out = r * cos(phi);
  y_out = r * sin(phi);
}

void make_orthonormals(vector N, output vector a, output vector b)
{
  if(N[0] != N[1] || N[0] != N[2]) a = cross(vector(1, 1, 1), N);
  else a = cross(vector(-1, 1, 1), N);
  
  a = normalize(a);
  b = cross(N, a);
}

vector sample_cos_hemisphere(vector N, float randu, float randv)
{
  vector T, B;
    
  make_orthonormals(N, T, B);
  to_unit_disk(randu, randv, randu, randv);
  float costheta = sqrt(max(1.0 - randu * randu - randv * randv, 0.0));

  return randu * T + randv * B + costheta * N;
}

shader node_occlusion2(
  color Effect = color(0),
  color No_Effect = color(1),
  int Mode = 0, /* 0: Concave (AO) 1:Convex (Wear) 2:Both */
  int InvertEffect = 0,
  float Distance = 0.2,
  int Samples = 1,
  output color Color = 0,
  output float Fac = 0,
  output normal outNormal = N
)
{
  int i, rng;
  float f, randu, randv, ray_t, hits = 0;
  vector ray_P, ray_R;
  normal hit_normal = N;
  float hit_dist;

  f = fmod(cellnoise(P*123456.0), 1.0);
  rng_seed(rng, int(f * 21374647));
  
  for(i = 0; i &lt; Samples; i++) 
  { randu = rng_uniform(rng);
    randv = rng_uniform(rng);
       
    ray_P = P;
    ray_R = sample_cos_hemisphere(-N, randu, randv);
    ray_t = Distance;
    
        
    if (!Mode)
    { if(trace(ray_P, -ray_R, &quot;maxdist&quot;, ray_t)) {
            hits += 1.0;
            int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
            outNormal = outNormal + (hit_normal);
        }
    }
    else if (Mode == 1)
    { if(trace(ray_P, ray_R, &quot;maxdist&quot;, ray_t)) {
           hits += 1.0;
           int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
           outNormal = outNormal - (hit_normal);
        }
    }
    else { 
        if(trace(ray_P, -ray_R, &quot;maxdist&quot;, ray_t)) {
            hits += 1.0;
            int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
            outNormal = outNormal + (hit_normal);
        }
        if(trace(ray_P, ray_R, &quot;maxdist&quot;, ray_t)) {
           hits += 1.0;
           int HitNormal = getmessage (&quot;trace&quot;, &quot;N&quot;, hit_normal);
           outNormal = outNormal - (hit_normal);
        }
    } 
  }
  Fac = 1.0 - (hits/Samples);
  if(InvertEffect) Color = mix(No_Effect, Effect, Fac);
  else Color = mix(Effect, No_Effect, Fac);
  outNormal = normalize(outNormal);
}

Here's a render of many spheres that are all separate objects.  It works well in some cases but not always.