  This article presents an optimized algorithm and implementation for
calculating resolution-of-the-identity Hartree-Fock (RI-HF) energies and
analytic gradients using multiple Graphics Processing Units (GPUs). The
algorithm is especially designed for high throughput \emph{ab initio} molecular
dynamics simulations of small and medium size molecules (10-100 atoms). Key
innovations of this work include the exploitation of multi-GPU parallelism and
a workload balancing scheme that efficiently distributes computational tasks
among GPUs. Our implementation also employs techniques for symmetry
utilization, integral screening and leveraging sparsity to optimize memory
usage and computational efficiency. Computational results show that the
implementation achieves significant performance improvements, including over
$3\times$ speedups in single GPU AIMD throughput compared to previous
GPU-accelerated RI-HF and traditional HF methods. Furthermore, utilizing
multiple GPUs can provide super-linear speedup when the additional aggregate
GPU memory allows for the storage of decompressed three-center integrals.
Additionally, we report strong scaling efficiencies for systems up to 1000
basis functions and demonstrate practical applications through extensive
performance benchmarks on up to quadruple-$\zeta$ primary basis sets, achieving
floating-point performance of up to 47\% of the theoretical peak on a
4$\times$A100 GPU node.
