Cavity size in reaction field theory

Supercomputer-Based Ensemble Docking Drug Discovery Pipeline with Application to Covid-19

We present a supercomputer-driven pipeline for in silico drug discovery using enhanced sampling molecular dynamics (MD) and ensemble docking. Ensemble docking makes use of MD results by docking compound databases into representative protein binding-site conformations, thus taking into account the dynamic properties of the binding sites. We also describe preliminary results obtained for 24 systems involving eight proteins of the proteome of SARS-CoV-2. The MD involves temperature replica exchange enhanced sampling, making use of massively parallel supercomputing to quickly sample the configurational space of protein drug targets. Using the Summit supercomputer at the Oak Ridge National Laboratory, more than 1 ms of enhanced sampling MD can be generated per day. We have ensemble docked repurposing databases to 10 configurations of each of the 24 SARS-CoV-2 systems using AutoDock Vina. Comparison to experiment demonstrates remarkably high hit rates for the top scoring tranches of compounds identified by our ensemble approach. We also demonstrate that, using Autodock-GPU on Summit, it is possible to perform exhaustive docking of one billion compounds in under 24 h. Finally, we discuss preliminary results and planned improvements to the pipeline, including the use of quantum mechanical (QM), machine learning, and artificial intelligence (AI) methods to cluster MD trajectories and rescore docking poses.

Variational Formulation of the Generalized Many-Body Expansion with Self-Consistent Charge Embedding: Simple and Correct Analytic Energy Gradient for Fragment-Based ab Initio Molecular Dynamics

The many-body expansion (MBE) and its extension to overlapping fragments, the generalized (G)MBE, constitute the theoretical basis for most fragment-based approaches for large-scale quantum chemistry. We reformulate the GMBE for use with embedding charges determined self-consistently from the fragment wave functions, in a manner that preserves the variational nature of the underlying self-consistent field method. As a result, the analytic gradient retains the simple "sum of fragment gradients" form that is often assumed in practice, sometimes incorrectly. This obviates (without approximation) the need to solve coupled-perturbed equations, and we demonstrate stable, fragment-based ab initio molecular dynamics simulations using this technique. Energy conservation fails when charge-response contributions to the Fock matrix are neglected, even while geometry optimizations and vibrational frequency calculations may yet be accurate. Stable simulations can be recovered by means of straightforward modifications introduced here, providing a general paradigm for fragment-based ab initio molecular dynamics.

Efficient Implicit Solvation Method for Full Potential DFT

With the advent of efficient electronic structure methods, effective continuum solvation methods have emerged as a way to, at least partially, include solvent effects into simulations without the need for expensive sampling over solvent degrees of freedom. The multipole moment expansion (MPE) model, while based on ideas initially put forward almost 100 years ago, has recently been updated for the needs of modern electronic structure calculations. Indeed, for an all-electron code relying on localized basis sets and-more importantly-a multipole moment expansion of the electrostatic potential, the MPE method presents a particularly cheap way of solving the macroscopic Poisson equation to determine the electrostatic response of a medium surrounding a solute. In addition to our implementation of the MPE model in the FHI-aims electronic structure theory code [ Blum , V. ; Comput. Phys. Commun. 2009 , 180 , 2175 - 2196 , DOI: 10.1016/j.cpc.2009.06.022 ], we describe novel algorithms for determining equidistributed points on the solvation cavity-defined as a charge density isosurface-and the determination of cavity surface and volume from just this collection of points and their local density gradients. We demonstrate the efficacy of our model on an analytically solvable test case, against high-accuracy finite-element calculations for a set of ≈140000 2D test cases, and finally against experimental solvation free energies of a number of neutral and singly charged molecular test sets [ Andreussi , O. ; J. Chem. Phys. 2012 , 136 , 064102 , DOI: 10.1063/1.3676407 ; Marenich , A. V. ; Minnesota Solvation Database , Version 2012; University of Minnesota : Minneapolis, MN, USA , 2012 . ]. In all test cases we find that our MPE approach compares very well with given references at computational overheads < 20% and sometimes much smaller compared to a plain self-consistency cycle.

A spherical cavity model for quadrupolar dielectrics

The dielectric properties of a fluid composed of molecules possessing both dipole and quadrupole moments are studied based on a model of the Onsager type (molecule in the centre of a spherical cavity). The dielectric permittivity ε and the macroscopic quadrupole polarizability αQ of the fluid are related to the basic molecular characteristics (molecular dipole, polarizability, quadrupole, quadrupolarizability). The effect of αQ is to increase the reaction field, to bring forth reaction field gradient, to decrease the cavity field, and to bring forth cavity field gradient. The effects from the quadrupole terms are significant in the case of small cavity size in a non-polar liquid. The quadrupoles in the medium are shown to have a small but measurable effect on the dielectric permittivity of several liquids (Ar, Kr, Xe, CH4, N2, CO2, CS2, C6H6, H2O, CH3OH). The theory is used to calculate the macroscopic quadrupolarizabilities of these fluids as functions of pressure and temperature. The cavity radii are also determined for these liquids, and it is shown that they are functions of density only. This extension of Onsager's theory will be important for non-polar solutions (fuel, crude oil, liquid CO2), especially at increased pressures.

Reaction Pathway and Free Energy Barrier for Urea Elimination in Aqueous Solution

