cclib 2.0: An updated architecture for interoperable computational chemistry.

Interoperability in computational chemistry is elusive, impeded by the independent development of software packages and idiosyncratic nature of their output files. The cclib library was introduced in 2006 as an attempt to improve this situation by providing a consistent interface to the results of various quantum chemistry programs. The shared API across programs enabled by cclib has allowed users to focus on results as opposed to output and to combine data from multiple programs or develop generic downstream tools. Initial development, however, did not anticipate the rapid progress of computational capabilities, novel methods, and new programs; nor did it foresee the growing need for customizability. Here, we recount this history and present cclib 2, focused on extensibility and modularity. We also introduce recent design pivots-the formalization of cclib's intermediate data representation as a tree-based structure, a new combinator-based parser organization, and parsed chemical properties as extensible objects.

Supramolecular Complexes of Plant Neurotoxin Veratridine with Cyclodextrins and Their Antidote-like Effect on Neuro-2a Cell Viability.

Veratridine (VTD) is a plant neurotoxin that acts by blocking the voltage-gated sodium channels (VGSC) of cell membranes. Symptoms of VTD intoxication include intense nausea, hypotension, arrhythmia, and loss of consciousness. The treatment for the intoxication is mainly focused on treating the symptoms, meaning there is no specific antidote against VTD. In this pursuit, we were interested in studying the molecular interactions of VTD with cyclodextrins (CDs). CDs are supramolecular macrocycles with the ability to form host-guest inclusion complexes (ICs) inside their hydrophobic cavity. Since VTD is a lipid-soluble alkaloid, we hypothesized that it could form stable inclusion complexes with different types of CDs, resulting in changes to its physicochemical properties. In this investigation, we studied the interaction of VTD with ß-CD, γ-CD and sulfobutyl ether ß-CD (SBCD) by isothermal titration calorimetry (ITC) and nuclear magnetic resonance (NMR) spectroscopy. Docking and molecular dynamics studies confirmed the most stable configuration for the inclusion complexes. Finally, with an interest in understanding the effects of the VTD/CD molecular interactions, we performed cell-based assays (CBAs) on Neuro-2a cells. Our findings reveal that the use of different amounts of CDs has an antidote-like concentration-dependent effect on the cells, significantly increasing cell viability and thus opening opportunities for novel research on applications of CDs and VTD.

The aug-cc-pVTZ-J basis set for the p-block fourth-row elements Ga, Ge, As, Se, and Br.

The aug-cc-pVTZ-J basis set family is extended to include the fourth-row p-block elements Ga, Ge, As, Se, and Br. We use the established approach outlined by Sauer and coworkers (J. Chem. Phys. 115, 1324 [2001], J. Chem. Phys. 133, 054308 [2010], J. Chem. Theory Comput. 7, 4070 [2011], and J. Chem. Theory Comput. 7, 4077 [2011]) where the completely uncontracted aug-cc-pVTZ basis set is saturated with tight s-, p-, d-, and f-functions to form the aug-cc-pVTZ-Juc basis set for the tested elements. The saturation is carried out on the simplest hydrides possible for the tested elements GaH, GeH4 , AsH3 , H2 Se, and HBr until an improvement is less than 0.01% for all s-, p-, and d-functions added. f-Functions are added to an improvement less than or equal to 1.0% due to the computational expense these functions add. The saturated aug-cc-pVTZ-Juc (26s16p12d5f) is then recontracted using the molecular orbital coefficients from self-consistent field calculations on the simple hydrides to improve computational efficiency. During contraction of the basis set, we observe that the linear hydrogen bromide molecule has a slower convergence than the other tested molecules which sets a limit on the accuracy obtained. All calculations with the contracted aug-cc-pVTZ-J [17s10p7d5f] gives results that are within 1.0% of the uncontracted results at considerable computational savings.

Cost-Effective Potential for Accurate Polarizable Embedding Calculations in Protein Environments.

