Intermolecular Potentials of the Methane Dimer Calculated with Møller-

Comment on "Intermolecular Interaction Potential of the Methane Dimer From the Local Density Approximation”

1. Introduction

To verify the recently calculated intermolecular interaction potentials of the methane dimer within the density functional theory using the (Perdew) local density approximation (LDA) [Chen et al., Phys. Rev. A 69, 034701 (2004)], we have performed a parallel series of calculations using the LDA/6-311++G (3df, 3pd) level of theory with selected exchange functionals (B, G96, MPW, O, PBE, PW91, S, and XA). None of the above calculated intermolecular interaction potentials from the local density approximation reproduce the results reported in the commented paper. In addition, we point out the inappropriateness of using the Lennard-Jones function to model the long-range parts of the calculated intermolecular interaction potentials, as suggested positively by Chen et al.

Chen et al. [1] recently reported the intermolecular interaction potentials of the methane dimer (CH4)2 within the standard density functional theory (DFT) scheme [2]

using the (Perdew) local density approximation (LDA) [3], the pseudopotential [4,5], and the plane-wave expansion [6]. Their results agreed surprisingly well with those

obtained by the correlation-corrected Møller-Plesset (MP2, MP3) [7] and coupled cluster (CCSD(T)) [8] methods using a large basis set [9]. Because it has been known for some time that the usual DFT based approaches, using either the LDA or the generalized gradient approximation (GGA), can not calculate the intermolecular interaction potentials of molecular dimers to such a high level of accuracy [10-12], it is important to perform a parallel series of calculations using the available implementations of commonly used exchange-correlation functionals to verify the proposed results by Chen et al.

Intermolecular interaction potentials, or van der Waals interactions, or non-covalent-bonded interactions, play an essential role in condensed matter physics, materials chemistry and structural biology. While these interactions are normally one or two orders of magnitude weaker than typical covalent bonds, they are crucial in determining the thermodynamic properties of molecular liquids and solids [13], the energy transfers among molecular complexes [14], and the conformational tertiary structures of macromolecules such as protein and DNA [15]. Unlike intramolecular covalent bonds, intermolecular bonds do not originate from sharing of electrons but rather arise from simultaneous electron correlation of the separated subsystems [16].

Different from stiff covalent bonds, they are relatively soft and non-rigid. Early studies of intermolecular interactions can be traced back to one century ago [17],

while measurements of these interactions are still challenging in the present time [18].

The main difficulty in determining intermolecular interactions experimentally resides at limited samplings of the potential energy surface. For example, experiments using the X-ray crystallography or the laser spectroscopy mainly explore the equilibrium regions of the potential, while thermodynamic measurements in the gas or liquid phase often yield isotropic potential data without the desired stereo-chemical responses. In addition, the measured potentials sensitively depend on the thermodynamic conditions. Usually two measurements carried out in different conditions cannot be compared directly but rely on auxiliary theoretical modeling.

Alternatively intermolecular potentials can be calculated in terms of correlation-corrected quantum chemistry methods [19-21] or density functional theory (DFT) [22-23]. These quantum mechanics based potentials are requested by ab initio molecular dynamics simulations [24] and by classical molecular simulations using force field constructions [25]. Among the components of an intermolecular interaction, the London dispersion force is the most difficult to calculate. The reason is that dispersion interactions arise from the non-local “dynamic” correlations [26]. This non-locality demands full exploration of the time-dependent perspective of quantum mechanics. Often an electron correlation-corrected method and a large basis set are

required to obtain accurate dispersion forces [27]. Also a technical note is in order.

Most present implementations of quantum chemistry programs utilize Gaussian type functions to fasten the calculations of Coulomb repulsion integrals. Because Gaussian type functions are local functions, a large basis set is indispensable to perform a correlation energy calculation. Moreover, these functions do not have the correct asymptotic behavior as the intermolecular separation becomes large. Therefore, the basis set limit of the calculated potential must be estimated so as to be consistent with the conventional perturbation theory.

For general polar molecular systems, the relatively weak dispersion energy is masked by the competing electrostatic energies and hydrogen bond interactions. Nonpolar atomic and molecular dimers are usually taken as a prototype case to study the dispersion energy. Many previous studies on dispersion forces have focused on atomic inert gas dimers and several important conclusions have been drawn from the calculations [28]. There are, however, comparatively fewer detailed studies on whether similar conclusions can be applied to molecular dimers. Because of the extra degrees of freedom and the stereo-chemical responses, the conclusions about atomic dimers may need to be extended or modified in dealing with molecular dimers.

Methane is a non-polar molecule with vanishing dipole and quadrupole moments and

the first nonvanishing electrostatic interaction is the octopole-octopole interaction.

