A Non-Orthogonal Block-Localized Effective Hamiltonian Approach for Chemical and Enzymatic Reactions

Abstract

The effective Hamiltonian-molecular orbital and valence bond (EH-MOVB) method based on non-orthogonal block-localized fragment orbitals has been implemented into the program CHARMM for molecular dynamics simulations of chemical and enzymatic reactions, making use of semiempirical quantum mechanical models. Building upon ab initio MOVB theory, we make use of two parameters in the EH-MOVB method to fit the barrier height and the relative energy between the reactant and product state for a given chemical reaction to be in agreement with experiment or high-level ab initio or density functional results. Consequently, the EH-MOVB method provides a highly accurate and computationally efficient QM/MM model for dynamics simulation of chemical reactions in solution. The EH-MOVB method is illustrated by examination of the potential energy surface of the hydride transfer reaction from trimethylamine to a flavin cofactor model in the gas phase. In the present study, we employed the semiempirical AM1 model, which yields a reaction barrier that is more than 5 kcal/mol too high. We use a parameter calibration procedure for the EH-MOVB method similar to that employed to adjust the results of semiempirical and empirical models. Thus, the relative energy of these two diabatic states can be shifted to reproduce the experimental energy of reaction, and the barrier height is optimized to reproduce the desired (accurate) value by adding a constant to the off-diagonal matrix element. The present EH-MOVB method offers a viable approach to characterizing solvent and protein-reorganization effects in the realm of combined QM/MM simulations.

Figure 1

Computed potential energy profile along the minimum energy path (Rc = R[C – H ] − R[H – N] ) for the hydride transfer reaction between trimethylamine and flavin model using EH-MOVB(AM1) in red, AM1 in light blue, B3LYP/6-31G(d) in navy blue, and M06-2X/6-31G(d) in green.

Figure 2

Optimized geometries for the reactant and product complexes and the transition state for the hydride transfer reaction depicted in Scheme 1.

Figure 3

Computed potential energy surfaces for the diabatic reactant state (blue), the diabatic product state (green), and the adiabatic ground state (red) along the minimum energy path.

Figure 4

Computed potential energy surfaces in Figure 3 for the diabatic reactant state (blue), the diabatic product state (green), and the adiabatic ground state (red) represented as a function of the energy difference between the reactant and product diabatic states (i.e., the energy-gap reaction coordinate).

Figure 5

Computed potential energy profiles for the reactant (blue) and product (green) diabatic states along with the adiabatic ground state (red) and the reference potential as a function of the coupling parameter linearly connecting the reactant and product potentials. This reaction path is called the reference minimum energy path, which has a different meaning from that of Figure 1.

Figure 6

The potential energy profile for the hydride transfer reaction between methylamine and the model flavin cofactor plotted against the geometrical reaction coordinate (eq 14) following the reference minimum energy path in Figure 5.

Figure 7

Computed potential energy profiles for the reactant (blue) and product (green) diabatic states and the adiabatic ground state (red) as a function of the energy gap-reaction coordinate for structures obtained along the reference minimum energy path in Figure 5.

Scheme 1

Schematic representation of the block-localization of molecular orbitals within individual molecular fragments for the reactant diabatic state (left) and the product diabatic state (right) for the hydride transfer reaction between trimethylamine (TMA-H) and a model for the flavin cofactor (Nf⁺). Atoms and charges in each rectangular specify the molecular block defined by the corresponding Lewis structure within which molecular molecular orbitals are localized. The antisymetric wave function constructed from the two blocks on the left-hand side of the arrow, TMA-H and Nf⁺, defines the reactant diabatic state, whereas that for the right-hand side blocks, TMA⁺ and H-Nf, define the product diabatic state.