Scientists’ searching for new and innovative materials with promising properties often finds their way back to the roots and fundamentals of nature. Porphyrins, for example, are engaged in many essential and complex processes of nature, such as photosynthesis or oxygen fixation in blood cells. Hence, porphyrins are excellent candidates as natural products to design novel compounds with a high adaptability to their environment (Kadish et al., 2000; Dolphin, 2012; Auwärter et al., 2015). Based on these tetrapyrrolic macrocycles, expanded porphyrins have been synthesized with an array of interesting properties, including a high conformational flexibility, novel coordination behaviours, versatile aromaticity and exceptional optical properties (Tanaka and Osuka, 2016). Expanded porphyrins are macrocycles containing more than four pyrrole rings (or analogous heterocyclic subunits), which are connected directly or by methine groups, keeping an internal pathway of at least 17 atoms (Sessler and Seidel, 2003). Different modifications to the core structure of expanded porphyrins can be introduced, either by 1) replacing one or more pyrrole rings with other five-membered heterocycles such as furan, thiophene or tellurophene (Chandrashekar and Venkatraman, 2003; Rath et al., 2005; Kumar et al., 2007; Chatterjee et al., 2017), or by 2) introducing different substituents at the meso-positions of the macrocycle (Marcos et al., 2014; Tanaka and Osuka, 2016; Woller et al., 2016; Torrent-Sucarrat et al., 2017b). These variations in the core structures have recently drawn a lot of attention as the modified structures exhibit outstanding properties for multidisciplinary applications such as bio-sensors, photodynamic therapy, photochemistry, catalysis and molecular electronics (Saito and Osuka, 2011; Sung et al., 2017; Cárdenas-Jirón et al., 2019).

The field of molecular electronics involves the use of individual molecules as active elements in the electronic circuit, such as diodes, rectifiers, transistors, or memory devices (Sun et al., 2014; Su et al., 2016; Xiang et al., 2016). Besides their potential to meet the ever-increasing demand for the miniaturization of electronics, more importantly, molecular electronics opens up the possibility of devices with novel functionalities beyond the complementary silicon-based technologies, such as molecular switches (Zhang et al., 2015). A considerable advantage over silicon-based electronics is that molecular electronics can be tailor-made using the concepts from self-assembly and thus becomes highly compatible with a large number of substrates (Su et al., 2016). As a consequence, the chemical and physical properties are easily tunable resulting in applications ranging from photopharmacology and drug delivery to imaging, spectroscopy and optoelectronics (Feringa and Browne, 2011; Feringa, 2017).

The switching devices consist of a molecular building block that can be reversibly shifted between two or more stable states upon application of an external stimulus, including light, temperature, bias voltage, chemical reactions and even mechanical forces (Zhang et al., 2015; Ladenthin et al., 2016). The prerequisite of a molecular switch is that the two or more states are separable and show a distinct difference in the property of interest (Zhang et al., 2015). Regarding optical switches, organic materials with commutable nonlinear optical responses are sought for optoelectronic applications, such as molecular-scale memory devices with multiple storage and non-destructive reading capacity (Kawata and Kawata, 2000). Interesting nonlinear (NLO) optical properties are second harmonic generation (SHG), two-photon absorption (TPA), and third harmonic generation (THG), among others. At the molecular level, the first hyperpolarizability (β) is a noteworthy quantity as this property is intimately connected to a variety of optical phenomena like SHG or the Hyper-Rayleigh Scattering (HRS). The potential of single molecules as nonlinear optical switches is mainly determined by the amplitude of the NLO contrast and thus large differences in first hyperpolarizability between the different states are required for high-performance optical switches (Castet et al., 2013; Kariduraganavar et al., 2021).