To accurately predict the free energy barrier for urea elimination in aqueous solution, we examined the reaction coordinates for the direct and water-assisted elimination pathways, and evaluated the corresponding free energy barriers by using the surface and volume polarization for electrostatics (SVPE) model-based first-principles electronic-structure calculations. Based on the computational results, the water-assisted elimination pathway is dominant for urea elimination in aqueous solution, and the corresponding free energy barrier is 25.3 kcal/mol. The free energy barrier of 25.3 kcal/mol predicted for the dominant reaction pathway of urea elimination in aqueous solution is in good agreement with available experimental kinetic data.

Reaction pathways and free energy profiles for spontaneous hydrolysis of urea and tetramethylurea: unexpected substituent effects

It has been difficult to directly measure the spontaneous hydrolysis rate of urea and, thus, 1,1,3,3-tetramethylurea (Me4U) was used as a model to determine the "experimental" rate constant for urea hydrolysis. The use of Me4U was based on an assumption that the rate of urea hydrolysis should be 2.8 times that of Me4U hydrolysis because the rate of acetamide hydrolysis is 2.8 times that of N,N-dimethyl-acetamide hydrolysis. The present first-principles electronic-structure calculations on the competing non-enzymatic hydrolysis pathways have demonstrated that the dominant pathway is the neutral hydrolysis via the CN addition for both urea (when pH < ~11.6) and Me4U (regardless of pH), unlike the non-enzymatic hydrolysis of amides where alkaline hydrolysis is dominant. Based on the computational data, the substituent shift of the free energy barrier calculated for the neutral hydrolysis is remarkably different from that for the alkaline hydrolysis, and the rate constant for the urea hydrolysis should be ~1.3 × 10(9)-fold lower than that (4.2 × 10(-12) s(-1)) measured for the Me4U hydrolysis. As a result, the rate enhancement and catalytic proficiency of urease should be 1.2 × 10(25) and 3 × 10(27) M(-1), respectively, suggesting that urease surpasses proteases and all other enzymes in its power to enhance the rate of reaction. All of the computational results are consistent with available experimental data for Me4U, suggesting that the computational prediction for urea is reliable.

Modeling in vitro inhibition of butyrylcholinesterase using molecular docking, multi-linear regression and artificial neural network approaches

Butyrylcholinesterase (BChE) has been an important protein used for development of anti-cocaine medication. Through computational design, BChE mutants with ∼2000-fold improved catalytic efficiency against cocaine have been discovered in our lab. To study drug-enzyme interaction it is important to build mathematical model to predict molecular inhibitory activity against BChE. This report presents a neural network (NN) QSAR study, compared with multi-linear regression (MLR) and molecular docking, on a set of 93 small molecules that act as inhibitors of BChE by use of the inhibitory activities (pIC₅₀ values) of the molecules as target values. The statistical results for the linear model built from docking generated energy descriptors were: r(2)=0.67, rmsd=0.87, q(2)=0.65 and loormsd=0.90; the statistical results for the ligand-based MLR model were: r(2)=0.89, rmsd=0.51, q(2)=0.85 and loormsd=0.58; the statistical results for the ligand-based NN model were the best: r(2)=0.95, rmsd=0.33, q(2)=0.90 and loormsd=0.48, demonstrating that the NN is powerful in analysis of a set of complicated data. As BChE is also an established drug target to develop new treatment for Alzheimer's disease (AD). The developed QSAR models provide tools for rationalizing identification of potential BChE inhibitors or selection of compounds for synthesis in the discovery of novel effective inhibitors of BChE in the future.

A density functional and quantum Monte Carlo study of glutamic acid in vacuo and in a dielectric continuum medium

We present density functional theory (DFT) and quantum Monte Carlo (QMC) calculations of the glutamic acid and glutamate ion in vacuo and in various dielectric continuum media within the polarizable continuum model (PCM). In DFT, we employ the integral equation formalism variant of PCM while, in QMC, we use a PCM scheme we have developed to include both surface and volume polarization. We investigate the gas-phase protonation thermochemistry of the glutamic acid using a large set of structural conformations, and find that QMC is in excellent agreement with the best available theoretical and experimental results. For the solvated glutamic acid and glutamate ion, we perform DFT calculations for dielectric constants, ε, between 4 and 78. We find that the glutamate ion in the zwitterionic form is more stable than the non-zwitterionic form over the whole range of dielectric constants, while the glutamic acid is more stable in its non-zwitterionic form at ε = 4. The dielectric constant at which the two glutamic acid species have the same energy depends on the cavity size and lies between 5 and 12.5. We validate these results with QMC for the two limiting values of the dielectric constant, and find qualitative agreement with DFT even though the solvent polarization is less pronounced at the QMC level.

Thermodynamics of chemical reactions with COSMO-RS: the extreme case of charge separation or recombination

Many technically relevant chemical processes in the condensed phase involve as elementary reactive steps the formation of ions from neutral species or, as the opposite, recombination of ions. Such reactions that generate or annihilate charge defy the standard gas phase quantum chemical treatment, and also continuum solvation models are only partially able to account for the right amount of stabilization in solution. In this work, for such types of reaction, a solvation treatment involving the COSMO-RS method is assessed, which leads to improved results, i.e., errors of only around 10 kJ/mol for both protic and aprotic solvents. The examples discussed here comprise protolysis reactions and organo halide heterolysis, for both of which a comparison with reliable experimental data is possible. It is observed that for protolysis, the quality of results does not strongly depend on the quantum chemical method used for energy calculation. In contrast, in the case of heterolytic carbon-chlorine bond cleavage, clearly better results are obtained for higher correlated (coupled cluster) methods or the density functional M06-2X, which is well known for its accuracy if applied to organic chemistry. This hints at least that the right answer is obtained for the right reason and not due to a compensation of errors from gas phase thermodynamics with those from the solvation treatment. Problems encountered with certain critical solvents or upon decomposing Gibbs free energies into heats or entropies of reaction are found to relate mostly to the parameterization of the H-bonding term within COSMO-RS.

