Minimally-corrected partial atomic charges for non-covalent electrostatic interactions

TL;DR

This paper introduces a new method for refining molecular partial atomic charges by making minimal corrections to existing reference PACs, resulting in improved electrostatic potential reproduction while maintaining symmetry invariance.

Contribution

The authors develop a minimally-corrected PAC scheme that enhances ESP accuracy and preserves molecular symmetry, applicable to various reference charge sets.

Findings

MC-PACs significantly improve ESP reproduction.

MC-PACs preserve reference PACs' symmetry invariance.

Comparable ESP quality to ChElPG method when using iterative-Hirshfeld PACs.

Abstract

We develop a new scheme for determining molecular partial atomic charges (PACs) with external electrostatic potential (ESP) closely mimicking that of the molecule. The PACs are the "minimal corrections" to a reference-set of PACs necessary for reproducing exactly the tensor components of the Cartesian zero- first- and second- molecular electrostatic multipoles. We evaluate the quality of ESP reproduction when "minimally correcting" (MC) Mulliken, Hirshfeld or iterated-Hirshfeld reference PACs. In all these cases the MC-PACs significantly improve the ESP while preserving the reference PACs' invariance under the molecular symmetry operations. When iterative-Hirshfeld PACs are used as reference the MC-PACs yield ESPs of comparable quality to those of the ChElPG charge fitting method.

Figures

Tables

Table 1. Table 1: The PAC change ‖ Δ ​ q ‖ ∞ = m ​ a ​ x 1 ≤ a ≤ A ​ | Δ ​ q a | subscript norm Δ 𝑞 𝑚 𝑎 subscript 𝑥 1 𝑎 𝐴 Δ subscript 𝑞 𝑎 \left\|\Delta q\right\|_{\infty}=max_{1\leq a\leq A}\left|\Delta q_{a}\right| induced by mcD and mcDQ for 0PA (where the reference PACs are all zero), MPA, HPA and iHPA for the set of molecules of Fig.  3 . Also shown the number of atoms A 𝐴 A the number of degrees of freedom F 𝐹 F and the number of constraints C 𝐶 C for each molecule.

Molecule	Sym	$A$	$F$	$C$	${\Vert \Delta q\Vert }_{\infty }$ (mcD)								${\Vert \Delta q\Vert }_{\infty }$ (mcDQ)
					0PA		MPA		HPA		iHPA		0PA		MPA		HPA		iHPA
Pyridazine	$C_{2v}$	10	7	6	0.	11	0.	02	0.	06	0.	04	0.	27	0.	29	0.	17	0.	12
Ethylamine	$C_s$	10	7	6	0.	04	0.	01	0.	03	0.	01	1.	17	1.	61	1.	22	1.	28
Acetoacetic acid	$C_1$	10	9	6	0.	07	0.	02	0.	03	0.	01	0.	54	0.	26	0.	25	0.	21
Acetone	$C_{2v}$	10	5	5	0.	17	0.	02	0.	05	0.	01	0.	52	0.	48	0.	34	0.	08
Ethylene-glycol	$C_1$	10	10	9	0.	09	0.	03	0.	05	0.	02	0.	67	0.	88	0.	64	0.	54
Cyclopentadienone	$C_{2v}$	10	6	6	0.	09	0.	04	0.	03	0.	00	0.	19	0.	22	0.	05	0.	02
Oxetane	$C_s$	10	7	6	0.	08	0.	05	0.	03	0.	01	1.	09	1.	29	1.	09	1.	09
3-Imino-2,3-dihydroisoxzole	$C_s$	10	8	6	0.	05	0.	02	0.	02	0.	01	0.	79	0.	83	0.	69	0.	41
Lithium-dimethylamine	$C_{2v}$	10	6	4	0.	31	0.	16	0.	05	0.	04	0.	31	0.	16	0.	05	0.	04
Methylacetate	$C_s$	11	9	6	0.	10	0.	02	0.	02	0.	01	0.	53	0.	22	0.	26	0.	01
Pyridine	$C_{2v}$	11	7	4	0.	06	0.	02	0.	03	0.	01	0.	15	0.	23	0.	09	0.	03
2-Cyanopyridine	$C_s$	12	12	6	0.	12	0.	05	0.	04	0.	01	0.	19	0.	19	0.	07	0.	03
1-fluoro-4-nitrobenzene	$C_{2v}$	14	9	6*	0.	05	0.	06	0.	01	0.	01	0.	13	0.	20	0.	05	0.	04
Morpholine	$C_s$	15	9	6	0.	03	0.	01	0.	02	0.	01	0.	47	0.	04	0.	33	0.	08
Quinoline	$C_s$	17	17	6	0.	03	0.	12	0.	02	0.	01	0.	10	0.	26	0.	05	0.	02
Fluorocyclohexane (A)	$C_s$	18	12	6	0.	05	0.	03	0.	02	0.	01	0.	17	0.	06	0.	05	0.	02
Fluorocyclohexane (E)	$C_s$	18	12	6	0.	04	0.	03	0.	01	0.	01	0.	17	0.	08	0.	07	0.	04
Median					0.	07	0.	03	0.	03	0.	01	0.	31	0.	23	0.	17	0.	04
Max					0.	31	0.	16	0.	06	0.	04	1.	17	1.	61	1.	22	1.	28

