Charge-Transfer States in Pentacene: Dimer versus Crystal - The

Article pubs.acs.org/JPCC

Charge-Transfer States in Pentacene: Dimer versus Crystal Piotr Petelenz,*,† Mateusz Snamina,† and Grzegorz Mazur‡ †

The K. Gumiński Department of Theoretical Chemistry, and ‡Department of Computational Methods in Chemistry, Faculty of Chemistry, Jagiellonian University, Ingardena 3, 30-060 Kraków, Poland S Supporting Information *

ABSTRACT: According to recent ab initio calculations, the energy gap between the two oppositely polarized charge-transfer (CT) states of a model pentacene dimer is anomalously large, attaining 0.8 eV. Here we introduce the self-consistent charge field approach to evaluate the electrostatic stabilization energies of the pertinent states in the pentacene crystal represented by a dedicated multiscale model system containing the model dimer as its core. We demonstrate that, contrary to common wisdom, the lower of the two CT states is barely affected by the crystalline environment whereas the upper one undergoes a large red shift. Effectively, embedding of the dimer in the crystal bulk reduces the pertinent splitting by an order of magnitude, because most of the intradimer charge−quadrupole interactions are compensated by similar interactions with surrounding molecules. This resolves the apparent contradiction between the ab initio result obtained for the dimer and the splitting of about 0.04 eV resulting from microelectrostatic calculations for the crystal.

■

INTRODUCTION The discovery of spontaneous singlet exciton fission (SF) in solid pentacene prompted widespread interest in the electronic properties of this system, in the recent years triggering a surge of theoretical papers that strive to unravel the underlying physical mechanisms.1−4 Various methodologies have been invoked for this purpose. No matter how vibronic coupling is treated and how the dynamic aspects of the problem are described,1−3 the process is ultimately controlled by the energetics of the electronic states involved. Yet, as the small energy gaps between the relevant eigenstates are inevitably sensitive to the details of the wave function, the requisite accuracy is just on the verge of the predictive power of the most sophisticated quantum chemistry tools available to date. By these standards, the pentacene molecule is large. Moreover, in order to describe the fission process, it is necessary to take into account at least the two molecules on which the two emerging triplets are located, which makes a dimer the minimum model system to be considered in this context. At the ab initio level, this creates a numerical problem of formidable complexity, practically ruling out inclusion of any elements of the crystal environment in which the pertinent pair of molecules is in reality embedded. In a number of papers1−3 focusing on diverse aspects of the fission phenomenon, various versions of the dimer model were adopted. In most cases, parametrization of the model Hamiltonian was in principle based on quantum chemistry calculations but fine-tuned on intuitive grounds to account for simplifications of the model and for potential intrinsic errors of the applied quantum chemistry methods. The paper by Zeng et al.4 stands out from this collection, representing a consistent ab initio approach at a highly © XXXX American Chemical Society

advanced level. However, in view of the accuracy it offers, even minor energy shifts resulting from the simplifications of the underlying model may significantly influence its final findings. Although in many respects one can agree with the authors’ belief that the central conclusions should not be affected, the disregard for the effect of crystalline environment calls for at least a tentative scrutiny. This is the objective of the present contribution. There are at least two specific issues that need to be addressed. First, one of the salient inferences of Zeng et al. is the small value of the charge-transfer (CT) admixture to the (crucial for fission) singlet component of the triplet-pair state (tt). Yet, CT configurations, owing to their large dipole moments, are usually substantially stabilized by immersion in a polarizable medium, such as solid pentacene. This might reduce their energy separation from the neutral tt state, endowed with no dipole moment and expected to be much less sensitive to environmental influence. This provokes the crucial question of how large the crystal-field-induced shifts of the CT configurations really are. Second, the intriguing observation4 that the energies of the diabatic CT configurations of opposite polarities (ac vs ca) dramatically differ, attributed to the charge−quadrupole (CQ) interaction, defies the early microelectrostatic estimates of the latter in solid pentacene,5 later superseded by a more sophisticated treatment6,7 and relied on in recent modeling.3 The question is whether this disagreement ensues from the superior quality of the Zeng et al. calculations, immensely Received: May 20, 2015 Revised: June 1, 2015

