Broadening access to small-molecule parameterization with the force field toolkit.

Access to accurate force-field parameters for small molecules is crucial for computational studies of their interactions with proteins. Although a number of general force fields for small molecules exist, e.g., CGenFF, GAFF, and OPLS, they do not cover all common chemical groups and their combinations. The Force Field Toolkit (ffTK) provides a comprehensive graphical interface that streamlines the development of classical parameters for small molecules directly from quantum mechanical (QM) calculations, allowing for force-field generation for almost any chemical group and validation of the fit relative to the target data. ffTK relies on supported external software for the QM calculations, but it can generate the necessary QM input files and parse and analyze the QM output. In previous ffTK versions, support for Gaussian and ORCA QM packages was implemented. Here, we add support for Psi4, an open-source QM package free for all users, thereby broadening user access to ffTK. We also compare the parameter sets obtained with the new ffTK version using Gaussian, ORCA, and Psi4 for three molecules: pyrrolidine, n-propylammonium cation, and chlorobenzene. Despite minor differences between the resulting parameter sets for each compound, most prominently in the dihedral and improper terms, we show that conformational distributions sampled in molecular dynamics simulations using these parameter sets are quite comparable.

Directional ΔG Neural Network (DrΔG-Net): A Modular Neural Network Approach to Binding Free Energy Prediction.

The protein-ligand binding free energy is a central quantity in structure-based computational drug discovery efforts. Although popular alchemical methods provide sound statistical means of computing the binding free energy of a large breadth of systems, they are generally too costly to be applied at the same frequency as end point or ligand-based methods. By contrast, these data-driven approaches are typically fast enough to address thousands of systems but with reduced transferability to unseen systems. We introduce DrΔG-Net (or simply Dragnet), an equivariant graph neural network that can blend ligand-based and protein-ligand data-driven approaches. It is based on a 3D fingerprint representation of the ligand alone and in complex with the protein target. Dragnet is a global scoring function to predict the binding affinity of arbitrary protein-ligand complexes, but can be easily tuned via transfer learning to specific systems or end points, performing similarly to common 2D ligand-based approaches in these tasks. Dragnet is evaluated on a total of 28 validation proteins with a set of congeneric ligands derived from the Binding DB and one custom set extracted from the ChEMBL Database. In general, a handful of experimental binding affinities are sufficient to optimize the scoring function for a particular protein and ligand scaffold. When not available, predictions from physics-based methods such as absolute free energy perturbation can be used for the transfer learning tuning of Dragnet. Furthermore, we use our data to illustrate the present limitations of data-driven modeling of binding free energy predictions.

Approximating large-basis coupled-cluster theory vibrational frequencies using focal-point approximations.

The focal-point approximation can be used to estimate a high-accuracy, slow quantum chemistry computation by combining several lower-accuracy, faster computations. We examine the performance of focal-point methods by combining second-order Møller-Plesset perturbation theory (MP2) with coupled-cluster singles, doubles, and perturbative triples [CCSD(T)] for the calculation of harmonic frequencies and that of fundamental frequencies using second-order vibrational perturbation theory (VPT2). In contrast to standard CCSD(T), the focal-point CCSD(T) method approaches the complete basis set (CBS) limit with only triple-ζ basis sets for the coupled-cluster portion of the computation. The predicted harmonic and fundamental frequencies were compared with the experimental values for a set of 20 molecules containing up to six atoms. The focal-point method combining CCSD(T)/aug-cc-pV(T + d)Z with CBS-extrapolated MP2 has mean absolute errors vs experiment of only 7.3 cm-1 for the fundamental frequencies, which are essentially the same as the mean absolute error for CCSD(T) extrapolated to the CBS limit using the aug-cc-pV(Q + d)Z and aug-cc-pV(5 + d)Z basis sets. However, for H2O, the focal-point procedure requires only 3% of the computation time as the extrapolated CCSD(T) result, and the cost savings will grow for larger molecules.

A quantum chemical interaction energy dataset for accurately modeling protein-ligand interactions.

