Polarizable charge equilibration model for predicting accurate electrostatic interactions in molecules and solids

For practical simulations of dynamical processes, such as ligands binding to proteins, nucleic acids, and polymers responding to externals fields and stresses, catalysts reacting with substrates, and external fields driving electrochemical reactions, it is necessary to go far beyond the time and length scales of QM through the use of a force field (FF) to describe the structures and forces as they evolve. A critical issue in all such multiscale models is how to accurately describe electrostatic interactions. One common approach is to break the system into fragments, perform QM calculations on each one, and then obtain partial charges from Electrostatic Potential fitting (ESP),¹ Mulliken Population Analysis (MPA),² or other QM charge assignment models.³ Additional discussion on these models can be found in the supplementary material. One disadvantage with these approaches is that the charges are fixed and not allowed to adjust to the changes in the electrostatic environment that occur during dynamics. It can also be burdensome to perform the QM calculations to obtain charges, for example, for the millions of ligands used in Virtual Screening (VS) applications.

The charge equilibration (QEq)⁴ method introduced by Rappé and Goddard in 1991 provides an alternative fast way to predict charges for systems too large for QM. Indeed, carrying out the QEq calculation along an MD trajectory takes into account some changes in polarization during dynamics. Advantages of QEq are that the 3 parameters per atom are derived from atomic ionization energies (valence averaged) and from covalent radii so that they are available for the whole periodic table (through Lr, atomic no. = 103). Also, because the charge on each atom is distributed over a Slater orbital having the size of the atom, QEq can be used to predict charges between bonded atoms to describe the changes during reactions. This made QEq useful for defining a general FF (UFF)⁵ for inorganic-organic systems and for reactive FF (ReaxFF).^6,7 However, it has not been demonstrated whether QEq is as accurate as ESP or MPA in reproducing QM energies nor that the predicted changes in polarization during dynamics agree with QM.

Describing the changes in polarization within a molecule or solid during dynamics or in response to an external electric field is crucial in many applications.^8,9 Consequently, many strategies have been proposed for including polarization into FFs particularly for liquid water and its interactions with ions,^10–12 for modeling of proteins,^13,14 DNA,¹⁵ enzymatic reactions,¹⁶ protein–ligand docking,¹⁷ peptides,¹⁸ and in small-molecule systems.^19,20 Polarization is also important in ion channels and aqueous solution,^11,21 superionic systems,²² piezoelectric and ferroelectrics materials,^23–25 lithium batteries,²⁶ crystal defects and surface energies,^27,28 lattice vibrational frequencies calculations,^28,29 dynamic dielectric response or Raman light scattering,^29,30 hydration energy calculations,³¹ carbon nanotubes,³² and predicting organic crystal structure.³³ This need has led to a number of approaches that have been discussed thoroughly in several reviews.^8,10,34 Our perspective about these methods is summarized in the supplementary material.

In this paper, we propose a new polarizable charge equilibration scheme that builds upon the success of QEq and includes polarization in a generic way that can be easily extended for the entire periodic table. This polarizable charge equilibration (PQEq) model allows charges and polarization to readjust dynamically to attain a constant chemical potential during the simulation. Here the polarization is described by an atomic sized Gaussian shaped electron density cloud that can polarize away from the atomic core in response to internal and external electric fields. The charges on the cores are also described by Gaussian functions and charge can flow from one atom to another based on the QEq scheme. The total electrostatic energy is expressed as a sum of internal atomic energy plus pairwise shielded Coulombic interactions. PQEq uses the same covalent bond radii and atomic ionization energies previously used in QEq. An additional atomic polarization parameter is based on the literature value for atomic polarizability. Thus, the parameters for PQEq are well defined for all elements up to Nobelium (atomic no. = 102).

For validation, we perform a series of high quality QM calculations for 30 structures using cyclohexane and benzene scaffolds containing H, C, N, O, F, Si, P, S, and Cl atoms. The interaction energy was computed by bringing a pair of ±1 point charges (probe dipole) towards each structure along several axes. We show that PQEq produces interaction energies in excellent agreement with QM. In addition, we optimize a parameter set (PQEq1) that increases the accuracy in the interaction energies for these particular compounds. Here, and for the rest of the paper, the total interaction energy between the dipole and target structure is referred to as the interaction energy. For fixed charge models, this interaction is just the electrostatic energy between the dipole and fixed charges (i.e., no polarization). For PQEq, with charge updates and shell polarization turned on, this interaction energy now also reflects the change in energy from polarization.

