Importance of the electronic kinetic energy in polarizable force fields

Traditionally, a polarizable force field (PFF) is a classical model for the local electronic linear response of -molecules and solids. The recent extension of popular biomolecular force fields (AMBER, CHARMM, OPLS) with polarization terms is driven by the large demand for more reliable nonbonding interactions in molecular mechanics simulations. In addition, in other fields, PFFs are used to improve the accuracy of noncovalent force fields, e.g., for the simulations of adsorption isotherms of small guest molecule in porous silica and metal-organic frameworks. Nearly all PFFs attach an inducible dipole to each atom, either in the form of a Drude oscillator or as an inducible point dipole. Some PFFs also include variable charges (e.g., CPE) or even use them exclusively (EEM, Qeq, CHEQ, FlucQ, or FQ). In any case, these PFFs always use classical electrostatic interactions between charges and/or dipoles at different sites. It is generally accepted that atomic partial charges are the consequence of differences in intrinsic chemical potential and an external field acting on the molecule. It should therefore be surprising that most PFFs only use atomic inducible dipoles and thus skip the first term in the atomic multipole expansion. However, when variable atomic charges are included in a PFF, the model obtains two major undesirable properties: (i) the dipole polarizability scales cubically with system size (unlike the linear scaling one expects for dielectric systems) and (ii) molecules dissociate into fragments with fractional partial charges. We recently proposed a new type of PFF, atom-condensed Kohn–Sham DFT approximation to second order (ACKS2), that overcomes both limitations. Compared to a traditional PFF, ACKS2 introduces a new energy term for the electronic kinetic energy. This extension implies an extra nonclassical interaction for charges and dipoles at different sites. Additional advantages of the ACKS2 model include the direct derivation (not calibration) of PFF parameters from a Kohn–Sham wavefunction and a modest increase in computational cost compared with a conventional PFF. In this study, we will briefly review the theory of the ACKS2 model and discuss the pragmatic aspects of obtaining transferable and accurate ACKS2 parameters for a set of organic molecules

Faster potential energy surfaces: the art of making force fields

Density Functional Theory (DFT) has reshaped the field of computational chemistry over the past decades, which is mainly the merit of the attractive trade-off between the efficiency of DFT computations and the relative accuracy of the DFT potential energy surface with modern functionals. Despite these benefits, DFT computations are not always fast enough, especially when facing sampling problems on systems with many degrees of freedom, e.g. docking of a ligand on a protein surface. It is not only the scaling of the computational cost with system size, but also the increase in complexity that requires more samples and hence more computing time. There are two major approaches to reduce the burden of such sampling problems: (i) smarter sampling algorithms that extract the same information from less samples and (ii) faster methods to compute the potential energy of a molecular system. The latter option is the topic of this paper. Force fields are the fastest models to evaluate the potential energy (and nuclear forces) for a given molecular geometry. They are unfortunately also known for their limited accuracy and the inability to make and break chemical bonds during a simulation. Another major roadblock is the determination of force-field parameters for new systems, due to the lack of systematic calibration procedures. In this paper we present our recent methodological advances to surmount the typical weaknesses of (polarizable) force fields

The Monomer Electron Density Force Field (MEDFF) : a physically inspired model for noncovalent interactions

