Isotopic Effect in Double Proton Transfer Process of Porphycene Investigated by Enhanced QM/MM Method

The proton transfer process in porphycenes holds potential applications in developing molecular switches, transistors and information storage devices^1,2. In particular, tautomerization in porphycenes through double proton transfer process has attracted wide interest in the fields of spectroscopy and photophysics². The inner hydrogen atoms of porphycene can migrate from one trans isomer to the other equivalent trans isomer through double proton transfer process as shown in Figure 1. Two mechanisms have been proposed for the double proton transfer process: the concerted and the stepwise mechanism^3,4. In the concerted double proton transfer process, both proton atoms move to the transition state synchronously in a symmetric way, whereas one proton completes the transfer before the other proton in a stepwise process. Two hydrogen atoms can transfer simultaneously or stepwise depending on the correlation strength between two hydrogen atoms⁵.

Isotopic substitution has been used to detect the structural properties of molecules and rate constants of reaction kinetics⁶. Single deuterium substitution in the inner hydrogen of porphycene leads to an asymmetric shape of the molecule. The hydrogen bond may expand or contract because of mass difference between the hydrogen and deuterium atoms. The isotopic substitution introduces a perturbation in the scaffold of porphycene. The question arises that whether asymmetric structure would affect the proton transfer process. Limbach and coworkers reported that the replacement of hydrogen with deuterium will compress both hydrogen bonds, and the cooperative coupling of two hydrogen bonds in porphycene may favor concerted mechanism⁷, whereas Yoshikawa stated the deuteration would make the stepwise mechanism contribute more than the concerted mechanism⁸. Experimental techniques, such as force spectroscopy, have been developed to capture tautomerization details in a single porphycene⁹. However, it is still challenging to determine the atomic details of proton transfer experimentally because of its transient nature.

Theoretical calculations and simulations can act as complementary tools in elucidating the reaction mechanisms of proton transfer. Among different theoretical methods, molecular dynamics (MD) simulations can monitor dynamic motions of each atom, and has been widely used to reveal complex mechanisms in chemical and enzymatic reactions. However, regular MD simulations tend to suffer from insufficient sampling issue, especially when high energy barrier exists in the process of interest. Therefore, enhanced sampling methods have been developed, which include transition path sampling^10,11, umbrella sampling (US)^12,13, and integrated tempering sampling (ITS)^14,15. Combination of different enhanced sampling methods can further increase sampling efficiency^16,17,18. To harness the enhanced sampling algorithms in simulating chemical reactions, we have implemented the selective integrated tempering sampling (SITS) method with quantum mechanical and molecular mechanical (QM/MM) potentials recently¹⁹. The proposed SITS-QM/MM method combines the advantages from both methods: the SITS method accelerates the sampling and can explore all possible reaction channels without prior knowledge of the reaction mechanism, and QM/MM provides more accurate description of the bond forming and bond breaking process, which cannot be simulated by MM methods solely. The implemented SITS-QM/MM approach has successfully uncovered concerted double proton transfer, uncorrelated and correlated stepwise double proton transfer mechanism in different systems, without pre-defining reaction coordinates¹⁹. For porphycene, the stepwise but correlated proton transfer character has been reported¹⁹. The hybrid SITS-QM/MM method was used to investigate the isotopic effect in porphycene in our study, and below are the detailed descriptions of the algorithm and protocol of our method.

We have implemented SITS method with hybrid QM/MM potentials. The effective potential of SITS was defined to include the potential energy at different temperatures with the weighting factors n_k to cover wider temperature ranges,

where, N is the number of canonical terms, β_k is the inverse temperature, and n_k is the corresponding weighting factor for each canonical component. U_E(R) and U_N(R) represent the enhanced and non-enhanced terms in SITS and are defined as,

U_s, U_se and U_e are the potential energy of sub-system, the interaction between sub-system and the environment, and the potential energy of environment. QM/MM potential is expressed as a hybrid summation of three components,

where U_qm, U_qm/mm, and U_mm are the internal energy term of the QM subsystem, the interaction energy between the QM and MM regions, and the interaction energy within the MM subsystem, respectively. The U_qm/mm term can be further divided into three components, which include the electrostatic, van der Waals, and covalent interaction energy terms between the QM and MM atoms,

We assign , and into one U_s term in SITS,

The full potential of the system was then decomposed into the energy of subsystem U_s, the interaction energy between the subsystem and the environment U_se, and the energy of environment U_e. For instance, in the system of the present work, the subsystem is the porphycence, and the environment the water.

The PMF profile along a collective variable τ(R) is derived as,

The generally used reaction coordinates for each hydrogen transfer of N1−H1··· N2 are q₁ = (r₁−r₂)/2 and q₂ = r₁ + r₂, where r₁ is the distance of N1-H1 and r₂ is the distance of H1-N2.

The method has been implemented in the QM/MM MD simulation package QM⁴D²⁰. The complete source code and documentation can be found here: http://www.qm4d.info/.

Generally, the SITS-QM/MM MD simulations involve four steps: pre-equilibrium (pre-sits); optimization n_k (opt-sits); production simulation and data analysis.