Simulations involving the Lennard-Jones potential usually employ a cutoff at *r* = 2.5σ. This communication investigates the possibility of reducing the cutoff. Two different cutoff implementations are compared, the standard shifted potential cutoff and the less commonly used shifted forces cutoff. The first has correct forces below the cutoff, whereas the shifted forces cutoff modifies Newton's equations at all distances. The latter is nevertheless superior; we find that for most purposes realistic simulations may be obtained using a shifted forces cutoff at *r* = 1.5σ, even though the pair force is here 30 times larger than at *r* = 2.5σ.

Molecular dynamics (MD) simulations solve Newton's equations of motion by discretizing the time coordinate. The time-consuming part of a simulation is the force calculation. For a system of *N* particles this is an *O*(*N*^{2}) process whenever all particles interact. In practice the interactions are negligible at long distances, however, and for this reason one always introduces a cutoff at some interparticle distance *r* = *r*_{c} beyond which interactions are ignored.^{1}

The standard Lennard-Jones (LJ) pair potential is given by

Usually, a cutoff at *r*_{c} = 2.5σ is employed; at this point the potential is merely 1.6% of its value at the minimum (−ε). Although a cutoff makes the force calculation an *O*(*N*) process, this calculation remains the most demanding in terms of computer time.

The present communication investigates the possibility of reducing the LJ cutoff below 2.5σ without compromising accuracy to any significant extent. Before presenting evidence that this is possible, it is important to recall that quantities depending explicitly on the free energy are generally quite sensitive to how large is the cutoff. Examples include the location of the critical point,^{2} the surface tension,^{2,3} and the solid–liquid coexistence line.^{4,5} For such quantities even a cutoff at 2.5σ gives inaccurate results, and in some cases the cutoff must be larger than 6σ to get reliable results.^{3} Note, however, that if a simulation gives virtually correct particle distribution, the thermodynamics can be accurately calculated by first-order perturbation theory.^{6}

We compared two cutoff implementations at varying cutoffs with the “true” LJ system, the latter being defined here by the cutoff *r*_{c} = 4.5σ. One cutoff is the standard “truncated and shifted potential” (SP for shifted potential), for which the radial force is given^{1} by [*f*_{LJ}(*r*) = −*u*′_{LJ}(*r*) is the LJ radial force]

This is referred to as a SP cutoff because it corresponds to shifting the potential below the cutoff and putting it to zero above, which ensures continuity of the potential at *r*_{c} and avoids an infinite force here.

The “truncated and shifted forces” cutoff (SF for shifted forces)^{1,7} has the force go continuously to zero at *r*_{c}, which is obtained by subtracting a constant term:

This corresponds to the following modification of the potential: *u*_{SF}(*r*) = *u*_{LJ}(*r*) − (*r* − *r*_{c})*u*′_{LJ}(*r*_{c}) − *u*_{LJ}(*r*_{c}) for

*u*

_{SF}(

*r*) = 0 for

^{8}