We propose a methodology to derive pairwise-additive noncovalent force fields from monomer electron densities without any empirical input. Energy expressions are based on the symmetry-adapted perturbation theory (SAPT) decomposition of interaction energies. This ensures a physically motivated force field featuring an electrostatic, exchange repulsion, dispersion, and induction contribution, which contains two types of parameters. First, each contribution depends on several fixed atomic parameters, resulting from a partitioning of the monomer electron density. Second, each of the last three contributions (exchange-repulsion, dispersion, and induction) contains exactly one linear fitting parameter. These three so-called interaction parameters in the model are initially estimated separately using SAPT reference calculations for the S66x8 database of noncovalent dimers. In a second step, the three interaction parameters are further refined simultaneously to reproduce CCSD(T)/CBS interaction energies for the same database. The limited number of parameters that are fitted to dimer interaction energies (only three) avoids ill-conditioned fits that plague conventional parameter optimizations. For the exchange repulsion and dispersion component, good results are obtained for all dimers in the S66x8 database using one single value for the associated interaction parameters. The values of those parameters can be considered universal and can also be used for dimers not present in the original database used for fitting. For the induction component such an approach is only viable for the dispersion dominated dimers in the S66x8 database. For other dimers (such as hydrogen-bonded complexes), we show that our methodology remains applicable. However, the interaction parameter needs to be determined on a case-specific basis. As an external validation:, the force field predicts interaction energies in good agreement with CCSD(T)/CBS values for dispersion dominated dimers extracted from an HIV-II protease crystal structure with a bound ligand (indinavir). Furthermore, experimental second virial coefficients of small alkanes and alkenes are well reproduced

Molecular Modeling of the Synthesis of Zeolites and Related Nanoporous Materials