A

DOI: 10.1021/acs.jpcc.5b04824 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

molecules or molecular ions separately. The obtained electron densities are attributed to individual atoms using the Merz− Singh−Kollman (MK9) population analysis. This scheme is favored because it was originally designed to reproduce correctly the electrical potential of the pertinent entity, which is the salient quantity in the present context. Other schemes have also been tested and found to yield very similar results. Then, with the relative orientation of the molecules/ions following the desired geometry of the model cluster, one of the standard options of the Gaussian program8 is applied to find the eigenstates of the following: (1) the cation in the electric field of the point atomic charges of the anion and of all surrounding neutral molecules; (2) the anion in the field of the point atomic charges of the cation and of all surrounding neutral molecules; and (3) every neutral molecule of the cluster in the field of both ions and of all surrounding molecules. Subsequently, the obtained electron densities of all the entities contained in the cluster (ions and molecules alike) in the field of their surrounding entities are used to generate the improved sets of atomic charges, corrected to first order by intermolecular interaction. These charges, in turn, are adopted as the environment for each entity in the next step of DFT eigenstate calculations. The procedure is continued until the atomic charges from two consecutive iterations agree. The above scheme is in fact a version of the familiar meanfield approximation. Accordingly, the electrostatic stabilization energy of the cluster is approximated as the difference between the sum of the eigenenergies of all molecules/ions in the field of their entire respective environments and the sum of all pairwise electrostatic interactions between the atomic charges of the molecules/ions constituting the pairs. As is any other version of the mean-field approximation, this recipe is adopted to prevent double counting. By subtracting the corresponding value calculated for the cluster with all molecules in their ground states and neutral (no ions), the energies are referred to the cluster ground state. In analogy to the self-consistent polarization field method (SCPF10−12), which is routine in microelectrostatics, the above procedure will be hereafter referred to as self-consistent charge field (SCCF) approach. Charge Pair Energies in the Crystal. Following our earlier test, in the way of benchmark we applied the SCCF method to the model dimer considered by Zeng et al.,4 using different basis sets. The results indicate that the intermoiety interaction in the pentacene ionic pair is very sensitive to the details of the electronic distribution. The energy difference between the oppositely polarized (ac and ca) CT states predicted by the SCCF approach strikingly changes with the quality of the basis set used in the calculations (cf. Table 1, top panel, columns 3−5). Reasonable convergence is reached at the 6-311G level; the minimum basis turns out to be blatantly insufficient in this context. Therefore, in view of the primitive representation of atomic orbitals in the ZINDO method, it is not surprising that the ZINDO-based TS approach6 grossly underestimated the pertinent energy difference. This should apply a fortiori to the earlier microelectrostatic submolecule approximation (SA5) which was still cruder in describing the electron density. For our present approach in the most advanced basis set, the pertinent difference (0.65 eV) is close to the value obtained from the sophisticated wave functionbased calculations (0.8 eV),4 confirming the validity of the SCCF approach.

improved with respect to earlier attempts, or is an artifact of the restrictive isolated-dimer model. The answer is vital for fission modeling.

■

