Streamlining and Optimizing Strategies of Electrostatic Parameterization

Abstract

Accurate characterization of electrostatic interactions is crucial in molecular simulation. Various methods and programs have been developed to obtain electrostatic parameters for additive or polarizable models to replicate electrostatic properties obtained from experimental measurements or theoretical calculations. Electrostatic potentials (ESPs), a set of physically well-defined observables from quantum mechanical (QM) calculations, are well suited for optimization efforts due to the ease of collecting a large amount of conformation-dependent data. However, a reliable set of QM ESP computed at an appropriate level of theory and atomic basis set is necessary. In addition, despite the recent development of the PyRESP program for electrostatic parameterizations of induced dipole-polarizable models, the time-consuming and error-prone input file preparation process has limited the widespread use of these protocols. This work aims to comprehensively evaluate the quality of QM ESPs derived by eight methods, including wave function methods such as Hartree-Fock (HF), second-order Møller-Plesset (MP2), and coupled cluster-singles and doubles (CCSD), as well as five hybrid density functional theory (DFT) methods, used in conjunction with 13 different basis sets. The highest theory levels CCSD/aug-cc-pV5Z (a5z) and MP2/aug-cc-pV5Z (a5z) were selected as benchmark data over two homemade data sets. The results show that the hybrid DFT method, ωB97X-D, combined with the aug-cc-pVTZ (a3z) basis set, performs well in reproducing ESPs while taking both accuracy and efficiency into consideration. Moreover, a flexible and user-friendly program called PyRESP_GEN was developed to streamline input file preparation. The restraining strengths, along with strategies for polarizable Gaussian multipole (pGM) model parameterizations, were also optimized. These findings and the program presented in this work facilitate the development and application of induced dipole-polarizable models, such as pGM models, for molecular simulations of both chemical and biological significance.

Figure 1.

Molecules utilized in this work. The sets of SMALL and LARGE molecules are shown in the upper and lower panels, respectively. For each molecule, the chemical structure together with its corresponding chemical formula is shown. The charged molecule is shaded in red; polar and nonpolar species are in blue and gray, respectively.

Figure 2.

Relative root-mean-square (RRMS) error of the electrostatic potentials (ESPs) between the given and reference methods. The left panel is the result from the SMALL set with the reference method CCSD/aug-cc-pV5Z (a5z), while the right panel refers to the LARGE set with the reference method MP2/aug-cc-pV5Z (a5z). Outliers are marked with a star sign, and median values are shown in orange. The outliers are indicated by “*”, which indicates when the values are outside the range of (Q1−1.5 × IQR) to (Q3 + 1.5 × IQR), where Q1 and Q3 represent the 25 and 75% quartiles, respectively, and IQR denotes the interquartile range between Q3 and Q1. Abbreviations of the basis sets can be found in Section 2.2.

Figure 3.

Relative root-mean-square (RRMS) error of ESPs derived from different QM methods with respect to different reference methods. Here, reference methods were CCSD/aug-cc-pV5Z and MP2/aug-cc-pV5Z for the SMALL (left) and LARGE (right) sets, respectively. The RRMS is shown here on a logarithmic scale; a linear scale plot can be found in Figure S1.

Figure 4.

Relative RMS per molecule (〈RRMS〉) as a function of the restraining strength at (a) the permanent point charge of the pGM-ind model (a_q in eq 4) and (b) the permanent dipole of the pGM-perm model (a_p in eq 7). Reference method and basis sets, namely, CCSD/aug-cc-pV5Z (a5z) and MP2/aug-cc-pV5Z (a5z), are utilized for SMALL (circle point) and LARGE (square point) sets, respectively.

Figure 5.

Comparison of the relative root-mean-square (RRMS) errors for multipole fitting to QM ESPs for three models: RESP, pGM-ind, and pGM-perm. For small set (a), QM ESPs are calculated at the CCSD/aug-cc-pV5Z (a5z) level, while for LARGE set (b), QM ESPs are calculated at the MP2/aug-cc-pV5Z (a5z) level. The star signs indicate the outliers, and their chemical structures are shown by the insets.