A double exponential potential for van der Waals interaction

Van der Waals (vdw) interaction is an important force between atoms and molecules. Many potential functions have been proposed to model vdw interaction such as the Lennard-Jones (L-J) potential. To overcome certain drawbacks of existing function forms, this work proposes a double exponential (DE) potential that contains a repulsive exponential term and an attractive exponential term. This potential decays faster than the L-J potential and has a soft core. The DE potential is very flexible and its two exponential parameters can be adjusted continuously to mimic many potential functions. Combined with the isotropic periodic sum (IPS) method, the DE potential can efficiently and accurately describe non-bonded interactions and is convenient for alchemical free energy calculation.

The function have a minimum of −ε 0 at r = rm. σ is the distance where the potential is zero. The L-J potential has only two parameters and is therefore limited in how accurately it can be fitted to the properties of any real material. The sixth-power term in the L-J potential models effectively the dipole-dipole interactions due to electron dispersion in noble gases (London dispersion forces), but it does not represent other kinds of interactions well. The twelfth-power term appearing in the potential is chosen for its ease of calculation for simulations (by squaring the sixth-power term) and is not theoretically or physically based.
Several other forms have been used to describe vdw interaction, 2 such as the Buckingham potential: 3 which keeps the sixth-power term of the L-J potential but uses an exponential term to describe the repulsion. Study from the interacting quantum atoms (IQA) method has demonstrated that exponential relationships can better capture atomistic short-range repulsion. 4 It has been found that the buffered 14−7 potential can better describe vdw interaction. 2 The AMEOBA force field uses the buffered 14−7 potential and has found success in many applications. 5 The L-J potential diverges when two atoms approach one another, which is also a problem for the 1 r 6 term of the Buckingham potential. It is found the dispersion interaction can be better described with a damped 1 r 6 term, indicating a new way to construct the dispersion component is needed for a new-generation force field. 6,7 With above considerations, we propose a double exponential (DE) potential of the following form:

ARTICLE scitation.org/journal/adv
The DE potential has a repulsive part in the form of exp(−α r r m ) and an attraction part in the form of exp(−β r r m ). The exponential parameters, α and β, define the steepness of the repulsive interaction and the decay of the attraction, respectively. Like the L-J potential, the DE potential also has a minimum of −ε 0 at r = rm.
The DE potential has many advantages. One advantage is that the DE potential has a soft core, ε(0) = ε 0 βe α −αe β α−β , which make it convenient for free energy calculation through chemical perturbation. Another advantage of DE is that it converges faster than the L-J potential because the exponential functions decay faster than power functions. Third, the exponential parameters, α and β, can be When na = 12, nc = 6 eq. (4) reduces to the L-J potential, eq. (1).
For efficient calculation of long-range nonbonded interactions, we have developed the isotropic periodic sum (IPS) method. 9 where rc is the radius of the IPS local region, which is also called the cutoff distance. The IPS potential power m can be any integer. We pick m = 10 and 4 for n = 12 and 6, respectively, because they make eq. (5) best fit the analytically solved IPS potentials. 8 Therefore, the L-J IPS potential has the following form: ε IPS LJ (r, rc) = ε 0 r 12 m 1 r 12 + φ IPS p (r, rc, 12, 10) − 2r 6 m 1 Similarly, to efficiently handle long-range DE interaction, we can derive a simple IPS potential for the exponential function, e −κr , based on the homogeneity condition: 8 The DE IPS potential has the following form: One major concern about the DE potential, eq. (3), is the expensive calculation of the two exponential functions, exp(−α r r m ) and exp(−β r r m ). This burden has been greatly reduced by the fast approximation functions proposed to replace exponential function calculation. 10 On GPU, CUDA's implementation of the singleprecision exponential function, expf(), makes use of the fast hardware approximation and significantly speed up the calculation. The DE IPS potential, eq. (8), does not increase the calculation cost very much as compared with the DE potential, eq. (3), because the first exponential function in eq. (7) is related to a reciprocal of the exponential function in the DE potential and the other two exponential functions in eq. (7) are constants that can be precalculated. Therefore, the DE IPS potential, eq. (8), has only two exponentials to calculate for each pair of particles, the same as the DE potential, eq. (3).
The DE potential can mimic many types of repulsive-attractive potentials like Lennard-Jones potential. As shown in Table I, the DE potential that mimic the L-J potentials has α = 17.470 and β = 4.099. For convenience, we use DE (α, β) to represent a DE potential with the given α, β values. Fig. 1 (bottom) 1, 1.0) shows an flat energy surface with an accessible core region around r=0. Small α, β values result in flat energy surface, which corresponds to a slow decay, and a large cutoff distance would be needed for the cutoff method to have converged ensemble averages. To avoid large cutoff distances, one can use the IPS method to obtain converged results with reasonable small cutoff distances. The middle panel of Fig. 1 shows their corresponding IPS potentials, eq. (6) and eq. (8), which are short ranged within the cutoff distance, here, rc = 3rm. The top panel shows the IPS pair potentials, ε IPS pair (r, rc) = ε IPS (r, rc) − ε IPS (rc, rc), which are continuous and smooth at the cutoff boundary.