The fragment-based polarizable embedding (PE) model combined with an appropriate electronic structure method constitutes a highly efficient and accurate multiscale approach for computing spectroscopic properties of a central moiety including effects from its molecular environment through an embedding potential. There is, however, a comparatively high computational overhead associated with the computation of the embedding potential, which is derived from first-principles calculations on individual fragments of the environment. To reduce the computational cost associated with the calculation of embedding potential parameters, we developed a set of amino acid-specific transferable parameters tailored for large-scale PE-based calculations that include proteins. The amino acid-based parameters are obtained by simultaneously fitting to a set of reference electric potentials based on structures derived from a backbone-dependent rotamer library. The developed cost-effective polarizable protein potential (CP3) consists of atom-centered charges and isotropic dipole-dipole polarizabilities of the standard amino acids. In terms of reproduction of electric potentials, the CP3 is shown to perform consistently and with acceptable accuracy across both small tripeptide test systems and larger proteins. We show, through applications on realistic protein systems, that acceptable accuracy can be obtained by using a pure CP3 representation of the protein environment, thus altogether omitting the cost associated with the calculation of embedding potential parameters. High accuracy comparable to that of the full fragment-based approach can be achieved through a mixed description where the CP3 is used only to describe amino acids beyond a threshold distance from the central quantum part.

Improving solvation energy predictions using the SMD solvation method and semiempirical electronic structure methods.

The PM6 implementation in the GAMESS program is extended to elements requiring d-integrals and interfaced with the conducter-like polarized continuum model of solvation, including gradients. The accuracy of aqueous solvation energies computed using AM1, PM3, PM6, and DFT tight binding (DFTB) and the Solvation Model Density (SMD) continuum solvation model is tested using the Minnesota Solvation Database data set. The errors in SMD solvation energies predicted using Neglect of Diatomic Differential Overlap (NDDO)-based methods are considerably larger than when using density functional theory (DFT) and HF, with root mean square error (RMSE) values of 3.4-5.9 (neutrals) and 6-15 kcal/mol (ions) compared to 2.4 and ∼5 kcal/mol for HF/6-31G(d). For the NDDO-based methods, the errors are especially large for cations and considerably higher than the corresponding conductor-like screening model results, which suggests that the NDDO/SMD results can be improved by re-parameterizing the SMD parameters focusing on ions. We found that the best results are obtained by changing only the radii for hydrogen, carbon, oxygen, nitrogen, and sulfur, and this leads to RMSE values for PM3 (neutrals: 2.8/ions: ∼5 kcal/mol), PM6 (4.7/∼5 kcal/mol), and DFTB (3.9/∼5 kcal/mol) that are more comparable to HF/6-31G(d) (2.4/∼5 kcal/mol). Although the radii are optimized to reproduce aqueous solvation energies, they also lead more accurate predictions for other polar solvents such as dimethyl sulfoxide, acetonitrile, and methanol, while the improvements for non-polar solvents are negligible.

Structure-guided approach identifies a novel class of HIV-1 ribonuclease H inhibitors: binding mode insights through magnesium complexation and site-directed mutagenesis studies.

Persistent HIV infection requires lifelong treatment and among the 2.1 million new HIV infections that occur every year there is an increased rate of transmitted drug-resistant mutations. This fact requires a constant and timely effort in order to identify and develop new HIV inhibitors with innovative mechanisms. The HIV-1 reverse transcriptase (RT) associated ribonuclease H (RNase H) is the only viral encoded enzyme that still lacks an efficient inhibitor despite the fact that it is a well-validated target whose functional abrogation compromises viral infectivity. Identification of new drugs is a long and expensive process that can be speeded up by in silico methods. In the present study, a structure-guided screening is coupled with a similarity-based search on the Specs database to identify a new class of HIV-1 RNase H inhibitors. Out of the 45 compounds selected for experimental testing, 15 inhibited the RNase H function below 100 µM with three hits exhibiting IC50 values <10 µM. The most active compound, AA, inhibits HIV-1 RNase H with an IC50 of 5.1 µM and exhibits a Mg-independent mode of inhibition. Site-directed mutagenesis studies provide valuable insight into the binding mode of newly identified compounds; for instance, compound AA involves extensive interactions with a lipophilic pocket formed by Ala502, Lys503, and Trp (406, 426 and 535) and polar interactions with Arg557 and the highly conserved RNase H primer-grip residue Asn474. The structural insights obtained from this work provide the bases for further lead optimization.

Relative Ligand-Binding Free Energies Calculated from Multiple Short QM/MM MD Simulations.

