Multiscale Models for Light-Driven Processes

A Polarizable QM/MM Model That Combines the State-Averaged CASSCF and AMOEBA Force Field for Photoreactions in Proteins

This study presents the polarizable quantum mechanics/molecular mechanics (QM/MM) embedding of the state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular applications (AMOEBA) force field for the purpose of studying photoreactions in protein environments. We describe two extensions of our previous work that combine SA-CASSCF with AMOEBA water models, allowing it to be generalized to AMOEBA models for proteins and other macromolecules. First, we discuss how our QM/MM model accounts for the discrepancy between the direct and polarization electric fields that arises in the AMOEBA description of intramolecular polarization. A second improvement is the incorporation of link atom schemes to treat instances in which the QM/MM boundary goes through covalent bonds. A single-link atom scheme and double-link atom scheme are considered in this work, and we will discuss how electrostatic interaction, van der Waals interaction, and various kinds of valence terms are treated across the boundary. To test the accuracy of the link atom scheme, we will compare QM/MM with full QM calculations and study how the errors in ground state properties, excited state properties, and excitation energies change when tuning the parameters in the link atom scheme. We will also test the new SA-CASSCF/AMOEBA method on an elementary reaction step in NanoLuc, an artificial bioluminescence luciferase. We will show how the reaction mechanism is different when calculated in the gas phase, in polarizable continuum medium (PCM), versus in protein AMOEBA models.

Importance of Polarizable Embedding for Computing Optical Rotation: The Case of Camphor in Ethanol

Solvation effects on optical rotation are notoriously challenging to model for computational chemistry, as the specific rotatory power of a molecule can vary wildly going from apolar to polar or even protic solvents. To address such a problem, we present a polarizable embedding implementation of an electric and magnetic response property based on density functional theory and the AMOEBA polarizable force field, and apply such an implementation to the study of the optical rotation of camphor in ethanol. By comparing a continuum model, and electrostatic and polarizable embedding QM/MM models, we observe that accounting for the environment's polarization gives rise to not only a different quantitative prediction, in very good agreement with experiments for the QM/AMOEBA model, but also to a very different qualitative picture, with the values of the optical rotation computed along a classical molecular dynamics trajectory with electrostatic embedding being statistically uncorrelated to the ones obtained with the polarizable description.

Multiscale biomolecular simulations in the exascale era

The complexity of biological systems and processes, spanning molecular to macroscopic scales, necessitates the use of multiscale simulations to get a comprehensive understanding. Quantum mechanics/molecular mechanics (QM/MM) molecular dynamics (MD) simulations are crucial for capturing processes beyond the reach of classical MD simulations. The advent of exascale computing offers unprecedented opportunities for scientific exploration, not least within life sciences, where simulations are essential to unravel intricate molecular mechanisms underlying biological processes. However, leveraging the immense computational power of exascale computing requires innovative algorithms and software designs. In this context, we discuss the current status and future prospects of multiscale biomolecular simulations on exascale supercomputers with a focus on QM/MM MD. We highlight our own efforts in developing a versatile and high-performance multiscale simulation framework with the aim of efficient utilization of state-of-the-art supercomputers. We showcase its application in uncovering complex biological mechanisms and its potential for leveraging exascale computing.

A Vision for the Future of Multiscale Modeling

The rise of modern computer science enabled physical chemistry to make enormous progresses in understanding and harnessing natural and artificial phenomena. Nevertheless, despite the advances achieved over past decades, computational resources are still insufficient to thoroughly simulate extended systems from first principles. Indeed, countless biological, catalytic and photophysical processes require ab initio treatments to be properly described, but the breadth of length and time scales involved makes it practically unfeasible. A way to address these issues is to couple theories and algorithms working at different scales by dividing the system into domains treated at different levels of approximation, ranging from quantum mechanics to classical molecular dynamics, even including continuum electrodynamics. This approach is known as multiscale modeling and its use over the past 60 years has led to remarkable results. Considering the rapid advances in theory, algorithm design, and computing power, we believe multiscale modeling will massively grow into a dominant research methodology in the forthcoming years. Hereby we describe the main approaches developed within its realm, highlighting their achievements and current drawbacks, eventually proposing a plausible direction for future developments considering also the emergence of new computational techniques such as machine learning and quantum computing. We then discuss how advanced multiscale modeling methods could be exploited to address critical scientific challenges, focusing on the simulation of complex light-harvesting processes, such as natural photosynthesis. While doing so, we suggest a cutting-edge computational paradigm consisting in performing simultaneous multiscale calculations on a system allowing the various domains, treated with appropriate accuracy, to move and extend while they properly interact with each other. Although this vision is very ambitious, we believe the quick development of computer science will lead to both massive improvements and widespread use of these techniques, resulting in enormous progresses in physical chemistry and, eventually, in our society.