Expanded porphyrins are highly desirable for near-infrared dyes and nonlinear optical applications (Rath et al., 2005; Pawlicki et al., 2009; Mori et al., 2013; Shimomura et al., 2020; Wang et al., 2020). Thanks to their extended π-system, expanded porphyrins experience a red-shift of absorption bands and exhibit a much larger TPA cross section compared to regular porphyrins. As expanded porphyrins can easily switch between different redox states and π-conjugation topologies induced by diverse external stimuli, they are an ideal test bed for molecular optical switches (Tanaka and Osuka, 2016). Indeed, recent computational studies have shown that the interconversion between different π-conjugation topologies, more precisely, Hückel and Möbius structures, induces dramatic changes in first and second hyperpolarizabilities, acting therefore as a novel type of efficient nonlinear optical switches (Torrent-Sucarrat et al., 2017b; Woller et al., 2019). Importantly, the performance of these optical switches can be improved by introducing different types of functionalizations (Torrent-Sucarrat et al., 2017b; Pino-Rios and Cárdenas-Jirón, 2019; Woller et al., 2019). For example, it has been demonstrated that the NLO contrast in [28]hexaphyrin topological switches can be enhanced by introducing push-pull substitution patterns, which incorporate electron-withdrawing and electron-donating substituents on meso-positions located on the opposite sides of the macrocycles (Torrent-Sucarrat et al., 2017b).

However, exploring the entire region of synthetically accessible expanded porphyrins in the chemical compound space (CCS) remains an immense challenge. Hence, a rational approach is needed to help in the design of new high-potential molecular switches. In traditional molecular design procedures, one directly introduces the desired functionalizations into the chemical structure and subsequently evaluates for the altered structure the NLO quantities (Torrent-Sucarrat et al., 2012; Torrent-Sucarrat et al., 2017b; Tonnelé et al., 2018; Woller et al., 2019). One thus mainly relies on the chemist’s intuition to design derivatives with favorable NLO properties, which can become a cumbersome task when multiple chemical modifications are demanded. Alternatively, inverse molecular design approaches target the optimization of the property of interest as a function of the chemical structure (Wang et al., 2006; Balamurugan et al., 2008; De Vleeschouwer et al., 2012; De Vleeschouwer et al., 2015; De Vleeschouwer et al., 2016; Teunissen et al., 2017; Freeze et al., 2019). Only modifications are made that effectuate property enhancement. A big advantage of these inverse design approaches is that only a small subset of all possible expanded porphyrin derivatives needs to be evaluated to devise high-performing chemical structures for a particular application. Recent advances in computational inverse design approaches have already impacted the pool of structures available to nanophotonics (Molesky et al., 2018; Su et al., 2020).

In this paper, we take up the challenge to design efficient molecular NLO switches by means of an inverse design algorithm. In this way, the discovery of new functionalized expanded porphyrins with improved NLO contrasts can be accelerated. In particular, we focus on hexaphyrins, which are very effective platforms to realize versatile electronic states including Möbius and Hückel aromatic and antiaromatic species (Sankar et al., 2008; Ahn et al., 2005; Ishida et al., 2011; Kim et al., 2009; Sun et al., 2016). Two types of hexaphyrin-based molecular switches are considered here, namely redox (28R → 26R) and topology (28R → 28M) switches, the latter involving Hückel-to-Möbius transformations while keeping the number of π-electrons constant (Scheme 1). Among the different porphyrinoids investigated so far in our research group, the redox-triggered aromaticity switches based on hexapyrrolic macrocycles exhibit the largest ON/OFF transmission ratios in conductance switching applications (Alonso et al., 2012, 2014; Alonso et al., 2013a; Alonso et al., 2013b; Woller et al., 2016; Woller et al., 2018; Stuyver et al., 2018). In addition, topology switches based on [28]hexaphyrin exhibit dramatic variations in the second- and third-order NLO properties, as shown by Torrent-Sucarrat et al. (2017b).

