Energy balance in a molecular dynamics simulation with ghost atoms and periodic boundary conditions?

Question

I am performing a simulation in LAMMPS with granules (DEM) and I had a question regarding the conservation of energy. If an interaction happens where one granule is in the simulation box and the other atom is on the other side (i.e it is a ghost atom), will its contribution to the potential energy be approximately halved?

I am guessing this because Newton's third law is not used for the ghost atoms? So even though the force is computed, it is only responsible for only half the resulting displacement? I see an overprediction in my total energy when I sum the constituent energies... It becomes greater than the total work done and I was wondering if this might be the cause?

Edit: I am using the default hookean potential $f=kd$ therefore the potential energy is $\frac{1}{2}kd^2$ where $k$ is a constant and $d$ is the overlap. I am assuming there isn't any friction right now. The total potential energy is therefore just summing this energy over all contacts and the kinetic energy is $\frac{1}{2}mv^2$.

The total work done is the energy input into the granular system: $$\int F_{piston}\cdot v_{piston} dt$$ Where the force term is the force exerted on the piston and the velocity is the velocity of the piston. The energy input should equal the summation of the kinetic and potential energies of the system.

The reason why I said mentioned Newton's 3rd law is because of

if (newton_pair || j < nlocal) {
          f[j][0] -= fx;
          f[j][1] -= fy;
          f[j][2] -= fz;
          torque[j][0] -= radj*tor1;
          torque[j][1] -= radj*tor2;
          torque[j][2] -= radj*tor3;
        }

in the source code. If I am interpreting this correctly, ghost atoms won't feel the f[j][:] force as they have j index greater than nlocal. I am running a serial simulation btw.

Although there exists a closed-form solution for the energy of a contact, if one looks at the incremental work done at a contact, it is given by the expression in (Asmar et al 2001 (Energy monitoring in distinct element models of particle systems):

If this is the case, shouldn't the j term be excluded in the above equation if one of the atoms is a ghost?

Answer 1

Newton's Third Law (N3L) is respected by force pairs in practically all molecular dynamics (MD) simulations, so where packages include a newton option (as LAMMPS does), it refers to the use of N3L in computing forces in parallel processes.

MD packages often parallelize calculations by assigning sets of atoms to different processes. They then have to decide how to compute forces for some pair ij where atoms i and j are "owned" by different processes ("local to", in LAMMPS code terminology) -- we'll call them processes i and j respectively. Since N3L is obeyed, the following two choices are equivalent:

1. Compute f_ij on (say) process i and send f_ij to process j, which can negate it and apply it to particle j. (LAMMPS's newton on)
2. Compute f_ij on process i and compute f_ji on process j. (LAMMPS's newton off)

The difference is that one send in choice 1 is replaced by one compute in choice 2, and thus which one is faster is determined by the exact problem and exact machine one is using.

For example, when simulating a small system on just one machine, sending information between processes is very cheap, so choice 1 is preferable. At the other extreme, some simulations contain millions of atoms spread across hundreds of machines, and all of that sending is much slower than just computing all forces locally on processes that need them (choice 2).

In either case, the final total force is the same (up to numerical errors). So when you say:

If I am interpreting this correctly, ghost atoms won't feel the f[j][:] force as they have j index greater than nlocal. I am running a serial simulation btw.

that's not correct. Later on in the LAMMPS code (including possibly outside the particular pair class) there's a "reverse communication" which sends those forces back from the ghost to the local atom.