SOLVENT REORGANIZATION ENERGY WITH DIELECTRIC GREEN FUNCTIONAL AND ITS APPLICATION TO RETURN ELECTRON TRANSFER IN TETRACYANOETHYLENE-HEXAMETHYLBENZENE SYSTEM

Based on classical electrodynamics, the proper integral formula for general cases is adopted to estimate the electrostatic free energy difference between two given states, e.g. equilibrium and non-equilibrium states in the charging process of solute molecule in solution. With adequate consideration of the change of dielectric property in establishing the non-equilibrium state, the electrostatic free energy for this state can be obtained in a form differing from the traditional theories. The general formula of the solvent reorganization energy is derived according to the definition of the free energy difference between the non-equilibrium and equilibrium states due to the same electronic state of the solute molecule. Applying it to the two-sphere model and single-sphere model, we obtain the new expressions of the solvent reorganization energy, which are convenient to be used in practice. In the investigation of the return electron transfer of tetracyanoethylene-hexmethylbenzene system, good results can be obtained by using our expression when they are compared with those from the experiment.

Catalytic mechanism of cytochrome P450 for 5'-hydroxylation of nicotine: fundamental reaction pathways and stereoselectivity

A series of computational methods were used to study how cytochrome P450 2A6 (CYP2A6) interacts with (S)-(-)-nicotine, demonstrating that the dominant molecular species of (S)-(-)-nicotine in CYP2A6 active site exists in the free base state (with two conformations, SR(t) and SR(c)), despite the fact that the protonated state is dominant for the free ligand in solution. The computational results reveal that the dominant pathway of nicotine metabolism in CYP2A6 is through nicotine free base oxidation. Further, first-principles quantum mechanical/molecular mechanical free energy (QM/MM-FE) calculations were carried out to uncover the detailed reaction pathways for the CYP2A6-catalyzed nicotine 5'-hydroxylation reaction. In the determined CYP2A6-(S)-(-)-nicotine binding structures, the oxygen of Compound I (Cpd I) can abstract a hydrogen from either the trans-5'- or the cis-5'-position of (S)-(-)-nicotine. CYP2A6-catalyzed (S)-(-)-nicotine 5'-hydroxylation consists of two reaction steps, that is, the hydrogen transfer from the 5'-position of (S)-(-)-nicotine to the oxygen of Cpd I (the H-transfer step), followed by the recombination of the (S)-(-)-nicotine moiety with the iron-bound hydroxyl group to generate the 5'-hydroxynicotine product (the O-rebound step). The H-transfer step is rate-determining. The 5'-hydroxylation proceeds mainly with the stereoselective loss of the trans-5'-hydrogen, that is, the 5'-hydrogen trans to the pyridine ring. The calculated overall stereoselectivity of ∼97% favoring the trans-5'-hydroxylation is close to the observed stereoselectivity of 89-94%. This is the first time it has been demonstrated that a CYP substrate exists dominantly in one protonation state (cationic species) in solution, but uses its less-favorable protonation state (neutral free base) to perform the enzymatic reaction.

An implicit solvent model for SCC-DFTB with Charge-Dependent Radii

Motivated by the need of rapidly exploring the potential energy surface of chemical reactions that involve highly charged species, we have developed an implicit solvent model for the approximate density functional theory, SCC-DFTB. The solvation free energy is calculated using the popular model that employs Poisson-Boltzmann for electrostatics and a surface-area term for non-polar contributions. To balance the treatment of species with different charge distributions, we make the atomic radii that define the dielectric boundary and solute cavity depend on the solute charge distribution. Specifically, the atomic radii are assumed to be linearly dependent on the Mulliken charges and solved self-consistently together with the solute electronic structure. Benchmark calculations indicate that the model leads to solvation free energies of comparable accuracy to the SM6 model (especially for ions), which requires much more expensive DFT calculations. With analytical first derivatives and favorable computational speed, the SCC-DFTB based solvation model can be effectively used, in conjunction with high-level QM calculations, to explore the mechanism of solution reactions. This is illustrated with a brief analysis of the hydrolysis of mono-methyl mono-phosphate ester (MMP) and tri-methyl mono-phosphate ester (TMP). Possible future improvements are also briefly discussed.

Structural assignment of 6-oxy purine derivatives through computational modeling, synthesis, X-ray diffraction, and spectroscopic analysis

6-Oxy purine derivatives have been considered as potential therapeutic agents in various drug discovery efforts reported in the literature. However, the structural assignment of this important class of compounds has been controversial concerning the specific position of a hydrogen atom in the structure. To theoretically determine the most favorable type of tautomeric form of 6-oxy purine derivatives, we have carried out first-principles electronic structure calculations on the possible tautomeric forms (A, B, and C) and their relative stability of four representative 6-oxy purine derivatives (compounds 1-4). The computational results in both the gas phase and aqueous solution clearly reveal that the most favorable type of tautomeric form of these compounds is A, in which a hydrogen atom bonds with the N1 atom on the purine ring. To examine the computational results, one of the 6-oxy purine derivatives (i.e., compound 4) has been synthesized and its structure has been characterized by X-ray diffraction and spectroscopic analysis. All of the obtained computational and experimental data are consistent with the conclusion that the 6-oxy purine derivative exists in tautomer A. The conclusive structural assignment reported here is expected to be valuable for future computational studies on 6-oxy purine derivative binding with proteins and for computational drug design involving this type of compounds.