Excitation landscape of the CP43 photosynthetic antenna complex from multiscale simulations

Photosystem II (PSII), the principal enzyme of oxygenic photosynthesis, contains two integral light harvesting proteins (CP43 and CP47) that bind chlorophylls and carotenoids. The two intrinsic antennae play crucial roles in excitation energy transfer and photoprotection. CP43 interacts most closely with the reaction center of PSII, specifically with the branch of the reaction center (D1) that is responsible for primary charge separation and electron transfer. Deciphering the function of CP43 requires detailed atomic-level insights into the properties of the embedded pigments. To advance this goal, we employ a range of multiscale computational approaches to determine the site energies and excitonic profile of CP43 chlorophylls, using large all-atom models of a membrane-bound PSII monomer. In addition to time-dependent density functional theory (TD-DFT) used in the context of a quantum-mechanics/molecular-mechanics setup (QM/MM), we present a thorough analysis using the perturbed matrix method (PMM), which enables us to utilize information from long-timescale molecular dynamics simulations of native PSII-complexed CP43. The excited state energetics and excitonic couplings have both similarities and differences compared with previous experimental fits and theoretical calculations. Both static TD-DFT and dynamic PMM results indicate a layered distribution of site energies and reveal specific groups of chlorophylls that have shared contributions to low-energy excitations. Importantly, the contribution to the lowest energy exciton does not arise from the same chlorophylls at each system configuration, but rather changes as a function of conformational dynamics. An unexpected finding is the identification of a low-energy charge-transfer excited state within CP43 that involves a lumenal (C2) and the central (C10) chlorophyll of the complex. The results provide a refined basis for structure-based interpretation of spectroscopic observations and for further deciphering excitation energy transfer in oxygenic photosynthesis.

Deciphering Photoreceptors Through Atomistic Modeling from Light Absorption to Conformational Response

In this review, we discuss the successes and challenges of the atomistic modeling of photoreceptors. Throughout our presentation, we integrate explanations of the primary methodological approaches, ranging from quantum mechanical descriptions to classical enhanced sampling methods, all while providing illustrative examples of their practical application to specific systems. To enhance the effectiveness of our analysis, our primary focus has been directed towards the examination of applications across three distinct photoreceptors. These include an example of Blue Light-Using Flavin (BLUF) domains, a bacteriophytochrome, and the orange carotenoid protein (OCP) employed by cyanobacteria for photoprotection. Particular emphasis will be placed on the pivotal role played by the protein matrix in fine-tuning the initial photochemical event within the embedded chromophore. Furthermore, we will investigate how this localized perturbation initiates a cascade of events propagating from the binding pocket throughout the entire protein structure, thanks to the intricate network of interactions between the chromophore and the protein.

Visible Light Induced Exciton Dynamics and trans-to-cis Isomerization in Azobenzene Aggregates: Insights from Surface Hopping/Semiempirical Configuration Interaction Molecular Dynamics Simulations

Assemblies of photochromic molecules feature exciton states, which govern photochemical and photophysical processes in multichromophoric systems. Understanding the photoinduced dynamics of the assemblies requires nonadiabatic treatment involving multiple exciton states and numerous nuclear degrees of freedom, thus posing a challenge for simulations. In this work, we address this challenge for aggregates of azobenzene, a prototypical molecular switch, performing on-the-fly surface hopping calculations combined with semiempirical configuration interaction electronic structure and augmented with transition density matrix analysis to characterize exciton evolution. Specifically, we consider excitation of azobenzene tetramers in the nπ* absorption band located in the visible (blue) part of the electromagnetic spectrum, thus extending our recent work on dynamics after ππ* excitation corresponding to the ultraviolet region [Titov, J. Phys. Chem. C2023, 127, 13678-13688]. We find that the nπ* excitons, which are initially strongly localized by ground-state conformational disorder, undergo further (very strong) localization during short-time photodynamics. This excited-state localization process is extremely ultrafast, occurring within the first 10 fs of photodynamics. We observe virtually no exciton transfer of the localized excitons in the nπ* manifold. However, the transfer may occur via secondary pathways involving ππ* states or the ground state. Moreover, we find that the nπ* quantum yields of the trans-to-cis isomerization are reduced in the aggregated state.