Our objectives are threefold. First, we aim to establish the structure-property relationships between molecular topology, aromaticity, and optoelectronic properties for different unsubstituted hexaphyrins, labelled as the parent structures. Second, we investigate the effect of selected electron-withdrawing and electron-donating meso-substituents on the macrocyclic aromaticity and nonlinear optical properties. Lastly, we perform the inverse molecular design procedure with the Best-First Search (BFS) algorithm (Wang et al., 2006; Balamurugan et al., 2008; De Vleeschouwer et al., 2012, 2015; De Vleeschouwer et al., 2016; Teunissen et al., 2017) to find meso-substitution functionalizations of hexaphyrins inducing a maximum NLO contrast. By applying the inverse design protocol, only a fraction of all possible structures is calculated. Different substitution patterns as previously described in literature are considered, involving multiple sites to be modified, from a library of possible substituents.

The Best-First Search (BFS) algorithm consists of a site-by-site optimization during which chemical modifications from a fragment library are introduced into a predefined molecular scaffold (Wang et al., 2006; Balamurugan et al., 2008; De Vleeschouwer et al., 2012; De Vleeschouwer et al., 2015; Teunissen et al., 2017; De Vleeschouwer et al., 2016). For each site, the modification’s impact on the target property of interest is evaluated and the most favourable substituent is selected based on the largest improvement in property value. The algorithm relies on the independent site approximation (ISA), which assumes that the sites can be optimized one by one, yielding the optimum after a single iteration over all sites and regardless of the site sequence used. Nevertheless, for most properties, the ISA is not fully valid as sites are affected by the type of functionalization on other positions in the structure. Therefore, several site iterations are carried out to ensure that intersite influences are evaluated as well. Figure 1 shows an example of the BFS algorithm on the 28R ⇆ 26R redox switch with 3 pairs of sites (R_1,4, R_2,5 and R_3,6) and 3 substituents (-F, -CN and -OH) available for functionalization. The three pairs of sites consist of combinations of 2 meso-positions that face each other on the macrocycle. First, the algorithm generates a random initial structure within the hexaphyrin CCS, with the CCS magnitude defined by the number of sites and the size of the substituent library. To this end, the algorithm randomly selects for each site a substituent from the fragment library. In the example of Figure 1, the randomly chosen substituents are a cyano group for the first (R_1,4), a hydroxyl group for the second (R_2,5), and a fluoro group for the third site (R_3,6). Next, the algorithm starts to optimize site 1 (R_1,4), while keeping all other sites (R_2,5 and R_3,6) fixed, i.e. the substituents on all sites, but site 1, remain unaltered. Optimizing site 1 involves placing each of the substituents from its fragment library on that position and evaluating the property for each generated structure. The substituent with the largest improvement in β _HRS contrast (here: -F) is then introduced into the initial structure. This altered structure is subsequently used as the starting point for the next step in the algorithm, corresponding to the optimization of site 2 while keeping the other two positions fixed. The same procedure is repeated until all sites have been visited. Convergence is reached when a single site modification does not result anymore in an improved property value. Even though there exist more recently introduced inverse design methods that may outperform the BFS method based on efficacy or larger CCS exploration (Elward and Rinderspacher, 2015; Strasser et al., 2016; Häse et al., 2018), BFS has the key advantage of systematically delivering essential information on how the type and location of each possible substituent influences our property of interest, since in every global iteration all substituents are tested on all available sites (Teunissen et al., 2017).

The figure-of-merit of the performance of an NLO switch is the NLO contrast as defined by Eq. 1. This contrast considers the ratio of the first HRS hyperpolarizabilities of both the ON and OFF state structures having the same meso-substitution pattern and has a high value for a well-performing NLO switch. c o n t r a s t = β H R S ( O N ) β H R S ( O F F ) (1)

In the literature, this contrast expression has been widely used for the evaluation of the performance of molecular switches (Plaquet et al., 2008; Castet et al., 2013; Tonnelé et al., 2018; Woller et al., 2019). The contrast is maximized during the BFS procedures, but the maximization of the ratio entails one problem. If the denominator becomes zero, the entire contrast reaches a value of infinity making it impossible to compare other structures facing the same issue. This widely used contrast definition may therefore not be an ideal target objective in the inverse design maximizations. To circumvent this complication, we implemented the following condition within the BFS procedure. If the β _HRS of the OFF state becomes lower than 10 a.u., the denominator is set to a value of 0.001 a.u.