We then compare the PQEq interaction energies with other common charge models such as Mulliken or ESP charges obtained from QM and those from standard force fields (OPLS³⁵ and AMBER³⁶). We find that the fixed charge methods do not describe the induced polarization in the system.

Based on these results, we believe that PQEq can be used to improve the description of electrostatic interactions for systems in which polarization is important. It can be incorporated into existing force fields such as OPLS, AMBER, CHARMM,³⁷ Dreiding, ReaxFF, and UFF, but it may be necessary to modify some parameters to account for the change in the charge model. This could be useful for a large range of applications including ligand docking to proteins, catalytic reactions, electrocatalysis, ferroelectrics, and growth of ceramics and films.

The PQEq model combines the QEq⁴ model with the shell (Drude oscillator) model.^13,20,24,38 The key difference from previous shell models is that PQEq does not use point charges. Rather the shell electron is described as a Gaussian function having the same size as the core charge. This leads to shielding as the shell electron interacts with its core and with other atoms so that the singularities in point charge descriptions are avoided. The polarization of the shell away from its core in response to the electrostatic field of all other charges and any external field accounts for polarization dynamically. Here we take the mass of the shell to be zero so that it responds adiabatically as the atoms move about in the MD.

For a system of N atoms, each atom, i, is partitioned into two charged sites (core and shell). The core (ρ_ic) consists of two parts: ρ_i with a variable total charge (q_i) and ρ_iZ with a fixed total charge (Z_i). The shell (ρ_is) has a fixed total charge of −Z_i. The shell and core of an atom are connected by an isotropic harmonic spring with force constant K_s (see Figure 1),

where δR is the distance between the core and shell. Equation (1) leads to an atomic or shell polarizability of

where Z is the shell charge and C_unit = 332.0637 is a unit conversion factor that expresses energy in kcal/mol, distance in angstroms (Å), and η as ų. This conversion constant is C_unit = 14.3994 for energies expressed in eV. The η values derived from the atomic polarizability can be computed using high quality ab initio calculations or measured for single atom polarization in response to an external electric field. We use the values tabulated by Miller.³⁹ 

Defining the total charge (core plus shell) on atom i as q_i, the individual charges on the core and shell are

The net atomic charge at the core (q_is + q_ic = q_i) is variable and adjusts to keep the chemical potential constant. There is a positive fixed charge (Z_i) at the core at the position $r→ic$ (i.e., $r→i$⁠) and a negative fixed charge (−Z_i) is at the shell position $r→is$⁠. The displacement of the shell i with respect to its core, $r→is,ic$⁠, is defined as $r→is−r→ic$⁠. The charge density of both the core and the shell is described by a 1s Gaussian charge distribution,

where α_ik, the width of the distribution, is given by

where R_ik is the covalent atomic radius in Å  units and λ is a parameter that converts the overlap of two Gaussian charges to the effective shielding. We determined the value of λ by comparing PQEq and QM electrostatic interaction energies (see below). The PQEq model uses equal atomic and shell radii (i.e., R_i = R_ic = R_is) for each atom i so that the above equations simplify to

This charge distribution has the shape of a spherical 1s Gaussian shape with a width determined by the atomic radius of the atom. The core and shell of each atom have a Columbic interaction with the cores and shells of every other atom in the system. We allow the atomic charges (q_i) to respond to the electrostatic environment based on the QEq scheme.^4,25 The position of the shell is then calculated by balancing the sum of all electrostatic forces on the shell with the spring force (see below). A shell charge of Z_i = 1 is used for all atoms.

The electrostatic energy between two Gaussian charges is given by $Cik,jl(r→)qikqjl$⁠, where

here i and j are the atomic indices, and k and l represent the core (c) and shell (s), respectively. In the case of $r→ik=r→jl$⁠, C_ik,jl is equal to

A well-known problem in other shell and induced dipole^16,40 models is that they suffer from a polarization catastrophe when the shells or dipoles are placed too close together. The Gaussian shielding present in PQEq addresses this issue as the Coulombic interaction energy remains finite, even in the limit of zero interatomic distance.

We describe the PQEq electrostatic energy, $E({r→ic,r→is,qi})$⁠, as a sum of an intra-atomic electrostatic energy $Ei({r→ic,r→is,qi})$ and interatomic pairwise Coulomb interactions,

