2866
J. Chem. Theory Comput. 2009, 5, 2866–2878
Electrostatic Potential Derived Atomic Charges for Periodic Systems Using a Modified Error Functional Carlos Campan˜a´, Bastien Mussard, and Tom K. Woo* Centre for Catalysis Research and InnoVation and Department of Chemistry, UniVersity of Ottawa, Ottawa, Canada, K1N 6N5 Received July 7, 2009
Abstract: A method to generate electrostatic potential (ESP) derived atomic charges in crystalline solids from periodic quantum mechanical calculations, termed the REPEAT method, is presented. Conventional ESP fitting procedures developed for molecular systems, in general, will not work for crystalline systems because the electrostatic potential in periodic systems is ill-defined up to a constant offset at each spatial position. In this work the problem is circumvented by introducing a new error functional which acts on the relative differences of the potential and not on its absolute values, as it is currently done with molecular ESP charge derivation methods. We formally demonstrate that the new functional reduces to the conventional error functional used in molecular ESP approaches when the simulation box of the periodic calculation becomes infinitely large. Several tests are presented to validate the new technique. For the periodic calculation of isolated molecules, the REPEAT charges are found to be in good agreement with those determined with established molecular ESP charge derivation methods. For siliceous sodalite, it is demonstrated that conventional molecular ESP approaches generate ‘unphysical’ charges, whereas the REPEAT method produces charges that are both chemically intuitive and consistent between different periodic electronic structure packages. The new approach is employed to generate partial atomic charges of various microporous materials and compared to both experimentally derived and molecular fragment ESP charges. This method can be used to generate partial atomic charges to be used in simulations of microporous and nanoporous materials, such as zeolites and metal organic framework materials.
I. Introduction Microporous materials are a fascinating class of condensed matter systems with a wide range of important applications. They are commonly employed as catalysts in the petrochemical industry, as adsorbents to trap impurities, as physical reservoirs for gas storage, as molecular sieves for size and shape selective separations, and as mediums for ion exchange processes.1 Perhaps the most well-known microporous materials are zeolites, which are crystalline solids with well-defined structures and pores that range in size from ∼1-20 Å.1 Recently, a new class of microporous materials known as metal organic frameworks (MOFs)2 have emerged that are composed of metal ions and bridging organic ligands. Compared to zeolites, MOFs have a much wider variety of * Corresponding author. E-mail:
[email protected].
structural and chemical motifs that promise a higher degree of tunability. As such, MOFs have attracted significant attention for their potential applications in hydrogen storage and CO2 sequestration.3-9 Atomistic modeling of microporous materials, particulary zeolites, have contributed enormously to increase our understanding on their properties and functionality.10,11 Most studies of the full periodic systems have employed empirical interatomic potentials (force fields), although periodic density functional theory (DFT) studies have recently emerged.11 Molecular dynamics and Monte Carlo simulations using empirical potentials have been successfully applied to investigate the dynamics of guest molecules as well as their adsorption characteristics within microporous materials.10-13 In some cases, atomistic simulations have led to the design of improved compounds with enhanced properties and
10.1021/ct9003405 CCC: $40.75 2009 American Chemical Society Published on Web 08/27/2009
Electrostatic Potential Derived Atomic Charges
performance.14,15 Although several generally transferable force fields have been developed for organic molecules, the treatment of electrostatic interactions still remains a challenge when employing empirical potentials. Most force fields use fixed partial atomic charges to treat the electrostatic interactions, although more sophisticated methods have been developed but are not in wide use. Since there are no ‘true’ atomic charges within polynuclear systems, an assortment of charge derivation methods have been created for various quantitative and descriptive purposes. The so-called electrostatic potential (ESP) derived charges are most commonly used for atomistic simulation of molecular systems.16-19 In order to compute ESP charges, a quantum chemical calculation is performed on a molecule, and partial atomic charges are fit to reproduce the quantum chemical electrostatic potential on a fine grid surrounding the molecule. The grid is chosen to lie outside of the van der Waals radius (or similar) of each atom in the molecule. The charges are ‘designed’ to reproduce the electrostatic potential in the region where it is most important when modeling intermolecular interactions with two-body, additive Coulomb potentials. In situations where there are deeply embedded (buried) atoms in the molecule, the ESP charges for these buried atoms can fluctuate widely, sometimes resulting in nonchemically intuitive values. To deal with this problem, Bayly and co-workers introduced the RESP method,18 where a penalty function is introduced to the fitting procedure to inhibit the ‘unphysical’ charges from arising. ESP charges have been shown to provide superior results compared to those using other schemes, such as charges derived from a Mulliken population analysis, and therefore, ESP charges are the norm for molecular simulation when employing point charge electrostatic models.16,18 For the atomistic simulation of periodic systems, ESP charges are not as universally applied as they are for molecular systems. In highly packed solids, where there is little volume outside of the van der Waals radii of the atoms to define valid fitting points, ESP charges may not be appropriate. However, even for the simulation of microporous materials where there are large pore volumes, Mulliken charges derived from periodic DFT calculations are still in wide usage despite the strong basis set dependence of the scheme.13,20 If ESP charges are used for periodic simulations, they are derived by extracting a fragment of the periodic lattice and performing a molecular ESP calculation.12,21-23 The assumption with this procedure is that the electrostatic potential surrounding the isolated fragment is going to be similar to that of the periodic structure. Considering the care taken to treat long-range electrostatic interactions when using periodic boundary conditions in molecular simulations, this assumption may not be generally valid. Additionally, special considerations also have to be made to ‘cap’ the electronic system to satisfy unfilled valences resulting from the extraction or account for the net charge of the retrieved fragment. Finally, the properties of a periodic wave function including the electrostatic potential are averages over the Brillouin zone. Molecular fragment calculations cannot capture the effect of proper K-point sampling, which is important for some periodic systems.
J. Chem. Theory Comput., Vol. 5, No. 10, 2009 2867
Figure 1. Electrostatic potential plotted as a function of the square of the distance from the cell origin for sodalite generated from three periodic DFT packages: VASP (bottom), CPMD (center), and SIESTA (top). The tabulation is performed on a cubic grid, and only points located outside the van der Waals radii centered on each atom are shown.
Although periodic ‘first principles’ (i.e., DFT) calculations of materials with large unit cells containing hundreds of atoms are now routinely performed, to the best of our knowledge, a method to derive charges from the electrostatic potential of periodic first principles calculations has not been reported.24,25 One of the reasons for the dearth of development of periodic ESP fitting techniques, as compared to molecular systems, may be due to the fact that the absolute energy of an atom is not an intrinsic bulk property, and therefore, the average electrostatic potential is ill-defined.26 In other words, the reference state of the electrostatic potential in a periodic system is arbitrary. This is numerically demonstrated in Figure 1, where the electrostatic potential for siliceous sodalite (a zeolite) is plotted as a function of the square distance to the cell origin on a mesh resulting from three different periodic DFT packages - VASP, CPMD, and SIESTA.27 Qualitatively, the general structure of the electrostatic potentials show the same features, except that the manifold derived from each periodic DFT package is shifted by a different but constant offset amount. Since the average of the electrostatic potential in an infinite system is ill-defined, there is no ‘correct’ offset value, and the straightforward application of the ESP charge fitting procedures developed for molecular systems will result in charges that are dependent on this arbitrary offset value. In this work, we introduce a modified ESP error functional that can be used in periodic systems as well as in molecular systems. Our functional allows for a robust computation of the ESP charges while overcoming the fundamental problem of the electrostatic potential having an ill-defined reference state. We also introduce a new RESP-like penalty function to prevent large fluctuations in the fitted charges of buried atoms that is based on the expansion of the energy of an atom as a function of charge. In the tradition of ESP charge methods, we also suggest an acronym for the charges derived from the method. We refer to these charges as REPEAT (Repeating Electrostatic Potential Extracted ATomic) charges. The remainder of the paper is as follows: The formalism of the method is introduced in detail with the appropriate equations in Section II. Section III describes the choice of parameters used in our implementation of the REPEAT
2868
J. Chem. Theory Comput., Vol. 5, No. 10, 2009
Campaña´ et al.
method as well as details on the periodic DFT calculations that were performed to generate the electrostatic potentials. Section IV contains the tests and the applications used to validate the REPEAT method, and Section V concludes the paper.
II. Methodology For molecular systems, the electrostatic potential at a point br that results from a set of Nq point charges, {qj}, is given in eq 1 (in atomic units), where b rj represents the position of each of the point charges: Nq
φq(rb) )
q
∑ |rb -j rbj|
(1)
j
In generating ESP charges for molecular systems, a quantum chemical calculation is first performed to generate a reference electrostatic potential, φQM(r b). The set of point charges {qj} is then adjusted to minimize the differences between the ‘quantum mechanical’ electrostatic potential, φQM(r b), and the electrostatic potential due to the atomic charges defined by eq 1, for grid points surrounding the atoms in the molecule. This is typically performed in a leastsquares manner where the following function is minimized: F({qj}) )
∑ (φQM(rbgrid) - φq(rbgrid))2 + λ( ∑ qj - qtot) grid
j
(2) The second term in eq 2 represents a Lagrange multiplier that ensures that the sum of atomic charges equals the total charge of the system, qtot. The grid points used for the fitting procedure are chosen to lie outside of the van der Waals (VDW) radii (or similar) of the atoms in the molecule. The actual algorithm for choosing the grid points varies between different methods, and several ESP charge fitting techniques are in wide usage, such as the CHELPG17 and the RESP18 approaches. When generating ESP derived atomic charges for periodic systems, there are additional complications that need to be addressed as compared to molecular ESP charges. First, the electrostatic potential due to the point charges has to account for the infinite periodic images of the atomic charges. The Ewald summation technique is a well-known mathematical procedure to deal with the long-range electrostatic interactions in periodic systems. Within the Ewald procedure, the electrostatic potential is factored in two terms, real and reciprocal space series. The real space series accounts for the effects of the short-range interactions. Conversely, the reciprocal space term includes the long-range effects arising from the presence of the infinite images. For a formal derivation of the mathematical expressions representing each one of the series, we refer the reader to the original and related works.28,29 Summarizing, the exact expression for the potential according to the Ewald formulation is given in eq 3 where the erfc stands for the complementary error function.
φq(rb) )
∑ j,T b
erfc(R|rb - rbj,Tb|) 4π qj + |rb - rbj,Tb| V
k2
ik b(r b-r bj) e 4R2
∑ qje j,k b
k2 (3)
The vectors b T ) n1b a1 + n2b a2 + n3b a3, appearing in the first double summation, map the real space lattice positions. Conversely, the vectors b k ) m1b b1 + m2b b2 + m3b b3, included within the second double summation, map the reciprocal space lattice. The variable R represents the width of the screening Gaussian charge distribution added to each atomic center. Equation 3 can be considered the periodic analogue of eq 1. With infinite terms, eq 3 is mathematically exact and independent of the value chosen for the parameter R. However, the value of R as well as the number of terms in the real and reciprocal space series determine the speed of convergence and accuracy of the numerical computation of the electrostatic potential. Techniques to select optimum values for these parameters have been discussed elsewhere.30 If the charges are to be used in a molecular dynamics simulation, then it can be argued that these Ewald summation parameters used to derive the ESP charges should ideally be the same as those used in the simulations. The second complication when deriving periodic ESP charges relates to the ill-defined nature of the reference state of the electrostatic potential, as depicted in Figure 1. Even if eq 3 is used for φq, minimization of the conventional error function (eq 2) developed for molecular systems will not be of any value when attempting to compute charges in periodic systems. For example, it is clear that the different electrostatic potentials shown in Figure 1 will yield drastically different fitted charges. (There are reports of conventional ESP charge fitting procedures yielding unsatisfactory results for periodic systems.)31,32 To overcome this fundamental problem, we introduce a modified error functional to perform the leastsquares fit. Our modified error functional can be expressed as the sum of two independent terms: F({qj, δφ}) )
∑ (φQM(rbgrid) - (φq(rbgrid) + δφ))2 + grid
λ(
∑ qj - qtot)
(4)
j
where the new parameter δφ is given by δφ )
∑
1 (φ (rb ) - φq(rbgrid)) Ngrid grid QM grid
(5)
As in molecular approaches, the second term in eq 4 represents a Lagrange multiplier that ensures that the total charge of the system remains equal to qtot, which is zero in the case of periodic systems. Examination of the definition of δφ in eq 5 reveals that the ESP fit is performed such that only the relative differences with respect to the average value of the potential on the grid points are considered. This automatically guarantees that the quality of the fit is not affected by a constant offset in the estimation of the potential values. It is important to note that the expression for δφ can be obtained in a simple and straightforward manner by assuming δφ in eq 4 is another independent variable in the error functional F({qj, δφ}). The definition of δφ in eq 5 can
Electrostatic Potential Derived Atomic Charges
J. Chem. Theory Comput., Vol. 5, No. 10, 2009 2869
then be derived by performing the minimization of F not only with respect to the charges but also with respect to δφ. Explicit consideration of δφ within the error functional form not only removes the artifact of an ill-defined reference state but implicitly introduces corrections to the dipole moment and the other higher-order multipoles, if our technique is to be employed to fit charges in molecular compounds. More specifically, it is easy to show that for a neutral molecular system if one expands δφ in the continuum representation up to the first-order (dipole) corrections one obtains δφ )
∑
1 (φ (rb ) - φq(rbgrid)) ≈ ∆φq-QM + Ngrid grid QM grid 1 V
∫
(p bQM - bpq)rˆ |rb | 2
d3r
(6)
∑ (φQM(rbgrid) - (φq(br grid) + δφ))2 +
∑ qj - qtot) + ∑ wj(Ej0 + χjqj + 21 Jj00qj2) j
∂qj)1,...,Nq ) ∂F ∂δφ
(
- δφ - φq(rbgrid)) ×
)
∂φq(rbgrid) ∂δφ + + λ + wj(χj + J00 j qj) ) 0, ∂qj ∂qj
∑ (φ
bgrid) QM(r
- δφ - φq(rbgrid)) ) 0
grid
(8)
The last equation above defines δφ as presented in eq 5. Using the definition of the derivative of the potential with respect to the charges: ∂φq(rbgrid) ) ∂qj
erfc(R|rbgrid - rbj,Tb |) + |rbgrid - rbj,Tb |
∑ b T
4π V
∑ bk
k2 ikb(rbgrid-rbj) e 4R2 e 2
k
) φj′(rbgrid)
(9)
combined with the first group of eq 8 and after some lengthy algebra one arrives at a matrix problem with Nq + 1 unknowns of the type: ˜ )B AQ
(10)
with matrix coefficients Ajm defined as
(7)
∑ grid
{(
(φm′ (rbgrid) -
j
As opposed to the original RESP method and derivative works,18 we have made no ad hoc assumptions regarding the shape of the penalty term. We use a physically motivated penalty multiplier that is an expansion of the energy of an atom (idem to the potential in atomic units) as a function of the atomic charges up to the second order.33,34 The parameters χj and J00 j are the electronegativity and the self-Coulomb interaction, respectively, of the individual chemical elements. The nature of the self-Coulomb interaction term J00 j and how to determine it is discussed in more detail by Rappe and Goddard.34 The free parameter included in the model is wj, which can be considered a weighting factor used to adjust the strength of the restraints. A single weighting factor w could be used for all atoms, or as expressed in eq 7, the weighting factors can be individually adjusted to selectively turn on penalty corrections only for atoms identified as buried prior to performing the charge calculations. Aside from being
bgrid) QM(r
grid
) -2
Ajm,(j,m)
