22 October 1999
Chemical Physics Letters 312 Ž1999. 319–324 www.elsevier.nlrlocatercplett
Pair interaction molecular orbital method: an approximate computational method for molecular interactions Kazuo Kitaura a
a,)
, Takuya Sawai a , Toshio Asada a , Tatsuya Nakano b, Masami Uebayasi c
Department of Chemistry, College of Integrated Arts and Sciences, Osaka Prefecture UniÕersity, Sakai-si, Osaka 599-8531, Japan DiÕision of Chem-Bio Informatic, National Institute of Health Sciences, 1-18-1 Kamiyoga, Setagaya-ku, Tokyo 158-8501, Japan c National Institute of Bioscience and Human Technology, Tukuba-si 305, Japan
b
Received 5 July 1999; in final form 19 August 1999
Abstract We propose an approximate method for calculating large molecular clusters using molecular orbital ŽMO. methods. The method allows one to obtain the total energy of a molecular cluster by performing MO calculations on the molecules and the molecular pairs in the system. The computational time is greatly reduced for a system containing a large number of molecules, since the supermolecule calculation of the whole system is avoided. Our method considers important terms of many-body interactions within a pair approximation and is expected to give results of reasonable accuracy even for molecular interactions having significant many-body effects. The numerical calculations of water trimers and tetramers revealed that our method included the major portion of the many-body energies. q 1999 Elsevier Science B.V. All rights reserved.
1. Introduction Theoretical studies on large molecular systems such as proteins and molecular clusters composed of many molecules by the molecular orbital ŽMO. method have received increasing interest as computer power has progressed. On the subject of molecular clusters, a supermolecule calculation has been utilized. The computational time required for a supermolecule calculation, however, increases linearly with the system size, even if the most efficient
) Corresponding author. Fax: q81-0722-54-9931; e-mail:
[email protected]
algorithm is employed w1x and may practically limit the size of molecular clusters to be handled. In this Letter, we present an approximate method for calculating molecular interactions. The computational time for large molecular clusters may be drastically reduced, since our method is a pair approximation and avoids a supermolecule calculation on the whole system. Many-body effects have been known to be important to intermolecular interactions and have been extensively studied in atomic and molecular systems w2x. Among them, the many-body effects in the interactions of water molecules have been paid special attention in relation to the accurate descriptions of the properties of bulk water and the analytical poten-
0009-2614r99r$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 9 . 0 0 9 3 7 - 9
K. Kitaura et al.r Chemical Physics Letters 312 (1999) 319–324
320
tial functions used in the molecular simulations w3,4x. The many-body polarization term has been shown to be the most important in the interactions of water molecules. This term, however, accounts for only about 60% of the total three-body energy of the water trimer w5x, indicating that an accurate model for many-body interactions should take into account other many-body terms along with the many-body polarization term. Our pair interaction molecular orbital ŽPIMO. method accounts for the full many-body polarization term and a portion of other many-body terms. The numerical calculations of water trimers and tetramers revealed that the method well reproduced the geometries and the total interaction energies of the water clusters obtained by the ab initio supermolecule calculations. It indicates that the PIMO method includes significant many-body effects.
Solving the following Schrodinger equations at an ¨ appropriate level of wavefunction such as the Hartree–Fock level, HICI s EIXCI HI JCI J s EIX JCI J one obtains the monomer electronic energy, EIX , and the dimer electronic energy, EIX J . The total energy of monomer, EI , and dimer, EI J , are obtained by adding the corresponding nuclear repulsion energies to EIX and EIX J , respectively. The total energy of N molecular system, E, is calculated using EI and EI J as follows: Es
Ý Ž E I J y EI y E J . q Ý E I , I)J
Es
I
Ý EI J y Ž N y 2 . Ý E I . I)J
2. PIMO method We consider a system composed of N molecules. The monomer Hamiltonian in the system, HI , is defined as follows: nl
HI s HI0 q Ý i
Ý VJesp Ž r i .
Ž 1.
J/I
where HI0 is the Hamiltonian of the isolated molecule I, n I number of electrons in the fragment I and VJesp Ž r . the electrostatic potential of the molecule J: VJesp Ž r . s H d rXr J Ž rX . r< r y rX < y
Ý Zsr< r y rs < sgJ
Ž 2. where r J Ž r . is the electron density of the molecule J, and Z s is the nuclear charge of atom s. The electron density is calculated by an iterative procedure Žvide infra.. Similarly, the dimer Hamiltonian, HI J , is defined as follows: n Iqn J
HI J s HI0J q
Ý
Ý
i
K/I , J
VKesp Ž r i .
Ž 3.
where HI0J is the Hamiltonian of the isolated molecular pair IJ.
Ž 4.
I
Our method gives the total energy of a molecular cluster, so the molecular interaction energy, D E, is calculated in the same manner as the supermolecule calculation: D E s E y Ý EI0 , I
D Es
Ý Ž EI J y EI y EJ . q Ý Ž EI y EI0 . , I)J
Ž 5.
I
where EI0 is the total energy of the free molecule I. Since the monomer energy, EI , is obtained by a self-consistent procedure, the second term of Eq. Ž5. contains the polarization energy up to the Nth order along with the electrostatic energy. The first term of Eq. Ž5. includes other quantum-mechanical interaction terms such as the exchange-repulsion and charge-transfer interaction terms of the dimer under the electrostatic potential from the surrounding Ž N y 2. molecules. The quantum-mechanical terms, hence, are not the same as those of the free dimer, since the exchange-repulsion and charge-transfer interactions are promoted or suppressed under the electrostatic environment. Thus our model includes extra manybody terms over a many-body polarization model and a simple pair interaction model. The many-body effects involved in our model will be demonstrated by numerical calculations on water clusters in Section 3.
K. Kitaura et al.r Chemical Physics Letters 312 (1999) 319–324
321
The computational procedure of the PIMO method is as follows: 1. Carry out the MO calculations of the free monomers to obtain the energies EI0 and the electron density distributions. The electron density distributions are used to construct the initial electrostatic potential in the monomer Hamiltonian ŽEq. Ž1... 2. Construct the monomer Hamiltonians using the given electron density distributions and solve the Schrodinger equations for all monomers in the ¨ system to obtain the monomer energies, EI and the density distributions. 3. Determine whether the electron density distributions of monomers are the same as the previous ones within a specified criterion. If the density has not converged, return to step Ž2. with the new density distributions. 4. Construct the dimer Hamiltonians using the converged density distributions r J Ž r ., and solve the Schrodinger equations for all pairs of molecules ¨ in the system to obtain EI J . 5. Calculate the total energy using EI and EI J ŽEq. Ž4.. and other properties such as electron populations.
3. Results and discussions It is known that many-body effects are significant in clusters of water molecules. Water clusters, therefore, may be appropriate systems to examine our PIMO method. Geometry optimizations were carried out on water trimers and tetramers at the HFr631GUU level. In all calculations the geometry of ˚ and water molecule was kept rigid ŽROH s 0.9572 A /HOHs 104.528 w6x. and the intermolecular degrees of freedom w ere optim ized. The GAUSSIAN 94 program w7x was used for the supermolecule calculations. The optimized geometries of water trimers are shown in Fig. 1. The geometries from the PIMO calculations are in good agreement with those from the ab initio supermolecule calculations for all water trimers. The largest error is found in 3-C; the O–O length and the O–O–O angle are overestimated by ˚ and 2.28, respectively. The interaction ener0.01 A gies of the trimers, D E, are given in Table 1. The
Fig. 1. The optimized geometries of water trimers. The lengths are ˚ and angles in degrees. The values from the ab initio in A supermolecule calculations are given in parentheses.
PIMO interaction energies agree well with the corresponding ab initio energies within an error of 0.1 kcalrmol. The interaction energy is expanded as sum of n-body energies, V Ž n. : D E s V Ž2. q V Ž3. q . . . qV Ž n. , V
Ž2.
s
ÝŽ
EI0J y EI0 y E J0
Ž 6.
.,
I)J
V Ž3. s
Ý Ž EI0JK y EI0 y E J0 y EK0 . y V Ž2. , I)J)K
where EI0J and EI0JK are the total energies of the free dimer IJ and the free trimer IJK, respectively. In
K. Kitaura et al.r Chemical Physics Letters 312 (1999) 319–324
322
Table 1 Energies of water trimers from PIMO calculationsa c
E D Ed V Ž3.e
3-Ab
3-B b
3-C b
y228.09747 Žy228.09716. y17.5 Žy17.4. y2.1 Žy1.9.
y228.08568 Žy228.08558. y10.2 Žy10.1. q0.3 Žq0.3.
Žy228.08156. Žy228.08157. y7.6 Žy7.6. q0.5 Žq0.6.
a
The energies from the ab initio MO supermolecule calculations are given in parentheses. See Fig. 1. c Total energy in a.u. d Interaction energy in kcalrmol. e Three-body energy in kcalrmol. b
order to show the many-body energies involved in the PIMO results, the three-body energies of the trimers, V Ž3., were calculated and shown in Table 1 along with those from the ab initio supermolecule calculations. The errors in the three-body energies are small for the trimers, proving that the PIMO method takes account of a major portion of threebody effects.
The polarization term of the many-body interaction energy has been known to be the most important for polar molecules w2x. From the formulation of the PIMO method, one can easily understand that the method includes the full many-body polarization term. The three-body polarization term, however, contributes about 45% of the total three-body energy in 3-A that has the largest three-body energy among
Fig. 2. The optimized geometries of water tetramers. See the caption of Fig. 1.
K. Kitaura et al.r Chemical Physics Letters 312 (1999) 319–324
323
Table 2 Energies of water tetramers from PIMO calculationsa a
E D Ea V Ž3.a V Ž4.c
4-Ab
4-B b
4-C b
4-D b
y304.14035 Žy304.13965. y29.9 Žy29.6. y5.5 Žy4.9. y0.4 Žy0.4.
y304.12987 Žy304.12944. y23.3 Žy23.1. y2.4 Žy2.2. y0.0 Žy0.2.
y304.12054 Žy304.12003. y17.5 Žy17.2. q1.2 Žq1.6. y0.1 Žy0.1.
y304.13310 Žy304.13190. y25.4 Žy24.6. y2.8 Žy2.0. q0.1 Žq0.1.
a
See Table 1. See Fig. 2. c Four-body energy in kcalrmol. b
the three trimers; the three-body polarization energy is estimated to be y0.9 kcalrmol from the energy decomposition analysis w8x at the HFr6-31GUU level. This suggests that the PIMO method includes not only the three-body polarization term but also other three-body terms within a pair approximation. A detailed analysis to reveal the many-body terms involved in the PIMO method will be reported elsewhere. The optimized geometries and the interaction energies of several water tetramers are shown in Fig. 2 and Table 2, respectively. The PIMO geometries are in almost perfect agreement with the ab initio ones; ˚ in the largest error in the O–O distance is 0.012 A 4-A. The PIMO interaction energies are also in good agreement with the ab initio ones. The errors in the interaction energy are q0.3, q0.2, q0.3 and q0.8 kcalrmol for 4-A, 4-B, 4-C and 4-D, respectively, which are of the same order of magnitude as those of the three-body energies. The results show that the PIMO method has a tendency to overestimate negative three-body energies and underestimate positive three-body energies. The PIMO method also gives reasonable four-body energies, V Ž4., though they are very small in magnitude.
major portion of many-body interactions into consideration within a pair approximation and to give reliable results for molecular clusters, as far as polar molecules are concerned. We have applied the method to other hydrogen bonding clusters and several hydrated ion clusters, in which many-body effects are known to be significant, and obtained the geometries and energies in reasonable accuracy. These works will be published elsewhere. The PIMO method enables one to calculate the total energy of a many molecular system without performing the supermolecule calculation of whole system. The PIMO method is, therefore, expected to reduce the computational time of large molecular clusters drastically. For practical applications, however, the use of the fast multipole expansion method w1x for Coulomb type two-electron integrals involved in the electrostatic potential terms ŽEq. Ž2.. may be critical for saving computational time. Another advantage of the method is its easiness for utilizing parallel processing, since the monomer terms and the dimer terms can be calculated independently. The PIMO method may provide not only a useful computational method but also a simple model to understand many-body effects in molecular interactions.
4. Summary A new pair approximation model has been proposed for calculating intermolecular interactions using the molecular orbital method: the PIMO method. The PIMO method has been successfully applied to calculations of geometries and the energies of water clusters which have significant many-body effects. Thus, the PIMO method has been proved to take a
Acknowledgements This work was supported in part by grant from the Ministry of Education ŽGrant No. 08221227.. The numerical calculations were partially carried out at the Computer Center of the Institute for Molecular Science.
324
K. Kitaura et al.r Chemical Physics Letters 312 (1999) 319–324
References w1x C.A. White, B.G. Johnson, P.M.W. Gill, M. Head-Gordon, Chem. Phys. Lett. 230 Ž1994. 8. w2x M.J. Rlrod, R.J. Saykally, Chem. Rev. 94 Ž1994. 1975. w3x J. Detrich, G. Corongiu, E. Clementi, Chem. Phys. Lett. 112 Ž1984. 426. w4x T.P. Lybrand, P.A. Kollman, J. Chem. Phys. 83 Ž1985. 2923. w5x G. Chalasinski, M.M. Szczesniak, P. Cieplak, S. Scheiner, J. Chem. Phys. 94 Ž1991. 2873. w6x W.S. Benedict, N. Gailar, E.K. Plyler, J. Chem. Phys. 24 Ž1956. 1139.
w7x GAUSSIAN 94, Revision D.4, M.J. Frisch, G.W. Trucks, H.B. Schlegel, P.M.W. Gill, B.G. Johnson, M.A. Robb, J.R. Cheeseman, T. Keith, G.A. Petersson, J.A. Montgomery, K. Raghavachari, M.A. Al-Laham, V.G. Zakrzewski, J.V. Ortiz, J.B. Foresman, J. Cioslowski, B.B. Stefanov, A. Nanayakkara, M. Challacombe, C.Y. Peng, P.Y. Ayala, W. Chen, M.W. Wong, J.L. Andres, E.S. Replogle, R. Gomperts, R.L. Martin, D.J. Fox, J.S. Binkley, D.J. Defrees, J. Baker, J.P. Stewart, M. Head-Gordon, C. Gonzalez, J.A. Pople, Gaussian, Pittsburgh, PA, 1995. w8x K. Kitaura, K. Morokuma, Int. J. Quantum Chem. 10 Ž1976. 325.