Development and Assessment of Pair Natural Orbitals Coupled Cluster Quadratic Response Properties

Abstract

The acceleration of quantum chemical methods has enabled the study of systems beyond small molecules, extending into the regime of large and chemically relevant molecules and materials. Through the incorporation of reduced-scaling techniques, such as local pair natural orbitals (LPNOs), these advancements allow us to match experimental results and probe the fundamental behavior of complex systems. However, while significant progress has been made for ground-state energetics, such as reaction pathways and binding energies, much less progress has been achieved for molecular properties associated with a system's interaction with electromagnetic fields. Advancing this area would enable accurate simulations of light–matter phenomena, such as nonlinear optical (NLO) processes, which are central to the development of optoelectronic and photonic devices. In the context of computing nonlinear optical properties, this work presents a study on the efficacy of LPNO-based coupled cluster (CC) methods for calculating static electric-field properties—including dipole moments, polarizabilities, and first hyperpolarizabilities via numerical differentiation. We find that constructing pair correlation domains independently for each electric field displacement introduces inconsistencies in the dimensionality and content of the CC wavefunction. This leads to significant deviations from conventional CC results, particularly for higher-order properties such as hyperpolarizabilities (third-order energy derivatives). These findings suggest that LPNO-type methods are fundamentally unsuitable for use with numerical differentiation when targeting nonlinear properties. As an alternative, we explore a response theory-based approach. By employing perturbation-aware PNOs (PNO++) and a combined variant (cPNO++) that includes the external perturbation in the construction of the pair correlation domains, we demonstrate that LPNO-CC response methods offer a promising route for reducing computational cost. Applying these methods to quadratic response properties—closely related to hyperpolarizabilities—we observe improved performance over standard PNOs, enabling more aggressive truncation of the wave function space while maintaining acceptable accuracy relative to canonical CC response theory. These results establish the viability of tailored PNO-based approximations for computing nonlinear optical properties, and we provide a Python-based implementation of LPNO-CCSD quadratic response theory as a foundation for future production-level development.