Numerical forces in lammps: implemented ?

LAMMPS LAMMPS Mailing List Mirror
1

Dear Mailing list,

I would like to test numerically the forces of a new potential that I am implementing in Lammps.

With that I mean to compare the formulas with the forces obtained by making the gradient of the total energy.

My questions are:

  1. Is this feature implemented already in LAMMPS ?

  2. If so, what keywords should I look at in the manual ? (I looked in the database of the mailing list but couldn’t find any clue)

  3. If not implemented: do you have any advice on how to that in a general way? What files should I handle ?

Thank you very much for your support

Best regards,
Daniele Scopece

2

Dear Mailing list,

I would like to test numerically the forces of a new potential that I am
implementing in Lammps.
With that I mean to compare the formulas with the forces obtained by making
the gradient of the total energy.

My questions are:
1) Is this feature implemented already in LAMMPS ?
2) If so, what keywords should I look at in the manual ? (I looked in the
database of the mailing list but couldn't find any clue)
3) If not implemented: do you have any advice on how to that in a general
way? What files should I handle ?

you can do potential energy scans using loops and fix move and then
just record the energy with all velocities set to 0.0.

for pairwise additive potentials, you can also use the pair_write
command (provided you implement your model also with the
Pair::single() method) and compare the resulting table against the
formula, and similarly, you could compare against the table pair
style.

axel.

3

If you have written a pair style with
energy and forces, and are asking how to
check that it is correct, one simple test

is to run a few timesteps. If it does not

conserve energy, that is often an indication
that the energy and forces are not self-consistent.

I.e. the force is not the negative gradient of the
energy. Ditto for a minimization. If you can’t converge

to zero force = minimum energy, then the

forces and energy are not consistent.

Steve

4

Dear Steve and Axel,

Thank you for your replies.

However I have another question
that has arisen while working on this,
that is actually a question on the generic way LAMMPS handles forces and the atoms’ replicas.

(Hoping that it will be useful also to others working on potential development, if it is not obvious).