CALCULATIONS Preliminary Test. As a preliminary step, it is necessary to assess the sensitivity of CT state energies to the details of the electron density distribution in pentacene molecular ions. In the present circumstances, the dimer considered by Zeng et al. is the obvious testing ground for the methodologies used in the past to calculate the CT state energies in molecular crystals.5 The most elaborate of them, introduced by Tsiper and Soos,6 is a self-consistent scheme based on atom−atom polarizabilities evaluated at semiempirical ZINDO level. As a benchmark, we applied it here to calculate the energy of the electrostatic interaction between the pentacene anion (a) and cation (c) in the relative orientation corresponding to the dimer invoked in ref 4. Conforming to the latter, for both ions we adopted the neutral-molecule ground-state geometries. The obtained value of 0.23 eV, when compared to 0.8 eV resulting from the stateof-the-art correlated wave function-based calculations,4 clearly demonstrates the poor performance of the semiempirical approach. Obviously, this should apply a fortiori to the earlier microelectrostatic submolecule approximation (SA5), which was still cruder in describing the electron density, being rooted in the multipole expansion of the molecular electrostatic potential. In the present case the expansion is possibly invalid and certainly impractical, because the discrepancy presented above is presumably due primarily to very large charge− quadrupole terms.4 Consequently, when applied to describe the interaction of the model dimer with the molecules surrounding it in the crystal lattice, neither of these standard methods can produce credible results. This creates demand for a new computational approach, combining the accuracy sufficient to emulate the ab initio results (at least semiquantitatively) with the ability to handle clusters of molecules (including at least the requisite charge pair with all of its nearest neighbors in the pentacene lattice). It would also be desirable that the interaction within the model dimer be treated on the same footing as the interaction between the dimer and its crystalline environment and that moderate computer resources be needed. Self-Consistent Charge Field Method. The approach we propose here is devised to calculate the electrostatic stabilization energies of charge-transfer (CT) states in molecular crystals. It is an extension of the Tsiper and Soos methodology,6 automatically including intramolecular polarization effects and allowing one to take full advantage of the benefits of extended basis sets offered by routine quantum chemistry software. We use the Gaussian package implementation of the density functional theory (DFT).8 The calculations may be performed for a cluster of the size limited by the available computer resources, composed of any kind of molecules and/or ions. In the following, the procedure is described for a single CT state (charge pair) embedded in an environment consisting of several molecules. In the context presented in the introduction, the charge pair is identified with one of the diabatic CT configurations (ac or ca) of the model dimer of ref 4. In the first step the standard DFT program with the B3LYP energy functional is employed to calculate the energy eigenvalues and electron densities for each kind of constituent B

DOI: 10.1021/acs.jpcc.5b04824 J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

■

DISCUSSION Contributions to Energy. In ref 4 the energy difference between the ac and ca configurations was attributed to nondipole contributions to energy, presumably dominated by the charge−quadrupole (CQ) term. As follows from our description of the SCCF method, in the present approach (just as in ref 4) the CQ terms are automatically included in the energies listed in Table 1. In the following they are explicitly invoked merely for interpretational purposes, the specific values being only approximate. In order to estimate the contributions from each of the individual molecules surrounding the central charge pair in the crystal, SCCF calculations were performed for the trimers where a pentacene ion [at the (0, 0, 0) position] was flanked by two neutral molecules either at the (1/2, 1/2, 0) and (−1/2, −1/2, 0) positions or at the (−1/2, 1/2, 0) and (1/ 2, −1/2, 0) positions. Then the corresponding approximate values of the charge−quadrupole interaction energy Q and polarization energy P could be derived from the following reasoning. Consider first a pentacene cation interacting with a pair of neutral pentacene molecules in the geometry appropriate for the (±1/2, ∓ 1/2, 0) neighboring sites; the rest of the crystal is ignored (Figure 1). Following the early literature on the

Table 1. Energies (in eV) of the Dimer CT Configurations (with Respect to Ground State) Calculated in Various Approximationsa Isolated Dimer ZINDOb Eac Eca Eca − Eac

Eac Eca Eca − Eac

I-A-2.47 I-A-2.24 0.23

B3LYP/ STO-3Gc

B3LYP/ 6-311c

I-A-2.44 I-A-2.68 I-A-2.14 I-A-2.01 0.30 0.67 Embedded Dimer

B3LYP/6311++G**c

ref 4

I-A-2.69 I-A-2.04 0.65

2.266 3.063 0.797

ZINDOd

B3LYP/ STO-3Ge

B3LYP/ 6-311e

B3LYP/6311++G**e

I-A-2.73 I-A-2.69 0.04

I-A-2.50 I-A-2.48 0.01

I-A-2.63 I-A-2.61 0.02

I-A-2.68 I-A-2.66 0.02

f

Article

SCCF/ ref 4. 2.28 2.45 0.17

a

I and A denote the molecular ionization potential and electron affinity, respectively. Their specific values are irrelevant in the present context. bCalculated using TS self-consistent approach.6 cSCCF approach. dTS, quoted from ref 7. eSCCF within 10-molecule cluster, TS approach for more distant molecules. fDimer-to-crystal shifts from column 5 applied to the corresponding diabatic energies from ref 4.

