Re: [lammps-users] pair_style nb3b/harmonic implementation vs documentation
Axel Kohlmeyer <akohlmey@...24...>
Fri, 30 Jun 2017 11:18:06 -0400
On Fri, Jun 30, 2017 at 11:02 AM, Mordechai Kornbluth <mkornbl@...24...> wrote:
> 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
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.  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
> The documentation  and a subsequent mailing-list discussion  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. 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  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.
> 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
> 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 .
> - Yet the commit message indicates only the changing of > to >=, not the
> deletion of the ~25 lines of code from pair_nb3b_harmonic.cpp .
> - Also, the deleted lines  don’t seem to be doing anything; perhaps
> in an earlier form they ensured the angle was not double-counted (point 2
> 4. 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
> - if (rsq1 > params[ijkparam].cutsq) continue;
> - if (rsq2 > params[ijkparam].cutsq) continue;
> Thanks in advance!
>  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.
>  http://lammps.sandia.gov/doc/pair_nb3b_harmonic.html
>  http://lammps.sandia.gov/threads/msg56857.html
>  Although theta_0 is given in degrees, it is apparent from the source
> code  l.386,440,483 that K is in energy/radian^2.
>  lines 117-130 here:
> Check out the vibrant tech community on one of the world's most
> engaging tech sites, Slashdot.org! http://sdm.link/slashdot
> lammps-users mailing list
Dr. Axel Kohlmeyer akohlmey@...24... http://goo.gl/1wk0
College of Science & Technology, Temple University, Philadelphia PA, USA
International Centre for Theoretical Physics, Trieste. Italy.