X. W. Sheng^{a} and
K. T. Tang*^{b}
^{a}Department of Physics, Anhui Normal University, Anhui, Wuhu 24100, China
^{b}Department of Physics, Pacific Lutheran University, Tacoma, Washington 98447, USA. E-mail: tangka@plu.edu
First published on 9th September 2020
A chronological account is given to the development of a full range interatomic potential. Starting with a simple phenomenological model, the terms in the model are gradually modified, so that they can carry some definite physical meaning. To gain insight, a systematic, order by order interaction potential theory is developed. Conversely, this theory suggests the functional form for the potential model. At present, we have a simple interaction model that is capable of describing the van der Waals potentials of many systems from R = 0 to R → ∞.
If, in some cataclysm, all of scientific knowledge were to be destroyed, and only one sentence passed on to the next generation of creatures, what statement would contain the most information in the fewest words? I believe it is the atomic hypothesis that all things are made of atoms, attracting each other when they are a little distance apart, but repelling upon being squeezed into each other.
In the second part of that one sentence, he is talking about interatomic potential.
Almost all physical phenomena, outside the atomic nucleus, can be attributed directly or indirectly to the forces between atoms. Basic concepts such as temperature, pressure, viscosity of liquid and strength of solid are all intimately related to the interatomic potentials.
With modern quantum chemical methods and the advent of powerful computers, the interatomic potentials can in principle be calculated. Indeed for some rare gas systems, such as He_{2} and Ar_{2}, very accurate ab initio calculations have been carried out.^{2–6} Such calculations are very important in that they provide the benchmark references. However, such calculations are very complicated and expensive and potentials are obtained in numerical form only. Furthermore, in the highly repulsive region, quantum chemical calculations still have difficulties.^{7} Therefore there is a long standing interest in simple yet reliable analytic models of interaction potentials, which should ideally consist of physically identifiable terms and which should be based on a limited number of available theoretical or empirical data.
The foundation of the theory of interatomic potentials was laid in the year of 1927. In that year the chemical bond of the singlet state of H_{2} molecule was successfully described by Heitler and London^{8} using a symmetrized wave function. The calculation showed that the stabilization of the bond is provided by the exchange integral. The same calculation also demonstrated that the exchanged integral is responsible for the repulsive potential of the triplet state. In the same year Wang^{9} pointed out that there is always a long range attraction between any pair of atoms, which is now known as the dispersion energy, and is due to the interaction of a dipole with an instantaneous induced dipole. This long range attraction together with the exchange repulsion are the origin of the so called van der Waals potential between systems which do not form a chemical bond.
In this paper we will describe a simple model for van der Waals potentials which can accurately describe the full range potentials of many systems. We will first provide a chronological account of the development of this model and then give some examples of its applications.
Atomic units will be used throughout this paper. To facilitate comparison with experiments other units will from time to time also be used. When this is the case, it will be explicitly stated.
In 1973, Toennies published a paper^{14} in which he pointed out that the Buckingham potential model modified by including the R^{−8} and R^{−10} dispersion terms
(1) |
This observation is the stimulus and the starting point of our journey of searching for physically transparent models of interatomic potentials.
It became clear, upon close examination of eqn (1) that the good agreement with experiment must be due to some fortuitous cancellation of errors.^{14} In particular, the C_{10} term in the equation was found to have a strong influence on well depth, contributing between 25–75% with the larger value applying to the heavier systems. While it is well known that the long range attractive potential can be described by the dispersion series obtained from the second order perturbation theory,^{15,16} but this is an asymptotic series and has more than three terms.^{17} This raises the question as to what effect the rest of the series will have. Unfortunately coefficients higher than C_{10} were not available at the time. Furthermore, the divergent nature of the asymptotic series was not accounted for. In other words, the dispersion series must be “damped” as R decreases.
Our first attempts to answer these questions were to use the fact that the dispersion series are asymptotic series. It is well known that the best estimate of an asymptotic series at any given R is obtained by terminating the series at the term before the smallest term and including a fraction of the smallest term.^{17} Thus, with increasing R the number of terms retained in the series should be increased. Based on this property, we were also able to derive a recurrence relation:^{18}
(2) |
With the knowledge of the first three dispersion coefficients, higher ones can be estimated within a few percent. Although much more complicated recurrence relations were subsequently derived,^{19} tests against the accurate coefficients of He–He interaction, which became available later,^{20} show that eqn (2) works just as well.
With these corrections, the Toennies' model was able to predict the potentials of a large number of atom–atom,^{21} atom–molecule^{22–24} interactions.
It should be mentioned that the well known HFD (Hartree–Fock Dispersion) potential was also based on the Toennies model, except the damping of the dispersion series is modeled after the H_{2} interaction.^{25,26} There are a great many different HFD potentials in literature. Either the repulsive functional form fitting the SCF data is somewhat different, or the attractive series is changed through twiddling the dispersion coefficients and/or the empirical damping function. For example, at least six different HFD potentials, from HFD-I to HFD-D, are mentioned in ref. 27 alone. The purpose of these efforts to modify the HFD potential is to fit all experimental results into the same potential. The results of these attempts are summarized in the Aziz's report.^{28}
The dispersion coefficients C_{2n} were derived by assuming that there is no overlap of the electron clouds of the two interacting atoms. This of course is not the case. We invoked the semiclassical Drude model,^{29} similar to the one used earlier by Kim and Gordon,^{30} to estimate the effect of the overlap of the electron clouds on the dispersion coefficients. It turned out that this effect can be written in such a way that the dispersion terms C_{2n}/R^{2n} are reduced to f_{2n}(bR)C_{2n}/R^{2n}, with f_{2n}(bR) given by
f_{2n}(bR) = 1 − P_{2n}(bR)exp(−bR) |
f_{2n}(bR) → 1, R → ∞ |
f_{2n}(bR) → 0 + 0(R^{2n+1}), R → 0, |
(3) |
The most discriminating test was the accurate ab initio calculations of Koide, Meath and Allnatt for H(1S)–H(1S).^{32} They presented the damping function for 2n = 6…(2)…20. The agreement with eqn (3) in all cases is remarkably good. The agreement is best for the most important terms with 2n = 6, 8 and 10. For the higher order terms, at internuclear distances where f_{2n} = 0.5, eqn (3) are larger by 2% for 2n = 12, and in the worst case by 7% for 2n = 20. Later, when accurate ab initio damping functions for some other systems became available,^{33,34} they are found to be all in good agreement with eqn (3).
Thus, with the Born–Mayer repulsion, we have the model potential^{31}
(4) |
Compare to the accurate potentials available for systems such as Ar_{2}, the well depth D_{e} (the potential minimum) obtained from the TT potential of eqn (4) with the repulsive term Aexp(−bR) fitted to SCF ab initio data, are about 15% too deep. This is consistent with the observation from the ab initio perturbation theory that the inclusion of the correlation effects increases the magnitude of repulsion. Since this and other small effects are difficult to determine and the well depth D_{e} and its location R_{e} are available for many systems, we propose to determine A and b from D_{e} and R_{e}. This is most conveniently done with the reduced potential. If we define x = R/R_{e}, then eqn (4) can be written as V(R) = D_{e}U(x) where U(x) is the reduced potential and is defined as
U(1) = −1, |
These two equations are sufficient for us to find A and b. Now
Solving for A* from the second equation and substituting it into the first equation, we have
(5) |
This equation can be solved for the only unknown b*. This can be most expediently done on a computer. In the appendix of ref. 35, there is a simple program to solve this equation for the particular case of He_{2} with N_{max} = 8. With b* we can easily calculate A*, so A = A*D_{e} and b = b*R_{e}.
In the original paper,^{31} we showed that with the parameters in eqn (4) so determined, it was possible to predict the interaction potentials for not only the rare gases, but also H^{2}(^{3}Σ), NaK(^{3}Σ), and LiHg complexes. Since then, this model has been successfully applied to many more systems.^{36–42}
(6) |
In this expression an error of a factor 4.5 in the result reported by Duman and Smirnov^{48} has been corrected. The one dimensional integral can be evaluated and D/A^{4} can be fitted to the form^{49}
D/A^{4} = 0.0129 + 0.1297β − 0.0403β^{2} |
ε_{x} = 0.818R^{2.5}e^{−2R}. |
This is the asymptotic exchange energy of H_{2} by Herring and Flicker.^{50} They were the first to obtain the correct asymptotic energy. Although a number of approximations were used by Herring^{51,52} to obtain the effect of electron–electron correlation, he was able to show that this is the asymptotic exact result. It is interesting to note that when Herring undertook this study, the purpose was to correct the conceptual defect of the Heitler–London method. Even Herring himself thought this approach is only “good for the soul”. Later it turned out that this is in excellent agreement with the most elaborate ab initio calculation^{53} not only in the van der Waals region but all the way down to R ≈ 2 a.u.
The surface integral method is based on the physical process of electrons hopping between nuclei. Although this interpretation is readily accepted in the charge exchange scattering experiments, it apparently contradicts textbook wisdom. While the physical interpretation might still be controversial, the mathematical formulation is indisputable.
The theory was further extended by Smirnov and Chibisov^{54} to treat general atomic systems. It is based on the fundamental approximation that the exchange interaction between two multielectron atoms is dominated by the exchange of a single pair of electrons at any one time. However, in order to antisymmetrize the total wavefunction, the coupling of orbital and spin angular momenta of all equivalent electrons has to be taken into account. As a result, the exchange energy of a multielectron system is found to equal to the exchange energy of single pair multiplied by a rather complicated coupling constant k which was first calculated by Duman and Smirnov.^{48} Unfortunately their final results were incorrect. The correct result is given in ref. 55. Although these k factors are the result of a rather complicated angular momentum coupling, all of them are identical to the number of possible exchanges between the valence electrons with the same spin in the two atoms. For example, the k values for H–He, He–He, He–Ne and Ar–Ar interactions are respectively 1, 2, 6 and 18.
According to the GHL, the van der Waals potential is given by V = E_{t} − ε_{0} and for M = 2,
V = ε_{1} + ε_{2} + (1 +s_{0} + s_{1})ε_{x} |
To use this theory to study the potential curves of He_{2}, Ne_{2} and Ar_{2},^{56} various terms were evaluated with CI wave functions. It turn out that ε_{2} is very well approximated by the damped dispersion series, and ε_{1} is almost parallel to ε_{x} and is only about 10% of ε_{x} in magnitude. The overlaps s_{0} and s_{1} are negligibly small in the van der Waals region. Thus by decrease A in the ε_{x} by about 2.5%, the ε_{1} contribution can be accounted for because exchange energy depends on A^{4}. With this modification, the van der Waals potential can be written as V = ε_{x} + ε_{2}.
In a parameter free study of the van der Waals potential of He_{2},^{57} we used the expression of the exchange energy for repulsion and for attraction we used the damped dispersion series with b given by the logarithmic derivative of the repulsion. Since the ionization energy of He is 0.9036, so β = 1.3443.The asymptotic amplitude A calculated from SCF is determined to be 2.84. Therefore D = 7.449. With these parameters and the known dispersion coefficients, we found a potential curve that is in amazing agreement with the no-appproximation bruteforce Monte Carlo calculations.^{58} This agreement is fortuitous. The CI wave function determined A is equal to 2.91^{56} which is 2.46% larger than 2.84 determined by the SCF wave function. Therefore the neglected first order Coulomb energy is completely compensated.
The expression for the heterogeneous molecules is far more complicated, and is currently an interesting subject for research.^{46,59} Furthermore, even for homonuclear dimers, the amplitude D is very difficult to determine. Therefore in our search for a simple potential model for both heterogeneous and homogeneous dimers, we used what we call the Tang–Toennies–Yiu (TTY) model^{60}
(7) |
Thus, the TT potential can easily be transformed into the TTY potential. We think it is an improvement because its repulsive potential has a stronger theoretical foundation. Furthermore, this functional form is softer in the intermediate distances and is in better agreement with experiment.^{61} The obvious disadvantage of TTY is that it turns to zero as R goes to zero.
E_{el}(0) = E_{ua} − E_{A} − E_{B} = C. |
The electronic energy must also satisfy the condition^{62}
Following Buckingham,^{63} we write the short range potential as
The final expression for the extended potential (TT2) is given by joining the short range to the TTY potential,^{60}
V_{TT2}(R) = V_{short}(R) + (1 − e^{−αR})V_{TTY}(R), | (8) |
In addition to Z_{A} and Z_{B} (which are known), there are four additional parameters a_{1}, a_{2}, a_{3} and α which can be determined in the following way. By expanding e^{−αR} in terms of power series, we see that
(9) |
Similarly, from
(10) |
Furthermore, we require that the new potential V_{TT2}(R) is consistent with the potential parameters D_{e} and R_{e} at the van der Waals minimum. That is
V_{TT2}(R_{e}) = D_{e}, |
The first condition leads to
The second condition leads to
Z_{A}Z_{B}(a_{1} + 2a_{2}R_{e} + 3a_{3}R_{e}^{2}) = D_{e}. |
Combining the last two equations, we have
Putting in a_{1} and a_{2} of eqn (9) and (10), we have a quadratic equation of α which can be solved to give
Thus all parameters are determined.
For H–H, H–He, and He–He interactions, the input data and the calculated parameters of eqn (8) are shown in Table 1. This table is similar to the Table 1 in ref. 60, except in the present case we include 3 more dispersion terms. Because of that, some of the calculated parameters are different. This will facilitate our discussion of the upper limit N_{max} of the dispersion series.
H–H | H–He | He–He | |
---|---|---|---|
Input parameters | |||
C_{6} | 6.49903 | 2.82134 | 1.46212 |
C_{8} | 124.399 | 41.8364 | 14.1258 |
C_{10} | 3285.83 | 871.540 | 183.781 |
C_{12} | 1.1977 × 10^{5} | 2.5507 × 10^{4} | 3.2199 × 10^{3} |
C_{14} | 6.0239 × 10^{6} | 1.0487 × 10^{6} | 7.5971 × 10^{4} |
C_{16} | 4.1811 × 10^{8} | 6.0575 × 10^{7} | 2.4138 × 10^{6} |
R_{e} | 7.83 | 6.66 | 5.608 |
D_{e} | −2.048 × 10^{−5} | −2.256 × 10^{−5} | −3.482 × 10^{−5} |
Z_{A} | 1 | 1 | 2 |
Z_{B} | 1 | 2 | 2 |
E_{A} | −0.49973 | −0.49973 | −2.90339 |
E_{B} | −0.49973 | −2.90339 | −2.90339 |
E_{ua} | −2.17503 | −7.47798 | −14.66844 |
Derived parameters | |||
TT parameters | |||
A | 10.51789 | 5.80043 | 25.18501 |
b | 1.68113 | 1.85898 | 2.41630 |
TTY parameters | |||
Λ | 0.74737 | 0.80930 | 7.94920 |
λ | 2.00034 | 2.18899 | 2.70024 |
TT2 parameters | |||
α_{1} | 0.90618 | 1.72603 | 1.84438 |
α_{2} | −0.28040 | −0.58597 | −0.75316 |
α_{3} | 0.01895 | 0.04568 | 0.06999 |
α | 2.08175 | 3.76346 | 4.05980 |
(11) |
In this potential, both U_{short}(x) and U_{long}(x) are dimensionless. Other than the Z_{A} and Z_{B} which are the well known number of nuclear charges, the potential depends only on the well depth D_{e} and its location R_{e} which should both be in atomic units. Note that Z_{A}Z_{B}/R_{e} is also energy in atomic unit. We call the potential of eqn (11), the TTS potential.
We will not repeat these calculations here, except in this paper we carry three more dispersion terms, in other words, in ref. 64 N_{max} = 5, and in this paper N_{max} = 8. This enables us to discuss the effects of more dispersion terms.
Theoretically, the damped dispersion series should include all the terms necessary for convergence. This usually means that a large number of terms must be included in the series to get the correct overall potential. However, if the repulsion parameters (A and b in the TT potential, Λ and λ in the TTY potential) are determined from R_{e} and D_{e} (like what we are doing here in solving eqn (5)) with only the first three leading dispersion terms, then this truncated model is usually adequate to predict an accurate potential in the van der Waals well region. This can be understood by noting that if R_{e} and D_{e} are correctly reproduced, then the small error caused by dropping terms higher than the third term in the dispersion series must be compensated by the adjustment of the repulsive potential near the minimum of the potential. Moreover, for R ≫ R_{e} the contribution from the neglected dispersion terms and the repulsive potential are all very small, and in this region the potential is dominated by the three leading dispersion terms. So there may be very little difference between the potential with three (N_{max} = 5) dispersion terms and the potential with six (N_{max} = 8) dispersion terms. On the other hand, for R ≪ R_{e}, since the compensation may not be complete in this region, and more dispersion terms contribute more to the negative component of the potential, therefore the potential with N_{max} = 8 may be softer than the potential with N_{max} = 5. Normally we should not make N_{max} much higher than 8, since higher dispersion coefficients are generally not available and the recurrence relation of eqn (2) is not precise. Repeated use of eqn (2) may produce unwarranted errors.
We illustrate this situation in Fig. 1 for the case of Ar_{2}. In Fig. 1a, R < 6 a.u., this is the region where the potential raises rapidly, therefore it is plotted in logarithmic scale. First we see how the accuracy is improved in going from TT potential to TTY potential, then from TTY potential to TTS potential. Secondly, we see that N_{max} = 8 curves are always smaller than the corresponding N_{max} = 5 curves. In Fig. 1b, R > 6 a.u. This is where the bowl of the van der Waals potential of Ar_{2} is located. In this region, the potential is plotted in linear scale. It is seen that in this region there is not much difference between these potentials. Therefore we can say that in the van der Waals region, the potential with N_{max} = 5 is a good enough approximation.
Fig. 1 The interaction potential of ground state Ar_{2}. The lines are the present result. Circles are the ab initio results taken from ref. 4 and 5. |
R_{e} | −D_{e} × 10^{5} | Z_{A} | Z_{B} | |
---|---|---|---|---|
a Taken from ref. 65 and 67.b Taken from ref. 66 and 67. | ||||
NeH | 6.425^{a} | 5.365^{a} | 10 | 1 |
ArH | 6.709^{b} | 17.052^{b} | 18 | 1 |
KrH | 6.898^{b} | 21.094^{b} | 36 | 1 |
XeH | 7.200^{b} | 25.357^{b} | 54 | 1 |
For the H–Ne system, we chose the results of Hishinuma^{65} who carefully measured the absolute integral cross sections between 1.8 and 330 meV. For H–Ar, H–Kr, and H–Xe systems, the orbiting resonance experiments of Toennies et al.^{66} still appear to be the most sensitive. These experiments were reanalysed with more accurate dispersion coefficients and the results are shown in ref. 67. The resulting potentials are shown in Fig. 2. Comparisons with the combining rule results of Tang and Toennies^{67} and with the MCSCF calculations of Das et al.^{68} show that these TTS potentials of eqn (11) are quite accurate.
Fig. 2 The interaction potential of the ground state rare-gas hydrides (Ne–H, Ar–H, Kr–H, Xe–H). The solid lines are the present results. Small circles are the combining rule results of Tang and Toennies.^{67} Crosses are the MCSCF results of Das et al.^{68} |
Fig. 3 The interaction potential of ground state HeLi. The red line is the present result. Circles are the ab initio results taken from ref. 69. |
Fig. 4 The interaction potential of ground state HeNa. The red line is the present result. Circles are the ab initio results taken from ref. 69. |
Fig. 5 The interaction potential of ground state HeK. The red line is the present result. Circles are the ab initio results taken from ref. 69. |
For Rb–He and Cs–He, accurate values of D_{e} and R_{e} are not known. We could estimate them using combining rules,^{35,70} but they probably can be more accurately calculated in the following way. The dispersion coefficients of alkali–helium listed in the table of Jiang el al.^{71} are rather accurate. Since these systems conform with each other, the reduced dispersion coefficients (C_{6}* = C_{6}/(R_{e}^{6}D_{e}), C_{8}* = C_{8}/(R_{e}^{8}D_{e}) should be nearly the same for all of them. These quantities are listed in Table 5.
R_{e} | −D_{e} × 10^{6} | C_{6} | C_{8} | C_{6}* | C_{8}* | |
---|---|---|---|---|---|---|
a Ref. 69.b Ref. 71.c Calculated in the present paper. | ||||||
HeLi | 11.47^{a} | 7.36^{a} | 22.535^{b} | 1084.2^{b} | 1.34462^{c} | 0.49173^{c} |
HeNa | 11.85^{a} | 6.96^{a} | 25.759^{b} | 1326.9^{b} | 1.33663^{c} | 0.49032^{c} |
HeK | 13.50^{a} | 5.00^{a} | 39.468^{b} | 2620.7^{b} | 1.30398^{c} | 0.47509^{c} |
It is seen that the reduced dispersion coefficients C_{6}* and C_{8}* are indeed nearly identical for all three interactions. Now it follows from the definitions of the reduced dispersion coefficients that
Using the average values of the reduced dispersion coefficients of C_{6}* and C_{8}*, and the C_{6} and C_{8} values from Jiang's Table, we calculate the R_{e} and D_{e} for all alkali–helium systems, their values are listed in Table 6. Together with Z_{A} and Z_{B}, The TTS potentials for these systems can be plotted.
C_{6} | C_{8} | R_{e} | −D_{e} × 10^{6} | Z_{A} | Z_{B} | |
---|---|---|---|---|---|---|
HeLi | 22.535 | 1084.2 | 11.47 | 7.45 | 2 | 2 |
HeNa | 25.759 | 1326.9 | 11.87 | 6.93 | 2 | 11 |
HeK | 39.468 | 2620.7 | 13.48 | 4.96 | 2 | 19 |
HeRb | 44.785 | 3144.9 | 13.86 | 4.76 | 2 | 33 |
CsHe | 52.679 | 4273.1 | 14.89 | 3.63 | 2 | 54 |
For Li–He, Na–He and K–He, they are almost indistinguishable from Fig. 3, 4, 5 respectively. Fig. 6 shows the potentials we predicted for Rb–He and Cs–He. The predicted R_{e} and D_{e} for Rb–He and Cs–He are compared with other determinations in Table 7. Clearly these predicted values are reasonable compared to some vastly more complicated calculations.
Fig. 6 The interaction potentials of the ground state He–Rb and He–Cs predicted by the TTS potential model with R_{e}, D_{e}, Z_{A}, Z_{B} listed in Table 4. |
RbHe | CsHe | |||
---|---|---|---|---|
R_{e} | −D_{e} × 10^{6} | R_{e} | −D_{e} × 10^{6} | |
Present | 13.86 | 4.76 | 14.89 | 3.63 |
Ayed et al.^{72} | 13.0 | 8.72 | 14.86 | 6.15 |
Pascale^{73} | 14.0 | 6.83 | 15.0 | 5.92 |
Baylis^{74} | 11.6 | 11.85 | 12.0 | 11.39 |
Kleinekathofer et al.^{49} | 14.24 | 3.71 | 15.02 | 3.17 |
Cvetko et al.^{75} | 13.86 | 4.92 | 14.61 | 4.14 |
Patil^{76} | 13.89 | 4.48 | 14.61 | 3.82 |
In the introduction, we mentioned that the foundation of the theory of interatomic potential was laid in 1927 by Heitler and London in their study of H_{2} molecule using symmetrized wavefunction. In the following year, 1928, Pauling^{77} applied the Heitler–London theory to study the H_{2}^{+} system. Although in this case, the approach is now commonly known as the linear combination of atomic orbitals, the formalism is the same. Thus, if the singlet state and triplet state wavefunctions of the H_{2} molecule are replaced by the gerade and ungerade state wavefunctions, respectively, of the H_{2}^{+} molecular ion, all the formulae derived for H_{2} are equally applicable to H_{2}^{+}. The specific results are of course different.
The most obvious difference is that for ion–atom systems, in addition to dispersion attraction, they also have induction attraction. Other than that, most of the present formalisms are equally applicable to ion–atom systems. For example, the alkali ion–rare gas and halogen ion–rare gas systems were studied with the TT potential model by Ahlrichs et al..^{78} It will be interesting to study these systems with TTY potential since the Herring type exchange terms for ion–atom systems are also availabe.^{79,80} In fact, we have used the Generalized Heitler–London theory to study the alkali dimer cations and obtained very good results.^{81,82}
Footnote |
† This paper is written for the celebration of the 90th birthday of Prof. J. P. Toennies who has been our inspiration and mentor for almost half a century. |
This journal is © the Owner Societies 2020 |