Why are per-atom forces in dump different from atom->f in pair::compute

LAMMPS LAMMPS Beginners
1

Dear LAMMPS users,

I am trying to run some simulations in LAMMPS and I want to look at per-atom forces at each timestep.

To do so, I included the following line in my input file:
dump 1 all custom 1 dump x y z fx fy fz

The values I got out from the dump command looked weird, so I tried to modify the compute function of the potential I’m using in a way such that before it returns, it loops over all atoms and prints the value of atom->f[i][0], atom->f[i][1], atom->f[i][2] for each i.

The values that get printed from the code inside the modified compute are wildly different from the ones I see in the file produced by the dump command.

I don’t understand why this happens. Why do the values differ?
Is there maybe some other function that gets called before the execution of the dump command that modifies the contents of atom->f? If so, what function is it and what is the purpose of such a function?

I hope you can help me clarify this situation.
Thank you in advance,
Lorenzo.

2

You are mistaken. If there was some kind of “force transformation” it would be documented and discussed in the publications about LAMMPS. It would also be extremely unusual for a classical MD code to do such a thing (which would be even more of a reason to discuss and publish it).

The probable issue here is that you are comparing forces for the wrong atoms.
The index i in atom->f[i][0] is different from the atom ID used when writing out dump files.
You can see/print the atom ID for atom at index i at atom->tag[i].
This differentiation between local index and global atom ID is required because of the data distribution in LAMMPS via domain decomposition. Each MPI process has a list of atoms from 0 to atom->nlocal, which are the “owned” atoms of an MPI process and then also atoms beyond that, which are “ghosts” or copies of atoms owned by other MPI ranks, so that interactions with atoms in neighboring subdomains can be computed.

3

Dear Axel, thank you so much for your response.
I apologize for wasting your time with such matters, because I’m sure I’m missing something trivial here, but even after I made sure I’m looking at the same atoms from both inside the code and in the dump file, I observe that the forces differ.

I am using the serial version of LAMMPS to avoid any problems related to domain decomposition and following your observations I modified the instructions that I added to the PairEAM::compute (I am using the EAM potential).
Now, for i ranging from 0 to ilist->inum - 1, I print the following values:
atom->tag[i], atom->x[i][0], atom->x[i][1], atom->x[i][2], atom->f[i][0], atom->f[i][1], atom->f[i][2].

I also modified my dump command to this one:
dump 1 all custom 1 dump id x y z fx fy fz

From what I understand, the output produced by my printf statements inside the potential and the one produced by the dump command should be the same, however this is not the case.
Specifically, both show the same values for the first 4 fields (atom id and position), confirming that I am looking at the same atoms in both outputs, but the values for forces differ: they are always basically 0 in the dump file (e.g. 2.10812e-15), and of course in the output produced by the printf statements I added they are not always 0.

What could be the reason for this?

4

Impossible to say without seeing your input.

Also, please provide some examples showing how different those values actually are.

5

I modified the PairEAM::compute function by adding the following code at the end of it:
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
printf("%d %f %f %f %f %f %f\n", atom->tag[i], xtmp, ytmp, ztmp, f[i][0], f[i][1], f[i][2]);
}

Here is a minimal input file which shows the issue, along with the LAMMPS log and the dump file produced by its execution. (sorry for the messy text but I’m a new user, so I can’t upload files)
And here you can find the potential I’m using, even though it shouldn’t matter.

Input file:

units metal
boundary p p p

lattice bcc 3.0
region box block 0 2 0 2 0 2 units lattice
create_box 1 box
create_atoms 1 box

pair_style eam/fs
pair_coeff * * /home/lorenzo/Desktop/lmp/potentials/Fe_Earth_core.eam.fs Fe

neighbor 2.0 bin
neigh_modify delay 0 check yes page 100000

dump 1 all custom 1 dump id x y z fx fy fz
run 0

Dump:

