Journal of Molecular Structure (Theochem), 122 (1985) 343-350 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
MNDO AND MNDO/H CALCULATIONS ON HYDROGEN BONDS -A COMPARISON WITH AB INITIO AND CNDO/B METHODS
J. KOLLER,
V. HARB, M. HODO%EK
Boris KidriE Znstitu te of Chemistry,
and D. HAD21
61000 Ljubljana,
Hajdrihova
19 (Yugoslavia)
(Received 13 December 1984)
ABSTRACT Semiempirical SCF MNDO and its MNDO/H modification were tested on a series of hydrogen bonded systems. H-bond lengths, H-bond energies, proton potential functions and charge transfer were calculated and compared with ab initio and CNDO/P results. It is shown that the MNDO method fails completely. Although the results obtained with the MNDO/H modification are, in general, superior, disturbing deviations from the ab initio results were obtained in some cases and thus MNDO/H cannot be relied upon. INTRODUCTION
There is still considerable scope for the application of semiempirical SCF MO techniques to hydrogen bonded systems that are too large for treatment by ab initio methods. These are particularly the molecules of biological importance and more generally, problems in which extensive variation of geometrical arrangements have to be treated and where only semiquantitative results are required. Despite its well known deficiencies the CNDO/B method is still being used in such cases. Although the theoretically better justified MINDO/3 and, even more so, MNDO techniques appeared to be successful in calculating the heats of formation, geometries, ionization potentials and dipole moments [l, 21 both of them failed when applied to simple hydrogen bonded test cases such as water, HCN dimers, etc. [ 4-91. The results were inferior even to those yielded by CNDO/B. The reason for failure of the methods appeared to be in the inadequate parameterization for intermolecular interactions rather than in the theoretical approximations. This prompted us to attempt a reparametrization of only those integrals where the atoms taking part in the H-bond, were involved. The limited reparameterization gave no satisfactory results. At this time a paper by Burstein and Isaev appeared, [3] describing a different approach. Only the core-core interaction function for the interaction between H and Y in the H-bond XH * * * Y have been modified by these authors and one empirical parameter has been implemented. The proposed interaction function is f(RYH) = exp (-c? YH), OLbeing the new empirical parameter. The best value for (Y reported is 2.0 A” for all pairs H * * * Y, Y = N: 0, F 0166-1280/85/$03.30
0 1985 Elsevier Science Publishers B.V.
344
The tests carried out on a set of small H-bonded complexes appeared promis_ ing and the authors claim the modified version MNDO/H to be at least equivalent in calculating the properties of H-bonded systems to the ab initio method with the 4-31G basis set. Since the test examples were limited both in number and the types of _H-bonding it appeared wise to extend it and include intramolecular H-bonds and, particularly, the calculation of proton potential functions. These are known to be very sensitive to the quality of calculation. On the other hand, situations with the proton transfer are often met in biological systems, e.g., enzyme mechanisms. The treatment of such systems by ab initio methods is usually limited to small basis sets and/or radical trimming of the number of atoms taken into consideration. A semiempirical method would be very useful for such problems, but its reliability has to be adequately tested. We have calculated by the MNDO/H method the distances R(X **. Y), H-bond energies and Mulliken populations for the following complexes: NH,-HOF, (HF),, HF-HOH, H,O-HF, (HCOOH)2, CO,-HF, (H,O), linear, (H,O), cyclic, CH3NHz-HCOOH, CH~NI&-COOH and benzoic acid dimer. As an example of proton transfer systems the benzoic acid dimer and the methylamine-formic acid complex were investigated. For a test on intramolecular hydrogen bonding the hydrogen maleate ion was treated. Since the calculation of the conformational energies of molecules with sp’ hybridized atoms is known to be a weak side of the MNDO techniques, we included in this work also malonic acid. For all these systems comparisons are made with results obtained by the original parameterization and with ab initio calculated quantities. CNDO/B calculations were also made for the purpose of comparison. COMPUTATIONAL
The computer program for the MNDO method with the standard parameterization was used. This version is hereafter named original MNDO. Modifications proposed by Burst&n and Isaev [3] were implemented to obtain another version denoted MNDO/H. The gradient method using the DFP (Davidon-Fletcher-Powel) algorithm [ll] for the geometry optimization was used with both methods. For the CNDO calculations standard CNDO/2 parameterization was used. Mulliken populations, H-bond energies and H-bond lengths were compared with the corresponding ones obtained by the ab initio technique. When possible the data were obtained from literature (preferably calculations with at least double zeta quality basis sets) otherwise the calculations were done with the split shell 4-31G basis. RESULTS
Figure 1 shows the interaction energy as the function of OH distance for the linear HBN*** HOF system. Two local minima were found by MNDO:
345 I
AM 6
5' b' 3' 2' 1' 0 -1
Fig. 1. Interaction energy as the function of OH distance for H,N . a*HOF calculated by are 4.04, 2.64 and 2.62 A for MNDO, MNDO/H MNDO, MNDO/H and CNDO/%. R,.... and CNDO/S, respectively.
R(OH) = 0.96 A, and R(OH) = 3.0 A, the latter being less deep and corresponds to the situation described as the ionic state. No ionic states are usually detected by the semiempirical methods. Neither CNDO/B produced the ionic state although the equilibrium bond lengths and bond energies were reasonable. MNDO/H and CNDO/B curves are very similar. The largest system treated was the benzoic acid dimer with 92 valence electrons out of the total of 128 electrons. Figure 2 presents the potential energy functions for this system. All the geometric parameters except for the geometry of the phenyl ring were optimized during the calculation of each point. It has been shown by Nagaoka and coworkers [17, see also refs. 18-211 that extensive optimization is needed in order to obtain reliable results. Because of the possibility of direct comparison with the
-1.01’ 0.0
’ ’ 1.0
’ ’ 1.2
’ ’ 1.b
’
1.6
’
yir
Fig. 2. Proton potential energy functions for benzoic acid dimer (double proton transfer) by MNDO and MNDO/H methods.
346
results obtained by these authors the Cih symmetry was retained during the calculations (double proton transfer). The proton potential energy functions have their minima at the O-H distance of 0.95 A and 0.99 R for the MNDO and MNDO/H scheme, respectively. For the transition state we report only the O* * 0 distances: 2.38 A (MNDO) and 2.49 A (MNDO/H). No comparison can be made in this case because the corresponding figures cannot be extracted from the data in ref. 17. The potential energy barrier heights are 64.61 and 27.86 kcal mol-’ for MNDO and MNDO/H, respectively. These values are too high compared [17] to the ab initio STO-3G minimal set value (4.4 kcal mol-‘) or to the 4-31G extended basis set value (12.7 kcal mol-‘). The interest in malonic acid lies in the reproduction of the overall geometry which is strongly influenced by the non bonded interactions, and in the weak intramolecular hydrogen bonds in the free carboxylic group. Both MNDO and MNDO/H overestimated the 0 * * - 0 distance and also the COH angle in comparison with the ab initio calculation. However, the overall geometry was reproduced quite well. This concerns particular the tilt of both carboxylic groups, the dihedral angles for MNDO/H turned out to be 87.7” and 12.3” (ab initio [28] : 90”; 0”; crystal structure [29] : 90”; 13”). Table 1 presents the H-bond energies for the considered systems. Coniplexes with intramolecular H-bonds are not included because of the well known difficulties with the reference structures. It can be seen that MNDO underestimates the H-bond energies by roughly one order of magnitude. Moreover, it even predicts, in one case, that the structure is unstable. l
TABLE 1 H-bond energies E,, calculated by MNDO and MNDO/H compared with the corresponding ab initio values (energies in kcal mol-’ ) Systems NH,-HOF (HF)2 HF-HOH H,O-HF ( HCOOH,)a CO,-HF (H,O), linear (H,O), cyclic CH,NH,-HCOOH (neutral form) CH,NH,-COOH (ionic form) Benzoic acida Dimer
E,(MNDO) -0.32 -1.25 -0.45 -1.39 -0.47 -0.71 -0.595 -0.465 + 2.23 -102.21 -1.11
E, (MNDO/H) -10.50 -8.97 -4.53 -9.375 -5.25 -5.97 -4.37 -6.07 -131.63 -13.647
E, (ab initio) -9.80b
-3.35 --5.05c -2.66d -9.29d -7.9---8.1e -2.38--4.95f -3.7--12.6g -4.od -14.47 129.59 -10.6
aEnergy per hydrogen bond. bReference 10. CReference 22, p. 65, and references therein; refs. 4, 23. dReference 4. eReference 22, p. 117, and references therein. fReferences 25, 26. gReference 22, p. 62, and references therein; ref. 4.
341
H-bond energies are better reproduced by MNDO/H and they are all of the same order of magnitude as the ab initio calculated ones. A short remark is necessary for the ionic form of methylamine-formic acid complex: the geometry of the complex is not at the energy minimum. The NH distance has to be kept fixed lest the proton be not pulled away by the oxygen. It has been shown by Kollman and Allen [ 131 that the MNDO method predicts equal stability for both the linear and cyclic water dimers. According to the MNDO approximation the linear conformation is preferred by 0.13 kcal mol.‘. Unfortunately, the bifurcated H-bond of the cyclic water dimer can not be fitted into the MNDO/H scheme and this appeared to be a major defect of the method. Optimized X **- Y distances are given in Table 2. All the internal coordinates except for the phenyl rings were optimized. As expected MNDO overestimates the intermolecular X *- - Y distances by a factor of 1.5 to 2.0. The trend of changing these distances by MNDO/H is correct. For the hydrogen maleate ion (intramolecular H-bond) the MNDO 0 - - - 0 distance is 2.36 R which happens to be equal to the distance predicted by the ab initio calculation with the large 6-311G basis set [16]. The planar geometry of the ion has been assumed. A symmetric proton potential function with single minimum is predicted by both the MNDO and MNDO/H methods, the proton being pushed slightly (0.15 R) out of the 0 - - - 0 line. MNDO/H enlarges the 0 - - - 0 distance by 0.1 A. Split valence basis (4-31G and 6-311G) yielded double minimum potential functions with a central barrier. TABLE 2 Hydrogen bond lengths R,. ..y calculated (distances in A ). System
E,
NH,HOF (HF), HF-HOH H,O-HF (HCOOH), CO,-HF (H,O), linear (H,O), cyclic CH,NH,-HCOOH (neutral form) CH,NH,-COOH (ionic form) Hydrogen maleate ion Benzoic acid (dimer) Malonic acid (monomer)
on MNDO, MNDO/H and ab initio level
R,. ..Y(MNDO/H)
R,.
4.04 3.93 3.89 3.92 4.43 4.14 4.225 4.22 3.98
2.53 2.55 2.62 2.53 2.55 2.60 2.615 2.525
2.77a 2.69-2.91b 2.94-3.15C*d 2.71c 2.73e 2.86f 2.73-3.01s 2.90c 2.72
2.40
2.39
2.55
2.36 3.78 2.23
2.46 2.51 2.21
2.36n 2.72’ 2.3+-2.42’
....(MNDO)
. . Y(ab initio)
*Reference 10. bReference 4, 22, 23. CReference 4. dReference 24. eReference 22 and references therein. fReference 25, 26. sReference 22 and references therein; ref. 4. hReference 16. ‘Reference 17. iReference 28.
348
The height of the barrier depends on the level of computation [16]. Lack of definitness is notable in the MNDO/H scheme when an assymetric H-bond is under consideration and the proton is moved along the bond: where to replace the initial core-core interaction by the new one? We tried to overcome this by changing both interaction functions: f(Rv,) as well as P(RXH). For H3N--- HOF only a slight enlargement of the N - -a0 distance is observed when compared with the calculation with only one core-core interaction function. It changes from 2.53 A (Table 2) to 2.59 A. The potential energy tends to have the minimum near the middle of the bond (Fig. 3). Further, we consider the ionization potentials. The differences between the first ionization potential of water monomer and dimer is compared with the nonempirical large basis set calculation of Diercksen [ 141. For the linear conformation the first ionization potential changes from 12.19 eV (monomer) to 11.90 eV (dimer, MNDO) and 11.81 eV (dirner, MNDO/H). The difference is 0.29 and 0.38 eV for MNDO and MNDO/H, respectively, compared to 0.78 eV by Diercksen (monomer: I. = 13.84 eV, dimer: IP = 13.06eV). Photoelectron spectroscopy data are also available concerning the two lowest ionization potentials of the water dimer [ 151. Table 3 provides the comparison between our MNDO and MNDO/H calculations, the ab initio calculations without and with configuration interaction and experimental data (position of the first two bands in the photoelectron spectrum). The agreement between our calculations and the experiment happens to be good. To demonstrate the charge rearrangement provided by the MNDO methods obtained by the Mulliken population analysis the linear water dirner was chosen. In Table 4 the MNDO and MNDO/H atomic charge
I
0.9
I
I
1.1
I
1
1.3
8
I
1.5
1
L
IQ1
Fig. 3. Proton potential energy curve for H,N . * * HOF by MNDO/H (RN.. new core--core interactionfunctions on both sides.
.
0
= 2.59 A ),
349
TABLE 3 The first two ionization potentials for the water dimers: comparison culated and experimental values (all values in eV).
between the cal-
SCF-MOa*b
CIb
MNDO (linear)
MNDO/H (linear)
MNDO (cyclic)
Experimentb
12.74 14.08
11.07 12.44
11.90 12.37
11.81 12.32
12.14’ 12.14
12.1 i: 0.1 13.2 * 0.2
‘6-31G basis set. b Reference 12. cThe two highest filled molecular orbitais are degenerate and of n symmetry. TABLE 4 Comparison between atomic charge density changes when H-bond linear form of water dimer (H,H,O, . . . H,O,H,). H 1.2 MNDO MNDO/H STO-3G CND0/2
-0.0172 +0.0147 +0.018 + 0.017
0, -I-0.0360 -0.0087 -0.007 + 0.002
H, -0.0013 +0.0315 + 0.034 + 0.042
0, + 0.0265 -0.0396 -0.040 -0.052
is formed for the
H6 -0.0259 -0.0129 -0.024 -0.026
Charge transfer 0.0016 0.0207 0.029 0.036
density changes on H-bond formation are compared by the corresponding quantities obtained by the semiempirical CNDO/B method and by the ab initio calculations with the minimal STO-3G basis set [22]. The overall charge transfer has the correct sign for both methods tested (charge transfered from proton acceptor to proton donor), the MNDO/H value being also of the same order of magnitude as those with the CND0/2 and ab initio calculations. Considering the individual changes of charge it becomes immediately obvious that MNDO fails completely. At best it reproduces the correct sign. MNDO/H charge differences are surprisingly good and are definitely better than those calculated by the CND0/2 method CONCLUSIONS
The MNDO semiempirical method with the original parameterization does not yield acceptable results for H-bonded systems. It has similar deficiencies to MIND0/3: the inadequate description of core-core repulsion, the incorrect p-S proportionality for interactions between neighbours (0,” = -k,, S,,) and it neglects all important higher order corrections for nonneighbour p’s [27]. In some cases the compensation of unwarranted errors takes place and the calculation is more accurate. The MNDO/H version seems to correct the first error from the above list but leaves all the others
350
untouched. This is the reason why the results of MNDO as well as of MNDO/H calculations are unreliable, although the latter modification yields superior results in most cases treated in this work. REFERENCES 1 2 3 4 5 6
M. J. S. Dewar and W. Thiel, J. Am. Chem. Sot., 99 (1977) 4899. M. J. S. Dewar and W. Thiel, J. Am. Chem. Sot., 99 (1977) 4907. K. Y. Burstein and A. N. Isaev, Theor. Chim. Acta, 64 (1984) 397. P. A. Kollman and L. C. Allen, J. Am. Chem. Sot., 92 (1970) 753. T. J. Zielinski, D. L. Breen and R. Rein, J. Am. Chem. Sot., 100 (1978) 6266. G. Klopman, P. Andreozzi, A. J. Hopfinger, 0. Kikuchi, and M. J. S. Dewar, J. Am. Chem. Sot., 100 (1978) 6267. 7 M. J. S. Dewar and G. P. Ford, J. Am. Chem. Sot., 101 (1979) 5558. 8 S. Scheiner, Theor. Chim. Acta, 57 (1980) 71. 9 V. Harb, J. Koller and D. HadZi, unpublished results. 10 J. E. DelBene, J. Am. Chem. Sot., 95 (1973) 5460. 11 R. Fletcher, Comput. J., 13 (1970) 317. 12 J. S. BinkIey, J. A. Pople, W. J. Hehre, J. Am. Chem. Sot., 102 (1980) 939. 13 P. A. Kollman and L. C. Allen, Chem. Rev., 72 (1972) 283. 14 G. H. F. Diercksen, Theor. Chim. Acta, 21 (1971) 335. 15 S. Tomoda, Y. Achiba and K. Kimura, Chem. Phys Lett., 87 (1982) 197. 16 M. Hodowek, F. Avbelj and D. HadZi, to be published. 17 S. Nagaoka, N. Hirota, T. Matsushita and K. Nishimoto, Chem. Phys. Lett., 92 (1982) 498. 18 B. H. Meier, R. Mayer, R. R. Ernst, P. Zolliker, A. Furrer and W. Halg, Chem. Phys. Lett., 103 (1983) 169. 19 K. Fur%, Chem. Phys. Lett., 108 (1984) 518. 20 B. H. Meier, R. Mayer, R. R. Ernst, A. Stockli, A. Furrer, W. HaIg and I. Anderson, Chem. Phys. Lett., 108 (1984) 522. 21 S. Nagaoka, T. Terao, F. Imashiro, N. Hirota and S. Hayashi, Chem. Phys. Lett., 108 (1984) 524. 22 P. Schuster, Energy surfaces for hydrogen bonded systems, in P. Schuster, G. Zundel and C. Sandorfy (Eds.), The Hydrogen Bond/I, North Holland, 1976. 23 G. Alagona, E. Scrocco and J. Tomasi, Theor. Chim. Acta, 47 (1978) 133. 24 P. Kollman, J. McKelvey, A. Johanson and S. Rothenberg, J. Am. Chem. Sot., 97 (1975) 955. 25 A.-M. Sapse, J. Am. Chem. Sot., 78 (1983) 5733. 26 A.-M. Sapse, J. Am. Chem. Sot., 78 (1983) 5738. 27 S. DeBruijn, Int. J. Quantum Chem., 25 (1984) 367. 28 M. Merchan, F. Tomas and I. Nebot-Gil, J. Mol. Struct., Theochem, 109 (1984) 51. 29 J. D. Goedkoop and H. C. MacGillavry, Acta Crystallogr., 10 (1957) 125.