Free Energies of Solvation with Surface, Volume, and Local Electrostatic Effects and Atomic Surface Tensions to Represent the First Solvation Shell

Building on the SVPE (surface and volume polarization for electrostatics) model for electrostatic contributions to the free energy of solvation with explicit consideration of both surface and volume polarization effects, on the SMx approach to including first-solvation-shell contributions, and on the linear relationship between the electric field and short-range electrostatic contributions found by Chipman, we have developed a new method for computing absolute aqueous solvation free energies by combining the SVPE method with semiempirical terms that account for effects beyond bulk electrostatics. The new method is called SMVLE, and the elements it contains are denoted by SVPE-CDSL where SVPE denotes accounting for bulk electrostatic interactions between solute and solvent with both surface and volume contributions, CDS denotes the inclusion of solvent cavitation, changes in dispersion energy, and possible changes in local solvent structure by a semiempirical term utilizing geometry-dependent atomic surface tensions as implemented in SMx models, and L represents the local electrostatic effect derived from the outward-directed normal electric field on the cavity surface. The semiempirical CDS and L terms together represent the deviation of short-range contributions to the free energy of solvation from those accounted for by the SVPE term based on the bulk solvent dielectric constant. A solute training set containing a broad range of molecules used previously in the development of SM6 is used here for SMVLE model calibration. The aqueous solvation free energies predicted by the parameterized SMVLE model correlate exceedingly well with experimental values. The square of the correlation coefficient is 0.9949 and the slope is 1.0079. Comparison of the final SMVLE model against the earlier SMx solvation model shows that the parameterized SMVLE model not only yields good accuracy for neutrals but also significantly increases the accuracy for ions, making it the best implicit solvation model to date for aqueous solvation free energies of ions. The semiempirical terms associated with the outward-directed electric field account in a physical way for the improvement in the predictive accuracy for ions. The SMVLE method greatly decreases the need to include explicit water molecules for accurate modeling of solvation free energies of ions.

First-Principles Determination of Molecular Conformations of Indolizidine (-)-235B' in Solution

Indolizidine (-)-235B' is a particularly interesting natural product, as it is the currently known, most potent and subtype-selective open-channel blocker of the alpha4beta2 nicotinic acetylcholine receptor (nAChR). In the current study, extensive first-principles electronic structure calculations have been carried out in order to determine the stable molecular conformations and their relative free energies of the protonated and deprotonated states of (-)-235B' in the gas phase, in chloroform, and in aqueous solution. The (1)H and (13)C NMR chemical shifts calculated using the computationally determined dominant molecular conformation of the deprotonated state are all consistent with available experimental NMR spectra of (-)-235B' in chloroform, which suggests that the computationally determined molecular conformations are reasonable. Our computational results reveal for the first time that two geminal H atoms on carbon-3 (C3) of (-)-235B' have remarkably different chemical shifts (i.e. 3.24 and 2.03 ppm). The computational results help one to better understand and analyze the experimental (1)H NMR spectra of (-)-235B'. The finding of remarkably different chemical shifts of two C3 geminal H atoms in a certain molecular conformation of (-)-235B' may also be valuable in analysis of NMR spectra of other related ring-containing compounds. In addition, the pK(a) of (-)-235B' in aqueous solution is predicted to be ~9.7. All of the computational results provide a solid basis for future studies of the microscopic and phenomenological binding of various receptor proteins with the protonated and deprotonated structures of this unique open-channel blocker of alpha4beta2 nAChRs. This computational study also demonstrates how one can appropriately use computational modeling and spectroscopic analysis to address the structural and spectroscopic problems that cannot be addressed by experiments alone.

First-principles determination of molecular conformations of cyclic adenosine 3',5'-monophosphate in gas phase and aqueous solution

Extensive first-principles electronic structure calculations were performed in this study to explore the possible molecular structures and their concentration distribution of an intracellular second messenger, that is, cyclic adenosine 3',5'-monophosphate (cAMP), and its protonated form (cAMPH) in the gas phase and aqueous solution. The calculations resulted in prediction of four different stable conformers of cAMP and eight different stable conformers of cAMPH and their relative Gibbs free energies in the gas phase and aqueous solution. All of the computational results consistently demonstrate that the predominant conformers of cAMP and cAMPH are always the cAMP-chair-anti and cAMPH-chair2-syn conformers, respectively, in both the gas phase and aqueous solution. It has been demonstrated that the free energy barriers calculated for the intertransformation reactions between different conformers are very low (below approximately 6 kcal/mol) such that the intertransformation reactions between different conformers are very fast so that the concentration distribution of the system can quickly reach the thermodynamic equilibration during the process of binding with a protein. The calculated phenomenological pKa of 3.66 is in good agreement with the experimental pKa of 3.9 reported in literature, suggesting that the computational predictions resulted from this study are reasonable.

An efficient implementation for determining volume polarization in self-consistent reaction field theory

