Neighbor Tables for Short-Ranged Potentials

We have already discussed how to cutoff short-ranged potentials at the box edge, necessitated by the need to have a continuous force for molecular dynamics. But what happens as the number of particles grows? Clearly the code fragment the operation count and hence the CPU scales as N ², because the number of pairs of to be summed is N(N-1)/2. . This is its "complexity". (The cutoff actually reduces the required number of pairs to be summed to get force in box of side length L by 4pi(r_c/L)³.)

But as the system gets bigger we do not have to keep increasing the cutoff radius. Once the cutoff radius gets outside the second neighbor distance, or more rigorously we can fix r_c by the requirement that V(r_c) < kT. If this is satisfied, we can use a more efficient scheme to reduce the complexity to N.

• The neighbor list (Verlet scheme). Here we go through the particles and identify which are within an expanded cutoff radius equal to rc plus a "skin". Those are stored in an array linked to each particle. Then in computing the pair interaction and force, only pairs on the list are examined. We need to refresh the list when the two largest particle displacements sum to be larger than the skin depth. One chooses the skin depth to minimize the CPU time/step. If it is too large the neighbor lists get excessive. If it is too small, one continually has to refresh the list. The total CPU time will be proportional to N ² for building the list and N for evaluating the potential and forces. Clearly if the particles do not diffuse such as in a crystal or glass, one will only have to rarely or never refresh the neighbor table.
• One can speed up the generation of the table by use of the linked cell method. This works best in periodic boundary conditions. The whole cell is divided into a Cartesian product of small cells. Particles are then put into the cells (using a linked list for example.) To compute the force or potential on a given particle, first one computes the label of its cell. Also the names of the cells within the cutoff radius are identified. The neighbor list is constructed from that and we proceed as in 1). Now the construction of the neighbor list is order (N). What we need to optimize is the number of cells. Too small and one will have to loop over an excessive number of cells to find the few ones that have particles; too large and many extra particles will be on the neighbor list. This is the method that is used in CLAMPS. (see output)

Long-Ranged Potentials

Long-ranged potentials can not be handled as we discussed with the Lennard-Jones system. By long-ranged we mean that the limit ( v(r) r ³ ) does not tend to zero at large r. (Why is this important? Consider trying to make the box bigger and bigger. How quickly does the potential have to drop off before we can simply ignore distant particles?) Surprisingly we can still treat these types of interactions in PBCs and in some cases with the same complexity as for short-ranged potentials.

Suppose we have a potential between pairs of particles in vacuo which goes like v(r). There are two distinct issues which sometimes get confused:

• How should we define the potential in PBCs?
• How do we go about doing it efficiently?

V_image (r) = sum over array of images v(r+image) - background.

A particle interacting with another particle will see an infinite array of the other particle. One has to decide if this is physically correct. Inhomogeneous systems like a single charged impurity are problematical. With the image potential, one is essentially computing the interaction at a low but finite concentration. The background is put in to make the sum converge. If one sums spherical shells, the series will be conditionally convergent. However one has a choice about surface charges. The Ewald method has the great advantage that it respects the Poisson equation. This implies for example that the long-wavelength limits are correct.

To answer the second question see the more detailed description . The image expansion is implemented by splitting the interaction into a part in real space and a part in k-space in such a way that both sums are quickly convergent. The complexity of the method described is N ^3/2 worse than the short-ranged potential but still better than the simplest direct method for any pairwise sum. The coefficient isn't so bad because the new features are Fourier transforms which are fast.

What tests can you make for your Ewald code? You can put the charges on a perfect lattice and make sure that you get the known Madelung constants. For example in a 3D cubic lattice with lattice spacing "a" the energy per atom should be exactly -1.41864874 e²/a .

(See exercises in CLAMPS)

In the last 10 years methods have been developed to do this computation quicker. Some of the methods go by the name Fast Multipole Methods, others are called particle-mesh methods. They are advantageous if the number of particles is greater than about 1000.