We simulated the standard single-component LJ liquid at the state point that in dimensionless units has density ρ = 0.85 and temperature *T* = 1.0.^{9} This is a typical moderate-pressure liquid state point.^{1,10} Other state points were also examined—including state points of the fcc crystal, at the liquid–gas interface, at the solid–liquid interface, and for a supercooled system—leading in all cases to the same overall conclusion. For this reason we report below results for just one state point of the LJ liquid and one of the Kob–Andersen binary LJ (KABLJ) liquid.^{11} 2000 LJ particles were simulated using the standard central-difference constant temperature/energy (*NVT* /*NVE*) algorithms (Figs. 2–4 and 6, respectively); 1000 particles of the KABLJ liquid were simulated using the *NVT* algorithm (Fig. 5).

Figure 1 shows the basics of the LJ system. In the upper figure the black curve gives the LJ pair potential *u*_{LJ}(*r*) and the black dashed curve the radial distribution function *g*(*r*), which has its maximum close to *u*'s minimum. In the lower figure the black (lower) curve shows the LJ pair force *f*_{LJ}(*r*). The red (upper) curve gives *f*_{SF}(*r*) when a cutoff at 1.5σ is introduced; note that the shifted force differs significantly from the true force.

Figure 2 shows the true pair-distribution function (black) and the simulated *g*(*r*) for three *r*_{c} = 1.5σ cutoffs: SF (red, barely visible), SP (green, slightly lower at *r* = 1.5σ), and a smoothed SP cutoff ensuring the force and its first derivative go continuously to zero at the cutoff^{12} [dashed (green) curve]. The curves deviate little, except near the cutoff where the smallest errors are found for a SF cutoff (inset).

In order to systematically compare the SP and SF cutoffs we studied the LJ liquid for a range of cutoffs. Figure 3 quantifies the difference between the computed *g*(*r*) and the true radial distribution function, *g*_{0}(*r*), by evaluating

*r*

_{c}above the Weeks–Chandler–Andersen (WCA) cutoff at the potential energy minimum

^{6}where SF = SP (

*r*

_{c}= 2

^{1/6}σ = 1.12σ). Smoothing a SP cutoff has only a marginal effect compared to not smoothing it (results not shown). Applying first-order perturbation theory with the

*g*(

*r*) obtained in a simulation with SF cutoff at

*r*

_{c}= 1.5σ leads to a pressure that deviates only 1% from the correct value.

Figure 4 studies energy drift in long *NVE* simulations for *r*_{c} = 1.5σ. The SF cutoff (red, horizontal) exhibits no energy drift, whereas SP (green, diverging) does. Figure 4 also gives results when the force of a SP cutoff is smoothed^{12} (green dashed, horizontal curve). This leads to much better energy conservation,^{1} but the energy fluctuations are somewhat larger than for a SF cutoff. The simulations indicate the existence of a hidden invariance in the central-difference algorithm for a continuous force function, deriving from a “shadow Hamiltonian.”^{13}

Not only static quantities, but also the dynamics are affected little by replacing a 2.5σ SP cutoff with a 1.5σ SF cutoff. This is demonstrated in Fig. 5, which shows simulations of the incoherent intermediate scattering function of the supercooled KABLJ liquid.^{11} For reference a WCA cutoff simulation is included (blue dashed curve, fastest relaxation), which was recently shown to be inaccurate despite the fact that the WCA radial distribution function is reasonably good for this system.^{14} A SP cutoff at *r*_{c} = 1.5σ_{AA} gives too slow dynamics (purple dotted curve). Within the numerical uncertainties incoherent scattering functions are identical for the “true” system, a SP cutoff at *r*_{c} = 2.5σ_{AA}, and a SF cutoff at *r*_{c} = 1.5σ_{AA}. Similar results were found for the single-component LJ liquid's dynamics. We conclude that a SF cutoff at *r*_{c} = 1.5σ generally works well for both statics and dynamics of LJ systems.

Why does a cutoff, for which the forces are modified at all distances (SF), work better than when the forces are correct below the cutoff (SP)? A SF cutoff modifies the pair force by adding a constant force for all distances below *r*_{c}; at the same time SF ensures that the pair force goes continuously to zero at *r* = *r*_{c}. Apparently, ensuring continuity of the force—and thereby that *u*″(*r*) does not spike artificially at the cutoff—is more important than maintaining the correct pair force below the cutoff. How large is the change induced by the added constant force of the SF cutoff? Figure 6 shows the x component of the force on a typical particle as a function of time (*r*_{c} = 1.5σ). The black curve gives the true force, the red, barely visible curve the SF force, and the blue, fluctuating curve the SF correction term. Although the true and SF individual pair forces differ significantly (Fig. 1), the difference between true and SF total forces is small and stochastic (∼3%). This reflects an almost cancellation of the correction terms deriving from the fact that the nearest neighbors are more or less uniformly spread around the particle in question. It was recently discussed why adding a linear term (∝*r*) to a pair potential hardly affects dynamics^{15} and statistical mechanics:^{16} For a given particle's interactions with its neighbors the linear terms sum to almost a constant because, if the particle is moved, some nearest-neighbor distances increase and others decrease, and their sum is almost unaffected.

Figure 6(b) shows details of Fig. 6(a); we here added the SP force for the same cutoff (green, upper curve). Both SP and correction terms are discontinuous; they jump whenever a particle pair distance passes the cutoff. Altogether, Fig. 6 shows that not only does the sum of the constant forces on a given particle from its neighbors cancel to a high degree, so do the interactions with particles beyond the cutoff. The result is that the particle distribution is affected little by the long-range attractive forces, a fact that lies behind the success of perturbation theory.^{6,17,18}

In summary, when a SF cutoff is used instead of the standard SP cutoff, errors are significantly reduced. Our simulations suggest that a SF cutoff at 1.5σ may be used whenever the standard SP cutoff at 2.5σ gives reliable results; this applies even though the pair force at *r* = 1.5σ is 30 times larger than at *r* = 2.5σ. A cutoff at 1.5σ is large enough to ensure that all interactions within the first coordination shell are taken into account (Fig. 2). Use of a 1.5σ SF cutoff instead of a SP cutoff at 2.5σ leads potentially to a factor of (2.5/1.5)^{3} = 4.7 shorter simulation time for LJ systems.

The center for viscous liquid dynamics “Glass and Time” is sponsored by the Danish National Research Foundation (DNRF).

## REFERENCES

*f*

_{SP}(

*r*) was smoothed in the interval

*r*

_{c}± δ by replacing it with the function

*f*(

*r*) =

*f*(

*r*

_{c}− δ)

*h*(

*x*) with

*h*(

*x*) = 1 +

*cx*− (3 + 2

*c*)

*x*

^{2}+ (2 +

*c*)

*x*

^{3}, where

*x*= (

*r*−

*r*

_{c}+ δ)/2δ and

*c*= 2δ

*u*″(

*r*

_{c}− δ)/

*u*′(

*r*

_{c}− δ). This ensures that

*f*

_{SP}and

*f*′

_{SP}go smoothly to zero in the cut interval. In the simulations δ = 0.025 (other values lead to similar conclusions).