Fast and accurate calculation of intermolecular interaction energies is desirable for understanding many chemical and biological processes, including the binding of small molecules to proteins. The Splinter ["Symmetry-adapted perturbation theory (SAPT0) protein-ligand interaction"] dataset has been created to facilitate the development and improvement of methods for performing such calculations. Molecular fragments representing commonly found substructures in proteins and small-molecule ligands were paired into >9000 unique dimers, assembled into numerous configurations using an approach designed to adequately cover the breadth of the dimers' potential energy surfaces while enhancing sampling in favorable regions. ~1.5 million configurations of these dimers were randomly generated, and a structurally diverse subset of these were minimized to obtain an additional ~80 thousand local and global minima. For all >1.6 million configurations, SAPT0 calculations were performed with two basis sets to complete the dataset. It is expected that Splinter will be a useful benchmark dataset for training and testing various methods for the calculation of intermolecular interaction energies.

Benchmark coupled-cluster lattice energy of crystalline benzene and assessment of multi-level approximations in the many-body expansion.

The many-body expansion (MBE) is promising for the efficient, parallel computation of lattice energies in organic crystals. Very high accuracy should be achievable by employing coupled-cluster singles, doubles, and perturbative triples at the complete basis set limit [CCSD(T)/CBS] for the dimers, trimers, and potentially tetramers resulting from the MBE, but such a brute-force approach seems impractical for crystals of all but the smallest molecules. Here, we investigate hybrid or multi-level approaches that employ CCSD(T)/CBS only for the closest dimers and trimers and utilize much faster methods like Møller-Plesset perturbation theory (MP2) for more distant dimers and trimers. For trimers, MP2 is supplemented with the Axilrod-Teller-Muto (ATM) model of three-body dispersion. MP2(+ATM) is shown to be a very effective replacement for CCSD(T)/CBS for all but the closest dimers and trimers. A limited investigation of tetramers using CCSD(T)/CBS suggests that the four-body contribution is entirely negligible. The large set of CCSD(T)/CBS dimer and trimer data should be valuable in benchmarking approximate methods for molecular crystals and allows us to see that a literature estimate of the core-valence contribution of the closest dimers to the lattice energy using just MP2 was overbinding by 0.5 kJ mol-1, and an estimate of the three-body contribution from the closest trimers using the T0 approximation in local CCSD(T) was underbinding by 0.7 kJ mol-1. Our CCSD(T)/CBS best estimate of the 0 K lattice energy is -54.01 kJ mol-1, compared to an estimated experimental value of -55.3 ± 2.2 kJ mol-1.

A quantitative assessment of deformation energy in intermolecular interactions: How important is it?

Dimer interaction energies have been well studied in computational chemistry, but they can offer an incomplete understanding of molecular binding depending on the system. In the current study, we present a dataset of focal-point coupled-cluster interaction and deformation energies (summing to binding energies, De) of 28 organic molecular dimers. We use these highly accurate energies to evaluate ten density functional approximations for their accuracy. The best performing method (with a double-ζ basis set), B97M-D3BJ, is then used to calculate the binding energies of 104 organic dimers, and we analyze the influence of the nature and strength of interaction on deformation energies. Deformation energies can be as large as 50% of the dimer interaction energy, especially when hydrogen bonding is present. In most cases, two or more hydrogen bonds present in a dimer correspond to an interaction energy of -10 to -25 kcal mol-1, allowing a deformation energy above 1 kcal mol-1 (and up to 9.5 kcal mol-1). A lack of hydrogen bonding usually restricts the deformation energy to below 1 kcal mol-1 due to the weaker interaction energy.

Assessment of three-body dispersion models against coupled-cluster benchmarks for crystalline benzene, carbon dioxide, and triazine.