Another measure of NLO contrast has been presented, among others, in the work of Torrent-Sucarrat et al. (2017b), taking the difference of the NLO responses between the ON and OFF state (Eq. 2). A high difference is connected to a high and low β _HRS response for the ON and OFF state, respectively, and reflects similar behaviour as our previous contrast definition. Note that a value of zero for the OFF state does not pose any optimization issues. The main difference between both contrast definitions is that high β _HRS differences can also appear from molecular switches in which the OFF state’s β _HRS response is far away from zero a.u. In this work, we carried out additional BFS searches in which the β _HRS difference is maximized. The best performing switches are then compared with the optimal structures connected to the ratio in β _HRS . c o n t r a s t = β H R S ( O N ) − β H R S ( O F F ) (2)

All quantum-chemical calculations were executed with the Gaussian16 software package (Frisch et al., 2016). The geometry optimizations and frequency calculations were performed using the Minnesota hybrid functional M06-2X (Zhao and Truhlar, 2008) with the cc-pVDZ basis set (Dunning, 1989). In recent benchmark studies (Marcos et al., 2012; Torrent-Sucarrat et al., 2017a; Sylvetsky et al., 2020; Woller et al., 2020), the M06-2X functional emerges among the best contemporary density functionals to describe the relative energies of Hückel and Möbius expanded porphyrins. Although singly-twisted Möbius topologies are a challenging test for electronic structure methods, M06-2X provides errors close to chemical accuracy relative to the golden-standard canonical CCSD(T) calculations. Harmonic vibrational analyses show that all optimized structures are exclusively characterized by positive eigenvalues and thus represent minima on the potential energy surface. To obtain more accurate electronic energies, subsequent single-point energy calculations on the optimized geometries were performed at the M06-2X/cc-pVTZ level of theory.

For all unsubstituted hexaphyrins, several aromaticity descriptors were computed, including the magnetic susceptibility exaltation (Λ) (Dauben et al., 1968), the nucleus-independent chemical shifts (NICS) (Lazzeretti, 2004; Chen et al., 2005; Fallah-Bagher-Shaidaei et al., 2006), the relative chemical hardness (Δη) (De Proft and Geerlings, 2004), and the recently introduced AV1245 and AV_min (Casademont-Reig et al., 2018; Matito, 2016). This set of aromaticity descriptors are rooted in magnetic, reactivity and electronic criteria in order to take the multidimensional character of the aromaticity into account (Solà, 2017; Woller et al., 2016; Wu et al., 2013). For the evaluation of the aromaticity descriptors, we employed the long-range corrected functional CAM-B3LYP (Yanai et al., 2004), in combination with 6–311+G(d,p) (Hehre et al., 1986), due to the relevance of the delocalization error in these measures. It has been recently found that the amount of Hartree-Fock exchange in DFT functionals largely influences the value of the aromaticity descriptors in expanded porphyrin systems, although aromaticity trends remain (Casademont-Reig et al., 2018, 2020; Szczepanik et al., 2017). The global aromaticity indices, Λ and Δη, were computed on the basis of the isomerization scheme (Figure 2), which involves the comparison of the magnetic susceptibilities and hardness values of a methyl and methylene isomer of each hexaphyrin structure (Schleyer and Pühlhofer, 2002; Alonso et al., 2012).

