We report the first ever non-Born-Oppenheimer variational calculations on the ground state of a four electron molecular system. The basis set used in the calculations consists of explicitly correlated Gaussians multiplied by powers of the internuclear distance. To accelerate the optimization of the many nonlinear variational parameters involved in the variational wave function, we performed the calculations on a cluster of Linux workstations using MPI and a parallel implementation of the formulas. Results for the nonadiabatic ground state energy of LiH, as well as expectation values for the kinetic and potential energies, the internuclear and square of the internuclear distance, the virial coefficient, and the square of the energy gradient norm are reported. (C) 2001 American Institute of Physics.