In the rgeos package, gNearestPoints takes two geometries and returns the two closest points of the two geometries.
If there are multiple closest points, it seems to return one pair of them (for example if you do gNearestPoints(g1,g1) from the example in the help). I suspect if you test it on two line geometries that intersect you'll get one of the intersection points (two polygon geometries that overlap have a set of nearest points in the overlap region and you'll get one of those, probably on one of the boundaries).
The rgeos packages calles GEOSNearestPoints_r in GEOS to do the work.
Questions: "Should it, like in rgeos, only work for a pair of sfg (simple feature geometry) objects, not for lists of them"
See here for the examples.
Most helpful comment
See here for the examples.