Magnetic susceptibilities were evaluated with the Continuous Set of Gauge Transformations (CGST) method (Keith and Bader, 1993), whereas the hardness of the methyl and methylene isomers was computed as the HOMO and LUMO energy difference. The Gauge-Independent Atomic Orbital (GIAO) method (Wolinski et al., 1990) was applied to obtain the NICS values at three different positions: the geometrical center of the macrocycle’s heavy atoms and at 1Å above and below the ring center, together with their respective out-of-plane tensor components. The anisotropy of the induced current density (AICD) method (Herges and Geuenich, 2001; Geuenich et al., 2005) was applied to visualize the induced delocalization of π-electrons. In both the NICS and AICD methods, an external magnetic field is applied perpendicular to the molecular plane for expanded porphyrins with a Hückel or Möbius topology. Recently, extensive studies have been performed on the aromaticity of both unsubstituted and substituted expanded porphyrins by means of the ring current criterion and current density analysis (Fliegl et al., 2010, 2011; Fliegl and Sundholm, 2012; Valiev et al., 2013; Valiev et al., 2017; Valiev et al., 2018; Foroutan-Nejad et al., 2018; Foroutan-Nejad and Ghosh, 2018; Ghosh et al., 2018; Conradie et al., 2019). The electronic delocalization indices AV1245 and AV_min indices were computed with the ESI-3D program (Matito, 2006; Casademont-Reig et al., 2018) at the CAM-B3LYP/6-311+G(d,p) level of theory. The calculation of these electronic indices relies on Quantum Theory of Atoms in Molecules (QTAIM) atomic partition performed with the AIMAll software (Keith, 2017). AV1245 consists of the average of the 4-center multicenter index (MCI) values along the ring that keeps a positional relationship of 1, 2, 4, 5, whereas AV_min corresponds to the minimal absolute value of the aforementioned 4-center MCI values along the conjugation pathway (Matito, 2016; Casademont-Reig et al., 2018).

Next to the aromaticity descriptors, two geometrical descriptors, i.e., the torsional ring strain (Φ _p ) (Alonso et al., 2012) and the π-conjugation index (Π) (Stępień et al., 2011; Alonso et al., 2012), were computed to gain insight into the factors governing the relative energies of Hückel and Möbius conformations. Φ _p corresponds to the average dihedral angle between neighboring pyrrole rings, whereas Π measures the effectiveness of the overlap of neighboring p-orbitals. Π is positive for Hückel conformations and negative for Möbius structures. Macrocyclic aromaticity in porphyrinoids corresponds to values of Π above a threshold of 0.30 (Stępień et al., 2011).

Hyperpolarizabilities related to the nonlinear optical responses in the static and dynamic regimes were obtained at the coupled-perturbed Kohn-Sham and TD-DFT levels of theory, respectively, with the CAM-B3LYP functional and a Pople basis set 6-311G+(d,p). Different benchmark studies highlight the dependency of extended organic dyes’ NLO properties on the exchange-correlation functional and basis set. However, recent research has shown that range-separated hybrid functionals, such as CAM-B3LYP, are the best choice for the calculation of first order hyperpolarizabilities (Garza et al., 2015; Lescos et al., 2020). For the particular case of Hückel and Möbius π-conjugation topologies, Torrent-Sucarrat et al. showed that the hybrid M05-2X and the long-range corrected CAM-B3LYP provide a semi-quantitative description of the NLO properties at a reasonable computational cost (Torrent-Sucarrat et al., 2011; Torrent-Sucarrat et al., 2012; Torrent-Sucarrat et al., 2017b). In Torrent-Sucarrat’s work on designing topological switches based on meso-substituted [28]hexaphyrins with high NLO properties, both functionals together with the 6-31G basis set were employed (Torrent-Sucarrat et al., 2017b). While the good performance of M05-2X on the NLO properties can be ascribed to the large percentage of exact HF exchange (56%), the success of CAM-B3LYP is based on the introduction of a growing fraction of exact HF exchange at larger inter-electronic distances. Regarding the basis set dependency, split-valence double- or triple-ζ basis sets including one set of diffuse and polarization functions were shown to be generally sufficient to describe the dominant β tensor components and the depolarization ratios (Plaquet et al., 2008; de Wergifosse and Champagne, 2011; Castet et al., 2013). In particular, for the evaluation of the nonlinear optical properties of a Hückel-Möbius aromaticity switch of a closely related porphyrinoid A,D-di-para-benzi[28]hexaphyrin(1.1.1.1.1.1), Torrent-Sucarrat et al. showed that the increment of basis set from 6-31G to 6-311G(d) does not provoke a large variation of α and γ values (Torrent-Sucarrat et al., 2012). As such, the 6-311+G(d,p) basis set employed in our study is adequate to compute the β _HRS values of the Hückel and Möbius conformers of hexaphyrins. Indeed, in a very recent benchmark study on the influence of the amount of exact Hartree–Fock exchange included in the DFT functional on the magnitude of the static HRS responses, a similar basis set 6-311+G(d) was employed (Lescos et al., 2020).

