Efficient Evaluation of Molecular Electrostatic Potential in Large Systems

An algorithm for the efficient computation of molecular electrostatic potential is reported. It is based on the partition/expansion of density into (pseudo) atomic fragments with the method of Deformed Atoms in Molecules, which allows to compute the potential as a sum of atomic contributions. These contributions are expressed as a series of irregular spherical harmonics times effective multipole moments and inverse multipole moments, including short-range terms. The problem is split into two steps. The first one consists of the partition/expansion of density accompanied by the computation of multipole moments, and its cost depends on the size of the basis set used in the computation of electron density within the Linear Combination of Atomic Orbitals framework. The second one is the actual computation of the electrostatic potential from the quantities calculated in the first step, and its cost depends on the number of computation points. For a precision in the electrostatic potential of six decimal figures, the algorithm leads to a dramatic reduction of the computation time with respect to the calculation from electron density matrix and integrals involving basis set functions.


Introduction
Despite of its undeniable success in the description of chemical structure and reactivity, the language of chemistry is largely grounded on concepts with loose theoretical foundations. Concepts like bond, atomic charges, orbitals, hybridization, resonance, delocalization, hyperconjugation and many others cannot be directly measured in physical experiments or expressed more precisely in the jargon of quantum mechanics, they are not observables.
The proposal of G.N Lewis [1], in the early twenties of last century, to explain all the previous work made on the study of the molecular structure was recognized from the very beginning as an outstanding contribution to the field. Because of the weak arguments that supported his model, his work was promptly followed by attempts to overcome its theoretical shortcomings by providing a more rigorous foundation based on the recently discovered quantum mechanics. In these attempts, many concepts were proposed to underpin Lewis ideas by pioneers like Coulson, Pauling, Hückel and Mulliken, just to mention some of the most conspicuous figures. Their early efforts to accommodate the complex and sometimes elusive chemical facts in the frame provided by the new emerging physics were driven in part by the available computational capabilities, that were very reduced at the time.
Thus, what at first was an admirable exercise of imagination in its origin has led to a burden of new concepts in an attempt to grasp the increasing knowledge on the complexity of the chemical Several results on performance and accuracy are presented in the fourth section and, finally, significant conclusions are drawn from these results.

Method
The algorithm for the efficient evaluation of MESP reported here is based on DAM [26] partition of electron density. In this partition the molecular density is expressed as a sum of atomic contributions: where the densities of the (pseudo)atomic fragments, ρ A (r A ), are obtained with a least-deformation criterion based on the fast convergence of the long-range multipole expansion of the electrostatic potential [28]. In practice, in the Linear Combination of Atomic Orbitals framework (LCAO), in which DAM is formulated, this fast convergence is achieved by assigning the one-center distributions, i.e., products of pairs of basis functions centered at the same point χ A a (r A ) χ A a (r A ), to their pertaining atoms, A, and partitioning the two-center distributions, χ A a (r A ) χ B b (r B ), between their respective centers: In this work, we will deal with basis functions, χ, consisting of Gaussian contracted functions (CGTO), which are most used in molecular calculations.
Contracted Gaussian functions, χ , are linear combinations of primitive Gaussian functions, g LM (ξ, r): where N G is the number of primitive functions (length of the contraction), index i runs over the primitive functions in the contraction, and the expansion coefficients, c i , are chosen so that the radial part of the contraction remains normalized. Primitive Gaussians on their side are defined as: and the same angular part is taken for all the primitives in a given contraction. Thus, the contracted functions can be written as: where N r ξ i and N Ω LM are radial and angular normalization constants defined by: and z M L (θ, φ) are unnormalized real spherical harmonics given by: where P |M| L (z) are the corresponding associated Legendre functions (see Equation 8.751.1 at [29]).
For two-center densities, g A 00 (ξ A , r A ) g B 00 (ξ B , r B ), consisting of pairs of spherical primitive Gaussians (L = 0) centered at two different points, A and B, it has been proved [28] that the best convergence in the long-range potential is achieved by assigning the whole density to the center with higher exponent ξ or, in case of equal exponents, by assigning one half to each center. This result reminds the nearest-site algorithm reported by A. Stone and M. Alderton [30], but without considering expansions at the bond center.
When applying this criterion to densities made of pairs of spherical contracted functions, the products of primitives become assigned to one center or another depending on their exponents, yielding a partition that can be regarded as a more realistic version of the Mulliken partition (see Figure 1). For distributions consisting of nonspherical functions (L = 0), the partition is applied to the products of brackets of Equation (5), and the remaining terms are translated to the centers as described below.
The next step in the procedure is to expand the densities of the atomic fragments thus obtained as a series of radial factors times unnormalized real regular solid harmonics (hereafter regular harmonics, for short), namely: where the regular harmonics, z m l (r), are related to unnormalized real spherical harmonics by: The partition of density given in Equation (1) combined with the expansion (9) allows us to write the electrostatic potential as [27]: in terms of atomic nuclei charges, ζ A , effective multipoles, Q lm (r), and inverse multipoles, q lm (r): q A lm (r) = (l − |m|)! (l + |m|)! r >r dr z m l (r ) r 2l+1 ρ A (r ) = 4 π 2l + 1 ∞ r dr r ρ A lm (r ) (13) In this way, the molecular electrostatic potential results in a sum of atomic contributions and the short/long-range separation can be carried out at the atomic level. Thus in the long-range region, the effective multipoles can be accurately replaced by point multipoles: (14) and the inverse multipoles do vanish: lim r→∞ q A lm (r) = 0. As it will be shown below, this is a most useful feature when dealing with large systems because a huge fraction of the atomic contributions to MESP can be computed from only the long range terms even in the regions where molecular short-range potential is necessary. This apparent paradox comes from the fact that, in these systems, the molecular short-range at a given point usually involves a reduced number of atomic short-range contributions coming from atoms in the neighborhood of the point, the remaining contributions being of long-range type.

