To evaluate the Coulomb interaction matrix in density functional theory (or Hartree-Fock theory) within a finite Gaussian type basis set, explicit evaluation of the two-electron integrals is not necessary. Instead, modified expressions can be developed directly for J matrix contributions from entire shell quartets of two-electron integrals. This can substantially reduce the computational cost. In this work we explore the use of the well-known McMurchie-Davidson recurrence relations to directly form J matrix contributions. Floating point operation counts and computer timings indicate that this approach yields significant improvements (typically a factor of two) over the previous J matrix engine.