The computation of polarization integrals proceeds in two steps. First, an array of all pairs of basis functions has to be generated. The upper bound may be used to skip negligible pairs. Each pair (instance of class PrimitivePair) contains the information about the position, total angular momentum, contraction coefficient and orbital exponent of the bra and ket primitive shells (both instances of the Primitive class). In the second step, the pair list is transferred to the device memory, and the kernel function polarization_prim_pairs_kernel <double, k, mx, my, mz, q> (pairs, npair, buffer, alpha); processes all npair pairs in parallel putting the output into a buffer array. Since each primitive shell with angular momentum l contains N L (l) = (l + 1)(l + 2)/2 Cartesian angular momentum components (all combinations of powers n i such that n x + n y + n z = l), each pair of primitive shells i and j generates an integral block of size N L (l i ) × N L (l j ) in the output buffer. For more information refer to the source code of the test cases and examples.