fromage: A library for the study of molecular crystal excited states at the aggregate scale

Abstract The study of photoexcitations in molecular aggregates faces the twofold problem of the increased computational cost associated with excited states and the complexity of the interactions among the constituent monomers. A mechanistic investigation of these processes requires the analysis of the intermolecular interactions, the effect of the environment, and 3D arrangements or crystal packing on the excited states. A considerable number of techniques have been tailored to navigate these obstacles; however, they are usually restricted to in‐house codes and thus require a disproportionate effort to adopt by researchers approaching the field. Herein, we present the FRamewOrk for Molecular AGgregate Excitations (fromage), which implements a collection of such techniques in a Python library complemented with ready‐to‐use scripts. The program structure is presented and the principal features available to the user are described: geometrical analysis, exciton characterization, and a variety of ONIOM schemes. Each is illustrated by examples of diverse organic molecules in condensed phase settings. The program is available at https://github.com/Crespo-Otero-group/fromage.


| INTRODUCTION
The computational study of photochemistry in molecular condensed phases represents a notoriously difficult problem. The increased computational cost associated with accurate excited state methods is compounded by the need for larger system sizes in these materials.
In the case of molecular crystals, bound by van der Waals (vdW) forces, photochemical phenomena are thought not to delocalize beyond the immediate vicinity of the absorbing molecule. [1] This renders full-unit cell's periodic excited-state property calculations less pertinent than they would initially seem. In fact, the structural features, which are relevant to these systems, are of the scale of the aggregate, ranging from monomer conformation and dimer packing all the way to molecular cluster shape.
There is therefore a potential overlap in methodology between finite size molecular aggregates and periodic molecular crystals when dealing with excited state processes. The traditional approach to investigating molecular crystals-using periodic electronic structure codes-should thus be complemented by adequately extending molecular aggregate methods for the excited state.
The generation and analysis of these aggregate nuclear configurations is often separated from the codes that produce their corresponding electronic structure. While the former tasks are typically less computationally demanding than the latter one and therefore might be added as an auxiliary feature to a larger code (e.g., Crystal17 [2] ), they often require a degree of flexibility which situates them more comfortably in the realm of the programming library.
Indeed the bridging of scales between the periodic crystal, finite cluster, dimer, and monomer are challenging to generalize due to the formidable conformational variety of organic molecules.
An optimal strategy is therefore to provide modular tools for the investigator to tailor to their system, accompanied by ready-made scripts that compile those tools for use in ubiquitous cases. In this way, nonexpert users are able to use the program's principal features with a reliable degree of robustness, while more comfortable users can repurpose and extend the code to better suit fringe cases.
There exists a variety of computational chemistry scripting libraries for different tasks, with the ecosystem in rapid development. Their promise is to optimize and standardize the research workflow for increasingly specialized operations thanks to robust, tailored tools.
The Cambridge Structural Database Python API [3] focuses on crystallographic property analysis, accessing its associated database. The Atomic Simulation Environment (ASE) [4] specializes in interfacing with numerous electronic structure codes and communicates with them via Python scripting or a graphical user interface. RDKit [5] provides programming tools for general-purpose chemoinformatics and is itself used in numerous child libraries. Chemshell provides additive QM: MM interfaces between electronic structure and force field codes. [6] Libra is a library designed for the development of quantum and classical dynamics. [7] To our knowledge, there still lacks a library dedicated to the examination of photochemistry in molecular aggregates and crystals, exploiting the overlap in methodology between the two materials. To address this, we offer the FRamewOrk for Molecular AGgregate Excitations (fromage). fromage is a standalone Python library, accompanied by ready-to-use command line scripts destined to facilitate the study of molecular aggregates in the excited state.
They are summarized in Scheme 1. The program is tested for Python 2.7 and 3.6, though the authors do not guarantee backwards compatibility with Python 2 in future releases. For geometry manipulation routines, the program only relies on the unit cell's Cartesian structure and lattice vectors. From this information, the user can obtain unique dimer configurations, molecular clusters and general structural information.
All of the electronic structure calculations are performed by popular quantum chemistry programs. Currently, interfaces are provided to run calculations in DFTB+, [8] Gaussian, [9] Molcas, [10] and Turbomole. [11]  and Molecular mechanics (ONIOM) calculation while taking advantage of the diversity of modeling methods of several programs instead of only one. The electrostatic embedding in these calculated can be extended to include the Coulomb interaction of the whole crystal by fitting point charges to an Ewald potential. [12] fromage has been used to study various Quantum Mechanics in Quantum Mechanics (QM:QM 0 ) electrostatic embedding schemes for applications in photoactive molecular crystals, [12] the aggregationinduced emission process in propeller-shaped molecules, [13] and the design principles for proton transfer luminescent materials. [14] In this article, we present an overview of fromage and its capabilities. First, we describe the program structure, illustrated by basic usage examples. Then, we enumerate its principal features, namely, geometric analysis, exciton coupling evaluation, and ONIOM QM:QM 0 models.

