An automatic code generated C++/HIP/CUDA implementation of the (auxiliary) Fock, or Kohn-Sham, matrix construction for execution in GPU-accelerated hardware environments is presented. The module is developed as part of the quantum chemistry software package VeloxChem, employing localized Gaussian atomic orbitals. The performance and scaling characteristics are analyzed in view of the specific requirements for self-consistent field optimization and response theory calculations. As an example, the electronic circular dichroism spectrum of a G-quadruplex is calculated at the level of time-dependent density functional theory in conjunction with the range-separated CAM-B3LYP exchange-correlation functional. Computational issues due to the high-density of states following the adoption of large-scale model systems are here bypassed with use of the complex polarization propagator approach. The origin of the negative Cotton effect in the long-wavelength onset of the experimental spectrum is elucidated in large-scale modeling and shown to be associated with the TTA nucleobase linkers in the G-quadruplex.