Spin glasses are spin-lattice models with quenched disorder and frustration. The mean field long-range Sherrington-Kirkpatrick (SK) model was solved by Parisi and displays replica symmetry breaking (RSB), but the more realistic short-range Edwards-Anderson (EA) model is still not solved. Whether the EA spin glass phase has many pairs of pure states as described by the RSB scenario or a single pair of pure states as described by two-state scenarios such as the droplet/scaling picture is not known yet. Rigorous analytical calculations of the EA model are not available at present and efficient numerical simulations of spin glasses are crucial in making progresses in the field. In addition to being a prototypical example of a classical disordered system with many interesting equilibrium as well as nonequilibrium properties, spin glasses are of great importance across multiple fields from neural networks, various combinatorial optimization problems to benchmark tests of quantum annealing machines. Therefore, it is important to gain a better understanding of the spin glass models. In an effort to do so, our work has two main parts, one is to develop an efficient algorithm called population annealing Monte Carlo and the other is to explore the physics of spin glasses using thermal boundary conditions. We present a full characterization of the population annealing algorithm focusing on its equilibration properties and make a systematic comparison of population annealing with two well established simulation methods, parallel tempering Monte Carlo and simulated annealing Monte Carlo. We show numerically that population annealing is similar in performance to parallel tempering, each has its own strengths and weaknesses and both algorithms outperform simulated annealing in combinatorial optimization problems. In thermal boundary conditions, all eight combinations of periodic vs antiperiodic boundary conditions in the three spatial directions appear in the ensemble with their respective Boltzmann weights, thus minimizing finite-size effects due to domain walls. With thermal boundary conditions and sample stiffness extrapolation, we show that our data is consistent with a two-state picture, not the RSB picture for the EA model. Thermal boundary conditions also provides an elegant way to study the phenomena of temperature chaos and bond chaos, and our results are again in agreement with the droplet/scaling scenario.