- Email: [email protected]
Comparison of ab initio and multipole determinations of the electrostatic interaction of acetamide dimers
Journal of Molecular Structure (Theochem), 282 (1993) 19-31 0166-1280/93/%06.00 0 1993 - Elsevier Science Publishers B.V., Amsterdam
19
Comparison of ab initio and multipole determinations of the electrostatic interaction of acetamide dimers Julia C. White, Ernest R. Davidson* Department
of Chemistry,
Indiana
University, Bloomington,
IN 47405,
USA
(Received 1 November 1991)
Abstract In this study, we used a constrained least-squares method to determine parameters for use in an atomic multipole description of the electrostatic energy. To test the accuracy of the parameters, we compared acetamide dimer binding energies determined from the multipole scheme with those from a Morokuma analysis. Interestingly, it was found that a model employing only atomic point charges was able to accurately reproduce the ab initio electrostatic energies at the van der Waals separation of the dimer. The least-squares charges were closer to Mulliken charges than to other definitions which have been suggested. At closer distances, penetration of the charge distributions invalidates the multipole expansion.
Introduction Because the size of many compounds exceeds our ability to perform sophisticated ab initio calculations, considerable effort is being put into the development and improvement of semiempirical formulae and parameters [l]; MMZ [2], AMBER and CHARMM [4] are among the most commonly used semiempirical programs. The quality and usefulness of each is determined by the force field it employs. For example, MM~ and its related programs include more sophisticated algorithms for the description of bond stretching and bending, torsional motion, and so forth, than AMBER or CHARMM. Thus MM~ is more effective in treating medium-sized systems, whereas the strength of AMBER and CHARMM comes from their ability to examine much larger systems. Goddard and coworkers have recently proposed a force field with * Corresponding
author.
parameters defined by atomic bonding or hybridization (e.g. sp’, sp3) [5]. In this paper, we deal with the problem of modeling the electrostatic part of the interaction between hydrogen bonded dimers. For this term, Weiner and Kollman [3] and Karplus and coworkers [4] use a chargesharge formula qAqB/ ERAB9 where E is the dielectric constant of the medium, qA and qB are atom centered point charges and RAB is the distance between the atoms. Allinger [2] employs bond dipoles rather than charges on atoms. All three programs treat hydrogen bonding as a separate type of interaction, often modeling it with a Lennard-Jones (10,12)-potential. Method We chose acetamide (1) as a simple compound, capable of forming hydrogen bonds like those in complex biological molecules and examined the electrostatic interaction within the dimer.
20
J.C.
White and E.R. DavicLFonlJ. Mol. Struct. (Theochem)
H3 0
Hs
,‘
/ Cl
H4
c2
N
Hl
>i-
Hi
Basis set
Restrictions on computer time and space forced us to limit the size of the basis set, and experience with standard basis sets (e.g. 6-31G*) [6] led us to conclude that we should create our own in order to describe the hydrogen bond more accurately. Van Duijneveldt [7] hydrogen (10s) and Partridge [8] carbon, nitrogen, and oxygen (13s,8p) primitives were used to calculate the corresponding atomic wavefunctions. The three most diffuse Partridge s primitives were modified so that the ratio of exponents L/L+ I was 2.3 and a fourth function was added. The same thing was done to the four most diffuse Partridge p primitives, and a fifth was added. In this way, the initial (13s,Sp) set was increased to (14s,9p). Atomic orbital coefficients were used in contractions to yield a triple-zeta valence basis set [4s3p] for each heavy atom. The 1s orbital of the neutral atom formed the 1s basis function. The remaining three s and p basis functions were obtained from the 2s and 2p atomic orbitals of the cationic, anionic, and neutral wavefunctions respectively. The cationic and anionic self-consistent field (SCF) calculations were done with frozen 1s cores from the neutral atom. The [2s] basis set for hydrogen was obtained from the 0 orbital of molecular hydrogen and the 1s orbital of atomic hydrogen. Polarization functions were added to give a final basis set of [4s3pld/2slp]. The energy
282 (1993)
19-3I
lost in the contractions is illustrated in Table 1. The SCF energies of the cation, anion, and neutral states calculated with the [4s,3p] basis agree well with the uncontracted basis set values (within 0.03 kcal mall’), and fairly well with the numerical Hartree-Fock values (within 0.50 kcal mall’ ). Additional states were examined for comparison (e.g. ‘S of neutral carbon). Another motivation for developing our own basis set was the reduction of basis set superposition error, an artificial lowering of the energy because the dimer has more basis functions to utilize than either monomer alone. To check the quality of the final basis set, the SCF energy of the monomer in its own basis (- 208.053785 a.u.) was compared with the energy of the monomer in the full dimer basis (- 208.05385 a.u.) using a doubly hydrogen bonded dimer configuration at an N * * *0 distance of 3.6A. We estimate that the basis set superposition error will be less than 0.1 kcal mall’ for any dimer configuration. The geometry of acetamide was taken from the crystal structure reported by Kitano and Kuchitsu [91. Linear least-squares jit
Potentials were calculated at about 1600 random points around the acetamide molecule using MELDF. The electron charge density for molecule d is defined, in the SCF approximation, as N/2 P.&) = 2 1 y?(r) i=I
(1)
where ‘Pi is an occupied molecular orbital. The electrostatic potential at the point rol due to the charge distribution of molecule d is IL&,) = - 1 P_&)/lr - r,ldr + 1 z,JJL
(2)
where R,, is the distance between point rrr and nucleus A. Buckingham [lo] proposed that the potential at some point in space outside the charge distribution
J.C. White and E.R. Davidson/J. Mol. Strut.
(Theochem)
282 (1993)
19-31
21
Table 1 Atomic SCF energies” [4%3Pl
MPb
Carbon ‘P(C)
- 37.688539
- 37.688557
2P(c+) %(C - ) ‘S(C) ‘D(C)
-
- 37.292166 - 37.708539
Nitrogen %N)
- 54.400802
- 54.40083 1
‘P(N+) ‘P(N-) *D(N) 2P(N)
-
53.887870 54.321322 54.295927 54.227684
- 53.887904 - 54.321345
- 53.887896 - 54.321344
Oxygen 3w)
- 74.809185
%(o+) 2P(O-) ‘S(O) ‘D(O)
-
- 74.809232 - 74.372442 - 74.789030
- 74.372440 - 74.789029
37.292141 37.708527 37.549110 37.631166
74.372392 74.788991 74.610601 73.727459
Frozen corec
Numerical Hartree-Fock
- 37.688619
- 37.292154 - 37.708539
- 37.292224 - 37.708844 - 37.549611 -37.631331 -
54.400934 53.888005 54.321959 54.296169 54.228102
- 74.809398 - 74.372606 - 74.789746 -74.611021 - 74.729264
Hydrogen
H H*
- 0.499998 - 1.128528 4
- 0.499999 - 1.128528
-0.5 - 1.13363d
“All energies are in atomic units. .,_’ b Modified Partridge primitive set, as described in text. ‘Calculations performed with the modified Partridge basis, freezing the 1s core. ‘At 1.4 a,,.
of a set of atoms can be approximated
as
certain matrix elements, (+)oA
+
1 1 &~a m n
%n~ &~a
I&
S=
’’
‘1
=
&ii
C$‘tfl
(3)
where qA, ,uA , and 6Aare the atom centered charges, dipoles, and quadrupoles, respectively, and m, n take the values x, y, z. Using a method previously explored by other authors [l 1,121 we performed a linear least-squares fit of the electrostatic potential at selected points about acetamide to a set of atomic properties. To ensure correct long range behavior of the potential, the set of atomic properties was constrained to reproduce certain moments of the SCF charge density. As an example, consider the simplest case using only atomic charges subject to the constraint of zero total charge. After defining
(V,
@A
=
=
t&,)
qA
and
ce)A
=
1
the least-squares error with a Lagrange multiplier for the total charge constraint can be written as E = (44 - VY(f#q - V) + Iq’e
(4)
We can write the equation for optimum charges derived from this as q = Sm’##2 v-
(1/2)IS-‘e
(5)
so that the final total charge Q is Q = eTq = (t~~S-‘f$~ V) - (1/2)A(eTSP’e)
(6)
22
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
which can be used to find 3, by setting Q to zero. These equations are easily generalized to include atomic dipoles and quadrupoles. The quality of the fit can be measured by the R2 and standard deviation (SD) values. The value of R2 is given by (7)
culate the electrostatic interaction energy of acetamide dimers in a variety of conformations with charge densities frozen at the free molecule SCF value. This electrostatic energy (ES), to which the multipole series converges for large intermolecular separations, is given by ES = 1 G K&)
- 1 p&)Y&)dz
(10)
B
where
where d and $? label the interacting molecules. By expanding each term in a Taylor series in Rik,the ES values may be rewritten in terms of atomic multipoles as
P= N,_’ CK CI
D,=v,-P
i,=N,-‘~+, K* = AA
16A
(11)
and the variance o2 is 2 = (44 - VT’&
- IY(N, - Nr)
(8)
Similarly, the variance in the value of an individual parameter qA can be estimated as a; = (S-‘),$
(9)
In Eq. (8), N, is the number of points at which the potential is sampled and Nr is the number of degrees of freedom in the calculation (the total number of atomic properties minus the constraints). For example, the charge-only model of a&amide with molecular charge and dipole constraints would use an Nr value of 5. One important factor that has been neglected so far is that the sum in the definition of the matrix S will grow like Ni’3 with most of the weight on distant points, if the points CIare spread uniformly over all space. We tried incorporating a distancedependent weight eeLR in the least-squares fit program in order to reduce the contribution of points that were far away. For the selection of points we used, this made no appreciable difference. Multiple determination of electrostatic energy Once the atomic properties were generated from the linear least-squares fit, they were used to cal-
To test the accuracy of the ES values calculated for the various dimer geometries with the atomic multipole expansion, we compared these values (Eq. (11)) with the corresponding ones obtained from ab initio calculations (Eq. (10)). Numerous methods have been proposed for decomposing the interaction energy between non-covalently bonded species [13]. The method of Umeyama and Morokuma [ 141has already been incorporated into the program MELDF, so it was a simple matter to generate a list of Morokuma ES values for a randomly generated set of configurations. These values for ES are identical to those defined by Eq. (10). Results and Discussion Potentials at points distributed about acetamide were calculated using the [4s3pld/2slp] basis set and were symmetric about the mirror plane of the molecule. This set includes points just at the molecular van der Waals radius and up to 15 A away. Atomic charges and dipoles were generated from the set of 1597 points, with the number of constants, NCNST, equal to 4 in order to exactly reproduce the SCF molecular charge and dipole and with L = 0 to set the distance dependent weight to unity. To test the accuracy of the fit, the
J.C. White and E.R. DaviakonlJ. Mol. Struct. (Theochem) 282 (1993) 19-31 30 *
15
10
5
0
-5
10’ -100
-50 Potentials
0
50
100
(159 7) in koal/mol
Fig. 1. Error in recalculated electrostatic potential about acetamide versus the potential for the 1597 point set. Potentials were recalculated with the [I,41 atomic properties (see text).
potentials were recalculated from the atomic properties and a plot made of the magnitude of the error (Fig. 1). (The [m. n] notation is defined in Table 2. In all figures and tables, the atomic properties were required to reproduce the four constraints of correct molecular charge and dipole moment components, and the distance-dependent weight was not applied to the potentials.) Another set of atomic properties was created from a subset of the 1597 potentials, for which only the points at distances greater than 1.2A from the molecular van der Waals radius were used (henceforth called the 589 point set). This distance corresponds to the closest approach of other nuclei to the charge centers of this molecule in a dimer (i.e. the smallest RAB in Eq. (11)). A plot of the error in the fitted potential was made for the model containing only atomic charges (Fig. 2) and one with charges and dipoles on the heavy atoms (Fig. 3). The reduction in the error from 20 kcal mol-’ to
23
less than 1 kcalmol-’ between point sets is dramatic. We believe that the large errors in the 1597 point set arose from the inappropriate use of a multipole-fit program to reproduce the potential at points which have penetrated the molecular charge density. The points in the 589 point set do not approach the molecule as closely and the fit is much better. The 589 point set will be used for all further calculations. Comparing Figs. 2 and 3, it may be seen that a charge-dipole model is only slightly superior to a simple charge model. Further, the R2 and SD values of the two models differ by very little. Acetamide atomic charges and dipoles from various types of calculations and the statistical analysis of each fit are listed in Table 2. Experimentation with various values for the distancedependent weighting mentioned above produced no significant changes in the atomic properties. All calculations were then carried out giving all points equal weight. Various acetamide dimer configurations were examined and the electrostatic energy was determined using both the multipole analysis and the Morokuma method. Some dimers were positioned so that the molecular van der Waals radii were just touching, and some configurations were made by moving the dimers apart another 2A. Plots of the error in the electrostatic energy (calculated ES Morokuma ES) against the SCF binding energy (INT) were made for the charge and charge-dipole models generated from the 589 point set (Figs. 4 and 5). Because we saw little difference between the two models, further analysis of the binding energy was based on the smaller set of atomic properties. The interaction energy is predominantly electrostatic and the error seemingly random but predominantly positive (Fig. 6). Plots of the errors in the electrostatic energies against the Morokuma exchange energy, which is related to the overlap of the molecular charge distributions, showed a linear correlation (Fig. 7). This result suggests that the fraction of the electrostatic energy not being accounted for in the multipole analysis is due to charge penetration, which is a function of the overlap of the molecular orbitals. Because the
24
J.C. White and E.R. Daviakon/J. Mol. Strut.
(Theochem) 282 (1993) 19-31
Table 2 Acetamide atomic properties” Property
Atom Hl
H2
H3
H4
H5
Cl
c2
N
0
1597 Point set [l,l]; R2 = 0.9746; max = 0.0481; SD = 0.0103 4 0.33 0.34 0.01 0.00 (0.01) (0.01) (0.01) (0.01)
0.00 (0.01)
- 0.03 (0.03)
0.87 (0.02)
-0.82 (0.02)
- 0.70 (0.01)
1597 Point set [1,4]; R2 = 0.9794; max = 0.0431; SD = 0.0093 4 0.29 0.35 0.05 - 0.06 (0.02) (0.01) (0.01) (0.01)
- 0.06 (0.01)
0.23 (0.04) 0.16 (0.04) - 0.23 (0.03) 0.00 (0.02)
0.59 (0.05) 0.21 (0.10) - 0.48 (0.09) 0.00 (0.02)
- 0.97 (0.04) 0.19 (0.06) -0.15 (0.05) 0.00 (0.01)
- 0.42 (0.05) -0.18 (0.03) -0.31 (0.04) 0.00 (0.01)
- 0.23 (0.04) - 0.06 (0.02) -0.10 (0.02) -0.11 (0.03)
0.64 (0.13) -0.13 (0.06) -0.17 (0.10) 0.00 (0.07)
0.81 (0.06) 0.22 (0.12) -0.51 (0.11) 0.00 (0.02)
- 1.91 (0.07) - 0.43 (0.13) 0.01 (0.11) 0.00 (0.02)
- 0.52 (0.05) -0.11 (0.03) -0.19 (0.04) 0.00 (0.01)
589 Point set [l,l]; R’ = 0.9998; max = 0.0014; SD = 0.0003 4 0.48 0.46 0.18 0.17 (0.00) (0.00) (0.00) (0.00)
0.17 (0.00)
- 0.62 (0.02)
1.03 (0.01)
- 1.18 (0.01)
- 0.68 (0.00)
589 Point set [1,4]; R2 = 0.9999; max = 0.0007; SD = 0.0001 4 0.50 0.49 0.22 0.16 (0.01) (0.00) (0.01) (0.00)
0.16 (0.00)
- 0.59 (0.03) 0.07 (0.02) - 0.08 (0.01) 0.00 (0.00)
1.68 (0.02) 0.64 (0.07) - 0.28 (0.06) 0.00 (0.00)
- 1.65 (0.02) 0.20 (0.04) - 0.42 (0.02) 0.00 (0.00)
- 0.96 (0.03) 0.11 (0.02) 0.27 (0.02) 0.00 (0.00)
1.50 (0.20) 0.16 (0.07) - 0.70 (0.18) 0.00 (0.03)
1.20 (0.07) 1.41 (0.14) -0.15 (0.14) 0.00 (0.01)
- 1.84 (0.05) - 0.30 (0.18) - 0.27 (0.08) 0.00 (0.01)
- 1.04 (0.04) 0.11 (0.02) 0.35 (0.03) 0.00 (0.00)
IJLW
1597 Point set [4,4]; R2 = 0.9840; max = 0.0367; SD = 0.0084 4 0.93 0.59 - 0.08 - 0.23 (0.06) (0.05) (0.06) (0.04) p’x - 0.55 0.19 0.03 - 0.06 (0.03) (0.02) (0.05) (0.02) PY 0.06 0.16 0.08 -0.10 (0.03) (0.02) (0.02) (0.04) p’z 0.00 0.00 0.00 0.11 (0.01) (0.01) (0.01) (0.03)
KY
ikz 589 Point set [4,4]; R’ = 0.9999; max = 0.0005; SD = 0.0001 -0.11 - 0.47 4 0.72 0.49 (0.06) (0.06) (0.07) (0.04) p’x -0.15 0.09 - 0.06 -0.16 (0.04) (0.02) (0.01) (0.02) p”y - 0.03 0.00 0.21 -0.19 (0.03) (0.01) (0.04) (0.01) & 0.00 0.00 0.00 0.33 (0.04) (0.00) (0.00) (0.00)
- 0.47 (0.07) -0.16 (0.02) -0.19 (0.01) 0.33 (0.04)
J.C. White and E.R. Davi&on/J. Mol. Struct. (Theochem) 282 (1993) 19-31
25
Table 2 (continued) Property
Atom Hl
Mulliken 4
H2
H3
H4
H5
Cl
c2
N
0
0.40
0.39
0.24
0.21
0.21
- 0.60
0.53
- 0.84
-0.54
0.26 0.17 0.00 0.00
0.25 - 0.08 -0.14 0.00
0.03 - 0.05 0.15 0.00
-0.01 - 0.07 - 0.08 0.11
-0.01 - 0.07 - 0.08 -0.11
- 0.09 -0.10 0.06 0.00
1.01 0.10 0.00 0.00
- 0.67 0.03 - 0.02 0.00
- 0.75 - 0.02 - 0.03 0.00
0.01 - 0.04 0.15 0.00
0.13 0.03 0.02 0.25
0.13 0.03 0.02 - 0.25
-0.16 0.12 0.11 0.00
0.52 0.12 - 0.04 0.00
- 0.43 0.06 - 0.04 0.00
- 0.54 0.00 - 0.05 0.00
CADPAC 9
p’x p’y p’z
Hirshfeld population analysis 4 0.15 PLX 0.23 IJY - 0.03 k 0.00
0.19 0.07 0.24 0.00
a The first six sets are properties in atomic units deduced from a multipole analysis program. They are constrained to reproduce the total molecular charge and dipole of acetamide. The 1597 set is composed of points about a&amide, approaching the van der Waals radius of the molecule. The 589 point set is a subset of 1597. All points less than the van der Waals radius plus 1.2A of the nearest atom were eliminated. The numbers in the brackets indicate how many properties will be determined for each light and heavy atom, respectively: [ 1,l], charges only on light (H) and heavy (C,N,O) atoms; [ 1,4], charges on light atoms, charges and dipoles on heavy atoms;
[4,4], charges and dipoles on all atoms. SD is defined in Eq. (8); max is the maximum error (in atomic units); the number in parentheses below each value is the uncertainty
in that value given by Eq. (9).
dimer configurations were generated at random, there is no simple correlation between relative orientation and the degree of charge penetration. We would expect a similar relationship between the Morokuma values for charge transfer and exchange energies, because they both vary with molecular overlap. Figure 8 shows that this is a fairly linear relationship. The polarization energy is orientation-specific so we did not expect any simple relationship with the other energy quantities, although it might vary quadratically with field strength, and hence quadratically with the electrostatic energy. Figure 9 is included mainly to give an indication of the magnitude of the polarization energy. One of the assumptions in molecular modeling is that atomic charges are transferable between similar molecules. We decided to explore this in a little more detail with acetone. The multipolefitting scheme was applied to a set of 2014 poten-
tials around
acetone.
Results
are given in Table 3.
Potentials were recalculated when the acetone was replaced with acetamide, positioned so that the
We can see that group does not significantly change the overall charge distribution. Interestingly, the heavy atoms obtained from a charge-only fit of the potentials are nearly ionic. For comparison, we also list the charges produced from the program MMX [151. The size of these charges may be understood if it is remembered that they are derived from Allinger’s MM3 force field, which is based on a C-H methyl group dipole of zero. Also included are charges and atomic moments using the Mulliken definition [ 161, the Stone-Buckingham definition [17] incorporated in CADPAC, and the Hirshfeld definition [18], as incorporated in MELD. The Mulliken definition agrees better than the others with the charges found from the electrostatic potential. molecules
replacing
nearly
overlapped.
one methyl
group
with an -NH,
J.C. White
and
E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
2
T .
.
1
. .
. 0
8
.
.
+ 5
.
l
.
..
.
l
l .
.
. .
.
, l .
. .i
0
. -1
-1
-20
-30
-10
Potentials
0 (589)
10
20
-10
0
-5
in kcal/mol
Morokuma
Fig. 2. Error in recalculated electrostatic potential about acetamide versus the potential for the 589 point set. Potentials were recalculated with the [l,l] atomic properties (see text).
5
10
INT (kcal/mol)
Fig. 4. Difference in electrostatic energy (multipole calculated ES - Morokuma ES) versus the Morokuma interaction energy. The [ 1, l] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
2 1
. . . 1
.
. . 8. 0
.
l
I
.
..
0
.:
. .
r
. .
e-
.
-1 -10
-1 -30
-20
-10
Potentials
0 (559)
10
20
in kcal/mol
Fig. 3. Error in recalculated electrostatic potential about acetamide versus the potential for the 589 point set. Potentials were recalculated with the [1,4] atomic properties (see text).
-5 Morokuma
0
5
10
INT Ikcallmol)
Fig. 5. Difference in electrostatic energy (multipole calculated ES - Morokuma ES) versus the Morokuma interaction energy. The [ 1,4] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
27
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
0.00
nln
-0.10
m
E
# o
i-
0
o
=E
-0.20
¢I
E -,i
I
o•
Inn
,
"0 @
"F
P
,.,,}
o
-0.30
-v -,
"5 "0i
-0.40
o
-0.50 0.00
-1 -10
-5
0
5
0.50
1.00
1.50
2.00
2.50
3.00
3.50
10 Morokuma EX (kcallmol)
Morokuma ES (kcallmol}
Fig. 6. Difference in electrostatic energy (multipole calculated ES - Morokuma ES) versus the Morokuma electrostatic energy. The [1,1] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
Fig. 8. Morokuma charge transfer energy (kcalmol -l) versus Morokuma exchange energy (kcalmol -~) for the acctamide dimer configurations under investigation. (589 point set [1,1]; N C N S T = 4; L = 0.)
!
0.00 111 -0.20
E 0 J a.
!
-0.40
m
E
uJ "o
-0.60
o 0
X "5
-0.80
0
o -1.00 -1
0
i
i
~
I
2
3
4
-1.20 -10
-5
o
5
10
Morokuma EX (kcal/mol)
Fig. 7. Difference in electrostatic energy (multipole calculated E S - Morokuma ES) versus the Morokuma exchange energy. The [1,1] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
Morokuma ES (kcallmol)
Fig. 9. Morokuma polarization energy (kcal mol -]) versus Morokuma electrostatic energy (kcal tool-') for the acetamide dimer configurations under investigation. (589 point set [1,1]; NCNST = 4; L = 0.)
28
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
Table 3 Acetamide and acetone atomic charges q” Molecule
Atom H2
Hl
H3
H4
2014 Point set [I,l]; R’ = 0.9999; max = 0.0005; SD = 0.0001 Acetamide 0.48 0.46 0.16 0.15 (0.00) (0.00) (0.00) (OW MMX (dipole = 3.63 D) Acetamide 0.15 2014 Point set [I,lJ; Acetone
0.15
0.00
0.00
R2 =0.9998; max = 0.0013; SD = 0.0002 0.16 0.15 (0.00) (0.00)
MMX (dipole = 2.80 D) Acetone
0.00
0.00
H5
0.15 (0.00)
Cl
c2
N
0
- 0.56 (0.01)
1.oo (0.00)
- 1.17 (0.00)
- 0.68 (0.00)
-0.17
- 0.39
0.00
0.04
0.15
- 0.57 (0.01)
(0.00) 0.00
0.04
0.22 0.87 (0.00)
- 0.63 (0.00)
0.32
- 0.39
a Using the procedure described in the text, charges were obtained for methane using a 2038 point set. The [ 1, l] model yielded charges of 0.14 (0.0) and - 0.56 (0.0) for the hydrogen and carbon atoms, respectively; R2 was 0.9891; max was 0.0004; SD was 0.00004. (See also footnote to Table 2.)
Electron correlation
The binding energy of various dimers has been dissected into meaningful components, but one important term cannot be obtained from the SCF energy. Configuration interaction (CI) calculations are required to estimate the dispersion energy of the system. We used a doubly hydrogen bonded dimer in which the O..*H bond lengths were about 2.6A to estimate the dispersion energy at one point on the potential surface. Results are listed in Tables 4 and 5. Table 4 illustrates the difficulty in obtaining a precise value for the correlation contribution to the binding energy by conventional methods. The total correlation energy of the dimer estimated by these calculations is about 800 kcal mol-’ and the expected change is less than 2 kcal mol-' . CI calculations using only single and double excitations are well known not to be size-consistent - i.e. the sum of the energies of two monomers calculated separately differs from the energy of the dimer at large internuclear separation. Quadruple corrections, an empirical method justified by perturbation theory, nearly solves this problem so that in Table 4 the super molecule CI energy differs from the sum of the monomers by only 1.1 kcal mol-’ .
Unfortunately, this error is the same size as the effect being calculated. Another source of error in calculations is that the basis set of the dimer, when used for the monomer alone, lowers the monomer energy. With our basis set, this effect is small in the SCF energy, but is 0.6 kcalmol-’ for each monomer in the CI energy. In the CI column of Table 4, this leads to a binding energy of - 7.66 kcal mall’ compared with the SCF binding of - 7.91 kcalmoll’, i.e., the straight CI result gives 0.25 kcal mall’ antibonding from correlation effects. The row labeled CP (basis set counterpoise corrections) includes corrections for the overbinding of the dimer due to the enlarged basis set. The row labeled EXT includes the correction for the non-additivity of monomer energies by using a supermolecule energy computed in the same way as the dimer is calculated. Finally, the line labeled CP + EXT is found assuming the previous two corrections are additive. This still shows net repulsion relative to the SCF result. The linearized couple cluster (LCC) method, and this method with finite electron pair corrections [ 193 are alternative methods for correcting a SDCI calculation for higher excitations. These methods are exactly size-consistent, so the supermolecule energy
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31 Table 4 Acetamide correlation
29
energies”
SCF
HFSDCS’
LC@
FEPCd
Dimer (ABlAB) -416.120184
-417.4696
-417.46993
-417.44651
Monomer (A/A) - 208.053785
- 208.7287
- 208.72822
-208.71651
Monomer (AIAB) - 208.053858
- 208.7297
-208.72913
- 208.71742
- 417.4557
-417.45645
- 417.43302
Supermolecule [(A/A) -416.107571
+ (BIB) /
Interaction energy (kcalmol-‘)’ - 7.91 CPf - 7.82 EXT* -7.91 CP + EXT - 7.82
-
7.66 6.40 8.72 7.46
-
8.46 7.32 8.46 7.32
a All energies in atomic units unless otherwise specified. b Hartree-Fock singles and doubles extrapolated energy with quadrupoles correction. ‘Linearized coupled-cluster energy. d Finite electron pair corrected, LCC. ‘The first row of interaction energies was obtained by subtracting twice the monomer ‘Calculated by_ using- the monomer in the dimer basis numbers, i.e. (A/AB). 8Determined by using the supermolecule values.
is just the sum of the dimer energies. Both of these methods show a net stabilization by electron correlation of about 0.5 kcal mol-’ before correction for basis set superposition errors. This correction, unfortunately, is larger than the effect being computed, so reverses the sign (and probably produces a less accurate result). Another approach to finding the dispersion energy is to partition the configurations into categories. Using orbitals localized on each molecule [6], the configurations can be classed as double excitations of d, double excitations of a’, and simultaneous single excitations of d and 2#. The latter configurations are the traditional dispersion terms. In Table 5 we give the result of calculations using each of these categories of configurations separately. The dispersion terms give a net binding contribution of -0.98 kcalmol-‘. Each of the monomers, however, has the magnitude of its correlation energy decreased by 0.75 kcal mol-' , which results in a net repulsive effect from electron correlation of 0.6 kcal mol-’ . Although we have not
-
8.46 7.32 8.46 7.32
energy from the dimer energy.
examined the question here, for the water dimer this monomer effect was shown to be due to the Pauli exclusion principle, which changes the correlation energy of d by removing from its virtual space the occupied orbitals of ~3. Conclusion
The goal of this study was to investigate the possibility of producing a set of atomic properties from ab initio wavefunctions which could then be used in a semiempirical determination of the electrostatic energy for dimers, and to test their accuracy by comparing them with the corresponding exact values. Results show that properties derived from a least-squares fit do a credable job of reproducing the electrostatic energy for acetamide dimers in a variety of conformations. In fact, a charge-only model did as well as a mixed chargedipole model, with the dominant error apparently due to interpenetration of the charge distributions rather than truncation of the series. This is in
J.C.
30 Table 5 Acetamide dispersion
White and E.R. Davihon/J.
Mol. Struck (Theochem)
282 (1993)
19-31
energy” HFSDCIb
SCF
LCC
FEPC”
-416.12152
-416.12152
-416.79427
- 416.78258
Dimer (ABlAB) Excite one electron from A and one electron from B -416.1215 -416.120184 Excite two electrons from one monomer -416.120184
(A)
into the virtual space
-416.7948
Dispersion energy terms’ Adilllcr A 1110011131~ ’
- 0.6746
- 0.6741
- 0.6624
- 0.6758
- 0.6753
- 0.6627
AE,=
- 0.0012
Dispersion”
- 0.0014
Net correlation
effect
0.0012 - 0.0014
0.0010
0.0010
0.0003 - 0.0014
- 0.0008
“All energies in atomic units unless otherwise specified. Quadruple corrections have been included in the CI energies. bHartree-Fock singles and doubles extrapolated enrgy. c Linearized coupled clusters. d Finite electron pair correction. ‘The CI energy of the dimer with two electrons excited from monomer A into the virtual space minus the SCF energy of the dimer. ‘The CI energy of (A/AB) minus its SCF energy. r The difference of the previous two terms (= AEB). hThe Cl energy of the dimer with one electron excited from each monomer into the virtual space minus the dimer SCF energy.
contrast with other studies, which have found that a multipole expansion through octopoles is required in order to produce a converged energy [20]. The least squares charges were closer to Mulliken charges than to any other commonly used definition. Although the electrostatic interaction dominates at distances as large as the van der Waals radius, it is less dominant at shorter distances. At the typical H-bond distance of 2.0 A for H***O,the ES term no longer dominates. At this shorter distance for the doubly hydrogen bonded dimer, the atomic moments method, including charges, dipoles, and quadrupole moments, gives - 21 kcal mol-’ for ES, compared with - 19 kcalmol-’ from integration of the charge density. However, the high penetration of the charges causes the exchange repulsion to be + 13 kcal mol-’ . The polarization contribution was -4 and the charge transfer was - 2 for a total SCF interaction of - 11 kcal mol-’ . Even at this shorter distance, the correlation contribution to the energy, as well as we can compute it, remains small.
Acknowledgment
The authors gratefully acknowledge financial support of this research by the National Institutes of Health, Grant number PHS ROl GM 34081.
References K.B. Lipkowitz and D.B. Boyd, Reviews in Computational Chemistry, VCH, New York, 1990. N.L. Allinger, J. Am. Chem. Sot., 99 (1977) 8127; 111 (1989) 8551; 111 (1989) 8566; 111 (1989) 8576. P.K. Weiner and P.A. Kollman, J. Comput. Chem., 2 (1981) 287. B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan and M. Karplus, J. Comput. Chem., 4 (1983) 187. S.L. Mayo, B.D. Olafson and W.A. Goddard III, J. Phys. Chem., 94 (1990) 8897. J.C. White and E.R. Davidson, J. Chem. Phys., 93 (1990) 8029. F.B. van Duijneveldt, IBM Research Report RJ 945, 1971.
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31 8 9 10 11 12 13
H. Partridge, J. Chem. Phys., 90 (1989) 1043. M. Kitano and K. Kuchitsu, Bull. Chem. Sot. Jpn., 46 (1973) 3048. A.D. Buckingham, Q. Rev. Chem. Sot., 13 (1959) 183. S.R. Cox and D.E. Williams, J. Comput. Chem., 2 (1981) 304. R.J. Woods, M. Khalil, W. Pell, S.H. Moffat and V.H. Smith, Jr., J. Comput. Chem., 3 (1989) 297. U. Dinur and A.T. Hagler, J. Chem. Phys., 91 (1989) 2949. M. Gutowski and L. Piela, Mol. Phys., 64 (1988) 337. W.A. Sokalalski, S. Roszak and K. Pecul, Chem. Phys. Lett., 153 (1988) 153. I.C. Hayes and A.J. Stone, Mol. Phys., 53 (1984) 83. N. Gresh, P. Claverie and A. Pullman, Int. J. Quantum Chem., 13 (1979) 243. P.A. Kollman and L.C. Allen, Theor. Chim. Acta, 18 (1978) 399. M. Dreyfus and A. Pullman, Theor. Chim. Acta, 19 (1970) 20. J.N. Murrell and G. Shaw, J. Chem. Phys., 46 (1967) 1768.
14 15
16 17
18
19 20
31
H. Umeyama and K. Morokuma, J. Am. Chem. Sot., 99 (1977) 1316. J.J. Gajewski, K.E. Gilbert and J. McKelvey, Advances in Molecular Modeling, Vol. 2, JAI Press, Greenwich, CT, 1990, p. 65. R.S. Mulliken, J. Chem. Phys., 23 (1955) 1833. R.D. Amos and J.E. Rice, CADPAC: The Cambridge Analytic Derivatives Package, Cambridge, 1987. A.D. Buckingham and P.W. Fowler, Can. J. Chem., 63 (1985) 2018. A.J. Stone, Chem. Phys. Lett., 83 (1991) 233. A.J. Stone and M. Alderton, Mol. Phys., 56 (1986) 1047. F.L. Hirshfeld, Theor. Chim. Acta, 44 (1977) 129. E.R. Davidson and S. Chakravorty, Theor. Chim. Acta, 83 (1992) 319. J. Gdanitz and R. Ahlrichs, Chem. Phys. Lett., 143 (1988) 413. J.T. Egan, T.J. Swissler and R. Rein, Int. J. Quantum Chem., 1 (1974) 71. R.L. Ornstein and R. Rein, Biopolymers, 17 (1978) 2341.
19
Comparison of ab initio and multipole determinations of the electrostatic interaction of acetamide dimers Julia C. White, Ernest R. Davidson* Department
of Chemistry,
Indiana
University, Bloomington,
IN 47405,
USA
(Received 1 November 1991)
Abstract In this study, we used a constrained least-squares method to determine parameters for use in an atomic multipole description of the electrostatic energy. To test the accuracy of the parameters, we compared acetamide dimer binding energies determined from the multipole scheme with those from a Morokuma analysis. Interestingly, it was found that a model employing only atomic point charges was able to accurately reproduce the ab initio electrostatic energies at the van der Waals separation of the dimer. The least-squares charges were closer to Mulliken charges than to other definitions which have been suggested. At closer distances, penetration of the charge distributions invalidates the multipole expansion.
Introduction Because the size of many compounds exceeds our ability to perform sophisticated ab initio calculations, considerable effort is being put into the development and improvement of semiempirical formulae and parameters [l]; MMZ [2], AMBER and CHARMM [4] are among the most commonly used semiempirical programs. The quality and usefulness of each is determined by the force field it employs. For example, MM~ and its related programs include more sophisticated algorithms for the description of bond stretching and bending, torsional motion, and so forth, than AMBER or CHARMM. Thus MM~ is more effective in treating medium-sized systems, whereas the strength of AMBER and CHARMM comes from their ability to examine much larger systems. Goddard and coworkers have recently proposed a force field with * Corresponding
author.
parameters defined by atomic bonding or hybridization (e.g. sp’, sp3) [5]. In this paper, we deal with the problem of modeling the electrostatic part of the interaction between hydrogen bonded dimers. For this term, Weiner and Kollman [3] and Karplus and coworkers [4] use a chargesharge formula qAqB/ ERAB9 where E is the dielectric constant of the medium, qA and qB are atom centered point charges and RAB is the distance between the atoms. Allinger [2] employs bond dipoles rather than charges on atoms. All three programs treat hydrogen bonding as a separate type of interaction, often modeling it with a Lennard-Jones (10,12)-potential. Method We chose acetamide (1) as a simple compound, capable of forming hydrogen bonds like those in complex biological molecules and examined the electrostatic interaction within the dimer.
20
J.C.
White and E.R. DavicLFonlJ. Mol. Struct. (Theochem)
H3 0
Hs
,‘
/ Cl
H4
c2
N
Hl
>i-
Hi
Basis set
Restrictions on computer time and space forced us to limit the size of the basis set, and experience with standard basis sets (e.g. 6-31G*) [6] led us to conclude that we should create our own in order to describe the hydrogen bond more accurately. Van Duijneveldt [7] hydrogen (10s) and Partridge [8] carbon, nitrogen, and oxygen (13s,8p) primitives were used to calculate the corresponding atomic wavefunctions. The three most diffuse Partridge s primitives were modified so that the ratio of exponents L/L+ I was 2.3 and a fourth function was added. The same thing was done to the four most diffuse Partridge p primitives, and a fifth was added. In this way, the initial (13s,Sp) set was increased to (14s,9p). Atomic orbital coefficients were used in contractions to yield a triple-zeta valence basis set [4s3p] for each heavy atom. The 1s orbital of the neutral atom formed the 1s basis function. The remaining three s and p basis functions were obtained from the 2s and 2p atomic orbitals of the cationic, anionic, and neutral wavefunctions respectively. The cationic and anionic self-consistent field (SCF) calculations were done with frozen 1s cores from the neutral atom. The [2s] basis set for hydrogen was obtained from the 0 orbital of molecular hydrogen and the 1s orbital of atomic hydrogen. Polarization functions were added to give a final basis set of [4s3pld/2slp]. The energy
282 (1993)
19-3I
lost in the contractions is illustrated in Table 1. The SCF energies of the cation, anion, and neutral states calculated with the [4s,3p] basis agree well with the uncontracted basis set values (within 0.03 kcal mall’), and fairly well with the numerical Hartree-Fock values (within 0.50 kcal mall’ ). Additional states were examined for comparison (e.g. ‘S of neutral carbon). Another motivation for developing our own basis set was the reduction of basis set superposition error, an artificial lowering of the energy because the dimer has more basis functions to utilize than either monomer alone. To check the quality of the final basis set, the SCF energy of the monomer in its own basis (- 208.053785 a.u.) was compared with the energy of the monomer in the full dimer basis (- 208.05385 a.u.) using a doubly hydrogen bonded dimer configuration at an N * * *0 distance of 3.6A. We estimate that the basis set superposition error will be less than 0.1 kcal mall’ for any dimer configuration. The geometry of acetamide was taken from the crystal structure reported by Kitano and Kuchitsu [91. Linear least-squares jit
Potentials were calculated at about 1600 random points around the acetamide molecule using MELDF. The electron charge density for molecule d is defined, in the SCF approximation, as N/2 P.&) = 2 1 y?(r) i=I
(1)
where ‘Pi is an occupied molecular orbital. The electrostatic potential at the point rol due to the charge distribution of molecule d is IL&,) = - 1 P_&)/lr - r,ldr + 1 z,JJL
(2)
where R,, is the distance between point rrr and nucleus A. Buckingham [lo] proposed that the potential at some point in space outside the charge distribution
J.C. White and E.R. Davidson/J. Mol. Strut.
(Theochem)
282 (1993)
19-31
21
Table 1 Atomic SCF energies” [4%3Pl
MPb
Carbon ‘P(C)
- 37.688539
- 37.688557
2P(c+) %(C - ) ‘S(C) ‘D(C)
-
- 37.292166 - 37.708539
Nitrogen %N)
- 54.400802
- 54.40083 1
‘P(N+) ‘P(N-) *D(N) 2P(N)
-
53.887870 54.321322 54.295927 54.227684
- 53.887904 - 54.321345
- 53.887896 - 54.321344
Oxygen 3w)
- 74.809185
%(o+) 2P(O-) ‘S(O) ‘D(O)
-
- 74.809232 - 74.372442 - 74.789030
- 74.372440 - 74.789029
37.292141 37.708527 37.549110 37.631166
74.372392 74.788991 74.610601 73.727459
Frozen corec
Numerical Hartree-Fock
- 37.688619
- 37.292154 - 37.708539
- 37.292224 - 37.708844 - 37.549611 -37.631331 -
54.400934 53.888005 54.321959 54.296169 54.228102
- 74.809398 - 74.372606 - 74.789746 -74.611021 - 74.729264
Hydrogen
H H*
- 0.499998 - 1.128528 4
- 0.499999 - 1.128528
-0.5 - 1.13363d
“All energies are in atomic units. .,_’ b Modified Partridge primitive set, as described in text. ‘Calculations performed with the modified Partridge basis, freezing the 1s core. ‘At 1.4 a,,.
of a set of atoms can be approximated
as
certain matrix elements, (+)oA
+
1 1 &~a m n
%n~ &~a
I&
S=
’’
‘1
=
&ii
C$‘tfl
(3)
where qA, ,uA , and 6Aare the atom centered charges, dipoles, and quadrupoles, respectively, and m, n take the values x, y, z. Using a method previously explored by other authors [l 1,121 we performed a linear least-squares fit of the electrostatic potential at selected points about acetamide to a set of atomic properties. To ensure correct long range behavior of the potential, the set of atomic properties was constrained to reproduce certain moments of the SCF charge density. As an example, consider the simplest case using only atomic charges subject to the constraint of zero total charge. After defining
(V,
@A
=
=
t&,)
qA
and
ce)A
=
1
the least-squares error with a Lagrange multiplier for the total charge constraint can be written as E = (44 - VY(f#q - V) + Iq’e
(4)
We can write the equation for optimum charges derived from this as q = Sm’##2 v-
(1/2)IS-‘e
(5)
so that the final total charge Q is Q = eTq = (t~~S-‘f$~ V) - (1/2)A(eTSP’e)
(6)
22
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
which can be used to find 3, by setting Q to zero. These equations are easily generalized to include atomic dipoles and quadrupoles. The quality of the fit can be measured by the R2 and standard deviation (SD) values. The value of R2 is given by (7)
culate the electrostatic interaction energy of acetamide dimers in a variety of conformations with charge densities frozen at the free molecule SCF value. This electrostatic energy (ES), to which the multipole series converges for large intermolecular separations, is given by ES = 1 G K&)
- 1 p&)Y&)dz
(10)
B
where
where d and $? label the interacting molecules. By expanding each term in a Taylor series in Rik,the ES values may be rewritten in terms of atomic multipoles as
P= N,_’ CK CI
D,=v,-P
i,=N,-‘~+, K* = AA
16A
(11)
and the variance o2 is 2 = (44 - VT’&
- IY(N, - Nr)
(8)
Similarly, the variance in the value of an individual parameter qA can be estimated as a; = (S-‘),$
(9)
In Eq. (8), N, is the number of points at which the potential is sampled and Nr is the number of degrees of freedom in the calculation (the total number of atomic properties minus the constraints). For example, the charge-only model of a&amide with molecular charge and dipole constraints would use an Nr value of 5. One important factor that has been neglected so far is that the sum in the definition of the matrix S will grow like Ni’3 with most of the weight on distant points, if the points CIare spread uniformly over all space. We tried incorporating a distancedependent weight eeLR in the least-squares fit program in order to reduce the contribution of points that were far away. For the selection of points we used, this made no appreciable difference. Multiple determination of electrostatic energy Once the atomic properties were generated from the linear least-squares fit, they were used to cal-
To test the accuracy of the ES values calculated for the various dimer geometries with the atomic multipole expansion, we compared these values (Eq. (11)) with the corresponding ones obtained from ab initio calculations (Eq. (10)). Numerous methods have been proposed for decomposing the interaction energy between non-covalently bonded species [13]. The method of Umeyama and Morokuma [ 141has already been incorporated into the program MELDF, so it was a simple matter to generate a list of Morokuma ES values for a randomly generated set of configurations. These values for ES are identical to those defined by Eq. (10). Results and Discussion Potentials at points distributed about acetamide were calculated using the [4s3pld/2slp] basis set and were symmetric about the mirror plane of the molecule. This set includes points just at the molecular van der Waals radius and up to 15 A away. Atomic charges and dipoles were generated from the set of 1597 points, with the number of constants, NCNST, equal to 4 in order to exactly reproduce the SCF molecular charge and dipole and with L = 0 to set the distance dependent weight to unity. To test the accuracy of the fit, the
J.C. White and E.R. DaviakonlJ. Mol. Struct. (Theochem) 282 (1993) 19-31 30 *
15
10
5
0
-5
10’ -100
-50 Potentials
0
50
100
(159 7) in koal/mol
Fig. 1. Error in recalculated electrostatic potential about acetamide versus the potential for the 1597 point set. Potentials were recalculated with the [I,41 atomic properties (see text).
potentials were recalculated from the atomic properties and a plot made of the magnitude of the error (Fig. 1). (The [m. n] notation is defined in Table 2. In all figures and tables, the atomic properties were required to reproduce the four constraints of correct molecular charge and dipole moment components, and the distance-dependent weight was not applied to the potentials.) Another set of atomic properties was created from a subset of the 1597 potentials, for which only the points at distances greater than 1.2A from the molecular van der Waals radius were used (henceforth called the 589 point set). This distance corresponds to the closest approach of other nuclei to the charge centers of this molecule in a dimer (i.e. the smallest RAB in Eq. (11)). A plot of the error in the fitted potential was made for the model containing only atomic charges (Fig. 2) and one with charges and dipoles on the heavy atoms (Fig. 3). The reduction in the error from 20 kcal mol-’ to
23
less than 1 kcalmol-’ between point sets is dramatic. We believe that the large errors in the 1597 point set arose from the inappropriate use of a multipole-fit program to reproduce the potential at points which have penetrated the molecular charge density. The points in the 589 point set do not approach the molecule as closely and the fit is much better. The 589 point set will be used for all further calculations. Comparing Figs. 2 and 3, it may be seen that a charge-dipole model is only slightly superior to a simple charge model. Further, the R2 and SD values of the two models differ by very little. Acetamide atomic charges and dipoles from various types of calculations and the statistical analysis of each fit are listed in Table 2. Experimentation with various values for the distancedependent weighting mentioned above produced no significant changes in the atomic properties. All calculations were then carried out giving all points equal weight. Various acetamide dimer configurations were examined and the electrostatic energy was determined using both the multipole analysis and the Morokuma method. Some dimers were positioned so that the molecular van der Waals radii were just touching, and some configurations were made by moving the dimers apart another 2A. Plots of the error in the electrostatic energy (calculated ES Morokuma ES) against the SCF binding energy (INT) were made for the charge and charge-dipole models generated from the 589 point set (Figs. 4 and 5). Because we saw little difference between the two models, further analysis of the binding energy was based on the smaller set of atomic properties. The interaction energy is predominantly electrostatic and the error seemingly random but predominantly positive (Fig. 6). Plots of the errors in the electrostatic energies against the Morokuma exchange energy, which is related to the overlap of the molecular charge distributions, showed a linear correlation (Fig. 7). This result suggests that the fraction of the electrostatic energy not being accounted for in the multipole analysis is due to charge penetration, which is a function of the overlap of the molecular orbitals. Because the
24
J.C. White and E.R. Daviakon/J. Mol. Strut.
(Theochem) 282 (1993) 19-31
Table 2 Acetamide atomic properties” Property
Atom Hl
H2
H3
H4
H5
Cl
c2
N
0
1597 Point set [l,l]; R2 = 0.9746; max = 0.0481; SD = 0.0103 4 0.33 0.34 0.01 0.00 (0.01) (0.01) (0.01) (0.01)
0.00 (0.01)
- 0.03 (0.03)
0.87 (0.02)
-0.82 (0.02)
- 0.70 (0.01)
1597 Point set [1,4]; R2 = 0.9794; max = 0.0431; SD = 0.0093 4 0.29 0.35 0.05 - 0.06 (0.02) (0.01) (0.01) (0.01)
- 0.06 (0.01)
0.23 (0.04) 0.16 (0.04) - 0.23 (0.03) 0.00 (0.02)
0.59 (0.05) 0.21 (0.10) - 0.48 (0.09) 0.00 (0.02)
- 0.97 (0.04) 0.19 (0.06) -0.15 (0.05) 0.00 (0.01)
- 0.42 (0.05) -0.18 (0.03) -0.31 (0.04) 0.00 (0.01)
- 0.23 (0.04) - 0.06 (0.02) -0.10 (0.02) -0.11 (0.03)
0.64 (0.13) -0.13 (0.06) -0.17 (0.10) 0.00 (0.07)
0.81 (0.06) 0.22 (0.12) -0.51 (0.11) 0.00 (0.02)
- 1.91 (0.07) - 0.43 (0.13) 0.01 (0.11) 0.00 (0.02)
- 0.52 (0.05) -0.11 (0.03) -0.19 (0.04) 0.00 (0.01)
589 Point set [l,l]; R’ = 0.9998; max = 0.0014; SD = 0.0003 4 0.48 0.46 0.18 0.17 (0.00) (0.00) (0.00) (0.00)
0.17 (0.00)
- 0.62 (0.02)
1.03 (0.01)
- 1.18 (0.01)
- 0.68 (0.00)
589 Point set [1,4]; R2 = 0.9999; max = 0.0007; SD = 0.0001 4 0.50 0.49 0.22 0.16 (0.01) (0.00) (0.01) (0.00)
0.16 (0.00)
- 0.59 (0.03) 0.07 (0.02) - 0.08 (0.01) 0.00 (0.00)
1.68 (0.02) 0.64 (0.07) - 0.28 (0.06) 0.00 (0.00)
- 1.65 (0.02) 0.20 (0.04) - 0.42 (0.02) 0.00 (0.00)
- 0.96 (0.03) 0.11 (0.02) 0.27 (0.02) 0.00 (0.00)
1.50 (0.20) 0.16 (0.07) - 0.70 (0.18) 0.00 (0.03)
1.20 (0.07) 1.41 (0.14) -0.15 (0.14) 0.00 (0.01)
- 1.84 (0.05) - 0.30 (0.18) - 0.27 (0.08) 0.00 (0.01)
- 1.04 (0.04) 0.11 (0.02) 0.35 (0.03) 0.00 (0.00)
IJLW
1597 Point set [4,4]; R2 = 0.9840; max = 0.0367; SD = 0.0084 4 0.93 0.59 - 0.08 - 0.23 (0.06) (0.05) (0.06) (0.04) p’x - 0.55 0.19 0.03 - 0.06 (0.03) (0.02) (0.05) (0.02) PY 0.06 0.16 0.08 -0.10 (0.03) (0.02) (0.02) (0.04) p’z 0.00 0.00 0.00 0.11 (0.01) (0.01) (0.01) (0.03)
KY
ikz 589 Point set [4,4]; R’ = 0.9999; max = 0.0005; SD = 0.0001 -0.11 - 0.47 4 0.72 0.49 (0.06) (0.06) (0.07) (0.04) p’x -0.15 0.09 - 0.06 -0.16 (0.04) (0.02) (0.01) (0.02) p”y - 0.03 0.00 0.21 -0.19 (0.03) (0.01) (0.04) (0.01) & 0.00 0.00 0.00 0.33 (0.04) (0.00) (0.00) (0.00)
- 0.47 (0.07) -0.16 (0.02) -0.19 (0.01) 0.33 (0.04)
J.C. White and E.R. Davi&on/J. Mol. Struct. (Theochem) 282 (1993) 19-31
25
Table 2 (continued) Property
Atom Hl
Mulliken 4
H2
H3
H4
H5
Cl
c2
N
0
0.40
0.39
0.24
0.21
0.21
- 0.60
0.53
- 0.84
-0.54
0.26 0.17 0.00 0.00
0.25 - 0.08 -0.14 0.00
0.03 - 0.05 0.15 0.00
-0.01 - 0.07 - 0.08 0.11
-0.01 - 0.07 - 0.08 -0.11
- 0.09 -0.10 0.06 0.00
1.01 0.10 0.00 0.00
- 0.67 0.03 - 0.02 0.00
- 0.75 - 0.02 - 0.03 0.00
0.01 - 0.04 0.15 0.00
0.13 0.03 0.02 0.25
0.13 0.03 0.02 - 0.25
-0.16 0.12 0.11 0.00
0.52 0.12 - 0.04 0.00
- 0.43 0.06 - 0.04 0.00
- 0.54 0.00 - 0.05 0.00
CADPAC 9
p’x p’y p’z
Hirshfeld population analysis 4 0.15 PLX 0.23 IJY - 0.03 k 0.00
0.19 0.07 0.24 0.00
a The first six sets are properties in atomic units deduced from a multipole analysis program. They are constrained to reproduce the total molecular charge and dipole of acetamide. The 1597 set is composed of points about a&amide, approaching the van der Waals radius of the molecule. The 589 point set is a subset of 1597. All points less than the van der Waals radius plus 1.2A of the nearest atom were eliminated. The numbers in the brackets indicate how many properties will be determined for each light and heavy atom, respectively: [ 1,l], charges only on light (H) and heavy (C,N,O) atoms; [ 1,4], charges on light atoms, charges and dipoles on heavy atoms;
[4,4], charges and dipoles on all atoms. SD is defined in Eq. (8); max is the maximum error (in atomic units); the number in parentheses below each value is the uncertainty
in that value given by Eq. (9).
dimer configurations were generated at random, there is no simple correlation between relative orientation and the degree of charge penetration. We would expect a similar relationship between the Morokuma values for charge transfer and exchange energies, because they both vary with molecular overlap. Figure 8 shows that this is a fairly linear relationship. The polarization energy is orientation-specific so we did not expect any simple relationship with the other energy quantities, although it might vary quadratically with field strength, and hence quadratically with the electrostatic energy. Figure 9 is included mainly to give an indication of the magnitude of the polarization energy. One of the assumptions in molecular modeling is that atomic charges are transferable between similar molecules. We decided to explore this in a little more detail with acetone. The multipolefitting scheme was applied to a set of 2014 poten-
tials around
acetone.
Results
are given in Table 3.
Potentials were recalculated when the acetone was replaced with acetamide, positioned so that the
We can see that group does not significantly change the overall charge distribution. Interestingly, the heavy atoms obtained from a charge-only fit of the potentials are nearly ionic. For comparison, we also list the charges produced from the program MMX [151. The size of these charges may be understood if it is remembered that they are derived from Allinger’s MM3 force field, which is based on a C-H methyl group dipole of zero. Also included are charges and atomic moments using the Mulliken definition [ 161, the Stone-Buckingham definition [17] incorporated in CADPAC, and the Hirshfeld definition [18], as incorporated in MELD. The Mulliken definition agrees better than the others with the charges found from the electrostatic potential. molecules
replacing
nearly
overlapped.
one methyl
group
with an -NH,
J.C. White
and
E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
2
T .
.
1
. .
. 0
8
.
.
+ 5
.
l
.
..
.
l
l .
.
. .
.
, l .
. .i
0
. -1
-1
-20
-30
-10
Potentials
0 (589)
10
20
-10
0
-5
in kcal/mol
Morokuma
Fig. 2. Error in recalculated electrostatic potential about acetamide versus the potential for the 589 point set. Potentials were recalculated with the [l,l] atomic properties (see text).
5
10
INT (kcal/mol)
Fig. 4. Difference in electrostatic energy (multipole calculated ES - Morokuma ES) versus the Morokuma interaction energy. The [ 1, l] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
2 1
. . . 1
.
. . 8. 0
.
l
I
.
..
0
.:
. .
r
. .
e-
.
-1 -10
-1 -30
-20
-10
Potentials
0 (559)
10
20
in kcal/mol
Fig. 3. Error in recalculated electrostatic potential about acetamide versus the potential for the 589 point set. Potentials were recalculated with the [1,4] atomic properties (see text).
-5 Morokuma
0
5
10
INT Ikcallmol)
Fig. 5. Difference in electrostatic energy (multipole calculated ES - Morokuma ES) versus the Morokuma interaction energy. The [ 1,4] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
27
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
0.00
nln
-0.10
m
E
# o
i-
0
o
=E
-0.20
¢I
E -,i
I
o•
Inn
,
"0 @
"F
P
,.,,}
o
-0.30
-v -,
"5 "0i
-0.40
o
-0.50 0.00
-1 -10
-5
0
5
0.50
1.00
1.50
2.00
2.50
3.00
3.50
10 Morokuma EX (kcallmol)
Morokuma ES (kcallmol}
Fig. 6. Difference in electrostatic energy (multipole calculated ES - Morokuma ES) versus the Morokuma electrostatic energy. The [1,1] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
Fig. 8. Morokuma charge transfer energy (kcalmol -l) versus Morokuma exchange energy (kcalmol -~) for the acctamide dimer configurations under investigation. (589 point set [1,1]; N C N S T = 4; L = 0.)
!
0.00 111 -0.20
E 0 J a.
!
-0.40
m
E
uJ "o
-0.60
o 0
X "5
-0.80
0
o -1.00 -1
0
i
i
~
I
2
3
4
-1.20 -10
-5
o
5
10
Morokuma EX (kcal/mol)
Fig. 7. Difference in electrostatic energy (multipole calculated E S - Morokuma ES) versus the Morokuma exchange energy. The [1,1] atomic properties obtained from the 589 point set were used in the multipole determination of ES (see text). (NCNST = 4; L = 0.)
Morokuma ES (kcallmol)
Fig. 9. Morokuma polarization energy (kcal mol -]) versus Morokuma electrostatic energy (kcal tool-') for the acetamide dimer configurations under investigation. (589 point set [1,1]; NCNST = 4; L = 0.)
28
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31
Table 3 Acetamide and acetone atomic charges q” Molecule
Atom H2
Hl
H3
H4
2014 Point set [I,l]; R’ = 0.9999; max = 0.0005; SD = 0.0001 Acetamide 0.48 0.46 0.16 0.15 (0.00) (0.00) (0.00) (OW MMX (dipole = 3.63 D) Acetamide 0.15 2014 Point set [I,lJ; Acetone
0.15
0.00
0.00
R2 =0.9998; max = 0.0013; SD = 0.0002 0.16 0.15 (0.00) (0.00)
MMX (dipole = 2.80 D) Acetone
0.00
0.00
H5
0.15 (0.00)
Cl
c2
N
0
- 0.56 (0.01)
1.oo (0.00)
- 1.17 (0.00)
- 0.68 (0.00)
-0.17
- 0.39
0.00
0.04
0.15
- 0.57 (0.01)
(0.00) 0.00
0.04
0.22 0.87 (0.00)
- 0.63 (0.00)
0.32
- 0.39
a Using the procedure described in the text, charges were obtained for methane using a 2038 point set. The [ 1, l] model yielded charges of 0.14 (0.0) and - 0.56 (0.0) for the hydrogen and carbon atoms, respectively; R2 was 0.9891; max was 0.0004; SD was 0.00004. (See also footnote to Table 2.)
Electron correlation
The binding energy of various dimers has been dissected into meaningful components, but one important term cannot be obtained from the SCF energy. Configuration interaction (CI) calculations are required to estimate the dispersion energy of the system. We used a doubly hydrogen bonded dimer in which the O..*H bond lengths were about 2.6A to estimate the dispersion energy at one point on the potential surface. Results are listed in Tables 4 and 5. Table 4 illustrates the difficulty in obtaining a precise value for the correlation contribution to the binding energy by conventional methods. The total correlation energy of the dimer estimated by these calculations is about 800 kcal mol-’ and the expected change is less than 2 kcal mol-' . CI calculations using only single and double excitations are well known not to be size-consistent - i.e. the sum of the energies of two monomers calculated separately differs from the energy of the dimer at large internuclear separation. Quadruple corrections, an empirical method justified by perturbation theory, nearly solves this problem so that in Table 4 the super molecule CI energy differs from the sum of the monomers by only 1.1 kcal mol-’ .
Unfortunately, this error is the same size as the effect being calculated. Another source of error in calculations is that the basis set of the dimer, when used for the monomer alone, lowers the monomer energy. With our basis set, this effect is small in the SCF energy, but is 0.6 kcalmol-’ for each monomer in the CI energy. In the CI column of Table 4, this leads to a binding energy of - 7.66 kcal mall’ compared with the SCF binding of - 7.91 kcalmoll’, i.e., the straight CI result gives 0.25 kcal mall’ antibonding from correlation effects. The row labeled CP (basis set counterpoise corrections) includes corrections for the overbinding of the dimer due to the enlarged basis set. The row labeled EXT includes the correction for the non-additivity of monomer energies by using a supermolecule energy computed in the same way as the dimer is calculated. Finally, the line labeled CP + EXT is found assuming the previous two corrections are additive. This still shows net repulsion relative to the SCF result. The linearized couple cluster (LCC) method, and this method with finite electron pair corrections [ 193 are alternative methods for correcting a SDCI calculation for higher excitations. These methods are exactly size-consistent, so the supermolecule energy
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31 Table 4 Acetamide correlation
29
energies”
SCF
HFSDCS’
LC@
FEPCd
Dimer (ABlAB) -416.120184
-417.4696
-417.46993
-417.44651
Monomer (A/A) - 208.053785
- 208.7287
- 208.72822
-208.71651
Monomer (AIAB) - 208.053858
- 208.7297
-208.72913
- 208.71742
- 417.4557
-417.45645
- 417.43302
Supermolecule [(A/A) -416.107571
+ (BIB) /
Interaction energy (kcalmol-‘)’ - 7.91 CPf - 7.82 EXT* -7.91 CP + EXT - 7.82
-
7.66 6.40 8.72 7.46
-
8.46 7.32 8.46 7.32
a All energies in atomic units unless otherwise specified. b Hartree-Fock singles and doubles extrapolated energy with quadrupoles correction. ‘Linearized coupled-cluster energy. d Finite electron pair corrected, LCC. ‘The first row of interaction energies was obtained by subtracting twice the monomer ‘Calculated by_ using- the monomer in the dimer basis numbers, i.e. (A/AB). 8Determined by using the supermolecule values.
is just the sum of the dimer energies. Both of these methods show a net stabilization by electron correlation of about 0.5 kcal mol-’ before correction for basis set superposition errors. This correction, unfortunately, is larger than the effect being computed, so reverses the sign (and probably produces a less accurate result). Another approach to finding the dispersion energy is to partition the configurations into categories. Using orbitals localized on each molecule [6], the configurations can be classed as double excitations of d, double excitations of a’, and simultaneous single excitations of d and 2#. The latter configurations are the traditional dispersion terms. In Table 5 we give the result of calculations using each of these categories of configurations separately. The dispersion terms give a net binding contribution of -0.98 kcalmol-‘. Each of the monomers, however, has the magnitude of its correlation energy decreased by 0.75 kcal mol-' , which results in a net repulsive effect from electron correlation of 0.6 kcal mol-’ . Although we have not
-
8.46 7.32 8.46 7.32
energy from the dimer energy.
examined the question here, for the water dimer this monomer effect was shown to be due to the Pauli exclusion principle, which changes the correlation energy of d by removing from its virtual space the occupied orbitals of ~3. Conclusion
The goal of this study was to investigate the possibility of producing a set of atomic properties from ab initio wavefunctions which could then be used in a semiempirical determination of the electrostatic energy for dimers, and to test their accuracy by comparing them with the corresponding exact values. Results show that properties derived from a least-squares fit do a credable job of reproducing the electrostatic energy for acetamide dimers in a variety of conformations. In fact, a charge-only model did as well as a mixed chargedipole model, with the dominant error apparently due to interpenetration of the charge distributions rather than truncation of the series. This is in
J.C.
30 Table 5 Acetamide dispersion
White and E.R. Davihon/J.
Mol. Struck (Theochem)
282 (1993)
19-31
energy” HFSDCIb
SCF
LCC
FEPC”
-416.12152
-416.12152
-416.79427
- 416.78258
Dimer (ABlAB) Excite one electron from A and one electron from B -416.1215 -416.120184 Excite two electrons from one monomer -416.120184
(A)
into the virtual space
-416.7948
Dispersion energy terms’ Adilllcr A 1110011131~ ’
- 0.6746
- 0.6741
- 0.6624
- 0.6758
- 0.6753
- 0.6627
AE,=
- 0.0012
Dispersion”
- 0.0014
Net correlation
effect
0.0012 - 0.0014
0.0010
0.0010
0.0003 - 0.0014
- 0.0008
“All energies in atomic units unless otherwise specified. Quadruple corrections have been included in the CI energies. bHartree-Fock singles and doubles extrapolated enrgy. c Linearized coupled clusters. d Finite electron pair correction. ‘The CI energy of the dimer with two electrons excited from monomer A into the virtual space minus the SCF energy of the dimer. ‘The CI energy of (A/AB) minus its SCF energy. r The difference of the previous two terms (= AEB). hThe Cl energy of the dimer with one electron excited from each monomer into the virtual space minus the dimer SCF energy.
contrast with other studies, which have found that a multipole expansion through octopoles is required in order to produce a converged energy [20]. The least squares charges were closer to Mulliken charges than to any other commonly used definition. Although the electrostatic interaction dominates at distances as large as the van der Waals radius, it is less dominant at shorter distances. At the typical H-bond distance of 2.0 A for H***O,the ES term no longer dominates. At this shorter distance for the doubly hydrogen bonded dimer, the atomic moments method, including charges, dipoles, and quadrupole moments, gives - 21 kcal mol-’ for ES, compared with - 19 kcalmol-’ from integration of the charge density. However, the high penetration of the charges causes the exchange repulsion to be + 13 kcal mol-’ . The polarization contribution was -4 and the charge transfer was - 2 for a total SCF interaction of - 11 kcal mol-’ . Even at this shorter distance, the correlation contribution to the energy, as well as we can compute it, remains small.
Acknowledgment
The authors gratefully acknowledge financial support of this research by the National Institutes of Health, Grant number PHS ROl GM 34081.
References K.B. Lipkowitz and D.B. Boyd, Reviews in Computational Chemistry, VCH, New York, 1990. N.L. Allinger, J. Am. Chem. Sot., 99 (1977) 8127; 111 (1989) 8551; 111 (1989) 8566; 111 (1989) 8576. P.K. Weiner and P.A. Kollman, J. Comput. Chem., 2 (1981) 287. B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan and M. Karplus, J. Comput. Chem., 4 (1983) 187. S.L. Mayo, B.D. Olafson and W.A. Goddard III, J. Phys. Chem., 94 (1990) 8897. J.C. White and E.R. Davidson, J. Chem. Phys., 93 (1990) 8029. F.B. van Duijneveldt, IBM Research Report RJ 945, 1971.
J.C. White and E.R. Davidson/J. Mol. Struct. (Theochem) 282 (1993) 19-31 8 9 10 11 12 13
H. Partridge, J. Chem. Phys., 90 (1989) 1043. M. Kitano and K. Kuchitsu, Bull. Chem. Sot. Jpn., 46 (1973) 3048. A.D. Buckingham, Q. Rev. Chem. Sot., 13 (1959) 183. S.R. Cox and D.E. Williams, J. Comput. Chem., 2 (1981) 304. R.J. Woods, M. Khalil, W. Pell, S.H. Moffat and V.H. Smith, Jr., J. Comput. Chem., 3 (1989) 297. U. Dinur and A.T. Hagler, J. Chem. Phys., 91 (1989) 2949. M. Gutowski and L. Piela, Mol. Phys., 64 (1988) 337. W.A. Sokalalski, S. Roszak and K. Pecul, Chem. Phys. Lett., 153 (1988) 153. I.C. Hayes and A.J. Stone, Mol. Phys., 53 (1984) 83. N. Gresh, P. Claverie and A. Pullman, Int. J. Quantum Chem., 13 (1979) 243. P.A. Kollman and L.C. Allen, Theor. Chim. Acta, 18 (1978) 399. M. Dreyfus and A. Pullman, Theor. Chim. Acta, 19 (1970) 20. J.N. Murrell and G. Shaw, J. Chem. Phys., 46 (1967) 1768.
14 15
16 17
18
19 20
31
H. Umeyama and K. Morokuma, J. Am. Chem. Sot., 99 (1977) 1316. J.J. Gajewski, K.E. Gilbert and J. McKelvey, Advances in Molecular Modeling, Vol. 2, JAI Press, Greenwich, CT, 1990, p. 65. R.S. Mulliken, J. Chem. Phys., 23 (1955) 1833. R.D. Amos and J.E. Rice, CADPAC: The Cambridge Analytic Derivatives Package, Cambridge, 1987. A.D. Buckingham and P.W. Fowler, Can. J. Chem., 63 (1985) 2018. A.J. Stone, Chem. Phys. Lett., 83 (1991) 233. A.J. Stone and M. Alderton, Mol. Phys., 56 (1986) 1047. F.L. Hirshfeld, Theor. Chim. Acta, 44 (1977) 129. E.R. Davidson and S. Chakravorty, Theor. Chim. Acta, 83 (1992) 319. J. Gdanitz and R. Ahlrichs, Chem. Phys. Lett., 143 (1988) 413. J.T. Egan, T.J. Swissler and R. Rein, Int. J. Quantum Chem., 1 (1974) 71. R.L. Ornstein and R. Rein, Biopolymers, 17 (1978) 2341.