Variational calculations employing explicitly correlated Gaussian basis functions have been performed for the ground state of the boron monohydride molecule (BH) and for the boron atom (B). Up to 2000 Gaussians were used for each system. The calculations did not assume the Born-Oppenheimer (BO) approximation. In the optimization of the wave function, we employed the analytical energy gradient with respect to the Gaussian exponential parameters. In addition to the total nonrelativistic energies, we computed scalar relativistic corrections (mass-velocity and Darwin). With those added to the total energies, we estimated the dissociation energy of BH. The non-BO wave functions were also used to compute some expectation values involving operators dependent on the interparticle distances.