Letter pubs.acs.org/JPCL
An Extended Charge Equilibration Method Christopher E. Wilmer, Ki Chul Kim, and Randall Q. Snurr* Department of Chemical and Biological Engineering, Northwestern University, 2145 Sheridan Road, Evanston, Illinois 60208, United States S Supporting Information *
ABSTRACT: We present a method for estimating partial atomic charges that uses all of the measured ionization energies (first, second, third, etc.) for every atom in the periodic table. We build on the charge equilibration (Qeq) method of Rappé and Goddard (which used only the first ionization energies) but reduce the number of ad hoc parameters from at least one for every type of atom to just two global parameters: a dielectric strength and a modified parameter for hydrogen atoms. Periodic electrostatic interactions are calculated via Ewald sums, and the partial charges are determined by simultaneously solving a system of linear equations; no iteration is required. We compare the predicted partial atomic charges of this extended charge equilibration (EQeq) scheme against plane-wave density-functional theory derived charges determined via the REPEAT method for 12 diverse metal−organic frameworks (MOFs). We also compare EQeq charges against ChelpG charges calculated using nonperiodic MOF fragments, as well as against Qeq charges as implemented in Accelrys Materials Studio. We demonstrate that for the purpose of ranking MOFs from best to worst for carbon capture applications, EQeq charges perform as well as charges derived from electrostatic potentials, but EQeq requires only a tiny fraction of the computational cost (seconds vs days for the MOFs studied). The source code for the EQeq algorithm is provided. SECTION: Molecular Structure, Quantum Chemistry, and General Theory pattern recognition techniques23 or large tables of fitting parameters27,28 to rapidly estimate partial charges in new systems that are closely related to a “training” set of systems. In contrast to this are recently proposed techniques21,25 that improve in generality (e.g., the REPEAT method for periodic systems21 ) without additional parameters, but employ computationally costly quantum mechanics (but only to obtain partial charges as a precursor step to classical simulations). In this Letter, we propose a general method of estimating partial atomic charges that, in comparison to the charge equilibration (Qeq) method of Rappé and Goddard, has a reduced number of fitting parameters, leading to a simpler model that nonetheless maintains reasonable accuracy. We build on our preliminary work,22 where we first proposed extending the Qeq method to incorporate higher ionization energies. Here we present a complete version of that work that can be applied to periodic systems without adding significant computational expense (source code available in the Supporting Information). We dub this new scheme the extended charge equilibration method (EQeq). In a prior work, Rappé et al.31 described a periodic Qeq (PQeq) method that required iteratively minimizing the system energy in a self-consistent manner. This iterative approach, though computationally expensive and described over 15 years prior to this Letter, is still widely used today.30 By contrast, the EQeq method requires no iteration, which makes it better suited for high throughput screening of materials and for calculating fluctuating charges on
A
ccurately modeling the electrostatic environments of large molecular systems (hundreds to thousands of atoms) has been a long-standing challenge in physical chemistry1−5 of great significance to a diverse spectrum of problems including protein folding,6 DNA melting,7 lyotropic self-assembly of colloids and nanoparticles,8,9 and adsorption of polar species in porous sorbents.10−12 In lieu of using prohibitively computationally expensive quantum mechanics, the dynamics of large systems are often modeled using classical force fields where electrostatic forces are described by pairwise Coulombic interactions between atoms having partial atomic charges. Numerous schemes for estimating partial atomic charges have been proposed over the past four decades,5,13−28 with a few finding very wide adoption.5,14 However, as can be inferred by the number of schemes proposed in the past two years alone,21−23,25,27,28 the existing methods have certain drawbacks, and entirely new classes of materials (e.g., ionic liquids, metal− organic frameworks (MOFs)) reveal deficiencies in the old methods and create a need for new ones. For example, while the split-charge method16 has improved charge assignment in nonbonded organic molecules (or nonbonded fragments of large molecules) by only allowing charge transfer along covalent bonds, Verstraelen et al. found that this failed to account for locally charged functional groups formed by ionic bond dissociation reactions.29 Additionally, neither the splitcharge method nor variations thereof can be applied when the bond topology is unknown or ill-defined, as in the case of screening hundreds of MOFs stored in public databases.30 A defensible strategy of several recently proposed schemes has been to focus on a narrow subset of systems and use either © 2012 American Chemical Society
Received: June 29, 2012 Accepted: August 13, 2012 Published: August 13, 2012 2506
dx.doi.org/10.1021/jz3008485 | J. Phys. Chem. Lett. 2012, 3, 2506−2511
The Journal of Physical Chemistry Letters
Letter
the fly throughout a molecular simulation. For example, the charge-optimized many-body (COMB) approach,32 has shown promise for molecular dynamics (MD) simulations of periodic Si/ SiO2 systems; however, iterative convergence is required to determine partial charges at every MD time-step. In order to validate the EQeq method for the purpose of highthroughput screening of materials, we calculated the partial charges of a diverse set of 12 MOFs and subsequently simulated CO2 adsorption at 0.1 bar and 298 K via grand canonical Monte Carlo methods. These MOFs differ significantly in their chemical composition and pore geometry and, moreover, we were able in each case to compare against experimental CO2 observations. We repeated each CO2 adsorption simulation using REPEAT21 charges, ChelpG5 charges derived from representative MOF clusters by Yazaydin et al.,33 Qeq charges (as implemented in Accelrys Materials Studio34), as well as with no framework charges. We found that EQeq charges were almost as good at identifying high-performance CO2 capture materials as REPEAT and ChelpG charges but required only a tiny fraction of the computational cost. The EQeq method can be concisely expressed in the few equations shown below. A pedagogically-oriented derivation, along with the algorithm code, is provided in the Supporting Information. We begin by assuming a Taylor series form for the energy of an isolated (gas-phase) atom as a function of its charge, EA(Q), centered about some particular partial charge, Q* ∈ , (in the original Qeq method,14 Q* = 0),
where In is the nth ionization energy (the energy required to remove the nth electron after n − 1 have already been removed), and we call χn and Jn the nth electronegativity and nth hardness, respectively (note, I0 is defined here as the energy to go from a charge of −1 to 0, which is equal to the atom’s electron affinity). Since ionization energies (and electron affinities) have been precisely measured for nearly every element in the periodic table35,36 (up to Z = 86), this leads to a quadratic energy function for every atom based entirely on measured quantities, 1 EA (Q ) = EA (Q *) + χQ * (Q − Q *) + JQ *(Q − Q *)2 2 (4)
When atoms are expected to have large positive charges, such as the transition metals in MOFs, the atomic energy given by a secondorder Taylor series centered at Q* = 0 is unlikely to yield accurate results because the energy is based on an extrapolation far from the measured value. The electronegativity of positively charged atoms, as defined in eq 2, is significantly higher than that of neutral atoms (see Figure 1a−c). It is better to extrapolate from the measured atomic energies as little as possible, by choosing an appropriate center for the second-order Taylor series (see Figure 1d−f). If higher order Taylor series terms are incorporated in eq 1, then the atom energy function, now a higher order polynomial, can pass through more than two of the measured points. This would reduce the need to estimate an appropriate center at the cost of requiring iterative solution methods. The energy of the entire system can be expressed as a sum of individual atomic energies and pairwise Coulombic interactions,
⎛ ∂E ⎞ 1 ⎛ ∂ 2E ⎞ EA (Q ) = EA (Q*) + ⎜ ⎟ (Q − Q*) + ⎜ 2 ⎟ (Q − Q*)2 2 ⎝ ∂Q ⎠Q = Q* ⎝ ∂Q ⎠Q = Q*
(1)
+ ...
ESys(Q 1 , Q 2 , ..., Q N )
We neglect higher than quadratic terms so that the final set of equations are linear and can be solved simultaneously. By considering the energy of an atom with one extra electron, EA(Q* − 1), and one less electron, EA(Q* + 1), we obtain the following relations: ⎛ ∂E ⎞ =⎜ ⎟ = χQ * ⎝ ∂Q ⎠Q = Q *
(2)
⎛ ∂ 2E ⎞ IQ *+1 − IQ * = ⎜ 2 ⎟ = JQ * ⎝ ∂Q ⎠Q = Q *
(3)
IQ *+1 + IQ * 2
N
=
∑ ⎛⎝EAk(Q *k ) + χQ * (Q k − Q *k ) ⎜
k
k=1
+
⎞ 1 JQ *(Q k − Q *k )2 + ECk + EOk ⎟ ⎠ k 2
(5)
where ECk is the energy of the charge on the kth atom interacting with the charges of all other atoms in the system, and EOk is a pairwise damping term to prevent infinite charge separation when two atoms are brought arbitrarily close together (see Supporting Information for more details). The Coulombic energy term can be expressed as follows:
where K = 1/4πϵrϵ0, ϵr and ϵ0 are the dielectric strength and permittivity constant, L is the number of “shells” of repeating unit cells to consider, and rkm = |⎯→ rkm + ua⃗ + vb⃗ + wc|⃗ and hkm = ⎯→ ⎯ ⎯→ ⎯ ⎯→ ⎯ ⎯→ ⎯ | hkm + u a−1 + v b−1 + w c −1| are the radial distance between atoms in the real-space and frequency-space lattice, respectively (a⃗,b⃗,c ⃗ are the translational symmetry vectors of the real-space ⎯→ ⎯ ⎯→ ⎯ ⎯→ ⎯ lattice and a−1, b−1, c −1 are of the frequency-space lattice). Note that one may choose to consider a different number of shells in real-space than in frequency-space. The “m≠k*”
notation indicates that the index m skips over the value of k except when considering periodic images (i.e., whenever u, v, or w are not equal to zero). The terms EWr and EWh are the real-space and frequencyspace components of the Ewald sum, erfc EWr (rkm) = 2507
( ) rkm η
rkm
(9)
dx.doi.org/10.1021/jz3008485 | J. Phys. Chem. Lett. 2012, 3, 2506−2511
The Journal of Physical Chemistry Letters
Letter
Figure 1. Measured ionization energies of three elements, (a) carbon, (b) magnesium, and (c) zinc, as well as their associated atom energies (d−f), which are derived from summing the ionization energies (i.e., the energy to acquire a +2 partial charge is the sum of the first and second ionization potential energies). In panels a−c, the measured electron affinity is shown when −1 electrons have been removed. In panels d−f, the circles are measured values, and the red, blue, and green curves represent the approximate quadratic energy function centered around a charge of 0, +1, and +2, respectively. Ionization energies, atom energies, and quadratic fits for the three different centers for all atoms are included in the Supporting Information.
⎛ 4π ⎞ e−( EWh(hkm) = ⎜ ⎟ ⎝Ω⎠
hkmη 2 2
⎯ ⎯→ ) cos( ⎯→ hkm · rkm) 2 hkm
Combining eqs 13 and 11 provides N equations with N unknowns (the charges). Provided the equations in eq 13 depend only linearly on the partial charges, then they can be solved simultaneously, without iteration (in the Supporting Information we show that this is always the case). In simple terms, the electronegativity and Coulombic energy terms in eq 5 favor charge separation, while the hardness terms always resist separation. Since both the Coulombic and hardness terms are of quadratic order, it is possible in general for a certain choice of hardness parameters to result in an energy function that has no minimum (energy decreases as charge is separated without bound). In the Qeq method,14 the hardness term of the hydrogen atoms was made to be of cubic order, which strictly enforced the existence of a global energy minimum on the one hand but required iterative solving on the other. Instead, to maintain a noniterative solution, we use a global dielectric strength, ϵr, to prevent infinite charge separation, which is one of only two ad hoc parameters employed (ϵr = 1.67 used in this work). For high-density systems (e.g., alloys, nonporous solids), a larger dielectric strength may be required. The other ad hoc parameter applies to hydrogen atoms and concerns the unphysical tendency for hydrogen atoms to be assigned negative charges when using the measured ionization potentials. We found that setting I0 = −2 eV for hydrogen (the measured value is +0.754 eV), without changing the form of the energy, led to reasonable charges in all cases. We applied the EQeq method to 12 MOFs (see Table 1) and subsequently calculated CO2 adsorption at 298 K using grand canonical Monte Carlo (GCMC) simulations following standard methods reported elsewhere.10,33 In performing the EQeq calculations, we explicitly considered periodic images contained
(10)
where η is the so-called splitting parameter (with units of distance) that determines how much of the energy to sum in frequency-space versus real-space, and Ω is the volume of the real-space unit cell. The charges Q1,Q2...,QN are determined by finding the minimum energy in eq 5 in a noniterative manner (regardless of whether the system is periodic), subject to the constraint that N
∑ Q k = QT
(11)
k=1
where QT is the net charge of the system (which must be zero in the periodic case). The general approach to minimizing a function of many variables subject to constraints is via Lagrange multipliers, ⎛N ⎞ ∇ESys − λ∇⎜⎜∑ Q k − QT⎟⎟ = 0 ⎝k=1 ⎠
(12)
where the gradient vector ∇ spans the N dimensional space defined by the partial charges. Matching the vector components of eq 12, we get the following N − 1 equations: ∂ESys ∂Q 1
=
∂ESys ∂Q 2
= ··· =
∂ESys ∂Q N
(13) 2508
dx.doi.org/10.1021/jz3008485 | J. Phys. Chem. Lett. 2012, 3, 2506−2511
The Journal of Physical Chemistry Letters
Letter
Table 1. MOFs Investigated in This Study and the Time Required to Obtain the Partial Chargesa MOF name HKUST-1 Pd(2-pymo)2 IRMOF-1 IRMOF-3 Zn-MOF-74 Ni-MOF-74 Co-MOF-74 Mg-MOF-74 ZIF-8 MIL-47 UMCM-150 UMCM-150(N2)
atoms time to time to calculate per unit calculate EQeq ESP (VASP)37,38 cell Q’s (s) (h) 624 252 424 472 162 162 162 162 276 72 354 330
40.6 6 18.2 23 2.8 2.4 2.8 2.6 7.2
