Nonlinear optical properties of molecular materials are governed by aggregation and cooperative effects in the condensed phase, requiring first-principles simulations on large molecular assemblies to achieve a realistic description of experiment. While modern electronic-structure theory implementations mitigate key computational challenges through low-memory two-electron integral engines and multifrequency response solvers, practical calculations of third-order nonlinear optical response in extended systems remain prohibitively expensive. The dominant bottleneck arises because conventional approaches require evaluating many independent components of the rank-four hyperpolarizability tensor, each of which is costly, making orientational averaging prohibitively expensive. Here, we introduce an analytic tensor-averaged formulation of cubic response theory within time-dependent Kohn-Sham density functional theory that removes this spatial bottleneck. By performing the orientational average directly at the level of perturbed densities and transformed Fock matrices, the full Cartesian second-order hyperpolarizability tensor is never constructed. Instead, the isotropic second-order hyperpolarizability is obtained directly. Exploiting the linearity of the Fock-matrix construction, this approach eliminates redundant spatial components and reduces the number of required exchange-correlation kernel integrations in the electronic quartic- and cubic-Hessian contractions by 72 and 90%, respectively. The resulting methodology enables cubic-response simulations of third-harmonic generation in large molecular aggregates with exchange-correlation functionals spanning multiple rungs of Jacob's ladder. Applications to oligo(thiophene-benzothiadiazole) (OTBP) clusters reveal pronounced aggregation-induced dampening of the third-harmonic response and demonstrate scalability to systems with nearly 5800 contracted basis functions. The present framework provides a quantitative tool for disentangling intrinsic molecular effects from supramolecular and morphological contributions, and for assessing the potential for further performance gains through crystal-structure engineering. In this way, theory can directly guide experimental efforts by identifying packing motifs that maximize cooperative enhancement and those that are detrimental. By resolving the missing spatial-domain optimization in the cubic response theory, the present tensor-averaged formulation shifts the dominant computational cost back to the frequency-dependent response solver, enabling routine material-level modeling of cooperative third-order nonlinear optical phenomena.
Ahmadzadeh et al. (Tue,) studied this question.