In this study, we mainly focus on the first hyperpolarizability associated to the Hyper-Rayleigh scattering phenomenon, termed as β _HRS (Clays and Persoons, 1991; Hendrickx et al., 1998; Verbiest et al., 2009). The incoherent scattered light is at twice the optical frequency (2ω) of the incident laser pulse (ω) and its intensity is related to the first hyperpolarizability. The entire HRS intensity can be written as Eq. 3, when the incoherent scattered light is observed perpendicular to the laser’s propagation plane. β H R S ( − 2 ω ; ω , ω ) = 〈 β Z Z Z 2 〉 + 〈 β Z X X 2 〉 (3)

⟨ β Z Z Z 2 ⟩ and ⟨ β Z X X 2 ⟩ represent the orientational averages of β and describe the isotropic distribution of molecular orientations. The full descriptions of these tensor components are written in Eqs. 4, 5. ⟨ β Z Z Z 2 ⟩ = 1 7 ∑ i x , y , z β i i i 2 + 4 35 ∑ i ≠ j x , y , z β i i j 2 + 2 35 ∑ i ≠ j x , y , z β i i i β i j j + 4 35 ∑ i ≠ j x , y , z β j i i β i i j + 4 35 ∑ i ≠ j x , y , z β i i i β j j i + 1 35 ∑ i x , y , z β j i i 2 + 4 105 ∑ i ≠ j ≠ k x , y , z β i i j β j k k + 1 105 ∑ i ≠ j ≠ k x , y , z β j i i β j k k + 4 105 ∑ i ≠ j ≠ k x , y , z β i i j β k k j + 2 105 ∑ i ≠ j ≠ k x , y , z β i j k 2 + 4 105 ∑ i ≠ j ≠ k x , y , z β i j k β j i k (4) ⟨ β Z X X 2 ⟩ = 1 35 ∑ i x , y , z β i i i 2 + 4 105 ∑ i ≠ j x , y , z β i i i β i j j − 2 35 ∑ i ≠ j x , y , z β i i i β j j i + 8 105 ∑ i ≠ j x , y , z β i i j 2 + 3 35 ∑ i ≠ j x , y , z β i j j 2 − 2 35 ∑ i x , y , z β i j j β j i i + 1 35 ∑ i ≠ j ≠ k x , y , z β i j j β j k k − 2 105 ∑ i ≠ j ≠ k x , y , z β i i k β j j k − 2 105 ∑ i ≠ j ≠ k x , y , z β i i j β j k k + 2 35 ∑ i ≠ j ≠ k x , y , z β i j k 2 − 2 105 ∑ i ≠ j ≠ k x , y , z β i j k β j i k (5)

Lastly, an inverse design procedure was executed using the best-first search (BFS) algorithm (Pearl, 1984) as implemented in our in-house program CINDES (Teunissen, 2019; Teunissen, 2016). Both definitions for the contrast are considered for each substitution pattern. The fragment library consists of -NO₂, -CN, -F, -H, -CH₃, -OH and -NH₂. For all BFS optimizations, the computational approach and level of theory for the geometry optimization, frequency and NLO calculations are kept the same as aforementioned.