This justifies SCCF application for a cluster where the model dimer is surrounded by its nearest neighbors in the pentacene lattice. The cluster for which the calculations have been done consists of 10 molecules, the central 2 being the charge pair of interest. It is of some relevance that in this description the interaction of the model dimer with the molecules surrounding it in the crystal lattice is treated on the same footing as the intradimer interaction, which should prevent potential imbalance-induced artifacts. The interaction of the relevant charge pair with the molecules of the pentacene crystal that are not contained in the model 10-molecule cluster is certainly weaker but may affect the charge-pair energies to some extent. In this work the contribution to CT state energy upon embedding of the cluster in the macroscopic crystal was extracted from the results of Tsiper et al. To this end, the energy of the cluster was recalculated following the ZINDO-based TS scheme6 and subtracted from the CT state stabilization energy derivable from the published effective electron−hole interaction energy.7 The result has been added to the stabilization energy of the 10molecule cluster described above. In this way, diverse methodologies were combined to design a multiscale shell model. The original charge pair (equivalent to that considered in ref 4) became a part of the 10-molecule cluster for which the calculations were done by means of the SCCF approach, with intramolecular density relaxation described at the DFT level and the MK scheme used to generate atomic charges representing the local environment. The long-range contributions, less sensitive to the details of charge distribution, were adopted from the TS results7 that had been obtained for much larger clusters and extrapolated to infinite cluster size. The results are collected in Table 1. They demonstrate that the enormous ac−ca splitting of 0.65 eV noted for the dimer, consistent with the best ab initio result,4 is dramatically suppressed by the crystalline environment. As follows from Table 1 (bottom line), the two energies calculated for the crystal differ only by 0.02 eV. This result is only seemingly counterintuitive, as will become clear below.

Figure 1. Charge-induced-dipole (italics) and charge−quadrupole interactions (boldface) of a pentacene cation/anion with the adjacent neutral molecules in the (1/2, −1/2, 0) (cyan) and (1/2, 1/2, 0) (violet) positions. The ion is at the center.

subject,13 the pentacene molecule is now represented by the lowest terms of the multipole expansion of its charge density, so that the electrostatic stabilization energy of the above system (with respect to the no-charge situation) may be approximated by ΔE = Q − P

(1)

where P is the polarization energy (due to induced dipoles, and quadratic in the charge of the ion), and Q is the charge− quadrupole interaction energy (linear in the charge of the ion). Supposing for the sake of simplicity that in the anion the total charge merely changes sign but its distribution remains roughly the same as in the cation (which is roughly correct), the stabilization energy of the corresponding anionic cluster is the sum of the same polarization energy as for the cation and the charge−quadrupole term with reversed sign. The energies of the cationic and anionic trimer, independently calculated by the SCCF method, combine to provide a set of two equations, enabling one to evaluate individually the Q and P terms. The corresponding values for the neighboring (±1/2, ±1/2, 0) sites are readily obtainable from a similar inference. The results are visualized in Figure 1. C