First-Principles Calculations of Excited-State Decay Rate Constants in Organic Fluorophores

In this Perspective, we discuss recent advances made to evaluate from first-principles the excited-state decay rate constants of organic fluorophores, focusing on the so-called static strategy. In this strategy, one essentially takes advantage of Fermi's golden rule (FGR) to evaluate rate constants at key points of the potential energy surfaces, a procedure that can be refined in a variety of ways. In this way, the radiative rate constant can be straightforwardly obtained by integrating the fluorescence line shape, itself determined from vibronic calculations. Likewise, FGR allows for a consistent calculation of the internal conversion (related to the non-adiabatic couplings) in the weak-coupling regime and intersystem crossing rates, therefore giving access to estimates of the emission yields when no complex photophysical phenomenon is at play. Beyond outlining the underlying theories, we summarize here the results of benchmarks performed for various types of rates, highlighting that both the quality of the vibronic calculations and the accuracy of the relative energies are crucial to reaching semiquantitative estimates. Finally, we illustrate the successes and challenges in determining the fluorescence quantum yields using a series of organic fluorophores.

Interpolating Nonadiabatic Molecular Dynamics Hamiltonian with Bidirectional Long Short-Term Memory Networks

Essential for understanding far-from-equilibrium processes, nonadiabatic (NA) molecular dynamics (MD) requires expensive calculations of the excitation energies and NA couplings. Machine learning (ML) can simplify computation; however, the NA Hamiltonian requires complex ML models due to its intricate relationship to atomic geometry. Working directly in the time domain, we employ bidirectional long short-term memory networks (Bi-LSTM) to interpolate the Hamiltonian. Applying this multiscale approach to three metal-halide perovskite systems, we achieve two orders of magnitude computational savings compared to direct ab initio calculation. Reasonable charge trapping and recombination times are obtained with NA Hamiltonian sampling every half a picosecond. The Bi-LSTM-NAMD method outperforms earlier models and captures both slow and fast time scales. In combination with ML force fields, the methodology extends NAMD simulation times from picoseconds to nanoseconds, comparable to charge carrier lifetimes in many materials. Nanosecond sampling is particularly important in systems containing defects, boundaries, interfaces, etc. that can undergo slow rearrangements.

Surface hopping modeling of charge and energy transfer in active environments

An active environment is any atomic or molecular system changing a chromophore's nonadiabatic dynamics compared to the isolated molecule. The action of the environment on the chromophore occurs by changing the potential energy landscape and triggering new energy and charge flows unavailable in the vacuum. Surface hopping is a mixed quantum-classical approach whose extreme flexibility has made it the primary platform for implementing novel methodologies to investigate the nonadiabatic dynamics of a chromophore in active environments. This Perspective paper surveys the latest developments in the field, focusing on charge and energy transfer processes.

Machine Learning Exciton Hamiltonians in Light-Harvesting Complexes

We propose a machine learning (ML)-based strategy for an inexpensive calculation of excitonic properties of light-harvesting complexes (LHCs). The strategy uses classical molecular dynamics simulations of LHCs in their natural environment in combination with ML prediction of the excitonic Hamiltonian of the embedded aggregate of pigments. The proposed ML model can reproduce the effects of geometrical fluctuations together with those due to electrostatic and polarization interactions between the pigments and the protein. The training is performed on the chlorophylls of the major LHC of plants, but we demonstrate that the model is able to extrapolate well beyond the initial training set. Moreover, the accuracy in predicting the effects of the environment is tested on the simulation of the small changes observed in the absorption spectra of the wild-type and a mutant of a minor LHC.

State averaged CASSCF in AMOEBA polarizable water model for simulating nonadiabatic molecular dynamics with nonequilibrium solvation effects

