Atomic forces for geometry-dependent point multipole and gaussian multipole models

Refinement of Atomic Polarizabilities for a Polarizable Gaussian Multipole Force Field with Simultaneous Considerations of Both Molecular Polarizability Tensors and In-Solution Electrostatic Potentials

Atomic polarizabilities are considered to be fundamental parameters in polarizable molecular mechanical force fields that play pivotal roles in determining model transferability across different electrostatic environments. In an earlier work, the atomic polarizabilities were obtained by fitting them to the B3LYP/aug-cc-pvtz molecular polarizability tensors of mainly small molecules. Taking advantage of the recent PCMRESPPOL method, we refine the atomic polarizabilities for condensed-phase simulations using a polarizable Gaussian Multipole (pGM) force field. Departing from earlier works, in this work, we incorporated polarizability tensors of a large number of dimers and electrostatic potentials (ESPs) in multiple solvents. We calculated 1565 × 4 ESPs of small molecule monomers and dimers of noble gas and small molecules and 4742 × 4 ESPs of small molecule dimers in four solvents (diethyl ether, ε = 4.24, dichloroethane, ε = 10.13, acetone, ε = 20.49, and water, ε = 78.36). For the gas-phase polarizability tensors, we supplemented the molecule set that was used in our earlier work by adding both the 4252 monomer and dimer sets studied by Shaw and co-workers and the 7211 small molecule monomers listed in the QM7b database to a combined total of 13,523 molecular polarizability tensors of monomers and dimers. The QM7b polarizability set was obtained from quantum-machine.org and was calculated at the LR-CCSD/d-aug-cc-pVDZ level of theory. All other polarizability tensors and all ESPs were calculated at the ωB97X-D/aug-cc-pVTZ level of theory. The atomic polarizabilities were developed using all polarizability tensors and the 1565 × 4 ESPs of small molecule monomers and were then assessed by comparing them to the 4742 × 4 ab initio ESPs of small molecule dimers. The predicted dimer ESPs had an average relative root-mean-square error (RRMSE) of 9.30%, which was only slightly larger than the average fitting RRMSE of 9.15% of the monomer ESPs. The transferability of the polarizability set was further evaluated by comparing the ESPs calculated using parameters developed in another dielectric environment for both tetrapeptide and DES monomer data sets. It was observed that the polarizabilities of this work retained or slightly improved the transferability over the one discussed in earlier work even though the number of parameters in the present set is about half of that in the earlier set. Excluding the gas-phase data, for the DES monomer set, the average transfer RRMSEs were 16.25% and 10.83% for pGM-ind and pGM-perm methods, respectively, comparable to the average fitting RRMSEs of 16.03% and 10.54%; for tetrapeptides, the average transfer RRMSEs were 5.62% and 3.95% for pGM-ind and pGM-perm methods, respectively, slightly larger than 5.41% and 3.61% of the fitting RRMSEs. Therefore, we conclude that the pGM methods with updated polarizabilities achieved remarkable transferability from monomer to dimer and from one solvent to another.

PCMRESP: A Method for Polarizable Force Field Parameter Development and Transferability of the Polarizable Gaussian Multipole Models Across Multiple Solvents

The transferability of force field parameters is a crucial aspect of high-quality force fields. Previous investigations have affirmed the transferability of electrostatic parameters derived from polarizable Gaussian multipole models (pGMs) when applied to water oligomer clusters, polypeptides across various conformations, and different sequences. In this study, we introduce PCMRESP, a novel method for electrostatic parametrization in solution, intended for the development of polarizable force fields. We utilized this method to assess the transferability of three models: a fixed charge model and two variants of pGM models. Our analysis involved testing these models on 377 small molecules and 100 tetra-peptides in five representative dielectric environments: gas, diethyl ether, dichloroethane, acetone, and water. Our findings reveal that the inclusion of atomic polarization significantly enhances transferability and the incorporation of permanent atomic dipoles, in the form of covalent bond dipoles, leads to further improvements. Moreover, our tests on dual-solvent strategies demonstrate consistent transferability for all three models, underscoring the robustness of the dual-solvent approach. In contrast, an evaluation of the traditional HF/6-31G* method indicates poor transferability for the pGM-ind and pGM-perm models, suggesting the limitations of this conventional approach.

Assessment of Amino Acid Electrostatic Parametrizations of the Polarizable Gaussian Multipole Model

