Multi-replica biased sampling for photoswitchable pi-conjugated polymers
Mariagrazia Fortino, Concetta Cozza, Massimiliano Bonomi, Adriana Pietropaolo
MMulti-replica biased sampling for photoisomerization processes in conjugated polymers
Mariagrazia Fortino , Concetta Cozza , Massimiliano Bonomi , Adriana Pietropaolo AFFILIATIONS Dipartimento di Scienze della Salute, Università di Catanzaro, Viale Europa, 88100 Catanzaro, Italy. Structural Bioinformatics Unit, Department of Structural Biology and Chemistry; CNRS UMR 3528; C3BI, CNRS USR 3756; Institut Pasteur, Paris, France. ORCID: 0000-0002-7321-0004. Dipartimento di Scienze della Salute, Università di Catanzaro, Viale Europa, 88100 Catanzaro, Italy. ORCID: 0000-0003-0955-2058 * Authors to whom correspondence should be addressed: [email protected], [email protected]
ABSTRACT
In recent years, conjugated polymers are attracting considerable interest in view of their light-dependent torsional reorganization around the p -conjugated backbone, which determines peculiar light-emitting properties. Motivated by the interest in designing conjugated polymers with tunable photoisomerization pathways, we devised a computational framework to enhance the sampling of the polymer conformational space and at the same time estimate ground to excited-state free-energy differences. This scheme is based on a combination of Hamiltonian Replica Exchange (REM), Parallel Bias metadynamics, and free-energy perturbation theory. In our scheme, each REM replica samples different intermediate states connecting the ground to the first two excited states, which are characterized by TD-DFT simulations at the B3LYP/6-31G* level of theory. We applied the method on a 5-mer of poly(9,9-dioctylfluoren-2,7-diyl) and compared the results with the emission energies measured experimentally, showing a quantitative agreement with the prediction provided by our simulation framework. I. INTRODUCTION
Conjugated polymers are a class of organic frameworks that are widely used in an extensive range of electronic applications, thanks to their promising optical and electronic properties. These polymers offer an excellent alternative to inorganic materials since they are moderately low-cost. For these reasons, conjugated polymers have found remarkable applications in devices such as organic photovoltaics (OPV), organic light-emitting diodes (OLED), organic field-effect transistors (OFET) and a variety of sensors.
A peculiarity of these systems is the photo-induced torsional reorganization around the conjugated backbone, due to a trans-cis isomerization, which allows a specific molecular response after light-irradiation. everal computational efforts have been performed in the last decades in order to optimize their structural and optical properties, their electronic structure and their charge-transport character. For instance, excited-state ab initio dynamics within non-adiabatic schemes or surface hopping algorithms were used to predict the photo-induced dynamics in such materials in order to gain accurate description of excited-states formation determining a specific non-radiative decay. In 2012,
Clark et al. have shown through non-adiabatic excited-state dynamic simulations that the photo-induced torsion around the backbone of conjugated polymers in some conditions may occur in sub-100 fs timescale. A flattening mechanism was predicted to occur during the decay, showing the conversion of the electronic potential energy into torsional kinetic energy. More recently, ab initio excited-state molecular dynamic simulations have been used to identify the torsional motions generating the signal observed in Femtosecond Stimulated Raman spectra of conjugated polymers. Specifically, the spectral evolution reflecting the excited state relaxation towards a planar conformation was investigated by considering the first excited state of two different quaterthiophene derivatives, relaxing through stretching modes coupled with dihedral backbone angle. However, non-adiabatic excited state molecular dynamics simulations often require excited-state quantum simulations performed at every time step on an ensemble of trajectories, each propagated for a picosecond or more. Consequently, although this field has matured significantly with remarkable improvements in the simulation of very large molecular systems, it often requires to make compromises in order to find the best balance between precision and numerical cost. Recently, we have developed a simple approach that is able to estimate the free-energy gap between ground and excited states, while at the same time enhancing the sampling of different conformers. Within this simulation framework, several sets of independent free-energy simulations are performed to enable an exhaustive sampling of the conformational space, and free-energy perturbation theory (FEP) is then used to provide an accurate estimate of the free-energy gap between ground and first excited state. Herein, we extend the method biasing multiple torsions of polyfluorene derivatives substituted in the 9-position by dioctyl chains (Scheme 1), whose substitution is usually employed to improve the solubility in organic solvents. The Parallel Bias metadynamics approach combined with Hamiltonian Replica Exchange (REM) are used to increase the efficiency of our scheme and, in combination with FEP, to characterize the free-energy landscapes of the first singlet S ground state and the lowest singlet excited states S and S as well as the transitions between these states. Scheme 1.
Photo-induced trans-cis isomerization of poly(9,9-dioctylfluoren-2,7-diyl).
II. THEORETICAL FRAMEWORK
Derivation of the torsional potential
The potential energy scans (PES) as a function of the dihedral rotation were performed within the Density Functional Theory (DFT) framework for the torsions belonging of a 5-mer polyfluorene (Scheme 1) in the ground and in the first two singlet excited states. The Becke’s three-parameter hybrid functional (B3) with the Lee, Yang, and Parr (LYP) expression for the nonlocal correlation, B3LYP were used with the 6-31G* basis set using the Gaussian16 package. The torsional conformers were optimized varying the ! torsion from 0 to 180 degrees in the S , S and S states. While in the ground state the " !" values of the PES of the consecutive torsions does not show appreciable differences, in the excited states the external ( ! ) and the internal torsion ( ! $ ) markedly vary (Figure 1). In particular, the PES of the second torsion ( ! $ ) shows lower energy values for coplanar conformers and a higher potential energy barrier at 90 degree that increases owing to the loss of conjugation, with respect to the PES of the first torsion ( ! ) . As expected, a dihedral rotation triggered by light irradiation is most likely to start in the external polymer section, which acts as a seed facilitating the isomerization of the consecutive main-chain dihedrals, according to the “sergeants and soldiers” principle. Coplanar states are stabilized through light irradiation owing to the lowering of the potential energy barriers at 0 degree in S and S excited states and to the increase of the potential energy barrier at 90 degree (11 kJ/mol for S , 17 kJ/mol for S and 21 kJ/mol for S concerning the ! dihedral and 11 kJ/mol for S , 30 kJ/mol for S , 24 kJ/mol for S concerning the j dihedral). Consequently, we derived two sets of torsional potentials for the external and the internal torsions of the polymer. Non-linear least-square interpolations were performed by minimizing the square difference between the molecular mechanics energies ( " "" ) and the ab initio energies ( " !" ) for the $ points on the PES of the electronic state i , reported in Eq. 1. % = ∑ (" !"% (! ' ) − " ""% (! ' )) $(') (1) where: " ""% (!) = " *+,-"" + , . ! (!) + ∑ - / (1 + cos (2! − ! )) (2) and: , . ! (!) = ∑ 3 ! (cos (! − 180)) (3) in which " *+,-"" represents the force-field potential CHARMM CGenFF and , . ! (!) represents the torsional potential in a given electronic state i related to the torsions for which each PES is calculated ( j is defined as from Scheme 1). The term - / and ! represent the force constants and the phase of a cosine series describing the concurrent torsions (1235, 5326, 6234 from Scheme 1). The values of " "" were calculated from an energy minimization of the optimized coordinates of the QM PES. We also tested the fitting parameters through the torsional free-energy profiles calculated at 200 K, where the free energy approaches the potential energy (Figure 1). The - / and ! terms of the concurrent torsions were calculated with iterative fits starting from - / =0 and using as initial guess the fitting parameters from the QM PES interpolation for the ! torsion in a given electronic state. The - / and ! terms found for the ground-state combined with the fitting parameters of the QM PES in a given electronic state, well reproduced the potential energy along the ! torsion (Figure 1). The fitting parameters describing the , . ! (!) potential are reported in Table 1, the parameters describing the concurrent torsions are reported in Table 2. The ESP charges were calculated for each optimized electronic state of the 5-mer oligofluorene at the B3LYP/6-31G* and are reported in Tables S1 of Supplementary Material. We also calculated the ESP charges at the MP2 and CIS level of theory. Those are consistent with the ones derived at the B3LYP/6-31G* level of theory and are reported in Table S2 of the Supplementary Material. Table 1 . The coefficients ! (kJ/mol) of the Ryckaert-Bellemans function described in Eq. 3 and used in GROMACS, obtained through a non-linear least-square interpolation of the PES values obtained at the B3LYP/6-31G* level. C0 C1 C2 C3 C4 C5 S ( j ) S ( j ) S ( j ) S ( j ) S ( j ) able 2. Torsional parameters for the concurrent torsions and described in Eq. 2. - / (kJ/mol) j m - / (kJ/mol) j m -6.40 180 4 -2.20 0.0 2 2.20 180 2 -8.20 0.0 4 Figure 1.
Fluorene-fluorene torsional potential calculated from energy-minimized conformations at MM level and from free energy calculated with Well-tempered metadynamics simulations at 200 K. The total energy calculated at B3LYP/6-31G* level of theory is also shown. a) Calculations related to the torsion ( ! ) of the 5-mer polyfluorene in the S state, b) calculations related to the torsion ( ! $ ) of the 5-mer polyfluorene in the S state, c) calculations related to the torsion ( ! ) of the 5-mer polyfluorene in the S state, d) calculations related to the torsion ( ! $ ) of the 5-mer polyfluorene in the S state, e) calculations related to the torsion ( ! ) of the 5-mer polyfluorene in the S state, f) calculations related to the torsion ( ! $ ) of the 5-mer polyfluorene in the S state. amiltonian Replica Exchange combined with Free-energy perturbation The free-energy as a function of the fluorene-fluorene dihedral angles ! % is defined as follows: . ! (! % ) = − log 9 ∫ 5(7 ! ! (9)); " <9∫ ; " <9 : (4) In this equation, ,(;) represents the interatomic potential (the force field), ; the atomic coordinates, and < represents the inverse of = = > , where = = is the Boltzmann constant and > the temperature of the system. We now add to the interatomic potential the torsional potential reported , . ! (! % ) defined in Eq.3. Hence, the free-energy surface for a given electronic state . ! (! % ) becomes: . ! (! % ) = − log ? ∫ 5(7 ! ! (9)); " <9∫ ; " <9 @ (5) In order to evaluate the free-energy difference between the S , S and S electronic states, FEP was used in combination with REM, following the approach introduced in Ref. . Multiple intermediate windows between the S , S and S electronic states were introduced through two parameters A and A $ . The former is used to reach S from S ( A : 0 → 1 ; A $ = 0 ); the latter to reach S from S ( A =1 ; A $ : 0 → 1 ) (Scheme 2). Scheme 2.
Scheme of the simulation framework herein proposed based on Hamiltonian Replica Exchange with multiple replicas differing from the torsional potentials through A and A $ parameters. Each torsion is biased within the parallel bias metadynamics approach. The external potential acting on each replica is defined as: , > (!) = ∑ (1 − A ), . / (! / ) ?/) + A (1 − A $ ), . (! / ) + A $ , . (! / ) (6) Where represents the number of torsions belonging to a 5-mer polymer. In order to obtain an adequate overlap between the energy distributions of adjacent windows, a step sizes Dl equal to 0.04 and Dl equal to 0.1 were used for the S -S and S -S transitions, respectively. In this REM framework, $ non-interacting replicas of the system, each identified by a given pair A = (A , A $ ) of control parameters, were simulated in parallel. At fixed time intervals, a swap between two configurations belonging to neighboring replicas was attempted and accepted using a Metropolis criterion based on the overlap between energy distributions of the two replicas. In this way, instead of multiple independent FEP runs, the simulation can be executed in a single multi-replica run. Within his framework, the ground state replica with higher coplanar free-energy barriers can be exchanged with excited state replica having lower coplanar free-energy barriers, thus enhancing sampling. Furthermore, the overlap between energy distributions of neighboring replicas, which is ensured to obtain efficient diffusion in the A = (A , A $ ) space, guarantees also an accurate calculation of the free-energy difference between two neighboring replicas using the FEP approach in Eq. 7: Δ6 (@ ! AB@) = − log〈F ) 〉 > ! (7) Finally, the free-energy gap between the S and S and S and S electronic states can be estimated as the sum over the free-energy differences between neighboring replicas using Eq. 8: Δ6 . ! →. !(0 = ∑ Δ6 (> ! AB>)% (8)
Enhanced sampling with Parallel Bias metadynamics.
Parallel Bias metadynamics (PBMetaD) was used to enhance sampling in each REM replica. The four main-chain dihedrals belonging to the polymer ( ! , ! $ , ! E , ! ? ) were used as Collective Variables (CVs). Multiple mono-dimensional metadynamics bias potentials , F (! % ; J) were used to construct the PBMetaD bias potential , G= (! , ! $ , ! E , ! ? ; J) defined by: , G= (! , ! $ , ! E , ! ? ; J) = − logK∑ F (7 ! ;I))?%) L (9) The individual metadynamics potential , F (! % ; J) is built adaptively during the course of the simulation by depositing Gaussian functions along the system trajectory in the CVs space: , F (! % ; J) = ∫ NJ JKL
O(J J ) ∙ exp T− M7 ! (9)87 ! N9OI PQR $S !1 U (10) where J is simulation time, V % the Gaussian width of the i- th CV, and O the deposition rate of the bias potential. Typically, Gaussians are added to the simulation with a discrete and constant deposition stride W . Therefore, the deposition rate O is expressed as the ratio between the Gaussian height X and the deposition stride W . Similarly to well-tempered metadynamics, in PBMetaD the Gaussian height decreases with simulation time as: X(J) = X ∙ exp Y− C (7 ! ;I)T UV Z ⋅ ; " ∑ ; " (11) where X is the initial Gaussian height and an input parameter with the dimension of a temperature can be used to limit the exploration to relevant regions of the CV space and thus void visiting excessively high free-energy areas. This parameter is often expressed in terms of the so-called biasfactor: \]^_J = VAUVV (12) At convergence, one-dimensional free-energy profiles !" ! X can be reconstructed directly from the bias potentials , F (! % ; J) as in standard well-tempered metadynamics simulations: , F (! % ; $ → ∞ ) = − !""+!" ∙ ! (! % ) + , (13) while multi-dimensional free-energy surfaces can be reconstructed by reweighting. The presence of the PBMetaD bias potential , G= (! , ! $ , ! E , ! ? ; J) was accounted for when calculating the swap acceptance probability in the REM scheme as well as when computing the FEP ensemble average in Eq. 7. Simulation Details
The initial 5-mer of poly(9,9-dioctylfluoren-2,7-diyl) structures were built by selecting the polymeric chains from the simulated system reported in Ref. In order to derive the force field in the ground and excited states, potential energy scans as a function of the dihedral rotation of the 5-mer polyfluorene were carried out using the Density Functional Theory (DFT) framework with the B3LYP functionals and the 6-31G* basis set within the Gaussian16 package. The torsional conformers were optimized varying the torsion from 0 to 180 degrees in the S and S states. Subsequently, the initial coordinates were energy-minimized in vacuo using the Steepest Descent Algorithm and equilibrated in vacuo at 200 K and at 300 K using the ab initio derived torsional potential combined with the CHARMM based General Force Field. The equilibrated conformation at 200 K was used as starting coordinates for reconstructing the free-energy profile at 200 K. Specifically, two well-tempered metadynamics runs were carried out using the first (external) or the second (inner) polyfluorene dihedrals. Gaussians with initial height equal to 1.2 kJ/mol were deposited with a stride of 1 ps and biasfactor equal to 25, with a temperature of 200 K. The equilibrated conformation at 300 K was used as starting coordinates for the REM simulations. 36 replicas were used: 26 to cover the interval S , S and 10 for the interval S , S Each replica was simulated in the NVT ensemble at 300 K enforced by the velocity rescaling thermostat. swap between two neighboring replicas was attempted every 100 MD steps and accepted using a Metropolis criterion. For the PBMetaD approach, Gaussians with initial equal to 1.2 kJ/mol were deposited with a stride of 1 ps and biasfactor equal to 40. In all simulations, the time step was set at 2 fs. The LINCS algorithm was used to fix all bonds lengths. For the Lennard-Jones and electrostatic interactions, a cutoff of 2.0 nm was used; the Particle Mesh Ewald method was used to calculate electrostatic interactions. All simulations were performed with GROMACS version 2020.2 equipped with PLUMED 2.7.0. III. FREE-ENERGY PATHWAYS IN THE GROUND AND EXCITED STATES
The free-energy profiles for the S - S and S - S excitations of the investigated polyfluorene pentamer are shown in Figure 2 as a function of the external dihedral angle j . All the free-energy profiles for each dihedral angle are reported in Figure S1 of the Supplementary Material. In the ground state, no significant differences have been observed between free-energy profiles as a function of different dihedral angles. Small variations appear in the free-energy profiles of the excited states. The S state shows different orthogonal free-energy barrier at 90 degrees on the external and internal torsions (0.22 eV for ! , ! $ ,0.33 eV for ! E and 0.24 for ! ? ), whereas the S state shows variations in the coplanar free-energy barrier (0.03 eV for ! or ! ? ; 0.09 eV for ! $ and 0.11 eV for ! E ). The calculated S - S free-energy gap estimated by FEP is equal to 2.946 eV with a statistical error of 310 -4 eV, while the calculated S - S free energy gap is predicted as 0.442 eV with an error of 210 -5 eV. Notably, the calculated S - S free-energy gap shows a good agreement with the experimental photoluminescence spectra centered at 2.88 eV reported for a 5-mer of poly(9,9-dioctylfluoren-2,7-diyl). Figure 2. A.
Free-energy profiles for the S - S excitation, computed using the Parallel Bias metadynamics approach combined with Hamiltonian Replica Exchange, as a function of the outer dihedral angle ! . B. Free-energy profiles for the S - S excitation as a function of the outer dihedral angle ! . he free-energy profiles of S , S and S states and the relative gaps are highlighted in Figure 3A. At variance with an unsubstituted 5-mer oligofluorene in which the free-energy landscape is shaped through four degenerate free-energy minima, the dioctyl chains stabilize the dihedral values of 144 degree typical of a 5 helix, instead of the values of 72 degree featuring a 5 helix. The more favorable packing of a 5 helix instead of a 5 helix was previously suggested from electron diffraction patters combined with geometry optimizations at RHF/6-31G(d,p) level of theory. The 5 helix is further stabilized in the S and S excited states where a higher orthogonal free-energy barrier at 90 degrees is observed for the free-energy profiles reconstructed along the first ! dihedral of S and S states (0.23 eV and 0.28 eV, respectively) compared to the free-energy barrier of 0.18 eV in the one of S . A switch region is detected at values of angles of 44 degree, where a reduced free-energy coplanar barrier is observed for the ! transition towards a coplanar state ( ~ ~ S and S states are shifted towards higher values approaching 155 degree, whereas the dihedral angles of the switch free-energy basin are shifted towards lesser values approaching 27 degrees (Figure 3A), favoring in this way a molecular switch from positive to coplanar dihedrals reaching negative twist dihedrals. These results are in agreement with the typical twisted-coplanar transition of biphenyl-based rotors, herein represented in the molecular conformations sketched in Figure 3B for S , S and S electronic states. Figure 3. A.
Calculated free-energy profiles for S , S and S states. B. Representation of the molecular conformational pathway occurring during light irradiation, starting from a bent ground state and crossing the flat conformation of the S and S excited states. IV. CONCLUSIONS
In summary, we have introduced a simulation scheme that combines Hamiltonian Replica Exchange, Parallel Bias metadynamics and the free-energy perturbation theory to enhance sampling of the conformational space of a polymer and at the same time to estimate the free-energy gaps between ground and excited states. The method was tested on a 5-mer of poly(9,9-dioctylfluoren-2,7-diyl) nd succeeded in predicting the energetics and structural properties, such as the ground to excited states emission energies as well as the highest stability of a 5 helix conformation compared to a 5 helix one. The presented simulation framework is general and can be used to predict large-scale rearrangements occurring through light-matter interactions. SUPPLEMENTARY MATERIAL
See the supplementary material for ESP charges (Tables S1 and S2) and for the simulated free-energy profiles for each dihedral angles.
DATA AVAILABILITY STATEMENT
ACKNOWLEDGMENTS
The ISCRA supercomputing initiative is acknowledged for computational time provided. This research was funded by Italian Ministry of Research under the program PRIN2017 no. 2017WBZFHL_003 (AP).
References
1. Y. Tsuji, Y. Morisaki, Y. Chujo, Polymer Journal , 203 (2017). 2. R.M. Pankow, B. C. Thompson, Polymer , 122874 (2020). 3. R. Noriega, J. Rivnav, K. Vandewall, F.P.V. Koch, N. Stingelin, P. Smith, M.F. Toney, A. Salleo, Nat Mater.
12 (11) , 1038 (2013). 4. A. Facchetti, Chem. Mater.
23 (3) , 733 (2011). 5. A. Grimsdale, K.L. Chan, R.E. Martin, P.G. Jokisz, A. B. Holmes, Chem. Rev.
109 (3) , 897 (2009). 6. J. Clark, G. Lanzani, Nature Photon. , 438 (2010).
7. B. K. Yap, R. D. Xia, M. Campoy-Quiles, P. N. Stavrinou, D. D. C. Bradley, Nature Mater. , 376 (2008). 8. T. Virgili, D. Marinotto, C. Manzoni, G. Cerullo, G. Lanzani, Phys. Rev. Lett. , 117402 (2005). 9. K. C. Vishnubhatla, J. Clark, G. Lanzani, R. Ramponi, R. Osellame, T. Virgili, Appl. Phys. Lett. , 041123 (2009). 10. D. Kabra, L. P. Lu, M. H. Song, H. J. Snaith, R. H. Friend, Adv. Mater. , 3194 (2010). 11. I. A. Barlow, T. Kreouzis, T, D. G., Lidzey, Appl. Phys. Lett. , 243301 (2009).
12. C. R. McNeill, N. C. Greenham, Adv. Mater. , 3840 (2009). 13. Y. Wang, T. Sakamoto, T. Nakano, Chem. Commun., , 1871(2012). 14. Y.Wang, T. Harada, L. Q. Phuong, Y. Kanemitsu, T. Nakano, Macromolecules
51 (17), , 5509 (2013). 6. A. Pietropaolo, S. Tang, F. M. Raymo, Nanoscale, , 368 (2010). 18. L. Angiolini, T. Benelli, L. Giorgini, F. Mauriello, E. Salatelli, R. Bozio, A. Daurù, D.Pedron, Eur. Polym. J. , 3550 (2007). 19. V.Marcon, N. van der Vegt, G. Wegner, G. Raos, J. Phys. Chem. B , 5253 (2006). 20. W. Chunwaschirasiri, B. Tanto, D. L. Huber, M. J. Winokur, Phys. Rev. Lett. , 107402 (2005). 21. X. K. Chen, Y. Tsuchiya, Y. Ishikawa, C. Zhong, C. Adachi, J. L. Bredàs, Adv. Mater. , 1702767 (2017). 22. S. Kilina, E. R. Batista, P. Yang, S. Tretiak, A. Saxena, R. L. Martin, D. L. Smith, ACS Nano , 1381 (2008). 23. T. Körzdörfer, J. L. Bredàs, Acc. Chem. Res. , 3284 (2014). 24. E. Fron, A. Deres, S. Rocha, G. Zhou, K. Müllen, F. C. De Schryver, M. Sliwa, H. Uji-i, J. Hofkens, T. Vosch, J. Phys. Chem. B , 1277 (2010). 25. J. S. Kim, L. Lu, P. Sreearunothai, A. Seeley, K. H. Yim, A. Petrozza, C. E. Murphy, D. Beljonne, J. Cornil, R. H. Friend, J. Am. Chem. Soc. , 13120 (2008). 26. J. S. Kim, L. Lu, P. Sreearunothai, A. Seeley, K. H. Yim, A. Petrozza, C. E. Murphy, D. Beljonne, J. Cornil, R. H., J. Am. Chem. Soc. , 13120 (2008). 27. I. Tavernelli, U.F. Röhrig, U. Rothlisberger, Mol. Phys. , 963 (2005). 28. B. F. E. Curchod, T. J. Martìnez, Chem. Rev., , 3305 (2018). 29. N. L. Doltsinis, D. Marx, Phys. Rev. Lett. , 166402 (2002). 30. N.L. Doltsinis, D. Marx, J. Theor. Comput. Chem , 319 (2002). 31. D. Hollas, L. Šisťík, E. G. Hohenstein, T. J. Martínez, P. Slavícěk, J. Chem. Theory Comput., , 339 (2018). 32. E. Tapavicza, G. D. Bellchambers, J. C. Vincent, F. Furche, Phys. Chem. Chem. Phys. , 18336 (2013). 33. H. Song, S. A. Fischer, Y. Zhang, C. J. Cramer, S. Mukamel, N. Govind, S. Tretiak, J. Chem. Theory Comput
16 (10) , 6418 (2020). 34. J. C. Tully, R. K. Preston, J. Chem. Phys , 562 (1971). 35. J. C. Tully, J. Chem. Phys. , 1061 (1990). 36. J. R. Schmidt, P. V. Parandekar, J. C. Tully, J. Chem. Phys. , 044104 (2008). 37. E. Tapavicza, I. Tavernelli, U. Rothlisberger, Phys. Rev. Lett. , 023001(2007). 38. T. Vreven, F. Bernardi, M. Garavelli, M. Olivucci, M. A. Robb, H. B. Schlegel, J. Am. Chem. Soc. , 12687 (1997). 39. T. R. Nelson, A. J. White, J. A. Bjorgaard, A. E. Sifain, Y. Zhang, B. Nebgen, S. Fernandez-Alberti, D. Mozyrsky, A.E. Roitberg, S. Tretiak, Chem. Rev. , 2215 (2020). 40. T. Nelson, S. Fernandez-Alberti, V. Chernyak, A.E. Roitberg, S. Tretiak, J. Phys. Chem. B , 5402 (2011). 41. T. Nelson, S. Fernandez-Alberti, A. E. Roitberg, S. Tretiak, Acc. Chem. Res. , 1155 (2014). 42. J.Clark, T. Nelson, S. Tretiak, G. Cirmi, G. Lanzani, Nat. Phys., , 225 (2012). 43. S. Dasgupta, J. M. Herbert, J. Phys. Chem. A , 6356 (2020). 44. I. Conti, G. Cerullo, A. Nenov, M. Garavelli, J. Am. Chem. Soc., , 16117 (2020). 5. E. V. Stolyarov, A. J. White, D. Mozyrsky, J. Chem. Phys, , 074116 (2020). 46. C. Cozza, M. Bonomi, A. Pietropaolo, J. Chem. Theory Comput., , 5441 (2018). 47. W.L Jorgensen, L.L. Thomas, J Chem Theory Comput. , 869 (2008). 48. B. W. J. Irwin, J. Chem. Theory Comput. , 5062 (2015). 50. Y. Sugita, Y. Okamoto. Chem. Phys. Lett. , 141–151 (1999). 51. P.J. Stephens, F.J.Devlin. C.F. Chabalowski, M.J. Frisch, J.Phys.Chem. , 11623 (1994). 52. Gaussian 16, Revision C.01, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2016 53. M.M. Green, M.P. Reidy, R.D. Johnson, G. Darling, D.J. O'Leary, G. Willson, J. Am. Chem. Soc.
111 (16) , 6452 (1989). 54. K. Cobos, R, Rodríguez, E. Quiñoá, R. Riguera, F. Freire, Angew. Chem. Int. Ed. , 23724 (2020). 55. Y. Nagata Y, T. Yamada, T. Adachi, Y. Akai, T. Yamamoto, M. Suginome, J. Am Chem Soc. , 10104 (2013). 56. K. Vanommeslaeghe, E. Hatcher, C. Acharya, S. Kundu, S. Zhong, J. Shim, E. Darian, O. Guvench, P. Lopes, I. Vorobyov, A. D. Jr.MacKerell, J. Comput. Chem. , 671 (2010). 57. M. J. Frisch, M. Head-Gordon, J. A. Pople, “Chem. Phys. Lett., , 135-49 (1992). 59. Schrödinger Release 2020-4: FEP+, Schrödinger, LLC, New York, NY, 2020. 60. A. Laio, M. Parrinello, Proc Natl Acad Sci U S A., , 12562 (2002). 61. A. Barducci, G. Bussi, M. Parrinello, Phys Rev Lett., , 020603 (2008). 62. G. M. Torrie, J. P. Valleau, J. Comput. Phys. , 187 (1977). 63. A. Pietropaolo, Y. Wang, T. Nakano, Angew. Chem. Int. Ed. , 2688 (2015). 64. G. Bussi, D. Donadio, M. Parrinello, J. Chem. Phys. , 014101 (2007). 65. B. Hess, J. Chem. Theory Comput. , 116 (2008). 66. U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem Phys. , 8577 (1995). 67. B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, J. Chem. Theory Comp. , 435-447 (2008). 68. The PLUMED consortium, Nat. Methods, , 670 (2019). 69. J. Kang, J. Jo, Y. Jo, S. Y. Lee, P. E. Keivanidis, G. Wegner, D. Y. Yoon, Polymer, , 5700 (2008). 70. G. Lieser, M. Oda, T. Miteva, A. Meisel, H.G. Nothofer, U. Scherf, D. Neher, Macromolecules,33