This paper presents a state-averaged complete active space self-consistent field (SA-CASSCF) in the atomic multipole optimized energetics for biomolecular application (AMOEBA) polarizable water model, which enables rigorous simulation of non-adiabatic molecular dynamics with nonequilibrium solvation effects. The molecular orbital and configuration interaction coefficients of the solute wavefunction, and the induced dipoles on solvent atoms, are solved by minimizing the state averaged energy variationally. In particular, by formulating AMOEBA water models and the polarizable continuum model (PCM) in a unified way, the algorithms developed for computing SA-CASSCF/PCM energies, analytical gradients, and non-adiabatic couplings in our previous work can be generalized to SA-CASSCF/AMOEBA by properly substituting a specific list of variables. Implementation of this method will be discussed with the emphasis on how the calculations of different terms are partitioned between the quantum chemistry and molecular mechanics codes. We will present and discuss results that demonstrate the accuracy and performance of the implementation. Next, we will discuss results that compare three solvent models that work with SA-CASSCF, i.e., PCM, fixed-charge force fields, and the newly implemented AMOEBA. Finally, the new SA-CASSCF/AMOEBA method has been interfaced with the ab initio multiple spawning method to carry out non-adiabatic molecular dynamics simulations. This method is demonstrated by simulating the photodynamics of the model retinal protonated Schiff base molecule in water.

Assessing the quality of QM/MM approaches to describe vacuo-to-water solvatochromic shifts

The performance of different quantum mechanics/molecular mechanics embedding models to compute vacuo-to-water solvatochromic shifts is investigated. In particular, both nonpolarizable and polarizable approaches are analyzed and computed results are compared to reference experimental data. We show that none of the approaches outperform the others and that errors strongly depend on the nature of the molecular transition to be described. Thus, we prove that the best choice of embedding model highly depends on the molecular system and that the use of a specific approach as a black box can lead to significant errors and, sometimes, totally wrong predictions.

Structural Disorder in Higher-Temperature Phases Increases Charge Carrier Lifetimes in Metal Halide Perovskites

Solar cells and optoelectronic devices are exposed to heat that degrades performance. Therefore, elucidating temperature-dependent charge carrier dynamics is essential for device optimization. Charge carrier lifetimes decrease with temperature in conventional semiconductors. The opposite, anomalous trend is observed in some experiments performed with MAPbI₃ (MA = CH₃NH₃⁺) and other metal halide perovskites. Using ab initio quantum dynamics simulation, we establish the atomic mechanisms responsible for nonradiative electron-hole recombination in orthorhombic-, tetragonal-, and cubic MAPbI₃. We demonstrate that structural disorder arising from the phase transitions is as important as the disorder due to heating in the same phase. The carrier lifetimes grow both with increasing temperature in the same phase and upon transition to the higher-temperature phases. The increased lifetime is rationalized by structural disorder that induces partial charge localization, decreases nonadiabatic coupling, and shortens quantum coherence. Inelastic and elastic electron-vibrational interactions exhibit opposite dependence on temperature and phase. The partial disorder and localization arise from thermal motions of both the inorganic lattice and the organic cations and depend significantly on the phase. The structural deformations induced by thermal fluctuations and phase transitions are on the same order as deformations induced by defects, and hence, thermal disorder plays a very important role. Since charge localization increases carrier lifetimes but inhibits transport, an optimal regime maximizing carrier diffusion can be designed, depending on phase, temperature, material morphology, and device architecture. The atomistic mechanisms responsible for the enhanced carrier lifetimes at elevated temperatures provide guidelines for the design of improved solar energy and optoelectronic materials.

Trajectory Surface Hopping for a Polarizable Embedding QM/MM Formulation

We present the implementation of trajectory surface-hopping nonadiabatic dynamics for a polarizable embedding QM/MM formulation. Time-dependent density functional theory was used at the quantum mechanical level of theory, whereas the molecular mechanics description involved the polarizable AMOEBA force field. This implementation has been obtained by integrating the surface-hopping program Newton-X NS with an interface between the Gaussian 16 and the Tinker suites of codes to calculate QM/AMOEBA energies and forces. The implementation has been tested on a photoinduced electron-driven proton-transfer reaction involving pyrimidine and a hydrogen-bonded water surrounded by a small cluster of water molecules and within a large water droplet.

Electron-phonon relaxation at the Au/WSe2 interface is significantly accelerated by a Ti adhesion layer: time-domain ab initio analysis

