Effect of Quantum Delocalization on Temperature Dependent Double Proton Transfer in Molecular Crystals of Terephthalic Acid

Supplementary Information: Effect of Quantum Delocalization on Temperature Dependent Double Proton Transfer in Molecular Crystals of Terephthalic Acid Unmesh Mondal,†Ivan Girotto,‡Ali Hassanali,¶and Prasenjit Ghosh∗,§ †Department of Chemistry, Indian Institute of Science Education and Research, Pune 411008, India ‡Data Science and Scientific Computing, International Center for Theoretical Physics, 34151, Trieste, Italy ¶Condensed Matter and Statistical Physics, International Center for Theoretical Physics, 34151, Trieste, Italy §Department of Physics, Indian Institute of Science Education and Research, Pune 411008, India E-mail: p[email protected],[email protected] CONTENTS 1. 0K structure of the three forms of Terephthalic Acid. 2. Convergence of PIGLET simulation w.r.t the number of beads. 3. NQE and finite temperature effects on double proton transfers. S0 4. Weak improper C−H· · · Ohydrogen bonds and its connection to the asymmetry of the DPT. 1 0K structure of the three forms of Terephthalic Acid. As mentioned in the main paper, Terephthalic acid exists in three polymorphs, namely, form I, II and III. While the lowest energy one, i.e., form I has been described in the main paper, here we discuss the results of form II and III. Figure S2 and S3 show the structure of form II and III respectively. The main difference between form I and the other forms is the stacking of the two layers. While in form II, the -COOH groups of two molecules in a chain that are involved in H-bonding, are stacked above the benzene ring of TPA in the layer below it, in form I and III, the molecules are perfectly stacked on top of each other. Additionally, in form I and II there are improper weak interchain H-bonds. Compared to form I, in form III, two neighbouring chains are shifted with respect to one another thereby disrupting the interchain C-H· · · O H-bonds. This also results in doubling of the cell size in form III. Table S1 compares our computed values of lattice parameters with the experimental results. We found that they are in reasonably good agreement. Further, from the computation of the cohesive energy, we find that while they are similar in form I and II, the TPA molecules are weakly bound in form III. Table S1: Comparison of the lattice parameters and cohesive energies (∆Ecoh) of the three polymorphs of TPA. Form a b c α β γ Vol. ∆Ecoh (˚ A) (˚ A) (˚ A) (◦) (◦) (◦) (˚ A3) (eV/mol.) I 9.625 7.606 3.623 73.80 105.07 139.44 165.13 1.811 I-exp19.54 7.73 3.75 70.85 107.36 137.78 174.76 II 9.623 5.247 4.945 85.76 136.10 104.36 165.17 1.809 II-exp19.54 5.34 5.02 86.95 134.65 104.90 172.77 III 8.899 10.449 3.612 86.30 89.36 91.06 335.12 1.748 III-exp28.940 10.442 3.79 90.0 91.21 90.0 353.72 S1 Figure S1: Crystal structure representation of TPA (form I) using ball and stick representation. (a) Side view (along the chain axis) represents the stacking of molecular chains (b) Top view (along z axis) shows layer-layer overlap. Top layer atoms in solid balls and bottom layer atoms as points. Green, red and white balls represent carbon, oxygen and hydrogen respectively Figure S2: Crystal structure representation of TPA (form II) using ball and stick representation. (a) Side view (along the chain axis) represents the stacking of molecular chains (b) Top view (along z axis) shows layer-layer overlap. Top layer atoms in solid balls and bottom layer atoms as points. Green, red and white balls represent carbon, oxygen and hydrogen respectively. S2 Figure S3: Crystal structure representation of TPA (form III) using ball and stick representation. (a) Side view (along the chain axis) represents the stacking of two molecular chains (b) Top view (along z axis) shows layer-layer overlap. Top layer atoms in solid balls and bottom layer atoms as points. Green, red and white balls represent carbon, oxygen and hydrogen respectively. 2 Convergence of PIGLET simulation w.r.t the number of beads. PIGLET simulations uses path integral formulation of statistical mechanics to incorporate NQE. In this approach, the quantum nucleus is transformed into a classical ring polymer of “P” beads/imaginary time slices where each bead is connected to (two) neighbouring beads through harmonic potential. In the limit P→ ∞, the transformation is exact and gives accurate quantum statistics of the nuclear system. However, in reality only a finite number of beads are used. Hence, it is necessary to perform convergence tests of the relevant properties with respect to the number of beads. Since the number of beads is inversely proportional to the temperature, convergence tests need to be performed at the lowest simulation temperature. In our study 70 K is the lowest temperature. Therefore, we have performed convergence test of the following thermodynamic quantities with respect to the number of beads at 70 S3 K: (i) Potential Energy (Figure S4(a)) and (ii) quantum kinetic energy (Figure S4(b)). The convergence tests were conducted on 4, 12, 16, 24, 32 and 64 beads where 4 and 12 beads were sampled for for 5 ps and the others for 10 ps. The potential and quantum kinetic energy (per atom) converged at 16 beads. Figure S4: Potential energy (top panel) and quantum kinetic energy (bottom panel) as a function of the number of beads. The convergence obtained for the thermodynamic quantities does not guarantee the accurate description of system properties. Hence, additional convergence tests were undertaken for quantities related to the investigation of proton transfers along the hydrogen bond. These are (i) the heavy atom distance (dO−O, Figure S5(a)), (ii) the hydrogen bond angle (O-H-O, Figure S5(b)), (iii) the radius of gyration of the proton (rG, Figure S5(c)) and (iv) the the effective free energy as a function of proton transfer coordinates (δ1and δ2, Figure S6). The converged number of beads for dO−Oand O-H-O is at 12 whereas the rGneeds 16 S4 beads to achieve convergence. The effective free energy as a function of the proton transfer coordinates needs a higher value of 64 beads. The reason for the distribution to not converge has to do with the fact that the proportion of the transferring hydrogen atom decreases as the number of beads increases. This is observed because Pis directly proportional to the force constant of the inter bead harmonic potential. Consequently, an increase in the number of beads slows the sampling of the potential energy landscape and would require longer simulation time. Therefore, an overall comparison for all the beads would not be completely relevant to check the convergence as we focus mainly on the nature of the proton transfer. In order to compare the nature of a specific proton transfer, i.e., whether the proton hops or tunnels, the radius of gyration as a function of sampling time is compared with the time evolution of the centroid proton transfer coordinate. The graphs at the two values of “P” (16 and 64) are shown in Figure S7. It is observed that nature of transfer does not change when the number of beads is increased from 16 to 64. Hence, 16 beads were used for the PIGLET simulations at all temperatures. 3 NQE and finite temperature effects on double proton transfers. In the main paper we have discussed the results only for 70 and 300 K. Here we show how the different quantities associated with the DPT evolves with increase in temperature. Figure S8 shows the temperature dependent evolution of the free energy surface as a function of dO−O and the proton transfer coordinate δ1,2, as obtained from the BOMD and PIGLET simulations. δ1,2has been described in the main paper. While, in the BOMD simulations proton transfers were not obtained until the system attains room temperature, in the PIGLET simulations even at 70 K there are proton transfer events. We observe that as the temperature is increased the free energy difference between the two minima decreases. In order to investigate whether the double proton transfers are stepwise or concerted, we S5 Figure S5: Convergence with respect to the number of beads for (a) the heavy atom separation (dO−O), (b) the hydrogen bond angle,(O-H-O) and (c) the radius of gyration of the O-H hydrogen (rG). S6 Figure S6: The effective free energy as a function of the proton transfer coordinates δ1and δ2for (a-f) 4, 12, 16, 24, 32 and 64 beads. Figure S7: The time evolution of the centroid proton transfer coordinate (δc o) and the radius of gyration (rG) for (a) 16 and (b) 64 beads. S7 Figure S8: Temperature evolution of the effective free energy as a function of the proton transfer coordinate (δ1,2) and the heavy atom distance (dO−O) obtained from (a) BOMD and (b) PIGLET simulations. have plotted the free energy surface as a function of the proton transfer coordinates δ1and δ2. Figure S9 shows the results obtained from BOMD and PIGLET simulations. Figure S9: Temperature evolution of the effective free energy as a function of the proton transfer coordinates δ1and δ2obtained from (a) BOMD and (b) PIGLET simulations. Figure S10 displays the joint density plots of δc o(centroid proton transfer coordinate) and r∥ G(component of rGparallel to the O-H bond) of the transferring protons as a function of S8