The internal energy (E_i) is a function of the electronic polarization and total charge on the atom,

The first term on the right hand side of Equation (16), E_i(0,q_i), is the energy required to create a charge, q_i, assuming zero polarization and neutral atomic state as the reference point. We use a Taylor expansion to express this energy term as

and truncate it after second-order terms. The original QEq method⁴ included a radius dependent scaling parameter for hydrogen atom in order to better fit the dipole moment for both alkali hydrides (e.g., LiH) and halogen hydrides (e.g., HF). This leads to a nonharmonic dependence of energy, which in our experience can lead to unstable systems, so PQEq eschews this complication.

We use three conditions to determine the rest of the parameters as follows:

where E_i(0,+1) is the ionization potential (IP), which is the energy required to remove one electron from the atom, and E_i(0,−1) is the electron affinity (EA), which is the energy gained when atom receives one additional electron. Both IP and EA are well known experimentally for nearly all elements. We use the same values as determined by Rappé and Goddard⁴ in which the experimental IP and EA are averaged over the ground state atomic configuration in order to reflect the averaging introduced by bonding to other atoms. Thus, for nitrogen atom, the IP and EA are derived using the averages over the ⁴S, ²D, and ²P states associated with the (2s)²(2p)³ configuration and the ³P, ¹D, and ¹S states associated with the (2s)²(2p)² and (2s)²(2p)⁴ configurations of the ions. Solving Equations (18)–(20) for the unknowns yields

where $χi0$ is the Mulliken electronegativity⁴¹ of atom i and $Jii0$ is the idempotential (hardness) or electron capacity of atom i, which resists the electron flow to or from an atom. We replace the second term on the right hand side of Equation (16) with the Coulombic interaction between the core and shell of atom i plus a spring interaction between the core and shell of atom i. Therefore, Equation (15) can be written as

Ignoring the $O(qi3)$ term in Equation (23), the electrostatic energy of the system is given by

where the second sum is the pairwise shielded Coulomb interaction energy between all cores and shells, which can be expanded to give the total electrostatic energy as

A serious problem in most classical MD/MM applications is that fixed charges are assigned to each atom. Such fixed charges (even the most reliable ones from ab initio calculations) do not respond to the changes in electrostatic environment, which decreases the accuracy. This problem becomes paramount for reactive force fields (e.g., ReaxFF⁷) where the bond connectivities of atoms change during reactive MD simulations, requiring updates of the atomic charges (ideally at each time step). As in QEq, the PQEq model allows the charge distribution on the various atoms to change as the electrostatic environment changes during the dynamics. The optimum charge distribution is computed from the conditions that the chemical potentials (∂E/∂q_i) are equal for all of the atoms (which provides N − 1 conditions where N is the number of atoms) and that the total charge is conserved

where Q is the total charge of the system. We use Lagrange multipliers to guarantee this constraint as the charges are optimized. The energy expressions with the Lagrange multiplier, μ, is

Setting the derivative of Equation (27) equal to zero yields

where H_ij is an N by N matrix and δ_ij is the Kronecker delta function. The diagonal elements of H_ij matrix (δ_ij = 1) denote the idempotential of the atoms while the off-diagonal elements represent the Coulombic interactions between the variable charge part of the cores (i.e., q_i). A_i in Equation (30) is a vector of length N. The first term of A_i is the electronegativity of the atom. The second term is the sum of Coulombic interaction coefficient between the variable charge part of core i (i.e., q_i) and the fixed charge component of all other cores and shells (i.e., Z_j and −Z_j). Note that in Equation (30), A_i is a fixed quantity for each atom during the charge minimization, which reduces to $Ai=χi0$ if polarization is not included, as in the QEq model.^4,42 In Equations (28)–(30), the Lagrange multiplier μ is the chemical potential that constrains the sum of the atomic charges to be equal to the total charge of Q. Solving Equations (28)–(30) leads to

Applying Equation (26) we get

which is solved to obtain μ,

where $q∼i$ and $q^i$ are the fictitious charges. In practice, we solve Equation (33) by partitioning it into two sub equations. For a neutral system, Equation (33) can be divided into two subequations as

Finally, the instantaneous total charges on each atom (q_i = q_ic + q_is) can be written as

The above formulation for PQEq omits the presence of external electric fields, which is included in the supplementary material. A frequency-dependent response can be obtained from time-dependent fields.

