Lead optimization mapper: automating free energy calculations for lead optimization

Alchemical free energy calculations hold increasing promise as an aid to drug discovery efforts. However, applications of these techniques in discovery projects have been relatively few, partly because of the difficulty of planning and setting up calculations. Here, we introduce lead optimization mapper, LOMAP, an automated algorithm to plan efficient relative free energy calculations between potential ligands within a substantial library of perhaps hundreds of compounds. In this approach, ligands are first grouped by structural similarity primarily based on the size of a (loosely defined) maximal common substructure, and then calculations are planned within and between sets of structurally related compounds. An emphasis is placed on ensuring that relative free energies can be obtained between any pair of compounds without combining the results of too many different relative free energy calculations (to avoid accumulation of error) and by providing some redundancy to allow for the possibility of error and consistency checking and provide some insight into when results can be expected to be unreliable. The algorithm is discussed in detail and a Python implementation, based on both Schrödinger’s and OpenEye’s APIs, has been made available freely under the BSD license.     

Modeling reaction pathways of low energy particle deposition on polymer surfaces via first principle calculations

The chemical processes that lead to polystyrene surface modification via low energy deposition of C(2)H(+), C(2)F(+), CH(2), CH(2)(+), and H(+) radicals and ions are examined using first principles calculations. Specifically, the reaction mechanisms responsible for products identified in classical molecular dynamics with reactive empirical bond-order potentials are examined using density functional theory. In addition, these calculations consider how the presence of charges on the incident particles changes the result for the CH(2) system through the comparison of barriers, transition states, and final products for CH(2) and CH(2)(+). The structures of the reaction species and energy barriers are determined using the B3LYP hybrid functional. Finally, CCSD/6-31G(d,p) single point energy calculations are carried out to obtain optimized energy barriers. The results indicate that the large variety of reactions occurring on the polystyrene surface are a consequence of complex interactions between the substrate and the deposited particles, which can easily be identified and characterized using advanced computational methodologies, such as first principle calculations.     

Calculations of the optical spectra of hydrocarbon radical cations based on Koopmans' theorem

The first few bands in the optical spectra of radical cations can often be interpreted in terms of A-type transitions that involve electron promotions from doubly occupied to the singly occupied molecular orbital (SOMO) and/or B-type transition which involve electron promotion from the SOMO to virtual molecular orbitals. We had previously demonstrated that, by making use of Koopmans' theorem, the energies of A-type transitions can be related to orbital energy differences between lower occupied MOs and the highest occupied MO (HOMO) in the neutral molecule, calculated at the geometry of the radical cation. We now propose that the energies of B-type transitions can be related similarly to energy differences between the lowest unoccupied MO (LUMO) and higher virtual MOs in the dication, also calculated at the geometry of the radical cation, by way of an extension of Koopmans' theorem to virtual MOs similar to that used sometimes to model resonances in electron scattering experiments. The optical spectra of the radical cations of several polyenes and aromatic compounds, the matrix spectra of which are known (or presented here for the first time), and for which CASSCF/CASPT2 calculations are available, are discussed in terms of these Koopmans-based models. Then the spectra of five poly(bicycloalkyl)-protected systems and that of hexabenzocoronene, compounds not amenable to higher level calculations, are examined and it is found that the Koopmans-type calculations allow a satisfactory interpretation of most of the features in these spectra. These simple calculations therefore provide a computationally inexpensive yet effective way to assign optical transitions in radical ions. Limitations of the model are discussed.     

Modeling the adiabatic connection in H2