DOI: 10.1021/acs.jpcc.5b04824 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Strictly speaking, the above arguments apply merely to the ion−neutral−molecule interactions, not to the interion interaction that is operative in the model dimer. As can be proved, however (cf. Supporting Information), only that part of the quadrupole moment of the ions which is equal to the quadrupole moment of the neutral molecule contributes to the ac−ca splitting, so only the corresponding contribution to energy needs to be discussed here. To start with, it is worth noting that in the absence of the crystalline environment the ultimate value of the difference between the energies of the ac and ca states of the model dimer4 is amplified by orientation effects peculiar to its specific geometry. Some components of the quadrupole moment tensor are positive and some are negative; the relative position of the ions defines which components are exposed by one ion to the electric field of the other. In the configuration labeled ac in ref 4, the interaction of the electric field originating from the anion with the set of cation atomic charges is stabilizing. At the same time, the electric field of the cation, interacting with the atomic charges of the anion, also yields a stabilizing contribution, despite the reversed charge (Figure 2a). The two contributions add, giving a large overall stabilizing effect. The effect is reversed for the ca pair (Figure 2b), giving rise to the reported enormous splitting of the two diabatic configurations.4 Upon embedding of the dimer in the crystal, all of these intradimer charge−quadrupole terms are counterbalanced by similar interactions with the new neighbors of the two ions, as shown in the diagram (Figure 2a, b). On this account, the (upper) ca configuration is stabilized, while the (lower) ac configuration is destabilized by roughly the same amount (to a large extent, this coincides with the unpublished results recently obtained by the authors of ref 4). On the other hand, the polarization contributions (mediated by the charge-induceddipole term) due to the new neighbors are always negative (Figure 2c), so that they stabilize both CT states. This leads to enormous net stabilization of the (originally much higher) ca configuration (where the additional polarization concurs with the decrease of the net CQ contribution), whereas the (originally lower) energy of the ac configuration barely changes with respect to its position in the absence of crystalline surroundings, because the two contributions practically cancel each other. Overall, if all simplifying assumptions adopted above would be strictly valid, the two dimer configurations, when embedded in the crystal, should have very similar energies, just as the simplified treatments7,5 suggest. Actually, their energy separation of 0.02 eV obtained here (in complete SCCF calculations, without the simplifications needed above for interpretational purposes) is quite close to the TS result of 0.04 eV.7 Validity of Various Approximations. It is tempting to improve the above result by taking full advantage of the multiscale model described earlier, whereby the reasonably accurate SCCF description of crystal field effects can be combined with the excellent accuracy of the sophisticated calculations4 for the model dimer. This is done by directly applying the SCCF dimer-to-crystal shifts calculated here to the diabatic CT state energies of ref 4. The values obtained in this way are listed in the last column of Table 1 (bottom panel). Qualitatively, the result at this level of approximation is reminiscent of that in the isolated dimer: the energy of the ca CT state still perceptibly exceeds that of the ac CT state. However, the difference of about 0.17 eV (cf. Table 1) is now almost 5 times smaller. Admittedly, the large discrepancy

Figure 2. Leading electrostatic interactions in the model cluster encompassing the CT dimer. The black-and-gray shapes mimic the molecules. The wedges, representing the interactions, point from the charges (C, circles) to the quadrupoles (Q, squares) or induced dipoles (ID, rectangles) with which the charges interact. The widths of the wedges are proportional to the interaction energies; the green/red color codes negative/positive contributions. CQ terms: (a) ac configuration; (b) ca configuration. The intradimer CQ contributions cumulate, both having the same sign. In the crystal, they acquire the counterparts of the same size, but opposite sign. (c) Polarization (CID) contributions (all values negative; valid for both configurations).

between the splitting estimated in this way and our original result (0.02 eV) may reflect the inherent accuracy of the D

DOI: 10.1021/acs.jpcc.5b04824 J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

■

underlying quantum chemistry methods and should not be surprising in view of the quality of Zeng et al. wave functions, definitely superior to our DFT approximation. We are inclined to suppose, though, that it is more likely to be an artifact due to incompleteness of the compensation between the CQ interactions within the dimer (described by high-grade wave function-based theory) and the corresponding terms engaging the surrounding molecules (described within the less advanced DFT approach). In this regard, it is crucial that in the SCCF description the interaction of the model dimer with the surrounding molecules is treated on the same footing as the intradimer interaction.

■

Article

ASSOCIATED CONTENT

* Supporting Information S

Electrostatic characteristics of the ions and their effect on the ac−ca splitting; complete list of authors for ref 8. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b04824.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +48 12 6632212. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We express our gratitude to P. Kubisiak for his inestimable advice concerning the approximations and software used here, to M. Andrzejak for access to his unpublished results on pentacene relaxation energies, and to T. Zeng, R. Hoffmann, and N. Ananth for their comments and for information concerning their new results prior to their publication. This research was performed with equipment purchased thanks to the financial support of the European Regional Development Fund in the framework of the Polish Innovation Economy Operational Program (contract no. POIG.02.01.00-12-023/08).

CONCLUSIONS