Table 2. Table 2: The median/maximal MARD (for the set of molecules used above) determined for each PAC reference found for: non-corrected, and minimally-corrected schemes, mcD and mcDQ .

	0PA	MPA	HPA	iHPA
no-correction	NA	0.39/1.90	0.38/0.78	0.14/0.32
mcD	0.41/0.67	0.18/0.60	0.18/0.53	0.08/0.17
mcDQ	0.12/0.26	0.10/0.28	0.11/0.19	0.07/0.15

Equations

$$4\pi {ϵ}_0\phi \left(\mathbf{r}\right)$$

$$=\frac{Q}{r}+\frac{{\mu }_i{r}_i}{r^3}+\frac{1}{2}\frac{r_i{\Theta }_{ij}{r}_j}{r^5}+....$$

$$W_{es}=Q{\phi }^{other}+{\mu }_i{\phi }_i^{other}+\frac{1}{6}{\Theta }_{ij}{\phi }_{ij}^{other}+\dots$$

$${\phi }^{PAC}\left(\mathbf{r}\right)=\frac{e}{4\pi {ϵ}_0}\sum _a\frac{q_a}{\left|\mathbf{r}-{\mathbf{R}}_a\right|},$$

$$MARD\left({\phi }^{PAC},\phi \right)=⟨\left|\frac{{\phi }^{PAC}\left(\mathbf{r}\right)-\phi \left(\mathbf{r}\right)}{\phi \left(\mathbf{r}\right)}\right|⟩,$$

$$c$$

$$c_i$$

$$c_{ij}$$

$$L_{mcDQ}$$

$$-\sum _{xy,yz,zx}{\lambda }_{ij}^{OD}{c}_{ij}-\sum _{x,y,z}{\lambda }_{ii}^D{c}_{ii}$$

$$S$$