An efficient algorithm of the surface and volume polarization for electrostatics (SVPE) method in self-consistent reaction field (SCRF) theory, denoted by SV(1)PE, has been proposed to simulate direct volume polarization potential with a single layer of point charges outside the solute cavity while the indirect effects of volume polarization on surface polarization are still simulated with multiple layers of point charges. The free energies of solvation calculated using the SV(1)PE algorithm (implemented in GAUSSIAN03) reproduce the corresponding values calculated using the standard SVPE implementation within an error of only approximately 0.1% when the solute cavity is defined by the standard 0.001e/a(0) (3) solute charge isodensity contour. The SV(1)PE results are much less sensitive to the used cavity size in comparison with the well-established surface and simulated volume polarization for electrostatics [SS(V)PE] method which simulates volume polarization through an additional surface charge distribution on the cavity surface. The SCRF calculations using the SV(1)PE method are more efficient than those using the original SVPE method.

Variable atomic radii for continuum-solvent electrostatics calculation

We have developed a method to improve the description of solute cavity defined by the interlocking-sphere model for continuum-solvent electrostatics calculations. Many models choose atomic radii from a finite set of atom types or uses an even smaller set developed by Bondi [J. Phys. Chem. 68, 441 (1964)]. The new model presented here allowed each atom to adapt its radius according to its chemical environment. This was achieved by first approximating the electron density of a molecule by a superposition of atom-centered spherical Gaussian functions. The parameters of the Gaussian functions were then determined by optimizing a function that minimized the difference between the properties from the model and those from ab initio quantum calculations. These properties included the electrostatics potential on molecular surface and the electron density within the core of each atom. The size of each atom was then determined by finding the radius at which the electron density associated with the atom fell to a prechosen value. This value was different for different chemical elements and was chosen such that the averaged radius for each chemical element in a training set of molecules matched its Bondi radius. Thus, our model utilized only a few adjustable parameters-the above density cutoff values for different chemical elements-but had the flexibility of allowing every atom to adapt its radius according to its chemical environment. This variable-radii model gave better solvation energy for 31 small neutral molecules than the Bondi radii did, especially for a quantum mechanics/Poisson-Boltzmann approach we developed earlier. The improvement was most significant for molecules with large dipole moment. Future directions for further improvement are also discussed.

Modeling the catalysis of anti-cocaine catalytic antibody: competing reaction pathways and free energy barriers

The competing reaction pathways and the corresponding free energy barriers for cocaine hydrolysis catalyzed by an anti-cocaine catalytic antibody, mAb15A10, were studied by using a novel computational strategy based on the binding free energy calculations on the antibody binding with cocaine and transition states. The calculated binding free energies were used to evaluate the free energy barrier shift from the cocaine hydrolysis in water to the antibody-catalyzed cocaine hydrolysis for each reaction pathway. The free energy barriers for the antibody-catalyzed cocaine hydrolysis were predicted to be the corresponding free energy barriers for the cocaine hydrolysis in water plus the calculated free energy barrier shifts. The calculated free energy barrier shift of -6.87 kcal/mol from the dominant reaction pathway of the cocaine benzoyl ester hydrolysis in water to the dominant reaction pathway of the antibody-catalyzed cocaine hydrolysis is in good agreement with the experimentally derived free energy barrier shift of -5.93 kcal/mol. The calculated mutation-caused shifts of the free energy barrier are also reasonably close to the available experimental activity data. The good agreement suggests that the protocol for calculating the free energy barrier shift from the cocaine hydrolysis in water to the antibody-catalyzed cocaine hydrolysis may be used in future rational design of possible high-activity mutants of the antibody as anti-cocaine therapeutics. The general strategy of the free energy barrier shift calculation may also be valuable in studying a variety of chemical reactions catalyzed by other antibodies or proteins through noncovalent bonding interactions with the substrates.

First-principles calculation of pKa for cocaine, nicotine, neurotransmitters, and anilines in aqueous solution

The absolute pKa values of 24 representative amine compounds, including cocaine, nicotine, 10 neurotransmitters, and 12 anilines, in aqueous solution were calculated by performing first-principles electronic structure calculations that account for the solvent effects using four different solvation models, i.e., the surface and volume polarization for electrostatic interaction (SVPE) model, the standard polarizable continuum model (PCM), the integral equation formalism for the polarizable continuum model (IEFPCM), and the conductor-like screening solvation model (COSMO). Within the examined computational methods, the calculations using the SVPE model lead to the absolute pKa values with the smallest root-mean-square-deviation (rmsd) value (1.18). When the SVPE model was replaced by the PCM, IEFPCM, and COSMO, the rmsd value of the calculated absolute pKa values became 3.21, 2.72, and 3.08, respectively. All types of calculated pKa values linearly correlate with the experimental pKa values very well. With the empirical corrections using the linear correlation relationships, the theoretical pKa values are much closer to the corresponding experimental data and the rmsd values become 0.51-0.83. The smallest rmsd value (0.51) is also associated with the SVPE model. All of the results suggest that the first-principles electronic structure calculations using the SVPE model are a reliable approach to the pKa prediction for the amine compounds.

Computational study of the release of H2 from ammonia borane dimer (BH3NH3)2 and its ion pair isomers

