H. Lee Woodcock

Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations†

A hybrid quantum mechanical/molecular mechanical (QM/MM) potential energy function with Hartree‐Fock, density functional theory (DFT), and post‐HF (RIMP2, MP2, CCSD) capability has been implemented in the CHARMM and Q‐Chem software packages. In addition, we have modified CHARMM and Q‐Chem to take advantage of the newly introduced replica path and the nudged elastic band methods, which are powerful techniques for studying reaction pathways in a highly parallel (i.e., parallel/parallel) fashion, with each pathway point being distributed to a different node of a large cluster. To test our implementation, a series of systems were studied and comparisons were made to both full QM calculations and previous QM/MM studies and experiments. For instance, the differences between HF, DFT, MP2, and CCSD QM/MM calculations of H2O···H2O, H2O···Na+, and H2O···Cl− complexes have been explored. Furthermore, the recently implemented polarizable Drude water model was used to make comparisons to the popular TIP3P and TIP4P water models for doing QM/MM calculations. We have also computed the energetic profile of the chorismate mutase catalyzed Claisen rearrangement at various QM/MM levels of theory and have compared the results with previous studies. Our best estimate for the activation energy is 8.20 kcal/mol and for the reaction energy is −23.1 kcal/mol, both calculated at the MP2/6‐31+G(d)//MP2/6‐31+G(d)/C22 level of theory.

CHARMMing: A new, flexible, web portal for CHARMM