$$S=\left(\begin{array}[]{ccccccccccccc}\text{\textSFi}&\cdots&\text{\textSFiii}&-1&\text{\textSFi}&&\text{\textSFiii}&\text{\textSFi}&&\text{\textSFiii}&\text{\textSFi}&&\text{\textSFiii}\\ \vdots&I_{A\times A}&\vdots&\vdots&\vdots&-D_{A\times 3}&\vdots&\vdots&-T_{A\times 3}^{D}&\vdots&\vdots&-T_{A\times 3}^{OD}&\vdots\\ \text{\textSFii}&\cdots&\text{\textSFiv}&-1&\text{\textSFii}&&\text{\textSFiv}&\text{\textSFii}&&\text{\textSFiv}&\text{\textSFii}&&\text{\textSFiv}\\ 1&\cdots&1&0&0&0&0&0&0&0&0&0&0\\ \text{\textSFi}&\cdots&\text{\textSFiii}&0&0&0&0&0&0&0&0&0&0\\ &D_{3\times A}&&0&0&0&0&0&0&0&0&0&0\\ \text{\textSFii}&\cdots&\text{\textSFiv}&0&0&0&0&0&0&0&0&0&0\\ \text{\textSFi}&\cdots&\text{\textSFiii}&0&0&0&0&0&0&0&0&0&0\\ &T_{3\times A}^{D}&&0&0&0&0&0&0&0&0&0&0\\ \text{\textSFii}&\cdots&\text{\textSFiv}&0&0&0&0&0&0&0&0&0&0\\ \text{\textSFi}&\cdots&\text{\textSFiii}&0&0&0&0&0&0&0&0&0&0\\ &T_{3\times A}^{OD}&&0&0&0&0&0&0&0&0&0&0\\ \text{\textSFii}&\cdots&\text{\textSFiv}&0&0&0&0&0&0&0&0&0&0\end{array}\right)$$

Full text

\SectionNumbersOn

1 Introduction

Partial atomic charges (PACs), i.e. point charges placed on the nuclei position of a molecule are often used in large-scale molecular mechanics calculations to replace the detailed quantum mechanical charge distributions. 1, 2, 3, 4, 5, 6, 7 The model is extremely useful since by using them the long-range electrostatic forces acting *between *molecules can be expressed as a sum of pairwise interactions, enabling a fast computation, important especially as molecules jiggle around and rotate quite a lot during the course of the simulation. The question of just how to determine PACs for this purpose is critical. We argue that the most important constraint is the exact reproduction of the low-order electrostatic moments (ESM), the monopole $Q=e\int\rho\left(\boldsymbol{r}\right)d\boldsymbol{r}$, which is the total charge of the system, the dipole $\mu_{i}=e\int\rho\left(\boldsymbol{r}\right)r_{i}d\boldsymbol{r}$ ($i=x,y,z$) and the quadrupole moment $\Theta_{ij}=e\int\rho\left(\boldsymbol{r}\right)\left(3r_{i}r_{j}-\delta_{ij}r^{2}\right)d\boldsymbol{r}$, where $e\rho\left(\boldsymbol{r}\right)$ is the charge distribution within the molecule.111When defining the moments it is customary to take the origin in the center of the positive charge distribution. These moments are of critical importance as they determine the far-field potential produced by the molecule , as evident from the monopole expnasion:222See reference 8; we use the Einstein convention by which repeated Cartesian indices are summed over.

[TABLE]

These low-order ESMs also control the electrostatic interaction energy $W_{es}$ between the molecule (and through it the forces) with a weakly non-constant potential $\varphi^{other}\left(\boldsymbol{r}\right)$ resulting from the other molecules or distant charged sources8:

[TABLE]

where $\varphi_{i}^{other}=\frac{\partial\varphi^{other}}{\partial r_{i}}$ and $\varphi_{ij}^{other}=\frac{\partial^{2}\varphi^{other}}{\partial r_{i}\partial r_{j}}$ (estimated at a central point within the molecule) etc. This pivotal *dual *role of ESMs is what drives the requirement that the charge distribution of the PACs reproduce exactly low-lying molecular ESMs (MOL-ESMs). This point was discussed at length in ref. 9 where the importance of adherence to the ESMs was demonstrated. An efficient elegant method for achieving this in as many as possible moments has been developed 10 although inapplicable for large molecule charges due to numerical instabilities.11

Another source of PACs are the quantum mechanical population analysis (PA) techniques, such as the Mulliken (MPA),12, Loewdin (LPA),13, Hirshfeld (HPA)14, and natural population (NPA)15 analyses. These PAs reflect not only the charge distribution but also aspects of the quantum mechanical wave function. In Fig. 1 (top left) we show as bar-plot the PACs produced by these methods applied to the methyl acetate molecule. It is seen that the different methods produce sometimes significantly different sets of PACs, even PAC signs are not preserved! For example, the LPA assigns positive charges to oxygen atoms, which seems awkward given their high electronegativity. Furthermore, standard PAs do not reproduce the MOL-ESPs closely, as shown in the ESP correlation plot of Fig. 1 (bottom left), where several PAC-ESPs,

