Jingzhi Pu

CHARMM: The biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model‐building capabilities. The CHARMM program is applicable to problems involving a much broader class of many‐particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical‐molecular mechanical force fields, to all‐atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Tests of second-generation and third-generation density functionals for thermochemical kineticsElectronic supplementary information (ESI) available: Mean errors for pure and hybrid DFT methods. See http://www.rsc.org/suppdata/cp/b3/b316260e/

We report tests of second- and third-generation density functionals, for pure density functional theory (DFT) and hybrid DFT, against the BH6 representative barrier height database and the AE6 representative atomization energy database, with augmented, polarized double and triple zeta basis sets. The pure DFT methods tested are G96LYP, BB95, PBE, mPWPW91, VSXC, HCTH, OLYP, and OPW91 and the hybrid DFT methods tested are B1B95, PBE0, mPW1PW91, B97-1, B98, MPW1K, B97-2, and O3LYP. The performance of these methods is tested against each other as well as against first-generation methods (BP86, BLYP, PW91, B3PW91, and B3LYP). We conclude that the overall performance of the second-generation DFT methods is considerably better than the first-generation methods. The MPW1K method is very good for barrier height calculations, and none of the pure DFT methods outperforms any of the hybrid DFT methods for kinetics. The B1B95, VSXC, B98, OLYP and O3LYP methods perform best for atomization energies. Using a mean mean unsigned error criterion (MMUE) that involves two sizes of basis sets (both with polarization and diffuse functions) and averages mean unsigned errors in barrier heights and in atomization energy per bond, we find that VSXC has the best performance among pure functionals, and B97-2, MPW1K, and B1B95 have the best performance of all hybrid functionals tested.

How subunit coupling produces the γ-subunit rotary motion in F1-ATPase

FoF1-ATP synthase manufactures the energy “currency,” ATP, of living cells. The soluble F1 portion, called F1-ATPase, can act as a rotary motor, with ATP binding, hydrolysis, and product release, inducing a torque on the γ-subunit. A coarse-grained plastic network model is used to show at a residue level of detail how the conformational changes of the catalytic β-subunits act on the γ-subunit through repulsive van der Waals interactions to generate a torque that drives unidirectional rotation, as observed experimentally. The simulations suggest that the calculated 85° substep rotation is driven primarily by ATP binding and that the subsequent 35° substep rotation is produced by product release from one β-subunit and a concomitant binding pocket expansion of another β-subunit. The results of the simulation agree with single-molecule experiments [see, for example, Adachi K, et al. (2007) Cell 130:309–321] and support a tri-site rotary mechanism for F1-ATPase under physiological condition.

Parametrized direct dynamics study of rate constants of H with CH4 from 250 to 2400 K

Four implicit potential energy surfaces (PESs) with specific-reaction-parameters (SRP) are developed and tested for the reaction CH4+H→CH3+H2. The first is called MPW60 and is based on the modified Perdew–Wang hybrid density-functional method with the percentage of the Hartree–Fock exchange equal to 60%. The other three PESs are constructed with multi-coefficient correlation methods (MCCMs). The second is called MCOMP2-SRP, and the third is called MC-QCISD-SRP. Both of them are parametrized for this specific reaction by starting with their corresponding global parameters. The fourth is called MCG3-SRP and is based on the MCG3-CHO semiglobal parametrization with further refinement for this specific reaction. All four SRP surfaces have a classical forward barrier height of 14.8 kcal/mol, and all three MCCM SRP surfaces have a classical endoergicity of 3.3 kcal/mol. The stationary point geometries, vibrational frequencies, and zero-point-energies are reported for several standard single-level methods and MCC...

Validation of variational transition state theory with multidimensional tunneling contributions against accurate quantum mechanical dynamics for H+CH4→H2+CH3 in an extended temperature interval

Variational transition state theory with multidimensional tunneling contributions (VTST/MT) is tested against quantum mechanical rate constants for the reaction H+CH4→H2+CH3 at temperatures up to 1000 K. The VTST/MT method can be and has been applied to many reactions that cannot be treated by rigorous quantum dynamics methods. Studying the accuracy of VTST/MT by comparison with accurate quantal results that are becoming available for systems of increasing size is important for validating the theory. In the present study, covering a factor of five in temperature, the VTST/MT method is found to have a mean deviation from accurate quantal rate constants for a six-body reaction of only 13% and maximum deviation of only 23%.

Test of variational transition state theory with multidimensional tunneling contributions against an accurate full-dimensional rate constant calculation for a six-atom system