To study the contribution of three-body dispersion to crystal lattice energies, we compute the three-body contributions to the lattice energies for crystalline benzene, carbon dioxide, and triazine using various computational methods. We show that these contributions converge quickly as the intermolecular distances between the monomers grow. In particular, the smallest value among the three pairwise intermonomer closest-contact distances, Rmin, shows a strong correlation with the three-body contribution to the lattice energy, and, here, the largest of the closest-contact distances, Rmax, serves as a cutoff criterion to limit the number of trimers to be considered. We considered all trimers up to Rmax=15Å. The trimers with Rmin<4Å contribute 90.4%, 90.6%, and 93.9% of the total three-body contributions for crystalline benzene, carbon dioxide, and triazine, respectively, for the coupled-cluster singles, doubles, and perturbative triples [CCSD(T)] method. For trimers with Rmin>4Å, the second-order Møller-Plesset perturbation theory (MP2) supplemented with the Axilrod-Teller-Muto (ATM) three-body dispersion correction reproduces the CCSD(T) values for the cumulative three-body contributions with errors of less than 0.1 kJ mol-1. Moreover, three-body contributions are converged within 0.15 kJ mol-1 by Rmax=10Å. From these results, it appears that in molecular crystals where dispersion dominates the three-body contribution to the lattice energy, the trimers with Rmin>4Å can be computed with the MP2+ATM method to reduce the computational cost, and those with Rmax>10Å appear to be basically negligible.

Benchmarking two-body contributions to crystal lattice energies and a range-dependent assessment of approximate methods.

Using the many-body expansion to predict crystal lattice energies (CLEs), a pleasantly parallel process, allows for flexibility in the choice of theoretical methods. Benchmark-level two-body contributions to CLEs of 23 molecular crystals have been computed using interaction energies of dimers with minimum inter-monomer separations (i.e., closest contact distances) up to 30 Å. In a search for ways to reduce the computational expense of calculating accurate CLEs, we have computed these two-body contributions with 15 different quantum chemical levels of theory and compared these energies to those computed with coupled-cluster in the complete basis set (CBS) limit. Interaction energies of the more distant dimers are easier to compute accurately and several of the methods tested are suitable as replacements for coupled-cluster through perturbative triples for all but the closest dimers. For our dataset, sub-kJ mol-1 accuracy can be obtained when calculating two-body interaction energies of dimers with separations shorter than 4 Å with coupled-cluster with single, double, and perturbative triple excitations/CBS and dimers with separations longer than 4 Å with MP2.5/aug-cc-pVDZ, among other schemes, reducing the number of dimers to be computed with coupled-cluster by as much as 98%.

Range-dependence of two-body intermolecular interactions and their energy components in molecular crystals.

Routinely assessing the stability of molecular crystals with high accuracy remains an open challenge in the computational sciences. The many-body expansion decomposes computation of the crystal lattice energy into an embarrassingly parallel collection of computations over molecular dimers, trimers, and so forth, making quantum chemistry techniques tractable for many crystals of small organic molecules. By examining the range-dependence of different types of energetic contributions to the crystal lattice energy, we can glean qualitative understanding of solid-state intermolecular interactions as well as practical, exploitable reductions in the number of computations required for accurate energies. Here, we assess the range-dependent character of two-body interactions of 24 small organic molecular crystals by using the physically interpretable components from symmetry-adapted perturbation theory (electrostatics, exchange-repulsion, induction/polarization, and London dispersion). We also examine correlations between the convergence rates of electrostatics and London dispersion terms with molecular dipole moments and polarizabilities, to provide guidance for estimating convergence rates in other molecular crystals.

Quantum Chemistry Common Driver and Databases (QCDB) and Quantum Chemistry Engine (QCEngine): Automation and interoperability among computational chemistry programs.

Community efforts in the computational molecular sciences (CMS) are evolving toward modular, open, and interoperable interfaces that work with existing community codes to provide more functionality and composability than could be achieved with a single program. The Quantum Chemistry Common Driver and Databases (QCDB) project provides such capability through an application programming interface (API) that facilitates interoperability across multiple quantum chemistry software packages. In tandem with the Molecular Sciences Software Institute and their Quantum Chemistry Archive ecosystem, the unique functionalities of several CMS programs are integrated, including CFOUR, GAMESS, NWChem, OpenMM, Psi4, Qcore, TeraChem, and Turbomole, to provide common computational functions, i.e., energy, gradient, and Hessian computations as well as molecular properties such as atomic charges and vibrational frequency analysis. Both standard users and power users benefit from adopting these APIs as they lower the language barrier of input styles and enable a standard layout of variables and data. These designs allow end-to-end interoperable programming of complex computations and provide best practices options by default.