[TABLE]

are plotted vs. the MOL-ESP $\varphi\left(\boldsymbol{r}\right)$ calculated from the QM density (Eq. 1) at a grid point $\boldsymbol{r}$. The thin red line in the plot corresponds to the perfectly correlated condition $\varphi^{PAC}=\varphi^{QM}$. In order to quantify the quality of $\varphi^{PAC}\left(\boldsymbol{r}\right)$ we define the mean absolute relative deviation (MARD) from $\varphi\left(\boldsymbol{r}\right)$ as

[TABLE]

where an average is taken over all grid-points $\boldsymbol{r}$ for which: 1) $\boldsymbol{r}$ is “outside of the molecule”, i.e. its distance from any nucleus $a$ is larger than the atomic van-der-Waals radius $R_{a}^{vdW}$ 16) and 2) $\boldsymbol{r}$ is not too far from the molecule, so that its potential $\left|\varphi\left(\boldsymbol{r}\right)\right|$ is not smaller than the threshold value of $e\varphi_{thresh}=0.3eV$.333Note that the expression in Eq. 5 cannot become singular due to this requirement. PACs obtained by “standard” PAs have large MARDs: ranging from 0.37 for HPA up to a whopping 1.76 for NPA. On the right panel of the figure we show data concerning the same molecule, but using the iterated-Hirshfeld method (iHPA),17, 18, 19 the CM5 method 20, which is a parameterized database correction to HPA charges, and the ChElPG method21, which selects PACs that reconstruct the *ab initio *ESP on a set of grid points as close as possible. The latter approach is taken here as representative of a class of methods routinely used for PACs determination. Other members of this method class are the “charge from ESPs” (ChElP)22, the Merz-Kollman23, 16, 24, the charge-restraint ESPs (RESP)25, 26, atomic multipoles ESPs27, in combination with molecular multipoles28 (related to the method proposed here), the dynamical RESP (D-RESP)29 and Hu-Yang fitting30.The iHPA, CM5 and ChElPG methods yield much improved description of the ESP with MARD going from 0.3 for iHPA and CM5 down to 0.08 for ChElPG. Despite the close ESP fit, ChElPG produces PACs that are usually not invariant under transformations preserving the point symmetry of the molecule or under rotations or translations of the nuclei with respect to the real space grid used to perform the fit. Furthermore, in larger molecules the PACs of atoms distant from the molecular surface can become unwieldy large. Both of these issues are discussed in the literature31, 30This instability is likely linked to the fact that the number of parameters derivable from the ESP in a statistically significant way is considerably less than the number of atoms.32 Therefore, iHPA and CM5 are often considered preferred approaches for PACs, although as seen in the figure, both methods leave ample room for improvement. Note that the iHPA charges for this molecule are close to the ChElPG PACs.

Here, we study a new idea: take PACs which are as close as possible to a reference set, for example the MPA, HPA or iHPA PACs, but insist that they reproduce exactly the components of the lowest ESM tensors (dipole and quadrupole) characterizing the molecular charge distribution. We formulate a straightforward method to determine such “minimally-corrected PACs” (section 2) and then benchmark the results using a subset of molecules taken from the database of ref. 20 (section 3). Final conclusions are summarized in section 4. All MPA, HPA and iHPA PACs, as well as the associated MOL-ESPs and MOL-ESMs were computed using developer versions of Q-Chem 4.3 and 4.4 33 at the M06-L DFT level 34 and using the MG3 semi-diffuse (MG3S) basis set 35. This functional/basis set combination was used for developing of the CM5 approach. The CM5, NPA and LPA results were taken from ref. 20.

2 Method

