Discontinuity, anisotropy, and cutoffs

I am simulating a slab with a coarse-grained model based on ellipsoids (sorry, not yet in the official distribution), compiled on LAMMPS 21 Nov 2023. The model is based on bonded ellipsoids decorated with point charges: think of them as a rigid body composed of a spheroid plus massless point charges attached to it. The computation of the neighbour list is affected by the anisotropy of the ellipsoids: here I have a particle whose axes are: 13.5 5.8 3.5, so I typically use a skin distance of about 5 Å to avoid missing interactions with molecules outside the cutoff distance which may suddenly fall inside the cutoff, e.g. due to the rotation of the long ellipsoid.

I am using the CG model in a deposition protocol which goes like this:

  1. Use long-range electrostatics with slab correction.
  2. Insert a molecule with read_data in a random position.
  3. Thermalise the molecule and the slab with two different thermostats, both fix nvt/asphere.
  4. Add a velocity to the inserted molecule and integrate its trajectory with fix nve/asphere.
  5. Monitor the vertical position (Z coordinate) and total force acting on the centre-of-mass of the deposited molecule.

These are the relevant commands:

variable cutoff equal  14.
neighbor 5. bin
neigh_modify check yes
boundary p p f
# a CG pair_style is defined with the chosen ${cutoff}.
# ...
read_data npb_slab.data group base # previously equilibrated.
read_data npb_single.data add append group m1
displace_atoms m1 rotate 0 0 0 -0.715 0.311 0.730 82.4
displace_atoms m1 move 14.7 18.6 34.
variable zcom  equal xcm(m1,z)
variable vzcom equal vcm(m1,z)
variable fzcom equal fcm(m1,z)
fix DEP all ave/time 10 100 1000 v_zcom v_vzcom v_fzcom file output.dat mode scalar
# Thermalisation of the inserted molecule.
fix             NVT1  base nvt/asphere temp 300 300 1000
fix             NVT2  m1 nvt/asphere temp 300 300 1000
run     40000 post no
unfix NVT2
# Kicking the inserted molecule towards the surface.
velocity m1 set -1.23e-5 1.30e-5 -2.26e-6 sum yes
fix             NVE2  m1 nve/asphere
run     2000000

The plot from the above simulation shows that:

  • The molecule is inserted at a distance from the surface exceeding the cutoff of 14 Angs, and as a consequence, the sum of the forces acting on the molecule is zero.
  • At step 40000 a vertical velocity is added to push the molecule towards the surface.
  • At step 67000 and vertical position zcom=33.4, the distance between the molecule and the surface is roughly equal to the cutoff + skin, and the molecule begins interacting with the surface, as highlighted by the spike in the force acting on it.
  • The molecule bounces away from the surface.

The force acting on the inserted molecule does not go to zero after being reflected, even when it reaches the original insertion point (outside the cutoff distance). My understanding is that the check yes option forces the neighbour list not to be rebuilt until at least one atom has moved more than half the neighbor skin distance. But why is the force non-zero outside the cutoff distance?

I repeated the calculation by forcing the neighbour list to be rebuilt every 100 steps:

neigh_modify check no delay 0 every 100

and indeed the force acting on the incoming molecule goes to zero as the distance from the surface increases:

This is a snapshot of the system to give a feeling of its geometry:


Highlighted in red are the two beads closer to the inserted molecule at zcom=33.9, which are located at a distance of 18.9 and 22.1 Angs, well outside the cutoff distance.

What is the origin of the interaction with the surface when the molecule is outside its cutoff range but the neighbour list is not updated?
The force on the incoming molecule doesn’t depend on the spatial decomposition, as it is the same if the simulation runs on a single processor.

Ultimately, I am looking for a way to reduce the discontinuities due to the interplay between the cutoff and the particle’s anisotropy, without blowing the computational cost too much e.g. by using an insanely long cutoff radius.

Thank you so much for sharing your thoughts.

Using long-range electrostatics with the slab correction complicates the code path a lot. Can you try using a Wolf summation instead, see pair_style coul/cut command — LAMMPS documentation? This would tell us if the issue is in the neighbor list, or specific to using LRE.

Shooting from the hip, but would a pair-specific neighbor list inclusion criterion be helpful for cutting down costs? I.e. if your specific pair class had a method that told npair whether or not to include two particles in the nlist. I’m not familiar with typical neighboring algorithms for highly anisotropic interactions, but maybe there’s then some elliptic or bounding box criterion that’d be more stingy?

If so, maybe the changes in this PR could be generalized? Add npair classes for a general pairwise cutoff by jtclemm · Pull Request #4060 · lammps/lammps · GitHub

Hi Stan,
I cannot try the Wolf summation out of the box, as the pair_style molc currently supports only the /long and /cut variants. It is not difficult to add more variants, as it is just a blend of the conventional gayberne and coul styles, plus some extra functions to link the position of the charges to the parent ellipsoid. However, I tried the molc/cut variant with:

variable cutoff equal  14.
neighbor 5. bin
neigh_modify check yes

and found that the force acting on the incoming molecule is indeed non-zero only when there are molecules inside the cutoff radius:


Maybe the problem is that the neighbour list is very coarse, as it is computed on only 5 points per molecule. If I were to use an all-atom force field for the same molecule, I would have 78 atoms (and neighbour lists), yielding a much finer breakdown of intermolecular interactions.

@jtclemm: my example would benefit from your PR, as it uses particles with significant anisotropy. I also found that a larger skin distance is necessary to avoid lost particles, especially when the particle’s size is comparable with the cutoff distance.

Thanks, that’s good to know. I’m debating simplifying that PR and removing the new npair classes, but if they have other use cases maybe it’s best to keep it general.