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

Summary

A protocol that uses enhanced QM/MM method to investigate the isotopic effect on the double proton transfer process in porphycene is presented here.

Abstract

The single deuterium substitution in porphycene leads to an asymmetric molecular geometry, which may affect the double proton transfer process in the porphycene molecule. In this study, we applied an enhanced QM/MM method called SITS-QM/MM to investigate hydrogen/deuterium (H/D) isotope effects on the double proton transfer in porphycene. Distance changes in SITS-QM/MM molecular dynamics simulations suggested that the deuterium substituted porphycene adopted the stepwise double proton transfer mechanism. The structural analysis and the free energy shifts of double proton transfer process indicated that the asymmetric isotopic substitution subtly compressed the covalent hydrogen bonds and may alter the original transition state location.

Introduction

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.

Protocol

1. Building model

1. Build porphycene structure: Open GaussView software by double-click the mouse. Then click Element Fragment button in the menu of the GaussView to choose the needed elements. Construct porphycene. Then click File button to save as the pdb file.
2. Solvate the model: Solvate porphycene in a cubic TIP3P²¹ water box with an edge length of 38 Å by issuing the command in the linux operating system: genbox_d -cp prp-vac.pdb -cs spc216.gro -o solv.pdb -maxsol 1484 -box 3.8.
3. Build deuterium porphycene: Issue the following command to generate topology file: cns < ppi_solv.inp. Then open prp-wat.psf with vi command and change the mass of H1 from 1.00800 to 2.01600 to replace one intramolecular hydrogen atom of porphycene with deuterium to build single deuterium substituted porphycene.
4. Set normal MD simulation parameters: Input method scctb, integral 0.5 fs, and cutoff 12 in the MD input file by opening it with vi command.
   NOTE: Adopt a cutoff distance of 12 Å for calculating both vdW and electrostatic interactions. Simulate the porphycene molecule with DFTB/MIO method²². Set the integration time step as 0.5 fs for MD simulations. Maintain the temperature of the simulation system at 300 K with Langevin thermostat. Then perform the simulations with QM⁴D software following the steps below.

2. Pre-sits

1. Set up temperature parameters: Input templow 260, temphigh 1100 and ntemp 160 in the input file.
   NOTE: The temperature range from 260K to 1100K was spread out by QM4D software to 160 temperature points during MD simulations. The template input files are included in Supplementary Files.
2. Initiate the pre-sits: Set runtype 100 and step 120,000 in the input file. Then issue the following command: $PATH/qm4d $INPUTFILE > $OUTPUTFILE.
   NOTE: The total step is 120,000 but can be adjusted depending on the specific need. The suggested parameters of MD simulations are saved in the $INPUTFILE. The same command is also used in the following opt-sits and production simulation steps, with the input file modified accordingly.
3. Calculating the decomposed energies
   1. Extract energy changes: During the pre-sits stage, monitor the energy of each term to calculate the mean values, as shown in Figure 1. Use grep Linux command to extract the energy as follows:
      grep ‘SITS-ener0’ $INPUTFILE | awk ‘{a+=$3;b+=$4;c+=$5}END{print a/NR,b/NR,c/NR}.
   2. Modify average energies in the MD input file: Calculate the average energies based on the output of the command line above, and modify the line vshift0 -30801.95; vshift1 -26.88; vshift2 -13888.28 in the input file with the newly generated averages.
      NOTE: Numbers -30801.85, -26.88 and -13888.28 are the average energies in the current model system. Please modify the values based on the specific systems.

3. Opt-sits

1. Initiate opt-sits: Set runtype 0 in the input file. Then initiate the QM⁴D program by typing the command as shown in step 2.2 to start the optimization step.
2. Monitoring the energy changes and n_k values.
   1. Plot the energy propagation with “grace” program and make sure the energy fluctuation can cover the lowest and the highest ends of the temperature range.
   2. After optimization, save the final n_k values of the opt-sits step into a new file, which is named as nk.dat in this protocol.

4. Running production simulations