We have devised a new efficient approach to compute combined quantum mechanical (QM) and molecular mechanical (MM, i.e. QM/MM) ligand-binding relative free energies. Our method employs the reference-potential approach with free-energy perturbation both at the MM level (between the two ligands) and from MM to QM/MM (for each ligand). To ensure that converged results are obtained for the MM → QM/MM perturbations, explicit QM/MM molecular dynamics (MD) simulations are performed with two intermediate mixed states. To speed up the calculations, we utilize the fact that the phase space can be extensively sampled at the MM level. Therefore, we run many short QM/MM MD simulations started from snapshots of the MM simulations, instead of a single long simulation. As a test case, we study the binding of nine cyclic carboxylate ligands to the octa-acid deep cavitand. Only the ligand is in the QM system, treated with the semiempirical PM6-DH+ method. We show that for eight of the ligands, we obtain well converged results with short MD simulations (1-15 ps). However, in one case, the convergence is slower (∼50 ps) owing to a mismatch between the conformational preferences of the MM and QM/MM potentials. We test the effect of initial minimization, the need of equilibration, and how many independent simulations are needed to reach a certain precision. The results show that the present approach is about four times faster than using standard MM → QM/MM free-energy perturbations with the same accuracy and precision.

Modeling Electronic Circular Dichroism within the Polarizable Embedding Approach.

We present a systematic investigation of the key components needed to model single chromophore electronic circular dichroism (ECD) within the polarizable embedding (PE) approach. By relying on accurate forms of the embedding potential, where especially the inclusion of local field effects are in focus, we show that qualitative agreement between rotatory strength parameters calculated by full quantum mechanical calculations and the more efficient embedding calculations can be obtained. An important aspect in the computation of reliable absorption parameters is the need for conformational sampling. We show that a significant number of snapshots are needed to avoid artifacts in the calculated electronic circular dichroism parameters due to insufficient configurational sampling, thus highlighting the efficiency of the PE model.

An averaged polarizable potential for multiscale modeling in phospholipid membranes.

A set of average atom-centered charges and polarizabilities has been developed for three types of phospholipids for use in polarizable embedding calculations. The lipids investigated are 1,2-dimyristoyl-sn-glycero-3-phosphocholine, 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine, and 1-palmitoyl-2-oleoyl-sn-glycerol-3-phospho-L-serine given their common use both in experimental and computational studies. The charges, and to a lesser extent the polarizabilities, are found to depend strongly on the molecular conformation of the lipids. Furthermore, the importance of explicit polarization is underlined for the description of larger assemblies of lipids, that is, membranes. In conclusion, we find that specially developed polarizable parameters are needed for embedding calculations in membranes, while common non-polarizable point-charge force fields usually perform well enough for structural and dynamical studies. © 2017 Wiley Periodicals, Inc.

Mapping Interaction Energies in Chorismate Mutase with the Fragment Molecular Orbital Method.

The Claisen rearrangement of chorismate to prephenate is mapped across the entire reaction pathway using the fragment molecular orbital (FMO) method. Three basis sets (6-31G(d), cc-pVDZ, and pcseg-1) are studied to provide guidance toward obtaining high accuracy with the FMO method on such systems. Using a fragmentation scheme of one residue per fragment, the FMO method using the 6-31G(d) basis set and second-order Møller-Plesset perturbation theory (MP2) with the hybrid orbital projection fragmentation scheme provides the most reliable results across the entire reaction pathway. Calculations using the multilayer FMO method are performed and shown to be in agreement with single-layer calculations in all cases with differences of less than one kilocalorie per mole for all tested basis set combinations along the entire reaction path. The use of restricted Hartree-Fock for the lower-level layer and MP2 for the higher-level layer gives the most consistent results when using the same basis set for both layers. Pair interaction energy decomposition analysis calculations confirm that electrostatic interactions are the predominant force between three key arginine residues and chorismate and that dispersion and charge transfer interactions in the binding pocket also play a role in the local chemistry of the reaction.

Computational Approach to Evaluation of Optical Properties of Membrane Probes.