Consider a molecule having $A$ nuclei at given Cartesian positions $\boldsymbol{R}^{a}=\left(R_{x}^{a},\,R_{y}^{a},\,R_{z}^{a}\right)$ ($a=1,\dots,A$), for which a QM calculation has determined the charge density $\rho\left(\boldsymbol{r}\right)$ of the molecule and from it, low order moments the charge $Q$, the dipole $\mu_{i}$ and the symmetric traceless quadrupole moment tensor $\Theta_{ij}$. Note that below, we use the notation $\Theta_{i}^{D}\equiv\Theta_{ii}$ for the diagonal elements of $\Theta$ and $\Theta_{i}^{OD}\equiv\Theta_{jk}$ where $i=x,y,z$ and $ijk$ is a cyclic permutation of $xyz$. For any set of PACs $\boldsymbol{q}=\left(q_{1},\dots,q_{A}\right)$ we define the PAC-ESMs as: the monopole (total charge) $Q^{PAC}\equiv e\sum_{a}q_{a}$, the dipole $\mu_{i}^{PAC}\equiv e\sum_{a}q_{a}R_{i}^{a}$ and the quadrupole $\Theta_{ij}^{PAC}\equiv e\sum_{a}q_{a}\left(3R_{i}^{a}R_{j}^{a}-\delta_{ij}\left(R^{a}\right)^{2}\right)$, ($i,j=x,y,z$) . Given a set of reference PACs $\boldsymbol{q}^{ref}$ we seek to determine a “minimally-corrected” set of PACs $\boldsymbol{q}^{mc}=\boldsymbol{q}^{ref}+\Delta\boldsymbol{q}$ such that that the size of the correction $\left\|\Delta\boldsymbol{q}\right\|^{2}=\Delta\boldsymbol{q}\cdot\Delta\boldsymbol{q}$ is *minimal *but the multipoles are equal to the QM determined multipoles, i.e. the following constraints are satisfied:

[TABLE]

Note that the number of constraints (denoted $C$) in Eq. (6) is $9$ and not $13$ since the electric quadrupole tensor is symmetric and traceless. Point symmetries can reduce this number of constraints further. If, for example, both positive and negative charge densities are symmetric against the reflection through a plane (the x-y plane, for example) then there are 3 constraint less (one from the z component of the dipole and and 2 from XZ and YZ components of the quadrupole, which are zero by symmetry). Only when the number of atoms $A$ in the molecule is greater than the number of constraints $C$ can we hope to reproduce the constraints exactly. We therefore demand that $A>9$ and use the $A-C$ additional “degrees of freedom” to minimize the deviance $\Delta\boldsymbol{q}$. When $2\leq A\leq 9$ we avoid the quadrupole moment constraint and use only the dipole moment constraint.

We are led to consider the Lagrangian

[TABLE]

as a function of the $A$ $q^{\prime}s$ and the ten Lagrange multipliers: one $\lambda,$ three $\lambda_{i}$’s , three diagonal $\lambda_{ii}^{D}$ and three off-diagonal$\lambda_{jk}^{OD}$ where $i=x,y,z$ and $ijk$ is a cyclic permutation of $xyz$. Taking derivatives with respect to these variables and equating to zero leads to the following set of $\left(10+A\right)$ linear equations in $\left(10+A\right)$ unknowns, given here in block-matrix/vector form444Since the matrix is dominated by zero’s one can formulate the linear equation in a more concise way. However, this form is straightforward to derive and manipulate when there are instabilities, discussed later.:

[TABLE]

The $\left(A+10\right)\times\left(A+10\right)$ matrix $S$ is of the following form:

[TABLE]

