Elucidating the Nuclear Quantum Dynamics of Intramolecular Double Hydrogen Transfer in Porphycene

We address the double hydrogen transfer (DHT) dynamics of the porphycene molecule, a complex paradigmatic system in which the making and breaking of H-bonds in a highly anharmonic potential energy surface require a quantum mechanical treatment not only of the electrons but also of the nuclei. We combine density functional theory calculations, employing hybrid functionals and van der Waals corrections, with recently proposed and optimized path-integral ring-polymer methods for the approximation of quantum vibrational spectra and reaction rates. Our full-dimensional ring-polymer instanton simulations show that below 100 K the concerted DHT tunneling pathway dominates but between 100 and 300 K there is a competition between concerted and stepwise pathways when nuclear quantum effects are included. We obtain ground-state reaction rates of 2.19 × 1011 s–1 at 150 K and 0.63 × 1011 s–1 at 100 K, in good agreement with experiment. We also reproduce the puzzling N–H stretching band of porphycene with very good accuracy from thermostated ring-polymer molecular dynamics simulations. The position and line shape of this peak, centered at around 2600 cm–1 and spanning 750 cm–1, stem from a combination of very strong H-bonds, the coupling to low-frequency modes, and the access to cis-like isomeric conformations, which cannot be appropriately captured with classical-nuclei dynamics. These results verify the appropriateness of our general theoretical approach and provide a framework for a deeper physical understanding of hydrogen transfer dynamics in complex systems.


point interpolation.
Regarding DFT calculations, three exchange correlation (xc) functionals, namely PBE, 7 PBE0 8 and B3LYP, 9,10 and two types of dispersion interactions corrections (Tkatchenko-Scheffler (TS) 11 and many body dispersion (MBD) 12 were tested. All the density functional electronic structure calculations were performed using the all-electron FHI-aims program 13 with both light and tight settings. These default settings include both the size of the basis sets as well as other numerical parameters and are readily available with the program package.
In Table S1 and Table S2 we show the energy and relevant geometric parameters for all the stationary points obtained with the different xc functionals and CCSD(T). We present here only the results obtained with light settings. Their difference with respect to tight settings was, in all the cases, below 5 meV for energies differences and 0.01 for geometrical parameters. We note that for the three xc functionals tested, dispersion interactions play a minor role in the relative energies, in most cases representing differences below 5 meV.
Zero point energy (ZPE) contributions were only computed for the functionals including the TS dispersion correction. All xc functionals predict a potential energy surface (PES) which is characterized by two pairs of degenerate local minima, 4 first order saddle point and 1 second order in agreement with the literature. 14,15 B3LYP+TS is the only functional where both barriers are still present after the ZPE corrections and, therefore, does not predict a 'trivial' delocalization of the hydrogen between the two trans states.  In Table S3 we show the average absolute value of the dipole projected along the Hbonded nitrogen atoms.

Additional Information on Free Energy Profiles
In Fig. S1 we show the effective free energy profile for the hydrogen transfer coordinate obtained from PIMD simulations at 290K using PBE+MBD as exchange correlation functional.
These are qualitatively different from the rest of the simulations where the B3LYP+TS functional was used. In this case, the combination of the small proton transfer barrier height and ZPE produce a total delocalization of the protons.
S3 Figure S1: Effective free energy profile for the hydrogen transfer coordinate obtained from PIMD simulations at 290K using PBE+MBD exchange correlation functional. The contour lines are separated by 1 k B T .
In Fig. S2 and S3 we show the same free energy profile at 100K with B3LYP+TS exchange correlation functional and in Fig. S4 we present the 1D free energy cuts along the trans-trans, cis-cis and trans-cis paths. Finally in Fig. S5 we show the regions considered to compute the trans/cis population ratio.
. Structures that fall in the gray area are considered transition-state-like structures and therefore discarded.

Additional Information on Instanton Calculations
The ring-polymer instanton calculations were performed using the standard procedure described elsewhere. 17,18 In this theory, the thermal rate is given by the expression 19 where Q reac and Q inst refers to the reactant and instanton partition functions respectively, S is the instanton action obtained as β P U P (q) where U P is the potential energy (including the spring term) of the extended ring-polymer Hamiltonian andq is the geometry of the instanton. B N is a normalization factor, β P = β/P with β = 1/k B T and P is the number of replicas of the system used. A more detailed definition of these factors can be found in Ref. 18 In table S4 we show the convergence of the rates at 100 and 150 K with respect to S6 the number of beads. The correction factor that was used to scale the B3LYP+TS barrier to the CCSD(T) ones can be read as following: whereq (k) i represents the position of the i-th degree of freedom at the imaginary time slice k evaluated at the instanton geometry, N is the number of atoms, which is 38 for the porphycene molecule, and f is the correction factor defined as the ratio of the energy barriers obtained with CCSD(T) and B3LYP+TS for the considered path. It is assumed that the zero of energy is set the reactant.
In Fig. S6 we show the potential energy along the cumulative mass-weighted path-length and in Fig. S7 we present the 2D projection of the instanton pathways in the δ 1 and δ 2 plane. Figure S6: Potential energy along the cumulative mass-weighted path-length 20 at 150 K (red) and 100K (blue) for concerted (circles) and stepwise (squares) mechanisms. The zero of energy is set at the trans minimum energy geometry. This is the actual barrier which the system must tunnel through and clearly shows that the concerted path-length is longer.  Figure S7: 2D projection of the instanton geometry in the δ 1 and δ 2 plane at 100K. Each circle and square represents one replica of the discretized concerted and stepwise tunnelling pathway respectively. White dots are placed at the corresponding trans and cis minimumenergy geometry for reference.

Derivation of the Ring-Polymer Expansion of the Hessian used in the Instanton Calculations
In this section, we demonstrate the the ring-polymer expansion of the Hessian. We make use of the transformation matrix T (P, P ) defined in the Ref., 21 where P and P are the old (smaller) and new (larger) number of beads, defining the contraction operation as where r The P × P expansion matrix T (P , P ) can be defined as It can be shown that this matrix has the following properties: j (T (P , P )) ij = 1 i (T (P , P )) ij = P P .
We start our derivation with the following approximation where V (k) represents the potential energy corresponding to the k-th new replica and V (j) represents the potential energy corresponding to the k-th old replica. Summing over k and S9 using the properties of the T matrix we get the well known expression Applying the chain rule , the forces can be approximated as where F Since the force comes from the derivative of Eq. 7, one needs to sum again over k and use the properties of the T matrix to get Applying the chain rule a second time to Eq. 10, one can obtain the expression for the S10 ring polymer expansion of the Hessian following similar steps.
where H

Additional Information on IR Spectra
In Fig. S8 we present the vibrational density of states (290 K) of a water molecule calculated with thermostatted ring polymer molecular dynamics (TRPMD) coupled to colored noise thermostats on the Partridge-Schwenke potential. 22 We show that convergence is reached at 16 beads in this case. We thus simulated the IR spectrum of porphycene at 290 K with 16 beads in this paper.  Table S1) and therefore, the hydrogen transfer barrier disappears.
As already mentioned, this fact gives rise to a total delocalization of the hydrogen atoms which can be observed not only in the free energy profile but also in the IR spectrum. In

