This paper describes the extension of a recently developed numerical solver for the Landau-Lifshitz Navier-Stokes (LLNS) equations to binary mixtures in three dimensions. The LLNS equations incorporate thermal fluctuations into macroscopic hydrodynamics by using white-noise fluxes. These stochastic PDEs are more complicated in three dimensions due to the tensorial form of the correlations for the stochastic fluxes and in mixtures due to couplings of energy and concentration fluxes (e.g., Soret effect). We present various numerical tests of systems in and out of equilibrium, including time-dependent systems, and demonstrate good agreement with theoretical results and molecular simulation