What is the formula used to find out the chemical potential of a molecule or atom?
The change of the expectation value of energy with respect to electron number, ## \partial <E>/\partial n_e ##.
I think this might be a bit more complicated. Since actual atoms and molecules have integer numbers of electrons, this derivative cannot strictly be defined. Also, a single molecule is a microscopic object, so should not be described by a statistical/thermodynamic quantity like a chemical potential.
That being said, people still offer definitions like μ = 1/2 (IP + EA) (where IP and EA are the first ionization potential and first electron affinity, respectively), or in the context of mean-field theory, μ = 1/2 (ε_homo + ε_lumo) (where the εs are the orbital energies of the highest occupied and lowest unoccupied orbital).
 (there are ways to define them in so called conceptual density functional theory, but one has to jump through several hoops, and I personally think it is legitimate to disregard non-integer electrons...)
I think there is no conceptual problem in calculating the mean energy per atom in a statistical ensemble in which only the average number of electrons per atom is specified, e.g. a mixture of 90% neutral atoms and 10% of atoms carrying a single positive or negative charge. It should be quite obvious how this relates to ionisation energy or electron affinity. At T=0 and integer values of electrons per atoms mu is discontinuous. At finite temperatures one can justify
μ = 1/2 (IP + EA) as always both positive and negatively charged atoms will be in the ensemble.
I agree that 1/2(IP+EA) it is a sensible choice for thermodynamic ensembles at T>0. But even then one still as to consider that, technically, there are also other IPs and EAs than the first, and they should enter in the thermodynamic sum. And one might have to consider ensembles of molecular geometries and molecular interactions, too.
I'm not saying that it is wrong to take or approximate μ like this, just that it is not something one has to accept as the 100% certain definition as which it sometimes comes across. To me it looks more like a sensible approximation derived from single atoms or molecules, which might then carry over to certain kinds of ensembles under some conditions. But this latter step is not often made explicit---people actually do talk about μs of single atoms and molecules, and define them like above.
I think that at T=0, as you said, the only ``really true'' fact is that μ is discontinuous, and that for finite T, technically one should actually set up and calculate the thermodynamic ensemble, from which one then would get μ as the mean energy derivative.
 Imo, this is more problematic than it looks like, because even in an ensemble it might not actually be the first IP/EA which is relevant, but some other one which would be favored by scattering or dynamic processes. If, for example, one had a complex with massive ligands all around a metal core, and the first EA/IP would actually refer to ionizations of this metal core, then the IPs/EAs refering to some outer ligand orbitals might be more important in the thermodynamic sum, even if higher in energy, because they could be triggered by collisions much more easily. (Yes, I am just making stuff up, but I think it's thinkable)
How about defining the chemical potential from the good old thermodynamic point of view to be the Gibbs free energy per atom (molecule). And then for the original poster we can suggest the ideal "classical" monoatomic gas formula:
μ=μ° + kBT ln (P/P0)
where kB is Boltzmann constant, T is the temperature, P is the partial pressure, P0 is a reference pressure and μ° a reference chemical potential.
This easy triggering would be a kinetic effect and should play no role in thermodynamic considerations.
You are right in that higher ionisation potentials have to be included in principle, however I think in praxis they are suppressed by a very small Boltzmann factor ##\exp(-(IP2-IP1)/kT)##.
Separate names with a comma.