# A High-Throughput Computational Protocol for Tuning Molecular Properties: Application to ESIPT Chromophores

Isabella C D Merritt, Frédéric Castet

► To cite this version:

Isabella C D Merritt, Frédéric Castet. A High-Throughput Computational Protocol for Tuning Molecular Properties: Application to ESIPT Chromophores. Journal of Chemical Information and Modeling, 2025, 65 (10), pp.5035-5049. 10.1021/acs.jcim.5c00692. hal-05114398

HAL Id: hal-05114398

https://hal.science/hal-05114398v1

Submitted on 16 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A High-Throughput Computational Protocol for Tuning Molecular Properties: Application to ESIPT Chromophores

Isabella C. D. Merritt* and Frédéric Castet

Univ. Bordeaux, CNRS, Bordeaux INP, Institut des Sciences Moléculaires (ISM), F-33405, Talence, France

E-mail: isabella.merritt@u-bordeaux.fr

## Abstract

Over the last decade, improvements in computing power and theoretical approaches have enabled high-throughput computational investigations of systems. In this work we present the development of a simple automated computational protocol for the study of molecular substitutions to known molecules, which minimises human error and effort while capitalising on existing calculations to optimise computational cost. We demonstrate the use of our protocol on three test-cases of known chromophores undergoing intramolecular proton transfer: (1) a focused study of 12 molecular derivatives, where the protocol is run locally on a standard laptop, (2) a larger study of 169 derivatives which allows for investigation of trends influencing ground and excited state reactivity, and (3) a realistic study of how our protocol could be used to investigate a large (over 700) set of derivative molecules and select candidates fulfilling required criteria for a particular given application. Our protocol, available online, is designed to be user-friendly and lightweight, and enables efficient and straightforward screening of hundreds of molecular derivatives.

## 1. INTRODUCTION

Over the last 30 years, improvements in both computing power and in theoretical approaches have led to an incredible growth in the capabilities of computational chemistry. Theoretical modelling, often using the game-changing KS-DFT method first introduced by Nobel laureates Kohn and Sham in 1965,¹ is now routinely used in chemical research, both in collaboration with experiment and as a driving force for new discoveries. With the now widespread availability of computational resources, made possible by the rapid scaling of computational power over the last decades,² it is now possible to carry out high-throughput theoretical investigations of systems. This can allow for explorations over large ranges of chemical space, which experimentally would be challenging or even impossible due to time and/or cost constraints.

High-throughput automated computational systems have been developed in many research fields. In materials research, the HTES package introduced by Nepal et al.³ can be used to automate *ab initio* calculations based on existing material databases; while Fortunato et al.⁴ developed and used a high-throughput approach to screen candidates for novel magnetocaloric materials. In reaction kinetics, the AutoTST method of Bhoorasingh et al.⁵ can be used to automatically calculate large numbers of reaction rates by locating transition states in both uni and bi-molecular reactions. In medical research, high-throughput calculations are routinely used to predict protein stability, folding, and binding;⁶ as well as for prediction of properties of new potential drug candidate molecules.⁷ These examples represent only a small fraction of the uses of high-throughput automated calculations in scientific

research, especially as many such systems are developed and used locally within research environments and not made publicly available.

In molecular chemistry, especially in the field of light-driven processes, a common approach both experimentally and theoretically is to modulate the properties of existing molecules through substitution of functional groups, very often on aromatic rings.⁸⁻¹⁵ These substitutions, often chosen to either donate or withdraw electron density to/from the rest of the molecule, can be used to adjust the properties of molecules undergoing a desired reaction/processes, making the substituted derivatives more suited for the chosen application: for instance, adding substitutions to a chromophore in order to red-shift UV absorption wavelengths into the visible part of the spectrum. Such chemical approaches have been used to design improved molecules and materials for many applications, such as OLEDs,¹⁶ organic solar cells,¹⁷ drug-design and discovery,¹⁸ and more. The advantage of computational investigations over experiment is that expensive and time-consuming synthesis can be avoided, making it more economically and ecologically friendly to investigate large numbers of molecular derivatives (although CO₂ production due to calculation cannot be ignored). Theoretical investigations can be best used as a tool to select candidate derivatives from a large number of potential molecules, with experimental testing thus being carried out on a smaller more targeted group of molecules.

One field of chemistry of high interest in recent years is the study of light-responsive molecules, in particular molecules which undergo photochemical reactions after the absorption
of a photon. One common type of chemical reaction which
may occur after light-absorption is the transfer of a proton
from a donor to an acceptor within the same molecule, typi-
cally along existing hydrogen bonds in the system : known as
excited state intramolecular proton transfer (ESIPT).¹⁹ ES-
IPT reactions are of interest for many different applications,
such as in chemical and light sensors, and as broadband-
light emitters thanks to their possible dual fluorescence.²⁰⁻²³
The use of a high-throughput approach to investigate the im-
pact of a large number of substitution on different molecules
which undergo the ESIPT reaction is of interest, in order to
select best-candidate molecules for applications as well as
determine trends of molecular properties.

In this work, we have implemented a simple automated
computational protocol to easily calculate properties of sub-
stituted molecules. Our method allows for a high-throughput
investigation of hundreds of variants of known molecules
and/or reactions, while minimising computational cost. This
protocol benefits from both the advantages of an automated
input generator, vastly reducing the amount of human time
required to build molecular input files and removing the po-
tential for human error, while also minimising computational
time required by capitalising on existing calculations of the
unsubstituted molecule, resulting in simpler geometry opti-
misation convergence. This protocol is made available online
as lightweight and easily modifiable python scripts designed
to create the required inputs for computations using either
the Gaussian16²⁴ or Orca²⁵ quantum chemistry program (in
this work we have used Gaussian16). The scripts are written
in a nature which would also allow easy adaptation to any
standard quantum chemistry software package.