Cartesian message passing neural networks for directional properties: Fast and transferable atomic multipoles.

The message passing neural network (MPNN) framework is a promising tool for modeling atomic properties but is, until recently, incompatible with directional properties, such as Cartesian tensors. We propose a modified Cartesian MPNN (CMPNN) suitable for predicting atom-centered multipoles, an essential component of ab initio force fields. The efficacy of this model is demonstrated on a newly developed dataset consisting of 46 623 chemical structures and corresponding high-quality atomic multipoles, which was deposited into the publicly available Molecular Sciences Software Institute QCArchive server. We show that the CMPNN accurately predicts atom-centered charges, dipoles, and quadrupoles and that errors in the predicted atomic multipoles have a negligible effect on multipole-multipole electrostatic energies. The CMPNN is accurate enough to model conformational dependencies of a molecule's electronic structure. This opens up the possibility of recomputing atomic multipoles on the fly throughout a simulation in which they might exhibit strong conformational dependence.

AP-Net: An atomic-pairwise neural network for smooth and transferable interaction potentials.

Intermolecular interactions are critical to many chemical phenomena, but their accurate computation using ab initio methods is often limited by computational cost. The recent emergence of machine learning (ML) potentials may be a promising alternative. Useful ML models should not only estimate accurate interaction energies but also predict smooth and asymptotically correct potential energy surfaces. However, existing ML models are not guaranteed to obey these constraints. Indeed, systemic deficiencies are apparent in the predictions of our previous hydrogen-bond model as well as the popular ANI-1X model, which we attribute to the use of an atomic energy partition. As a solution, we propose an alternative atomic-pairwise framework specifically for intermolecular ML potentials, and we introduce AP-Net-a neural network model for interaction energies. The AP-Net model is developed using this physically motivated atomic-pairwise paradigm and also exploits the interpretability of symmetry adapted perturbation theory (SAPT). We show that in contrast to other models, AP-Net produces smooth, physically meaningful intermolecular potentials exhibiting correct asymptotic behavior. Initially trained on only a limited number of mostly hydrogen-bonded dimers, AP-Net makes accurate predictions across the chemically diverse S66x8 dataset, demonstrating significant transferability. On a test set including experimental hydrogen-bonded dimers, AP-Net predicts total interaction energies with a mean absolute error of 0.37 kcal mol-1, reducing errors by a factor of 2-5 across SAPT components from previous neural network potentials. The pairwise interaction energies of the model are physically interpretable, and an investigation of predicted electrostatic energies suggests that the model "learns" the physics of hydrogen-bonded interactions.

Psi4 1.4: Open-source software for high-throughput quantum chemistry.

PSI4 is a free and open-source ab initio electronic structure program providing implementations of Hartree-Fock, density functional theory, many-body perturbation theory, configuration interaction, density cumulant theory, symmetry-adapted perturbation theory, and coupled-cluster theory. Most of the methods are quite efficient, thanks to density fitting and multi-core parallelism. The program is a hybrid of C++ and Python, and calculations may be run with very simple text files or using the Python API, facilitating post-processing and complex workflows; method developers also have access to most of PSI4's core functionalities via Python. Job specification may be passed using The Molecular Sciences Software Institute (MolSSI) QCSCHEMA data format, facilitating interoperability. A rewrite of our top-level computation driver, and concomitant adoption of the MolSSI QCARCHIVE INFRASTRUCTURE project, makes the latest version of PSI4 well suited to distributed computation of large numbers of independent tasks. The project has fostered the development of independent software components that may be reused in other quantum chemistry programs.