- Email: [email protected]
The fragment molecular orbital method and understanding monomer polarization
Chemical Physics Letters 554 (2012) 185–189
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
The fragment molecular orbital method and understanding monomer polarization Cassandra D.M. Churchill ⇑ Department of Chemistry, University of Alberta, 11227 Saskatchewan Drive, Edmonton, Alberta, Canada T6G 2G2
a r t i c l e
i n f o
Article history: Received 16 August 2012 In final form 20 September 2012 Available online 1 October 2012
a b s t r a c t The magnitudes of select p–p interactions are studied with MP2, as well as the fragment molecular orbital (FMO) method, which evaluates the total energy as a sum of polarized monomer energies and interfragment interaction energies. The monomer polarization is of particular interest and is found to increase in magnitude with system size and the electrostatic contribution to the interaction. As a result, differences will exist in the interaction energy reported from FMO calculations, and those calculated using a traditional supermolecular approach. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Noncovalent interactions are important contributors to the structure and stability of molecular complexes, such as proteins and nanostructures. However, accurate determination of these interactions within large assemblies presents a computational challenge, particularly with ab initio methods [1,2]. Approaches such as the fragment molecular orbital (FMO) method have been developed to perform accurate computations on macromolecules [3]. This method structurally decomposes large assemblies into smaller, computationally-manageable components [4]. Subsequently, the total energy is calculated as a sum of the polarized fragment energies and interfragment interaction energies (IFIEs) [5]. This method has been successfully applied to a variety of problems [2,6], including the investigation of protein mutations in a system containing over 36 000 atoms [7]. Given its promise and applications in the literature, it is important to understand how the FMO method varies from other approaches of quantifying noncovalent interactions, before applying FMO to study noncovalent interactions in large molecular systems. The supermolecular approach is the most popular for examining the strength of noncovalent interactions in small systems, where the interaction energy is calculated as the difference in energy between the complex and its constituents [8]. This approach may even be extended to calculate many-body interactions. For example, in a simple tetramer not only can the total interaction energy be calculated, but also two-, three- and four-body contributions, allowing additivity to be investigated. However, with such an approach, effects like monomer polarization are implicitly included in the interaction energy, rather than explicitly calculated [5]. The FMO method has two advantageous features that are of interest: (1) the ability to study noncovalent interactions in large
⇑ Fax: +1 (780) 492 8231. E-mail address: [email protected] 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.09.047
molecular systems; and (2) the explicit calculation of monomer polarization. Therefore, the goal of this Letter is to understand the role monomer polarization plays in noncovalent interactions and how FMO results differ from those obtained using a supermolecular approach. This will allow the FMO method to be applied with confidence in studies that examine noncovalent interactions in protein-containing systems, or other macromolecules. The present Letter applies the FMO method to quantify p–p interactions in benzene systems. The benzene geometries of interest have been previously studied in the literature [9–11]. In addition, these molecular clusters are noncovalently bound and therefore do not require the fragmentation of bonds, which will allow the interactions to also be evaluated using a supermolecular approach. Furthermore, it is expected that the benzene monomers will experience little polarization from the surrounding benzene molecules since the electrostatic component of these interactions is small [12]. Results from these test systems provide a measure of the minimum polarization one may expect in noncovalentlybound systems. Ultimately, the present Letter will offer insight into the effect of monomer polarization within larger systems, as well as an enhanced understanding of the FMO method.
2. Methods Structures of the benzene-containing dimers, trimers and tetramers were obtained from the literature [9,10]. For brevity, the dimer geometries in this Letter have been renamed from those in Ref. [9] according to the following: 2S = dimer S; 2PD = dimer PD; 2T = dimer T. Similarly, the trimer and tetramer geometries in this Letter have been renamed from Ref. [10]: 3A = trimer S; 3B = trimer D; 3C = trimer T2; 3D = trimer T1; 3E = trimer C; 4A = tetramer S; 4B = tetramer PD; 4C = tetramer T. Binding energies were first evaluated at MP2/6-31G(d,p) using a supermolecular approach for an N-body system (DEsup ):
186
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
DEsup ¼ EðNÞ
N X Ei
ð1Þ
i¼1
where EðNÞ is the energy of the N-body system and Ei is the energy of an isolated monomer. Interaction energies were then evaluated with the FMO method at the MP2/6-31G(d,p) level of theory, designated FMO-MP2/ 6-31G(d,p). This comparative study examines the ability of FMO-MP2 to reproduce MP2 energies with the same basis set. Although it is well known that MP2/6-31G(d,p) is not the most accurate approach for calculating p–p interaction energies [1], this method was chosen with the ultimate goal of studying protein– ligand interactions, where a balance between computational efficiency and accuracy will be necessary. The FMO method calculates fragment energies self-consistently in a Coulomb field generated by surrounding fragments to create a converged electrostatic potential, such that each fragment is polarized by its environment [4]. The polarization is the sum of both destabilizing (due to mutual polarization) and stabilizing (due to electrostatic interactions between polarized charge distributions) components. However, in molecular clusters, the stabilizing component is equal to the negative of two times the destabilizing component, resulting in the total polarization being related to the negative of the destabilization component [13]. Therefore, polarization discussed in the present Letter refers to only the destabilization component. Subsequently, interfragment interactions are calculated and added to the polarized monomer energies to obtain the total energy of the system, EFMO [4]. Energy obtained via FMO2 ðEFMO2 Þ includes dimer interactions [5]:
EFMO2 ¼
N N N X X X E0i þ E0ij E0i E0j þ TrðDDij Vij Þ i¼1
i
ð2Þ
i
In Eq. (2), E0i is the energy of a polarized monomer, E0ij represents the energy of a polarized dimer and TrðDDij Vij Þ denotes the interaction of the relaxed density change with the embedding electrostatic potential [14]. DDij is the difference in the density matrix between a dimer and its constituent monomers [14]. Vij stands for the electrostatic potential from the surrounding fragments acting on a dimer [14]. If trimer interactions are included, EFMO3 is obtained [5]:
EFMO3 ¼ EFMO2 þ
N h i X DE0ijk þ DEDijk
ð3Þ
i
DE0ijk ¼ E0ijk E0i E0j E0k E0ij E0ik E0jk
ð4Þ
DEDijk ¼ TrðDDijk V ijk Þ TrðDDij Vij Þ TrðDDik Vik Þ TrðDDjk Vjk Þ
ð5Þ
In the equations above, DE0ijk is the interaction energy of a polarized trimer and DEDijk represents the interaction of the relaxed density with the embedding electrostatic potential [14]. Interaction energies were evaluated from the FMO-calculated energy of the system ðEFMO2 or EFMO3 Þ in two ways, which differ subtly in their interpretation. In one case, binding was evaluated with respect to the isolated monomers in their free states [2]:
DEFMO ¼ EFMO
N X Ei
This energy represents the interaction between monomers within their current environment, and it is this quantity that is reported by FMO in the GAMESS package. The FMO2 results were further probed by examining the various contributions to the interaction energies using the Pair Interaction Energy Decomposition Analysis (PIEDA) scheme. However, only DE0FMO2 is analyzed, rather than DEFMO2 , since this is the quantity obtained by FMO in GAMESS. The electrostatic ðDE0ES Þ, exchangerepulsion ðDE0EX Þ and dispersion ðDE0DI Þ contributions, as well as charge transfer and higher-order mixed terms ðDE0CTþmix Þ are reported [5,13].
DE0FMO2 ¼ DE0ES þ DE0EX þ DE0DI þ DE0CTþmix
Given that the default implementation of FMO employs approximations for calculating IFIEs [4,5,18,19], all FMO calculations were FMO performed both with ðDEFMO approx Þ and without ðDEfull Þ these approximations when evaluating interaction energies. Specifically, approximations corresponding to the RESPPC, RESDIM and RCORSD keywords were disabled, as well as RITRIM for FMO3 calculations [5]. All calculations were performed using the GAMESS-US (Version October 2010 R1) program [15]. The program Facio was used to prepare input files and analyze output [16,17]. 3. Results and discussion Initially, benzene-containing dimers were considered. The interactions calculated using the supermolecular approach and FMO2 (with and without approximations) are included in Table 1. It is expected FMO2 will provide exact results for dimer systems since pair-wise interactions are the only many-body interactions that exist in dimers. Indeed, binding energies calculated with reFMO2 spect to the isolated monomers ðDEsup ; DEFMO2 approx ; and DEfull Þ are identical. These results indicate FMO2 provides exact MP2/631G(d,p) energies for these systems. However, examination of FMO2 the output for DEFMO2 calculations shows that the approx and DEfull small interfragment distances for the benzene dimers are within the FMO cutoffs, and therefore the interaction is explicitly evaluated and the implementation of approximations has no effect for these systems. FMO treatment of the benzene dimers reveals small polarization of the monomers (<1 kJ mol1, Table S1). Therefore, the binding energies calculated with respect to the polarized monomers 0FMO2 ðDE0FMO2 Þ differ from those calculated with respect to approx and DEfull FMO2 the isolated monomers ðDEFMO2 approx and DEfull Þ. This originates from whether monomer polarization is included in the reported binding energy. For example, calculations for dimer 2T show DEFMO2= 13.7 kJ mol1 and DE0 FMO2 = 15.1 kJ mol1 with a difference of 1.4 kJ mol1. However, Fragment 1 is polarized (destabilized) by Table 1 MP2/6-31G(d,p) benzene dimer interaction energies (kJ mol1). Structurea
N X E0i i¼1
ð7Þ
2PD
2T
9.8 9.8
14.3 14.3
13.7 13.7
c,e DE0FMO2 approx
8.4
14.5
15.1
f,d DEFMO2 full f,e DE0FMO2 full
9.8
14.3
13.7
8.4
14.5
15.1
c,d DEFMO2 approx
i¼1
DEFMO ¼ EFMO
2S
DEsup b
ð6Þ
This energy represents the change in stability of the isolated monomers upon complex formation. Alternatively, interaction energies were evaluated with respect to the polarized monomers:
ð8Þ
a
Geometries shown in Figure 1. Interaction energy calculated according to Eq. (1) using a supermolecular approach. c FMO results obtained with (default) approximations. d Calculated with respect to isolated monomers according to Eq. (6). e Calculated with respect to polarized monomers according to Eq. (7). f FMO results obtained without approximations. b
187
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
0.55 kJ mol1 while Fragment 2 is polarized (destabilized) by 0.92 kJ mol1 (Table S1), for a total monomer destabilization of 1.47 kJ mol1. Therefore, binding energies calculated with respect to isolated monomers ðDEsup and DEFMO2 Þ inherently include the energy difference resulting from monomer polarization, whereas binding energies calculated with respect to the polarized monomers ðDE0FMO2 Þ do not. These results for dimers clearly show the relationship between binding energies calculated with respect to isolated ðDEsup and DEFMO2 Þ and polarized ðDE0FMO2 Þ monomers. It should be noted that it is DE0FMO2 that is reported in the output for FMO calculations in GAMESS. Next, benzene trimers were considered (Table 2). For such systems, FMO3 should provide exact energies since it includes up to three-body interactions. Binding results for FMO3 with respect to FMO3 isolated monomers ðDEFMO3 approx and DEfull Þ and the supermolecular approach ðDEsup Þ show that near-identical energies are obtained, where minor discrepancies (<1 kJ mol1) may be attributed to tighter convergence in FMO than in ab initio calculations. As seen in the FMO2 calculations of dimer systems, the application of FMO3 to trimer systems yields exact energies. Small differences exist between the FMO2 energies and the MP2/6-31G(d,p) energies (up to 0.6 kJ mol1), where the use of approximations further increases these deviations (up to 1.1 kJ mol1). Unlike what is found in the dimers systems, the distances between non-nearest neighbor fragments (Fragments 1 and 3) in trimers 3A to 3D (Figure 1) fall outside the implemented cutoffs. Therefore, approximations are included in the evaluation of the FMO energies for the benzene trimers. The use of approximations leads to minimal deviations FMO between DEFMO approx and DEfull , which are slightly larger for FMO2 than FMO3 (0.6 and 0.2 kJ mol1, respectively). Within the benzene trimers, monomer polarization is as large as 2 kJ mol1 in magnitude, which is greater than that observed for fragments within the dimers (Table S1). Larger polarization is observed for the middle fragment (Fragment 2, Table S1) in structures 3A–3D, which is engaged in two nearest-neighbor interactions. The largest polarizations are observed for 3E (Table S1), which is the only genuine three-body system (i.e. each monomer has two nearest-neighbor contacts). These results illustrate that the number of nearest-neighbor monomers influences the magnitude of the calculated monomer polarization. Different binding energies result depending on whether they are calculated with respect to isolated ðDEFMO Þ or polarized ðDE0FMO Þ monomers. As in the case of the dimers, this difference can be rationalized. For example, results for 3E show a difference in binding energy of 6.0 kJ mol1 depending
on which monomer energies are used (DEFMO3 = 43.4 kJ mol1 versus DE0 FMO3 = 49.4 kJ mol1, Table 2). This difference can be explained by the fragment polarization, which has a destabilizing effect of 2.0 kJ mol1 for each monomer (Table S1). Overall, this indicates that the discrepancy between binding energies calculated with respect to isolated and polarized monomers increases with system size.
Table 2 MP2/6-31G(d,p) benzene trimer interaction energies (kJ mol1).
Table 3 MP2/6-31G(d,p) benzene tetramer interaction energies (kJ mol1).
Figure 1. Geometry and naming of the benzene (a) dimers, (b) trimers, and (c) tetramers, where fragments are numbered from top to bottom.
Structurea 3B
3C
3D
3E
19.7 19.3
29.4 28.3
27.5 26.9
26.9 26.9
43.5 43.2
DEsup b
c,e DE0FMO2 approx
16.4
28.8
29.6
29.2
49.2
f,d DEFMO2 full f,e DE0FMO2 full c,d DEFMO3 approx 0FMO3 c,e DEapprox f,d DEFMO3 full f,e DE0FMO3 full
19.6
28.8
27.1
26.9
43.2
16.7
29.4
29.8
29.2
19.6
29.3
27.5
16.7
30.0
30.2
19.8
29.3
16.9
30.0
DEsup b c,d DEFMO2 approx
a
Structurea
3A
4A
4B
4C
29.6 28.7
44.6 42.3
40.9 40.2
c,e DE0FMO2 approx
24.3
43.2
43.8
f,d DEFMO2 full
29.4
43.5
40.4
49.2
f,e DE0FMO2 full
25.0
44.4
44.1
27.0
43.4
c,d DEFMO3 approx
29.3
44.4
40.9
29.3
49.4
c,e DE0FMO3 approx
24.9
45.4
44.5
27.5
27.0
43.4
f,d DEFMO3 full
29.7
44.5
40.9
30.2
29.3
49.4
f,e DE0FMO3 full
25.3
45.4
44.5
Geometries shown in Figure 1. Interaction energy calculated according to Eq. (1) using a supermolecular approach. c FMO calculation performed with (default) approximations. d Calculated with respect to isolated monomers according to Eq. (6). e Calculated with respect to polarized monomers according to Eq. (7). f FMO calculation performed without approximations. b
c,d DEFMO2 approx
a
Geometries shown in Figure 1. Interaction energy calculated according to Eq. (1) using a supermolecular approach. c FMO calculation performed with (default) approximations. d Calculated with respect to isolated monomers according to Eq. (6). e Calculated with respect to polarized monomers according to Eq. (7). f FMO calculation performed without approximations. b
188
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
Figure 2. A selection from the S66 database containing p–p interactions, including (a) benzene–benzene, (b) benzene–pyridine, (c) benzene–uracil, (d) pyridine–pyridine, (e) pyridine–uracil and (f) uracil–uracil dimers. The change in monomer energy (kJ mol1) upon polarization, as calculated by FMO, is indicated next to each monomer.
Table 4 MP2/6-31G(d,p) dimer interaction energies (kJ mol1) of a select subset of the S66 database. Structurea benzenebenzene
benzenepyridine
benzeneuracil
pyridinepyridine
pyridineuracil
uraciluracil
14.0 14.0
17.3 17.3
27.6 27.6
20.1 20.1
32.3 32.3
46.9 46.9
c,e DE0FMO2 approx
13.9
17.8
32.4
20.8
38.2
55.4
f,d DEFMO3 full
14.0
17.3
27.6
20.1
32.3
46.9
f,e DE0FMO3 full
13.9
17.8
32.4
20.8
38.2
55.4
DEsup b c,d DEFMO2 approx
a b c d e f
Geometries shown in Figure 2. Interaction energy calculated according to Eq. (1) using a supermolecular approach. FMO calculation performed with (default) approximations. Calculated with respect to isolated monomers according to Eq. (6). Calculated with respect to polarized monomers according to Eq. (7). FMO calculation performed without approximations.
The benzene-containing systems were further extended to tetramers (Table 3). However, neither FMO2 nor FMO3 will fully account for the four-body interactions in these systems. By comparing FMO results calculated with respect to the isolated monomers ðDEFMO2 and DEFMO3 full full ) to the supermolecular binding energies, DEsup , one can see that the FMO2 results (deviations up to 1.1 kJ mol1) can be improved by FMO3 (deviations up to 0.1 kJ mol1). As seen for the trimers, the use of approximations tends to worsen results, with greater disagreement between FMO DEFMO being observed for FMO2 (deviations up to approx and DEfull 1 1.2 kJ mol ) than FMO3 (deviations up to 0.4 kJ mol1). The results for the benzene tetramers again illustrate that fragment polarizations increase with system size. As observed for the trimers, the middle fragments (Fragments 2 and 3, Table S1) experience greater polarization, while the outer fragments are less affected. However, a comparison of monomer polarization between trimers and tetramers reveals a given fragment is most affected by the presence of a nearest neighbor, and less affected by nextnearest-neighbor fragments. For example, identical intermonomer geometries exist within trimer 3A and tetramer 4A. Within 3A, Fragment 2 is polarized by –1.4 kJ mol1 from two nearest-neighbor fragments (Fragments 1 and 3). In comparison, within 4A, Fragment 2 is also polarized by –1.4 kJ mol1, but in the presence of two nearest-neighbor fragments (Fragments 1 and 3) and a nextnearest-neighbor fragment (Fragment 4). Therefore, a next-nearest-neighbor monomer has little influence on monomer polarization within these systems. The tetramer binding energies differ depending on whether they are calculated with respect to isolated or polarized monomer energies. The differences reach 4.4 kJ mol1, reflecting the increased polarization of fragments within a tetramer. This discrepancy between binding energies calculated with respect to isolated and polarized monomers increases with system size for two reasons: (1) as the number of nearest-neighbor contacts increases, the magnitude of fragment polarization increases; and (2) as the number of polarizable fragments that contribute to the binding energy increases, the deviations can accumulate (Table 3).
The systematic study of the benzene-containing systems detailed above illustrates the effect of fragment polarization on binding energies for systems of increasing size in which the electrostatic components are small. However, similar effects on fragment polarization will also be observed for systems with larger electrostatic components, which will lead to discrepancies in the calculation of binding between either DEFMO or DEsup and DE0FMO . Therefore, in addition to the above systems, a subset of the S66 database was chosen to consider interactions with larger electrostatic contributions [20]. Specifically, interactions in p–p dimers containing: (1) benzene–benzene; (2) benzene–pyridine; (3) benzene–uracil; (4) pyridine–pyridine; (5) pyridine–uracil; and (6) uracil–uracil were evaluated at their equilibrium geometries (Figure 2) [20]. Although this test set also contains a stacked benzene–benzene dimer (Figure 2), its geometry differs from that of 2S (Figure 1). Differences between the electrostatic contribution to the interactions ðDE0ES Þ for the benzene systems and the S66 subset can be seen in the PIEDA results (Figures S2, S3 and Tables S2, S3). As found for the benzene dimers above, the calculated interaction energies of the S66 subset with respect to the isolated monomers (Table 4) show that the same energy is obtained from both FMO2-MP2 and MP2. Furthermore, approximations are not invoked for these small systems since interfragment distances are within the defined cutoffs. However, large discrepancies are found when interaction energies are evaluated with respect to isolated ðDEFMO2 Þ and polarized ðDE0FMO2 Þ monomers. These discrepancies range from 0.0 kJ mol1 for the benzene–benzene dimer, up to 8.4 kJ mol1 for the uracil–uracil dimer. Indeed, monomer polarization is greatest for the uracil homodimer (Figure 2f). These results illustrate that monomer polarization and discrepancies in binding energies will also increase in magnitude as the electrostatic contribution to the interaction increases. 4. Conclusions The present comparative study examines p–p interaction energies of benzene-containing dimers, trimers and tetramers, as well
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
as a subset of the S66 database using both a supermolecular approach and the fragment molecular orbital method. Results show that both MP2/6-31G(d,p) and FMO-MP2/6-31G(d,p) predict nearly identical absolute energies for the systems considered. Deviations can occur depending on whether the calculation of binding energies is performed with respect to isolated ðDEFMO Þ or polarized ðDE0FMO Þ monomers. Deviations between DE0FMO and DEFMO tend to increase as the number of monomers increases and the electrostatic contribution to the interactions increases. Therefore, caution is suggested in the direct comparison of binding energies reported in the FMO output ðDE0FMO Þ and those calculated with respect to isolated monomers ðDEFMO or DEsup Þ. This is especially true for systems that are large and/or have substantial electrostatic contributions to interaction energy, since the magnitude of fragment polarization is greater in both these instances. Acknowledgments This work was supported by the Natural Sciences and Engineering Research Council (NSERC CGS-D), Alberta Innovates-Technology Futures, and the University of Alberta. I am grateful to Dr. Mariusz Klobukowski for his support and input during the preparation of this manuscript. Thanks are also given to a Reviewer for his or her expertise and suggested improvements to this work. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2012. 09.047.
189
References [1] P. Hobza, Accounts of Chemical Research 45 (2012) 663. [2] D.G. Fedorov, T. Nagata, K. Kitaura, Physical Chemistry Chemical Physics: PCCP 14 (2012) 7562. [3] M.S. Gordon, D.G. Fedorov, S.R. Pruitt, L.V. Slipchenko, Chemical Reviews 112 (2012) 632. [4] D.G. Fedorov, K. Kitaura, J. Chem. Phys. 120 (2004) 6832. [5] D. Fedorov, K. Kitaura, in: D. Fedorov, K. Kitaura (Eds.), The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems, CRC Press, New York, 2009, p. 5. [6] D.G. Fedorov, K. Kitaura (Eds.), The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems, CRC Press, 2009. [7] A. Yoshioka et al., J. Mol. Graphics Modell. 30 (2011) 110. [8] P. Hobza, Annu. Rep. Prog. Chem. C: Phys. Chem. 107 (2011) 148. [9] M.O. Sinnokrot, E.F. Valeev, C.D. Sherrill, Journal of the American Chemical Society 124 (2002) 10887. [10] T.P. Tauer, C.D. Sherrill, Journal of Physical Chemistry A 109 (2005) 10475. [11] M.O. Sinnokrot, C.D. Sherrill, Journal of Physical Chemistry A 110 (2006) 10656. [12] S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, K. Tanabe, Journal of the American Chemical Society 124 (2001) 104. [13] D.G. Fedorov, K. Kitaura, Journal of Computational Chemistry 28 (2007) 222. [14] D.G. Fedorov, K. Kitaura, in: E.B. Starikov, J.P. Lewis, S. Tanaka (Eds.), Modern Methods for Theoretical Physical Chemistry of Biopolymers, Elsevier, 2006, p. 3. [15] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, Journal of Computational Chemistry 14 (1993) 1347. [16] M. Suenaga, Journal of Computer Chemistry, Japan 4 (2005) 25. [17] M. Suenaga, Journal of Computer Chemistry, Japan 7 (2008) 33. [18] T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa, Y. Akiyama, M. Uebayasi, K. Kitaura, Chemical Physics Letters 351 (2002) 475. [19] D.G. Fedorov, K. Ishimura, T. Ishida, K. Kitaura, P. Pulay, S. Nagase, Journal of Computational Chemistry 28 (2007) 1476. [20] J. Rezac, K.E. Riley, P. Hobza, Journal of Chemical Theory and Computation 7 (2011) 2427.
Contents lists available at SciVerse ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
The fragment molecular orbital method and understanding monomer polarization Cassandra D.M. Churchill ⇑ Department of Chemistry, University of Alberta, 11227 Saskatchewan Drive, Edmonton, Alberta, Canada T6G 2G2
a r t i c l e
i n f o
Article history: Received 16 August 2012 In final form 20 September 2012 Available online 1 October 2012
a b s t r a c t The magnitudes of select p–p interactions are studied with MP2, as well as the fragment molecular orbital (FMO) method, which evaluates the total energy as a sum of polarized monomer energies and interfragment interaction energies. The monomer polarization is of particular interest and is found to increase in magnitude with system size and the electrostatic contribution to the interaction. As a result, differences will exist in the interaction energy reported from FMO calculations, and those calculated using a traditional supermolecular approach. Ó 2012 Elsevier B.V. All rights reserved.
1. Introduction Noncovalent interactions are important contributors to the structure and stability of molecular complexes, such as proteins and nanostructures. However, accurate determination of these interactions within large assemblies presents a computational challenge, particularly with ab initio methods [1,2]. Approaches such as the fragment molecular orbital (FMO) method have been developed to perform accurate computations on macromolecules [3]. This method structurally decomposes large assemblies into smaller, computationally-manageable components [4]. Subsequently, the total energy is calculated as a sum of the polarized fragment energies and interfragment interaction energies (IFIEs) [5]. This method has been successfully applied to a variety of problems [2,6], including the investigation of protein mutations in a system containing over 36 000 atoms [7]. Given its promise and applications in the literature, it is important to understand how the FMO method varies from other approaches of quantifying noncovalent interactions, before applying FMO to study noncovalent interactions in large molecular systems. The supermolecular approach is the most popular for examining the strength of noncovalent interactions in small systems, where the interaction energy is calculated as the difference in energy between the complex and its constituents [8]. This approach may even be extended to calculate many-body interactions. For example, in a simple tetramer not only can the total interaction energy be calculated, but also two-, three- and four-body contributions, allowing additivity to be investigated. However, with such an approach, effects like monomer polarization are implicitly included in the interaction energy, rather than explicitly calculated [5]. The FMO method has two advantageous features that are of interest: (1) the ability to study noncovalent interactions in large
⇑ Fax: +1 (780) 492 8231. E-mail address: [email protected] 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.09.047
molecular systems; and (2) the explicit calculation of monomer polarization. Therefore, the goal of this Letter is to understand the role monomer polarization plays in noncovalent interactions and how FMO results differ from those obtained using a supermolecular approach. This will allow the FMO method to be applied with confidence in studies that examine noncovalent interactions in protein-containing systems, or other macromolecules. The present Letter applies the FMO method to quantify p–p interactions in benzene systems. The benzene geometries of interest have been previously studied in the literature [9–11]. In addition, these molecular clusters are noncovalently bound and therefore do not require the fragmentation of bonds, which will allow the interactions to also be evaluated using a supermolecular approach. Furthermore, it is expected that the benzene monomers will experience little polarization from the surrounding benzene molecules since the electrostatic component of these interactions is small [12]. Results from these test systems provide a measure of the minimum polarization one may expect in noncovalentlybound systems. Ultimately, the present Letter will offer insight into the effect of monomer polarization within larger systems, as well as an enhanced understanding of the FMO method.
2. Methods Structures of the benzene-containing dimers, trimers and tetramers were obtained from the literature [9,10]. For brevity, the dimer geometries in this Letter have been renamed from those in Ref. [9] according to the following: 2S = dimer S; 2PD = dimer PD; 2T = dimer T. Similarly, the trimer and tetramer geometries in this Letter have been renamed from Ref. [10]: 3A = trimer S; 3B = trimer D; 3C = trimer T2; 3D = trimer T1; 3E = trimer C; 4A = tetramer S; 4B = tetramer PD; 4C = tetramer T. Binding energies were first evaluated at MP2/6-31G(d,p) using a supermolecular approach for an N-body system (DEsup ):
186
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
DEsup ¼ EðNÞ
N X Ei
ð1Þ
i¼1
where EðNÞ is the energy of the N-body system and Ei is the energy of an isolated monomer. Interaction energies were then evaluated with the FMO method at the MP2/6-31G(d,p) level of theory, designated FMO-MP2/ 6-31G(d,p). This comparative study examines the ability of FMO-MP2 to reproduce MP2 energies with the same basis set. Although it is well known that MP2/6-31G(d,p) is not the most accurate approach for calculating p–p interaction energies [1], this method was chosen with the ultimate goal of studying protein– ligand interactions, where a balance between computational efficiency and accuracy will be necessary. The FMO method calculates fragment energies self-consistently in a Coulomb field generated by surrounding fragments to create a converged electrostatic potential, such that each fragment is polarized by its environment [4]. The polarization is the sum of both destabilizing (due to mutual polarization) and stabilizing (due to electrostatic interactions between polarized charge distributions) components. However, in molecular clusters, the stabilizing component is equal to the negative of two times the destabilizing component, resulting in the total polarization being related to the negative of the destabilization component [13]. Therefore, polarization discussed in the present Letter refers to only the destabilization component. Subsequently, interfragment interactions are calculated and added to the polarized monomer energies to obtain the total energy of the system, EFMO [4]. Energy obtained via FMO2 ðEFMO2 Þ includes dimer interactions [5]:
EFMO2 ¼
N N N X X X E0i þ E0ij E0i E0j þ TrðDDij Vij Þ i¼1
i
ð2Þ
i
In Eq. (2), E0i is the energy of a polarized monomer, E0ij represents the energy of a polarized dimer and TrðDDij Vij Þ denotes the interaction of the relaxed density change with the embedding electrostatic potential [14]. DDij is the difference in the density matrix between a dimer and its constituent monomers [14]. Vij stands for the electrostatic potential from the surrounding fragments acting on a dimer [14]. If trimer interactions are included, EFMO3 is obtained [5]:
EFMO3 ¼ EFMO2 þ
N h i X DE0ijk þ DEDijk
ð3Þ
i
DE0ijk ¼ E0ijk E0i E0j E0k E0ij E0ik E0jk
ð4Þ
DEDijk ¼ TrðDDijk V ijk Þ TrðDDij Vij Þ TrðDDik Vik Þ TrðDDjk Vjk Þ
ð5Þ
In the equations above, DE0ijk is the interaction energy of a polarized trimer and DEDijk represents the interaction of the relaxed density with the embedding electrostatic potential [14]. Interaction energies were evaluated from the FMO-calculated energy of the system ðEFMO2 or EFMO3 Þ in two ways, which differ subtly in their interpretation. In one case, binding was evaluated with respect to the isolated monomers in their free states [2]:
DEFMO ¼ EFMO
N X Ei
This energy represents the interaction between monomers within their current environment, and it is this quantity that is reported by FMO in the GAMESS package. The FMO2 results were further probed by examining the various contributions to the interaction energies using the Pair Interaction Energy Decomposition Analysis (PIEDA) scheme. However, only DE0FMO2 is analyzed, rather than DEFMO2 , since this is the quantity obtained by FMO in GAMESS. The electrostatic ðDE0ES Þ, exchangerepulsion ðDE0EX Þ and dispersion ðDE0DI Þ contributions, as well as charge transfer and higher-order mixed terms ðDE0CTþmix Þ are reported [5,13].
DE0FMO2 ¼ DE0ES þ DE0EX þ DE0DI þ DE0CTþmix
Given that the default implementation of FMO employs approximations for calculating IFIEs [4,5,18,19], all FMO calculations were FMO performed both with ðDEFMO approx Þ and without ðDEfull Þ these approximations when evaluating interaction energies. Specifically, approximations corresponding to the RESPPC, RESDIM and RCORSD keywords were disabled, as well as RITRIM for FMO3 calculations [5]. All calculations were performed using the GAMESS-US (Version October 2010 R1) program [15]. The program Facio was used to prepare input files and analyze output [16,17]. 3. Results and discussion Initially, benzene-containing dimers were considered. The interactions calculated using the supermolecular approach and FMO2 (with and without approximations) are included in Table 1. It is expected FMO2 will provide exact results for dimer systems since pair-wise interactions are the only many-body interactions that exist in dimers. Indeed, binding energies calculated with reFMO2 spect to the isolated monomers ðDEsup ; DEFMO2 approx ; and DEfull Þ are identical. These results indicate FMO2 provides exact MP2/631G(d,p) energies for these systems. However, examination of FMO2 the output for DEFMO2 calculations shows that the approx and DEfull small interfragment distances for the benzene dimers are within the FMO cutoffs, and therefore the interaction is explicitly evaluated and the implementation of approximations has no effect for these systems. FMO treatment of the benzene dimers reveals small polarization of the monomers (<1 kJ mol1, Table S1). Therefore, the binding energies calculated with respect to the polarized monomers 0FMO2 ðDE0FMO2 Þ differ from those calculated with respect to approx and DEfull FMO2 the isolated monomers ðDEFMO2 approx and DEfull Þ. This originates from whether monomer polarization is included in the reported binding energy. For example, calculations for dimer 2T show DEFMO2= 13.7 kJ mol1 and DE0 FMO2 = 15.1 kJ mol1 with a difference of 1.4 kJ mol1. However, Fragment 1 is polarized (destabilized) by Table 1 MP2/6-31G(d,p) benzene dimer interaction energies (kJ mol1). Structurea
N X E0i i¼1
ð7Þ
2PD
2T
9.8 9.8
14.3 14.3
13.7 13.7
c,e DE0FMO2 approx
8.4
14.5
15.1
f,d DEFMO2 full f,e DE0FMO2 full
9.8
14.3
13.7
8.4
14.5
15.1
c,d DEFMO2 approx
i¼1
DEFMO ¼ EFMO
2S
DEsup b
ð6Þ
This energy represents the change in stability of the isolated monomers upon complex formation. Alternatively, interaction energies were evaluated with respect to the polarized monomers:
ð8Þ
a
Geometries shown in Figure 1. Interaction energy calculated according to Eq. (1) using a supermolecular approach. c FMO results obtained with (default) approximations. d Calculated with respect to isolated monomers according to Eq. (6). e Calculated with respect to polarized monomers according to Eq. (7). f FMO results obtained without approximations. b
187
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
0.55 kJ mol1 while Fragment 2 is polarized (destabilized) by 0.92 kJ mol1 (Table S1), for a total monomer destabilization of 1.47 kJ mol1. Therefore, binding energies calculated with respect to isolated monomers ðDEsup and DEFMO2 Þ inherently include the energy difference resulting from monomer polarization, whereas binding energies calculated with respect to the polarized monomers ðDE0FMO2 Þ do not. These results for dimers clearly show the relationship between binding energies calculated with respect to isolated ðDEsup and DEFMO2 Þ and polarized ðDE0FMO2 Þ monomers. It should be noted that it is DE0FMO2 that is reported in the output for FMO calculations in GAMESS. Next, benzene trimers were considered (Table 2). For such systems, FMO3 should provide exact energies since it includes up to three-body interactions. Binding results for FMO3 with respect to FMO3 isolated monomers ðDEFMO3 approx and DEfull Þ and the supermolecular approach ðDEsup Þ show that near-identical energies are obtained, where minor discrepancies (<1 kJ mol1) may be attributed to tighter convergence in FMO than in ab initio calculations. As seen in the FMO2 calculations of dimer systems, the application of FMO3 to trimer systems yields exact energies. Small differences exist between the FMO2 energies and the MP2/6-31G(d,p) energies (up to 0.6 kJ mol1), where the use of approximations further increases these deviations (up to 1.1 kJ mol1). Unlike what is found in the dimers systems, the distances between non-nearest neighbor fragments (Fragments 1 and 3) in trimers 3A to 3D (Figure 1) fall outside the implemented cutoffs. Therefore, approximations are included in the evaluation of the FMO energies for the benzene trimers. The use of approximations leads to minimal deviations FMO between DEFMO approx and DEfull , which are slightly larger for FMO2 than FMO3 (0.6 and 0.2 kJ mol1, respectively). Within the benzene trimers, monomer polarization is as large as 2 kJ mol1 in magnitude, which is greater than that observed for fragments within the dimers (Table S1). Larger polarization is observed for the middle fragment (Fragment 2, Table S1) in structures 3A–3D, which is engaged in two nearest-neighbor interactions. The largest polarizations are observed for 3E (Table S1), which is the only genuine three-body system (i.e. each monomer has two nearest-neighbor contacts). These results illustrate that the number of nearest-neighbor monomers influences the magnitude of the calculated monomer polarization. Different binding energies result depending on whether they are calculated with respect to isolated ðDEFMO Þ or polarized ðDE0FMO Þ monomers. As in the case of the dimers, this difference can be rationalized. For example, results for 3E show a difference in binding energy of 6.0 kJ mol1 depending
on which monomer energies are used (DEFMO3 = 43.4 kJ mol1 versus DE0 FMO3 = 49.4 kJ mol1, Table 2). This difference can be explained by the fragment polarization, which has a destabilizing effect of 2.0 kJ mol1 for each monomer (Table S1). Overall, this indicates that the discrepancy between binding energies calculated with respect to isolated and polarized monomers increases with system size.
Table 2 MP2/6-31G(d,p) benzene trimer interaction energies (kJ mol1).
Table 3 MP2/6-31G(d,p) benzene tetramer interaction energies (kJ mol1).
Figure 1. Geometry and naming of the benzene (a) dimers, (b) trimers, and (c) tetramers, where fragments are numbered from top to bottom.
Structurea 3B
3C
3D
3E
19.7 19.3
29.4 28.3
27.5 26.9
26.9 26.9
43.5 43.2
DEsup b
c,e DE0FMO2 approx
16.4
28.8
29.6
29.2
49.2
f,d DEFMO2 full f,e DE0FMO2 full c,d DEFMO3 approx 0FMO3 c,e DEapprox f,d DEFMO3 full f,e DE0FMO3 full
19.6
28.8
27.1
26.9
43.2
16.7
29.4
29.8
29.2
19.6
29.3
27.5
16.7
30.0
30.2
19.8
29.3
16.9
30.0
DEsup b c,d DEFMO2 approx
a
Structurea
3A
4A
4B
4C
29.6 28.7
44.6 42.3
40.9 40.2
c,e DE0FMO2 approx
24.3
43.2
43.8
f,d DEFMO2 full
29.4
43.5
40.4
49.2
f,e DE0FMO2 full
25.0
44.4
44.1
27.0
43.4
c,d DEFMO3 approx
29.3
44.4
40.9
29.3
49.4
c,e DE0FMO3 approx
24.9
45.4
44.5
27.5
27.0
43.4
f,d DEFMO3 full
29.7
44.5
40.9
30.2
29.3
49.4
f,e DE0FMO3 full
25.3
45.4
44.5
Geometries shown in Figure 1. Interaction energy calculated according to Eq. (1) using a supermolecular approach. c FMO calculation performed with (default) approximations. d Calculated with respect to isolated monomers according to Eq. (6). e Calculated with respect to polarized monomers according to Eq. (7). f FMO calculation performed without approximations. b
c,d DEFMO2 approx
a
Geometries shown in Figure 1. Interaction energy calculated according to Eq. (1) using a supermolecular approach. c FMO calculation performed with (default) approximations. d Calculated with respect to isolated monomers according to Eq. (6). e Calculated with respect to polarized monomers according to Eq. (7). f FMO calculation performed without approximations. b
188
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
Figure 2. A selection from the S66 database containing p–p interactions, including (a) benzene–benzene, (b) benzene–pyridine, (c) benzene–uracil, (d) pyridine–pyridine, (e) pyridine–uracil and (f) uracil–uracil dimers. The change in monomer energy (kJ mol1) upon polarization, as calculated by FMO, is indicated next to each monomer.
Table 4 MP2/6-31G(d,p) dimer interaction energies (kJ mol1) of a select subset of the S66 database. Structurea benzenebenzene
benzenepyridine
benzeneuracil
pyridinepyridine
pyridineuracil
uraciluracil
14.0 14.0
17.3 17.3
27.6 27.6
20.1 20.1
32.3 32.3
46.9 46.9
c,e DE0FMO2 approx
13.9
17.8
32.4
20.8
38.2
55.4
f,d DEFMO3 full
14.0
17.3
27.6
20.1
32.3
46.9
f,e DE0FMO3 full
13.9
17.8
32.4
20.8
38.2
55.4
DEsup b c,d DEFMO2 approx
a b c d e f
Geometries shown in Figure 2. Interaction energy calculated according to Eq. (1) using a supermolecular approach. FMO calculation performed with (default) approximations. Calculated with respect to isolated monomers according to Eq. (6). Calculated with respect to polarized monomers according to Eq. (7). FMO calculation performed without approximations.
The benzene-containing systems were further extended to tetramers (Table 3). However, neither FMO2 nor FMO3 will fully account for the four-body interactions in these systems. By comparing FMO results calculated with respect to the isolated monomers ðDEFMO2 and DEFMO3 full full ) to the supermolecular binding energies, DEsup , one can see that the FMO2 results (deviations up to 1.1 kJ mol1) can be improved by FMO3 (deviations up to 0.1 kJ mol1). As seen for the trimers, the use of approximations tends to worsen results, with greater disagreement between FMO DEFMO being observed for FMO2 (deviations up to approx and DEfull 1 1.2 kJ mol ) than FMO3 (deviations up to 0.4 kJ mol1). The results for the benzene tetramers again illustrate that fragment polarizations increase with system size. As observed for the trimers, the middle fragments (Fragments 2 and 3, Table S1) experience greater polarization, while the outer fragments are less affected. However, a comparison of monomer polarization between trimers and tetramers reveals a given fragment is most affected by the presence of a nearest neighbor, and less affected by nextnearest-neighbor fragments. For example, identical intermonomer geometries exist within trimer 3A and tetramer 4A. Within 3A, Fragment 2 is polarized by –1.4 kJ mol1 from two nearest-neighbor fragments (Fragments 1 and 3). In comparison, within 4A, Fragment 2 is also polarized by –1.4 kJ mol1, but in the presence of two nearest-neighbor fragments (Fragments 1 and 3) and a nextnearest-neighbor fragment (Fragment 4). Therefore, a next-nearest-neighbor monomer has little influence on monomer polarization within these systems. The tetramer binding energies differ depending on whether they are calculated with respect to isolated or polarized monomer energies. The differences reach 4.4 kJ mol1, reflecting the increased polarization of fragments within a tetramer. This discrepancy between binding energies calculated with respect to isolated and polarized monomers increases with system size for two reasons: (1) as the number of nearest-neighbor contacts increases, the magnitude of fragment polarization increases; and (2) as the number of polarizable fragments that contribute to the binding energy increases, the deviations can accumulate (Table 3).
The systematic study of the benzene-containing systems detailed above illustrates the effect of fragment polarization on binding energies for systems of increasing size in which the electrostatic components are small. However, similar effects on fragment polarization will also be observed for systems with larger electrostatic components, which will lead to discrepancies in the calculation of binding between either DEFMO or DEsup and DE0FMO . Therefore, in addition to the above systems, a subset of the S66 database was chosen to consider interactions with larger electrostatic contributions [20]. Specifically, interactions in p–p dimers containing: (1) benzene–benzene; (2) benzene–pyridine; (3) benzene–uracil; (4) pyridine–pyridine; (5) pyridine–uracil; and (6) uracil–uracil were evaluated at their equilibrium geometries (Figure 2) [20]. Although this test set also contains a stacked benzene–benzene dimer (Figure 2), its geometry differs from that of 2S (Figure 1). Differences between the electrostatic contribution to the interactions ðDE0ES Þ for the benzene systems and the S66 subset can be seen in the PIEDA results (Figures S2, S3 and Tables S2, S3). As found for the benzene dimers above, the calculated interaction energies of the S66 subset with respect to the isolated monomers (Table 4) show that the same energy is obtained from both FMO2-MP2 and MP2. Furthermore, approximations are not invoked for these small systems since interfragment distances are within the defined cutoffs. However, large discrepancies are found when interaction energies are evaluated with respect to isolated ðDEFMO2 Þ and polarized ðDE0FMO2 Þ monomers. These discrepancies range from 0.0 kJ mol1 for the benzene–benzene dimer, up to 8.4 kJ mol1 for the uracil–uracil dimer. Indeed, monomer polarization is greatest for the uracil homodimer (Figure 2f). These results illustrate that monomer polarization and discrepancies in binding energies will also increase in magnitude as the electrostatic contribution to the interaction increases. 4. Conclusions The present comparative study examines p–p interaction energies of benzene-containing dimers, trimers and tetramers, as well
C.D.M. Churchill / Chemical Physics Letters 554 (2012) 185–189
as a subset of the S66 database using both a supermolecular approach and the fragment molecular orbital method. Results show that both MP2/6-31G(d,p) and FMO-MP2/6-31G(d,p) predict nearly identical absolute energies for the systems considered. Deviations can occur depending on whether the calculation of binding energies is performed with respect to isolated ðDEFMO Þ or polarized ðDE0FMO Þ monomers. Deviations between DE0FMO and DEFMO tend to increase as the number of monomers increases and the electrostatic contribution to the interactions increases. Therefore, caution is suggested in the direct comparison of binding energies reported in the FMO output ðDE0FMO Þ and those calculated with respect to isolated monomers ðDEFMO or DEsup Þ. This is especially true for systems that are large and/or have substantial electrostatic contributions to interaction energy, since the magnitude of fragment polarization is greater in both these instances. Acknowledgments This work was supported by the Natural Sciences and Engineering Research Council (NSERC CGS-D), Alberta Innovates-Technology Futures, and the University of Alberta. I am grateful to Dr. Mariusz Klobukowski for his support and input during the preparation of this manuscript. Thanks are also given to a Reviewer for his or her expertise and suggested improvements to this work. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2012. 09.047.
189
References [1] P. Hobza, Accounts of Chemical Research 45 (2012) 663. [2] D.G. Fedorov, T. Nagata, K. Kitaura, Physical Chemistry Chemical Physics: PCCP 14 (2012) 7562. [3] M.S. Gordon, D.G. Fedorov, S.R. Pruitt, L.V. Slipchenko, Chemical Reviews 112 (2012) 632. [4] D.G. Fedorov, K. Kitaura, J. Chem. Phys. 120 (2004) 6832. [5] D. Fedorov, K. Kitaura, in: D. Fedorov, K. Kitaura (Eds.), The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems, CRC Press, New York, 2009, p. 5. [6] D.G. Fedorov, K. Kitaura (Eds.), The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems, CRC Press, 2009. [7] A. Yoshioka et al., J. Mol. Graphics Modell. 30 (2011) 110. [8] P. Hobza, Annu. Rep. Prog. Chem. C: Phys. Chem. 107 (2011) 148. [9] M.O. Sinnokrot, E.F. Valeev, C.D. Sherrill, Journal of the American Chemical Society 124 (2002) 10887. [10] T.P. Tauer, C.D. Sherrill, Journal of Physical Chemistry A 109 (2005) 10475. [11] M.O. Sinnokrot, C.D. Sherrill, Journal of Physical Chemistry A 110 (2006) 10656. [12] S. Tsuzuki, K. Honda, T. Uchimaru, M. Mikami, K. Tanabe, Journal of the American Chemical Society 124 (2001) 104. [13] D.G. Fedorov, K. Kitaura, Journal of Computational Chemistry 28 (2007) 222. [14] D.G. Fedorov, K. Kitaura, in: E.B. Starikov, J.P. Lewis, S. Tanaka (Eds.), Modern Methods for Theoretical Physical Chemistry of Biopolymers, Elsevier, 2006, p. 3. [15] M.W. Schmidt, K.K. Baldridge, J.A. Boatz, S.T. Elbert, M.S. Gordon, J.H. Jensen, S. Koseki, N. Matsunaga, K.A. Nguyen, S. Su, T.L. Windus, M. Dupuis, J.A. Montgomery, Journal of Computational Chemistry 14 (1993) 1347. [16] M. Suenaga, Journal of Computer Chemistry, Japan 4 (2005) 25. [17] M. Suenaga, Journal of Computer Chemistry, Japan 7 (2008) 33. [18] T. Nakano, T. Kaminuma, T. Sato, K. Fukuzawa, Y. Akiyama, M. Uebayasi, K. Kitaura, Chemical Physics Letters 351 (2002) 475. [19] D.G. Fedorov, K. Ishimura, T. Ishida, K. Kitaura, P. Pulay, S. Nagase, Journal of Computational Chemistry 28 (2007) 1476. [20] J. Rezac, K.E. Riley, P. Hobza, Journal of Chemical Theory and Computation 7 (2011) 2427.