The higher order electrostatic interactions are rather weak and decay fast at large intermolecular separation. The dominant long-range attraction for the methane dimer is thus due to the London dispersion force. On the other hand, the strong repulsive force almost comes from the exchange-repulsion interaction due to the overlapping of electron clouds. Because the exchange-repulsion interactions have been incorporated in the Hatree-Fock (HF) self-consistent theory, post-HF methods such as the Møller-Plesset (MP) perturbation theory and the coupled cluster (CC) theory are often used to calculate the correlation effect. Contrasting both sets of calculation helps to delineate the relative importance of the dispersion energy in the overall intermolecular interaction. In addition, the interaction potentials of alkanes, among them methane being the smallest model, is very crucial in determining the packing morphology in solids and liquids and in lipid bilayers [29-31]. The potentials are also requested for mesoscale simulations for macromolecules [32] because many polymers contain alkyl groups as their backbone units or side chains. Therefore, the calculation of intermolecular interactions of the methane dimer is a “must” and serves as a prototype example to start to investigate the various factors affecting the calculations of these interactions.

There have been many quantum chemistry studies of intermolecular potentials of the methane dimer using correlation-corrected methods [33-44]. Szczesniak et al. [33]

have used the MP2 method with the Sadlej basis to calculate the interaction energies for the six conformers of the methane dimer. They found that the minimum-energy conformation is the D3d conformer and the dimer structure is determined by minimizing the steric repulsion between hydrogen atoms belonging to opposite subsystems. Tsuzuki et al. [36] have studied the basis set effects using basis sets from 6-31G* to 6-311++G (2d,2p) to calculate the interaction energies for two conformers of the methane dimer (but none of them is the D3d conformer). They observed little basis set effect on the HF calculations as long as one uses a basis set larger than the 6-31G* one. The basis set effect is significant for the MP2 interaction energies and the basis set superposition error (BSSE). The dispersion energy can be seriously underestimated if a smaller (than 6-31G*) basis set has been used. They found that the BSSE uncorrected interaction energies do not systematically converge to a destined value, in contrast to those with BSSE corrections. Tsuzuki et al. [39, 43] studied twelve conformers and verified that the D3d conformer corresponds to the minimum-energy geometry. They observed that a large basis set with multiple polarization functions is necessary to evaluate the dispersion energy accurately. They found that augmentation of the diffuse d and p functions to the 6-311G** basis set

more efficiently yields the dispersion energy. Tsuzuki et al. [40, 43] explored the effect of the choice of the correlation-correction methods on the calculations of intermolecular interactions. They demonstrated that the MP2 and MP3 energies are not too far away from the higher level MP4 (SDTQ) calculations, while the latter is not less expensive than the CCSD(T) calculation. They tested the DFT using the BLYP, BPW91 and B3LYP functionals but found unbound interactions while the PW91 functional underestimates only 8% of the potential well depth. They suggest that DFT with the PW91 functional could be an alternative to the ab initio methods.

Recently, Tsuzuki et al. [44] have estimated the MP2 and CCSD(T) interaction energies of the n-alkane dimers at the basis set limit using Dunning’s correlation consistent basis sets.

Many previous studies mainly focused on the equilibrium region of the potentials with relatively few discussions of the full potential curves. Nevertheless, to construct a reliable force field model for molecular simulations, the full intermolecular potential surfaces are required. In this paper we perform a comprehensive up-to-date study on interaction potentials of the prototype methane dimer in terms of the HF, MP2 and DFT methods to gain more understanding of this system. With current computational powers, a detailed editing of the potential data base can be obtained for small size

molecular clusters. It is thus so important to obtain general features of the calculations that we can follow to explore large scale molecular simulations via similar procedures.

The purpose of this paper is twofold. First, we use the state-of-the-art methodology to obtain accurate potential energies for the methane dimer. We would like to verify or modify the previous conclusions about the basis set effects and the effect of including the BSSE on the calculation details of the intermolecular interactions. The basis set effects on repulsion exponents, equilibrium bond lengths, binding energies, and asymptotic coefficients of the calculated intermolecular potentials are thoroughly studied. This is achieved using basis sets from STO-3G [45] to aug-cc-pVQZ [46].

The full potential curves are presented in order to see the overall scope of the potential. In particular, both the BSSE corrected and uncorrected results are presented to emphasize the importance of these corrections. Second, this paper attempts to re-assess the utilities of using the available implementation of the density functional theory in determining the intermolecular interactions. From the studies of atomic dimers, it has been found that conventional DFT based on the local density approximation (LDA) and generalized gradient approximation (GGA) cannot calculate the intermolecular interactions to a satisfying level of accuracy [47]. To address the title issue, we carry out a systematic DFT study of the equilibrium binding energies and bond lengths of the methane dimer using 90 functionals. Methane is a