A new web portal for the CHARMM macromolecular modeling package, CHARMMing (CHARMM interface and graphics, http://www.charmming.org), is presented. This tool provides a user-friendly interface for the preparation, submission, monitoring, and visualization of molecular simulations (i.e., energy minimization, solvation, and dynamics). The infrastructure used to implement the web application is described. Two additional programs have been developed and integrated with CHARMMing: GENRTF, which is employed to define structural features not supported by the standard CHARMM force field, and a job broker, which is used to provide a portable method for using grid and cluster computing with CHARMMing. The use of the program is described with three proteins: 1YJP , 1O1O , and 1UFY . Source code is provided allowing CHARMMing to be downloaded, installed, and used by supercomputing centers and research groups that have a CHARMM license. Although no software can replace a scientists own judgment and experience, CHARMMing eases the introduction of newcomers to the molecular modeling discipline by providing a graphical method for running simulations.

Artificial reaction coordinate "tunneling" in free-energy calculations: the catalytic reaction of RNase H.

We describe a method for the systematic improvement of reaction coordinates in quantum mechanical/molecular mechanical (QM/MM) calculations of reaction free‐energy profiles. In umbrella‐sampling free‐energy calculations, a biasing potential acting on a chosen reaction coordinate is used to sample the system in reactant, product, and transition states. Sharp, nearly discontinuous changes along the resulting reaction path are used to identify coordinates that are relevant for the reaction but not properly sampled. These degrees of freedom are then included in an extended reaction coordinate. The general formalism is illustrated for the catalytic cleavage of the RNA backbone of an RNA/DNA hybrid duplex by the RNase H enzyme of Bacillus halodurans. We find that in the initial attack of the phosphate diester by water, the oxygen‐phosphorus distances alone are not sufficient as reaction coordinates, resulting in substantial hysteresis in the proton degrees of freedom and a barrier that is too low (∼10 kcal/mol). If the proton degrees of freedom are included in an extended reaction coordinate, we obtain a barrier of 21.6 kcal/mol consistent with the experimental rates. As the barrier is approached, the attacking water molecule transfers one of its protons to the O1P oxygen of the phosphate group. At the barrier top, the resulting hydroxide ion forms a penta‐coordinated phosphate intermediate. The method used to identify important degrees of freedom, and the procedure to optimize the reaction coordinate are general and should be useful both in classical and in QM/MM free‐energy calculations.

Multiscale Free Energy Simulations: An Efficient Method for Connecting Classical MD Simulations to QM or QM/MM Free Energies Using Non-Boltzmann Bennett Reweighting Schemes.

The reliability of free energy simulations (FES) is limited by two factors: (a) the need for correct sampling and (b) the accuracy of the computational method employed. Classical methods (e.g., force fields) are typically used for FES and present a myriad of challenges, with parametrization being a principle one. On the other hand, parameter-free quantum mechanical (QM) methods tend to be too computationally expensive for adequate sampling. One widely used approach is a combination of methods, where the free energy difference between the two end states is computed by, e.g., molecular mechanics (MM), and the end states are corrected by more accurate methods, such as QM or hybrid QM/MM techniques. Here we report two new approaches that significantly improve the aforementioned scheme; with a focus on how to compute corrections between, e.g., the MM and the more accurate QM calculations. First, a molecular dynamics trajectory that properly samples relevant conformational degrees of freedom is generated. Next, potential energies of each trajectory frame are generated with a QM or QM/MM Hamiltonian. Free energy differences are then calculated based on the QM or QM/MM energies using either a non-Boltzmann Bennett approach (QM-NBB) or non-Boltzmann free energy perturbation (NB-FEP). Both approaches are applied to calculate relative and absolute solvation free energies in explicit and implicit solvent environments. Solvation free energy differences (relative and absolute) between ethane and methanol in explicit solvent are used as the initial test case for QM-NBB. Next, implicit solvent methods are employed in conjunction with both QM-NBB and NB-FEP to compute absolute solvation free energies for 21 compounds. These compounds range from small molecules such as ethane and methanol to fairly large, flexible solutes, such as triacetyl glycerol. Several technical aspects were investigated. Ultimately some best practices are suggested for improving methods that seek to connect MM to QM (or QM/MM) levels of theory in FES.

Characterizing the Mechanism of the Double Proton Transfer in the Formamide Dimer

The double proton transfer in the formamide dimer is characterized computationally by combining density functional theory and ab initio methods. The intrinsic reaction coordinate (IRC) is obtained at the B3LYP level of theory. Energies of several points along the IRC are treated by the more rigorous focal point method to test the validity of the B3LYP functional. The reaction mechanism is examined in terms of the energy profile, the reaction force, the chemical potential, and the reaction electronic flux. The energy profile for the activation process of the formamide dimer to the imino ether product obtained with the B3LYP functional is in agreement with the results of the focal point method. Together with the reaction force analysis and the reaction electronic flux a precise assignment of the structural and electronic contributions to the activation barrier becomes possible. The results show that the reaction starts with a structural rearrangement, where the two dimers approach each other, and is followed by electronic changes before the system reaches the transition state. This electronic contribution to the activation barrier steers the activation process. After the transition state is reached, deviations of the B3LYP functional from the more accurate focal point energies become apparent, where the errors may be rationalized in terms of the treatment of exchange. The inconsistency could be assigned to the incapacity of the functional to describe delocalization effects over the whole system.

Computation of Hydration Free Energies Using the Multiple Environment Single System Quantum Mechanical/Molecular Mechanical Method

A recently developed MESS-E-QM/MM method (multiple-environment single-system quantum mechanical molecular/mechanical calculations with a Roothaan-step extrapolation) is applied to the computation of hydration free energies for the blind SAMPL4 test set and for 12 small molecules. First, free energy simulations are performed with a classical molecular mechanics force field using fixed-geometry solute molecules and explicit TIP3P solvent, and then the non-Boltzmann-Bennett method is employed to compute the QM/MM correction (QM/MM-NBB) to the molecular mechanical hydration free energies. For the SAMPL4 set, MESS-E-QM/MM-NBB corrections to the hydration free energy can be obtained 2 or 3 orders of magnitude faster than fully converged QM/MM-NBB corrections, and, on average, the hydration free energies predicted with MESS-E-QM/MM-NBB fall within 0.10-0.20 kcal/mol of full-converged QM/MM-NBB results. Out of five density functionals (BLYP, B3LYP, PBE0, M06-2X, and ωB97X-D), the BLYP functional is found to be most compatible with the TIP3P solvent model and yields the most accurate hydration free energies against experimental values for solute molecules included in this study.

Efficient Calculation of QM/MM Frequencies with the Mobile Block Hessian

The calculation of the analytical second derivative matrix (Hessian) is the bottleneck for vibrational analysis in QM/MM systems when an electrostatic embedding scheme is employed. Even with a small number of QM atoms in the system, the presence of MM atoms increases the computational cost dramatically: the long-range Coulomb interactions require that additional coupled perturbed self-consistent field (CPSCF) equations need to be solved for each MM atom displacement. This paper presents an extension to the Mobile Block Hessian (MBH) formalism for QM/MM calculations with blocks in the MM region and its implementation in a parallel version of the Q-Chem/CHARMM interface. MBH reduces both the CPU time and the memory requirements compared to the standard full Hessian QM/MM analysis, without the need to use a cutoff distance for the electrostatic interactions. Special attention is given to the treatment of link atoms which are usually present when the QM/MM border cuts through a covalent bond. Computational efficiency improvements are highlighted using a reduced chorismate mutase enzyme system, consisting of 24 QM atoms and 306 MM atoms, as a test example. In addition, the drug bortezomib, used for cancer treatment of myeloma, has been studied as a test case with multiple MBH block choices and both a QM and QM/MM description. The accuracy of the calculated Hessians is quantified by imposing Eckart constraints, which allows for the assessment of numerical errors in second derivative procedures. The results show that MBH within the QM/MM description not only is a computationally attractive method but also produces accurate results.

The microwave and infrared spectroscopy of benzaldehyde: Conflict between theory and experimental deductions

Recently, it has been proposed that ab initio calculations cannot accurately treat molecules comprised of a benzene ring with a pi-conjugated substituent, for example, benzaldehyde. Theoretical predictions of the benzaldehyde barrier to internal rotation are typically a factor of 2 too high in comparison to the experimental values of 4.67 (infared) and 4.90 (microwave) kcal mol(-1). However, both experiments use Pitzers 1946 model to compute the reduced moment of inertia and employ the experimentally observed torsional frequency to deduce benzaldehydes rotational barrier. When Pitzers model is applied to a system with a nonconjugated functional group, such as phenol, the model and theoretical values are in close agreement. Therefore, we conclude the model may not account for conjugation between the substituent and the pi-system of benzene. The experimental values of the benzaldehyde rotational barrier are therefore misleading. The true rotational barrier lies closer to the theoretically extrapolated limit of 7.7 kcal mol(-1), based on coupled cluster theory.

ProBiS-CHARMMing: Web Interface for Prediction and Optimization of Ligands in Protein Binding Sites.

Proteins often exist only as apo structures (unligated) in the Protein Data Bank, with their corresponding holo structures (with ligands) unavailable. However, apoproteins may not represent the amino-acid residue arrangement upon ligand binding well, which is especially problematic for molecular docking. We developed the ProBiS-CHARMMing web interface by connecting the ProBiS ( http://probis.cmm.ki.si ) and CHARMMing ( http://www.charmming.org ) web servers into one functional unit that enables prediction of protein-ligand complexes and allows for their geometry optimization and interaction energy calculation. The ProBiS web server predicts ligands (small compounds, proteins, nucleic acids, and single-atom ligands) that may bind to a query protein. This is achieved by comparing its surface structure against a nonredundant database of protein structures and finding those that have binding sites similar to that of the query protein. Existing ligands found in the similar binding sites are then transposed to the query according to predictions from ProBiS. The CHARMMing web server enables, among other things, minimization and potential energy calculation for a wide variety of biomolecular systems, and it is used here to optimize the geometry of the predicted protein-ligand complex structures using the CHARMM force field and to calculate their interaction energies with the corresponding query proteins. We show how ProBiS-CHARMMing can be used to predict ligands and their poses for a particular binding site, and minimize the predicted protein-ligand complexes to obtain representations of holoproteins. The ProBiS-CHARMMing web interface is freely available for academic users at http://probis.nih.gov.

Unlocking the binding and reaction mechanism of hydroxyurea substrates as biological nitric oxide donors.

Hydroxyurea is the only FDA approved treatment of sickle cell disease. It is believed that the primary mechanism of action is associated with the pharmacological elevation of nitric oxide in the blood; however, the exact details of this are still unclear. In the current work, we investigate the atomic level details of this process using a combination of flexible-ligand/flexible-receptor virtual screening coupled with energetic analysis that decomposes interaction energies. Utilizing these methods, we were able to elucidate the previously unknown substrate binding modes of a series of hydroxyurea analogs to hemoglobin and the concomitant structural changes of the enzyme. We identify a backbone carbonyl that forms a hydrogen bond with bound substrates. Our results are consistent with kinetic and electron paramagnetic resonance (EPR) measurements of hydroxyurea-hemoglobin reactions, and a full mechanism is proposed that offers new insights into possibly improving substrate binding and/or reactivity.