Exactly solving for the charges that satisfy the QEq condition involves inverting an N by N matrix, which scales as O(N³). Since this is required every time step,⁴ this process is computationally too expensive to be practical for large systems. A practical solution to this problem is the PCG method implemented in the PuReMD^43,44 and LAMMPS⁴⁵ software packages.

We use PCG to solve Equations (34) and (35). The efficiency and convergence of this iterative conjugate-gradient (CG) method depend on the spectrum of the coefficient matrix. The PCG method uses a second matrix (preconditioner) to transform the coefficient matrix to obtain improved spectral properties. This preconditioner involves an incomplete factorization of the coefficient matrix. In particular, incomplete LU (ILU) factorization (where L and U are lower and upper triangular) can be used for solving these sparse linear systems.⁴⁶ For QEq, Aktulga et al. studied the performance, stability, and accuracy of the ILU-based preconditioners for various model systems.⁴⁴ They showed that ILU-based preconditioners dramatically reduce the number of iterations while allowing the same L and U factors to be used effectively as preconditioners over several steps, due to the slow changes in the simulation environment. We extended this ILU-based preconditioner method to PQEq and coupled it with shell relaxation (see Sec. II E) to calculate the PQEq charges while updating the shell position.

In our formulation of the PQEq dynamics, we choose to displace the core of each atom together with its shells as a rigid body during every time step of the dynamics. After moving core plus shell, the first step of the next iteration is to calculate the electrostatic field on every particle and to solve the PQEq equations (using the PCG method) for the new charges. Next, we fix the core positions and update the positions of all shells simultaneously using a one-step relaxation as follows. If necessary, this process of updating atomic charges and shell relaxation with fixed cores can be repeated for several iterations to attain self-consistency for troublesome geometries. However, we find that one cycle is normally sufficient to reach equilibrium for each time step after the first.

The shell position for each atom is obtained by balancing the effect of the electrostatic field due to all external atoms with intra-atomic interactions involving only the core and its shell. These forces are calculated by taking the derivative of the electrostatic energy (Equation (25)) with respect to the shell position,

We solve Equations (37) and (38) to find the optimal position of shells (r_is) using a single iteration of the Newton-Raphson method. Since, the shell is typically very close to its core (usually < 0.1 Å), we neglect the effect of the external core and shell charges on this second derivative to avoid inverting the Hessian. We assume here that the shell is massless, so that it relaxes instantaneously to its zero-force position, with no inertial delay. Therefore, we estimate the new position of the shell by assuming that F_intra and F_external are collinear. Thus, the new position of the shell is computed by

Although, this problem is not strictly one dimensional, we use the above second derivative allowing the shell to rotate in the direction in which the external field acts.

In order to validate the accuracy of PQEq and to perform optimization of the model (if needed), we must decide the criteria to use for comparison and optimization. The normal practice in most FFs is to use QM charges. We discuss in the supplementary material that QM charges are not reliable for this purpose. Instead, we use the QM interaction energy. We probe each of the 30 molecules in our validation set with a pair of ±1 point charges separated by 1 Å  to describe the interaction with dipole and higher order multipoles. For convenience, we refer to this pair of point charges as an electric dipole. The interaction energy from QM is shown as a function of distance and used as the reference energy.

We selected scan axes along a variety of symmetry directions to provide insight into about how the polarization depends on the elements. Care was taken to avoid close contacts with the nearby atoms. Figure 2 shows an example for a cyclohexane molecule of electric dipole scans along several different directions. These scans are performed towards a backbone atom (d1 and d2), along a bond (d3 and d4), perpendicular to a bond (d5), and toward the center of mass of the structure (d6). For all calculations, we use the standard B3LYP hybrid flavor of density functional theory (DFT), including both the generalized gradient approximation and a component of the exact Hartree–Fock (HF) exchange.^47,48 These calculations were performed with the 6-311G(d,p) (or 6-311G**++) basis set.⁴⁹ 