ARTICLE scitation.org/journal/adv
To illustrate the application of the DE potential, a cubic box of 500 particles are simulated with ε 0 = 119.8K = 0.238kcal mol and rm = 3.822Å, which are used as the Lennard-Jones parameters for argon. The cubic box side is 28.53Å and the periodic boundary condition (PBC) is applied. 10ns NVT molecular dynamics simulations are performed for each data point. Cutoff distances of 6, 8, 10, 12, 15, and 20Å are used to examine simulation convergency. Simulations with a straight cutoff method where energy goes to zero when r > rc are performed to compare with the IPS method.
NVT simulations are performed at T=100K for the L-J and the DE (17.470, 4.099) systems, and T=10 5 K for the DE (1.1, 1.0) systems to keep the systems in a fluid state. The average potential energies at different cutoff distances are shown in Fig. 2. Clearly, the average energies of the three systems are different. Even though the DE (17.470, 4.099) potential mimics the L-J potential, they are different in many ways: core softness, long-range decay, and IPS potentials. They are similar around the well and have exactly the same well depth, −ε 0 , and well location, rm (see Fig. 1), but the mentioned differences will result in different average energies. Therefore, to use the DE potential for vdw interaction, the parameters, ε 0 and rm, need be reparametrized and should not be taken directly from the L-J parameters.
Comparing the cutoff and the IPS results we can see clearly the IPS method produces converged results with a cutoff distance as small as 8 Å. The cutoff method needs much larger cutoff distance to reach the converged result. For the DE (17.470, 4.099) system, the cutoff results reach the IPS converged result at 12 Å (Fig. 2, middle panel), while for the L-J system the cutoff results do not reach the IPS converged result until 20 Å (Fig. 2, top panel), demonstrating that the DE potential decays faster than the L-J potential.
When α, β are small, the potential decay is slow and larger cutoff distances are needed for the cutoff method to obtain converged results. As can be seen from the bottom panel of Fig. 2, for the DE (1.1,1.0) system, even with a cutoff distance of 20 Å, the cutoff result is still far away from the IPS converged value. Again, the IPS results is almost independent of the cutoff distance. Therefore, for DE with small α, β values, the IPS method has significant advantage over the cutoff method. The soft core of the DE potential can be clearly seen from the radial distribution functions (RDF) shown in Fig. 3. The RDF of the L-J system and the DE (17.470, 4.099) system are almost identical, indicating the DE potential is a good mimic of the L-J potential. For the DE (1.1, 1.0) system, we can see all distance become accessible, indicating there is no hard core to stop particles from approaching. This is an important advantage of the DE potential that makes it suitable for alchemical free energy calculation.
With the four parameters, α, β, ε 0 , and rm, the DE potential is flexible enough to better fit into experiment properties of real matter. We will examine the DE potential in term of parameter fitting and property prediction in our future work. With all its advantages, we believe that the DE potential will find roles in new force field development and in molecular simulation.