The Algorithm
The algorithm for the evaluation of the electrostatic potential in large molecular systems, according to the method described in the previous section, consists of two main steps, which are executed in sequence. First, a partition of the molecular density must be carried out with the DAM procedure followed by the expansion of the atomic fragments and the computation of effective and inverse multipoles. Second, once the partition/expansion has been made, the electrostatic potential is computed in the desired points with the aid of Equation (11).
It is worth noticing that the first step depends on the size of the basis set used in the computation of molecular density, but it is independent of the number of points where the MESP has to be computed. On the other hand, the second step depends on the number of points for computation, but it is independent of the basis set size.
This decoupling of processes depending on the system size (number of basis functions) from those depending on the number of MESP computation points, combined with the short/long-range separation at atomic level, makes the procedure reported here most efficient to deal with large systems, provided that the MESP is to be evaluated in a not too reduced number of points.

Algorithm for Molecular Density Partition/Expansion
According to DAM partition criterion, the density of atomic fragment at A is given by: where, for two-center CGTO distributions: Θ(x) being the step function: As our purpose is to express the fragments d A ab in coordinates referred to center A and the second primitive in each term of Equation (16) is centered at B, it is necessary to translate the functions g b j to center A. The translation of the exponential factor in g b j can be made in terms of Bessel I functions as proposed by Kaufmann and Baumeister [31]. Working in an aligned frame, with A placed at the origin and B lying on the z axis, i.e., R A = (0, 0, 0), R B = (0, 0, R B ), the translation formula reads: where R B ≡ |R B |, I l+1/2 (z) are the corresponding Bessel functions (see [29] 8.467) and P l (z) are the Legendre polynomials (ibid 8.91).
On the other hand, the remaining factor, which is in essence the regular harmonic z M B L B (r) (apart from a constant factor), can be translated to center A by well known formulas [32]. In the aligned frame, the formula reads: Once the functions are referred to center A, the one-center products of regular harmonics, appearing both in the one-center and in the translated two-center distributions, are decomposed in terms of regular harmonics as described in section, Decomposition of products of regular harmonics, in SF. The final radial factors are identified as the quantities that multiply the corresponding regular harmonics resultant from the decomposition.
In practice, the algorithm proceeds as follows: 1.
The interval 0 ≤ r ≤ 20 bohr is partitioned in n interv subintervals with boundaries corresponding to previously selected values of r that will be noted as λ i , i = 0, 1, ... n interv (currently n interv = 30, see SF for details).

2.
For each interval, the variable r is mapped onto the interval [−1, 1] according to: and a set of values of t is chosen as the zeroes of the Chebyshev T polynomial of order n (currently n = 30) given by (see [33] 22.16.4): 3.
For each center, A, of the system, one-center distributions are expanded as: As described in section, Expansion of one-center distributions, of SF. The radial factors are evaluated in the tabulation points r j,i = (t j + 1) (λ i − λ i−1 )/2 + λ i−1 , multiplied by the ρ aa element of the density matrix, and accumulated.