non-polar molecule with vanishing dipole and quadrupole moments. The first nonvanishing electrostatic interaction is the octopole-octopole interaction and all higher order interactions are weak and decay fast at large intermolecular separation.

The dominant attraction for the methane dimer is thus due to the van der Waals force.

Therefore, the calculation of van der Waals interactions of the methane dimer serves as a prototype study to investigate the various factors affecting the calculations of these interactions.

2. Methods and Calculations

All the calculations were performed using the Gaussian03 package suite [48] and followed a theoretical procedure similar to that employed by Tsuzuki et al. [9] Fig.

1.1-1 shows the calculated interaction potentials of (CH4)2 with a set of selected exchange functionals (B, G96, MPW, O, PBE, PW91, S, and XA) [49], together with the Perdew correlation functional dubbed as PL (Perdew local) [3]. Although we have used a pretty large basis set, 6-311++G (3df, 3pd), which has been shown to lead to convergent results for (CH4)2 at a chemical accuracy [50], none of the calculated intermolecular interaction potentials reproduce the results of Chen et. al. A puzzling point in the commented paper is that there are two sets of data, one reported in their Fig. 1.2-1 and the other as numerical data in the text [1], while the latter is twice the former. Because there was no further clarification on this apparent inconsistency in

the commented paper, we present both data for comparison in Fig. 1.2-1 (open symbol-lines). Restated, neither of them can be reproduced in the present calculations.

In the case of methane dimers, a large part of the exchange-repulsion interactions can be calculated by the HF method. The calculation of electron correlation energies depends on the level of the correlation-corrected method, the size of the basis set, and the correction of the BSSE. The state-of-the-art choice of the correlation-corrected method is either the Møller-Plesset (MPx, x=2, 3, 4) perturbation method [51] or the coupled cluster method with iterative single and double substitutions and with noniterative triple excitations [CCSD(T)] method [52]. It has been found that the MP2 results for the methane dimer are not too much different from those calculated by the much more expensive CCSD(T) as long as a large basis set has been used [43]. To study the basis set effects, we have employed comprehensive basis sets from the Slater-type orbitals fitted with Gaussian functions (STO-nG, n=3~6) [45], Pople’s medium size basis sets [up to 6-311++G (3df, 3pd)] [53] to Dunning’s correlation consistent basis sets (aug-cc-pVXZ, X=D, T, Q) [46]. The basis set superposition error (BSSE) was corrected by the counterpoise (CP) method of Boys and Bernardi [54]. The MP2 interaction potentials at the basis set limit have been estimated using the methods of Helgaker et al. [55] and Feller [56] and a numerical extrapolation scheme based on the Lagrangian formula [57].

All the HF, MP2 and DFT calculations are performed using the Gaussian 03 program package [58] on a single node 2-processor AMD 250 PC cluster with distributed memory. The equilibrium geometry of a single methane molecule was first optimized at the MP2/6-31G* level of theory. To obtain the most stable intermolecular geometry, the methane dimer has been modeled by first fixing the carbon-carbon (C-C) distance while letting the two monomers to rotate freely. By approaching the monomers from the far side with several initial choices of mutual orientation, we found the minimum-energy conformation corresponds to the D3d symmetry conformer. This optimized conformer has been reached through the interplay of the steric stabilization of repulsive hydrogens in opposite monomers [33]. Subsequently the C-C distance was sampled in step 0.1 Å for a quite large range of intermolecular separation (normally 3~9 Å). During the scan we allow the individual methane molecule to be fully relaxed. This means that we do not fix the monomer geometry and the methane molecule is not assumed to be rigid. Although it is not expected to see much deviation from the rigid molecule approximation, in the real condensed phase environment, stretching, bending and torsional relaxations could be important for many subtle thermodynamic properties. The inclusion of intramolecular relaxation is especially relevant to the construction of force fields for use in molecular dynamics simulations

where flexible models often work better than rigid models [59].

The density functionals used in the present work include the 90 combinations chosen among 9 exchange (B88 [60], OPTX [61], MPW [62], PBE [63], PW91 [64], TPSS [65] , Slater [66], HCTH [67], XAlpha [68]) and 10 correlation (TPSS [65], PBE [63], PW91 [64], P86 [69], HCTH [67], VWN5 [70], PL [71], VWN [70], LYP [72]) functionals. We also consider several hybrid functionals of B971 [73], BB98 [74], BHandH [75], O3LYP [76], B3PW91 [77], PBE1PBE [63], MPW1PW91 [78]. The chosen functionals are selective representations of the most commonly used density functionals for van der Waals interactions in current literature. Recent studies showed that the PW91PW91 functional could yield reasonable binding energy of the methane dimer interaction [79] but the relative performance of the exchange and correlation functionals has not been systematically studied.