Taken at face value, our findings highlight the importance of accurate representation of the electron density in the calculations of electrostatic energies in crystals and are in acute conflict with common wisdom, where polarization of the surrounding molecules would be expected to stabilize all CT states of the model dimer to more or less the same extent. In contrast, the present accurate results reveal that in the pentacene case the crystal-field stabilization of the ca charge configuration is anomalously large, whereas the ac configuration is barely stabilized at all. This behavior is detectable only by applying relatively advanced quantum chemical methodology. The same evidence, however, may be looked upon from a different angle. As a matter of fact, the predicted huge crystalfield stabilization is merely undoing the effect of the very large ac−ca splitting predicted for the model dimer, which within the simplistic SA or TS description would not emerge at all. Effectively, the predicted ac−ca splitting in the crystal is almost the same, irrespective of the level of sophistication of the approach applied. This may be rationalized by taking into account the fact that the polarization term, stabilizing a CT configuration in the crystal, is negative-definite, and for this reason not very sensitive to the details of the wave function. In contrast, the charge−quadrupole term may have either sign, with the latter sensitively depending both on the details of the electron density distribution in the ions and on the relative orientation of the molecules with which the charges interact. Owing to crystal symmetry, the charge−quadrupole contributions of different signs to a large extent compensate each other, obliterating the energetic effects due to subtle details of the molecular electron density distribution. Effectively, at crystal scale the earlier microelectrostatic approaches5−7 may be viewed as a valid coarse-grained picture. By and large, our present results reconcile the enormous ac− ca splitting reported by Zeng et al. for the isolated model dimer with the modest value that emerges for the model dimer embedded in the crystal. This value is not an artifact of the simplistic microelectrostatic approach,6,5 as might be surmised, but is a genuine consequence of the interplay between the intradimer interactions and the crystal field effect. This implies that if the dimer model is to be used to study singlet fission in the crystal, the energies of its CT states have to be treated as effective values (as was done in refs 1−3), rather than evaluated by quantum chemistry calculations for an isolated pair of molecules.

■

REFERENCES

(1) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. Microscopic Theory of Singlet Exciton Fission. I. General Formulation. J. Chem. Phys. 2013, 138, 114102. (2) Berkelbach, T. C.; Hybertsen, M. S.; Reichman, D. R. Microscopic Theory of Singlet Exciton Fission. II. Application to Pentacene Dimers and the Role of Superexchange. J. Chem. Phys. 2013, 138, 114103. (3) Beljonne, D.; Yamagata, H.; Brédas, J. L.; Spano, F. C.; Olivier, Y. Charge-Transfer Excitations Steer the Davydov Splitting and Mediate Singlet Exciton Fission in Pentacene. Phys. Rev. Lett. 2013, 110, 226402. (4) Zeng, T.; Hoffmann, R.; Ananth, N. The Low-Lying Electronic States of Pentacene and Their Roles in Singlet Fission. J. Am. Chem. Soc. 2014, 136, 5755−5764. (5) Bounds, P. J.; Siebrand, W.; Eisenstein, I.; Munn, R. W.; Petelenz, P. Calculation and Spectroscopic Assignment of Charge-Transfer States in Solid Anthracene, Tetracene, and Pentacene. Chem. Phys. 1985, 95, 197−212. (6) Tsiper, E. V.; Soos, Z. G. Charge Redistribution and Polarization Energy of Organic Molecular Crystals. Phys. Rev. B 2001, 64, 195124. (7) Tsiper, E. V.; Soos, Z. G. Electronic Polarization in Pentacene Crystals and Thin Films. Phys. Rev. B 2003, 68, 085301. (8) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (9) Singh, U. C.; Kollman, P. A. An Approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5, 129− 145. (10) Jurgis, A.; Silinsh, E. A. On the Interaction of Electrons and Holes in a Molecular Crystal. Phys. Status Solidi B 1972, 58, 735−743. (11) Silinsh, E. A. Organic Molecular Crystals: Their Electronic States; Springer: Berlin, 1980. (12) Mazur, G. An Improved SCPF Scheme for Polarization Energy Calculations. J. Comput. Chem. 2008, 29, 988−993. (13) Bounds, P. J.; Munn, R. W. Polarization Energy of a Charge in a Molecular Crystal. II. Charge−Quadrupole Energy. Chem. Phys. 1981, 59, 41−45.

E

DOI: 10.1021/acs.jpcc.5b04824 J. Phys. Chem. C XXXX, XXX, XXX−XXX