The combinatorial explosion of empirical parameters in tens of thousands presents a tremendous challenge for extended Debye-Hückel models to calculate activity coefficients of aqueous mixtures of the most important salts in chemistry. The explosion of parameters originates from the phenomenological extension of the Debye-Hückel theory that does not take steric and correlation effects of ions and water into account. By 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.