Full configuration interaction (FCI) data are used to quantify the accuracy of approximate adiabatic connection (AC) forms in describing the ground state potential energy curve of H2, within spin-restricted density functional theory (DFT). For each internuclear separation R, accurate properties of the AC are determined from large basis set FCI calculations. The parameters in the approximate AC form are then determined so as to reproduce these FCI values exactly, yielding an exchange-correlation energy expressed entirely in terms of FCI-derived quantities. This is combined with other FCI-derived energy components to give the total electronic energy; comparison with the FCI energy quantifies the accuracy of the AC form. Initial calculations focus on a [1/1]-Padé-based form. The potential energy curve determined using the procedure is a notable improvement over those from existing DFT functionals. The accuracy near equilibrium is quantified by calculating the bond length and vibrational wave numbers; errors in the latter are below 0.5%. The molecule dissociates correctly, which can be traced to the use of virtual orbital eigenvalues in the slope in the noninteracting limit, capturing static correlation. At intermediate R, the potential energy curve exhibits an unphysical barrier, similar to that noted previously using the random phase approximation. Alternative forms of the AC are also considered, paying attention to size extensivity and the behavior in the strong-interaction limit; none provide an accurate potential energy curve for all R, although good accuracy can be achieved near equilibrium. The study demonstrates how data from correlated ab initio calculations can provide valuable information about AC forms and highlight areas where further theoretical progress is required.     

Computation of binding free energy with molecular dynamics and grand canonical Monte Carlo simulations

The binding of a ligand to a receptor is often associated with the displacement of a number of bound water molecules. When the binding site is exposed to the bulk region, this process may be sampled adequately by standard unbiased molecular dynamics trajectories. However, when the binding site is deeply buried and the exchange of water molecules with the bulk region may be difficult to sample, the convergence and accuracy in free energy perturbation (FEP) calculations can be severely compromised. These problems are further compounded when a reduced system including only the region surrounding the binding site is simulated. To address these issues, we couple molecular dynamics (MD) with grand canonical Monte Carlo (GCMC) simulations to allow the number of water to fluctuate during an alchemical FEP calculation. The atoms in a spherical inner region around the binding pocket are treated explicitly while the influence of the outer region is approximated using the generalized solvent boundary potential (GSBP). At each step during thermodynamic integration, the number of water in the inner region is equilibrated with GCMC and energy data generated with MD is collected. Free energy calculations on camphor binding to a deeply buried pocket in cytochrome P450cam, which causes about seven water molecules to be expelled, are used to test the method. It concluded that solvation free energy calculations with the GCMC/MD method can greatly improve the accuracy of the computed binding free energy compared to simulations with fixed number of water.     

Theoretical calculations: Can Gibbs free energy for intermolecular complexes be predicted efficiently and accurately?

The theoretical study has been performed to refine the procedure for calculations of Gibbs free energy with a relative accuracy of less than 1 kcal/mol. Three benchmark intermolecular complexes are examined via several quantum-chemical methods, including the second-order Moller-Plesset perturbation (MP2), coupled cluster (CCSD(T)), and density functional (BLYP, B3LYP) theories augmented by Dunnings correlation-consistent basis sets. The effects of electron correlation, basis set size, and anharmonicity are systematically analyzed, and the results are compared with available experimental data. The results of the calculations suggest that experimental accuracy can be reached only by extrapolation of MP2 and CCSD(T) total energies to the complete basis set. The contribution of anharmonicity to the zero point energy and TDeltaSint values is fairly small. The new, economic way to reach chemical accuracy in the calculations of the thermodynamic parameters of intermolecular interactions is proposed. In addition, interaction energy (De) and free energy change (DeltaA) for considered species have been evaluated by Carr-Parrinello molecular dynamics (CPMD) simulations and static BLYP-plane wave calculations. The free energy change along the reaction paths were determined by the thermodynamic integration/"Blue Moon Ensemble" technique. Comparison between obtained values, and available experimental and conventional ab initio results has been made. We found that the accuracy of CPMD simulations is affected by several factors, including statistical uncertainty and convergence of constrained forces (TD integration), and the nature of DFT (density functional theory) functional. The results show that CPMD technique is capable of reproducing interaction and free energy with an accuracy of 1 kcal/mol and 2-3 kcal/mol respectively.     

Accurate potential energy surface and quantum reaction rate calculations for the H+CH4-->H2+CH3 reaction