and depends *only *on the location of the atomic nuclei. The matrix is composed of blocks: the $I_{A\times A}$ block is a $A\times A$ unit matrix, $D_{A\times 3}$, $T_{A\times 3}^{D}$ and $T_{A\times 3}^{OD}$ are matrices of dimension $A\times 3$ ($3$ columns each of length $A$) of matrix elements: $D_{ai}=eR_{i}^{a}$ and $T_{ai}^{D}=e\left(3R_{i}^{a}R_{i}^{a}-\|R^{a}\|^{2}\right)$ for $a=1,\dots,A$ and $i=x,y,z$ and and $T_{ai}^{OD}=3eR_{j}^{a}R_{k}^{a}$ (where the ordered set $i,j,k$ is a cyclic permutation of $x,y,z$). The $D_{3\times A}$, $T_{3\times A}^{D}$, and $T_{3\times A}^{OD}$ blocks are respectively the transposed matrices. The $A+10$ column-vector on the left-hand-side of Eq. 22 includes the unknowns, the $A$ partial charges $q_{a}^{mc}$ and the $\lambda^{\prime}s$, the ten Lagrange multipliers for the ten constraints. The $10+A$ column-vector on the right-hand-side has $A$ values of the reference charges $q_{a}^{ref}$, followed by the total charge on the molecule $Q$, then the three values of the QM dipole moment $\mu_{i}$ followed by the three values of the diagonal elements of the given QM quadrupole tensor $\Theta_{i}^{D}=\Theta_{ii}$ and finally the QM values of the three off-diagonal elements $\Theta_{i}^{OD}=\Theta_{jk}$ where $i=x,y,z$ and $ijk$ is a cyclic permutation of $xyz$). A similar equation holds for the mcD method, where the six last rows are erased from $S$ and from the column vectors and the six right columns are erased from $S$ as well. This leaves us with a $\left(A+4\right)\times\left(A+4\right)$ system of equations.

The structural matrix $S$ may become singular or rank deficient. One trivial source for singularity the use of 3 diagonal constraints while their sum is composed to be zero. The use of the singular-value-decomposition pseudo-inverse 36 for solving Eq. 22 helps to bypass such a singularity. A more delicate source of singularities may arise from symmetry. For example, when the molecule is perfectly planar (or has a plane of symmetry) in the x-y plane then the row corresponding to the dipole in the z direction $D_{az}=eR_{z}^{a}$ must be identically zero and the matrix $S$ will be rank deficient. In this case the $T_{ai}^{OD}$ with $i=x$and $y$ must also be zero). In these cases the SVD pseudoinverse will automatically eliminate constraints that cannot be met due to this kind of symmetry. But for near-symmetrical configurations, instabilities may exist. In cases such as these we can still spot problems by examining the values of the Lagrange multipliers \lambda,$$\lambda_{i} and $\lambda_{ij}$in the solution vector of Eq. 22. The Lagrange multiplier is equal to the derivative of the minimal value of the Lagrangian $L$ with respect to the constraint value ($Q$, $\mu_{i}$ and $\Theta_{ij}$, respectively). Thus if the *ab initio *dipole moment $\mu_{x}$ is given to precision $\delta\mu_{x}$, the product $\left|\lambda_{x}\delta\mu_{x}\right|$is expected to be the error in the minimal value of $L$. Clearly, the minimizing procedure is meaningless unless this error is much smaller than 1. Hence, it is important to eliminate “offending” constraints from the matrix equation (the corresponding row and column in the matrix and the entry in the column vectors) for those having large Lagrange multipliers. We know that *ab initio *multipole properties are usually given to 3 digits hence we eliminate constraints corresponding to Lagrange multipliers large than 1000. The reduced equation is then solved and the remaining Lagrange multipliers are examined again. We repeat such elimination until all Lagrange multipliers have proper magnitudes. This pruning procedure helps avoid cases where small inaccuracies of the input data dominate the final result. Within the molecules studied here such a pruning procedure was used only for few cases of molecules having a near plane symmetry.

When symmetry is active, our procedures reduce the number of constraints and hence the number of *independent *$q_{a}$’s (called number of degrees of freedom (NDOFs)). For example, the water molecule has 3 nuclei but due to symmetry the two H nuclei will have the same PACs and so NDOF=2. Due to the symmetry only the dipole moment in the direction of the $C_{2}$ axis is a constraint (the components perpendicular to the C2 axis are zero by symmetry). Together with the charge of water (0) we already have 2 constraints so one must give up imposing the quadrupole moment for water.