| Principal classes
fromage makes use of two main classes in most of its operations, Atom and Mol. An Atom object is defined as a point in Cartesian space with associated physical properties. If the point is to represent an atom, supplying its element string provides standard properties such as covalent radius or atomic mass. The partial charge can also be specified, which can become useful in representing both atoms and point charges.
In practice, the user has little direct interaction with Atom objects.
Instead, they manipulate Mol objects, which contain a list of the former.
Mol extends common Python list methods via composition, allowing for intuitive appending, indexing, iterating, and so on, but also provides methods to manipulate molecular aggregate geometries and unit cells.
Mol objects can be created explicitly or generated from a typical geometry file such as .xyz. For instance, one may generate a molecular cluster and a supercell from unit cell information as illustrated in Figure 1.
Apart from the constituent atoms and lattice vectors, Mol also has attributes pertaining to the definition of a bond within the collection of atoms. Two atoms are said to be bonded when their distance falls below a certain threshold mol_thresh. This distance can be measured from nucleus to nucleus but also from the edge of the spheres of vdW radius or covalent radius. The method is selected by varying the mol_bonding attribute. Having this flexibility is required for highly distorted molecular geometries or diverse element combinations. Armed with the definition of a bond, the Mol class can single out covalently bonded complexes from an aggregate, generate molecular clusters from a single crystal, and detect atomic connectivity.
These tools are in and of themselves useful as a library for the Python literate user. However, several ready-made scripts are supplied for more complicated procedures and frequently required operations. Of the most practical use is perhaps fro_uc_tools.py, a command line script, which performs operations related to unit cells.
It is essential for comfortably communicating between periodic and finite systems which is a central concept in fromage. For instance, given a unit cell geometry file cell.xyz and a text file with the lattice vectors, the line: fro_uc_tools:py cell:xyz vectors− r 15 will produce a file cluster_out.xyz containing the geometry of a cluster of whole molecules where all atoms lie within a radius of 15 Å from the origin.
The other scripts are used to operate on unit cells, clusters, dimers, and monomers. They are discussed in Section 3 and represented in Figure 2 along with the core class structure.

| Dependencies
The most common calculation performed by the geometry manipulation routines is the evaluation of interatomic distances. Therefore, for an overall speedup, these distances are calculated in C++ and wrapped as a Python function using the Simplified Wrapper and Interface Generator (SWIG). Some other more involved operations F I G U R E 1 Example of use of fromage as a Python library. The three dimensional structures are shown in orthographic projection [Color figure can be viewed at wileyonlinelibrary.com] are sped up by using the numpy [15] library. Geometry optimization is carried out with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) implementation in scipy. [16] The program Ewald by Derenzo et al. [17,18] is used for the fitting of point charges to the Ewald potential. It is modified to allow the use of partial charges and redistribution with permission from the authors here https://github.com/Crespo-Otero-group/ Ewald.  hindrance. [19,20] It is, however, appropriate to begin such an investigation by obtaining computationally inexpensive and easily interpretable fea- In an aggregate of atoms a point belongs to the Voronoi volume of a given molecule if the atom it is nearest to belongs to that molecule. The application of Voroni cells [21] to molecular systems has successfully been used to characterize the geometry of condensed phases [22] and notably liquids. [23] This definition is refined by scaling the distance metric by the vdW radius of the atom; thus, for instance, assigning more space to oxygen than to hydrogen. For the calculation of distances, we employ a grid-based scheme, which makes it robust and allows us to choose an arbitrary resolution for the volume. [24,25] The resulting Voronoi volume V V can be compared to the sum of the An illustrative example is a derivative of 2 0 -hydroxychalcone (HC), which upon the addition of a methoxy group in para position with respect to the hydroxyl group turns off its emissive character in crystal form, bringing the fluorescence quantum yield from 0.32 to less than 0.01. [26,27] The new molecule is (E)-1-(2-hydroxy-5-methoxyphenyl)-

