Discovering a Transferable Charge Assignment Model using Machine Learning
Partial atomic charge assignment is of immense practical value to force field parametrization, molecular docking, and cheminformatics. Machine learning has emerged as a powerful tool for modeling chemistry at unprecedented computational speeds given ground-truth values, but for the task of charge assignment, the choice of ground-truth may not be obvious. In this letter, we use machine learning to discover a charge model by training a neural network to molecular dipole moments using a large, diverse set of CHNO molecular conformations. The new model, called Affordable Charge Assignment (ACA), is computationally inexpensive and predicts dipoles of out-of-sample molecules accurately. Furthermore, dipole-inferred ACA charges are transferable to dipole and even quadrupole moments of much larger molecules than those used for training. We apply ACA to long dynamical trajectories of biomolecules and successfully produce their infrared spectra. Additionally, we compare ACA with existing charge models and find that ACA assigns similar charges to Charge Model 5, but with a greatly reduced computational cost.
<!>Details of Charge Assignment
<p>Electrostatic interactions contribute strongly to the forces within and between molecules.</p><p>These interactions depend on the charge density field ⇢(r), which is computationally demanding to compute. Simplified models of the charge density, such as atom-centered monopoles, are commonly employed. These partial atomic charges result in faster computation as well as provide a qualitative understanding of the underlying chemistry. [1][2][3][4] However, the decomposition of charge density into atomic charges is, by itself, an ambiguous task. Additional principles are necessary to make the charge assignment task well-defined. Here we show that a Machine Learning model, trained only on the dipole moments of small molecules, discovers a charge model that is transferable to quadrupole predictions and extensible to much larger molecules.</p><p>Existing popular charge models have also been designed to reproduce observables of the electrostatic potential. The Merz-Singh-Kollman (MSK) [5,6] charge model exactly replicates the dipole moment and approximates the electrostatic potential on many points surrounding the molecule, resulting in high-quality electrostatic properties exterior to the molecule. However, MSK suffers from basis set sensitivity, particularly for "buried atoms" located inside large molecules. [7][8][9] Charge model 5 (CM5) [8] is an extension of Hirshfeld analysis, [10] with additional parametrization in order to approximately reproduce ab initio and experimental dipoles of 614 gas-phase dipoles. Unlike MSK, Hirshfeld and CM5 are nearly independent of basis set. [9] This insensitivity allows CM5 to use a single set of model parameters. The corresponding tradeoff is that its charges do not reproduce electrostatic fields as well as MSK.</p><p>A limitation of these conventional charge models is that they require expensive ab initio calculation, which can be computationally impractical, especially for large molecules, long time scales, or systems exhibiting great chemical diversity. Recent advances in machine learning (ML) have demonstrated great potential to build quantum chemistry models with ab initio-level accuracy while bypassing ab initio costs. [11] Trained to reference datasets, ML models can predict energies, forces, and other molecular properties. [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] They have been used to discover materials [28][29][30][31][32][33][34][35][36][37] and study dynamical processes such as charge and exciton transfer. [38][39][40][41] Most related to this work are ML models of existing charge models, [9,[42][43][44] which are orders of magnitude faster than ab initio calculation. Here, we show that ML is able to go beyond emulation and discover a charge model that closely reproduces electrostatic properties by training directly to the dipole moment.</p><p>In this letter, we use HIP-NN (Hierarchically Interacting Particle Neural Network) [45] a deep neural network for chemical property prediction-to train our charge model, called Affordable Charge Assignments (ACA). ACA is effective at predicting quadrupoles despite being trained only to dipoles, demonstrating the remarkable ability of ML to infer quantities not given in the training dataset. Furthermore, its predictions are extensible to molecules much larger than those used for training. We validate ACA by comparing it to other popular charge models, and find that it is similar to CM5. We then apply ACA to long-time dynamical trajectories of biomolecules, and produce infrared spectra that agree very well with ab initio calculations.</p><p>We briefly review HIP-NN's structure. A more complete description is reported elsewhere in Ref. [45]. HIP-NN takes a molecular conformation as input. The input representation consists of the atomic numbers of all atoms and the pairwise distances between atoms. This representation is simple and ensures that the network predictions satisfy translational, rotational, and reflection invariances. Figure 1 illustrates how HIP-NN processes molecules using a sequence of on-site and interaction layers. On-site layers generate information specific to each local atomic environment and interaction layers allow sharing of information between nearby atomic environments. HIP-NN has previously been successful in modeling energy [45] and pre-existing charge models. [9] Here, we extend the model for dipole prediction using where r i and q i are the position and charge of atom i. HIP-NN's learned charge assignment q i (the ACA charge) is decomposed as a sum over hierarchical corrections,</p><p>As depicted in Fig. 1, each q ì is calculated from the activations (i.e. outputs) of the `-th set of HIP-NN on-site layers. An equivalent decomposition is µ = P `µ`w here µ `= P i q ì r ì is the `-th hierarchical dipole correction. HIP-NN is designed such that higher-order corrections (i.e. µ `for larger `) tend to decay rapidly.</p><p>Training of HIP-NN proceeds by iterative optimization of the neural network model parameters using stochastic gradient descent. The goal of training is to maximize the accuracy of HIP-NN's dipole predictions (as quantified by the root-mean-square-error) subject to regularization. The full ACA model of this paper was generated by an ensemble of four networks.</p><p>More details about HIP-NN and its training process are provided in Ref. [45] and Supporting Information.</p><p>The HIP-NN training and testing data are drawn from the ANI-1x dataset, which includes non-equilibrium conformations of molecules with C, H, N, and O. [46] The ANI-1x dataset was constructed through an active learning procedure [47][48][49] that aims to sample chemical space with maximum diversity. Although ANI-1x was originally designed for potential energy modeling, its chemical diversity also enhances the transferability of ML predictions for other properties, such as the dipole moment. We restrict molecule sizes to 30 atoms or less, and randomly select 396k for training and 44k for testing. Dataset calculations were performed with Gaussian 09 using the !B97x density functional and 6-31G ⇤ basis set. [50] This level of theory will be referred to as the quantum-mechanical (QM) standard throughout this paper. We benchmark the ACA model according to the accuracy of its dipole and quadrupole predictions. To demonstrate extensibility, we test on the DrugBank (⇠ 13k structures) and Tripeptides (2k structures) subsets of the COMP6 benchmark, [46] which contain non-equilibrium conformations of drug molecules and tripeptides. Figure 2 shows the molecular size distribution of these datasets; the molecules in the extensibility sets are roughly four times larger on average than those of ANI-1x, which we used to train ACA.</p><p>Figure 3 shows 2D histograms comparing ACA predicted dipoles and quadrupoles to the QM reference, for all three datasets. We measure the root-mean-square-error (RMSE) and mean-absolute-error (MAE). Left panels of Fig. 3 compare Cartesian dipole components in units of Debye (D). The MAE of 0.078 D for predicting ANI-1x dipoles is comparable to the error between the QM level of theory and experimental dipole measurements. [51] The MAE of ⇡ 0.3 D for predicting DrugBank and Tripeptides dipoles demonstrates the strong extensibility of ACA. Right panels of Fig. 3 compare quadrupole Cartesian components in units of Buckingham (B). The agreement with QM is remarkable (MAE = 0.705 B for the ANI-1x tests) in light of the fact that ACA was trained only to dipoles. Furthermore, ACA continues to make good quadrupole predictions for the much larger COMP6 molecules. We conclude that the ACA charges are physically useful for reproducing electrostatic quantities.</p><p>Additional material quantifying the distributions depicted in Figs. 2 and 3, including error as a function of molecular size, are available in Supporting Information.</p><p>Next, we compare the dipole-inferred ACA model to some conventional charge models. This analysis uses a subset of GDB-11, denoted here as GDB-5, which contains up to 5 heavy atoms of types C, N, and O. [52] The dataset contains a total of 517,133 structures, including non-equilibrium conformations. Four charge models were included in the reference dataset: Hirshfeld, [10] MSK, [5,6] CM5, [8] and population analysis from natural bond orbitals [53] (NBO). Hirshfeld assigns atomic contributions to the electron density based on their relative weighting to the proto-density. MSK charges are constrained to reproduce the dipole moment while attempting to match the electrostatic potential at many points surrounding the molecule. CM5 is an extension of Hirshfeld, empirically parametrized to reproduce ab initio and experimental dipoles. NBO charges are computed as a sum of occupancies from all natural atomic orbitals on each atom. The NBO model is more popular for capturing features such as bond character.</p><p>Figure 4 shows the correlation between each pair of charge models and demonstrates the inconsistency between different approaches to charge partitioning. The strongest correspondence is between CM5 and ACA, with a mean-absolute-deviation of 0.031 e. Other model pairs have mean-absolute-deviations that range from three to eight times larger-a consequence of differing principles used to design these models.</p><p>Conceptually, MSK, CM5, and ACA are similar in that they attempt to partition charge such that the molecular dipole moment is preserved in the point charge representation. We note, however, that MSK differs significantly from CM5 and ACA (Fig. 4). MSK is constrained to match the QM dipole exactly for each given input molecular configuration. This constraint alone is under-determined, and so MSK therefore invokes additional principles for its charge assignment, attempting to fit the far-field electrostatic potential. However, the far-field potential is relatively insensitive to the partial charge assignments of internal atoms. [7][8][9] Because MSK performs its charge assignments according to global (rather than local) criteria, the assigned charges can deviate significantly from the local charge density field.</p><p>Another related difficulty of MSK is that it exhibits a noticeable basis set dependence. [7,9] CM5 was designed to address such drawbacks. [8] Like CM5, our ACA charge model is local-by-design, thus averting the problem of artificial long-range effects. Specifically, ACA seeks a local charge assignment model that best reproduces the QM dipoles over the whole training dataset. We remark that the ACA dipole predictions do not perfectly reproduce the QM dipoles. Allowing for this imperfection may actually be important; collapsing a charge density field into a relatively small number of monopoles while simultaneously forcing the molecular dipole to be exact may be incompatible with locality of the charge model. As we showed in Fig. 4, the CM5 and ACA charges are remarkably consistent, a result we did not anticipate. CM5 reproduces the molecular dipole well, but not as accurately as ACA (See Supporting Information). The reduced accuracy of CM5 dipoles may be due to the fact that it is a is fit to a hybrid of ab initio and experimental data. In contrast, ACA trains to a homogeneous database of QM dipoles. The ML approach has a conceptual advantage: it is fully automated and requires few design decisions (primarily, the specification of an error metric for training). As a consequence, the extension of ACA to new atomic species and to new classes of molecules should be straightforward.</p><p>A strong practical advantage of ACA is that assignment does not require any new QM calculations. We highlight this efficiency advantage by applying ACA to calculate an experimentally-relevant quantity. Inspired by the work of Ref. [26], we use ACA to calculate dynamic dipoles and subsequently infrared spectra for select molecules. Ground-state trajectories were generated from the ANI-1x potential [46] and were 100 ps in length with a 0.1 fs time-step-amounting to a total of 10 6 time-steps. Dipoles were predicted along these trajectories using ACA. Both the molecular dynamics and dipole prediction were performed using only ML, i.e., without any QM calculation. Spectra were made by Fourier transforming the dipole moment autocorrelation function. Harmonic spectra were calculated with the Gaussian 09 software. A comparison of time-domain ML spectra to QM harmonic spectra is shown Figure 5, left panels. Although time-domain and harmonic spectra are not one-toone, the comparison is reasonable since spectral features are harmonic to first order. ACA recovers the harmonic features across all molecules.</p><p>To further validate the ACA dipole predictions, QM calculations were performed at 10 3 subsampled time-steps throughout the trajectories. Fig. 5, right panels, shows that the ACA dipole predictions are in excellent agreement with QM, another validation of ACA's extensibility. The dipole errors are consistent with those observed in the datasets of Fig. 3.</p><p>Note that cholesterol and morphine have 74 and 40 atoms, respectively, whereas our training dataset has no molecules with more than 30 atoms. The quality of the ML-predicted spectra for cholesterol and morphine is similar to those of smaller molecules, such as aspirin.</p><p>We carried out an additional test with smaller molecules of sizes 6 to 15 atoms, using QM to calculate dipoles at all 10 6 time-steps. The resulting infrared spectra are shown in Supporting Information, and are in excellent agreement with our ML-based approach. For these smaller molecules, ACA yields a factor of greater than 10 4 computational speed-up.</p><p>The results are even more dramatic for large molecules.</p><p>In summary, the key contribution of this paper is the formulation of an electrostatically consistent charge model called Affordable Charge Assignments (ACA). We construct the ACA model using a deep neural network that outputs charges. The network is trained to DFT-computed molecular dipole moments over a diverse set of chemical structures. The fast and accurate predictive power of the model was evidenced with extensibility tests (Fig. 3) and infrared spectra (Fig. 5). Although ACA is only trained directly to the molecular dipole, we show that it also captures quadrupole moments, demonstrating transferability.</p><p>ACA is compared with four conventional charge models on a dataset containing over 500k molecules (Fig. 4). The rather poor correlation between most model pairs confirms the ambiguity in charge partitioning. The ACA model correlates well to Charge Model 5 (CM5). CM5 was designed to combine advantages of the Hirshfeld and MSK models. It is parameterized to reproduce a combination of ab initio and experimental dipoles. ACA, like CM5, is a local model that is designed to reproduce dipoles, but unlike CM5, is built entirely from ab initio data. In addition to fast charge assignments, a potential advantage of ACA is its applicability to a wide range of chemically diverse systems, assuming that appropriate training data is available. This work is also a testament to how physics-informed ML can be used to discover properties (here, charge assignment) not employed as an explicit target in the training process. Future work will focus on improving and utilizing ACA for quantum-chemical prediction.</p><p>Improvements to extensible dipole prediction may be made by engaging in dipole-driven active learning. Furthermore, ACA could be trained to higher-order multipole moments such as quadrupoles-this could be important for systems where the dipole does not provide enough of a constraint for charge assignments. Currently ACA is limited to CHNO atoms, but this could be overcome when more diverse datasets are available. Another important drawback of the current model is that charged systems, such as anionic and cationic species, cannot yet be treated. An application using ACA is underway to predict dynamic charges in neutral biomolecular systems to parametrize force fields for molecular dynamics.</p><!><p>The full charge assignments are given by an ensemble prediction using four different random initializations of HIP-NN, each separately trained to the same data. Figure S1 shows the correlation between charge predictions by the members of the ensemble; networks agree to approximately 0.01 e (Fig. S1).The charges produced by the ensemble are not exactly neutral, and so when predicting the charge on a molecule, excess total charge is redistributed evenly across atoms. This redistribution constitutes a very small change, typically 0.001 e or less per atom.</p><p>Dipoles for each datapoint µ are constructed as</p><p>and traceless quadrupoles are constructed as</p><p>where ⌦ is the outer product, and I is the unit dyadic (or Kronecker delta).</p>