3 Results

To demonstrate the efficacy of the method we show in Fig. 2 the MPA and iHPA PACs and their ESP correlation plots before and after applying the minimal corrections required for imposing dipole and quadrupole moments (denoted mcD/mcDQ-MPA and mcD/mcDQ-iHPA respectively).555Minimally-corrected PACs that reproduce only the dipole moment are designated mcD and those that reproduce the components of the dipole and the quadrupole ESMs are designated mcDQ.

Notice that the MPA-ESP has low correlation with the MOL-ESP, as can be evident visually and also by the reported MARD of 2. The mcD corrections improve the ESP but only mcDQ corrections show high quality ESP (with MARD of 0.05). In accordance with previous reports,19 the iHPA ESP already correlates nicely with the MOL-ESP (MARD of 0.16) but the mcDQ-iHPA improves the correlation significantly and the MARD reduces by a factor of 4. For this molecule, both mcDQ-MPA and mcDQ-iHPA have similar MARDs but this is not typical, for most molecules the mcDQ-iHPA MARDs are much smaller than those of mcDQ-MPA (see Fig.3). The mcDQ-MPA PACs are not drastically different from the MPA PACs yet their MARDs are considerably lower. This shows the power of the minimally-corrected PACs, where a small change in PACs can improve the PAC based ESP considerably.

In Fig. 3 we display a log-scale bar-plot of MARDs of several PAC-based potentials on selected molecules containing 10-18 atoms. Each PAC method can be characterized by a pair of numbers (shown in parenthesis within the legend box) indicating the median/maximal MARD taken over the given set of molecules. The PACs obtained by minimally-correcting the $q\equiv 0$ reference (called 0PA) are actually the minimal PACs that give the dipole and quadrupole of the molecules. It is seen that their correlation with the exact ESP is considerably higher than that of MPA and HPA, somewhat similar to that of mcDQ-MPA and mcDQ-HPA, close to that of CM5. This goes to show that the fit of just the dipole and quadrupole, keeping the charges as small as possible gives a reasonably behaved ESP, although in general, for very large molecules the mcDQ-0PA performance may degrade with size compared to the PA methods. We see that MPA and HPA have similar MARDs while iHPA seems to give considerably smaller MARDs (by a factor 2-3). The minimal corrected (mcDQ) to MPA and HPA yield smaller MARDs by a factor 4 and for iHPA by a factor 2. Altogether the mcDQ significantly improves the ESP. The mcDQ-iHPA median MARD is 7% is similar to that of ChElPG (5%).

It is worthwhile to examine the sensitivity of the MARD estimation with respect to the distance of grid points from the nearest nuclei. In Fig. 3 all sampling grid points were at a distance larger than $1.5\times r_{vdW}$ from any atom. When MARD is estimated using points further way (distance larger than a value of $2\times r_{vdW}$) the iHPA MARD dropped from 0.14 to 0.09 and mcDQ-iHPA MARD dropped from 0.07 to 0.03. ChElPG MARD also reduced, from 0.05 to 0.03. This finding is consistent with the fact that the MCDQ methods provide an asymptotically exact far-field ESP resulting from their reconstruction of the molecular dipole and quadrupole moments.

In Table 1 we show, for each set of PACs the magnitude of the charge correction $\left\|\Delta q\right\|_{\infty}$. For a given molecule the mcD correction is largest for 0PA and then for MPA and HPA and it is smallest for iHPA. mcDQ corrections are in general considerably larger than mcD but in both methods $\left\|\Delta q\right\|_{\infty}$ decreases as the number of atoms in the molecule grows. This is due to the fact that in large systems even small charge shifts have a large affect on the dipole and the quadrupole moments.