3-(4[dimethylamino]phenyl)prop-2-en-1-one (DAP). While HC exhibits
herringbone style packing, DAP has a complex unit cell structure whose steric constraints on the individual molecules are unclear at first glance. Figure 3 shows the difference in tightness of crystal packing between inequivalent monomers of DAP.

| Dimeric arrangement
Excitations in molecular aggregates are not guaranteed to remain confined to one absorbing monomer. Indeed, the electronic wavefunctions of the neighboring molecules may have enough overlap to produce intermolecular electronic interferences in the excited state. This can manifest in all sorts of photochemical processes and is central to exciton-governed mechanisms like charge transfer or singlet fission. In crystals, typical packing motifs like herringbone or sheet-like, produce a limited set of archetypal dimer arrangements such as edge-to-face or face-to-face which have generalizable excitonic behavior for chemically similar molecules. [14,28] For instance, in face-to-face aromatic systems, π-π interactions are a defining feature of the excitonic states.
Regardless of whether a researcher is investigating an amorphous cluster or a crystal structure, it is therefore informative to extract all of the significant dimers in the system and to quantify their geometrical arrangement in order to classify them under the principal dimer archetypes at a glance. In fromage, the user can extract the possible dimers whose distance falls below a given threshold. The intermolecular distance can be defined either as the centroid-to-centroid distance or as the nearest intermolecular atom pair distance. The latter can also be complemented by the vdW radii of the atoms. If the molecules are in lattice positions, the symmetry of the unit cell will produce groups of dimers identical up to a reflection or rotation. To filter out repeated configurations, all the intermolecular atomic distances are evaluated and sorted which provides a fingerprint for the dimer geometry.
Equivalent dimers are then defined as ones with the same fingerprint up to a root-mean-squared deviation threshold of 10 −4 Å by default.
The dimeric arrangement can then be characterized quantita- tively. An orthonormal set of principal, secondary, and tertiary axes is calculated for each constituent fragment, and the angles between same axes of two molecules can be associated to an archetypical dimer, effectively classifying the pair.
The generalized procedure to obtain characteristic vectors for a molecule is shown in Figure 4: 1. All atoms of the molecule are projected onto an averaged plane by singular value decomposition.
2. The two longest interatomic distances are identified, forming a quadrilateral ABCD such that AC > BD and AB is the longest side.
This imposes an arbitrary but consistent direction for the vectors.
This orthonormal set of axes presumes the significance of a plane which defines the shape of the molecule. Conjugated organic systems often have most atoms in one plane, and in other cases, the judicious elimination of extrenuous atoms from the analysis-a feature present in the code-can reduce the significant coordinates to a plane. To exclude a particular atom type from the calculation, only one atom need be specified, and others with the same identity can be automatically detected using the method outlined in Section 3.3.4.
If one same atom is involved in both largest interatomic distances, the procedure remains valid and the Points B and C become degenerate. In certain cases, the principal axis should simply be the vector connecting the largest interatomic distance. When this option is selected, the secondary axis becomes the perpendicular axis which lies on the averaged plane. In the general case where one wishes to extract the geometric information from a cluster of dimers clust.xyz, one should use the command: fro_dimer_tools:py clust:xyz, which will return the file dimers.dat containing all of the angles between dimers, the centroid distances and a proposed classification into different archetypical dimers. An additional slip angle is printed, to estimate the amount of face-to-face overlapping area giving rise to π-π interactions. It is defined as the smallest angle between the centroid-to-centroid axis and either tertiary axis of the constituent monomers, as is discussed by Dommett et al. [14].
The above method was used to investigate (2E)-3-(dimeth- This molecule exhibits lasing behavior with a fluorescence quantum yield of 0.77 in crystal form compared to 0.19 in polymethylmethacrylate film. [29] The crystal structure was optimized in Quan- because the process of charge separation and migration is fundamental to their mechanism. One approach is that of analyzing the one-electron transition density matrices associated with the excitation, which contains information about the migration of the charge during the process. [30] This method is implemented in TheoDORE. [31] Crespo-Otero and M. Barbatti [32] and Sen et al. [33] propose a different scheme to qualify the nature of the exciton by comparing the charge density distribution on fragments of a dimer before and after excitation. The original implementation is in a program named CALCDEN, which combines Perl and Fortran. The reimplementation in fromage is Python importable, making it an attractive alternative for use and extension depending on the user's preference. A Mulliken partition scheme can supply integrated orbital-specific densities located on one molecule (A): where S μν is the overlap integral between basis functions ϕ μ and ϕ ν , and c μk is the coefficient of basis function ϕ μ in molecular orbital k.
This density, ρ k A , can be used to produce two indices related to an excitation I: To illustrate the use of this feature, a dimer of perylene was extracted from its experimental crystal structure. [34] Its excited states were calculated in Gaussian16 using TD-ωB97X-D/

| Exciton coupling evaluation
The coupling associated with an exciton is a measure of the correlation between the individual isolated excitations within the exciton.
Evaluating it can give us a quantitative bearing on the importance of the many-body effects in the excited state process.
Kasha's exciton model is the initial approach, reducing the electronic densities of the fragments to point dipoles and comparing their relative geometry. [28] Under the point dipole approximation where μ n is the electronic transition dipole moment (TDM) for monomer n and R ij is the vector connecting the centroids of monomers i and j.
For a better spatial resolution of the electrostatic interactions, one can instead calculate the interaction between atomic transition charges (ATC): [35] where R k c is the position of atom c of monomer k with ATC q c and N k is the number of atoms of monomer k.
For more resolution, the transition electronic density itself can be used, yielding costly but exact Coulombic exciton coupling. Additional correction can be included, for example, by including the dielectric response of the environment as a polarizable continuum model as is done in the EXcitonic Analysis Tool (EXAT). [36] These models are purely Coulombic in nature and do not take into account the short-range interaction affecting the excited state behavior of adjacent molecules, which, for example, is dominant in chargetransfer states.
Aragó and Troisi have devised a procedure which evaluates the coupling J ij in a dimer by diabatising the Hamiltonian: [1] 1. Select an excited state property. For the sake of argument, we will use the TDM μ but anything bearing relation to the excited electronic density would be applicable. 6. Compute the diabatic Hamiltonian: The off-diagonal elements of the Hamiltonian the exciton coupling values J ij . More details can be found in ref. [1]. The diabatization method has already been employed in fromage to investigate the aggregate behavior of propeller-shaped emitters. [13] We have further extended it to calculate N-dimensional diabatic Hamiltonians, thus, for example, allowing for the calculation of pairwise exciton coupling within a trimer in the excited state, taking into account the influence of the third monomer.
An additional approximate method is that of exploiting the exciton energy splitting of the molecular excited state upon formation of a dimer. In the dimer, the S 1 and S 2 states are separated by twice the magnitude of the exciton coupling, provided that the individual constituent molecules are in perfect resonance. [37] Even in less symmetric cases, this approximation has been used with reasonable success. [14,38,39] fromage implements the half-gap method, the PDA Kasha This will print out the diabatic Hamiltonian where the three lower triangular off-diagonal elements correspond to the three exciton couplings of the trimer.
As an illustrative example, the crystal structure of 1,4-bis-(4-styryl-styryl)-benzene (4PV) was used to extract inequivalent dimers and evaluate their couplings. The structure of 4PV has been experimentally shown to contain six monomers, departing from the high symmetry of an ideal herringbone crystal due to variable slight rotation of the extreme phenyl rings. [40] First, the unit cell was optimized using PBE-D2 with a basis set cutoff of 50 Ry and a Monkhorst-Pack grid of 1 × 2 × 1 k-points as implemented in Quantum Espresso. [41] Then, the inequivalent dimers were detected by fromage using a centroid-centroid distance threshold of 7 Å. The TDMs of these dimers were calculated using TD-ωB97X-D/6-31G (d) in Gaussian16 [9] and processed in fromage in order to evaluate It is expected that in cofacial dimers with couplings of such magnitude, the short-range interactions should account for most of the coupling, making the use of a coupling scheme which accounts for exchange imperative. [42] Overall, the inclusion of the third molecule in the trimer diabatic Hamiltonian reduces the magnitude of the couplings by about 10 meV which is significant since it is in the order of the difference between certain edge-to-face and face-to-face values.
The irregularities in the herringbone packing produce two different face-to-face dimers with respective centroid-centroid distances of 6.05 and 6.15 Å. This relatively small increase in distance reduces the exciton coupling by 12 meV (14% of the average value), indicating that the natural packing of 4PV produces in dimers in an excitonically sensitive geometry.

| ONIOM flavors
Calculating an excited state electronic structure is often unavoidable when evaluating an absorption or emission energy, locating nonradiative pathways, or building a PES. This often proves to be challenging due to the environmental effects in such media.

| The ONIOM scheme
Many approaches have been developed and tailored specifically to reflect the condensed phase nature of these materials. For periodic systems, plane-wave basis ab initio calculations represent a tempting option due to their inherent treatment of the long-range interactions present in molecular crystals. However, due to the often local nature of excited state phenomena in molecular aggregates, it can be instructive to treat them as defects and not to periodically repeat the excitation throughout the material. [43,44] Furthermore, their treatment from a periodic perspective implies calculating the excited-state electronic structure of all explicit electrons in the unit cell, which, considering the relative asymmetry of molecular crystals, can represent many more electrons than are directly involved in the process. codes are not uncommon. GARLEEK [45] communicates between Gaussian and Tinker for subtractive QM:MM calculations. Chemshell [6,46] provides additive QM:MM formulations. The ASE [4] offers both additive ONIOM is a popular subtractive hybrid method scheme with the following energy equation: [47][48][49] where Region 1 is the active site, Region 2 is its vicinity and E i (n) is the energy of region n calculated at the i level of theory. In the usual ONIOM vernacular, Region 1 is the "model system" and Regions 1 and 2 are the "real system"; however, these terms are more appropriate in the general case where the inter-region boundary can cross bonds and link atoms are employed. [47] For excitations in molecular aggregates, 1 contains the molecules involved in the excitation and 2 a cluster of molecules encasing them. While ONIOM was originally formulated for QM:MM hybrid calculations, the extension to QM: QM 0 , where QM 0 is a lower accuracy method, soon followed. [48] This variant allows for a potentially more accurate descriptions of Region 2 and the inter-region interactions in exchange for an increased computational cost. The various quantum methods available and tested in fromage are listed in Table 1 along with their corresponding programs. For QM:QM 0 , the equation becomes: where EE stands for electrostatic embedding. For QM:MM, the point charges should simply correspond to the ones used in the force field; however, for QM:QM 0 , the choice is up for discussion. [12,50] In the embedding of E EE QM 0 1 ð Þ, the type of partial charges should be motivated by the cancellation of ground state inter-region electrostatic contributions introduced in the second term of Equation 6. Using that same set of charges for E EE QM 1 ð Þ has been argued to compensate overpolarization effects in ground state ONIOM calculations [51] ; however, this invites a low-resolution description of the environmental electrostatic interactions due to the stringent constraints on the computational cost of the QM 0 level of theory. In fact Restricted ElectroStatic Potential (RESP) charges originating from a high level of theory population analysis should by definition provide a point charge potential closest to that of the molecular electron density. fromage allows the user to choose between RESP and Mulliken charges for all of its electrostatic embedding formulations.

| Ewald embedding
Whatever the choice, a serious problem arises with the use of cluster models to represent the condensed phase. Indeed environmental effects on the excitation stemming from molecules beyond the immediate vicinity of Region 1 have been completely omitted. This is acceptable for short-range contributions to the interaction energy since they decay very quickly with distance and are therefore dominated by nearest neighbor contributions. However, the Coulomb interaction decays as an inverse law and therefore has significant long-range contributions to the excitation, compounded by the longrange order symmetry in the case of molecular crystals.
Simply increasing the size of Region 2 is not an adequate solution as the Madelung sum in crystals is known to be slowly and conditionally convergent, thus requiring an impractically large cluster with difficult to prove accuracy. The Ewald sum is an equation designed to closely approximate the Madelung potential with fast convergence by splitting the electrostatic potential terms into real and reciprocal space via the use of judicious Fourier transforms. The resulting potential is as follows: [52,53,67,68] V Ewals r ð Þ = X where L and G are the real and reciprocal space lattice vectors, q s are the charges of each lattice site s of the unit cell at positions R s , γ is the Ewald constant, and v c is the volume of the unit cell. A singularity arises at r = R i , which is corrected by replacing the L = 0 and s = i case of the first term of Equation (7) with − 2γq i ffiffi π p .
Klintenberg, Derenzo, and Weber have developed a method to fit an array of approximately 10,000 point charges to match the Ewald potential in a cluster region of a crystal. [17,18,54] The charges are distributed in a supercell and fitted with a system of linear equations to match the evaluated Ewald potential at the atomic sites of the central region.

| Self-consistent embedding
A further problem to address regarding electrostatic embedding in ONIOM is that of the mutual polarization between regions. Indeed if the embedding charge values are fixed, they are unable to react as the real charge density would to changes in the electronic structure of Region 1. Most notably, real excitations are liable to provoking an equilibration process with the environment. In order to reflect such processes, fromage combines the self-consistent Ewald embedding method by Wilbraham et al. [55] with the ONIOM QM:QM 0 formalism.
The method is as follows: 1. Region 1 is embedded in an Ewald point charge array.
2. An excited state population analysis is performed on Region 1.  Table 2.

| Charge value distribution
OEC, OEEC, and SC-OEEC calculations all require the distribution of point charges from single molecule population analyses to aggregate geometries or unit cell structures. This is a routine operation for preparing force field calculations, where sets of atomic types are associated to a potential or another with predetermined charge values. [56,57] The atomic type is usually dependent on its neighboring bonded atoms and the bond types. For QM:QM 0 calculations, it would be more attractive to use a broader definition of the atomic type which distinguishes between same function atoms at nonequivalent parts of the molecule.
To this end, fromage implements a connectivity detection tool which reads population analysis information from a single molecule or unit cell calculation and redistributes it onto any other finite or periodic ensemble of same molecules. The procedure builds a bond order matrix B where B ij is the shortest path connecting atoms i to j in number of bonds. The construction of the matrix is represented in Figure 8 and is performed as follows: will produce a file out char which contains the coordinates of clust.xyz with an added column stating the partial charge of each atom.

| Geometry optimization
Equation 6 provides the energy of the central molecule embedded in its environment. It can readily be extended to its energy gradients, thus allowing for the optimization of geometries. fromage uses the BFGS algorithm as implemented in scipy [16] to locate ground and excited state minima.
Other critical regions of the PES also play a role in determining the emissive behaviour of the material. When accessible, conical intersections can act as ultrafast funnels producing a nonradiative decay mechanism. It is therefore desirable to locate minimal energy conical intersections (MECIs) and assess their energetics and accessibility. In performing this search, the penalty function method [58] is employed to avoid the need for nonadiabatic coupling vectors which are not necessarily available depending on the quantum mechanical method employed. The quantity to be minimized is then: where E 1− 0 is the average of the S 1 and S 0 energies ΔE the energy gap, σ is a Lagrangian multiplier, and αis a parameter much smaller than ΔE. While multireference methods should be used in order to obtain a formal topology of the conical intersection, MECIs obtained with single-reference methods have proven to also yield informative results at a more modest computational cost. [59][60][61][62] 3.3.6 | Example of use  1-one (DMAP). [26,27] Both HC and DMAP can experience excited-state intramolecular proton transfer, splitting the excited state into two potential decay pathways. It was recently computationally shown that in HC, both the enol and the keto nonradiative decay pathways were rendered energetically inaccessible by a combination of molecular and crystalline factors. [19,63] OEEC was employed to investigate the critical points in the PES of DMAP. Geometry optimizations were carried out to locate the ground and excited state minima, as well as its MECI geometries. The protocols. [50,[64][65][66]

| CONCLUSIONS
We have detailed the principal capabilities of the Python library fromage, which aims to facilitate the computational investigation of the excited states of molecular aggregates. They include geometrical and exciton analysis as well as QM:QM 0 geometry optimization tools, which have already been successfully employed in three past publications. [12][13][14] The features were tested on a diverse array of molecules, in order to challenge their robustness. They are implemented with enough flexibility that Python literate researchers can employ fromage scripts as part of a larger workflow with little added effort. By virtue of being an open source, unit tested and documented piece of software, fromage represents an addition to the fast expanding pool of sustainable chemical software libraries. We hope that by enabling researchers to use the framework and manipulate the source code, the field of aggregate photochemistry will continue to mature towards modern reproducible workflows. Rachel Crespo-Otero https://orcid.org/0000-0002-8725-5350