LJ Potential

Empirical Lennard-Jones potential

Victor Anisimov

Introduction

The Lennard-Jones potential, hereinafter LJ, provides simple but effective means to model dispersion interaction between atoms. The importance of this function is in counter measuring the unsurpassed difficulty of first principle description of dispersion interactions in atomic systems which requires minimum CCSDT level of quantum mechanical theory and aug-cc-pVQZ basis set, the level which is not feasible for the systems of practical interest. On the other hand, the use of simple empirical LJ function, assuming the availability of highly optimized parameters, can provide even more reliable treatment of dispersion interactions than it can be obtained by using the respectful MP2 level of theory.

Similar to the popular non-empirical methods the modern NDDO-based semiempirical methods have limited ability to describe dispersion interactions. However, due to their parametric nature this limitation can be effectively fixed by adding empirical correction in the form of LJ potential.

LocalSCF program includes empirical LJ correction to the QM potential energy. The LJ term can be invoked with keyword "UseLJ". Following is the description of the program implementation.

Analytic expression of the LJ potential

The empirical LJ function, V_ij, describing the energy of dispersion interaction between atoms i and j separated by distance r, has the following form:

(1)

Here, the weight parameter w_ij is set to zero if the atoms i and j form a covalent bond (the so called 1,2-contact) or a valence angle (1,3-contact). For the atoms composing a dihedral angle (1,4-contact), the weight parameter varies in the range of [0.4 – 1.0]. For all other cases the weigh parameter is set to unit.

In our work we are going to use the CHARMM force field and its definition of the LJ function, which is slightly different from the above form within a freedom of choosing a linear constant.

(2)

In the CHARMM force field the available LJ parameters support proteins and DNA, which is minimally necessary for biomolecular simulations. The necessary parameters can be downloaded from the Web-site of Prof. Alexander D. MacKerell http://mackerell.umaryland.edu/CHARMM_ff_params.html

The LJ interaction energy has the following graphical form

Fig. 1. LJ interaction energy in diatomic system.

Here, r is the distance between atoms i and j, e_ij is well depth, and R_ij is generalized atomic radius computed from the radii of individual atoms:

(3)

In the CHARMM force field the division by 2 is implicitly included into the atomic radii R. The generalized well depth e_ij is computer from combination rule using the well depth parameters assigned to individual atoms

(4)

Correspondingly, each atom in the LJ potential is defined by two parameters, R_ii and e_ii. In the CHARMM force field the number of parameters is actually doubled; one set of parameters is used in the case of non-bond interactions, and the other set serves the case of 1,4-interactions. Because of that the weight parameter, w_ij, is set to zero for 1,2- and 1,3-contacts, and set to unit for all other cases.

Switching off the LJ potential

Since LJ energy calculation scales quadratically with the number of atoms in the system, direct calculation of all pair interactions for tens of thousand of atoms could easily become a computational bottleneck. Because of inverse power of 6 in the attractive part of the LJ function the interaction energy is expected to rapidly decrease with the distance between the particles. Therefore in order to avoid quadratic scaling, it is reasonable to cut-off the dispersion interactions starting from some distance. However, due to favorable character of the long-range dispersion interactions the cut-off distance should not be applied too early. In a large system encountering tens of thousands of atoms the dispersion interactions would quickly add up. For instance, in liquid water droplet having 60 Å in diameter and consisting of 4013 water molecules, the missing energy term would amount to 244 kcal/mol when treated with cut-off distance of 8 Å. If the LJ interactions are truncated at 6 Å, the interaction energy would be underestimated by 640 kcal/mol. The omission of such large energy value may impair the ability of the computational model to provide computational predictions with chemical accuracy. Therefore for large systems setting the LJ cut-off to distance less than 8 Å should be avoided. Suggested cut-off distance for the present system is 10 Å, where the energy difference over quadratic scaling LJ calculation is reduced to 113 kcal/mol.

Another problem, which requires serious consideration when dealing with systems encountering many thousands of atoms, is the property of the LJ energy in the vicinity of the cut-off distance. Since the value of the LJ energy drops to zero right after the cut-off distance the first derivative of the LJ function does not exist in that point. Such first derivatives are incompatible with molecular dynamics simulation. Therefore a common solution to avoid this problem is to incorporate a switching function which smoothly changes its value from 1 to 0 on a specified interval. Such property is provided by a polynomial function

s(r) = 0

r < r_on

s(r) = 1

r > r_off

Here we turn on the switching function at r_on= 8 Å, which starts gently approaching zero, and switches the LJ potential completely off at r_off = 10 Å. Unlike the original problematic case, the resulting LJ function has smooth first derivatives. In this example the switching occurs in the interval of 2 Å. Increasing the interval would improve the derivatives, though it may not be entirely necessary because of small magnitude of the LJ first derivatives. We will return to discussion of this issue after obtaining analytic derivatives of the LJ potential.

Figure 2 illustrates the work of switching function for dispersion interaction between two oxygen atoms in water dimer. Energy units are kcal/mol, the units of distance are Å. The LJ energy up to 8 Å is computed according to Eq. 2.

Fig. 2. Comparison of the original and switched (on distance from 8 to 10 Å) LJ interaction energy between two oxygen atoms.

Analytic first derivative of the original LJ function in the single-dimensional case of two-atomic system on the interval [0;8[ Å has the following form:

The LJ derivative for the switched interval [8;10] Å is

where

For a 3-dimensional system in Cartesian coordinates the differentiation of LJ potential can be done in the following way

The distance between interacting atoms is defined according to the following expression

Gradient of this function in respect to coordinates of atoms i and j is following

Correspondingly, the gradient of the LJ potential in Cartesian coordinates is

Now when we have the tool to calculate derivatives of the LJ function we can inspect the behavior of LJ gradient on the switched interval.

Fig. 3. First derivative of LJ potential on switched interval for water dimer (O-O interaction).

In Figure 3 we examined three choices of the cut-off interval. As one can see, choosing small switching interval introduces hump in the shape of gradients. Though, it is not very pronounced in absolute value on the 2 Å interval (cases 08-10 and 09-11). The general rule is the later one applies switching function the smaller the hump becomes (compare 09-11 vs. 08-10). However, the best solution is to extend the switching interval to 3 Å, see the case 08-11 where the hump is the smallest.

The visual inspection of LJ gradients shows a break at the point of application of switching function. It means the LJ gradient is not a smooth function now, i.e. the LJ second derivative will not be defined in that point. Luckily it has a minor numerical implication. Unless we are concerned with vibrational frequencies, which require second derivatives, the highest derivatives needed in typical molecular simulations, e.g. molecular dynamics, are of first order. Since we are dealing with the LJ potential in the QM framework, in order to improve accuracy of first derivatives in presence of switching function, we can afford applying very large cut-offs and still the LJ computation will be only a tiny fraction of the total computation time.

In semiempirical LocalSCF calculations the use of LJ correction proved its utility in calculation of protein-ligand absolute binding free energy, particularly when dealing with large and hydrophobic ligands.