We used a set of 30 molecular structures in our validation of the PQEq model. We designed this data set to cover H, C, N, O, F, Si, P, S, and Cl elements in a balanced manner. These structures are depicted in Figure S3 of the supplementary material. We use cyclohexane and benzene rings as the framework for these molecular structures, replacing C and H with the above atoms. This framework provides a reasonable number of atoms and bond types for studying charge transfer and polarization effects. The molecular structures are at their equilibrium geometries optimized using QM with same DFT method and basis set described above. Then, the electric dipole is scanned along various directions with respect to these molecular structures. The scan directions were selected after extensive preliminary calculations to probe properly the amount of polarization and electrostatic potential change during the scan. We excluded the cases that resulted in less than 2 kcal/mol change in the energy throughout the scan. We also avoided scanning directions that could lead to very close interaction of the dipole with nearby atoms. In addition, to avoid non-electrostatic interactions arising from Pauli principle repulsion at close distances, we scan only up to the inflection point (attractive forces) of the electrostatic potential curve. We find this distance to be near 2.5 Å  for most of the cases so that the electric dipole is scanned from 10.0 Å  up to 2.5 Å  with respect to the reference point for all cases.

The above considerations resulted in a total of 68 scans for the above molecular structures. The change of QM electric dipole energy with the distance for each case is shown in Figure S4 of the supplementary material and for several selected cases in Figure 3.

In this paper, we present two sets of PQEq parameters. The first set (denoted as PQEq) uses the same χ, J, and R_cs parameters as in the QEq method⁴ which were obtained from standard bond radii and experimental ionization energies. These are available for all the elements of the periodic table up to Lawrencium (Lr) (atomic no. = 103). Using Equation (3), we derived the K_s values based on experimental or high quality ab initio calculations of atomic polarizabilities in the presence of an external electric field. These values are available up to Nobelium (No) (atomic no. = 102).³⁹ The exception is for H atom where we use IP = 11.02 eV and EA = 1.96 eV to define χ and J as did Rappé and Goddard⁴ and we take the K_s value for H from the Karasawa and Goddard calculation for the polarizability of Polyvinylidene-fluoride crystal that they fitted to a shell-model.²⁴ Our results in Sec. V B show that this default PQEq parameter set, with no additional optimization predicts interaction energies from QM very well.

For the second series of parameters (PQEq1), we performed a constrained optimization of the atomic χ and J parameters with respect to the QM derived energy for all 68 cases. We used CG optimization to minimize the difference between the energies computed by PQEq1 and QM. The total error is defined as the weighted mean square error (MSE),

where ω_i is the weight and E_PQEq and E_QM are the energies computed by PQEq1 and QM, respectively. Constraints were applied to ensure that the parameters obey the general trends of the periodic table. That is, we require that within a row of the periodic table, the atoms should become more electronegative as we move to the right. Similarly, we require that atoms become more electronegative as we move down a column in the periodic table. This was enforced at each step of the optimization. These constraints defend against overfitting and ensure better transferability of the final parameter set. We find the total error in Equation (41) to decrease from 1219.29 to 471.55 during the PQEq1 optimization. Thus for this class of systems, the PQEq1 parameters should provide a more accurate description of the electrostatics. The changes in parameter values are generally small. The maximum change in the PQEq1 parameters is for fluorine (F) atom, which led to a 19.96 percent change in the χ parameter. The comparison between the parameters before and after optimization is shown in Table S3 in the supplementary material. The PQEq and PQEq1 parameters are tabulated in Tables S1 and S2 in the supplementary material, respectively. We also provide the electronic versions of these files.

Figure 3 shows the comparison (one for each atom type) between the interaction energies computed by QM, PQEq, and PQEq1 for the scan of the electric dipole at different distances. Here, the interaction energy includes the polarization effect during the scan. The dipole scan directions are shown with the dotted lines on the molecular structure schematics for each case. The comparisons for the rest of the cases are shown in Figure S4 in the supplementary material.

Based on these results, we choose the effective shielding parameter in Equation (10) to be λ = 0.4628. This value is close to the corresponding number in QEq model (0.4913) using Slater-type orbitals.⁴ The results show good agreement of PQEq with QM. This suggests that the PQEq general parameter set can describe accurately the electrostatic potential for a variety of molecular structures and environments. Thus, we expect good transferability of the PQEq model to new materials. This often has been a challenge for previous FFs and charge calculation models. As expected, the results from the PQEq1 parameter set show better agreement with QM and may be useful for other systems that contain similar structures as in our database. In particular, PQEq1 provides a dramatic improvement for molecules containing Fluorine element (Figure 3(e)).