ITEM: TIMESTEP
0
ITEM: NUMBER OF ATOMS
16
ITEM: BOX BOUNDS pp pp pp
0.0000000000000000e+00 6.0000000000000000e+00
0.0000000000000000e+00 6.0000000000000000e+00
0.0000000000000000e+00 6.0000000000000000e+00
ITEM: ATOMS id x y z fx fy fz
1 0 0 0 0 0 2.77556e-17
2 1.5 1.5 1.5 0 -1.11022e-16 -1.11022e-16
3 3 0 0 -2.22045e-16 -2.77556e-16 5.55112e-17
4 4.5 1.5 1.5 0 -1.11022e-16 -4.44089e-16
5 0 3 0 0 0 2.22045e-16
6 1.5 4.5 1.5 -1.11022e-16 0 -3.05311e-16
7 3 3 0 0 0 3.33067e-16
8 4.5 4.5 1.5 -5.55112e-17 5.55112e-17 -2.22045e-16
9 0 0 3 1.11022e-16 -1.11022e-16 -4.44089e-16
10 1.5 1.5 4.5 0 1.11022e-16 2.22045e-16
11 3 0 3 -2.22045e-16 0 -8.32667e-17
12 4.5 1.5 4.5 5.55112e-17 -1.11022e-16 3.33067e-16
13 0 3 3 4.44089e-16 0 0
14 1.5 4.5 4.5 0 -1.11022e-16 1.11022e-16
15 3 3 3 0 0 1.11022e-16
16 4.5 4.5 4.5 8.32667e-17 8.32667e-17 1.94289e-16

LAMMPS output:

LAMMPS (29 Oct 2020)
Lattice spacing in x,y,z = 3.0000000 3.0000000 3.0000000
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (6.0000000 6.0000000 6.0000000)
1 by 1 by 1 MPI processor grid
Created 16 atoms
create_atoms CPU = 0.000 seconds
WARNING: No fixes defined, atoms won’t move (…/verlet.cpp:54)
Neighbor list info …
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8
ghost atom cutoff = 8
binsize = 4, bins = 2 2 2
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair eam/fs, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d/newton
bin: standard
Setting up Verlet run …
Unit style : metal
Current step : 0
Time step : 0.001
1 0.000000 0.000000 0.000000 -0.273007 -0.273007 -0.273007
2 1.500000 1.500000 1.500000 0.594028 0.594028 0.594028
3 3.000000 0.000000 0.000000 0.167102 -0.320763 -0.320763
4 4.500000 1.500000 1.500000 0.000000 0.800601 0.800601
5 0.000000 3.000000 0.000000 -0.583036 0.540776 -0.216852
6 1.500000 4.500000 1.500000 0.259825 0.317974 0.848811
7 3.000000 3.000000 0.000000 0.334204 0.334204 -0.630338
8 4.500000 4.500000 1.500000 -0.000000 0.000000 1.055383
9 0.000000 0.000000 3.000000 -0.583036 -0.583036 1.129762
10 1.500000 1.500000 4.500000 0.259825 0.259825 0.684158
11 3.000000 0.000000 3.000000 0.334204 -0.996521 0.923190
12 4.500000 1.500000 4.500000 0.000000 0.466397 0.366184
13 0.000000 3.000000 3.000000 -1.314495 0.874980 0.874980
14 1.500000 4.500000 4.500000 0.092723 0.262273 0.262273
15 3.000000 3.000000 3.000000 0.668407 0.668407 0.668407
16 4.500000 4.500000 4.500000 0.000000 0.000000 0.000000
Per MPI rank memory allocation (min/avg/max) = 4.479 | 4.479 | 4.479 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -61.414831 0 -61.414831 -234319.22
Loop time of 0 on 1 procs for 0 steps with 16 atoms

0.0% CPU use with 1 MPI tasks x no OpenMP threads

MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total

Pair | 0 | 0 | 0 | 0.0 | 0.00
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0 | 0 | 0 | 0.0 | 0.00
Output | 0 | 0 | 0 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0 | | | 0.00

Nlocal: 16.0000 ave 16 max 16 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 839.000 ave 839 max 839 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 1344.00 ave 1344 max 1344 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 1344
Ave neighs/atom = 84.000000
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:00

Thank you again for taking the time to respond.

6

Please note this indicates that you are using a “half” neighbor list, which means for pairs of atoms that cross subdomain boundaries (which is equivalent to periodic boundaries in the case of a serial calculation) some of the forces are tallied to “ghost” atoms, i.e. atoms that “belong” to a different subdomain or are periodic copies (LAMMPS is not subject to minimum image conventions it has real copies of the periodic replica up to the ghost atom communication cutoff). So you don’t have the final forces on the “owned” atom until after you have a reverse communication. See the flow of control in the pseudocode here: 4.3. How a timestep works — LAMMPS documentation

You should compare the result to a run with the command “newton off” added at the beginning of the input (newton command — LAMMPS documentation) where those interactions are computed twice and the resulting force only tallied to “owned” atoms and then the reverse communication is avoided.

1 Like