  The combinatorial explosion of empirical parameters in tens of thousands
presents a tremendous challenge for extended Debye-H\"uckel models to calculate
activity coefficients of aqueous mixtures of most important salts in chemistry.
The explosion of parameters originates from the phenomenological extension of
the Debye-H\"uckel theory that does not take steric and correlation effects of
ions and water into account. In contrast, the Poisson-Fermi theory developed in
recent years treats ions and water molecules as nonuniform hard spheres of any
size with interstitial voids and includes ion-water and ion-ion correlations.
We present a Poisson-Fermi model and numerical methods for calculating the
individual or mean activity coefficient of electrolyte solutions with any
arbitrary number of ionic species in a large range of salt concentrations and
temperatures. For each activity-concentration curve, we show that the
Poisson-Fermi model requires only three unchanging parameters at most to well
fit the corresponding experimental data. The three parameters are associated
with the Born radius of the solvation energy of an ion in electrolyte solution
that changes with salt concentrations in a highly nonlinear manner.