Accurate parametrization of amino acids is pivotal for the development of reliable force fields for molecular modeling of biomolecules such as proteins. This study aims to assess amino acid electrostatic parametrizations with the polarizable Gaussian Multipole (pGM) model by evaluating the performance of the pGM-perm (with atomic permanent dipoles) and pGM-ind (without atomic permanent dipoles) variants compared to the traditional RESP model. The 100-conf-combterm fitting strategy on tetrapeptides was adopted, in which (1) all peptide bond atoms (-CO-NH-) share identical set of parameters and (2) the total charges of the two terminal N-acetyl (ACE) and N-methylamide (NME) groups were set to neutral. The accuracy and transferability of electrostatic parameters across peptides with varying lengths and real-world examples were examined. The results demonstrate the enhanced performance of the pGM-perm model in accurately representing the electrostatic properties of amino acids. This insight underscores the potential of the pGM-perm model and the 100-conf-combterm strategy for the future development of the pGM force field.

Streamlining and Optimizing Strategies of Electrostatic Parameterization

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.

Transferability of the Electrostatic Parameters of the Polarizable Gaussian Multipole Model

Accuracy and transferability are the two highly desirable properties of molecular mechanical force fields. Compared with the extensively used point-charge additive force fields that apply fixed atom-centered point partial charges to model electrostatic interactions, polarizable force fields are thought to have the advantage of modeling the atomic polarization effects. Previous works have demonstrated the accuracy of the recently developed polarizable Gaussian multipole (pGM) models. In this work, we assessed the transferability of the electrostatic parameters of the pGM models with (pGM-perm) and without (pGM-ind) atomic permanent dipoles in terms of reproducing the electrostatic potentials surrounding molecules/oligomers absent from electrostatic parameterizations. Encouragingly, both the pGM-perm and pGM-ind models show significantly improved transferability than the additive model in the tests (1) from water monomer to water oligomer clusters; (2) across different conformations of amino acid dipeptides and tetrapeptides; (3) from amino acid tetrapeptides to longer polypeptides; and (4) from nucleobase monomers to Watson-Crick base pair dimers and tetramers. Furthermore, we demonstrated that the double-conformation fittings using amino acid tetrapeptides in the αR and β conformations can result in good transferability not only across different tetrapeptide conformations but also from tetrapeptides to polypeptides with lengths ranging from 1 to 20 repetitive residues for both the pGM-ind and pGM-perm models. In addition, the observation that the pGM-ind model has significantly better accuracy and transferability than the point-charge additive model, even though they have an identical number of parameters, strongly suggest the importance of intramolecular polarization effects. In summary, this and previous works together show that the pGM models possess both accuracy and transferability, which are expected to serve as foundations for the development of next-generation polarizable force fields for modeling various polarization-sensitive biological systems and processes.

Accurate Reproduction of Quantum Mechanical Many-Body Interactions in Peptide Main-Chain Hydrogen-Bonding Oligomers by the Polarizable Gaussian Multipole Model

A key advantage of polarizable force fields is their ability to model the atomic polarization effects that play key roles in the atomic many-body interactions. In this work, we assessed the accuracy of the recently developed polarizable Gaussian Multipole (pGM) models in reproducing quantum mechanical (QM) interaction energies, many-body interaction energies, as well as the nonadditive and additive contributions to the many-body interactions for peptide main-chain hydrogen-bonding conformers, using glycine dipeptide oligomers as the model systems. Two types of pGM models were considered, including that with (pGM-perm) and without (pGM-ind) permanent atomic dipoles. The performances of the pGM models were compared with several widely used force fields, including two polarizable (Amoeba13 and ff12pol) and three additive (ff19SB, ff15ipq, and ff03) force fields. Encouragingly, the pGM models outperform all other force fields in terms of reproducing QM interaction energies, many-body interaction energies, as well as the nonadditive and additive contributions to the many-body interactions, as measured by the root-mean-square errors (RMSEs) and mean absolute errors (MAEs). Furthermore, we tested the robustness of the pGM models against polarizability parameterization errors by employing alternative polarizabilities that are either scaled or obtained from other force fields. The results show that the pGM models with alternative polarizabilities exhibit improved accuracy in reproducing QM many-body interaction energies as well as the nonadditive and additive contributions compared with other polarizable force fields, suggesting that the pGM models are robust against the errors in polarizability parameterizations. This work shows that the pGM models are capable of accurately modeling polarization effects and have the potential to serve as templates for developing next-generation polarizable force fields for modeling various biological systems.

