Curtis E.
Colwell
a,
Tavis W.
Price
a,
Tim
Stauch
bc and
Ramesh
Jasti
*a
aDepartment of Chemistry & Biochemistry, Materials Science Institute, Knight Campus for Accelerating Scientific Impact, University of Oregon, Eugene, Oregon 97403, USA. E-mail: rjasti@uoregon.edu
bUniversity of Bremen, Institute for Physical and Theoretical Chemistry, Leobener Str. NW2, D-28359 Bremen, Germany
cMAPEX Center for Materials and Processes, University of Bremen, Bibliothekstraße 1, D-28359 Bremen, Germany
First published on 23rd March 2020
Strain has a unique and sometimes unpredictable impact on the properties and reactivity of molecules. To thoroughly describe strain in molecules, a computational tool that relates strain energy to reactivity by localizing and quantifying strain was developed. Strain energy is calculated local to every coordinate in the molecule and areas of higher strain are shown experimentally to be more reactive. Not only does this tool directly compare strain energy in parts of the same molecule, but it also computes total strain to give a full picture of molecular strain energy. It is freely available to the public on GitHub under the name StrainViz and much of the workflow is automated to simplify use for non-experts. Unique insight into the reactivity of curved aromatic molecules and strained alkyne bioorthogonal reagents is described within.
Total strain energies are commonly reported and the strain in specific parts of the molecule cannot be discerned. While the total strain energy does provide some information, it does not correlate perfectly to reactivity. For example, when the same amount of strain energy is spread over more atoms, the molecule is more stable than when it is concentrated in fewer atoms. If this is corrected for by dividing by the total atoms, non-participating atoms artificially lower the strain energy per atom determined. Local strain can sometimes be inferred in highly symmetric molecules, such as strain per phenylene in a cycloparaphenylene (CPP). However, unsymmetric molecules have unevenly distributed strain energy leading to locations of higher reactivity that may be unintuitive. Some alternative metrics have been devised to measure local strain. For example, the deplanarization of an aromatic ring, torsional angle in a biphenyl segment, and bond lengths can be compared to an unstrained comparative molecule.15,22 For non-planar π-systems, a measure of pyramidalization was developed that estimates relative strain in non-planar aromatics such as corannulene and fullerene.23,24 However, these measurements are not quantitative and, therefore, cannot be compared across molecules which limits their utility. Even when combining total strain calculation with these other metrics, it results in an incomplete depiction of molecular strain energy.
A method that determines strain both quantitatively and locally is quite useful. Therefore, a computational method was developed that identifies the quantity of strain energy local to every coordinate (bond, angle, and torsional angle) in a molecule. This strain visualization software is called StrainViz and has been made freely available. A similar method was previously reported for mechanochemistry where unstrained molecules are stretched and the tension that appears upon stretching is analyzed.25 Attempts were made to apply this method to macrocycles by comparing stretched and unstretched macrocycles as well as conformational changes that induce stretching,26,27 however, the inherent strain was not addressed. StrainViz can find this elusive strain energy. Our new method was evaluated to establish its accuracy using prior calculated strain energies and experimental reaction results from the literature. It is freely available on GitHub.28 The resulting computational method provides an interactive and insightful strain map. Knowledge of specific strain location facilitates and enhances synthetic efforts towards strained macrocycles by providing the exact and specific impact on strain of structural changes in a molecule. As is demonstrated herein with StrainViz, it is now possible to make inferences about the local properties and reactivity in strained molecules.
Estrain = (Emacrocycle + Ecap) − Elinear | (1) |
This is an example of how a homodesmotic reaction may be used for estimating strain energy. The homodesmotic reaction has been rigorously defined elsewhere.30
It is also possible to disassemble the molecule into an infinite polymer and compare the repeating unit energies in both the macrocycle and the unstrained polymer.31,32 If it were possible to connect these geometries by creating a trajectory between the strained and unstrained states, one could comment on how the energy of each atom changes to release strain energy as the trajectory proceeds. Ideally, the trajectory would begin in the optimized geometry of the strained molecule and descend to an unstrained infinite polymer (Fig. 2a). The macrocycle cannot be broken without changing the atomic environment and introducing additional strain into the analysis. A theoretical trajectory between the strained macrocycle and unstrained infinite polymer that isolates macrocyclic strain energy is therefore impossible.
It is, however, possible to fragment the molecule so that it may descend into an unstrained state without introducing new strain. By deleting certain atoms, the trajectory shown in Fig. 2b becomes possible and allows the initial geometry to share the location of its atoms (highlighted in the inset of Fig. 3) with the strained molecule while still relaxing to an unstrained state upon geometry optimization. The trajectory of each atom accurately represents the trajectory of atoms in the strained molecule to atoms in an unstrained state. This approximates an ideal strain energy determining experiment by averaging these trajectories for multiple fragments.
Fig. 3 Workflow for strain analysis. The coordinates in each molecule fragment relax to release strain energy that is quantified per coordinate (r: bond length, θ: angle, ϕ: torsional angle). |
In practice, a segment of the molecule is removed, such as a phenylene or ethylene, to create a fragment as shown in steps 1 and 2 of Fig. 3. The choice of fragments does not appear to significantly impact overall strain energy determination, however, it can impact strain distribution. Therefore, obtaining accurate results requires as many symmetrically created fragments as possible. For example, when analyzing CPPs, there should be as many fragments as there are phenylenes to be removed. Once a fragment is removed the ends are capped with hydrogen atoms. This requires the segment removed to be at least two atoms (e.g. ethylene) to accommodate replacement with hydrogens. These capping hydrogen atoms are optimized by freezing all atoms in the fragment that match the initial geometry. This ensures they do not add additional strain to the fragment when the fragment is optimized. All fragments created for the analyses in this paper are in the ESI† to remove any ambiguity.
The trajectory of each individual atom during this process follows the optimization algorithm given by the program used. In these studies, Gaussian09 (ref. 33) was used with the quasi-Newton rational function optimization (RFO) method that is the default for Gaussian03 due to it converging more smoothly than the newer direct inversion in the iterative subspace (DIIS) method.34 StrainViz can also be used with Orca, delivering similar results and being free for academic users.35,36 The energy of each atom is given by its relationship to other atoms via internal coordinates. The internal coordinates describe the distances and angles between atoms shown in a zoomed in image in Fig. 3. There are three coordinates that together describe the position of every atom relative to each other: the distance between two atoms, the angle between three atoms, and the torsional angle between four atoms. The optimization algorithm minimizes the energy of the geometry by interrogating these internal coordinates and adjusting them to release energy. The algorithm estimates a force for each coordinate (F) and, depending on step size, assigns a displacement (Δx). It is also possible to define “strain” as the force and “strain energy” as the total energy associated with that force. Multiplying the force by the displacement, as shown in eqn (2), identifies the change in energy (ΔEcoord est) each coordinate experiences in each step.
FΔx = ΔEcoord est | (2) |
The energy determined from this specific calculation is only an estimate. The algorithm overestimates the total energy released for each displacement due to the necessary use of redundant internal coordinates.37 Therefore, each step is scaled relative to the actual change in the single point energy calculated (ΔEstep actual in eqn (3)) at each step.
ΔEcoord actual ≈ ΔEcoord est(ΔEstep actual/ΔEstep est) | (3) |
After symmetrically fragmenting the molecule, optimizing the fragments, and analyzing the trajectory of the internal coordinates, this data must be displayed in a way that effectively communicates the information gathered. The effective color mapping scheme used for analyzing mechanical force was adopted.25 The bond strain energy associated is simple to display because there is a single value per bond and the bonds are colored accordingly. For the energy associated with the angle between three atoms, the energy is divided in half among the two contributing bonds. Finally, for the torsional angle between four atoms, the energy is split evenly among the three bonds connecting the four contributing atoms. These maps are produced for bond, angle, and torsional strain energy in each fragment. Then the energies per bond are averaged among every fragment containing that bond and a single map for each type of strain energy. Finally, the three types are summed and a total strain map is produced.
It is important to note that almost none of the above mentioned processes are done manually. A package of freely available scripts on GitHub automate Gaussian input file creation, job submission, and VMD script generation.28 Each analysis only requires the manual generation of an optimized geometry and appropriate fragments, StrainViz does the rest.
When molecule fragments are being analyzed, the fragment size must be judiciously chosen to reduce the impact of edge effects where the molecule is cut. It has been previously seen, in the strain-induced retro-Huisgen cycloaddition of triazoles,38 that when the edge atoms are connected to the triazole by coordinates (torsional angle across an ethyl or propyl group), results do not match expectations. Therefore, the program does not include any coordinates that contain the end capping atoms or the atoms attached to them. This trims away forces at the ends of the geometry that are most susceptible to these edge effects. By doing so, the results become more relevant regardless of fragment size, but limits how small the fragments may be made. Despite this consideration, the chosen fragment size does still have an impact on the accuracy. Within each fragment, the variability of strain energy measurement for each bond from the analysis also increases with decreasing fragment size. The four largest[8]CPP fragment sizes, shown in Fig. 4, all determine strain energies within 3% of their mean. These fragments also are internally consistent. Each individual strain energy determined for each bond is also within 3% of the mean. This consistency shows that when the fragments used are at least half the original molecule edge effects are minimal.
A variety of molecules with macrocyclic strain energy were used in this analysis (Fig. 5). Given that CPPs are well studied in this respect, CPPs having six to ten phenylenes were analyzed and compared to an analysis using homodesmotic reactions (Fig. 5a).17 Comparing these two analyses, we can see that the results are most similar at larger CPP sizes. This is consistent with the aforementioned accuracy of the StrainViz analysis where larger CPP fragments result in more accurate strain energy determinations. Analyzing Itami's carbon nanobelt resulted in even better matching with previous efforts (Fig. 5b).41 The high accuracy may be owed to the fragment optimization trajectory quality. See ESI† for further comments. A recently synthesized highly strained small cyclophane from the Bodwell group42 and [2.2]paracyclophane43 were also analyzed and confirm that StrainViz is consistent with prior computational efforts (Fig. 5c).
Fig. 5 Literature examples of strain energy determinations compared to StrainViz analysis. (a) CPP strain determined by homodesmotic reactions.17 (b) Carbon nanobelt strain energy extrapolated from increasing size belts.41 (c) [2.2]paracyclophane strain and Bodwell's more strained analogue determined by isodesmotic reaction B3LYP/6-31G(d,p)43 and M06-2X/Def2TZVP42 respectively. Please see the ESI† for links to the 3D models. |
The literature report of [2.2]paracyclophane does attempt to quantify the strain energy present in the phenylene and ethylene segments.43 This was done in a similar manner to our method. The molecule was broken up and the strain energy in each fragment was determined by comparing single point energies of the strained and unstrained states. Their analysis, however, found different amounts of strain than their homodesmotic reaction; 10.2 kcal mol−1 per phenylene and 5.6 kcal mol−1 per ethylene summing to 31.6 kcal mol−1, but 30.8 kcal mol−1 from a homodesmotic reaction using the ωB97X-D functional. This begs the question which analysis is more accurate. With our method, the local strain adds up to the total strain by definition.
By corroborating results found in the literature, we establish that StrainViz determines total strain energies that are reasonable. This shows that generating a map of local strain does not compromise the total strain analysis quality. More importantly, we see exactly where in the structure the strain is distributed. We hope to confirm that local strain is more instructive than total strain for reactivity.
Fig. 6 Torsional angles and total strain of [6]CPP and m[6]CPP. Changing connectivity from para to meta decreases total strain, but increases local strain energy. Please see the ESI† for links to the 3D models. |
Analysis using StrainViz in Fig. 6 shows that strain is concentrated across from the meta phenylene. Changing a phenylene in a CPP from being para to meta connected relieves strain at that end of the molecule, but adds strain at the opposite end. If this program provides a meaningful molecular strain analysis, then this high strain area should react faster in a strain relieving reaction. For example, bromination of a meta-CPP should occur exactly across from the meta phenylene where the majority of the strain is located. Indeed, upon bromination of a m[6]CPP, bromination occurs exactly where predicted by our calculations (Scheme 1).
Even more striking from the analysis is that despite meta-CPPs being less strained in total, they should be more reactive due to a higher amount of local strain relative to a CPP where strain is spread equally over the molecule. This phenomenon is seen clearly in the aforementioned publications by Yamago. In the case of either bromination10 or C–C bond activation by platinum,11 a reaction at one phenylene is followed by a faster second reaction. This is not consistent with the total quantity of strain present in each reacting molecule. A homodesmotic reaction of the starting CPP and singly brominated CPP shows that the second has much less strain energy. However, when analyzed using StrainViz as shown in Fig. 7, it is clear that the singly brominated intermediate has more strain across from the first site of bromination and that the molecule has been activated to brominate the second time at a faster rate.
Fig. 7 Strain release during bromination of [6]CPP. First bromination activates molecule to be more reactive in the second step. Please see the ESI† for links to the 3D models. |
Despite these specific examples of reactivity correlating extremely well to the strain energy present at certain locations in the molecule, it is important to note that the specific reaction taking place is being aided by relief of strain to effect the transformation. This is why it is not the bond with the highest determined strain energy that is broken during the reaction, but that the chemical reactivity occurring is preferred when in proximity to higher strain energy.
Fig. 8 Möbius molecules have more strain due to an internal twist when compared to a non-Möbius belt. (a) The Möbius molecule synthesized in the Tanaka group is more strained than the non-Möbius despite being a larger size. (b) In a symmetric Möbius molecule, strain is centered at the twist entry point. Please see the ESI† for links to the 3D models. |
In addition to macrocyclic molecules, multimacrocyclic molecules can be analyzed as long as each fragment fully releases all strain present. The Yamago group has reported a highly symmetric nanoball with multiple macrocyclic connections (Fig. 9).50 Analyzing a single panel shows that strain is spread relatively evenly around the periphery aside from some anisotropy induced by the C2 symmetry of this lowest energy conformation. In the ball, however, it appears that there is more strain at the corners as opposed to the edges. This indicates that the three additional macrocycles in the ball add strain where they attach at the corners of a panel. As they do not apply force directly in the direction of any edges, the force is split between the edges and concentrates at the corners. This unique multimacrocyclic strain contribution is easily apparent using this analysis.
Fig. 9 Strain energy present in Yamago's nanoball. More strain at the ball corners relative to the edge. Please see the ESI† for links to the 3D models. |
While this analysis was designed for analyzing curved aromatic molecules, it applies widely in the analysis of strained molecules. For example, strained hydrocarbons play a very important role as bioorthogonal reagents. Specifically, cyclooctynes and trans-cyclooctenes are used as reactive reagents for copper-free click reactivity in biological media. A previous Houk group analysis51 noted that a strain energy based analysis fails to accurately predict molecular reactivity and instead used a distortion/interaction model52 to accurately predict reactivity. However, our novel analysis can correctly order the reactivity of these strained reagents using a local strain analysis. The analysis of cyclooctyne in Fig. 10 reveals 13.7 kcal mol−1 of strain whereas trans-cyclooctene has 17.4 kcal mol−1 both have nearly the same proportion (42% and 41%) located at the reactive site. This is in agreement with the relative rate of reaction with a tetrazine of 30 and 13000 M−1 s−1 respectively. Furthermore, when comparing the more reactive trans-bicyclo[6.1.0]nonene to trans-cyclooctene in Fig. 10, the total strain energy increases only slightly to 18.9 kcal mol−1, however, the reaction rate increases 160 fold.53 By increasing the macrocycle rigidity, the strain is shifted to the reactive alkene. The trans-cyclooctene has 7.2 kcal mol−1 of strain energy located in the alkene, whereas trans-bicyclo[6.1.0]nonene has 9.3 kcal mol−1. This increases the reactivity more than would be predicted by total strain energy. Again, the StrainViz analysis provides the missing information that previous strain analyses lack. While this analysis does not provide accurate prediction of rates as the distortion/interaction model of the Houk group does,52 it may be put to good use in informing the design of new strained bioorthogonal reagents by evaluating strain local to the reactive site.
Fig. 10 Strain energy in copper-free click reagents. Trans-cyclooctene is more strained than cyclooctyne, but has similar strain distribution. Increasing rigidity by addition of a cyclopropyl fusion increases strain energy and shifts it to the reactive site. Please see the ESI† for links to the 3D models. |
Footnote |
† Electronic supplementary information (ESI) available: Additional details on computations and synthesis (PDF). Raw data for all calculations (ZIP). See DOI: 10.1039/d0sc00629g |
This journal is © The Royal Society of Chemistry 2020 |