Computed optical properties of membrane probes are typically evaluated in the gas phase, i.e. neglecting the influence of the membrane. In this study, we examine how and to what extent a membrane influences the one- and two-photon absorption (1PA and 2PA, respectively) properties for a number of cholesterol analogs and thereby also evaluate the validity of the common gas phase approach. The membrane is modeled using the polarizable embedding scheme both with and without the effective external field extension of the polarizable embedding model. The shifts in excitation energies and 1PA oscillator strengths compared to the gas phase are relatively small, while the 2PA cross section is more affected. The electric field inside the membrane induces a larger change in the permanent electric dipole moment upon excitation of the analogs compared to the gas phase, which leads to an almost 2-fold increase in the 2PA cross section for one cholesterol analog. The relative trends observed in the membrane are the same as in the gas phase, and the use of gas phase calculations for qualitative comparison and design of cholesterol membrane probes is thus a useful and computationally efficient strategy.

Automated Fragmentation Polarizable Embedding Density Functional Theory (PE-DFT) Calculations of Nuclear Magnetic Resonance (NMR) Shielding Constants of Proteins with Application to Chemical Shift Predictions.

Full-protein nuclear magnetic resonance (NMR) shielding constants based on ab initio calculations are desirable, because they can assist in elucidating protein structures from NMR experiments. In this work, we present NMR shielding constants computed using a new automated fragmentation (J. Phys. Chem. B 2009, 113, 10380-10388) approach in the framework of polarizable embedding density functional theory. We extend our previous work to give both basis set recommendations and comment on how large the quantum mechanical region should be to successfully compute 13C NMR shielding constants that are comparable with experiment. The introduction of a probabilistic linear regression model allows us to substantially reduce the number of snapshots that are needed to make comparisons with experiment. This approach is further improved by augmenting snapshot selection with chemical shift predictions by which we can obtain a representative subset of snapshots that gives the smallest predicted error, compared to experiment. Finally, we use this subset of snapshots to calculate the NMR shielding constants at the PE-KT3/pcSseg-2 level of theory for all atoms in the protein GB3.

Computational Approach for Studying Optical Properties of DNA Systems in Solution.

In this paper we present a study of the methodological aspects regarding calculations of optical properties for DNA systems in solution. Our computational approach will be built upon a fully polarizable QM/MM/Continuum model within a damped linear response theory framework. In this approach the environment is given a highly advanced description in terms of the electrostatic potential through the polarizable embedding model. Furthermore, bulk solvent effects are included in an efficient manner through a conductor-like screening model. With the aim of reducing the computational cost we develop a set of averaged partial charges and distributed isotropic dipole-dipole polarizabilities for DNA suitable for describing the classical region in ground-state and excited-state calculations. Calculations of the UV-spectrum of the 2-aminopurine optical probe embedded in a DNA double helical structure are presented. We show that inclusion of polarizabilities in the embedding potential stemming from the DNA double helix is of crucial importance, while the water cluster surrounding the DNA system is well represented using a continuum approach.

Electronic Energy Transfer in Polarizable Heterogeneous Environments: A Systematic Investigation of Different Quantum Chemical Approaches.

Theoretical prediction of transport and optical properties of protein-pigment complexes is of significant importance when aiming at understanding the structure-function relationship in such systems. Electronic energy transfer (EET) couplings represent a key property in this respect since such couplings provide important insight into the strength of interaction between photoactive pigments in protein-pigment complexes. Recently, attention has been payed to how the environment modifies or even controls the electronic couplings. To enable such theoretical predictions, a fully polarizable embedding model has been suggested (Curutchet, C., et al. J. Chem. Theory Comput., 2009, 5, 1838-1848). In this work, we further develop this computational model by extending it with an ab initio derived polarizable force field including higher-order multipole moments. We use this extended model to systematically examine three different ways of obtaining EET couplings in a heterogeneous medium ranging from use of the exact transition density to a point-dipole approximation. Several interesting observations are made, for example, the explicit use of transition densities in the calculation of the electronic couplings, and also when including the explicit environment contribution, can be replaced by a much simpler transition point charge description without comprising the quality of the model predictions.

Polarizable density embedding: a new QM/QM/MM-based computational strategy.