We present calculations of the H+CH4 reaction rate on the Jordan–Gilbert surface using canonical variational transition state theory with microcanonical optimized multidimensional tunneling contributions (CVT/μOMT). The purpose of the calculation is to compare the results to the recent accurate dynamical calculations of Bowman, Wang, Huang, Huarte-Larranaga, and Manthe for this potential energy surface. Over the full 200–500 K range for which accurate results are available we find a mean absolute deviation of only 17% and a maximum absolute deviation of 23%. This provides a rigorous validation of this popular method for a larger system than has previously been possible and indicates that previous validations for atom–diatom reactions were indeed indicative of the kind of accuracy one can obtain for larger systems.

Tests of potential energy surfaces for H+CH4↔CH3+H2: Deuterium and muonium kinetic isotope effects for the forward and reverse reaction

In previous work, three implicit potential energy surfaces with specific reaction parameters (SRP), namely MPW60, MC-QCISD-SRP, and MCG3-SRP, were developed for the reaction CH4+H→CH3+H2. Forward reaction rate constants obtained by variational transition state theory with multidimensional tunneling (VTST/MT) dynamics calculations on these surfaces give good agreement with recently re-analyzed experimental results. In the present work, again employing VTST/MT, kinetic isotope effects (KIEs) for isotopic variants of the title reaction in both the forward and reverse directions are examined on these SRP surfaces. Various primary and secondary deuterium (D) kinetic isotope reactions are studied; we also calculated the KIE for the reaction between methane and muonium (Mu), which is an ultralight isotope of protium with the Mu/H mass ratio being 0.113. The results are compared with several sets of experimental studies. With the VTST/MT dynamical method and harmonic vibrations, the proposed surfaces predict the ...

Isotropic Periodic Sum Treatment of Long-Range Electrostatic Interactions in Combined Quantum Mechanical and Molecular Mechanical Calculations

The isotropic periodic sum (IPS) method was extended to describe long-range electrostatic interactions in combined quantum mechanical and molecular mechanical (QM/MM) calculations. The resulting method, designated QM/MM-IPS, was tested for two ion association processes and a model SN2 reaction in aqueous solution. Potential of mean force (PMF) profiles and radial distribution functions computed from the QM/MM-IPS simulations were compared with those obtained by using the existing QM/MM-Ewald sum and cutoff (QM/MM-Cutoff) methods. In contrast to the QM/MM-Cutoff method, with which PMFs of ion separation tend to display a spurious linear drift, the QM/MM-IPS method successfully eliminates such artifacts, in excellent agreement with the QM/MM-Ewald results. The PMF obtained with the QM/MM-IPS method for the SN2 reaction that transfers an NH3 group between two chloride anions closely resembles that from the QM/MM-Ewald simulations. Compared with QM/MM-Ewald, the QM/MM-IPS method reduces the computational cost by 60-70% when a local region of 12 to 14 Å is used. These results suggest that the QM/MM-IPS method can be used as a reliable and efficient alternative to the QM/MM-Ewald method to incorporate long-range electrostatic effects in simulating solution-phase chemical reactions.

Lateral confinement of image electron wave function by an interfacial dipole lattice

Image-potential states on Cu(111) surfaces covered by thin films of C60 fullerene have been characterized by angle-resolved two-photon photoemission spectroscopy. Metal-to-molecule electron transfer within the first layer creates a 4×4 superlattice of surface dipoles. We show that such a surface dipole lattice provides lateral confinement of image-electron wave functions. Measurements of parallel dispersion indicate that the n=1 image state is localized in the presence of one monolayer of C60 but becomes delocalized by the addition of a second layer. Quantum mechanical calculations explain this in terms of the screening of the dipole potential, thus, restoring the free-electron behavior parallel to the surface. These results show that a surface dipole lattice can effectively control the interfacial electronic structure.

Searching for Saddle Points by Using the Nudged Elastic Band Method: An Implementation for Gas-Phase Systems.

A new implementation of the Nudged Elastic Band (NEB) optimization method is presented. This approach uses a global procedure that yields the whole reaction path, and thus it provides an alternative to the sequential optimization of the transition state and consequent calculation of the minimum energy path. Furthermore the algorithm is very useful when one is not sure if a saddle point exists, because it can be used to eliminate the possibility of a saddle point when one does not exist. Three different versions of the NEB algorithm have been implemented. The influences of various parameters and methodological choices on the performance of the method have been studied, and the quality of the results is assessed by comparison with the saddle point and minimum energy path calculations sequential method. Recommendations are made for algorithmic choices and default parameters.