pair_style nb3b/harmonic implementation vs documentation

LAMMPS LAMMPS Mailing List Mirror
1

I apologize if this is not the correct place to post this.

The nb3b/harmonic pair style doesn’t seem to behave as expected or documented with respect to the input-file interpretation. [0] I would greatly appreciate it if someone would enlighten me whether I am misreading the code & documentation, or if this should be submitted as an issue to be fixed.

The documentation [1] and a subsequent mailing-list discussion [2] indicate a few major features of this pair_style syntax:

  1. It is OK if not all triplets (i,j,k) are listed in the potential file.
  2. The value of K corresponds to the energy of a displacement of the angle of 1 radian. [3]
  3. If atoms j and k are the same, the cutoff is used, and the physical parameters (K, theta_0) are totally ignored.
  4. If atoms j and k are different, the physical parameters (K, theta_0) are used and the cutoff is totally ignored

However, the source code [4] indicates the following:

  1. In setup_params, line 425, an error is raised if the potential file is missing an entry, against #1 above.
  • This was added in the commit of 30 Jul 2016, and seems to make the documentation now out-of-date.
  1. The i-j-k loop iterates over all triplets, meaning that the (i,j,k) angle will be counted twice: Once as (i,j,k) and once as (i,k,j). If so, the energy of a displacement of 1 radian would correspond to 2K due to this double-counting.
  2. Even if atoms j and k are the same,
  • The K, theta_0 parameters are read (in read_file, lines 386-387)
  • The i,j,k triplet is part of the loop (in compute, lines 103-125, over all atoms)
  • The parameters (K,theta_0) are applied (in threebody, lines 482-483)
  • This seems to be an undocumented effect of the commit on 15 Oct 2015, a mere 3 days after the mailing list discussion [2].
  • Yet the commit message indicates only the changing of > to >=, not the deletion of the ~25 lines of code from pair_nb3b_harmonic.cpp [4].
  • Also, the deleted lines [4] don’t seem to be doing anything; perhaps in an earlier form they ensured the angle was not double-counted (point 2 above).
  1. The implementation is in agreement with the documentation. However, perhaps the developers should think about implementing a cutoff for the (i,j,k) triplet; it entails merely adding the following two lines before the call to threebody on line 138, and of course modifying the documentation appropriately:
  • if (rsq1 > params[ijkparam].cutsq) continue;
  • if (rsq2 > params[ijkparam].cutsq) continue;

Thanks in advance!

[0] I am assuming a convention where the input “i j k” identifies an angle where “i” is the central atom, i.e. the j-k-i angle. This is the convention of the nb3b.harmonic file.
[1] http://lammps.sandia.gov/doc/pair_nb3b_harmonic.html
[2] http://lammps.sandia.gov/threads/msg56857.html
[3] Although theta_0 is given in degrees, it is apparent from the source code [3] l.386,440,483 that K is in energy/radian^2.
[4] https://github.com/lammps/lammps/blob/master/src/MANYBODY/pair_nb3b_harmonic.cpp
[5] lines 117-130 here: https://github.com/lammps/lammps/blob/aa6624c029fd5eb32aeec35ab41e79574ddf6a29/src/MANYBODY/pair_nb3b_harmonic.cpp

2

I apologize if this is not the correct place to post this.

it is OK to post your observations here.
an alternative would be to file an issue at
https://github.com/lammps/lammps/issues
since the latter is quite new and we are still refining our guidelines
which should be sent where.
thus for the time being we accept both.

i am copying the author of the code so you have an authoritative
source for the items, that i don't know about.

The nb3b/harmonic pair style doesn’t seem to behave as expected or
documented with respect to the input-file interpretation. [0] I would
greatly appreciate it if someone would enlighten me whether I am misreading
the code & documentation, or if this should be submitted as an issue to be
fixed.

The documentation [1] and a subsequent mailing-list discussion [2] indicate
a few major features of this pair_style syntax:
1. It is OK if not all triplets (i,j,k) are listed in the potential file.
2. The value of K corresponds to the energy of a displacement of the angle
of 1 radian. [3]
3. If atoms j and k are the same, the cutoff is used, and the physical
parameters (K, theta_0) are totally ignored.
4. If atoms j and k are different, the physical parameters (K, theta_0) are
used and the cutoff is totally ignored

However, the source code [4] indicates the following:
1. In setup_params, line 425, an error is raised if the potential file is
missing an entry, against #1 above.
    - This was added in the commit of 30 Jul 2016, and seems to make the
documentation now out-of-date.

you are correct. i observed an issue, where out-of-bounds array
accesses would happen and suggested two options to resolve this. the
chosen resolution was to require a complete set of entries in the
potential file by adding dummy entries. this should indeed by
documented. i will add this.

i have to defer to answers to the rest of your questions to others.

axel.

3

I apologize if this is not the correct place to post this.

The nb3b/harmonic pair style doesn’t seem to behave as expected or
documented with respect to the input-file interpretation. [0] I would
greatly appreciate it if someone would enlighten me whether I am
misreading the code & documentation, or if this should be submitted as
an issue to be fixed.

The documentation [1] and a subsequent mailing-list discussion [2]
indicate a few major features of this pair_style syntax:
1. It is OK if not all triplets (i,j,k) are listed in the potential file.
2. The value of K corresponds to the energy of a displacement of the
angle of 1 radian. [3]
3. If atoms j and k are the same, the cutoff is used, and the physical
parameters (K, theta_0) are totally ignored.
4. If atoms j and k are different, the physical parameters (K,
theta_0) are used and the cutoff is totally ignored

However, the source code [4] indicates the following:
1. In setup_params, line 425, an error is raised if the potential file
is missing an entry, against #1 above.
    - This was added in the commit of 30 Jul 2016, and seems to make
the documentation now out-of-date.
2. The i-j-k loop iterates over all triplets, meaning that the (i,j,k)
angle will be counted twice: Once as (i,j,k) and once as (i,k,j). If
so, the energy of a displacement of 1 radian would correspond to 2K
due to this double-counting.
3. Even if atoms j and k are the same,
    - The K, theta_0 parameters are read (in read_file, lines 386-387)
    - The i,j,k triplet is part of the loop (in compute, lines
103-125, over all atoms)
    - The parameters (K,theta_0) are applied (in threebody, lines
482-483)
    - This seems to be an undocumented effect of the commit on 15 Oct
2015, a mere 3 days after the mailing list discussion [2].
    - Yet the commit message indicates only the changing of > to >=,
not the deletion of the ~25 lines of code from pair_nb3b_harmonic.cpp
[4].
    - Also, the deleted lines [4] don’t seem to be doing anything;
perhaps in an earlier form they ensured the angle was not
double-counted (point 2 above).

I started to work with this potential (in an hybrid scheme) and I
confirm that the 3b interaction is effectively taken into account even
if atom j and k are the same.
I'm however not sure about many other aspect, like the units or how the
cutoff is considered... I need to look more carefully at the code.

Do you have any news on your side on the discrepancies between the doc
and the code? (except the comments from Axel on the 30.06)

4

Also, the author of pair nb3b/harmonic is listed at the top

of the source file. I suggest you contact him if you

have specific Qs: tzeitle at sandia.gov

Steve