4.
Likewise, for each center, A, a loop over all the remaining centers B = A is performed. In this loop, for each center B, all the fragments d A ab coming from two-center distributions with one function at A and the other at B, and attributted to center A, are expanded in an aligned frame as a series of regular solid harmonics times radial factors as described in section, One-center expansion of two-center fragments, of SF. The radial factors are evaluated in the tabulation points r j,i = (t j + 1) (λ i − λ i−1 )/2 + λ i−1 , multiplied by theρ ab element of the density matrix, which has been previously rotated to the aligned frame (that is what the tilde means) and accumulated in the aligned frame. Next, the locally accumulated radial factors (i.e., for fixed B) are rotated back to the molecular frame, and the resultant radial factors are further accumulated together with those coming from the one-center distributions and with the radial factors of other pairs of centers to yield the full radial factor ρ A lm (r A ) of Equation (9). Details on rotations of both density matrix and radial factors are given in section, Rotations, in SF.

5.
The tabulations of the ρ A lm (r A ) radial factors are used to decide whether they are negligible or not and, for non-negligible factors, to carry out a numerical projection on Chebyshev T polynomials of variable t in each interval r ∈ [λ i−1 , λ i ]. This projection yields the corresponding piecewise expansion of ρ A lm (r A ). Details of this expansion are given in section Expansion of atomic radial factors in SF. Thus, the final expansion in the i-th interval takes the form: where t is a function of r A , as defined in (20), and the exponential factor is introduced when ρ A 00 (r A ) (leading term in expansion (9)) decays steeply in the interval (see SF for details); otherwise, ζ i = 0 is taken. The number of polynomials taken in the expansion at the i-th interval, n i , is determined on the fly by analyzing the convergence of the projections. The expansion coefficients, c (i) k (l, m), of non-negligible factors are stored in a buffer. An array with a set of suitable pointers to address the coefficients is also generated and stored. 6.
Once the radial factors of expansion have been piecewise expanded, they are used to compute the auxiliary partial integrals: in the same tabulation points, r j,i , as used for the density, as well as the auxiliary constants: and Details are given in section, Effective multipoles from density expansion, in SF. 7.
The tabulations of Q A lm (r j,i ) and q A lm (r j,i ) are used to project these partial integrals onto Chebyshev T polynomials in the same intervals as used for the radial factors of density. In this case, no exponential factor is necessary: The numbers of polynomials in the intervals, n i , are the same as in the corresponding radial factors. In this way, the pointers defined for addressing density expansion coefficients, c  Atomic point multipoles of Equation (14) are obtained by: where the sum runs over the intervals. 9.
Molecular geometry and data corresponding to the tabulation of radial factors are stored in an external binary file with extension damqt, ready to be used for computation of DAM expansion of density. In particular, the following information is stored: number of atoms, number of basis functions and number of shell functions, atomic number and Cartesian coordinates of nuclei, basis set, length of expansion (9) (l max ), and for each center A: pointers to expansion coefficients of radial factors, fitting exponents, ζ i , and expansion coefficients, c (i) k (l, m).

10.
Atomic multipole moments Q A lm , auxiliary quantities Q A lm (λ i ) and q A lm (λ i ), and expansion coefficients b k (l, m) are stored in another external binary file with extension dmqtv. Since the pointers to b k and d k are the same as those used for a k by construction, they do not require to be stored again.

Algorithm for Electrostatic Potential Expansion
Once the files containing the information on MED partition/expansion and the auxiliary quantities for MESP have been generated, MESP can be computed at any desired points using this information. This step is completely independent of the first one, so that the computation can be made as many times as necessary and in different sets of points without requiring repetition of the partition/expansion process.
To compute the MESP, the following algorithm is employed: 1. MED partition/expansion data stored in file damqt are read and stored in memory.

2.
MESP auxiliary data are read from dmqtv, stored in memory and used for computing further auxiliary quantities. In particular, partial accumulated sums: and are computed and stored too.

3.
A double loop over atoms (outer) and tabulation intervals (inner) is performed to determine the length of expansion (11) in each interval and the long-range radius for the atom. This radius is chosen as the lower limit of the interval i, λ i−1 , for which q i,A lm is lower than a user defined long-range threshold.