High-level electronic structure calculations have been used to map out the relevant portions of the potential energy surfaces for the release of H2 from dimers of ammonia borane, BH3NH3 (AB). Using the correlation-consistent aug-cc-pVTZ basis set at the second-order perturbation MP2 level, geometries of stationary points were optimized. Relative energies were computed at these points using coupled-cluster CCSD(T) theory with the correlation-consistent basis sets at least up to the aug-cc-pVTZ level and in some cases extrapolated to the complete basis set limit. The results show that there are a number of possible dimers involving different types of hydrogen-bonded interactions. The most stable gaseous phase (AB)2 dimer results from a head-to-tail cyclic conformation and is stabilized by 14.0 kcal/mol with respect to two AB monomers. (AB)2 can generate one or two H2 molecules via several direct pathways with energy barriers ranging from 44 to 50 kcal/mol. The diammoniate of diborane ion pair isomer, [BH4-][NH3BH2NH3+] (DADB), is 10.6 kcal/mol less stable than (AB)2 and can be formed from two AB monomers by overcoming an energy barrier of approximately 26 kcal/mol. DADB can also be generated from successive additions of two NH3 molecules to B2H6 and from condensation of AB with separated BH3 and NH3 molecules. The pathway for H2 elimination from DADB is characterized by a smaller energy barrier of 20.1 kcal/mol. The alternative ion pair [NH4+][BH3NH2BH3-] is calculated to be 16.4 kcal/mol above (AB)2 and undergoes H2 release with an energy barrier of 17.7 kcal/mol. H2 elimination from both ion pair isomers yields the chain BH3NH2BH2NH3 as product. Our results suggest that the neutral dimer will play a minor role in the release of H2 from ammonia borane, with a dominant role from the ion pairs as observed experimentally in ionic liquids and the solid state.

Heats of Formation of Boron Hydride Anions and Dianions and Their Ammonium Salts [BnHmy-][NH4+]y with y = 1−2

Thermochemical parameters of the closo boron hydride BnHn2- dianions, with n=5-12, the B3H8- and B11H14- anions, and the B5H9 and B10H14 neutral species were predicted by high-level ab initio electronic structure calculations. Total atomization energies obtained from coupled-cluster CCSD(T)/complete basis set (CBS) extrapolated energies, plus additional corrections were used to predict the heats of formation of the simplest BnHmy- species in the gas phase in kcal/mol at 298 K: DeltaHf(B3H8-)=-23.1+/-1.0; DeltaHf(B5H52-)=119.4+/-1.5; DeltaHf(B6H62-)=64.1+/-1.5; and DeltaHf(B5H9)=24.1+/-1.5. The heats of formation of the larger species were evaluated by the G3 method from hydrogenation reactions (values at 298 K, in kcal/mol with estimated error bars of+/-3 kcal/mol): DeltaHf(B7H72-)=51.8; DeltaHf(B8H82-)=46.1; DeltaHf(B9H92-)=24.4; DeltaHf(B10H102-)=-12.5; DeltaHf(B11H112-)=-11.8; DeltaHf(B12H122-)=-86.3; DeltaHf(B11H14-)=-57.3; and DeltaHf(B10H14)=18.7. A linear correlation between atomization energies of the dianions and energies of the BH units was found. The heats of formation of the ammonium salts of the anions and dianions were predicted using lattice energies (UL) calculated from an empirical expression based on ionic volumes. The UL values (0 K) of the BnHn2- dianions range from 319 to 372 kcal/mol. The values of UL for the B3H8- and B11H14- anions are 113 and 135 kcal/mol, respectively. The calculated lattice energies and gas-phase heats of formation of the constituent ions were used to predict the heats of formation of the ammonium crystal salts [BnHmy-][NH4+]y. These results were used to evaluate the thermodynamics of the H2 release reactions from the ammonium hydro-borate salts.

On the Temperature Dependence of Amide I Frequencies of Peptides in Solution

The temperature dependence of the amide I vibrational frequencies of peptides in solution was investigated. In D2O, the amide I' bands of both an alpha-helical oligopeptide, the random-coil poly(L-lysine), and the simplest amide, N-methyl acetamide (NMA), exhibit linear frequency shifts of approximately 0.07 cm(-1)/degrees C with increasing temperature. Similar amide I frequency shifts are also observed for NMA in both polar (acetonitrile and DMSO) and nonpolar (1,4-dioxane) organic solvents, thus ruling out hydrogen-bonding strength as the cause of these effects. The experimental NMA amide I frequencies in the organic solvents can be accurately described by a simple theory based on the Onsager reaction field with temperature-dependent solvent dielectric properties and a solute molecular cavity. DFT-level calculations (BPW91/cc-pVDZ) for NMA with an Onsager reaction field confirm the significant contribution of the molecular cavity to the predicted amide I frequencies. Comparison of the computations to experimental data shows that the frequency-dependent response of the reaction field, taken into account by the index of refraction, is crucial for describing the amide I frequencies in polar solvents. The poor predictions of the model for the NMA amide I band in D2O might be due, in part, to the unknown temperature dependence of the refractive index of D2O in the mid-IR range, which was approximated by the available values in the visible region.

Theoretical studies of the transition-state structures and free energy barriers for base-catalyzed hydrolysis of amides

