Peter Pulay of the University of Arkansas and Michel Dupuis of Washington State University are supported by an award from the Theoretical and Computational Chemistry program to develop methods that extend the current limits of applicability of electronic structure methods in two currently difficult research areas. They are developing algorithms to reduce the computational cost of ab initio based QM/MM (quantum/molecular mechanics) Monte Carlo or molecular dynamics free energy calculations of molecules electronically embedded in a medium (solvent or protein matrix) by several orders of magnitude. Additionally, they are working on a method that replaces atom-centered basis functions for the diffuse part of the electron density by plane waves in accurate quantum calculations.
Software and algorithms originated by these scientists have become widely adopted by the community. Pulay and Dupuis are exploiting the educational possibilities of computer modeling through graphical computer simulations. These ideas are implemented in a three week course being presented to undergraduates at several colleges and universities including an Historically Black University.
Background and objectives. The number of ways atoms can combine to molecules and larger structures is practically infinite (just the number of possible small organic molecules is estimated to be 10200 !). It is not possible to explore all, or a significant part of this vast universe, experimentally or computationally. However, computer simulations can be potentially much faster than experiments. Chemical behavior is mainly determined by the motion of the electrons around the atomic nuclei. The Schrödinger equation of quantum mechanics describes his motion and gives in principle a complete description of all chemical systems. However, solving these equations numerically is a huge computational task. The main goal of theoretical chemistry is to develop methods which are both accurate and efficient enough to be can be applied to practical chemical and biochemical problems. Currently, we can calculate in vacuum (i.e., by themselves) small molecules (~20 atoms) very accurately, and larger systems (a few hundred atoms) at a fair accuracy. This is often sufficient. However, the vast majority of chemistry and biochemistry takes place in solution, most importantly in water. Water is highly polar, and its electric field can have a profound effect on dissolved molecules. Modeling solutions is much more difficult than modeling molecules because the solvent molecules are continuously moving and turning, and one has to calculate the average over a large number (millions) of solvent arrangements. This can be done by fast Molecular Mechanics (MM) methods. However, MM cannot model chemical reactions, or unusual systems. For these, quantum mechanical (QM) calculations are need for the solute while the solvent can be modeled by MM. The combined QM/MM procedure is generally employed for accurate atomistic modeling of solutions. We have worked on two different projects under this award. The main project was to speed up the QM/MM simulation of solutions without sacrificing accuracy. A second project aimed at improving the description of the intermolecular forces (forces between molecules). These are weaker than chemical bonds but they accumulate, and for large molecules and materials, for instance for biomolecules, they become very important. Paradoxically, they are harder to calculate than the stronger chemical bonds. Our proposed method is somewhat technical but it aimed at removing two serious problems in the calculation of intermolecular forces. Accomplishments. We have developed a computer program which is improves the speed of the most widely used method for accurate QM/MM computer simulation of molecules in solutions by about four orders of magnitude (10,000 times or more) by eliminating quantum calculations in the simulation. The main idea is to precalculate the response of the QM system (the solute) to the electric field generated by the solvent, instead of carrying out a quantum calculation (other interactions, like Pauli repulsion and dispersion attraction, are highly localized and can be included at low cost, in contrast with the long-range Coulomb force). The interaction between the solute and the solvent can be calculated by perturbation theory. Calculating the first order energy (interaction between fixed charge distributions) is widely used and fast because it, like our method, does not require quantum calculations. However, its accuracy is insufficient, as shown in Figs. 1 and 2 which show the distribution of errors relative to full QM/MM calculations in the first-order treatment (broad curve) and our method (narrow peaks). Without second-order (polarization) contributions, absolute energies have rather large errors (~5-6 kcal/mol) with a broad distribution. Our method has negligible errors. Achieving such efficiency and accuracy will undoubtedly have a significant effect on QM/MM simulations, particularly for enzyme reactions. We have tested our method first on a well-known reaction, the SN2 chloride exchange reaction in substituted benzyl chlorides. Polarizability (second order) contributions are significant: they lower the reaction barrier by about 5 kcal/mol (Fig. 3). We have created the first genuine mixed basis (Gaussians and plane waves) ab initio program. This is a fundamentally new technique for accurate wavefunction-based quantum chemistry. (Plane waves are much easier to use in Density Functional Theory, DFT. Our method is quite different and more difficult and ambitious.) Using plane waves instead of diffuse basis functions eliminates two significant technical problems with accurate large basis set calculations: basis set overcompleteness and basis set superposition error. It will be most important for weak intermolecular interactions which attracted a lot of interest lately because of their importance in biochemistry and supramolecular chemistry. We have carried out benchmark calculations for intermolecular interactions, using our older (2005) Coupled Cluster program which can perform calculations on larger systems than other programs. The new plane wave program is restricted to small systems as we work on technical improvements). Fig. 4 shows some supramolecular systems for which we have determined benchmark interaction energies, in collaboration with a leading group in this field (P. Hobza).