4.
Next, MESP is computed, running over the atoms, with Equation (11). For points placed in the long-range region, lr, of atom A, the contribution to MESP, V A lr (r A ), is computed in terms of the corresponding atomic point multipoles Q A lm as: For points placed in the short-range region, sr, the contribution is computed by means of: (35) and the quantities Q A lm (r A ) and q A lm (r A ) are obtained in terms of Q i,A lm and q i+1,A lm of Equations (32) and (33), i being the index of the interval such that λ i ≤ r A < λ i+1 , plus the expansions (29) and (30) for the integrals in the interval [λ i , r A ] and [r A , λ i+1 ], respectively.
In all cases, the regular solid harmonics are fast and accurately computed by recursion, as described in section, Recurrence relations of regular solid harmonics, of SF. In the short-range case, eqs (29) and (30) are evaluated with the coefficients b k and d k previously retrieved from file dmqtv and stored in memory, and with the Chebyshev polynomials computed by recursion, as shown in section, Recurrence relations of Chebyshev polynomials, of SF.

5.
If MESP derivatives are wanted, they can be computed together with MESP and using the same auxiliary quantities [34]. The procedure is quoted in section, Computing MESP derivatives, in SF.

6.
Data on basis set and density matrix are only necessary if computation of MESP in terms of nuclear attraction integrals and density matrix, without DAM partition/expansion, is required. As this is an expensive procedure, its usage should be restricted to those cases in which a reference is necessary for testing the accuracy of the algorithm reported here.