Zeolites are microporous inorganic crystals with a framework that mainly contains silicon, oxygen and aluminum. The chemical composition is comparable to quartz, but the crystal structure is fundamentally different. Zeolites contain regular cavities (channels and cages) with a maximum diameter of 1.2 nanometer, large enough to contain small organic molecules. The walls of the cavities represent a huge internal surface (of the order of 1000 m2/g). Due to the well-defined shape of these channels and cages, only the molecules with a compatible structure can enter the internal cavities of a zeolite. This principle is called shape selectivity and it enables the separation of chemically similar but structurally different molecules. A Si/Al substitution in the framework introduces a catalytically active site on the walls of the zeolite cavities. During the past 60 years, all these specific properties have driven the development of new zeolites and the implementation of zeolites in industrial applications. The world-wide production of zeolites amounts to 4.2 million tons per annum, with applications in diverse markets such as petrochemistry, agriculture, animal husbandry, pharmacy, and so on. The most important applications in the western world include catalytic cracking in the petrochemical industry, ion exchange (softening and purification of water), and the separation and extraction of gases and solvents. A recent development in this field is the search for micro- and mesoporous materials. The IUPAC defines a microchannel as a pore with a diameter ranging from 0.25 nm to 2 nm, while a mesopore has a diameter between 2 nm and 50 nm. The limited diffusion rate of guest molecules in the micropores of a conventional zeolite is a well-known obstacle for the transport of reagents and reaction products. Another issue is the formation of side-products that reduce the diffusion rates and block the access to catalytically active sites. The regeneration of a zeolite catalyst, in which side products are burned, is an expensive and impractical procedure. One can avoid or at least reduce these limitations by introducing mesopores in the zeolite that have a diameter of the order of 1 to 50 nanometer. These mesopores act as highways for the transport of guest molecules. The Center for Surface Chemistry and Catalysis (COK) of the KULeuven has developed a new route for the synthesis of such micro- and mesoporous zeolites with well-controlled sizes of both types of pores. This revolutionary concept was the onset for the strategic basic research (SBO) project on biporous materials (BIPOM's). This PhD is mainly carried out in the context of the SBO-BIPOM project. A detailed description of the molecular mechanisms that give rise to the nucleation and growth of zeolite crystals is not yet available. Such insights are however of great interest for the tailor-made synthesis of zeolite catalysts in industrial applications. One of the major questions is the exact role of organic template molecules. Templates direct the synthesis process in a selective way to the formation of zeolites whose channels and cavities are structurally complementary to the shape of the template molecule. Another controversy is the exact structure of the precursors that precede the formation of zeolite crystals. A few (incompatible) models are postulated in the literature. Although several fragments of the synthesis mechanism have been unraveled through an impressive history of experimental studies, the remaining unknowns are not easily disentangled. Molecular modeling is a research field that plays a complementary role in this context. With the aid of computer simulations, one gains insights into molecular interactions that are not easily accessible to the experiment due to the small scale. There are two large categories of models to describe molecular interactions on a theoretical basis. In first instance there are the quantum-mechanical methods that solve the electronic many-body problem (approximately) to obtain the potential energy of a molecular system. An important subcategory are the ab initio methods, which rely only on the elementary quantum-mechanical postulates and do not depend on any empirical input. The second large category consists of the molecular mechanics methods that approximate molecular interactions with mainly empirical models. These models are computationally very efficient, but the downside is that molecular mechanics methods have only a limited accuracy and can not (or at least not correctly) describe important chemical processes such as the formation and the breaking of chemical bonds. Methods from both categories are applied to obtain microscopic properties of molecular systems such as optimal geometries, vibrational modes, reaction mechanisms, and so on. Only the quantum mechanical methods are reliable for the description of chemical reactions. Statistical mechanics provides the theoretical foundations to translate these microscopic data into relevant macroscopic observables. Molecular dynamics is a very versatile technique to apply the laws of statistical physics. By integrating the equations of motion of the molecular system, one runs through all the relevant microscopic states at a given temperature and pressure. Such simulations take into account the complete molecular environment such as a solvent or a zeolite framework. By using the proper statistical techniques, one can derive macroscopic parameters as averages over the ensemble of microscopic states. This thesis covers all the aspects of molecular modeling that we have addressed to gain more insights into the synthesis of zeolites. The larger part of the research deliverables in this PhD were only possible through the development of new software. The most visible example is ZEOBUILDER, a computer program to construct atomic models of biporous zeolites. The geometric models made with ZEOBUILDER are an essential starting point for molecular simulations. Another important aspect of this work is the development of molecular mechanics models that allow efficient simulations on zeolite precursors. In particular, we proposed new methods to obtain reliable force-field parameters. These technological and methodological tools are applied in theoretical studies on distinct intermediate steps of the synthesis of zeolites. Infrared spectra of zeolite precursors and zeolite nanocrystals have been derived from molecular dynamics simulations. This investigation confirms the experimental observation that a shift of the so-called MFI-fingerprint in the infrared spectrum can be associated with the formation of nanoscopic zeolite crystals. The tetrapropylammonium (TPA) template is known for the synthesis of MFI-structured zeolites such as ZSM-5 and Silicalite-I. We have studied the interactions between TPA and zeolite precursors proposed by the COK, using a broad scala of modeling techniques. The principal conclusion is that the preferential position of TPA with respect to the precursors corresponds to the place where the crossing of two channels will form in a later phase of the synthesis. This is in agreement with NMR measurements on macroscopic MFI-structured zeolite crystals from which the TPA has not been removed yet. Without the presence of the TPA, the investigated precursors would collapse, which is a confirmation of the function of the template molecule: providing a structural support for the initial zeolitic species. The new software and theoretical models in this thesis are the cornerstone for the modeling of zeolite synthesis processes. The development of polarizable force fields must be continued and their implementation in molecular dynamics software is indispensable to gain additional insights in the molecular mechanisms that lead to the formation of zeolites

Extended Hirshfeld: atomic charges that combine accurate electrostatics with transferability

C

Robust methods for predicting the transition states of chemical reactions: new approaches that focus on key coordinates

A new method for optimizing transition state and minima structures using redundant internal coordinates is presented. The new method is innovative because it allows the user to select a few key reduced coordinates, whose Hessian components will be accurately computed by finite differencing; the remaining elements of the Hessian are approximated with a quasi-Newton method. Usually the reduced coordinates are the coordinates that are involved in bond breaking/forming. In order to develop this method, several other innovations were made, including ways to (a) select the key reduced coordinates automatically, (b) guess the transition state quickly and efficiently, (c) choose dihedrals so that the “linear angle problem” is avoided, (d) robustly convert redundant internal coordinates to Cartesian coordinates. These, and other technical developments (e.g., new quasi-Newton Hessians, new trust-radius updates), were validated using a database of 7000 initial transition-state guesses for a diverse set of 140 chemical reactions