I am running a simulation where a large portion of my object is fixed (i.e. fix setforce 0 0 0). To speed up my calculations, I attempted to exclude the fixed atoms from the pair list (i.e. neigh_modify exclude). However, when I added this command, the simulation did not arrive at the same result as before. After a bit of experimentation, I discovered that the forces on the unfixed atoms were incorrect after adding the neigh_modify command. However, the energy was still correct after the neigh_modify command was added.
After seeing that the neigh_modify command had this unexpected behavior with an EAM potential, I tested the same script with an LJ potential and both the forces and the energy were correct. This leads me to suspect that there is something I don’t understand about how LAMMPS interacts with EAM potentials or about EAM potentials in general. Does anyone know what would cause this?
Here is an MWE that reproduces what I am seeing. Attached at the bottom is the LAMMPS atom data file.
# LAMMPS input file that keeps all atoms but one fixed without neigh_modify exclude
# Nothing really moves here, I am just checking the force
units metal
dimension 3
boundary s s s
atom_style atomic
atom_modify map array
read_data lammps.data
pair_style eam/alloy
pair_coeff * * Au_Zhou04.eam.alloy Au
neighbor 100.0 bin
neigh_modify one 50000 delay 10 check yes page 500000
# I have 5122 atoms, generated semi-randomly with atomsk
group fixed id <= 5121
group unfixed id > 5121
fix noForce fixed setforce 0.0 0.0 0.0
compute eng unfixed pe/atom
compute x unfixed property/atom x
compute y unfixed property/atom y
compute z unfixed property/atom z
compute fx unfixed property/atom fx
compute fy unfixed property/atom fy
compute fz unfixed property/atom fz
compute id all property/atom id
variable eng equal c_eng[5122]
variable x equal c_x[5122]
variable y equal c_y[5122]
variable z equal c_z[5122]
variable fx equal c_fx[5122]
variable fy equal c_fy[5122]
variable fz equal c_fz[5122]
variable id equal c_id[5122]
thermo_style custom step v_eng v_x v_y v_z v_fx v_fy v_fz v_id
run 0
# LAMMPS input file that keeps all atoms but one fixed with neigh_modify exclude
# Nothing really moves here, I am just checking the force
units metal
dimension 3
boundary s s s
atom_style atomic
atom_modify map array
read_data lammps.data
pair_style eam/alloy
pair_coeff * * Au_Zhou04.eam.alloy Au
neighbor 2.0 bin
neigh_modify one 50000 delay 10 check yes page 500000
# I have 5122 atoms, generated semi-randomly with atomsk
group fixed id <= 5121
group unfixed id > 5121
fix noForce fixed setforce 0.0 0.0 0.0
neigh_modify exclude group fixed fixed
compute eng unfixed pe/atom
compute x unfixed property/atom x
compute y unfixed property/atom y
compute z unfixed property/atom z
compute fx unfixed property/atom fx
compute fy unfixed property/atom fy
compute fz unfixed property/atom fz
compute id all property/atom id
variable eng equal c_eng[5122]
variable x equal c_x[5122]
variable y equal c_y[5122]
variable z equal c_z[5122]
variable fx equal c_fx[5122]
variable fy equal c_fy[5122]
variable fz equal c_fz[5122]
variable id equal c_id[5122]
thermo_style custom step v_eng v_x v_y v_z v_fx v_fy v_fz v_id
run 0
Attached is the input LAMMPS data file. lammps.data (400.5 KB)
Why such a gigantic neighbor list skin value? You are massively blowing up your neighbor list with mostly entries that do not interact and thus slow down your calculation. Typical for metal units would be something in the range between 1.0 and 2.0.
The extreme neighbor list skin value is probably also the reason for requiring the modified values for “one” and “delay”. That is usually a sign for a very dense system. Also, delay 10, is a very aggressive value. You get about the same performance gain (after you lost much, much more due to the neighbor skin choice) with delay 2. Since you have “check yes” the neighbor list rebuilt will still be delayed if not needed, but the check done by check yes is very fast.
Yes. The EAM model contains a pairwise term and an “embedding term”. For the computation of the embedding term, you need to loop over all neighbors and collect a “density” contribution from the neighbors. So in EAM you loop over neighbors twice: first to collect the embedding density and derive from it the strength of the embedding term, and a second time to compute the pairwise forces and the forces from the embedding term.
Sorry about that, I was attempting to do a little debugging and thought that if all (or most) of the atoms were included in the neighbor list that it would give me more information about what was happening.
That is usually a sign for a very dense system.
Yes, in my case the “system” can be very dense. What I am actually doing is using LAMMPS in an optimization scheme, with LAMMPS keeping the system physically reasonable. The other end of the optimizer can often push atoms too close together. However, I am probably being overly cautious with my parameter choices here.
So in EAM you loop over neighbors twice:
I apologize if I seem a bit dense here, but I don’t understand how this would result in the forces being different when I have added in the neigh_modify command. Even though the EAM must loop over the neighbors twice, as long as they are the same neighbors, the EAM should spit out the same forces. Or do you mean that there are “neighbors” that affect it outside of its local atomic environment (i.e. something that is not in the neighbor list)?
In general you should expect three-body forces to depend on neighbors of neighbors unless proven otherwise. Here’s how this works for EAM specifically, based on my two-minute skim of the EAM source code:
Let’s say we’re calculating the force between neighbor atoms 1 and 2.
Let’s say some atom 3 would normally be atom 2’s neighbor, but is not thanks to neigh_modify delete.
The EAM potential includes an embedding term due to the electron density on atom 2, “rho[2]”.
This rho[2] would be calculated by looping over all neighbors of atom 2. So normally atom 3 would contribute to rho[2]. But with your neigh_modify setting, it does not.
Therefore the neigh_modify setting changes rho[2]. And since the atom1-atom2 force depends on rho[2], it is also changed.
As an exercise, work through the same logic to see that a pairwise Lennard-Jones potential would not exhibit similar behaviour, and the pairwise forces (wherever evaluated) would not be modified by neigh_modify.