The energy comparison is the crucial criterion to test the accuracy of the PQEq model, but we are also concerned to determine if the computed charges are consistent with chemical intuition. This is particularly important for using PQEq partial atomic charges in the electrostatic potential term of FFs that have been developed with different charge models. We compute the partial atomic charges for all of the molecular structures using PQEq and PQEq1 parameter sets and compare them with ESP and MPA charges.

For the MPA and ESP charge calculations, we use several flavors of DFT including B3LYP,⁴⁸ M06,⁵⁰ and Perdew-Burke-Ernzerhof (PBE)⁵¹ with several Gaussian basis sets including 6-311G, polarizable 6-311G**, polarizable and diffusive 6-311G**++.^49,52 The results for two selected cases are shown in Figure 4 and for the other structures in Figure S5 in the supplementary material. We find for all cases that the PQEq and PQEq1 charges are in the range of ESP and MPA charges. It is well known that ESP and MPA charges sometime lead to unintuitive charge assignments (see Section 2 of the supplementary material), but we have not found such cases with PQEq and PQEq1.

To test the stability of the PQEq model for MD/MM simulations, we examined the reactive MD simulations of the hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX)⁵³ crystal at high temperatures using ReaxFF-lg⁵⁴ reactive force field. These calculations use the LAMMPS⁴⁵ MD simulation package with our implementation of the PQEq methodology. First, we minimized the total energy of the crystal (168 atoms) using the CG method. Then, we equilibrated this structure using the (NVT) ensemble at 50 K for 2 ps. Then, we carried out MD-NVT simulations using a heating rate of 0.7 K/fs, during which the temperature increased from 50 to 3500 K. Finally, the structure was maintained at 3500 K for ∼50 ps using MD-NVT simulations. See Section 9 of the supplementary material for more details of the simulations. Under these conditions, bonds are broken with the fragments interacting to form new bonds. We consider this as a good test case for PQEq. Indeed, we find that PQEq provides a stable description of the complex evolution of the dynamics as bonds break and rearrange. There are smooth changes of the temperature (T), potential energy (E_p), and electrostatic energy (E_PQEq) of the RDX crystal during the simulation (Figure S8 in the supplementary material). We note that at 3500 K, both E_p and E_PQEq decrease for several ps due to fast chemical reactions at this high temperature and then reach equilibrium. We find that the atomic charges and shell positions fluctuate in response to the changes in the electrostatic environment as they were updated every time step. The changes with time of the charges and shell positions are shown in Figure 5. The shell positions remain stable with respect to the core with up to 0.05 Å  displacement from the core. This shows that the K_s values derived from the literature atomic polarizabilities are useful for simulation of these elements at high temperatures. The value of K_s should be tested for other elements prior to dynamics, particularly at high temperatures.

In addition, we performed a series of MD simulations to demonstrate the stability of PQEq during dynamics. For this purpose, we utilized the MD-NVT simulations above to heat the system from 50 to 3500 K and maintained the temperature at 3500 K for 5 ps. After this step, we performed a MD-NVE simulation for 20 ps. We observed reasonable dynamics throughout the simulation. One point that requires attention is the small drift in energy during the MD-NVE simulation, a known problem with the ReaxFF force field used here. To determine if shell polarization contributed to the change in energy, we repeated the simulation without shell polarization (QEq) and found a very similar change in energy, as shown in Figure S9 in the supplementary material. We conclude that the shell model does not introduce any new errors or instabilities in the conservation of the total energy of the system.

In this section, we compare the dipole interaction energies from QM, PQEq, and PQEq1 with the results from ESP, MPA, OPLS, AMBER, PQEq0, QEq, and QEq0. Here

• PQEq0 refers to PQEq with the charges fixed prior to the introduction of the dipole. Here, some part of the polarization energy is included via the shell polarization.

• QEq keeps the shell fixed to the core and equilibrates the charge as the dipole is scanned. Here, the charge updates capture part of the polarization energy.

• QEq0 keeps the shell fixed to the core and the charges fixed prior to the introduction of the dipole. In this case, no polarization is included.

For ESP, MPA, OPLS, and AMBER, PQEq0, and QEq0, we first compute the charges for each molecular structure in the absence of the electric dipole and then fix the charges to calculate the interaction energy at different distances of the electric dipole from the molecule.