1. Prepare MD input file: Set runtype 1 in the new input file to start the production simulation step. Specify the file name of stored n_k file as nkfile nk.dat in the input file. The number of time steps was set as 6,400,000 in the present systems.
2. Initiate production MD simulation: Issue the following command to start MD simulations: $PATH/qm4d $INPUTFILE > $OUTPUTFILE.
   NOTE: Make sure the n_k values can be read in by the QM⁴D software. The simulation time is system dependent so change the simulation step based on specific demands. Select a proper number of time steps to ensure enough simulation time for your own system. This step is likely time consuming, so save restart files to avoid restarting the production from the very beginning once being disrupted.

5 Data analysis

1. Monitoring the distance changes
   1. Monitor the bond forming and breaking process during the production phase, use the grep command to check the distance changes of H1-N1 and H1-N2 along the simulation time. The same operation can be conducted for H2-N3 and H2-N4. Then plot the distance propagation using the accumulated distance value during the production simulations.
2. Extracting reaction coordinates
   1. Extract the reaction coordinates and the energy terms from the production output file generated by QM⁴D by grep command:
      grep ‘dist 1’ $OUTPUTFILE | awk ‘{print $5}’ > distance1;
      grep ‘ener0’ $OUTPUTFILE > ener0.
   2. Organize data in four columns: q₁, q₂, U₀ and U’ (U₀ and U’ are the output normal energy and weighted energy), and write them into the data file at each time frame.
3. Calculating Free energy
   1. Calculate the free energy by issuing the following command:
      sits-pmf 300 $INPUTFILE PMF2 [hist_minx hist_maxx num_binsx] [hist_miny hist_maxy num_binsy] > $OUTPUTFILE.
      NOTE: sits-pmf is the histogram-based analysis method. [hist_minx hist_maxx num_binsx] defines the range and number of bins for the first reaction coordinate. The second reaction coordinates can be set by [hist_miny hist_maxy num_binsy].
   2. To project the free energy on the two-dimensional landscape, type the following command:
      sits-pmf 300 h1-2d.dat PMF2 -0.6 0.6 24 2.45 4.25 36 > sits-pmf.out.
      NOTE: Use a total of 24 bins and 36 bins to cover the distance changes in two selected reaction coordinates, q₁ and q₂, respectively. Save the 2D PMF data for each hydrogen/deuterium to the sits-pmf.out file.

Representative Results

The single deuterium substitution effect on double proton transfer process in porphycene was examined in the current protocol (Figure 1). The potential energy of QM sub-system and the water during pre-equilibrium and optimization step were checked to make sure the energy has been broadened to a wider energy range (Figure 2). The representative distance and angle changes (Figure 3 and Figure 4), and the projected free energy changes (Figure 5) were used to characterize the deuterium substitution effect on the geometry and proton transfer process of porphycene.

Figure 1. The structures of investigated molecules.
The structures of porphycene (A) and deuterated porphycene(B). Please click here to view a larger version of this figure.

Figure 2. The potential energy changes during MD simulations.
The potential energy changes of QM region (A) and environment (B) in pre-sits and opt-sits steps. Please click here to view a larger version of this figure.

Figure 3. The characteristic distance changes.
(A) The distance changes of H1-N1 and H2-N3 for porphycene, and (B) the distance changes of D1-N1 and H2-N3 for deuterated porphycene during SIRTS-QM/MM simulations; (C) the distribution of distance changes for H1-N1 and H2-N3 for porphycene, and (D) D1-N1 and H2-N3 for deuterated porphycene. Please click here to view a larger version of this figure.

Figure 4. The hydrogen bond angles during the production MD simulations.
The hydrogen bond angels for (A) prophycene and (B) deuterated porphycene. Please click here to view a larger version of this figure.

Figure 5. The free energy landscape of each hydrogen transfer process that was projected on two reaction coordinates (q₁, q₂).
(A) and (B) are 2D free energy landscapes of H1 and H2 transfer in porphycene; (C) and (D) 2D free energy landscape of D1 and H2 transfer in deuterated porphycene. Please click here to view a larger version of this figure.