Model of coupled quantum harmonic oscillators
We used the methodology developed by Henri-Rousseau and Blaise. 23 The full Hamiltonian in the absence of damping may be written as H fast and H slow are given by and where q, Q, p, P , m ,M , ω and Ω correspond to the position, conjugate momenta, mass and frequency of the fast and slow mode respectively. The angular frequency of the fast mode is expanded up to first order obtaining where ω • is the fast mode frequency in absence of coupling. One may define a convenient dimensionless coupling parameter α as Consequently, using the adiabatic approximation to decouple the fast and slow degrees of motion and with the introduction of direct (γ) and indirect (γ • ) damping, one arrives to the following expression of the spectral density (SD): where n = 1/(e Ω/k B T ), k B is the Boltzmann constant and T is the absolute temperature.
Finally the SD in the frequency space is obtained by the Fourier transform of Eq. 18.
In table S5 we show the range of parameters used to reproduce the simulated TRPMD spectrum of porphycene. In Fig. S10 two illustrative examples of the spectra obtained with those parameters are presented. S14 Table S5: Range of parameters that yield a good agreement for the TRPMD spectrum of porphycene. We could also compute b values by finite difference displacements of the harmonic normal modes of porphycene (Eq. 13) and convert them to the corresponding α using the Eq. 16.
In table we report the computed α coupling parameter of all modes with both symmetric and antymmetric N-H stretch. We used the following definition of M i (the generalized mass of the i-th mode) needed to convert b into α where m j is the mass of the j-th degree of freedom, and |a j i | j is the j-th element of the i-th   Table S6: Dimensionless coupling parameter between harmonic modes with symmetric (α 1 , ω = 2902cm −1 ) and antisymmetric (α 2 ,ω = 2905cm −1 ) N-H stretching motion.