In table 2 we summarize the MARD statistics (median and maximal) for for four sets of reference charges: 0PA (reference charges are equal to zero) and MPA, HPA, iHPA. The efficiency of the mc procedure is apparent for MPA, HPA and iHPA, where the mcD reduces the median/maximal MARD by about a factor of 2. mcDQ reduces the MARD further, by a factor of 3 for 0PA and ~2 for MPA and HPA and only 1.1 for iHPA. We thus see that iHPA reconstruction of the ESP strongly benefits from a dipole correction and, interestingly, much less a quadrupole correction.

PACs are sometimes used when molecules distort. In this case, it is important that the they remain continuous under the distortion, so as to enable force calculations. The MPA/HPA/iHPA do not show non-smooth behavior and the mcDQ which is a minimization procedure does not show it as well.666We cannot rule out possible issues if the matrix $S$ of Eq. 23 becomes rank deficient. However, we believe this is an unlikely or quite rare event. In 4 we show the MPA, mcDQ-MPA and ChElPG PACs of the oxygen atom in N-methylethanamide30 the as a function of the dihedral angle $\phi$. It is seen that as the angle increases from 0 the mcDQ-MPA PAC slightly decreases and then increase rapidly followed by a rapid yet continuous drop near $\phi_{c}=0.25$ from a value of $q_{O}=-0.25$ to $q_{O}\approx-0.64$. An additional very sharp feature is seen near $\phi=\pi$. We have checked that this sharp feature is not discontinuous (see inset in Fig. 4) and that the matrix $S$ of Eq. 23 does not become rank deficient. Similar behavior is seen for the PACs of other atoms. We thus conclude that the charges change continuously although sometimes very rapid charge fluctuations can occur.

4 Summary and conclusions

We have studied a new scheme for minimally correcting reference PACs so that they reproduce the exact dipole and quadrupole moments of a molecule and we found that such a minimal correction greatly improves the correlation of the PAC-ESP with respect to the MOL-ESP. The minimal correction scheme does not alter symmetry properties of the reference PACs. Hence, minimally-corrected PACs (mc-PACs) based on MPA, HPA, iHPA fully respect the point-symmetry and rotational/translational symmetries of the molecule.

An additional benefit of the mc-PACs is their stability for inner (or buried) atoms of large molecules. This rises from the stability of the standard population schemes themselves and the fact that mc-PACs involve rather small corrections. As an example, consider the 2-(Dimethylamino)-2-propanol molecule:

for which the ChElPG, HPA and mcDQ-HPA PACs are shown in Fig. 5. Here, ChElPG tends to polarize the molecule: the oxygen and nitrogen share between them a negative unit charge and this is counteracted by the positive unit charge of the central carbon atom C3. On the other hand, MPA assigns a low charge for C3 and spreads rather evenly the remaining positive charge on the 12 terminal hydrogen atoms. mc-MPA charges are very close to those of MPA and thus yet they improve significantly the ESP description for this molecule: the MPA MARD is 0.35 while that of the mc-MPA is 0.1. It is worthwhile to note that the PACs assigned by ChElPG also have a MARD of 0.1.

When the underlying reference is the iHPA set of PACs the resulting ESP is of similar quality to that of the ChElPG set of PACs resulting from a best-fit to ESPs. The dependence of the PACs on the molecular distortion was demonstrated to have sometimes very sharp features however all the changes were smooth, hence forces can be calculated on the atoms of the molecule.

The method here bears a similarity to the optimal point-charge model of Ref. 10 which determines PACs that reproduce as many low-order moments as possible. The crucial difference is best seen when systems grow, model of Ref. 10 would target increasingly higher electrostatic moments as more atoms are included while the present method targets multipoles up to second order and not beyond, thereby avoiding the numerical instabilities described in see Ref. 11. On the other hand. the optimal point-charge model treats the multipole constraints in a more systematic way by minimizing the error over unused moments in the last incomplete spherical shell.

Acknowledgments

Authors express special thanks to Dr. Yihan Shao from Q-CHEM Inc. for his advice and critical assistance in performing the iHPA calculations. We also gratefully acknowledge the support of the Israel Science Foundation Grant No. 189/14.