This is really a difficult problem, but possibly not for the reason you imagine. The following naif criterion seems at first highly appropriate:
take any pair of stars, subtract the center-of-mass motion, compute the kinetic ($T$) and gravitational ($W$) energies, and check whether
$T + W < 0$.
If so, the pair is bound, otherwise it is not.
How can this criterion fail? I can think of three ways:
Even in an uncrowded field, there may be three or more stars bound. It is well-known that multiple systems are unstable, so establishing which of the stars will remain bounded requires more analysis. Normally, the two heaviest stars are those that remain bounded, because the lighter ones gain energy at the expense of the heavy double, to escape the bound system.
A supposedly unbound pair may later turn out to be bound after all, because of loss of angular momentum. It is well-know that the centrifugal barrier, under angular momentum conservation from the original cloud in which the stars form, is so high that nearly no stars should form. Since gravity conserves the overall amount of angular momentum, there must be non-gravitational forces at work that remove it from the collapsing gas. No one knows for sure what is responsible for these torques, even though everyone claims that they must be magnetic in origin.
Thus it can happen that a given pair has, initially, so much angular momentum as to violate the above condition, but it can also happen that angular momentum removal is still occurring (at the time your simulations end) so that more binary pairs will form than is apparent from your simulations. Since this involves processes which are, as of now, largely speculative, no one can claim to know the exact time-evolution of these torques. This is a fundamental element of ignorance, on our part.
Lastly, and least likely, in crowded fields, interactions between the (transient) pair and passers-by may unbind the binary, at the expense of the kinetic energy of the passer-by. In the limit of fully-developed stars, this has been studied especially in connection with the evolution of binaries in globular clusters (GCs, see Lyman Spitzer's Dynamical Evolution of Globular Clusters, especially Ch. 6). But you should keep in mind that in GCs two circumstances enhance the relevance of these interactions: the high stellar densities (higher than in Open Clusters), and the long time scales available, of the order of a Hubble time as compared to the few million years available before newly-formed stars walk away from their formation sites.
Still, one circumstance makes many-body encounters more effective in Open Clusters, the large size of forming stars: when a binary encounters a passer-by, energy may be dumped from the enterloper's kinetic motion into tidal heating of either star in the binary, to be later dissipated thru radiative mechanisms. This tidal heating slows down the interloper, hence it makes the contact last longer, and the impulse transferred to the stars in the binary larger. In other words, it facilitates unbinding the binary.
Of the above, objections 1 and 3 are easily kept under control if your simulations do not involve unrealistically crowded fields (after all, Open Clusters are barely bound, and do no last long!). But objection 2 is more fundamental, and harder to keep under control. Keep in mind that, in Nature, 1 out of two stellar systems observed from Earth is a binary, so just about any simulation resulting in many fewer binaries than this is likely to be tainted by an improper accounting of electromagnetic torques.