For simplicity I started to compare the “analytical” forces (computed in my 2compute" method)
with the “numerical” ones (coming from the gradient of the total energy)
for an isolated system (i.e. with such a large cell that the atoms are not interacting with any of their replica).

The agreement was fine.

Then I reduced the size of the cells so to have a 3D bulk.

If I compare the forces as computed by the “compute” method in my potential with the numerical,
the former turned out not to be zero (that is unphysical).

I compared the same also with pair_tersoff.cpp and, with my great surprise, the same results came out.

But, in applying a CG the Fmax registered there is effectively zero (idem in my potential).

From this and other inspections,
I suspect that the total forces experienced by the atoms are not only the ones computed in the “compute” method in the pair_XXX.cpp file,
but that the interactions with the images are taken into account also somewhere else in the code (I haven’t studied where yet).

My questions are:

  1. Is my understanding correct ?

  2. Where else the forces are modified in LAMMPS (to take the correct interactions with the images into account) ?

  3. Is there a command I can use to plot the forces after this “correction” in the output file, so to compare the “complete” forces with my numerical ones ?

I’d like to plot the forces on all atoms.

Thank you very much for the clarifications you will give me

(and also for your patience)

Best regards,

Daniele Scopece

5

Dear Steve and Axel,

Thank you for your replies.

However I have another question
that has arisen while working on this,
that is actually a question on the generic way LAMMPS handles forces and the atoms’ replicas.
(Hoping that it will be useful also to others working on potential development, if it is not obvious).

For simplicity I started to compare the “analytical” forces (computed in my 2compute" method)
with the “numerical” ones (coming from the gradient of the total energy)
for an isolated system (i.e. with such a large cell that the atoms are not interacting with any of their replica).
The agreement was fine.

Then I reduced the size of the cells so to have a 3D bulk.
If I compare the forces as computed by the “compute” method in my potential with the numerical,
the former turned out not to be zero (that is unphysical).
I compared the same also with pair_tersoff.cpp and, with my great surprise, the same results came out.
But, in applying a CG the Fmax registered there is effectively zero (idem in my potential).

From this and other inspections,
I suspect that the total forces experienced by the atoms are not only the ones computed in the “compute” method in the pair_XXX.cpp file,
but that the interactions with the images are taken into account also somewhere else in the code (I haven’t studied where yet).

My questions are:

  1. Is my understanding correct ?

No.

  1. Where else the forces are modified in LAMMPS (to take the correct interactions with the images into account) ?

Nowhere.

  1. Is there a command I can use to plot the forces after this “correction” in the output file, so to compare the “complete” forces with my numerical ones ?
    I’d like to plot the forces on all atoms.

You cannot. Buy you make claims without backing them up with tangible data. One cannot start a discussion this way. So please provide proper evidence that things are not the way they should be.

Axel.

6

Dear all,

I apologize if I didn’t give any detail, but I was frankly expecting a “YES”…

But probably I was wrong …

Here I provide all the details of my test, so that everything can be checked and looking for my possible error.

I am using lammps-1Feb14.

Attached you can find a folder with the test on Tersoff.

I am testing the forces on a simple and very symmetric system:

A simple cubic structure with lattice constant = 2.5 Angstrom.

I multiplied the cell with many atoms (512) to be sure that an atom was not interacting with its own images,
nor with the images of their own neighbors (to be 100% sure).

Sub-folder 01-* is the test on the forces:

  1. pair_tersoff.cpp = this is the original file where I simply added a print of the forces at lines 239 and 243-249 (with a loop equal to the one on ii in the same file)
  2. Test-cube.0.x4.y4.z4.data = is my cell data
  3. in.testcube.Ters = is my input file (maybe I should have inserted something else here? Did I make a mess?)
  4. RUN-Lammps-pavilion-here.v01.sh = my script to run the calculation on my laptop
  5. out-DS.testcube.Ters = the output obtained by lammps

The lines starting with ^F-FINAL are the final forces as printed from the compute method in the format:

F-FINAL: i_atom = f[i_atom][0] f[i_atom][1] f[i_atom][2]

From here we see that a force of ~1.60 is registered for those atoms at the borders of the cells
(the ones that are interacting with some image atoms).

In this same file at line 1037-1038 the output of the zero-th step of the CG are printed:

Fmax = 4E-16.

This same results can be found in sub-folder 02-* (where I run the full CG) in the file
out-DS.testcube.Ters at lines 1037-1038 and 4114.

I admit that I might be doing something wrong, but I cannot really understand what …

Maybe I put the print of the forces in pair_tersoff.cpp in the wrong position ? But it is at the end of compute() just before the computation of the virial (that shouldn’t modify forces, or does it?)

Maybe the cell is wrong ?
Maybe I should have added another command in my input file to consider something that was necessary ?

But why this inconsistency between the forces that I print from compute() and the one revealed by CG ?

I apologize in advance if this is annoying: I am just trying to understand …

I made also some tests with a Morse-like potential (written by me) in order to delucidate better what compute() was doing.

In this other simple case I had a 8-atoms simple cubic cell with PBC just along z (subfolder 3).

In this case the cutoff was such that the only pair of atoms that were interacting were the one on the same edge of the cube.

The force of atom 0 (lower corner of the cube, number 1 in the data file)
should be zero due to the compensation of the interactions due to atom 4
(just above atom 0, number 3 in the data file)
and the replica of 4 below 0 (labelled atom 68 by LAMMPS).

The cycles in the compute() of the potential are made with the standard for sum_{i != j}:
inum = list->inum;
ilist = list->ilist;
for (int inow = 0; inow<inum ; inow++){
i_atom = ilist[inow]; // local index of the atom
num_neigh_i = numneigh[i_atom]; // number of neighbors of atom i_atom (NOT inow)
jlist = firstneigh[i_atom]; // list of neighbors
for (int jnow=0; jnow<num_neigh_i ; jnow++){

// here j_atom aumatically different from i_atom since it is the list of i_atom

j_atom = jlist[jnow]; // local index of the atom
j_atom &= NEIGHMASK; // needed (?)

---- here computation of the morse-like potential and of forces on i_atom and j_atom ----

}
}

It turns out that the final forces along z on atom 0 is NOT null.

The interactions registered in this case are:

0-4 (i.e. i_atom = 0 ; j_atom = 4) + 0-68 (68 is in the list of 0),

these 2 cancels.

BUT also
4-0

that is not cancelled by
68-0

since 68 is apparently not in the ilist (since it is a replica ?)

Should I modify the loops or put a condition in them in order to get the right results ?

If it is needed I can provide the file, but probably in order to understand focusing on Tersoff might be enough.

Thank you in advance for your help and my anticipated apologize if I messed up anything.

Best regards,

Daniele Scopece

Test-forces-Tersoff-x-mailing-list.tar.bz2 (21 KB)

7

Dear all,

I apologize if I didn't give any detail, but I was frankly expecting a
"YES"...
But probably I was wrong ...

Here I provide all the details of my test, so that everything can be checked
and looking for my possible error.

I am using lammps-1Feb14.

Attached you can find a folder with the test on Tersoff.

I am testing the forces on a simple and very symmetric system:
A simple cubic structure with lattice constant = 2.5 Angstrom.
I multiplied the cell with many atoms (512) to be sure that an atom was not
interacting with its own images,
nor with the images of their own neighbors (to be 100% sure).

Sub-folder 01-* is the test on the forces:
1) pair_tersoff.cpp = this is the original file where I simply added a
print of the forces at lines 239 and 243-249 (with a loop equal to the one
on ii in the same file)
2) Test-cube.0.x4.y4.z4.data = is my cell data
3) in.testcube.Ters = is my input file (maybe I should have inserted
something else here? Did I make a mess?)
4) RUN-Lammps-pavilion-here.v01.sh = my script to run the calculation on my
laptop
5) out-DS.testcube.Ters = the output obtained by lammps
    The lines starting with ^F-FINAL are the final forces as printed from