In this work, we discovered efficient nonlinear optical (NLO) switches with high β _HRS contrasts by applying the Best-First Search (BFS) algorithm. These innovative switches are based on meso-substituted hexaphyrins with different oxidation states and π-conjugation topologies and contain various electron-donating and electron-withdrawing substituents. Initially, we unveiled the structure-property relationships involving their aromaticity and NLO properties. As a proof of concept, we investigated the substituent effect on the geometry and (nonlinear) optical properties of different meso-substituted hexaphyrins previously reported in the literature. With the help of an inverse molecular design algorithm, the NLO contrasts (for which two definitions were taken into account, namely the ratio and the differences of the β _HRS responses of the ON and OFF states) of the most promising redox and topological switches were maximized by chemically modifying the ON and OFF states.

Unsubstituted [26] and [28] hexaphyrins prefer planar conformations driven by the low ring strain and high π-conjugation. A close relationship is found between the number of π-electrons, the topology and aromaticity as revealed by magnetic, electronic and reactivity aromaticity descriptors. While 26D, 26R and 28M are aromatic macrocycles, 28R is highly antiaromatic. The aromaticity is also reflected in their molecular orbital diagrams. Symmetry is the dominating factor in tweaking the β _HRS response of the parent structures. Other, contributing factors are the planarity and the ring strain of the macrocycle. The highest NLO contrasts are associated to redox or topological switches consisting of an OFF state (28R) and ON state (26R and 28M) with low and high β _HRS values, respectively. Changes in topology, aromaticity and symmetry are responsible for the high NLO contrasts.

The effect of meso-substitution on the β _HRS response of hexaphyrins was explored in detail. A proof of concept was included in which fully substituted macrocycles and 4 additional substitution patterns were evaluated, carefully selected from currently synthesizable hexaphyrins and a previous computational study. In general, meso-substituents distort the planarity of the macrocycle. As a consequence, the ring strain increases leading to an enhanced β _HRS response. The degree of the distortion depends on the number of introduced functional groups as well as their electronic character. Nonlinear optical properties (β _HRS ) can be enhanced by the introduction of different types of meso-substituents, especially by EDGs and push-pull combinations of substituents. A high difference in symmetry between the ON and OFF state of the switch significantly increases the NLO ratio. The results for the difference-based contrast are more subtle.

Multiple inverse molecular design procedures were performed, starting from three distinct substitution patterns, to design redox (28R → 26R) and topological (28R → 28M) switches with maximal NLO contrasts. Depending on the number, type and position of meso-substituents, ratio-based β _HRS contrasts up to 15 times larger for the redox switch and 2 times larger for the topological switch, compared to the unsubstituted versions, were discovered. For the difference-based contrast definition, the improvement is even more pronounced, reaching values 25 times larger for the redox switch and 6 times for the topological switch. Hence, the efficiency of the 28R → 26R redox switch can be enhanced upon meso-substitution to a much greater extent than for the 28R → 28M Hückel-Möbius switch. When the NLO response ratio is considered, both the optimal meso-substituted redox and topology switches show a preference for EDGs in combination with a centrosymmetric OFF state. Symmetry, more precisely the presence or absence of an inversion center, plays a pivotal role in tweaking the OFF state’s β _HRS as it can vary to a larger extent than the ON state’s response, from zero when in C _i symmetry to the order of 10³ a.u. Meso-substitution can induce profound changes in the symmetry resulting in a wide spectrum of different NLO responses. When centrosymmetry is far less likely due to the selected substitution pattern, e.g. A₂BC₂D, push-pull functionalizations take the upper hand. Nonetheless, the absence of low-response OFF states significantly reduces the ratio, despite higher ON state values. The difference-based contrast maximizations tend to almost always converge to push-pull configurations with OFF states showing a minimal NLO response. The main difference between the redox and topology switches is that EDGs generally outnumber EWGs for the former, while the reverse is true for the latter. On a final note, we would like to point out that replacing functional groups by substituents of the opposite nature can in some cases lead to a reversal of ON/OFF switching behaviour of the redox and especially topological switch, displaying more or less equivalent NLO contrasts.