Development of Polarizable Gaussian Model for Molecular Mechanical Calculations I: Atomic Polarizability Parameterization To Reproduce ab Initio Anisotropy

A set of atomic polarizability parameters for a new polarizable Gaussian model (pGM) has been developed with the goal to accurately reproduce the polarizability anisotropy, taking advantage of its ability to attenuate all short-range electrostatic interactions, by fitting the ab initio molecular polarizability tensors ( A _pq) calculated at the B3LYP/aug-cc-pVTZ level. For comparison, we also rederived the parameters for three Thole models in which the 1-2 (bonded), 1-3 (separated by two bonds), and 1-4 (separated by three bonds) interactions are fully included. The average percent errors (APEs) of molecular polarizability tensors for 4842 molecules or dimers are 2.98, 3.76, 3.28, and 3.82% for the pGM, Thole linear, Thole exponential, and Thole Amoeba models, respectively, with atom-type independent, universal screening factors (USF). The APEs are reduced further to 2.30, 2.69, 2.25, and 2.48% for the four corresponding polarizable models with atom-type dependent, variable screening factors (VSF). It is encouraging that the pGM with variable screening factors achieved APEs of 1.83 for 1155 amino acid analogs, dipeptides, and tetrapeptides, 1.39 for 28 nucleic acid bases, 0.708 for 1464 water clusters, and 1.99 for 85 dimers of water and biological building blocks. Compared to the new set of models, the APEs of the old Thole models that were fitted to isotropic molecular polarizabilities are 8.7% for set A (without the 1-2 and 1-3 interactions) and 6.3% for set D (with the 1-2 and 1-3 interactions) models, respectively. MPAD, a metric of molecular polarization anisotropy difference based on the diagonal terms of molecular polarizability tensors was defined and applied to assess the polarizable models in reproducing the ab initio molecular polarization anisotropy. The MPADs are 3.71, 4.70, 4.11, and 4.77% for the pGM, Thole linear, Thole exponential, and Thole Amoeba USF models, respectively. The APEs are reduced further to 2.85, 3.58, 2.90, and 3.15% for the four corresponding VSF models. Thus, the new pGM and Thole models notably improve molecular polarization anisotropy. Since pGM attenuates all short-range electrostatic interactions, its application is expected to improve stability in charge fitting, energy, and force calculations and the accuracy of multibody polarization.

Ultrafast dynamics induced by the interaction of molecules with electromagnetic fields: Several quantum, semiclassical, and classical approaches

Several strategies for simulating the ultrafast dynamics of molecules induced by interactions with electromagnetic fields are presented. After a brief overview of the theory of molecule-field interaction, we present several representative examples of quantum, semiclassical, and classical approaches to describe the ultrafast molecular dynamics, including the multiconfiguration time-dependent Hartree method, Bohmian dynamics, local control theory, semiclassical thawed Gaussian approximation, phase averaging, dephasing representation, molecular mechanics with proton transfer, and multipolar force fields. In addition to the general overview, some focus is given to the description of nuclear quantum effects and to the direct dynamics, in which the ab initio energies and forces acting on the nuclei are evaluated on the fly. Several practical applications, performed within the framework of the Swiss National Center of Competence in Research "Molecular Ultrafast Science and Technology," are presented: These include Bohmian dynamics description of the collision of H with H2, local control theory applied to the photoinduced ultrafast intramolecular proton transfer, semiclassical evaluation of vibrationally resolved electronic absorption, emission, photoelectron, and time-resolved stimulated emission spectra, infrared spectroscopy of H-bonding systems, and multipolar force fields applications in the condensed phase.

Torque and atomic forces for Cartesian tensor atomic multipoles with an application to crystal unit cell optimization

New equations for torque and atomic force are derived for use in flexible molecule force fields with atomic multipoles. The expressions are based on Cartesian tensors with arbitrary multipole rank. The standard method for rotating Cartesian tensor multipoles and calculating torque is to first represent the tensor with n indexes and 3(n) redundant components. In this work, new expressions for directly rotating the unique (n + 1)(n + 2)/2 Cartesian tensor multipole components Θpqr are given by introducing Cartesian tensor rotation matrix elements X(R). A polynomial expression and a recursion relation for X(R) are derived. For comparison, the analogous rotation matrix for spherical tensor multipoles are the Wigner functions D(R). The expressions for X(R) are used to derive simple equations for torque and atomic force. The torque and atomic force equations are applied to the geometry optimization of small molecule crystal unit cells. In addition, a discussion of computational efficiency as a function of increasing multipole rank is given for Cartesian tensors. © 2016 Wiley Periodicals, Inc.