3. Results and Discussions

(a) Chen et al. also concluded that through a nonlinear fitting their calculated intermolecular potentials can be well described by the Lennard-Jones (L-J) potential function Because this conclusion is contrary to what has been believed that results based on the

LDA can not be used to model the long-range dispersion interaction well [16-18], the determined accuracy from their calculations remains to be verified. For the sake of comparison, in Fig. 1.2-2 we present the calculated raw data and the claimed fitting curve by Chen et al. using their fitting values of a and b [1]. To our great surprise, the fitting curve is anything but like the calculated raw data. To clarify this point, we perform a nonlinear fitting of their calculated data to the L-J function and obtain a=2.0910⁶ Ź² kcal/mol, b=1.8410³ Å⁶kcal/mol, and the fitting is shown in Fig.

1.1-2 As expected, although the L-J function can model the strong repulsive part quite well, there is a significant discrepancy from the calculated data for the long-range interaction part (R > 4 Å). The calculated data using the LDA often decays faster than -1/R⁶ for the long-range part, due to the local nature of the functionals used. To demonstrate this point, we perform another nonlinear fitting using the exponential function can be seen in Fig. 1.2-2, the long-range part of the calculated data is well modeled by the fast-decaying exponential function, but not the L-J function.

(b) The intermolecular interaction potentials of the D3d conformer of the methane

dimer have been calculated with the HF, the MP2 and the DFT methods. We present the respective results along with discussions and make comparisons among the results.

A. Hartree-Fock self-consistent field calculations

The BSSE corrected HF interaction potentials of the methane dimer using several basis sets are shown in Fig. 1.2-3. All the HF calculations yield purely repulsive potentials without minima for all the basis sets used. This can be attributed to the rather weak electrostatic interaction for the methane dimer. In the short range, the strong exchange-repulsion interaction dominates with little alternation from electrostatic and induction attractions. The HF potential is insensitive to the basis size as long as the 6-31G** basis set has been used. We can model the HF potential using the repulsive Buckingham function [80]

R

HF R Ae

V ( ) ^ (1) where R is the C-C distance, A and α (the repulsion exponent) are the fitting parameters. The dependence of the repulsion exponent on the basis size is shown in Table 1.2-1. It is seen that the repulsion exponent converges quickly after the 6-31G**

basis set being used.

B. MP2 calculations

Unlike the HF potentials, the MP2 potentials shown in Fig. 1.2-4 display clear minima and long range attractive potential tails. Because the contributions from the octopole-octople interactions are small, the dispersion energy is mainly responsible for the attractions. Sharp differences between the HF calculations and the MP2 calculations indicate the importance of including the correlation corrections in the wave function calculations. The HF method in principle does not include the correlation effect so the attraction forces are due exclusively to the correlation effect.

In Fig. 1.2-4, we compare the MP2 potentials with and without the BSSE corrections (denoted as CP and NCP, respectively). We see very strong dependence of the interaction potentials on the BSSE corrections. The potentials without the BSSE corrections fluctuate with increasing basis size and do not systematically converge to the destined curve at the basis set limit. On the contrary, the BSSE corrected potentials systematically approach to the destined curve with increasing basis size.

Therefore, it is important to consider the BSSE correction in calculating the intermolecular interactions, in particular for small basis sets.

As shown in Table 1.2-1, the basis set effect on the BSSE corrected interaction

potentials is significant. The STO-3G basis set yields very small binding energy. The interaction energy becomes more accurate as one adds polarization functions and augments diffuse functions in the Pople’s basis sets. Small cc-pVDZ and cc-pVTZ basis sets lead to underestimated binding energies and it requires cc-pVQZ to saturate the result. Augmentation of the diffuse functions has significant effect on optimizing the binding energy. The cc-pVTZ basis set underestimates the energy by 30%, while the aug-cc-pVTZ basis set underestimates only 5% of the binding energy. Some subtle basis set features can also be observed. For small basis sets, adding polarization functions to the basis set does not significantly change the potential. On the other hand, augmentation of the diffuse functions has pretty significant effect. For example,

potentials is significant. The STO-3G basis set yields very small binding energy. The interaction energy becomes more accurate as one adds polarization functions and augments diffuse functions in the Pople’s basis sets. Small cc-pVDZ and cc-pVTZ basis sets lead to underestimated binding energies and it requires cc-pVQZ to saturate the result. Augmentation of the diffuse functions has significant effect on optimizing the binding energy. The cc-pVTZ basis set underestimates the energy by 30%, while the aug-cc-pVTZ basis set underestimates only 5% of the binding energy. Some subtle basis set features can also be observed. For small basis sets, adding polarization functions to the basis set does not significantly change the potential. On the other hand, augmentation of the diffuse functions has pretty significant effect. For example,