Supplementary files. Topology file, force field parameter file, coordinates file and input file. Please click here to download the file.

Supplementary movie 1. Porphycene. Please click here to download the video.

Discussion

The structure of porphycene was shown in Figure 1. The electrostatic embedding QM/MM hybrid potential with SITS method was used to describe the chemical reactions in water^23,24. The proton transfer occurs within porphycene³ and thus porphycene is set as QM region and the reminding water is set as MM region. Herein we adopted DFTB/MIO as our QM method to treat the porphycene by balancing the efficiency and accuracy^22,25. As a sampling enhancement technique, SITS simulation was shown to broaden the distribution of U_s to high-energy regions and meantime preserve sufficient sampling around the energy region at the temperature of interest. For current case, the energy of U_s in the “opt-sits” step was broadened to wider ranges encompassing the energy of the standard MD simulations in the “pre-sits” step as shown in Figure 2. Meantime, the smooth energy changes of U_e indicated that the higher temperature in the QM subsystem would not bring perturbation into the environment. SITS-QM/MM method realized enhanced sampling at the interested QM region without affecting the potential energy of water.

From distance changes in Figure 3, we noticed that H1 transferred from N1 to N2 to form a transit cis state, and then initiated a consecutively rapid H2 transfer to arrive the other trans state again; and vice versa. The dynamics proton transfer process is shown in Supplementary Movie 1. Deuterium D1 transfer between N1 and N2 in the single deuterated porphycene invoked the transfer of H2 between N3 and N4. Asynchronous distance changes indicated the stepwise double proton transfer process for both of porphycene and single deuterium substituted porphycene. The similar distance distributions of D1-N1 and H2-N3 suggested that the cooperative effect on two hydrogen bonds²⁶. Consistent with previously reported primary geometric isotope effect²⁶, the distance of D1-N1 is shorter than the distance of H1-N1 (1.048 Å vs. 1.051 Å). As shown in Figure 3, we observed around 135, and 65 times of transfer of H or D for porphycene and its isotopomers within 3.2 ns MD simulations, respectively. Deuteriation may exert less effect on the hydrogen bond angles as shown in Figure 4. The sufficient sampling on the two reaction channels enabled us to calculate the free energy changes of each proton transfer. The obvious isotopic effect was observed in the 2D free energy landscape. The transition state has been shifted from (0.01 Å, 2.52 Å) to (-0.01 Å, 2.76 Å) as revealed from reaction coordinates (q₁, q₂) (see Figure 5). The higher q₂ value means the nonbonded hydrogen bond were expanded. This may come from the asymmetric scaffold of the deuterated porphycene.

Both proton transfer processes in porphycene and deuterated porphycene can be captured by SITS-QM/MM MD simulations without pre-defining the reaction coordinates. Moreover, SITS-QM/MM MD simulations revealed the structural difference that was introduced by isotopic effect. The hydrogen bond D1-N1 was shortened in comparison with H1-N1. The transition state has been shifted toward higher q₂ value because of an asymmetric shape caused by deuteration. Though only subtle difference was detected in covalent hydrogen bond, the distance difference may invoke bigger energy difference around the equilibrium bond distance. We are planning to further validate this observation at higher level QM method in the future studies.

The feasibility of SITS-QM/MM has been well validated in the double channel of reactions without pre-defining reaction coordinates in the present study. This method holds potential of searching reaction products from known reactant states if no prior reaction mechanism is provided. We have adopted DFTB/MIO method in our current implementation of SITS-QM/MM approach, and gained a better understanding of the isotopic effect. It is worth noting that the implemented approach is able to capture the free energy changes, but may not capture dynamic properties without considering the quantum tunneling effect. Still, this protocol acts as a starting point to investigate the chemical reaction mechanisms in the condensed environment. We expect SITS-QM/MM method to be extended to higher level QM methods and thus can exploit more complex systems in the future.

Materials

Name	Company	Catalog Number	Comments
operating system	CentOS Linux release 6.0
QM4D software	http://www.qm4d.info/		in-house program
Computer desktop	HP