We present a new QM/QM/MM-based model for calculating molecular properties and excited states of solute-solvent systems. We denote this new approach the polarizable density embedding (PDE) model, and it represents an extension of our previously developed polarizable embedding (PE) strategy. The PDE model is a focused computational approach in which a core region of the system studied is represented by a quantum-chemical method, whereas the environment is divided into two other regions: an inner and an outer region. Molecules belonging to the inner region are described by their exact densities, whereas molecules in the outer region are treated using a multipole expansion. In addition, all molecules in the environment are assigned distributed polarizabilities in order to account for induction effects. The joint effects of the inner and outer regions on the quantum-mechanical core part of the system is formulated using an embedding potential. The PDE model is illustrated for a set of dimers (interaction energy calculations) as well as for the calculation of electronic excitation energies, showing promising results.

A third-generation dispersion and third-generation hydrogen bonding corrected PM6 method: PM6-D3H+.

We present new dispersion and hydrogen bond corrections to the PM6 method, PM6-D3H+, and its implementation in the GAMESS program. The method combines the DFT-D3 dispersion correction by Grimme et al. with a modified version of the H+ hydrogen bond correction by Korth. Overall, the interaction energy of PM6-D3H+ is very similar to PM6-DH2 and PM6-DH+, with RMSD and MAD values within 0.02 kcal/mol of one another. The main difference is that the geometry optimizations of 88 complexes result in 82, 6, 0, and 0 geometries with 0, 1, 2, and 3 or more imaginary frequencies using PM6-D3H+ implemented in GAMESS, while the corresponding numbers for PM6-DH+ implemented in MOPAC are 54, 17, 15, and 2. The PM6-D3H+ method as implemented in GAMESS offers an attractive alternative to PM6-DH+ in MOPAC in cases where the LBFGS optimizer must be used and a vibrational analysis is needed, e.g., when computing vibrational free energies. While the GAMESS implementation is up to 10 times slower for geometry optimizations of proteins in bulk solvent, compared to MOPAC, it is sufficiently fast to make geometry optimizations of small proteins practically feasible.

Inhibitor ranking through QM based chelation calculations for virtual screening of HIV-1 RNase H inhibition.

Quantum mechanical (QM) calculations have been used to predict the binding affinity of a set of ligands towards HIV-1 RT associated RNase H (RNH). The QM based chelation calculations show improved binding affinity prediction for the inhibitors compared to using an empirical scoring function. Furthermore, full protein fragment molecular orbital (FMO) calculations were conducted and subsequently analysed for individual residue stabilization/destabilization energy contributions to the overall binding affinity in order to better understand the true and false predictions. After a successful assessment of the methods based on the use of a training set of molecules, QM based chelation calculations were used as filter in virtual screening of compounds in the ZINC database. By this, we find, compared to regular docking, QM based chelation calculations to significantly reduce the large number of false positives. Thus, the computational models tested in this study could be useful as high throughput filters for searching HIV-1 RNase H active-site molecules in the virtual screening process.

Hybrid RHF/MP2 geometry optimizations with the effective fragment molecular orbital method.

The frozen domain effective fragment molecular orbital method is extended to allow for the treatment of a single fragment at the MP2 level of theory. The approach is applied to the conversion of chorismate to prephenate by Chorismate Mutase, where the substrate is treated at the MP2 level of theory while the rest of the system is treated at the RHF level. MP2 geometry optimization is found to lower the barrier by up to 3.5 kcal/mol compared to RHF optimzations and ONIOM energy refinement and leads to a smoother convergence with respect to the basis set for the reaction profile. For double zeta basis sets the increase in CPU time relative to RHF is roughly a factor of two.

Nuclear Magnetic Shielding Constants from Quantum Mechanical/Molecular Mechanical Calculations Using Polarizable Embedding: Role of the Embedding Potential.

We present NMR shielding constants obtained through quantum mechanical/molecular mechanical (QM/MM) embedding calculations. Contrary to previous reports, we show that a relatively small QM region is sufficient, provided that a high-quality embedding potential is used. The calculated averaged NMR shielding constants of both acrolein and acetone solvated in water are based on a number of snapshots extracted from classical molecular dynamics simulations. We focus on the carbonyl chromophore in both molecules, which shows large solvation effects, and we study the convergence of shielding constants with respect to total system size and size of the QM region. By using a high-quality embedding potential over standard point charge potentials, we show that the QM region can be made at least 2 Å smaller without any loss of quality, which makes calculations on ensembles tractable by conventional density functional theory calculations.