Thermal transport at nanoscale metal-semiconductor interfaces via electron-phonon coupling is crucial for applications of modern microelectronic, electro-optic and thermoelectric devices. To enhance the device performance, the heat flow can be regulated by modifying the interfacial atomic interactions. We use ab initio time-dependent density functional theory combined with non-adiabatic molecular dynamics to study how the hot electron and hole relaxation rates change on incorporating a thin Ti adhesion layer at the Au/WSe₂ interface. The excited charge carrier relaxation is much faster in Au/Ti/WSe₂ due to the enhanced electron-phonon coupling, rationalized by the following reasons: (1) Ti atoms are lighter than Au, W and Se atoms and move faster. (2) Ti has a significant contribution to the electronic properties in the relevant energy range. (3) Ti interacts strongly with WSe₂ and promotes its bond-scissoring which causes Fermi-level pinning, making WSe₂ contribute to electronic properties around the Fermi level. The changes in the relaxation rates are more pronounced for excited electrons compared to holes because both relative and absolute Ti contributions to the electronic properties are larger above than below the Fermi level. The results provide guidance for improving the design of novel and robust materials by optimizing the heat dissipation at metal-semiconductor interfaces.

Evaluation of molecular photophysical and photochemical properties using linear response time-dependent density functional theory with classical embedding: Successes and challenges

Time-dependent density functional theory (TDDFT) based approaches have been developed in recent years to model the excited-state properties and transition processes of the molecules in the gas-phase and in a condensed medium, such as in a solution and protein microenvironment or near semiconductor and metal surfaces. In the latter case, usually, classical embedding models have been adopted to account for the molecular environmental effects, leading to the multi-scale approaches of TDDFT/polarizable continuum model (PCM) and TDDFT/molecular mechanics (MM), where a molecular system of interest is designated as the quantum mechanical region and treated with TDDFT, while the environment is usually described using either a PCM or (non-polarizable or polarizable) MM force fields. In this Perspective, we briefly review these TDDFT-related multi-scale models with a specific emphasis on the implementation of analytical energy derivatives, such as the energy gradient and Hessian, the nonadiabatic coupling, the spin–orbit coupling, and the transition dipole moment as well as their nuclear derivatives for various radiative and radiativeless transition processes among electronic states. Three variations of the TDDFT method, the Tamm–Dancoff approximation to TDDFT, spin–flip DFT, and spin-adiabatic TDDFT, are discussed. Moreover, using a model system (pyridine–Ag20 complex), we emphasize that caution is needed to properly account for system–environment interactions within the TDDFT/MM models. Specifically, one should appropriately damp the electrostatic embedding potential from MM atoms and carefully tune the van der Waals interaction potential between the system and the environment. We also highlight the lack of proper treatment of charge transfer between the quantum mechanics and MM regions as well as the need for accelerated TDDFT modelings and interpretability, which calls for new method developments.

The atomistic modeling of light-harvesting complexes from the physical models to the computational protocol

The function of light-harvesting complexes is determined by a complex network of dynamic interactions among all the different components: the aggregate of pigments, the protein, and the surrounding environment. Complete and reliable predictions on these types of composite systems can be only achieved with an atomistic description. In the last few decades, there have been important advances in the atomistic modeling of light-harvesting complexes. These advances have involved both the completeness of the physical models and the accuracy and effectiveness of the computational protocols. In this Perspective, we present an overview of the main theoretical and computational breakthroughs attained so far in the field, with particular focus on the important role played by the protein and its dynamics. We then discuss the open problems in their accurate modeling that still need to be addressed. To illustrate an effective computational workflow for the modeling of light harvesting complexes, we take as an example the plant antenna complex CP29 and its H111N mutant.

Interpolating Nonadiabatic Molecular Dynamics Hamiltonian with Inverse Fast Fourier Transform

Nonadiabatic (NA) molecular dynamics (MD) allows one to investigate far-from-equilibrium processes in nanoscale and molecular materials at the atomistic level and in the time domain, mimicking time-resolved spectroscopic experiments. Ab initio NAMD is limited to about 100 atoms and a few picoseconds, due to computational cost of excitation energies and NA couplings. We develop a straightforward methodology that can extend ab initio quality NAMD to nanoseconds and thousands of atoms. The ab initio NAMD Hamiltonian is sampled and interpolated along a trajectory using a Fourier transform, and then, it is used to perform NAMD with known algorithms. The methodology relies on the classical path approximation, which holds for many materials and processes. To achieve a complete ab initio quality description, the trajectory can be obtained using an ab initio trained machine learning force field. The method is demonstrated with charge carrier trapping and relaxation in hybrid organic-inorganic and all-inorganic metal halide perovskites that exhibit complex dynamics and are actively studied for optoelectronic applications.

