To make efficient use of the GPUs, different kernels need to be implemented to handle different combinations of angular momenta in the ERIs. For the J matrix, where the [µν|λσ] integrals need to be evaluated, we need ss, sp, pp, etc. for the µν indices and same goes for the λσ indices. This results in nine Coulomb kernels for up to p-functions, and in our case 36 for up to d-functions. For the K matrix, where the [µλ|νσ] integrals need to be evaluated, we need ss, sp, pp, etc. for the µν indices, while for the λσ indices we need to take into account all combinations including ss, sp, ps, pp, etc. This leads to twelve exchange kernels for up to p-functions, and in our case 54 for up to d-functions. Note that we only take into account the upper triangular combinations for the µν indices when evaluating J µν and K µν , since we only need to deal with symmetric and anti-symmetric density matrices. Also, for anti-symmetric density matrices we skip the formation of the J matrix. After forming the J and K matrices in the primitive basis, the transformation to contracted basis is carried out on the node CPU.