Calculations for the cumulative reaction probability N(E) (for J=0) and the thermal rate constant k(T) of the H+CH(4)-->H(2)+CH(3) reaction are presented. Accurate electronic structure calculations and a converged Shepard-interpolation approach are used to construct a potential energy surface which is specifically designed to allow the precise calculation of k(T) and N(E). Accurate quantum dynamics calculations employing flux correlation functions and multiconfigurational time-dependent Hartree wave packet propagation compute N(E) and k(T) based on this potential energy surface. The present work describes in detail the various convergence test performed to investigate the accuracy of the calculations at each step. These tests demonstrate the predictive power of the present calculations. In addition, approximate approaches for reaction rate calculations are discussed. A quite accurate approximation can be obtained from a potential energy surface which includes only interpolation points on the minimum energy path.     

Theoretical investigation of the structure and coverage of the Si(111)-OCH3 surface

The surface structure, strain energy, and charge profile of the methoxylated Si(111) surface, Si(111)-OCH3, has been studied using quantum mechanics, and the results are compared to those obtained previously for Si(111)-CH3 and Si(111)-C2H5. The calculations indicate that 100% coverage is feasible for Si(111)-OCH3 (similar to the methylated surface), as compared to only approximately 80% coverage for the ethylated surface. These differences can be understood in terms of nearest-neighbor steric and electrostatic interactions. Enthalpy and free energy calculations indicate that the formation of the Si(111)-OCH3 surface from Si(111)-H and methanol is favorable at 300 K. The calculations have also indicated the conditions under which stacking faults can emerge on Si(111)-OCH3, and such conditions are contrasted with the behavior of Si(111)-CH3 and Si(111)-CH2CH3 surfaces, for which stacking faults are calculated to be energetically feasible when etch pits with sufficiently long edges are present on the surface.     

Stability,metastability, and unstability of three-electron-bonded radical anions. A model ab initio theoretical study

The stability of O therefore O, N therefore N, S therefore S, P therefore P, and Si therefore Si three-electron bonds in anionic radicals isoelectronic to dihalogen radical anions is studied by means of ab initio calculations on model systems. The difficulty of generating the dissociation energy profiles of such anions and their rearrangement to neutral species is solved by a practical method which consists of calculating the neutral and anionic energy profiles separately and shifting the curves with respect to each other to match the experimental energy gap between the asymptotes. Here the neutral and anionic reaction profiles are calculated at the CASPT2 and MP2 levels, respectively. The calculations predict that the O therefore O bond is likely to be observed in anions of the type [RO therefore OR](*-), where R is any alkyl substituent or carbon chain. The anion Si(2)H(6)(*-) is found to be a metastable species, with a fair barrier to electron detachment. The barrier is much smaller for N(2)H(4)(*-) and P(2)H(4)(*-), thus precluding experimental observation. However, these species can be stabilized by electron-attractor substituents, the effect of which can be quantitatively estimated by means of the parent anion's diagrams and some fast complementary calculations. An example is given with the [CF(3)HN therefore NHCF(3)](*-) anionic complex.     

Novel properties from experimental charge densities: an application to the zwitterionic neurotransmitter taurine

The charge distribution of taurine (2-aminoethane-sulfonic acid) is revisited by using an orbital-based method that describes the density in a fixed molecular orbital basis with variable orbital occupation numbers. A new neutron data set is also employed to explore whether this improves the deconvolution of thermal motion and charge density. A range of molecular properties that are novel for experimentally determined charge densities are computed, including Weinhold population analysis, Mayer bond orders, and local kinetic energy densities, in addition to charge topological analysis and quantum theory of atoms-in-molecules (QTAIM) integrated properties. The ease with which a distributed multipole analysis can be performed on the fitted density matrix makes it straightforward to compute molecular moments, the lattice energy, and the electrostatic interaction energies of molecules removed from the crystal. Results are compared with high-level (QCISD) gas-phase calculations and band structure calculations employing density functional theory. Finally, the avenues available for extending the range of molecular properties that can be calculated from experimental charge densities still further using this approach are discussed.     

On the use of exact state counting methods in RRKM rate calculations