Surface Hopping Dynamics with the Frenkel Exciton Model in a Semiempirical Framework

We present an implementation of the Frenkel exciton model in the framework of the semiempirical floating occupation molecular orbitals-configuration interaction (FOMO-CI) electronic structure method, aimed at simulating the dynamics of multichromophoric systems, in which excitation energy transfer can occur, by a very efficient approach. The nonadiabatic molecular dynamics is here dealt with by the surface hopping method, but the implementation we proposed is compatible with other dynamical approaches. The exciton coupling is computed either exactly, within the semiempirical approximation considered, or by resorting to transition atomic charges. The validation of our implementation is carried out on the trans-azobenzeno-2S-phane (2S-TTABP), formed by two azobenzene units held together by sulfur bridges, taken as a minimal model of multichromophoric systems, in which both strong and weak exciton couplings are present.

Protein Response Effects on Cofactor Excitation Energies from First Principles: Augmenting Subsystem Time-Dependent Density-Functional Theory with Many-Body Expansion Techniques

We investigate the possibility of describing protein response effects on a chromophore excitation by means of subsystem time-dependent density-functional theory (sTDDFT) in combination with a many-body expansion (MBE) approach. While sTDDFT is in principle intrinsically able to include such contributions, addressing cofactor excitations in protein models or entire proteins with full environment-response treatments is currently out of reach. Taking different model structures of the green fluorescent protein (GFP) and bovine rhodopsin as examples, we demonstrate that an embedded-MBE approach based on sTDDFT in its simplest version leads to a good agreement of the predicted protein response effect already at second order. To reproduce reference response effects from nonsubsystem TDDFT calculations quantitatively (error ≤ 5%), however, a third- or even fourth-order MBE may be required. For the latter case, we explore a selective inclusion of fourth-order terms that drastically reduces the computational burden. In addition, we demonstrate how this sTDDFT-MBE treatment can be utilized as an analysis tool to identify residues with dominant response contributions. This, in turn, can be employed to arrive at smaller structural models for light-absorbing proteins, which still feature all of the main characteristics in terms of photoresponse properties.

Ultrafast Transient Infrared Spectroscopy of Photoreceptors with Polarizable QM/MM Dynamics

Ultrafast transient infrared (TRIR) spectroscopy is widely used to measure the excitation-induced structural changes of protein-bound chromophores. Here, we design a novel and general strategy to compute TRIR spectra of photoreceptors by combining μs-long MM molecular dynamics with ps-long QM/AMOEBA Born-Oppenheimer molecular dynamics (BOMD) trajectories for both ground and excited electronic states. As a proof of concept, the strategy is here applied to AppA, a blue-light-utilizing flavin (BLUF) protein, found in bacteria. We first analyzed the short-time evolution of the embedded flavin upon excitation revealing that its dynamic Stokes shift is ultrafast and mainly driven by the internal reorganization of the chromophore. A different normal-mode representation was needed to describe ground- and excited-state IR spectra. In this way, we could assign all of the bands observed in the measured transient spectrum. In particular, we could characterize the flavin isoalloxazine-ring region of the spectrum, for which a full and clear description was missing.

Efficient Open-Source Implementations of Linear-Scaling Polarizable Embedding: Use Octrees to Save the Trees

We present open-source implementations of the linear-scaling fast multipole method (FMM) within the polarizable embedding (PE) model for efficient treatment of large polarizable environments in computational spectroscopy simulations. The implementations are tested for accuracy, efficiency, and usability on model systems as well as more realistic biomolecular systems. We explain how FMM parameters affect the calculation of molecular properties and show that PE calculations employing FMM can be carried out in a black-box manner. The efficiency of the linear-scaling approach is demonstrated by simulating the UV/vis spectrum of a chromophore in an environment of more than 1 million polarizable sites. Our implementations are interfaced to several open-source quantum chemistry programs, making computational spectroscopy simulations within the PE model and FMM available to a large variety of methods and a broad user base.