<figure><img src="image_2.png" /><figcaption>Figure 1: The three molecules investigated using the computational protocol: (a) 2-(2'-hydroxyphenyl) benzotriazole (BZT) (b) 2,11-dihydro-1H-cyclopent[de]indeno[1,2-b]quinolin-7-amine (CPIQ), (c) salicylidene aniline (SA). Labelled atoms correspond to positions where substitutions were tested, guided by experimental knowledge of synthesizable derivatives of such compounds.</figcaption></figure>

Here, we demonstrate the advantages of this protocol by
investigating the effects of substitutions on three molecules
known to undergo excited-state intramolecular proton trans-
fer (ESIPT), shown in Figure 1, in order to demonstrate
our protocol for different use-cases. The first molecule se-
lected is 2-(2'-hydroxyphenyl) benzotriazole (BZT), the par-
ent species of a wide range of molecules commonly used as
UV-absorbers in industrial applications to prevent UV-driven

degradation.²⁶,²⁷ We use this molecule as an example of using
our protocol to carry out a simple investigation of the impact
of a small number of different substituents on the absorp-
tion wavelength. The second molecule investigated (Figure
1(b)), which we name CPIQ in line with previous literature,
is a recently designed and synthesised⁸ rigid chromophore,
which undergoes ultrafast excited state intramolecular proton
transfer (ESIPT) along the N–H···N intramolecular hydro-
gen bond upon excitation in the near-UV. This molecule is
taken as an example of the use of our protocol to investigate a
large number of different substitutions at a limited number of
given positions, where experimental substitution is expected
to be viable. The final molecule is salicylidene aniline (SA),
a photochrome which upon excitation also undergoes ES-
IPT, converting from the enol to the keto form. This ESIPT
is then followed by *cis-trans* photo-isomerisation about the
C–C phenolic bond in the keto form,²⁸⁻³⁰ generating a short-
lived *trans*-keto photochromic form with measured lifetimes
on the order of microseconds.³¹,³² This final molecule is
taken as an example of the use of our protocol to evaluate
large volumes of substitutions at many possible substitution
positions, the synthesis of such molecules being relatively
straightforward and allowing for a large number of synthe-
sizable derivatives.³³

We investigate using our protocol the impact of chemical
substitution on several properties relating to possible appli-
cations of these photoactive molecules. These include; exci-
tation energies and oscillator strengths (i.e. predicted UV-vis
absorption peaks), tautomeric shifts of these absorptions,
relative stability of different tautomers, barriers to proton
transfer on the ground and excited state and/or *trans-cis* iso-
merisation, and relative dipole moments between different
forms of the photoswitches.

## 2. METHODOLOGY

### 2.1 High-Throughput Protocol

An overview of the protocol used for high throughout com-
putational screening is outlined in Figure 2, summarised as
a fundamental series of steps. The first three steps are those
carried out by the user and involve the definition and pre-
investigation of the system, before the succeeding steps are
then carried out in an automated fashion by the protocol, gen-
erating input files which must then be run using the chosen
quantum chemistry program.

#### 2.1.1 Initial system selection and investigation

The first step when using the protocol is to select the sys-
tem to be studied. As an example, we shall take the second
molecule defined in Figure 1, CPIQ. We shall suppose in this
case we are interested in the relative stability and ground state
intramolecular proton transfer (GSIPT) conversion between
the two tautomeric forms of this molecule. The second step
is then to identify and select key geometries relevant to this
GSIPT reaction. We can identify three key geometries to
be investigated depending on the location of the transferring
proton: the two tautomers and the GSIPT transition state con-
necting the two. We must then optimise these three geome-
tries using our chosen level of theory, to correctly describe
the GSIPT in the non-substituted molecule. These optimised
geometries, in standard xyz format, must then be collected
<figure><img src="image_3.png" /><figcaption>Figure 2: Protocol for high-throughput calculation of properties of substituted molecules.</figcaption></figure>

into a single location - importantly, the atomic ordering in the XYZ files must be consistent across all geometries.

The third step is to select the desired substitution positions and which functional groups we wish to investigate. Provided as part of the protocol code is a library of common electron donating and electron withdrawing substituents in the required format (ESI, Section S2), but any additional substituent can be added to this library with minimal effort. These substituents are defined in standard xyz format, and are optimised at the DFT level (see Section 2.2). Importantly, they are defined with an additional H atom in place of the bond where they would connect to the molecule - e.g., for the –OH functional group, an xyz file for H₂O is given in the substituent library. The comment line of the xyz file must follow a particular format, containing the identification of 3 key atoms in the substituent : the H atom which should be removed when joining this substituent to a molecule, the connecting atom from which a bond should be formed to the molecule, and the next adjacent atom to this “connecting” atom to allow for testing of substituent rotation. The full list of all functional groups available in the “substituent library” can be found in the ESI.

It is from these first three steps that the input file(s) to the protocol should be constructed. The input file contains all necessary information for the program to generate the inputs of all substitution combinations requested - an annotated example of such an input is given in the ESI. It follows a modular structure, starting with a [GENERAL] section defining the system name, computational program and paths to necessary folders, as well some additional options for the program. Following this is a section defining parameters relating to the computational method (e.g. DFT functional) ([PROG_PARAMS]). Finally, a block is given for each key geometry to be investigated where the name, nature (minimum, transition state), and electronic state are given, followed by a block for each substitution position within which the desired functional groups and substitution position are defined. We note that a separate input file is required for each substitution pattern to be investigated, such that if one wishes

to investigate substitutions at position A or position B on a molecule an individual input file will be required for each position. Multiple substitution blocks in the input file instead generate multiply-substituted molecules, for instance if one was interested in the combination of different substitutions on positions A and B.

### 2.1.2 Automated substitution connection protocol

Taking the input file given by the user, along with the optimised geometries of the unsubstituted molecule, the program automatically carries out the next steps given in Figure 2. All possible combinations of substitutions (when multiple-functionalised molecules are considered) are defined and constructed. For each particular key geometry, in order to build each substituted molecule, the steps outlined in Figure 3 are carried out. We note that for each substituted molecule, we only consider the most stable conformer and do not account for alternative conformers due to the addition of a substituent.

In step 1, the program starts with the optimised geometry of the unfunctionalised molecule, henceforth referred to as the “core”. A unit vector ($V_1$) is defined along the bond-coordinate between the substitution position $C_C$ and the atom $C_R$ (typically hydrogen) to be replaced by the functional group (as defined in the input file). The functional group geometry, henceforth referred to as the “fragment”, is then also read, and a second unit vector ($V_2$) is defined along the bond coordinate for the connection atom $F_C$ and the atom to be removed $F_R$ (as defined in the substitution library file). Then, in step 2, the fragment and core are translated in space such that the connecting atom of the core ($C_C$) and the atom to remove of the fragment ($F_R$) are superposed at the origin. A transformation matrix to rotate vector $V_2$ onto $V_1$ is then defined and used ($T_{rot}$), aligning these two vectors, and the two hydrogens to be removed are deleted (from Step 2 to 3). Through this approach, an approximation to the substituted geometry is generated, exploiting the pre-optimised core geometry and fragment geometries to be as close as possible to the true minimum or transition state of the substituted form.

However, there is still one key coordinate associated with
Step 1

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Step 2

<figure><img src="image_6.png" /></figure>

Step 3

<figure><img src="image_7.png" /></figure>

Figure 3: Procedure for automatically joining a substituent (fragment) to a particular molecular core. Vectors $V_1$ and $V_2$ are defined along the bond between the atom to be deleted and the "connecting" atom for both the core and the fragment geometries. Vector $V_3$ is defined as that between the "connecting" atom and an adjacent atom for the fragment.

the connection of the fragment to the core: the relative rotation of the fragment about the connecting bond (i.e. the dihedral angle between the core and the fragment). In order to start with a reasonable guess to this relative orientation, and avoid possible clashes and overlaps of atoms on the core and fragment, a scan of energies along this rotation is carried out at the Hartree-Fock level with a minimal basis set using the pyscf package,³⁴ and the lowest energy orientation is selected as the final geometry. This rotation of the fragment relative to the core is defined using the Rodrigues Formula³⁵ for the rotation of a vector about an axis, rotating the vector $V_3$ (Figure 3, Step 3) between the fragment connecting $F_C$ and adjacent $F_A$ atoms about the previously defined vector $V_1(\equiv V_2)$.

This connecting procedure as defined in Figure 3 is carried out in succession for each substitution added to the core geometry. Finally, once all substitutions have been added, input files for the full optimisation are written starting from this "initial best guess" generated geometry, and these can be submitted and run individually in parallel.

### 2.1.3 Installing and using the protocol

The required scripts along with a pre-optimised substituent library and example input files are publicly available at github.com/issymerritt/generate_subs. Instructions for installation and usage are given in the README.md file as

well as in the SI, including the installation of required packages ( numpy, pyscf, periodictable).

## 2.2 Computational Details

All calculations in this work were carried out using the Gaussian16 program, using the ultrafine integration option (pruned (99,590) grid) and standard optimisation convergence criteria (RMS force threshold of $3 \times 10^{-4}$).²⁴ Geometry optimizations were systematically followed by vibrational frequency calculations, to ensure the correct number of imaginary frequencies (0 for minima, 1 for transition states). For the BZT molecule, geometry optimisations were carried out at the DFT level and vertical excitations at the TD-DFT level, using the PBE1PBE³⁶ functional and 6-31+G(d,p) basis set,³⁷⁻⁴⁰ in acetonitrile using the polarizable continuum model (integral equation formalism variant) for solvent effects.⁴¹,⁴² The same basis set, functional, and solvent treatment were used for the CPIQ photochrome, with geometry optimisations carried out at both the DFT and TD-DFT level. Ground state Hirshfeld charge analysis was carried out on the optimised ground state, and excited state charge analysis on the optimised $S_1$ state, for each substituted CPIQ initial tautomer. The same basis set and functional were used for the SA molecule, but the solvent was instead defined manually with values of relative permittivity $\epsilon = 5.5$ and square of index of refraction $n^2 = 1.56$, values matching experimental measurements for chitosan.⁴³,⁴⁴ All orbital visualisations and molecular images were produced using the VMD software package.⁴⁵ Initial geometries used in the "core library" for each molecule were optimised at the same level as given above, while the "substitution library" was built by optimising each functional group, with a hydrogen in place of the bond which connects the substitution to the main molecule, at the B3LYP/6-31G(d,p) level⁴⁶ in gas phase.

# 3. RESULTS AND DISCUSSION

## 3.1 Hydroxyphenyl benzotriazole (BZT)

We first demonstrate the advantages of our protocol for the simple investigation of excitation energy of the 2-(2'-hydroxyphenyl) benzotriazole (BZT) molecule (Figure 1(a)). We have investigated the impact of substitution at different positions on the hydroxyphenyl ring; ortho ($o$), meta ($m$) and para ($p$) to the hydroxy group which takes part in the ES-IPT reaction after excitation. Four substitutions are chosen: -NO₂ and -Cl as examples of strong and weak electron withdrawing groups (EWGs) respectively, and -OH and -CH₃ as examples of strong and weak electron donating groups (EDGs) respectively. We highlight that such an investigation of the BZT molecule is not totally novel, with similar results having been previously outlined in literature for related derivatives.⁴⁷⁻⁵⁰ We do not aim in this section to discover new chemical or physical properties of BZT, instead to demonstrate how a typical small-scale investigation of such a molecule could be carried out thanks to our protocol, while also providing a clear description of trends which can be hard to find in literature.

Using our automated system, we have calculated the vertical excitation energies (VEE) of the 12 possible combinations of substitution and position (Figure 4), re-optimising the ground state Franck-Condon geometry for each compound
<figure><img src="image_8.png" /><figcaption>Figure 4: Vertical excitation energies (VEE) for the first (bright) excited state of variants of 2-(2'-hydroxyphenyl) benzotriazole, substituted on the hydroxyphenyl ring *ortho*, *meta* and *para* to the -OH group. The dotted line corresponds to the VEE of unsubstituted BZT.</figcaption></figure>

before calculating VEE. In all cases apart from -OH in the *ortho* position, the first excited state corresponds to the first optically allowed transition of type $π → π^*$. For *o*-OH, the first excited state is dark ($f = 0.09$) and comes from excitation from a stabilised localised orbital on the hydroxyphenyl ring (see ESI, Section 3). We take instead the second excited state for this molecule, being the first bright state and corresponding to the same transition as the $S_1$ in the other molecules.

The VEE as calculated for these molecules correlates closely, as expected, with the HOMO-LUMO (HOMO-1 for *o*-OH) energy gap, with a Pearson correlation coefficient squared of $R^2 = 0.95$ (Pearson's coefficient is henceforth referred to as R). However, if we exclude the -NO₂ substitutions, this correlation increases to $R^2 = 0.98$ (see ESI, Section 3), matching with the fact that the first excited state for all substitutions is well described by a single orbital transition, while for all three -NO₂ substitutions there is an important (weight > 0.1) contribution from one other excitation. We can thus look at the relative HOMO (HOMO-1 for *o*-OH) and LUMO orbital energies and forms to interpret the calculated VEE values.

Figure 4 shows the excitation energies of the different substitution patterns, compared to the unsubstituted molecule (dashed line, VEE = 3.75 eV), while Figure 5 shows the orbitals associated with the electronic transition and their relative energies. Three molecules show significant red-shifted VEE below 3.5 eV - *o*-NO₂, *m*-NO₂ and *p*-OH, while only *m*-OH shows significantly blue-shifted VEE at 4.00 eV. The weak EDG and EWG (-CH₃ and -Cl) groups have only a minor impact on the excitation energy as expected.

With the exception of -OH, we can see that EWGs stabilise both the HOMO and LUMO, while the EDG -CH₃ destabilises them. The differences in VEE arise from which orbital (HOMO or LUMO) is more (de)stabilised relative to the other. The weak EDG (-CH₃) always leads to red-shifted excitation compared to the unsubstituted molecule, as the HOMO is more destabilised than the LUMO, leading to a HOMO-LUMO gap smaller than that of the unsubstituted (4.50 eV). Interestingly, para substitution leads to red-shifted VEEs no matter the type of substitution, with increasing red-shift with higher electron donating character. In general we see also that no matter the substitution position,

<figure><img src="image_9.png" /><figcaption>Figure 5: Orbitals and energies associated with the first (bright) electronic transitions of variants of 2-(2'-hydroxyphenyl) benzotriazole, substituted on the hydroxyphenyl ring *ortho*, *meta* and *para* to the -OH group. Dotted lines correspond to the orbital energies of the unsubstituted BZT, and the numbers in grey give the HOMO-LUMO energy gaps.</figcaption></figure>
VEE(−CH₃) < VEE(−Cl), indicating that red-shifting of the VEE of these molecules is more favoured by weak electron donating character than weak electron withdrawing character.

The outliers which are o-NO₂, o-OH, m-NO₂ and p-OH are explained through resonance effects leading to significant orbital stabilisation. In both o-NO₂ and m-NO₂, the LUMO is highly stabilised thanks to the impact of π* orbitals coming from the −NO₂ group mixing with the rest of the π* framework. In o-OH, the presence of the two adjacent −OH groups leads to stabilisation of the HOMO-1 orbital, with the π contribution from the o-OH group leading to mixing with a different π orbital coming from the hydroxyphenyl ring, with little-to-no π contribution on the original −OH group taking part in the intramolecular hydrogen bond.

Finally for the BZT molecule we can take a look at the relative lengths of the O–H covalent bond and H···N hydrogen bond, as a measure of the “proton transfer” character in the ground state minimum, with the view that if the hydrogen is situated relatively closer to the N atom the proton transfer is likely to be more favourable. Figure 6 shows the value for each of the 12 molecules of

$$H_{\text{trans}} = \frac{OH_{\text{dist}}}{OH_{\text{dist}} + NH_{\text{dist}}} \quad (1)$$

$H_{\text{trans}}$ thus being a measure of how close the proton is located relative to the O and N atoms.

<figure><img src="image_10.png" /><figcaption>Figure 6: $H_{\text{trans}}$, a measure of the relative position of the proton with reference to the O and N atoms, for each of the 12 substituted BZT molecules, compared to the unsubstituted molecule (dashed line).</figcaption></figure>

We can see that for both *ortho* and *para* substitution, there is a clear trend for a smaller value of $H_{\text{trans}}$ with increasing EDG character - meaning the proton is located relatively closer to the O atom with a relatively longer hydrogen bond to the N atom, potentially disfavouring proton transfer. Comparing *ortho* and *para* substitution, $H_{\text{trans}}$ is systematically lower for *para* substitution than *ortho*, suggesting *ortho* substitution disfavours the OH bond and favours the H···N hydrogen bond compared to *para*. Finally, *meta* substitution appears to have little effect on the relative position of the H atom, and indeed on closer inspection has little impact on both the OH and H···N bond lengths.

For this first molecule, we have chosen to use our protocol

to demonstrate a small and focused study of a small number of substitutions on a single key geometry (the Frank-Condon point), aiming to simply describe VEE energies and geometrical differences between the molecules. With a small number of results, we have been able to enter into details of why certain VEE values and trends appear. Importantly, as only 12 different geometries needed to be built by our protocol, we were able to run the geometry building locally on a laptop, taking only 3 minutes. Indeed, for such a small investigation the time taken to run the protocol is roughly the same as how long it takes to set up the required files and write the input for our protocol, which we also estimate at 3 minutes, taking advantage of template input files and scripts to write xyz files. For such a simple case of substitutions on a single geometry this is not expected to be a large time-saver - optimistically assuming construction by hand of an input file takes ca. 1 minute per molecule for an experienced user, this would mean generation by hand of all 12 input files would take 12 minutes, meaning our protocol only saves 6 minutes in this case. However, the main advantage in this example of a small and focused study is the minimisation of possible human errors - as the geometry is built automatically, it is impossible to accidentally delete atoms or place new atoms incorrectly.

## 3.2 CPIQ

Next, we demonstrate the use of our protocol for the investigation of both ground-state and excited-state reactivity in a photochrome known to undergo an ESIPT reaction, the 2,11-dihydro-1H-cyclopent[de]indeno[1,2-b]quinolin-7-amine (CPIQ) molecule, first designed and explored in Ref. 8. After light absorption, the ESIPT reaction of the molecule occurs over a barrier of 33.7 kJ/mol (0.35 eV), this kinetic barrier to ESIPT resulting in dual-emission from both the normal (N) and the tautomeric (T) forms.

For this molecule, we will investigate properties relating to both the ground state reactivity (GSIPT) as well as the excited state using our automated approach. Figure 7 shows an overview of the six key geometries to be investigated as identified in Ref. 8, which describe five key properties of this photochrome, which can be divided into two groups. The first group consists of three properties relating to the ground state reactivity, and are described by DFT optimised geometries (Figure 7, left):

1. Relative stability of initial (N) and tautomeric (T) forms.

2. Energy barrier to the (N → T) GSIPT reaction.

3. Energy barrier to the (T → N) reverse GSIPT reaction.

The second consists of two properties relating to the excited-state reactivity, requiring TD-DFT, being:

1. Stokes shift (N absorption, T* emission).

2. Energy barrier to the (N* → T*) ESIPT reaction.

We have used our protocol to investigate the impact of substituting a range of EDG and EWG (SI, Section 4) at three positions on the CPIQ molecule (Figure 1), identified as experimentally feasible sites for substitution. We have also investigated the impact of dual-substitution of type push-pull,
<figure><img src="image_11.png" /><figcaption>Figure 7: Overview of the key geometries describing the excited-state (ESIPT) and ground state (GSIPT) intramolecular proton transfer reactions, as well as the absorption and emission leading to the large stokes shift, in the unsubstituted CPIQ photochrome.</figcaption></figure>

with one EWG and one EDG group on opposing sides of the molecule.

We shall begin by discussing the first group of properties, those relating to the forward and backwards GSIPT reactions determined from ground state DFT calculations. Figure 8 shows the spread of these three properties with different types of substitution. We can see that the trends for the relative N-T stability and the energy barrier to ($N \to T$) GSIPT reaction are almost identical, which is expected since the GSIPT transition state (TS) is very close in energy and similar in form to the ground state optimised tautomer. The ($T \to N$) reverse GSIPT barrier shows different trends to the other two.

Starting with EWG substitution, we can see that all three properties can be notably tuned, with the most impact coming from substitutions on positions 3 and 13, in line-with the fact that these substitution positions are closer to either the proton donor or acceptor nitrogens and might thus be expected to have more influence on this reaction compared to position 24. Substituting an EWG on position 3, *ortho* to the –NH₂ group, has a tendency to lower the GSIPT barrier and the relative N-T stability, by up to 20 and 28 kJ/mol respectively; i.e. substitution here favours more the formation of the tautomer. This makes sense, as an EWG at this position will lead to reduced electron-density on the –NH₂ group, increasing the acidity of the proton. The impact is much less clear on the barrier to reverse GSIPT, with some substitutions increasing and others decreasing the barrier. On the other hand, substituting an EWG at position 13 has a varied impact on the relative N-T stability and GSIPT barrier, both increasing and decreasing these properties by up to 15 kJ/mol depending on the specific substitution. While an EWG at position 13 might be expected to have the inverse effect to position 3, by reduc-

ing the basicity of the acceptor-N, the real impact is harder to predict. This is likely due to the in-space positioning of the EWG, which depending on the EWG size and shape can also interact with the proton on the –NH₂ group. The impact on the reverse GSIPT barrier is also varied, lowering or raising it to as high as 16 kJ/mol and as low as 5 kJ/mol. Finally, substituting an EWG at position 24 has a much smaller impact on the “ground state” properties, only modulating them by a maximum of 7 kJ/mol, possibly due to the further distance from the acceptor nitrogen and the greater physical separation from the –NH₂ group. However, at this position we see the expected impact of placing an EWG near to the proton-acceptor nitrogen, with higher N-T relative stability and barrier to ESIPT, and lower barrier to GSIPT.

Next, looking at EDG substitution, we see that with the exception of a small group of substitutions at position 13, there is a much smaller impact on these “ground-state” properties no matter which position is chosen. While substitution at position 3 tends to again decrease relative N-T stability, lower the GSIPT energy barrier, and increase the reverse GSIPT barrier; at position 24 the opposite effect is seen. For these two positions, variations of only ca. 3-4 kJ/mol are seen, much less significant than with EWG substitution. At position 13, some more variation is produced, but still less than with EWGs, with the exception of a small group of three molecules with notably increased relative N-T stability and GSIPT barriers. These three molecules correspond to substitutions of acylamido groups (–NHCOR, $R=CH_3,CH_2CH_3,C_6H_6$), and show this particular behaviour due to the positioning of the acylamido group, which is in plane with the $\pi$-framework in the N tautomer, but is pushed slightly out of plane in the GSIPT TS geometry (CCNC=24°)
<figure><img src="image_12.png" /><figcaption>Figure 8: Impact of substitution of electron withdrawing (EWG) and donating (EDG) groups at different positions of the CPIQ molecule on the “ground-state” properties: (a) relative stability of the initial (N) isomer compared to the tautomer generated by proton transfer (T), (b) energy barrier to the GSIPT reaction (N → T), and (c) energy barrier to the reverse GSIPT reaction (T → N). Orange lines correspond to the values of each property for the unsubstituted molecule. Figure (d) shows the positions at which Hirshfeld charges are best predictors of these three thermal properties, with highest R² correlations (see Table 1).</figcaption></figure>

and even further out of plane in the T tautomer (CCNC=53°) (see SI, Section 4), due to a clash between the transferred proton and the proton on the acylamido group. This leads to destabilisation of these two geometries, and thus the notable stabilisation of the N tautomer in comparison, leading to the highest N-T stabilities and GSIPT barriers for any substitution type.

Finally, we can look at the impact of “push-pull” substitution, referring to the dual addition of an EWG at position 3(13) and an EDG at position 13(3), labelled as 3-13 (13-3). We can see that push-pull substitution has the potential to vary significantly all three of these properties, with the opposing positions of EWG vs EDG leading to inverted trends. EWG(3)-EDG(13) leads generally to lower relative N-T stability and lower GSIPT barrier similarly to the effect of a single EWG on position 3, with the additional EDG not seeming to have a significant effect. However, we see a much larger ability to increase the reverse GSIPT barrier using push-pull substitution pattern compared to single EWG (or EDG) substitution, with barriers as high as 20 kJ/mol, corresponding to a harmonic transition state theory based unimolecular half-life⁵¹,⁵² of 350 ps for the photochemically generated tautomeric ground state, significantly longer than the predicted half-life of the unsubstituted CPIQ of 4 ps. The opposite substitution pattern, EWG(13)-EDG(3) has a less clear impact on the relative N-T stability and GSIPT barrier, but mostly tends to increase these two properties. The impact on the reverse GSIPT barrier is varied, with some compounds having a lower barrier and thus a predicted faster recovery of the initial tautomer, while others show a lower barrier. We can also note that, for the push-pull substitution, we see a clear range of values between extremes for all three properties - suggesting that, by carefully selecting different push-pull substitutions, we could tune the properties as desired to specific values.

To complete the investigation of these three properties relating to ground state reactivity, we have investigated if the results of these properties could be predicted without calculation of the tautomer or transition state, instead looking only at the ground state initial form of the molecule. In a previous theoretical study, Luchini et al.⁵³ found strong correlation between Hirshfeld charges⁵⁴ on adjacent atoms, and experimental Hammett parameters of substituents, used to predict their influence on reactivity of aromatic systems.⁵⁵ We decided therefore to investigate whether correlations existed between the charges at different atomic sites for substituted CPIQ, and the three parameters in Figure 8, taking each substitution position individually. The correlations between the properties and calculated charges on every atom are given in the SI (Section 4). Interestingly, we find that while Hirshfeld charges can be used to predict the values of all three properties, the atom on which the charge is the best predictor varies depending on the substitution position and type. Figure 8(d) shows the atoms which are the best predictors (the highest R² for all three properties), for each substitution position and type, with Table 1 showing their respective R² values.

No atomic charges are found to be moderate predictors (R² > 0.7) of the reactivity for EDG 3 substitution, consistent with that fact that this position has only weak impact on the
<table><caption>Table 1: Selected R² correlations between the values of the relative N-T stability (ΔE<sub>N</sub> – ΔE<sub>T</sub>), GSIPT energy barrier ΔE<sub>N→T</sub> and reverse GSIPT energy barrier ΔE<sub>T→N</sub>, with the Hirshfeld charges calculated at different atomic positions on the optimised Frank-Condon geometry (N) of substituted molecules. Charge position refers to the atom numbering as given in Figure 8(d).</caption><thead><tr><th rowspan="2">Type</th><th rowspan="2">Charge Position</th><th colspan="3">R<sup>2</sup> correlations</th></tr><tr><th>ΔE<sub>N</sub> – ΔE<sub>T</sub></th><th>ΔE<sub>N→T</sub></th><th>ΔE<sub>T→N</sub></th></tr></thead><tbody><tr><td>EWG 3</td><td>19</td><td>0.85</td><td>0.77</td><td>0.70</td></tr><tr><td>EWG 13</td><td>22</td><td>0.88</td><td>0.84</td><td>0.85</td></tr><tr><td rowspan="3">EWG 24</td><td>10</td><td>0.93</td><td>0.94</td><td>0.83</td></tr><tr><td>31</td><td>0.97</td><td>0.96</td><td>0.86</td></tr><tr><td>32</td><td>0.95</td><td>0.95</td><td>0.85</td></tr><tr><td rowspan="3">EDG 13</td><td>3</td><td>0.90</td><td>0.89</td><td>0.96</td></tr><tr><td>7</td><td>0.96</td><td>0.94</td><td>0.98</td></tr><tr><td>11</td><td>0.98</td><td>0.94</td><td>0.97</td></tr><tr><td rowspan="5">EDG 24</td><td>10</td><td>0.96</td><td>0.93</td><td>0.96</td></tr><tr><td>15</td><td>0.96</td><td>0.95</td><td>0.93</td></tr><tr><td>19</td><td>0.92</td><td>0.89</td><td>0.91</td></tr><tr><td>31</td><td>0.93</td><td>0.90</td><td>0.94</td></tr><tr><td>31</td><td>0.95</td><td>0.93</td><td>0.95</td></tr></tbody></table>

properties, while only one atom shows moderate prediction for EWG 3 substitution (0.7 < R² < 0.85). This difficulty to predict reactivity when substituting at atom 3, ortho to the –NH₂ group, is likely due to the same reason Hammett parameters cannot be established for ortho-substitution, with steric hindrance outweighing the electronic effect differently depending on the substituent.⁵⁶ The best predictor of the reactivity with EWG substitution at position 13 is the –NH₂ donor proton, which could be expected to predict well GSIPT reactions. Similarly, the charge at the acceptor nitrogen is also an excellent predictor (R² > 0.9) for EDG substitution at position 24. However, the other atoms which best predict reactivity are not those which might be expected. Unusually, the aliphatic carbon 10 and its attached hydrogens are excellent predictors of reactivity for substitution at position 24, for both EWG and EDG, with R² > 0.95 for all properties aside from reverse GSIPT for EWG substitution (where R² is still ca. 0.85). Carbon 19 is another excellent predictor for EDG 24 substitution, while EDG 13 reactivity is best predicted by atomic charges at aromatic carbons ortho and para to the –NH₂ group. Overall, we can see that while it would be possible to predict GS reactivity for substitution at positions 13 and 24 using Hirshfeld charges, investigation is necessary to determine at which atom these charges should be taken.

We can now move on to discussing the second group of properties, those relating to the excited-state reactivity and determined through TD-DFT energies and optimisations. Figure 9 shows the variation of these two properties, Stokes shift and ESIPT barrier, with different types of substitution. The Stokes shift here is determined as the difference between absorption of the initial tautomer (N → N*), and emission of the alternative tautomer (T* → T) formed after proton transfer. We only consider the Stokes shift if the calculated first excited state of N has an oscillator strength from TD-DFT greater than f = 0.1 - this excludes the substitution

of the –NO group at all positions, and the –NO₂ group at position 3.

<figure><img src="image_13.png" /><figcaption>Figure 9: Impact of substitution of electron withdrawing (EWG) and donating (EDG) groups at different positions of the CPIQ molecule on the “excited-state” properties: (a) difference in energy between N absorption and T* emission i.e., Stokes shift, (b) energy barrier to the ESIPT reaction.</figcaption></figure>

Again starting with EWG substitution, we can see that this tends to lower the Stokes shift by up to 0.4 eV, especially when placed at positions 3 and 24. The majority of substitutions however only modulate the Stokes shift by a small amount, less than 0.2 eV. If we look at the barrier to ESIPT, we see that certain EWG substitutions can greatly increase this barrier, but the majority have an impact of less than 12 kJ/mol. Substitution on position 3 and 13 has no systematic impact, while position 24 always increases the barrier. The three outliers with the highest ESIPT barriers for each position correspond to –NO group substitution, where the first excited state corresponds to a different type of excited state which is lower in energy, being n → π* type arising from an excitation of the lone pair on the –NO group as opposed to π → π* type as is the case for the standard CPIQ molecule (and other substituted versions). We still consider the ESIPT barrier on this state, as we assume that the photochemistry will occur after rapid relaxation to this dark state, and the ESIPT process is still both favourable and occurs with a barrier on this state, with a transition state corresponding to the proton transfer having been located. The barrier is however much higher to ESIPT on this state, between 72 and 96 kJ/mol depending on the substitution position.

Next looking at EDG substitution, the Stokes shift is modulated by even smaller amounts than with EWG substitution, matching the thermal properties in that EWGs seem to have a smaller input than EDGs. In general, EDG substitution
tends to slightly increase the Stokes shift, with the largest impact seen at position 3. If we look at the barrier to ESIPT, the majority of substitutions again have very small impact again on this barrier, although a few specific molecules have significantly larger barriers. The highest barriers for positions 3 and 13 arise from substitution of an acylamidophenyl group (–NHCOC<sub>6</sub>H<sub>6</sub>), where the higher barrier is due to a stabilisation of the initial tautomer optimised in the excited state arising from a twist around the N–C bond, found to be unfavourable for the ESIPT TS and an alternative tautomer in the excited state. The other high value of the ESIPT barrier at position 13, around 65 kJ/mol, arises from substitution of an acryloyloxy group (–OCOCHCH<sub>2</sub>).

Moving onto push-pull substitution, the Stokes shift is in general even less affected than for single EDG or EWG substitution, with no changes larger than 0.2 eV, with the exception of a group of molecules with extremely high Stokes shifts between 2.2 and 2.5 eV. These high Stokes shifts arise from combinations of any EDG at position 13 (excepting alkoxy, –OR, which have *f* < 0.1) with the –NO<sub>2</sub> EWG group at position 3. This is due to a significant additional geometrical change which occurs in the excited state - in the ground state the favourable orientation of the –NO<sub>2</sub> is in plane with the rest of the molecule, but in the excited state a twist occurs of this –NO<sub>2</sub> group, with the most favourable orientation being perpendicular to the plane of the aromatic system. This leads to a larger stabilisation of the excited state and destabilisation of the ground state of the tautomeric form at this excited state geometry, and thus a smaller energy gap and lower-energy emission, resulting in a larger Stokes shift. While this twist is also seen in the single substitution of –NO<sub>2</sub> at position 3, the first excited state only has a high enough oscillator strength (in this case, *f* > 0.2) when this EWG substitution is combined with certain EDG at position 13, thus this high Stokes shift is characteristic of these push-pull substitutions.

Looking now at the barrier to ESIPT, we can also note that while single EWG or EDG substitution has little effect aside from a few outliers, combining both types of substitution can allow for a large variety of ESIPT barriers ranging between 27 to 105 kJ/mol - especially for EWG(3)-EDG(13) substitution where a good spread of intermediate ESIPT barriers is observed. The opposite substitution pattern, EDG(3)-EWG(13) leads to a lower spread of energy barriers, although still larger than single-substitution alone, apart from a group of molecules with very high ESIPT barriers ca. 100 kJ/mol. These correspond to the combination of any EDG at position 3 and substitution of the –NO EWG in position 13, arising again due to differences in the excited state character, as for single –NO substitution.

Finally, we also tried to correlate Hirshfeld charges for the Franck-Condon geometry (N) of the molecule with these properties, the results of all atoms again being found in the SI (Section 4). However, we found that, unlike for ground-state reactivity, there is in general no correlation between charges and Stokes shift or ESIPT barrier, investigating both the atomic charges for the ground state and the first excited state. Excited-state reactivity appears to be much harder to predict using a simple approach, compared to the ground-state.

For this second molecule, we have chosen to demonstrate the use of our protocol for a high volume study of substituted compounds with few substitution positions but many possible EWG and EDG groups. Thanks to our protocol, we are able to investigate both ground and excited-state intramolecular proton transfer occurring in the CPIQ molecule, with results for 169 different substituted versions of CPIQ. This allowed us to have enough data to investigate trends which arise with different substitution types at different positions, as well as establish correlations between electronic properties and ground-state reactivity, which can now allow us to predict the ground-state reactivity of new CPIQ derivatives. This investigation entailed optimisation of 6 geometries for each substituted molecule, and thus over 1000 calculations. The creation by hand of inputs for this volume of calculations is clearly impractical - taking the same assumption of ca. 1 minute to generate an input file by hand, this would take over 16 hours. While our protocol still took several hours to generate all the initial geometries automatically, the main advantage is that this can be run without user-input, running remotely on a standard high-performance computing cluster. In addition, the optimisation of transition states is often complex without a good initial guess, limiting the use of molecule generation softwares such as OpenBABEL.<sup>57</sup> which can generate molecular geometries from connectivity information (SMILES), but are limited to "reasonable 3D structures". Our protocol uses the transition state of the unsubstituted molecule as a starting point allowing for successful transition state optimisations in both the ground and excited states. When optimising transition states, it is common to use one of many existing interpolation methods, which can often provide a good guess to these structures, starting from geometries of the reactants and the products.<sup>58–62</sup> However, these methods can often fail - starting from a previously existing optimised transition state geometry as in our protocol can allow for a more robust identification of transition states in substituted molecules.

## 3.3 Salicylidene Aniline (SA)

<figure><img src="image_14.png" /><figcaption>Figure 10: Overview of the photochemistry of the salicylidene aniline photochrome. The key geometries required to model the identified properties are circled in green.</figcaption></figure>

Finally, we shall demonstrate the use of our protocol for an example of a theoretical investigation to select best-candidate molecules with particular properties, based on substitutions to the well-known salicylidene aniline (SA) photochrome (Figure 1(c)). Substituted derivatives of this molecule are easily accessible at any position on either the phenolic or
<figure><img src="image_15.png" /><figcaption>Figure 11: Impact of substitution of electron withdrawing (EWG) and donating (EDG) groups at different positions of the SA molecule on the selected target properties: (a) Enol Excitation Energy (VEE) (b) ratio in polarity (dipole moments) between the Enol and Keto-*trans* (DM<sub>trans-Keto</sub>/DM<sub>Enol</sub>), and (c) energy barrier to *trans*-cis isomerisation of Keto-SA (Δ*E*<sub>isom</sub>). The x-axis labels correspond to the substitution positions of either the EWG, EDG, or both; numbering convention follows that of Figure 1, with O, M and P referring to *ortho*, *meta*, and *para* substitution respectively. For push-pull compounds, consisting of both EWG and EDG substitution on opposite sides of the molecule, the EWG position is given first, followed by that of the EDG.</figcaption></figure>

anilic ring, due to the nature of the coupling reaction through which SA can be synthesised.³³

In this case, the three properties which we would like to
optimise are:

1. Absorption wavelength ($\lambda_{abs}$)

2. Ratio of polarity between photoswitch forms (DM<sub>trans-Keto</sub>/DM<sub>Enol</sub>, DM = dipole moment)

3. Barrier to reverse isomerisation (Δ*E*<sub>isom</sub>)

These three properties can be modelled theoretically using
the three key geometries highlighted in green in Figure 10,
which describes the generally accepted photochemistry of
the SA molecule, where initial ESIPT is followed by iso-
merisation about the phenolic CCCN dihedral coordinate to
form the *trans*-keto photochromic product. We highlight that
these geometries which we use in our investigation are based
on prior literature of the SA molecule. ²⁸⁻³⁰ We will suppose
that, for this example investigation, we wish to identify pho-
tochromes with: (1) a longer wavelength of absorption $\lambda_{abs}$

than unsubstituted SA, calculated at 332 nm (3.73 eV) : we
shall choose a desired threshold of $\lambda_{abs} > 350$ nm (VEE <
3.54 eV), (2) a good ratio of polarity between the two forms
DM<sub>trans-Keto</sub>/DM<sub>Enol</sub> > 2, and (3) a lower barrier to re-
verse *trans*-keto → *cis*-keto isomerisation : we shall choose
Δ*E*<sub>isom</sub> < 75 kJ/mol, corresponding to a harmonic transition
state theory based unimolecular half-life of 1.6 seconds.

We use our protocol to investigate the impact of EWG and
EDG substitution at multiple different positions, as shown in
Figure 1(c), substituting EWG on the phenolic ring (positions
1, 2, 3 and 6) and EDG on the anilic ring (21, 23, 25). The
combination of multiple positions and EWG/EDG types led
us to an investigation of a total of 717 different derivatives of
SA, an investigation clearly impossible without an automated
approach.

Figure 11 shows the impact on the three selected properties
of different substitution positions and types. Beginning with
the excitation energy of the photochrome, we see that all
modifications tend to lower this value, by as much as ca.
1 eV. For EWG, substitution on positions 3 and 2 (ortho and para to the CCNC bridge) have a larger effect to lower this energy, while meta substitutions (1 and 6) have a much smaller impact. For EDG substitution on the opposing anilic ring, substitutions para (21) and ortho (25) to the central bridge again tend to lower the excitation energy slightly more than at the meta (23) position, although this trend is less evident than for EWG. For all substitution positions we also note a group of significantly lower excitation energies - these correspond to -NR₂ substitutions.

For push-pull substitution (EWG on the phenolic ring, EDG on the anilic ring), the combination of EWG and EDG appears to behave as expected from the single-substitution results. *Meta* EWG (1,6) substitution, especially on position 1, leads to the lowest variations of excitation energy, with dual-*meta* substitution (M(1)-M(23)) leading to the smallest impact. *Para* (2) and *ortho* (3) EWG lead to larger variations, especially when combined with *para* and *ortho* EDG (21, 25), where we see the largest impacts on the electronic excitation energy. The EWG clearly has a larger impact on the absorption energy than the EDG, although for each EWG position the general trend of changing EDG clearly follows that of EDG alone (P,O > M).

We can see that through substitution, we can modulate the excitation energy by a reasonable amount, to as low as 2.8 eV (443 nm). However, it remains to be seen how many of these low excitation energies are also associated with favourable values of the other two target properties. Finally, one final set of data is given for push-pull substitution which is an example of the inverse EWG and EDG positioning, with the EWG instead substituted on the anilic ring (at position 25, *ortho*) while the EDG is placed on the phenolic ring (at position 6, *meta*). The excitation energy of all of these molecules is systematically decreased, but by only ca. 0.25 eV, with a lower spread of values than for the opposite substitution of EWG and EDG (M(6)-O(25)).

Now, we can look at the ratio of dipole moments, where we see a spread from values lower than 1 (swapped enol and *trans*-keto polarities) to larger than 10. We can see that substitutions, especially EWG and push-pull, have a general tendency to decrease this ratio, although there are some notable exceptions. Within EWG substitutions, only a few molecules have a ratio higher than the unsubstituted molecule (orange line, 2.42). These correspond to -COOR substitutions, apart from the two highest *ortho* (3) points which arise from -CCl₃ and -CF₃ substitution. In contrast, EDG groups on the anilic ring actually increase the dipole moment ratio for the majority of cases, up to values as high as 5, with *meta* (23) substitution being the best choice here to increase this ratio.

For push-pull substitutions, the majority of combinations again decrease the dipole moment ratio, apart from certain substitutions with *ortho* (3) and *meta* (1) EWG (although not for *meta* (6) EWG) combined with any EDG substitution but particularly *meta* (23) and *ortho* (25). In these cases, there are a number of substitutions with ratio of dipole moments significantly increased, to as high as 11. These points all arise from substitutions including halide atoms (-F, -Cl, -CF₃, -CCl₃) as the EWG. Finally, we can again look at the

results for swapped EWG and EDG (O(25)-M(6)). For the dipole moment ratio, swapping the EWG and EDG groups to the opposite rings leads to a systematic decrease in this ratio, with the values centred around 1, meaning the initial enol and photochemically formed *trans*-keto forms have identical polarity - this would prevent their use in applications which rely on a measurable polarity difference.

Finally, we can look at the results for the energy barrier to the *trans*-keto → *cis*-keto isomerisation. We see a general tendency to decrease this barrier when substituting EWG on the phenolic ring and EDG on the anilic - the only systematic increase in energy barrier we see is for (O(25)-M(6)) push-pull substitution, where the EWG and EDG are swapped to the opposite rings. Substitutions can have a significant impact on this energy barrier, lowering it by as much as 43 kJ/mol, from 97 kJ/mol in the unsubstituted molecule to the lowest calculated barrier of 53.7 kJ/mol. As for the dipole moment ratio, we see that the trends in push-pull substitution are mostly defined by those of the EWG position, in this case lowering the barrier heights more for *ortho* (3) and *meta* (1, 6) substitutions; EDG substitution on the other hand has a less significant impact, leading to no clear systematic decrease in barriers aside from a small shift downwards for *para* EDG (21) substitutions, less clearly visible in push-pull compounds where the trends are dominated by the EWG. The largest decreases in barrier clearly arise from combinations of *ortho* EWG (3) substitutions, with the next largest arising from *meta* EWG (1, 6), followed by *para* EWG (2). The lowest energy barriers in general arise for EWG substitutions of -NO, -NO₂, -CCl₃ and -SO₂CH₃ - the stronger EWGs chosen.

<table><caption>Table 2: Best-candidate molecules selected from those investigated according to the three thresholds - $\lambda_{abs} > 350$ (VEE < 3.54 eV), $DM_{trans-Keto}/DM_{Enol} > 2$, $\Delta E_{isom} < 75$ kJ/mol. Substitution family notes the type of substitution pattern, matching the naming convention from Figure 11. Each row corresponds to one selected compound, noting the EWG group and position (Pos.), along with those of the EDG substitution for Push-Pull compounds.</caption><thead><tr><th>Substitution Family</th><th>EWG Group</th><th>EDG Group</th></tr></thead><tbody><tr><td rowspan="4">EWG O(3)</td><td>-CCl<sub>3</sub></td><td></td></tr><tr><td>-COOCH<sub>3</sub></td><td></td></tr><tr><td>-COOCHCH<sub>2</sub></td><td></td></tr><tr><td>-COOH</td><td></td></tr><tr><td rowspan="4">O(3)-O(25)</td><td>-CF<sub>3</sub></td><td>-OCH<sub>2</sub>CH<sub>3</sub></td></tr><tr><td>-CF<sub>3</sub></td><td>-C<sub>6</sub>H<sub>6</sub></td></tr><tr><td>-SO<sub>2</sub>CH<sub>3</sub></td><td>-OCH<sub>3</sub></td></tr><tr><td>-SO<sub>2</sub>CH<sub>3</sub></td><td>-OCH<sub>2</sub>CH<sub>3</sub></td></tr><tr><td rowspan="3">O(3)-P(21)</td><td>-CF<sub>3</sub></td><td>-CH<sub>2</sub>CH<sub>3</sub></td></tr><tr><td>-CF<sub>3</sub></td><td>-CH<sub>2</sub>CH<sub>2</sub>CH<sub>3</sub></td></tr><tr><td>-CF<sub>3</sub></td><td>-OCHCH<sub>2</sub></td></tr><tr><td rowspan="2">M(6)-M(23)</td><td>-NO<sub>2</sub></td><td>-OCH<sub>3</sub></td></tr><tr><td>-NO<sub>2</sub></td><td>-OCH<sub>2</sub>CH<sub>3</sub></td></tr></tbody></table>

We can now suppose that, from our calculated range of compounds, we wish to select some best-candidates which have favourable values of all three properties: excitation en-
ergy VEE < 3.54 eV, dipole moment ratio > 2, and energy barrier < 75 kJ/mol. We can simply set these as threshold values, and search our database of calculated molecules for those which fit these criteria. With these thresholds, we identify 13 molecules which fulfil these desired criteria, given in Table 2 (see also SI, Section 5).

Four EWG substitutions are selected all on position 3 (*ortho* to the bridging CCNC), the position with the largest possible impact on all three properties. No EDG compounds fulfil the desired criteria, as single EDG substitution has a lesser impact on all three desired properties as seen in Figure 11 - in particular, no EDG substituted compound has a low enough energy barrier ($ΔE_{isom} < 75$ kJ/mol). Finally, the identified push-pull molecules include 4 compounds with *ortho-ortho* (O(3)-O(25)) substitution pattern, 2 with *meta-meta* (M(6)-M(23)) pattern, and 3 with *ortho-para* (O(3)-P(21)) pattern.

## 4. CONCLUSIONS

In this work, we have demonstrated the use of a simple computational protocol for the high-throughput calculation of large numbers of molecules based on substitution to known chromophores which undergo excited-state intramolecular proton transfer. We selected three molecules for this investigation, demonstrating the use of our automated scheme at three levels. For the BZT molecule, we demonstrate how automation of calculation input files can both save time for a small focused investigation of a small number of substitutions, and more importantly avoid potential for human error when building geometries. For a small investigation of 12 molecules such as this, we show our protocol can be run locally on a laptop within only a few minutes. For the CPIQ molecule, we demonstrate how our protocol can be used to investigate trends in properties of molecules and their reactions, investigating both ground state and excited state reaction barriers as well as absorption and relative tautomer stabilities, allowing us to evaluate large volumes of molecules and extract trends for the influence of electron withdrawing and donating groups substituted onto this molecular core. The automation of a large volume of calculations also allowed us to investigate the potential for the prediction of molecular properties in these molecules using simple atomic charges, which would allow the prediction of ground-state reactivity for additional substituted molecules without need to optimise reaction barriers. Finally, for the SA molecule we demonstrate how our protocol could be used in a high-throughput investigation of a molecule with numerous potential substitution positions, in order to select best-candidate molecules for a particular application based on required properties. For the second and third applications, our protocol is run remotely on a high-throughput computational cluster.

For BZT, thanks to our automated investigation we demonstrate that substitution of both the strong and weak electron donating groups (–CH₃, –OH) lead to red-shifted absorption (lower VEE) for the BZT molecule, with the exception of *ortho*-OH, while electron withdrawing groups either lower or increase the VEE depending on the strength and position. The most red-shifted absorptions come from the strong EWG substitution of –NO₂ at either *ortho* or *meta* position, while weak EWG –Cl substitution at these positions leads to a

slight blue-shift in absorption. Finally, for *ortho* and *para* substitution we find a trend in the location of the proton which undergoes ESIPT in this molecule, being relatively closer to the donating oxygen than the accepting nitrogen with increasing electron-donating character. The proton is also systematically relatively closer to the accepting nitrogen when *para* substituted compared to *ortho*, suggesting *para* substitution may favour proton transfer relative to *ortho*. *Meta* substitution appears to have little impact on the relative positioning of the proton between the acceptor and donor atoms.

For CPIQ, thanks to the large volumes of calculations made possible by our automation protocol, we are able to establish trends between the substitution position/type and properties in the ground and excited state of this photochrome. EWG and EDG, especially combined in push-pull substitution, can significantly modulate the reactivity on the ground and excited states, decreasing and increasing barriers to GSIPT and ESIPT depending on the substitution position and type. Our dataset allows us to extract strong correlations between atomic charge of the initial tautomer and ground state properties, with different atoms being the best predictors depending on the position and type of the substitution, although substitution at position 3 (i.e. *ortho* to the –NH₂ group) is generally less predictable due to steric effects. The best predictive atomic charges are also often unexpected, sometimes being aliphatic carbons or hydrogens. Only in two cases are the charges on the proton undergoing GSIPT/ESIPT and the acceptor nitrogen good predictors of reactivity, and the charge on the donor nitrogen is never a good predictor. We also find that, while ground state reactivity often correlates well with the atomic charges, excited state reactivity cannot be predicted in the same manner, no matter whether ground state or excited state atomic charges are investigated.

For SA, the large volumes of calculations made possible by our automation protocol again allow for the establishment of trends between molecular properties and reactivity with substitution position and type. Substitution tends to lower the VEE of the SA molecule, no matter the position or type, leading to red-shifted absorption by as much as 0.9 eV. The ratio of dipole moments between the initial tautomer and the photochromic *trans*-keto product tends to be reduced when substituting EWGs on the phenolic ring and increased with EDGs on the anilic ring, and push-pull substitutions generally have lower ratio aside from EWG including halide atoms where the ratio can be as high as 11 (compared to 2.42 in unsubstituted SA). The barrier to *trans-cis*-keto isomerisation is generally lowered with EWG on the phenolic ring, with EDG substitution having a lesser impact on this barrier, and the lowest barriers arise from EWG substitution *ortho* to the CCNC bridge which can reduce the barrier by as much as 43 kJ/mol. Swapping the substitutions, such that EWG is on the anilic ring and EDG on the phenolic, leads to the opposite impact, increasing this isomerisation barrier.

Finally, we must discuss the limitations of this automated approach. First, the automated creation of large volumes of inputs necessitates either cheap calculations or a good availability of computational resources. In addition, it is still necessary to ensure the calculations converge correctly, es-
pecially for transition states (TS) where in some cases in
this study an incorrect TS geometry was optimised despite
starting from the guess of the unsubstituted molecule. Us-
ing automated data analysis, we were able to locate when
incorrect TS were found, based on key geometrical coordi-
nates - for instance, the SA isomerisation TS should have a
C3C10C13N14 dihedral angle close to 90°. These calculations
then had to be modified and rerun by hand - although this only
occurred in a small percentage of cases (for SA, in ca. 6% of
molecules). When automating high-throughput calculations
in this way, it is crucial to establish and run automatisable
checks on final results for this reason, since manually check-
ing outputs for such a volume of calculations is unfeasible. In
the same way, we had to ensure that for substitutions which
lead to multiple conformers, the same conformer was taken
when this is correct to do so. For this reason, we intend to
extend the automatisation protocol to include an option to
automatically check and adjust certain internal coordinates,
to ensure geometry calculations for a given substitution do
not optimise to different and thus incomparable conformers,
or to incorrect TS geometries. Another potential improve-
ment would be the extension of our protocol to account for
the potential coexistence of multiple conformers, through in-
terfacing with a search algorithm such as CREST,⁶³ which
could be used to generate and account for thermally acces-
sible conformers. Finally, the most important drawback of
this automated approach is that it assumes that the reactiv-
ity of substituted molecules will be the same as that of the
unsubstituted one, in that the same critical geometries are
taken. This is generally a reasonable approximation to take,
but when using this approach to select candidate molecules,
it is recommended to subsequently carry out investigations
by hand on these selected candidates, to ensure other reac-
tion pathways are not made possible by the addition of the
substitution.

The automated protocol, along with a short descrip-
tion of installation and use, is available for download at
github.com/issymerritt/generate_subs. We note also that the
way the protocol is written would allow for a simple adapta-
tion to other quantum chemistry programs - information on
how this would be possible is also clearly given. Future work
will aim to include automated checking of optimisations di-
rectly into the protocol, allowing for additional calculations
in case of optimisation failure to be also generated automati-
cally by the protocol. This will be relatively trivial to include.
When calculations fail to optimise at all, the lowest energy
point can then be taken as a starting point for follow-on cal-
culation, while if the optimisation converges to an incorrect
stationary point with too many imaginary frequencies, a mod-
ified geometry can be taken as a start point by automatically
reading and deforming the geometry along the nuclear mode
associated with the imaginary frequency.

DATA AND SOFTWARE AVAILABILITY

The complete set of scripts required to run the protocol is
available on GitHub: github.com/issymerritt/generate_subs.
Instructions for running the protocol are provided in the
README.md file, as well as in the electronic supporting
information. Installation of the protocol is also outlined in

the README.md file, as well as in the ESI. Raw data for
the results of the three investigations, on BZT CPIQ and SA,
are also given in .json format within a zip folder, provided as
supporting information.

SUPPORTING INFORMATION AVAILABLE

Annotated example of input file used to run the substitution
generator protocol. Description of the code layout and ad-
ditional scripts for analysis included with the protocol. De-
scription of the use of Harmonic (Eyring) Transition State
Theory to obtain unimolecular half-lives. Additional infor-
mation for BZT: orbitals describing the first excited state for
BZT and the first two excited states of *ortho*-OH substituted
BZT, correlations between HOMO-LUMO energy gap and
vertical excitation for substituted BZT both with and without
NO₂ group substitution. Additional information for CPIQ:
list of EDG/EWG groups tested, images of the “out-of-plane”
motion of the acylamido group when substituted at position
13, *R*² correlation values for atomic charges on every atom
with the five properties studied for both ground and excited
state Hirshfeld charges. Additional information for SA: list
of EDG/EWG tested, diagrams of the 13 selected molecules
using the thresholds outlined in the text. Database0s (in
json format) containing calculated properties for all modified
CPIQ and SA molecules.

ACKNOWLEDGEMENTS

This work was performed using computational resources
from GENCI-TGCC and GENCI-IDRIS (Grant 2024-
103754) and the MCIA (Mésocentre de Calcul Intensif
Aquitain) at the University of Bordeaux. The authors ac-
knowledge the support of the French Agence Nationale de
la Recherche (ANR), under grant ANR-23-CE24-0001-02
(project De-TECT) as well as the GPR (grand programme
de recherche) LIGHT of the University of Bordeaux. The
authors would like to thank J.-L. Pozzo and D. Bassani for
discussions on experimentally viable substitutions of CPIQ
and SA.

Author contributions

I.C.D.M.: conceptualization, investigation, software and
script writing, manuscript writing, editing, review. F.C: re-
view and editing, supervision, funding acquisition, direction
and planning.

Notes

The authors declare no competing financial interests.

REFERENCES

(1) Kohn, W.; Sham, L. J. Self-Consistent Equations In-
cluding Exchange and Correlation Effects. *Phys. Rev.*
**1965**, *140*, 1133–1138.

(2) Thompson, N.; Manso, G.; Ge, S. The Importance of (Exponentially More) Computing Power. *Academy of Management Proceedings* **2023**, *1*.

(3) Nepal, N. K.; Canfield, P. C.; Wang, L.-L. HT-ESP (High-throughput electronic structure package): A package for high-throughput ab initio calculations. *CComput. Mater. Sci.* **2024**, *244*, 113247.
(4) Fortunato, N. M.; Taubel, A.; Marmodoro, A.; Pfeifer, L.; Ophale, I.; Ebert, H.; Gutfleisch, O.; Zhang, H. High-Throughput Design of Magnetocaloric Materials for Energy Applications: MM'X alloys. *Adv. Sci.* **2023**, *10*, 2206772.

(5) Bhoorasingh, P. L.; Slakman, B. L.; Seyedzadeh Khandshan, F.; Cain, J. Y.; West, R. H. Automated Transition State Theory Calculations for High-Throughput Kinetics. *J. Phys. Chem. A* **2017**, *121*, 6896–6904.

(6) Sora, V.; Laspiur, A. O.; Degn, K.; Arnaudi, M.; Utichi, M.; Beltrame, L.; De Menezes, D.; Orlandi, M.; Stoltze, U. K.; Rigina, O.; Sackett, P. W.; Wadt, K.; Schmiegelow, K.; Tiberti, M.; Papaleo, E. RosettaDDGPrediction for high-throughput mutational scans: From stability to binding. *Protein Sci.* **2022**, *32*, e4527.

(7) Ahmad, I.; Jagatap, V.; Patel, H. *Phytochemistry, Computational Tools and Databases in Drug Discovery*; Elsevier, 2023; p 371–392.

(8) Liu, Z.-Y.; Hu, J.-W.; Huang, T.-H.; Chen, K.-Y.; Chou, P.-T. Excited-state intramolecular proton transfer in the kinetic-control regime. *Phys. Chem. Chem. Phys.* **2020**, *22*, 22271–22278.

(9) Bañuelos, J. BODIPY Dye, the Most Versatile Fluorophore Ever? *Chem. Rec.* **2016**, *16*, 335–348.

(10) Jelínková, V.; Dellai, A.; Verwaerde, L.; Rodriguez, V.; Fecková, M.; Vachtlová, M.; Podlesný, J.; Klikar, M.; Sempé, B.; Hugget, M.; Pařík, P.; Pytela, O.; Castet, F.; Růžičková, Z.; Bureš, F. Property Tuning in N-Methylpyrrole Azo-Photoswitches via Modification of the Peripheral Substituents. *Chem. Eur. J.* **2025**, *31*, e202404221.

(11) Meeto, W.; Suramitr, S.; Lukeš, V.; Wolschann, P.; Hannongbua, S. Effects of the CN and NH2 substituents on the geometrical and optical properties of model vinylfluorenes, based on DFT calculations. *Journal of Molecular Structure: THEOCHEM* **2010**, *939*, 75–81.

(12) Khan, R. K.; Turlapati, S.; Begum, N.; Mohiuddin, G.; Rao, N. V. S.; Ghosh, S. Impact of terminal polar substitution on elastic, electro-optic and dielectric properties of four-ring bent-core nematic liquid crystals. *RSC Advances* **2018**, *8*, 11509–11516.

(13) Ying, P.; Zhang, J.; Zhang, X.; Zhong, Z. Impacts of Functional Group Substitution and Pressure on the Thermal Conductivity of ZIF-8. *The Journal of Physical Chemistry C* **2020**, *124*, 6274–6283.

(14) Turan, H. T.; Yavuz, I.; Aviyente, V. Understanding the Impact of Thiophene/Furan Substitution on Intrinsic Charge-Carrier Mobility. *The Journal of Physical Chemistry C* **2017**, *121*, 25682–25690.

(15) Muhammad, S.; Al-Sehemi, A. G.; Irfan, A.; Chaudhry, A. R.; Gharni, H.; AlFaify, S.; Shkir, M.; Asiri, A. M. The impact of position and number of methoxy group(s) to tune the nonlinear optical properties of chalcone derivatives: a dual substitution strategy. *Journal of Molecular Modeling* **2016**, *22*, 73.

(16) Mi, S.; Lu, N. Impact of substitution on the reorganization energy of bis-triarylamine derivatives. *Journal of Molecular Structure: THEOCHEM* **2010**, *940*, 1–5.

(17) Yang, H.; Cui, C.; Li, Y. Effects of Heteroatom Substitution on the Photovoltaic Performance of Donor Materials in Organic Solar Cells. *Accounts of Materials Research* **2021**, *2*, 986–997.

(18) Ertl, P.; Altmann, E.; Racine, S.; Lewis, R. Ring replacement recommender: Ring modifications for improving biological activity. *European Journal of Medicinal Chemistry* **2022**, *238*, 114483.

(19) Joshi, H. C.; Antonov, L. Excited-State Intramolecular Proton Transfer: A Short Introductory Review. *Molecules* **2021**, *26*, 1475.

(20) Sedgwick, A. C.; Wu, L.; Han, H.-H.; Bull, S. D.; He, X.-P.; James, T. D.; Sessler, J. L.; Tang, B. Z.; Tian, H.; Yoon, J. Excited-state intramolecular proton-transfer (ESIPT) based fluorescence sensors and imaging agents. *Chem. Soc. Rev.* **2018**, *47*, 8842–8880.

(21) Fu, P.-Y.; Yi, S.-Z.; Pan, M.; Su, C.-Y. Excited-State Intramolecular Proton Transfer (ESIPT) Based Metal–Organic Supramolecular Optical Materials: Energy Transfer Mechanism and Luminescence Regulation Strategy. *Acc. Mater. Res.* **2023**, *4*, 939–952.

(22) Yao, D.; Zhao, S.; Guo, J.; Zhang, Z.; Zhang, H.; Liu, Y.; Wang, Y. Hydroxyphenyl-benzothiazole based full color organic emitting materials generated by facile molecular modification. *J. Mater. Chem.* **2011**, *21*, 3568.

(23) Azarias, C.; Budzák, S.; Laurent, A. D.; Ulrich, G.; Jacquemin, D. Tuning ESIPT fluorophores into dual emitters. *Chem. Sci.* **2016**, *7*, 3763–3774.

(24) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseeman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams-Young, D.; Ding, F.; Lipparini, F.; Egidi, F.; Gings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian-16 Revision C.01. 2016; Gaussian Inc. Wallingford CT.

(25) Neese, F. The ORCA program system. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2012**, *2*, 73–78.

(26) Cantwell, M. G.; Sullivan, J. C.; Burgess, R. M. Persistent Organic Pollutants (POPs): Analytical Techniques, Environmental Fate and Biological Effects; Elsevier, 2015; p 513–545.

(27) Avagyan, R.; Luongo, G.; Thorsén, G.; Östman, C. Benzothiazole, benzotriazole, and their derivates in cloth-
ing textiles—a potential source of environmental pollutants and human exposure. *Environ. Sci. Pollut. Res.* **2014**, *22*, 5842–5849.

(28) Pijeau, S.; Foster, D.; Hohenstein, E. G. Effect of Non-planarity on Excited-State Proton Transfer and Internal Conversion in Salicylideneaniline. *J. Phys. Chem. A* **2018**, *122*, 5555–5562.

(29) Ortiz-Sánchez, J. M.; Gelabert, R.; Moreno, M.; Lluch, J. M. Electronic-structure and quantum dynamical study of the photochromism of the aromatic Schiff base salicylideneaniline. *J. Chem. Phys.* **2008**, *129*, 214308.

(30) Silfies, M. C.; Mehmood, A.; Kowzan, G.; Hohenstein, E. G.; Levine, B. G.; Allison, T. K. Ultrafast internal conversion and photochromism in gas-phase salicylideneaniline. *J. Chem. Phys.* **2023**, *159*, 104304.

(31) Ziółek, M.; Kubicki, J.; Maciejewski, A.; Naskręcki, R.; Grabowska, A. An ultrafast excited state intramolecular proton transfer (ESIPT) and photochromism of salicylideneaniline (SA) and its “double” analogue salicylaldehyde azine (SAA). A controversial case. *Phys. Chem. Chem. Phys.* **2004**, *6*, 4682–4689.

(32) Ziółek, M.; Burdziński, G.; Douhal, A. Long-living structures of photochromic salicylaldehyde azine: polarity and viscosity effects from nanoseconds to hours. *Photochem. Photobiol. Sci.* **2012**, *11*, 1389–1400.

(33) Schiff, H. Mittheilungen aus dem Universitätslaboratorium in Pisa: Eine neue Reihe organischer Basen. *Liebigs Ann.* **1864**, *131*, 118–119.

(34) Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K. PySCF: the Python-based simulations of chemistry framework. *Wiley Interdiscip. Rev. Comput. Mol. Sci.* **2017**, *8*, e1340.

(35) Rodrigues, O. Des lois géométriques qui régissent les déplacements d'un système solide dans l'espace, et de la variation des coordonnées provenant de ces déplacements considérés indépendamment des causes qui peuvent les produire. *J. Math. Pures Appl.* **1840**, 380–440.

(36) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. *J. Chem. Phys.* **1999**, *110*, 6158–6170.

(37) Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. *The Journal of Chemical Physics* **1982**, *77*, 3654–3665.

(38) Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self—Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian—Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. *The Journal of Chemical Physics* **1972**, *56*, 2257–2261.

(39) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. *The Journal of Chemical Physics* **1984**, *80*, 3265–3269.

(40) Clark, T.; Chandrasekhar, J.; Spitznagel, G. W.; Schleyer, P. V. R. Efficient diffuse function-augmented

basis sets for anion calculations. III. The 3-21+G basis set for first-row elements, Li–F. *Journal of Computational Chemistry* **1983**, *4*, 294–301.

(41) Mennucci, B.; Cances, E.; Tomasi, J. Evaluation of Solvent Effects in Isotropic and Anisotropic Dielectrics and in Ionic Solutions with a Unified Integral Equation Method: Theoretical Bases, Computational Implementation, and Numerical Applications. *The Journal of Physical Chemistry B* **1997**, *101*, 10506–10517.

(42) Cossi, M.; Barone, V. Time-dependent density functional theory for molecules in liquid solutions. *The Journal of Chemical Physics* **2001**, *115*, 4708–4717.

(43) Bonardd, S.; Robles, E.; Barandiaran, I.; Saldias, C.; Leiva, A.; Kortaberria, G. Biocomposites with increased dielectric constant based on chitosan and nitrile-modified cellulose nanocrystals. *Carbohydr. Polym.* **2018**, *199*, 20–30.

(44) Azofeifa, D. E.; Arguedas, H. J.; Vargas, W. E. Optical properties of chitin and chitosan biopolymers with application to structural color analysis. *Opt. Mater.* **2012**, *35*, 175–183.

(45) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual molecular dynamics. *J. Mol. Graphics* **1996**, *14*, 33–38.

(46) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. *The Journal of Chemical Physics* **1993**, *98*, 5648–5652.

(47) Bojinov, V. B.; Grabchev, I. K. Novel functionalized 2-(2-hydroxyphenyl)-benzotriazole – benzo[de]isoquinoline-1, 3-dione fluorescent UV absorbers. *Journal of Photochemistry and Photobiology A: Chemistry* **2005**, *172*, 308–315.

(48) Crawford, J. 2(2-hydroxyphenyl)2H-benzotriazole ultraviolet stabilizers. *Progress in Polymer Science* **1999**, *24*, 7–43.

(49) Panahi, H. A.; Rabbani, M.; Zabarjad-Shiraz, N.; Mofavvaz, S.; Moniri, E.; Kanghari, S.; Entezari, M.; Hasanzadeh, A. Synthesis and application of 5-amino-2-benzotriazol-2-yl-phenol for preconcentration and determination of zinc (II) in water samples by flame atomic absorption spectrometry. *Desalination and Water Treatment* **2010**, *22*, 330–339.

(50) DeBellis, A. D.; Iyengar, R.; Kaprinidis, N. A.; Rodebaugh, R. K.; Suhadolnik, J. Service Life Prediction; American Chemical Society, 2001; p 453–467.

(51) Evans, M. G.; Polanyi, M. Some applications of the transition state method to the calculation of reaction velocities, especially in solution. *Trans. Faraday Soc.* **1935**, *31*, 875.

(52) Eyring, H. The Activated Complex in Chemical Reactions. *J. Chem. Phys.* **1935**, *3*, 107–115.

(53) Luchini, G.; Paton, R. S. Bottom-Up Atomistic Descriptions of Top-Down Macroscopic Measurements: Computational Benchmarks for Hammett Electronic Parameters. *ACS Phys. Chem. Au* **2024**, *4*, 259–267.

(54) Hirshfeld, F. L. Bonded-atom fragments for describing molecular charge densities. *Theor. Chim. Acta* **1977**, *44*, 129–138.

(55) Hammett, L. P. Some Relations between Reaction
Rates and Equilibrium Constants. *Chem. Rev.* **1935**, *17*, 125–136.

(56) Clayden, J.; Greeves, N.; Warren, S. *Organic Chemistry*, Second Edition; Oxford University Press, 2012.

(57) O'Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. *J. Cheminf.* **2011**, *3*, 33.

(58) E, W.; Ren, W.; Vanden-Eijnden, E. String method for the study of rare events. *Phys. Rev. B* **2002**, *66*, 052301.

(59) Burger, S. K.; Ayers, P. W. Methods for finding transition states on reduced potential energy surfaces. *J. Chem. Phys.* **2010**, *132*, 234110.

(60) Bonfanti, S.; Kob, W. Methods to locate saddle points in complex landscapes. *J. Chem. Phys.* **2017**, *147*, 204104.

(61) Ionova, I. V.; Carter, E. A. Ridge method for finding saddle points on potential energy surfaces. *J. Chem. Phys.* **1993**, *98*, 6377–6386.

(62) Peng, C.; Bernhard Schlegel, H. Combining Synchronous Transit and Quasi-Newton Methods to Find Transition States. *Isr. J. Chem.* **1993**, *33*, 449–454.

(63) Pracht, P.; Bohle, F.; Grimme, S. Automated exploration of the low-energy chemical space with fast quantum chemical methods. *Physical Chemistry Chemical Physics* **2020**, *22*, 7169–7192.
TOC Graphic

<figure><img src="image_16.png" /></figure>