The transition-state structures and free energy barriers for the rate-determining step (i.e. the formation of a tetrahedral intermediate) of base-catalyzed hydrolysis of a series of amides in aqueous solution have been studied by performing first-principle electronic structure calculations using a hybrid supermolecule-polarizable continuum approach. The calculated results and a revisit of recently reported experimental proton inventory data reveal that the favorable transition-state structure optimized for the tetrahedral intermediate formation of hydroxide ion-catalyzed hydrolysis of formamide may have three solvating water molecules remaining on the attacking hydroxide oxygen and two additional water molecules attached to the carbonyl oxygen of formamide. The calculated results have also demonstrated interesting substituent effects on the optimized transition-state geometries, on the transition-state stabilization, and on the calculated free energy barriers for the base-catalyzed hydrolysis of amides. When some or all of the hydrogen atoms of formamide are replaced by methyl groups, the total number of water molecules hydrogen-bonding with the attacking hydroxide in the transition state decreases from three for formamide to two for N-methylacetamide, N,N-dimethylformamide (DMF), and N,N-dimethylacetamide (DMA). The larger substituents of the amide hinder the solvent water molecules approaching the attacking hydroxide oxygen in the transition state and, therefore, destabilize the transition-state structure and increase the free energy barrier. By using the optimized most favorable transition-state structures, the calculated free energy barriers, i.e., 21.6 (or 21.7), 22.7, 23.1, and 26.0 kcal/mol for formamide, N-methylacetamide, DMF, and DMA, respectively, are in good agreement with the available experimental free energy barriers, i.e., 21.2, 21.5, 22.6, and 24.1 kcal/mol for formamide, N-methylacetamide, DMF, and DMA, respectively.

New formulation and implementation for volume polarization in dielectric continuum theory

In the use of dielectric continuum theory to model bulk solvation effects on the electronic structure and properties of a solute, volume polarization contributions due to quantum mechanical penetration of the solute charge density outside the cavity nominally enclosing it are known to be significant. This work provides a new formulation and implementation of methods for solution of the requisite Poisson equation. In previous formulations the determination of the surface polarization contribution required evaluation of the difficult to calculate electric field generated by the volume polarization. It is shown that this problematic quantity can be eliminated in favor of other more easily evaluated quantities. That formal advance also opens the way for a more efficient apparatus to be implemented for calculation of the direct contribution of volume polarization to the solvation energy. The new formulation and its practical implementation are described, and illustrative numerical results are given for several neutral and ionic solutes to study the convergence and precision in practice.

Cation electric field is related to hydration energy

Electronic structure calculations on ions that use dielectric continuum theory to mimic solvent around the bare ionic solutes are often prone to make large errors in the hydration energies. It is found for cations that much of the error can be accounted for by a simple linear correlation with the maximum value found anywhere on the dielectric cavity surface of the solute potential or, even better, the outgoing normal electric field, thus mirroring analogous results previously obtained for anions. This correlation allows for significantly improved estimates of cation hydration energies while still retaining cavities of physically reasonable size in determination of the bulk dielectric contributions.

Modeling multiple species of nicotine and deschloroepibatidine interacting with alpha4beta2 nicotinic acetylcholine receptor: from microscopic binding to phenomenological binding affinity

A variety of molecular modeling, molecular docking, and first-principles electronic structure calculations were performed to study how the alpha4beta2 nicotinic acetylcholine receptor (nAChR) binds with different species of two typical agonists, (S)-(-)-nicotine and (R)-(-)-deschloroepibatidine, each of which is distinguished by different free bases and protonation states. On the basis of these results, predictions were made regarding the corresponding microscopic binding free energies. Hydrogen-bonding and cation-pi interactions between the receptor and the respective ligands were found to be the dominant factors differentiating the binding strengths of different microscopic binding species. The calculated results and analyses demonstrate that, for each agonist, all the species are interchangeable and can quickly achieve a thermodynamic equilibrium in solution and at the nAChR binding site. This allows quantitation of the equilibrium concentration distributions of the free ligand species and the corresponding microscopic ligand-receptor binding species, their pH dependence, and their contributions to the phenomenological binding affinity. The predicted equilibrium concentration distributions, pK(a) values, absolute phenomenological binding affinities, and their pH dependence are all in good agreement with available experimental data, suggesting that the computational strategy from the microscopic binding species and affinities to the phenomenological binding affinity is reliable for studying alpha4beta2 nAChR-ligand binding. This should provide valuable information for future rational design of drugs targeting nAChRs. The general strategy of the "from-microscopic-to-phenomenological" approach for studying interactions of alpha4beta2 nAChRs with (S)-(-)-nicotine and (R)-(-)-deschloroepibatidine may also be useful in studying other types of ligand-protein interactions involving multiple molecular species of a ligand and in associated rational drug design.

First-principle studies of intermolecular and intramolecular catalysis of protonated cocaine

We have performed a series of first-principles electronic structure calculations to examine the reaction pathways and the corresponding free energy barriers for the ester hydrolysis of protonated cocaine in its chair and boat conformations. The calculated free energy barriers for the benzoyl ester hydrolysis of protonated chair cocaine are close to the corresponding barriers calculated for the benzoyl ester hydrolysis of neutral cocaine. However, the free energy barrier calculated for the methyl ester hydrolysis of protonated cocaine in its chair conformation is significantly lower than for the methyl ester hydrolysis of neutral cocaine and for the dominant pathway of the benzoyl ester hydrolysis of protonated cocaine. The significant decrease of the free energy barrier, approximately 4 kcal/mol, is attributed to the intramolecular acid catalysis of the methyl ester hydrolysis of protonated cocaine, because the transition state structure is stabilized by the strong hydrogen bond between the carbonyl oxygen of the methyl ester moiety and the protonated tropane N. The relative magnitudes of the free energy barriers calculated for different pathways of the ester hydrolysis of protonated chair cocaine are consistent with the experimental kinetic data for cocaine hydrolysis under physiologic conditions. Similar intramolecular acid catalysis also occurs for the benzoyl ester hydrolysis of (protonated) boat cocaine in the physiologic condition, although the contribution of the intramolecular hydrogen bonding to transition state stabilization is negligible. Nonetheless, the predictability of the intramolecular hydrogen bonding could be useful in generating antibody-based catalysts that recruit cocaine to the boat conformation and an analog that elicited antibodies to approximate the protonated tropane N and the benzoyl O more closely than the natural boat conformer might increase the contribution from hydrogen bonding. Such a stable analog of the transition state for intramolecular catalysis of cocaine benzoyl-ester hydrolysis was synthesized and used to successfully elicit a number of anticocaine catalytic antibodies.