The use of direct counts of molecular energy level degeneracies, P(E), t internal energy, E, RRKM unimolecular reaction rate calculations is clarified. Highly accurate densities of states for use in practical rate calculations can be obtained by utilizing an energy level “grain” size of 1–5 cm^?1 and by suitable averaging of P(E) over small energy intervals. An easily programmed algorithm for rounding energy levels to multiples of this “grain” size is given and applied to sum and density calculations.     

Computing Alchemical Free Energy Differences with Hamiltonian Replica Exchange Molecular Dynamics (H-REMD) Simulations

Alchemical free energy calculations play a very important role in the field of molecular modeling. Efforts have been made to improve the accuracy and precision of those calculations. One of the efforts is to employ a Hamiltonian replica exchange molecular dynamics (H-REMD) method to enhance conformational sampling. In this paper, we demonstrated that HREMD method not only improves convergence in alchemical free energy calculations but also can be used to compute free energy differences directly via the Free Energy Perturbation (FEP)algorithm. We show a direct mapping between the H-REMD and the usual FEP equations, which are then used directly to compute free energies. The H-REMD alchemical free energy calculation (Replica exchange Free Energy Perturbation, REFEP) was tested on predicting the pK(a) value of the buried Asp26 in thioredoxin. We compare the results of REFEP with TI and regular FEP simulations. REFEP calculations converged faster than those from TI and regular FEP simulations. The final predicted pK(a) value from the H-REMD simulation was also very accurate, only 0.4 pK(a) unit above the experimental value. Utilizing the REFEP algorithm significantly improves conformational sampling, and this in turn improves the convergence of alchemical free energy simulations.     

Reduction of the virtual space for coupled-cluster excitation energies of large molecules and embedded systems

We investigate how the reduction of the virtual space affects coupled-cluster excitation energies at the approximate singles and doubles coupled-cluster level (CC2). In this reduced-virtual-space (RVS) approach, all virtual orbitals above a certain energy threshold are omitted in the correlation calculation. The effects of the RVS approach are assessed by calculations on the two lowest excitation energies of 11 biochromophores using different sizes of the virtual space. Our set of biochromophores consists of common model systems for the chromophores of the photoactive yellow protein, the green fluorescent protein, and rhodopsin. The RVS calculations show that most of the high-lying virtual orbitals can be neglected without significantly affecting the accuracy of the obtained excitation energies. Omitting all virtual orbitals above 50 eV in the correlation calculation introduces errors in the excitation energies that are smaller than 0.1 eV. By using a RVS energy threshold of 50 eV, the CC2 calculations using triple-ζ basis sets (TZVP) on protonated Schiff base retinal are accelerated by a factor of 6. We demonstrate the applicability of the RVS approach by performing CC2/TZVP calculations on the lowest singlet excitation energy of a rhodopsin model consisting of 165 atoms using RVS thresholds between 20 eV and 120 eV. The calculations on the rhodopsin model show that the RVS errors determined in the gas-phase are a very good approximation to the RVS errors in the protein environment. The RVS approach thus renders purely quantum mechanical treatments of chromophores in protein environments feasible and offers an ab initio alternative to quantum mechanics/molecular mechanics separation schemes.     

"In silico" mechanistic studies as predictive tools in microwave-assisted organic synthesis

Computational calculations can be used as a predictive tool in Microwave-Assisted Organic Synthesis (MAOS). A DFT study on Intramolecular Diels-Alder reactions (IMDA) indicated that the activation energy of the reaction and the polarity of the stationary points are two fundamental parameters to determine "a priori" if a reaction can be improved by using microwave irradiation.     

Various energy minima and corresponding fragmentation processes: Alkylsilanes

Fragmentation of alkylsilanes, in particular trimethylethylsilane, were studied by mass-analysed ion kinetic energy (MIKE) and collision-induced decomposition MIKE techniques. Ab initio and semi-empirical molecular orbital calculations were applied to explain the main fragmentation processes. These calculations indicate that more than one minimum can be located on the potential energy surface of a given ground-state molecule ion. These differ from each other mainly in the length of the silicon–carbon bonds. The structures can be adequately described as complexes of a trivalent silyl ion and an alkyl radical. Each of these complexes fragments by the loss of the weakly bound alkyl radical. The calculated energetics of these reactions were found to be in good agreement with the observed energy dependence of the mass spectra.     