Crystal structure prediction of rigid molecules

A non-polarizable force field based on atomic multipoles fit to reproduce experimental crystal properties andab initiogas-phase dimers is described. The Ewald method is used to calculate both long-range electrostatic and 1/r6dispersion energies of crystals. The dispersion energy of a crystal calculated by a cutoff method is shown to converge slowly to the exact Ewald result. A method for constraining space-group symmetry during unit-cell optimization is derived. Results for locally optimizing 4427 unit cells including volume, cell parameters, unit-cell r.m.s.d. and CPU timings are given for both flexible and rigid molecule optimization. An algorithm for randomly generating rigid molecule crystals is described. Using the correct experimentally determined space group, the average and maximum number of random crystals needed to find the correct experimental structure is given for 2440 rigid single component crystals. The force field energy rank of the correct experimental structure is presented for the same set of 2440 rigid single component crystals assuming the correct space group. A complete crystal prediction is performed for two rigid molecules by searching over the 32 most probable space groups.

Limiting assumptions in molecular modeling: electrostatics

Molecular mechanics attempts to represent intermolecular interactions in terms of classical physics. Initial efforts assumed a point charge located at the atom center and coulombic interactions. It is been recognized over multiple decades that simply representing electrostatics with a charge on each atom failed to reproduce the electrostatic potential surrounding a molecule as estimated by quantum mechanics. Molecular orbitals are not spherically symmetrical, an implicit assumption of monopole electrostatics. This perspective reviews recent evidence that requires use of multipole electrostatics and polarizability in molecular modeling.

HPAM: Hirshfeld Partitioned Atomic Multipoles

An implementation of the Hirshfeld (HD) and Hirshfeld-Iterated (HD-I) atomic charge density partitioning schemes is described. Atomic charges and atomic multipoles are calculated from the HD and HD-I atomic charge densities for arbitrary atomic multipole rank l(max) on molecules of arbitrary shape and size. The HD and HD-I atomic charges/multipoles are tested by comparing molecular multipole moments and the electrostatic potential (ESP) surrounding a molecule with their reference ab initio values. In general, the HD-I atomic charges/multipoles are found to better reproduce ab initio electrostatic properties over HD atomic charges/multipoles. A systematic increase in precision for reproducing ab initio electrostatic properties is demonstrated by increasing the atomic multipole rank from l(max) = 0 (atomic charges) to l(max) = 4 (atomic hexadecapoles). Both HD and HD-I atomic multipoles up to rank l(max) are shown to exactly reproduce ab initio molecular multipole moments of rank L for L ≤ l(max). In addition, molecular dipole moments calculated by HD, HD-I, and ChelpG atomic charges only (l(max) = 0) are compared with reference ab initio values. Significant errors in reproducing ab initio molecular dipole moments are found if only HD or HD-I atomic charges used.

A finite field method for calculating molecular polarizability tensors for arbitrary multipole rank

A finite field method for calculating spherical tensor molecular polarizability tensors α(lm;l'm') = ∂Δ(lm)/∂ϕ(l'm')* by numerical derivatives of induced molecular multipole Δ(lm) with respect to gradients of electrostatic potential ϕ(l'm')* is described for arbitrary multipole ranks l and l'. Interconversion formulae for transforming multipole moments and polarizability tensors between spherical and traceless Cartesian tensor conventions are derived. As an example, molecular polarizability tensors up to the hexadecapole-hexadecapole level are calculated for water using the following ab initio methods: Hartree-Fock (HF), Becke three-parameter Lee-Yang-Parr exchange-correlation functional (B3LYP), Møller-Plesset perturbation theory up to second order (MP2), and Coupled Cluster theory with single and double excitations (CCSD). In addition, intermolecular electrostatic and polarization energies calculated by molecular multipoles and polarizability tensors are compared with ab initio reference values calculated by the Reduced Variation Space method for several randomly oriented small molecule dimers separated by a large distance. It is discussed how higher order molecular polarizability tensors can be used as a tool for testing and developing new polarization models for future force fields.