Polarizable continuum model study of solvent effects on electronic circular dichroism parameters

We present an implementation of the polarizable continuum model for the calculation of solvent effects on electronic circular dichroism spectra. The computational model used is density functional theory in the length-gauge formulation, and gauge-origin independence is ensured through the use of London atomic orbitals. Results of calculations carried out for methyloxirane and bicyclic ketones, camphor, norcamphor, norbornenone, and fenchone are presented, and the theoretically obtained solvent effects are compared with experimental observations.

Continuous medium theory for nonequilibrium solvation: III. Solvation shift by monopole approximation and multipole expansion in spherical cavity

According to the classical electrodynamics, a new and reasonable method about electrostatic energy decomposition of the solute-solvent system has been proposed in this work by introducing the concept of spring energy. This decomposition in equilibrium solvation gives the clear comprehension for different parts of total electrostatic free energy. Logically extending this cognition to nonequilibrium leads to the new formula of electrostatic free energy of nonequilibrium state. Furthermore, the general solvation shift for light absorption/emission has been reformulated and applied to the ideal sphere case with the monopole approximation and multipole expansion. Solvation shifts in vertical ionizations of atomic ions of some series of main group elements have been investigated with monopole approximation, and the variation tendency of the solvation shift versus atomic number has been discussed. Moreover, the solvation shift in photoionization of nitrate anion in glycol has been investigated by the multipole expansion method.

Reaction Pathways and Free Energy Barriers for Alkaline Hydrolysis of Insecticide 2-Trimethylammonioethyl Methylphosphonofluoridate and Related Organophosphorus Compounds: Electrostatic and Steric Effects

Reaction pathways and free energy barriers for alkaline hydrolysis of the highly neurotoxic insecticide 2-trimethylammonioethyl methylphosphonofluoridate and related organophosphorus compounds were studied by performing first-principles electronic structure calculations on representative methylphosphonofluoridates, (RO)CH3P(O)F, in which R = CH2CH2N+(CH3)3, CH3, CH2CH2C(CH3)3, CH2CH2CH(CH3)2, CH(CH3)CH2N+(CH3)3, and CH(CH3)CH2N(CH3)2. The dominant reaction pathway was found to be associated with a transition state in which the attacking nucleophile OH- and the leaving group F- are positioned on opposite sides of the plane formed by the three remaining atoms attached to the phosphorus in order to minimize the electrostatic repulsion between these two groups. The free energy barriers calculated for the rate-determining step of the dominant pathway are 12.5 kcal/mol when R = CH2CH2N+(CH3)3, 15.5 kcal/mol when R = CH3, 17.9 kcal/mol when R = CH2CH2C(CH3)3, 16.5 kcal/mol when R = CH2CH2CH(CH3)2, 13.4 kcal/mol when R = CH(CH3)CH2N+(CH3)3, and 18.7 kcal/mol when R = CH(CH(3))CH(2)N(CH(3))(2). The calculated free energy barriers are in good agreement with available experimentally derived activation free energies, i.e. 14.7 kcal/mol when R = CH(3), 13.4 kcal/mol when R = CH2CH2N+(CH3)3, and 13.9 kcal/mol when R = CH(CH3)CH2N+(CH3)3. A detailed analysis of the calculated energetic results and available experimental data suggests that the net charge of the molecule (M) being hydrolyzed is a prominent factor affecting the free energy barrier (DeltaG) for the alkaline hydrolysis of phosphodiesters, phosphonofluoridates, and related organophosphorus compounds. The electrostatic interactions between the attacking nucleophile OH- and the molecule M being hydrolyzed favor such an order of the free energy barrier: DeltaG(M(+)+OH-) < DeltaG(M0+OH-) < DeltaG(M(-)+OH-), where M+, M0, and M- represent the cationic, neutral, and anionic molecules, respectively. The change of the substituent R in (RO)CH(3)P(O)F from CH3 to CH2CH2N+(CH3)3 is associated with both the electrostatic and steric effects on the free energy barrier, but the electrostatic effect dominates the substituent shift of the free energy barrier. This helps to better understand why the alkaline hydrolysis of (RO)CH3P(O)F with R = CH2CH2N+(CH3)3 and CH(CH3)CH2N+(CH3)3 is significantly faster than that with R = CH3. The effect of electrostatic interaction also helps to understand why the rate constants for the alkaline hydrolysis of phosphodiesters, such as intramolecular second messenger adenosine 3',5'-phosphate (cAMP), are generally smaller than those for the alkaline hydrolysis of the phosphonofluoridates and related phosphotriesters.