Dissociative photoionization mechanism of 1,8-dihydroxyanthraquinone: an experimental and theoretical study

The photoionization and dissociative photoionization mechanism of 1,8-dihydroxyanthraquinone (1,8-DHAQ) have been investigated by infrared laser desorption/tunable synchrotron vacuum ultraviolet photoionization mass spectrometry (IR LD/VUV PIMS) technique and theoretical calculations. Consecutive losses of two carbon monoxides and elimination of hydroxyl group are found to be the major fragmentation channels in low photon energy range. Photoionization efficiency (PIE) spectrum of 1,8-DHAQ was measured in the photon energy range of 8.2-15.0 eV. Adiabatic ionization energy (IE) of 1,8-DAHQ (M) and appearance energies (AEs) of the major fragments (M-CO) (+), (M-C 2O 2) (+), and (M-OH) (+) are determined to be 8.54 +/- 0.05, 10.8 +/- 0.1, 11.0 +/- 0.1, and 13.1 +/- 0.1 eV, respectively, which are in fair agreement with calculated results. The B3LYP method with the 6-31+G(d) basis set was used to study fragmentation of 1,8-DHAQ. Theoretical calculations indicate that five lowest-energy isomers of 1,8-DHAQ cations can coexist by virtue of bond rotation and intramolecular proton transfer. A number of decarbonylation and dehydroxylation processes of 1,8-DHAQ cations are well established.     

Kernel energy method illustrated with peptides

We describe a kernel energy method (KEM) for applying quantum crystallography to large molecules, with an emphasis on the calculation of the molecular energy of peptides. The computational difficulty of representing the system increases only modestly with the number of atoms. The calculations are carried out on modern parallel supercomputers. By adopting the approximation that a full biological molecule can be represented by smaller “kernels” of atoms, the calculations are greatly simplified. Moreover, collections of kernels are, from a computational point of view, well suited for parallel computation. The result is a modest increase in computational time as the number of atoms increases, while retaining the ab initio character of the calculations. We describe a test of our method, and establish its accuracy using 15 different peptides of biological interest. © 2005 Wiley Periodicals, Inc. Int J Quantum Chem, 2005     

Intermolecular interactions in nitrogen-containing aromatic systems

π–π and CH···N interactions are vital in biological systems. In this study, stacking and hydrogen-bonded interactions in pyrazine and triazine dimers were investigated by density functional theory combined with symmetry-adapted perturbation theory (DFT-SAPT) and counterpoise (CP)-corrected supermolecular MP2, SCS-MP2, B3LYP-D and CCSD(T) calculations. All interaction energies were computed using the optimized structures at the CP-corrected SCS/aug-cc-pVDZ level, which gave 1–2 kJ/mol lower interaction energies than the ones computed at the MP2 level. For both dimers, doubly hydrogen-bonded and cross-(displaced) stacked orientations were found to be the lowest energy ones. The reference CCSD(T) calculations favored the former structure in both dimer systems, whereas MP2 and SCS-MP2 located the latter as the lowest energy isomer. In particular, the former was found to be lower in energy than the latter by 2.28 and 1.01 kJ/mol at the CCSD(T)/aug-cc-pVDZ level for pyrazine and triazine, respectively. B3LYP-D produced interaction energies in agreement with the CCSD(T) at the equilibrium geometries, but it overestimates them at the short range and underestimates at the long intermonomer separations. Furthermore, it tends to give smaller equilibrium distances compared to the CCSD(T). DFT-SAPT method was in a good agreement with the reference CCSD(T) calculations. This suggests that DFT-SAPT can be employed to compute the full potential energy surface of these dimers. Moreover, DFT-SAPT calculations showed that the electrostatic and dispersion contributions are the most important energy components stabilizing these dimers. The present study aims to show which theoretical method is the most promising one for the investigation of intermolecular interactions dominated by π–π and CH···N. Therefore, the findings obtained in this study can be used to unravel the structures of nucleic acid bases and other systems stabilized by π–π and CH···N interactions.