Focused ultrasound is a non-invasive, non-ionizing technology with great potential for various clinical applications, including thermal ablation of tumors, targeted drug delivery, and neuromodulation. Focused ultrasound uses ultrasound energy to treat tissue deep in the body. Optimizing treatment parameters to achieve desired clinical outcomes while minimizing adverse effects remains a significant challenge. Computational simulations are powerful tools to address this challenge, develop patient-specific treatment plans and general safety guidelines, and optimize ultrasound transducers. This study presents the development of an open-source Python library, named OptimUS, for calculating ultrasound wave propagation in large computational domains in 3D using the boundary element method, specifically for focused ultrasound applications. The numerical calculations only require surface meshes at the scatterers' interfaces to define the model's geometry. Also, the computations are fast and accurate for high-frequency waves through materials with high contrast in density and speed of sound. An intercomparison exercise supports the fidelity of the simulations. Finally, simulations using anatomical models for abdominal applications of focused ultrasound reliably show the aberration of the focus from reflections by ribs and the presence of prefocal hotspots due to the lensing effect of fat layers.