The OPLS and AMBER FFs are often used for simulations of large organic and protein systems. These charges are fixed and assigned based on the type of the atoms and its bonding type. AMBER and CHARMM have standard charges for standard amino acids and nucleic acid bases, but for other molecules the charges are assigned from QM using MPA or ESP. Thus we also include the ESP and MPA charges computed using the B3LYP flavor of DFT and 6-311G** basis set. The results for six selected cases are shown in Figure 6 while the remaining cases are in Figure S6 in the supplementary material.

The importance of polarization is clearly shown in Figures 6(a)–6(c) where the scans are performed towards the backbone C atom in cyclohexane (Figure 6(a)), toward the H atom and perpendicular to the benzene ring (Figure 6(b)), and toward the C–Si bond middle point in a cyclohexane-based molecule (Figure 6(c)). Here only PQEq, PQEq1, and QEq predict interaction energies in a good agreement with QM. Fixed charge methods sometimes fail to predict the correct sign of the interaction energy as shown in Figures 6(a) and 6(c). For the remaining cases in Figure 6, the scans are performed towards the N (Figure 6(d)) and O (Figure 6(e)) atoms in cyclohexane-based molecules and towards O (Figure 6(f)) atom in the plane of nitrobenzene molecule. For these polar systems involving N and O atoms, the fixed charge models account for some of the polarization that occurs along the bonds of polar to nonpolar atoms. We see here that PQEq1 does an excellent job of fitting QM, whereas PQEq is accurate for N but overestimates the polarization for O. This suggests that the reference polarizability for O may be too large.

We note here that QEq0 leads to an accuracy similar to the ESP or MPA obtained from QM. Thus for assigning fixed charges for use in MD, there is no longer a need to do QM, which can save considerable expense for applications such as virtual screening over millions of molecules or simulations on very large molecules.

In existing software codes, such as NAMD,⁵⁵ CHARMM,⁵⁶ and DESMOND,⁵⁷ major changes would be needed to recalculate the charges along the MD trajectory. However, including just the shell polarization would be fairly simple to add to current software packages. This would allow the accuracy of PQEq0, which captures 21.8%, 55.6%, and 62.6% of the total polarization energies in Figures 6(a)–6(c), respectively. Therefore, PQEq0 could provide dramatically improved descriptions of the polarization in very large systems (the shell polarization requires only a one step update in shell position each iteration). However, some reoptimization of the force field parameters might be needed when PQEq methodology is used to replace the charge model in other force fields.

PQEq and PQEq0 should be particularly interesting for MD simulations of highly polarizable systems such as ferroelectrics and electrochemical systems with solvents and applied fields.

We show that the PQEq polarizable charge equilibration method provides accurate descriptions of the electrostatic interactions for MD simulations. This PQEq model uses atomic sized Gaussian shaped core and shell densities connected with an isotropic harmonic spring. The atomic parameters of PQEq are obtained from standard atomic ionization energies, standard covalent radii, and literature atomic polarizabilities, which we provide here up to Nobelium (atomic no. = 102). Thus, no parameters have been optimized.

We validated the accuracy of PQEq by comparing the electrostatic polarization energies as an electric dipole is brought up to the molecule for 30 molecules (68 cases) involving H, C, N, O, F, Si, P, S, and Cl atoms. We find that PQEq is in good agreement with QM. We also considered the PQEq1 model in which the atomic parameters (χ and J) are optimized against QM polarization energy. This led to improvements especially for fluorine element.

We also presented the results for various fixed charge models: ESP, MPA, AMBER, OPLS, QEq0, and PQEq0. These methods are generally similar and much less accurate than the polarized models. However, we see that PQEq0 is capable of capturing significant parts of the polarization with just adjustments of the shell polarization while keeping the charges fixed. Thus, PQEq0 can offer significantly improved accuracy compared to other fixed charge models. We expect that PQEq and PQEq0 will be useful for many applications including ligand docking to proteins, catalytic reactions, electrocatalysis, ferroelectrics, fuel cells, lithium ion batteries, and growth of ceramics and films.

See supplementary material for the following:

• PQEq database of test molecules, full set of PQEq parameters, PQEq1 optimized parameters, polarization energy comparison, charge comparison, and reactive simulation of RDX;

• PQEq and PQEq1 electronic parameter files.

This work was supported by the Joint Center for Artificial Photosynthesis, a DOE Energy Innovation Hub, supported through the Office of Science of the U.S. DOE under Award No. DE-SC0004993. We thank Dr. Qingsong Zhang, Dr. Sergey Zybin, and Dr. Andres Jaramillo-Botero for useful discussions.