A computational scheme has been developed within the framework of Rayleigh-Schrodinger perturbation theory to evaluate nonlinear interaction energy contributions for a molecule in the presence of an external spatially uniform, time-independent magnetic field. Terms connected with the fourth power of the perturbing field, representing the fourth-rank hypermagnetizabilities of five small molecules, have been evaluated at the coupled Hartree-Fock level of accuracy within the conventional common-origin approach. Gaugeless basis sets of increasing size and flexibility have been employed in a numerical test, adopting two different coordinate systems to estimate the degree of convergence of theoretical tensor components.