the compute method in the format:
     F-FINAL: i_atom = f[i_atom][0] f[i_atom][1] f[i_atom][2]
From here we see that a force of ~1.60 is registered for those atoms at the
borders of the cells
(the ones that are interacting with some image atoms).
In this same file at line 1037-1038 the output of the zero-th step of the CG
are printed:
Fmax = 4E-16.

This same results can be found in sub-folder 02-* (where I run the full CG)
in the file
out-DS.testcube.Ters at lines 1037-1038 and 4114.

I admit that I might be doing something wrong, but I cannot really
understand what ...
Maybe I put the print of the forces in pair_tersoff.cpp in the wrong
position ? But it is at the end of compute() just before the computation of
the virial (that shouldn't modify forces, or does it?)

if you have newton's 3rd law for ghost atoms enabled (which is the
default), then the forces at the end of compute are not yet complete.
you are still missing the contribution from atom pairs that are
computed on other processors, which is communicated in the following
reverse communication step by the integrator class.

in any case, modifying the source code to output forces without
knowing LAMMPS overall flow of control is not a good idea to begin
with. the clean way of outputting forces is via the dump custom
command, which also doesn't require any modification of the code.

axel.

8

I’m reasonably confident you can ask LAMMPS
to build a single unit cell of Si, e.g. with the diamond
lattice, run it for 0 steps with Tersoff and dump the
forces on each atom. And they will all be 0.0.
That will involve lots of ghost atom communication and
force summation. Axel is correct that the forces won’t
be zero internal to the code until after all that communication
has occurred. They should be if you dump them,
since output is at the very end of the timestep.

Steve