Results
To test the performance of the method, we have started by analizing the accuracy of the results attained with the current algorithm. For this purpose, we have computed MESP values for benzene molecule in a set of equally spaced points corresponding to a 129 × 129 × 129 grid in the octant defined by: x, y, z ∈ [0, 20] (length in bohr). These results have been compared with the MESP exact values, V exact , computed using the electron density matrix and the integrals involving basis set functions: where ρ rs are the elements of the density matrix, and nbasis stands for the number of basis functions, χ r , in the LCAO calculation. The results of this comparison are reported in Table 1 for four different lengths in expansion (9), namely: lmax = 5, 10, 15, 20. The highest absolute error, ∆ max , and the root mean square error, rmse, in the grid points significantly decrease with the length of multipole expansion, and they suggest that an expansion with l max = 10 is sufficiently precise for most practical applications. Nevertheless, higher precision can be attained for more demanding applications like topological analysis at a very moderate cost.  Remarkably, for expansions with l max = 10 or higher, a great amount of points are computed with a high number of accurate decimal figures (twelve or more) as is shown in the last row of Table 1. This can be explained because short-range contributions in most points of the grid (those lying outside the molecular volume) are very small or even negligible, and precision is greatly determined by the convergence of the long-range expansion, i.e., the number of terms in (9), the accuracy in the radial factors playing a minor role.
In this regard, it is interesting to check how the algorithm performs in those points in which short-range terms are important. For this purpose, the number of points with a number of decimal accurate figures is plotted in Figure 2 for the four expansions. As it can be seen, the number of points addressed in the curve corresponding to l max = 5 is significantly greater than in the remaining curves. This is due to the insufficient convergence of the long-range MESP in this case, as mentioned before. In the remaining curves, the convergence in the long-range MESP has been achieved to a great extent, and the precision is mainly determined by the short-range terms. As the curves show, also in these cases the precision readily increases with the length of the expansion in two senses: a steady augment of the number of correct decimal figures in the least precise points is observed, and a significant decrease in the number of points for a given precision occurs. Raw data used for the figure can be found in section, Precision of MESP calculation, in SF. The computation wall-clock time with Equation (36) was 5600 seconds, and 2.9, 6.6, 16.2 and 28.0 s for the respective expansions. On the other hand, DAM partition/expansion step with l max = 20 took only 2.2 s for benzene density computed with Dunning's cc-pVDZ basis set [35]. These times were measured on a laptop with processor Intel(R) Core(TM) i7-6700HQ CPU @ 2.60 GHz, running an MPI parallel version of the codes compiled with gfortran, and using 4 processors. Although the calculation was made at HF level, it is important to stress that the computational cost of the full algorithm is independent of the computational level of the molecular density. This is so because the partition/expansion step only depends on the number of elements of the density matrix (i.e., the square of the basis set size), and the MESP computation only depends on the length of expansion (9) and the number of atoms in the system.
It is fairly evident that the algorithm based on DAM yields results that can be sufficient for most practical purposes at a computational time that, in this case, is between two and three orders of magnitude lower than the time required for computation from density matrix and integrals. This is so in a rather small system in which the short-range atomic contributions in the selected grid are about 20% of the total contributions, for a long-range threshold equal to 10 −7 taken in cases of l max = 5, 10, and about 40% for a threshold of 10 −8 taken in cases of l max = 15, 20. As it will be shown below, for really large systems, the fraction of short-range contributions in equivalent grids is much smaller, and the gain in the computational cost increases with the system size with respect to the computation by means of Equation (36). A further test on MESP surface extrema on a density isosurface has been included in the supplementary files.
Once the validity of the results has been established from the point of view of accuracy, we have analyzed the performance of the algorithm in large systems. In Table 2 we collect the results of MESP calculations in a set of molecules ranging from a small system like benzene (12 atoms, 222 basis functions) to two large ones: CC-MMIM BF 4 , consisting of a three circumcoronene slices with two pairs of MMIM BF 4 ionic liquid molecules embedded between the circumcoronene sheets (617 atoms) and a DNA fragment with 750 atoms. In the last two, the electron densities correspond to PM7 calculations [36] (only basis valence functions). With these results, we have analyzed the dependence of the computational cost of the two steps involved. Three different expansions have been used (l max = 10,15,20), except in CC-MMIM FB 4 and DNA fragment. In these cases, due to the ZDO aproximation involved in PM7 method, only valence one-center distributions contribute to densities. Consequenty, the partition/expansion is a very fast process and terms in MED and MESP expansions with l max higher than twice the highest L value in the basis set are zero.
As expected, the computational time of the partition/expansion step in small systems increases with the square of the number of basis functions, and at a lower rate when dealing with large systems, in which cases the dependence tends to be linear. The increment of the cost with respect to the length of expansion (11) is also smaller than the predicted l 2 max dependence. In the supplementary files, a little movie showing the MESP of the DNA fragment over its molecular surface defined as MED isosurface with ρ = 0.001 au is included. Red color means positive MESP values, blue color, negative values.
Finally, in Table 3 we analyze the performance of MPI parallelization of the codes corresponding to both steps of the algorithm. Calculations have been carried out on a polyhidroxilated circumcoronene system with 360 centers, with a basis set consisting of 7560 contracted GTOs. The expansion of the atomic fragments has been made up to l max = 20, and the MESP has been computed with this expansion on a 129 × 129 × 129 cubic grid (2,146,689 points) within an interval x, y, z ∈ [−38, 38]. The wall clock time, the average time per processor and the standard deviation are provided for the DAM partition/expansion and for the MESP tabulation.
In both steps, the computational time scales very well with the number of processors. In the partition/expansion step, a linear fit of the clock wall time vs. the reciprocal of the number of processors gives a regression parameter of R 2 = 0.9998, and the same value is obtained for the fit of the average time. In case of MESP tabulation, the scaling is likewise, with R 2 = 0.9982 and R 2 = 0.9990, respectively. The standard deviations show that the time is more evenly shared among processors in the partition/expansion than in the tabulation, but in both cases the time distribution is satisfactory. Furthermore, in the second step the standard deviation is affected by the fact that the main processor spends more time than the remaining ones, due to the workload associated to the gathering of tabulations accomplished by ancillary processors, which are stored in external files, and to tidying tasks.

Conclusions
A method for the efficient computation of molecular electrostatic potential has been reported. The method splits the problem in two steps. The first one consists of the partition of the molecular electron density in terms of (pseudo)atomic densities, which are further expanded as a truncated series of spherical harmonics times radial factors. The second step corresponds to the computation of the electrostatic potential from these expansions. The computational cost of the first step depends on the size of the basis set used in the LCAO computation and the length of the expansions, and it is independent of the number of points in which the potential is computed. The cost of the second step is independent of the basis set size and only depends on the expansions length and the number of tabulation points. The result is an algorithm that reduces the computational time around two orders of magnitude when compared with the calculation from electron density matrix and integrals involving the basis set, for a precision of 6 decimal figures. The algorithm can be very efficiently parallelized with MPI, showing a good sharing of computational cost among processors in both steps. The algorithm has been implemented in DAMQT package for the analysis of molecular electron density and related properties. A recent version of the package can be downloaded freely from: https://data.mendeley.com/datasets/7mwfftd2x4/2.